\section{Introduction}

Ultraluminous X-ray sources (ULXs) are point-like objects with X-ray luminosities above the Eddington luminosity of stellar-mass objects.
The origin of ULXs was believed to be either stellar-mass black holes with mass accretion above the Eddington limit or intermediate-mass black holes with mass accretion rates below the Eddington limit \citep{kaaret2017}.
Nevertheless, X-ray pulses have been detected associated with the ULX, M82 X-2, indicating that neutron stars could also be a source of ULX radiation \citep{bachetti2014}.
In recent years, several ULXs showing X-ray pulses have been identified and it has been widely recognised that many ULXs likely originate from neutron stars \citep{furst2016, israel2017a, israel2017b, tsygankov2017, carpano2018, doroshenko2018, rodriguez2020}.
Additionally, there is a ULX that is thought to have a neutron star origin because of the presence of a strong magnetic field, although no associated X-ray pulse has been detected \citep{brightman2018}.
Even though there are an increasing number of observations of neutron stars with luminosities above the Eddington limit, the mechanism of the X-ray luminosity and the magnetic field strength of such luminous neutron stars remain under debate \citep{eksi2015, tong2015, mushtukov2015, king2019, erkut2020, mushtukov2021}.

The fraction of ULXs for which neutron stars are the X-ray source is also an interesting issue \citep{fragos2015, shao2015, misra2020, kuranov2020}.
To estimate the fraction of neutron star ULXs, it is important to understand the binary nature of ULXs with accreting neutron stars.
From the current limited number of observations, it appears that neutron star ULXs are a type of high mass X-ray binary (HMXB) with a massive companion \citep{karino2016, king2019}.

HMXBs are traditionally classified as either OB-type or Be-type, depending on the type of donor.
If the system has an OB-type star and a short orbital period, the mass transfer from the donor to the neutron star will be via Roche lobe overflow (RLOF), which provides a large mass transfer rate \citep{bildsten1997, frank2002}.
In this case, mass accretion above the Eddington limit occurs in the neutron star and the X-ray luminosity will exceed the Eddington luminosity if the outflow around the neutron star does not overload it too much \citep{lipunova1999, kosec2018, chashkina2019}.
Note that nearby OB-type HMXBs often accrete via strong stellar wind from the OB donors; nevertheless, the stellar wind accretion is unlikely to exceed the Eddington magnitude \citep{karino2016}.

Conversely, if the donor is a Be-type star, the neutron star accretes a large amount of mass as it passes through the circumstellar disc of the Be star during its orbital motion and emits X-rays for a short duration.
Consequently, Be-type high mass X-ray binaries (BeHMXBs) produce X-ray bursts (Type I bursts) with a period related to their orbital period.
Sometimes, regardless of the orbital period, BeHMXBs also produce giant bursts (Type II bursts) that are brighter than normal bursts \citep{reig2011, rivinius2013}.
Some neutron star ULXs appear to have a Be-type companion; thus, it is an interesting question as to whether BeHMXB systems can achieve X-ray luminosities above the Eddington luminosity \citep{okazaki2013, karino2016, karino2021}.
Table 1 lists the neutron star ULXs that are suggested to be accreted from Be-type donors.

Generally, the X-ray luminosity of BeHMXBs is smaller than that of X-ray binary systems with OB-type donors (the X-ray luminosity in a typical Type I burst is approximately $10^{36} \rm{erg \, s}^{-1}$).
However, BeHMXBs constitute a large fraction of HMXBs \citep{bildsten1997, coe2015}.
Thus, if neutron star ULXs are part of the evolutionary process of BeHMXBs, there may be a significant number of neutron star ULXs (even if no pulses are detected) \citep{inoue2020, king2020, kuranov2020, mushtukov2021}.
The neutron star/black hole fraction in ULXs is not only important concerning the initial mass function of massive stars but is also interesting for estimating the population of gravitational wave sources such as double neutron star binaries.

Given this motivation, to investigate the characteristics of neutron star ULXs, \citet{karino2021} considered two conditions under which a BeHMXB can become a ULX: (1) the mass transfer rate from the Be donor to the neutron star exceeds the Eddington rate and (2) the transferred material accretes to the neutron star at a rate exceeding the Eddington rate.
In this study, we adopt the above conditions to the evolution of binary systems and investigate the conditions under which BeHMXBs can become ULXs.
Particularly, we consider the length of time that binary systems can be super-Eddington sources in the evolutionary process and investigate their observability.

In the next section, we discuss the conditions under which a BeHMXB can become a ULX.
In Section 3, we introduce an evolutionary model of binary systems.
In Section 4, as a result of numerical calculations, we determine the period during which the X-ray luminosity of a BeHMXB exceeds the Eddington luminosity.
We discuss the characteristics of binary systems that become ULXs on the basis of the conditions presented in Section 2, and we consider their observability.
In Section 5, we compare the observed binary systems with the results of the model calculations.
Additionally, we discuss the limitations of the models and future modelling prospects.
The final section is devoted to our conclusions.

\begin{table*}
	\centering
	\caption{\label{table1}
List of neutron star ultraluminous X-ray sources (ULXs) suggested to have Be-type donors.
For each systems, maximum luminosities, spin periods, orbital periods, donor masses, and orbital eccentricities are shown.
	}
	\label{tab:Table}
	\begin{tabular}{lcccccl}
		\hline
		Object & $L_{\rm{X, max}}\ \rm{[erg \, s^{-1}]}$ & $P_{\rm{s}}\ \rm{[s]}$ & $P_{\rm{orb}}\ \rm{[d]}$ & $M_{\rm{d}} \ \rm{[M_{\odot}]}$ & $e$	& Ref. \\
		\hline

		SMC X-3 & $3.5 \times 10^{39}$ & 7.8 & 45 & ${>}3.7$ & $0.23-0.26$ & \citet{weng2017,tsygankov2017} \\
		Swift J0243.6+6124 & $2.6 \times 10^{39}$ & 9.9 & $27.6/28.3$ &  & $0.1$ & \citet{wilson2018, doroshenko2018}\\

		RXJ0209.6-7427 & $1.6 \times 10^{39}$ & 9.3 & $39/47$ & & & \citet{chandra2020}\\

		\hline
	\end{tabular}
	\\

\end{table*}

\section{Conditions under which a BeHMXB burst exceeds the Eddington luminosity}

In a previous study, \citet{okazaki2013} investigated the conditions under which a BeHMXB causes a giant burst.
The conditions under which the luminosity of such giant bursts exceeds the Eddington luminosity can also be studied to determine the conditions that allow a BeHMXB to become a ULX \citep{karino2021}.

The following has been shown by \citet{okazaki2013} in the case of a massive burst of a neutron star supplied with mass via Bondi–Hoyle–Littleton (BHL)-type accretion from a warped disc around a Be donor.
According to the BHL theory, a neutron star captures material inside the accretion radius,
\begin{equation}
R_{\rm{acc}} = \frac{2 G M_{\rm{NS}}}{v_{\rm{rel}}^{2}},
\label{eq:racc}
\end{equation}
when it interacts with the disc material.
The mass transfer rate from the Be disc to the neutron star is given by
\begin{equation}
\dot{M}_{\rm{T}} = \min \left(1, \frac{R^2_{\rm{RL}}}{R^2_{\rm{acc}}} \right) \times \dot{M}_{0}.
\label{eq:MdotT}
\end{equation}
This equation implies that only material inside the Roche radius,
\begin{equation}
R_{\rm{RL}} = D \frac{0.49 q^{2/3}}{0.69 q^{2/3} + \ln (1 + q^{1/3}) },
\label{eq:Rrl}
\end{equation}
where the gravity of the neutron star is dominant, can be accreted.
Here, $D = a(1-e)$ is the separation of the binary at the periastron passage, where $a$ is the semi-major axis of the orbit, which can be obtained from Kepler's law by specifying the donor mass, $M_{\rm{d}}$, the mass of neutron star, $M_{\rm{NS}}$, the orbital period, $P_{\rm{orb}}$, .and the orbital eccentricity, $e$.
$\dot{M}_{0}$ is the mass transfer rate in the BHL accretion theory given by \citet{hoyle1939}, \citet{bondi1944} and \citet{edgar2004},
\begin{equation}
\dot{M}_{0} = \rho v_{\rm{rel}} \pi R_{\rm{acc}}^{2},
\label{eq:Mdot0}
\end{equation}
where $v_{\rm{rel}}$ is the relative velocity between the orbital motion of the neutron star and the rotational velocity of the Be disc material, which depends on the inclination of the disc surface and the orbital plane of the neutron star.
Here, the inclination angle between the disc and the orbital plane is represented by the angle between the orbital velocity vector and the Keplerian velocity vector of the disc material, $\delta$.
Then, the relative velocity can be written as
\begin{equation}
v_{\rm{rel}}^{2} = v_{\rm{orb}}^{2} + v_{\rm{Kep}}^{2} -2 v_{\rm{orb}} v_{\rm{Kep}} \cos \delta .
\label{eq:vrel}
\end{equation}
Here, $v_{\rm{orb}}$ is the orbital velocity at the periastron point and $v_{\rm{Kep}}$ is the Keplerian velocity of the disc material around the Be star at the periastron passage of the neutron star orbit.
If the Be disc and the orbital plane are aligned, the disc is truncated and becomes smaller.
Thus, a giant outburst may be more likely to occur if the disc is inclined to the orbital plane \citep{okazaki2013}.
Conversely, the inclination angle between the disc and the orbital plane also affects the mass transfer rate via the relative velocity between the orbital velocity of the neutron star and the rotational motion of the Be disc.
A smaller $\delta$ also implies a smaller relative velocity, which allows a larger mass transport rate as a result of BHL accretion and therefore a larger X-ray luminosity \citep{okazaki2013, karino2021}.
The dependence of the X-ray luminosity on the inclination is determined by a combination of these factors.

The density of the Be disc, $\rho$, is modelled here using the density of the inner edge of the disc, $\rho_{0}$, as follows \citep{rivinius2013}:
\begin{equation}
\rho = \rho_{0} \left(\frac{D}{R_{\rm{d}}} \right)^{-7/2},
\label{eq:rho}
\end{equation}
The density of the disc and the rotational velocity of the disc material vary with distance from the Be star; however, here all the physical quantities of the disc are represented by their values at the position $D$ of the periastron point.

The material transferred by the BHL process forms an accretion disc around the neutron star.
If the accretion disc is a standard disc, the accretion falls off slowly on a viscous timescale and does not produce a bright and sharp X-ray light curve as in Type II BeHMXB bursts.
Thus, the disc must be an advection-dominant, slim disc type that accretes quickly \citep{okazaki2013}.
To achieve this, the magnetic radius, $R_{\rm{mag}}$, of the neutron star must be smaller than the photon trapping radius, $R_{\rm{trap}}$.
This condition can be replaced by the condition that the mass transport rate $\dot{M}_{\rm{T}}$ be larger than a certain limit, $\dot{M}_{\rm{sc}}$:
\begin{equation}
\dot{M}_{\rm{T}} > \dot{M}_{\rm{sc}} = 1.7 \times 10^{19} k^{7/9} \mu_{30}^{4/9}.
\label{eq:Mdotsc}
\end{equation}
Here, the normalised magnetic moment of the neutron star is used:
\begin{equation}
\mu_{30} = \frac{B R_{\rm{NS}}^{3}}{10^{30}\ \rm{G \, \rm{cm}}^3},
\label{eq:mu}
\end{equation}
where $B$ is the surface field strength and $R_{\rm{NS}}$ denotes the radius of the neutron star.
$k$ is a parameter that depends on the accretionary geometry; here, we use $k = 0.5$ for the disc accretion.
The above condition, Eq.~(\ref{eq:Mdotsc}), must be satisfied for a BeHMXB to undergo a giant burst.

When Eq.~(\ref{eq:Mdotsc}) is satisfied, the accreting material transferred inwards through the slim disc is eventually trapped by the magnetic field of the neutron star forming an accretion column near the magnetic pole, which emits strong X-rays \citep{pringle1972}.
When the mass transfer exceeds the Eddington rate, the radiation pressure becomes dominant somewhere in the accretion disc and an outflow occurs \citep{lipunov1982, lipunova1999, poutanen2007}.
The radius at which the radiation pressure becomes dominant, $R_{\rm{sph}}$, can be written as
\begin{equation}
R_{\rm{sph}} = 1.4 \times 10^{6} \frac{\dot{M}_{T}}{\dot{M}_{\rm{Edd}}}\ \rm{cm},
\label{eq:rsph}
\end{equation}
where $\dot{M}_{\rm{Edd}}$ denotes the Eddington mass accretion rate \citep{king2017}.
Inside $R_{\rm{sph}}$, the accretion rate of the neutron star eventually decreases to
\begin{equation}
\dot{M}_{\rm{acc}} = \frac{R_{\rm{mag}}}{R_{\rm{sph}}} \dot{M}_{T}
\label{eq:MdotM}
\end{equation}
because of the loss of material via outflow.
Here, the magnetic radius of the neutron star is
\begin{equation}
R_{\rm{mag}} = \left[ \frac{\mu_{30}^{4} }{ 8 k^{2} G M_{\rm{NS}} \dot{M}_{\rm{acc}}^{2} } \right]^{1/7}.
\label{eq:rmag}
\end{equation}
Using the mass accretion rate, $\dot{M}_{\rm{acc}}$, given above, the luminosity of a giant burst of a BeHMXB can be obtained such that
\begin{equation}
L_{\rm{X}} = \dot{M}_{\rm{acc}} c^2,
\label{eq:LX}
\end{equation}
and, if this luminosity is greater than the Eddington luminosity, the system could be a candidate for a ULX.

Note that, in Eq.~(\ref{eq:MdotM}), because $R_{\rm{sph}} \propto \dot{M}_{0}$, the dependence of the mass accretion rate on the mass transfer rate vanishes and ultimately becomes
\begin{equation}
\dot{M}_{\rm{acc}} \propto B^{4/9}.
\label{eq:MB}
\end{equation}
This means that the luminosity of the neutron star-powered ULX only depends on the strength of the magnetic field of the accreting neutron star \citep{lipunov1982}.

Furthermore, previous studies suggest that a giant burst is more likely to occur when the Be disc, which is inclined to the orbital plane of the binary system, is warped by the tidal action of the neutron star \citep{negueruela2001, reig2007, moritani2011, okazaki2013}.
This is because, when the neutron star passes through the warped region, the mass transport rate increases, leading to a large outburst.

The conditions for the Be disc to warp can be expressed as \citep{karino2021}
\begin{equation}
R_{\rm{warp}} < R_{\rm{tr}}
\label{eq:warpcondition2}
\end{equation}
based on a combination of the typical radius of warp ($R_{\rm{warp}}$) induced by tidal action \citep{martin2011, rivinius2013} and the disc truncation radius ($R_{\rm{tr}}$) \citep{paczynski1977}.
These radii are expressed as
\begin{equation}
R_{\rm{tr}} = 0.5 a (1-e)
\label{eq:rtr}
\end{equation}
and
\begin{equation}
\begin{split}
R_{\rm{warp}} & = 4.91 \times 10^{11} \left(1- e^2 \right)^{3/4} \alpha_2^{1/2} h
\left(\frac{P_{\rm{orb}}}{1 \rm{d}} \right) \\
& \times
\left(\frac{M_{\rm{d}}}{M_{\odot}} \right)^{1/2}
\left(\frac{M_{\rm{NS}}}{M_{\odot}} \right)^{-1/2}
\left(\frac{M_{\rm{d}} + M_{\rm{NS}}}{M_{\odot}} \right)^{1/2}
\left(\frac{R_{\rm{d}}}{R_{\odot}} \right)^{-1/2}\ \rm{[cm]},
\end{split}
\label{eq:rwarp}
\end{equation}
respectively.
Here, $\alpha_{2}$ is the coefficient of the shear viscosity of the Be disc in the vertical direction and the value $\alpha_{2} = 2.73$ is adopted in \citet{martin2011}.
Conversely, \citet{cheng2014}, who investigated the presence of giant bursts in BeHMXBs depending on the warp conditions of the Be disc, uses $\alpha_{2} = 0.5 - 1$ to explain the bimodality of BeHMXBs on the Corbet diagram.
For large values of $\alpha_{2}$, most systems with long orbital periods and large eccentricities are rejected as ULX candidates \citep{karino2021}.
However, observations of BeHMXBs indicate that giant bursts occur even in systems with relatively large eccentricities.
For example, GROJ1008-57, which is a giant burst source, has a highly distorted orbit with an orbital period of 248 d and an orbital eccentricity of 0.65 \citep{coe2007, reynolds2020}.
To avoid rejecting such a system, $\alpha_{2}$ should be a small value.
Because $\alpha_{2}$ is highly indeterminate, it is given here for reference only; the value $\alpha_{2} = 0.5$ is used as a reference value.

In summary, a neutron star undergoing mass accretion from a Be donor to be a ULX must simultaneously satisfy Eq.~(\ref{eq:Mdotsc}) and
\begin{equation}
\dot{M}_{\rm{acc}} > \dot{M}_{\rm{Edd}}.
\label{eq:cond2}
\end{equation}
Additionally, we consider Eq.~(\ref{eq:warpcondition2}) as the warp condition for reference.

\section{Evolution of BeHMXBs and ULXs}

Here, we investigate how long a binary system could produce super-Eddington bursts during its evolutionary process by varying each BeHMXB parameter.
We fix the neutron star mass, $M_{\rm{NS}}$ and radius, $R_{\rm{NS}}$, at 1.4 solar masses and 10 km, respectively.
The neutron star magnetic field, $B$, is typically approximately $10^{12} \rm{G}$ in most of the systems observed via cyclotron resonant scattering and other methods \citep{christodoulou2017}; nevertheless, there is debate concerning the neutron star magnetic fields in ULXs \citep{eksi2015, tong2015, koliopanos2017, king2019, mushtukov2015, mushtukov2021}.
Here, we consider cases from $10^{11} \rm{G}$ to $10^{13} \rm{G}$.

B-type stars typically have masses close to 10 solar masses.
Here, we consider a Be donor mass, $M_{\rm{d}}$, of $8 - 16$ solar masses at the zero age main sequence stage.
The Be donor has already evolved by some amount at the moment of the supernova of the initial primary star and the birth of the BeHMXB.
Although the evolutionary stage of the donor depends on the mass of the initial primary star, here, we uniformly assume that the donor age is 5 Myr at the time of the birth of the BeHMXB.

Be stars eject their mass via stellar winds and/or deccretion discs.
The Be disc is considered to have a base density between $10^{-11} \rm{g \, cm}^{-3}$ and $10^{-10} \rm{g \, cm}^{-3}$, and the density distribution follows the power law given in Eq.~(\ref{eq:rho}) \citep{rivinius2013}.

The observed orbital periods of BeHMXBs, $P_{\rm{orb}}$, range from ten to several hundreds of days.
Additionally, the orbital eccentricity, $e$, has been observed to range from near circular to highly eccentric \citep{townsend2011}.
In this study, we consider the orbital period and the orbital eccentricity as being the main binary parameters and try various combinations of these parameters over a wide range.
There is little observational information concerning the parameter $\delta$, describing the inclination between the orbital plane and the Be disc plane.
Here, we compare three possible values of $\delta$: $\pi/16$, $\pi/8$ and $\pi/4$.

Under the initial conditions described above, we study the time evolution of a binary system comprising a Be donor and a neutron star.
The calculation of the binary evolution is based on the method of \citet{hurley2002}, with modifications to the code used in \citet{karino2020}.
For the donor Be star, however, we use a model with an age of 5 Myr at the start of the evolutionary calculation (using Eqs.(1)-(30) in \citet{hurley2000}).
The donor is assumed to be experiencing mass loss according to the mass loss rate of \citet{vink2001}, including ejection as a disc and loss via stellar wind.
Because of the short lifetime of binary systems, however, this mass loss rate does not have a significant effect on the subsequent evolution.
The binary period and the orbital eccentricity evolve as a result of the mass loss, mass transfer, gravitational wave emission and tidal action; however, they do not change significantly during the short lifetimes of HMXBs because of their large separation ($P_{\rm{orb}} {>}$ several days) \citep{reig2011}.
The base density of the disc and the inclination between the orbital plane and the disc plane are assumed to remain unchanged.
Under these conditions, the evolutionary effects are exclusively due to the evolution of the donor radius.

At the orbital periastron, if the Roche radius of the Be star is less than the radius of the Be star itself, mass transfer via RLOF from the donor occurs instead of BHL accretion from the disc.
In this case, a massive mass transfer may occur and the orbit may rapidly approach a circular orbit or, in some cases, a common envelope may be formed \citep{ivanova2013}.
When the binary mass ratio is large and the donor spin angular momentum is dominant, the binary separation shrinks rapidly due to tidal instability \citep{darwin1879,counselman1973,hut1980,taam2010,eggleton2006}.
Particularly, if the RLOF occurs after the donor has significantly evolved, it is unlikely to be observable as a stationary ULX.
In this study, we focus on the evolutionary stage when the BeHMXB becomes a (transient) ULX; accordingly, we terminate the calculation when the donor fills its Roche lobe.
If the donor does not fill its Roche lobe, the calculation is terminated when the donor evolves to the point where it causes a supernova explosion.
We solve for both the binary and donor evolution and, at the same time, consider the ULX conditions to identify the period when the BeHMXB would be considered a ULX.

\section{Results}

Choosing a set of BeHMXB initial parameters ($P_{\rm{orb}}$, $e$, $\rho_{0}$, $\delta$ and $B_{\rm{NS}}$), we can solve for the binary evolution.
At the same time, we consider the conditions under which the luminosity of the bursts exceeds the Eddington luminosity and calculate the length of time during which the BeHMXB system becomes a ULX.
Figure 1 shows colour contours indicating the length of time during which the burst luminosity exceeds the Eddington luminosity during the binary evolution.
In this case, the initial (zero age) donor mass is assumed to be 12 solar masses.
The horizontal axis indicates the initial orbital period of the binary, and the vertical axis indicates the initial orbital eccentricity.
The magnetic field of the neutron star, the base density of the Be disc and the inclination angle of the Be disc plane to the orbital plane are fixed at $B = 10^{12} \rm{G}$, $\rho_{0} = 5 \times 10^{-11} \rm{g \, cm}^{-3}$ and $\delta = \pi/8$, respectively.
Results are shown for evolutionary calculations of approximately 7500 binary systems with varying initial orbital periods and eccentricities.
The band-like area from the lower-left to the upper right of the figure is the region that satisfies the ULX condition during the binary evolution.
The parameter regimes that do not exceed the Eddington luminosity are shown in black.
Particularly, the black region in the upper left of the plot corresponds to the parameter regime in which the Roche lobe of the donor is filled at the beginning of the binary evolution.
Conversely, the lower right region of the plot mainly corresponds to a parameter regime in which the neutron star orbit does not pass through the high-density region of the Be disc and, therefore, sufficient mass transfer is not possible.

The figure indicates that a binary system with a period of less than approximately 20 d can be observed as a ULX for a certain period during its evolution regardless of the orbital eccentricity.
This result is consistent with the observational indication that Type II BeHMXB bursts are more likely to occur with shorter orbital periods \citep{reig2011}.
The length of time that a BeHMXB system could be a ULX during its evolutionary process depends on the orbital parameters but can be up to 10 Myr and is typically on the order of 1 Myr.
Because the lifetime of a 12 solar mass donor is roughly 18 Myr
, this time is not short in comparison with the lifetime of the HMXB.
This means that BeHMXBs with short orbital periods or large orbital eccentricities are expected to produce giant bursts that exceed the Eddington luminosity and could be observed as transient ULXs for a significant period of their lifetimes.
The duration of a ULX in the presence of stable mass transfer has been calculated by \citet{kuranov2020} to be from approximately 1 kyr–1 Myr.
Thus, the lifetimes of ULXs with Be donors could be longer than those of ULXs with super-giant donors.

Nevertheless, none of the observed BeHMXBs has binary periods shorter than 10 d.
This may be because such an orbital period is too short for the donor to supply enough mass to the disc or because the tidal action is too strong for the Be disc to form.
Hence, it is also expected that ULXs with Be donors will not have orbital periods of less than 10 d.
In Figure 1, the limit of $P_{\rm{orb}} = 10 \rm{d}$ is shown as a vertical line.
If we consider only the right side of this vertical line as a feasible region, we can see that the possible ULXs with Be donors are limited to systems with relatively large orbital eccentricities.

In the figure, the warp conditions of the Be disc, shown in Eq.~(\ref{eq:warpcondition2}), are indicated by two curves for reference.
This condition varies with the evolution of the donor.
The lower curve corresponds to the condition at the initial stage of the binary, i.e. when the donor age is 5 Myr.
The upper curve corresponds to the parameters when the donor reaches the terminal main sequence stage.
The lower parts of these curves indicate the parameter regime in which the Be disc warps inside the truncation radius.
Note that this condition is subject to the indefiniteness of the viscosity coefficient.
Because the viscosity parameter $\alpha_{2}$ is highly indeterminate, we have adopted $\alpha_{2} = 0.5$ as a likely value for warping
\footnote{
Even with the most optimistic viscous parameter ($\alpha_{2} = 0.5$), a Type II burst source with a long orbital period, i.e. GROJ1008-57, still does not satisfy the warp condition \citep{reynolds2020}.
Although it is necessary to account for the fact that the actual disc is more dynamically warped \cite[for example]{franchini2021}, we will consider the analytical warp condition as a reference.
}.
According to this line, systems with high orbital eccentricities and/or long orbital periods will not produce bright bursts.
Thus, possible ULXs with Be donors are restricted to a limited region in which the orbital period is less than roughly 100 d and the orbital eccentricity is smaller than approximately 0.8.
The observable lifetime of a system with these characteristics as a ULX is from approximately 0.1–10 Myr; however, this depends on the parameters.
Such a lifetime remains sufficiently long compared with the lifetimes of binary systems, such that the observability of objects in these parameter regimes as ULXs is ensured.
Additionally, as we will see later, this region contains many observed BeHMXBs.

\begin{figure}
\includegraphics[width=\columnwidth]{modified_m12_b12_r511_d8_warp_new.eps}
\caption{
Figure 1: Parameter regime in which a Be-type high mass X-ray binary (BeHMXB) undergoes a burst of super-Eddington luminosity.
The colour contours show the possible duration of the ultraluminous X-ray source (ULX) in years.
The horizontal axis indicates the initial orbital period of the binary, and the vertical axis indicates the initial orbital eccentricity.
The initial mass of the donor is 12 solar masses.
The parameters of the Be disc are fixed to $\rho_{0} = 5 \times 10^{-11} \rm{g \, cm}^{-3}$ and $\delta = \pi / 8$.
The vertical line indicates the lower limit of BeHMXB observations, a period of 10 d.
The two curves show the warp conditions at the beginning of the binary system (upper curve) and at the end of the donor main sequence (lower curve) for reference.
}
\label{fig:result1}
\end{figure}

The maximum possible burst luminosity of a BeHMXB is primarily determined via Eqs.~(\ref{eq:LX}) and (\ref{eq:MB}), which in the case of Figure 1 ($B = 10^{12} \rm{G}$) is approximately $L_{\rm{X}} = 16 L_{\rm{Edd}}$.
The magnetic field affects the maximum X-ray intensity but as only a weak contribution to the ULX condition via Eq.~(\ref{eq:Mdotsc}).
Thus, a variation in the magnetic field of one order of magnitude has little effect on the possible ULX duration.

Conversely, the parameters of the Be disc play an important role.
The density of the disc, $\rho_{0}$, directly affects the mass transfer rate, and the inclination of the disc, $\delta$, affects the mass capture rate in the BHL process.
Namely, a binary system is likely to be observed as a ULX when the disc density is high.
Additionally, a system is likely to be observed as a ULX when the inclination of the disc plane is small.
At the same time, the donor radius changes the size of the disc changes and affects ULX condition accordingly.

Figure 2 shows the possible super-Eddington luminosity bursting duration for different base densities, $\rho_{0}$, of the Be discs.
The other parameters are the same as those in Figure 1.
Figure 2 indicates that the density of the Be disc has a significant effect on ULX formation.
This implies that a Be disc with a low density of less than $10^{-11} \rm{g \, cm}^{-3}$ (see left panel) would not be observed as a ULX because it would not be able to achieve a sufficiently high mass transfer rate to become a ULX.
Conversely, when the base density of the disc increases to $10^{-10} \rm{g \, cm}^{-3}$ (right panel), super-Eddington bursts can be produced even in systems with longer orbital periods.
Based on observations of single Be star discs, Be stars with high-density discs of $5 \times 10^{-11} \rm{g \, cm}^{-3}$ are not rare \citep{touhami2011, rivinius2013, draper2014}.
Thus, it does not seem surprising that even a BeHMXB with typical parameters could become a ULX during its evolution.
When the density of the Be disc is high, the possible ULX duration is extended for systems with large eccentricities.
However, it is not clear whether this is actually observed because the warp condition is more difficult to satisfy for such systems.

\begin{figure*}
\includegraphics[width=6.0cm]{modified_m12_b12_r111_d8_warp_line_new.eps}
\includegraphics[width=5.2cm]{modified_m12_b12_r511_d8_warp_line_new.eps}
\includegraphics[width=6.0cm]{modified_m12_b12_r110_d8_warp_line_new.eps}
\caption{
Figure 2: Possible super-Eddington luminosity bursting duration for different base densities of Be discs.
The other parameters are the same as those in Figure 1.
The base densities are, from left to right, $1 \times 10^{-11} \rm{g \, cm}^{-3}$, $5 \times 10^{-11} \rm{g \, cm}^{-3}$ and $10 \times 10^{-11} \rm{g \, cm}^{-3}$.
}
\label{fig:result2}
\end{figure*}

Figure 3 shows the results when the inclination angle $\delta$ is varied from $\pi/16$ to $\pi/4$ to determine the effect of the misalignment between the disc plane and the orbital plane.
The other parameters are the same as those in Figure 1.
The figure indicates that the angle of the orbit also affects the formation of ULXs.
Additionally, the duration of the possible ULX increases when the orbital inclination angle decreases (left panel).
However, this change is primarily observed in the short-period region with a period of less than 10 d and may not be observed in real systems.
Conversely, when the orbital inclination angle is small, the possible ULX region is somewhat wider.
When the disc plane and the orbital plane are close to each other, relatively long-period binary systems with orbital periods of approximately 50 d and not so large orbital eccentricities ($e < 0.5$) can become ULXs for several millions of years.

\begin{figure*}
\includegraphics[width=6.0cm]{modified_m12_b12_r511_d16_warp_line_new.eps}
\includegraphics[width=5.2cm]{modified_m12_b12_r511_d8_warp_line_new.eps}
\includegraphics[width=6.0cm]{modified_m12_b12_r511_d4_warp_line_new.eps}
\caption{
Figure 3: Possible super-Eddington luminosity bursting duration for different inclination angles of the Be disc to the orbital plane.
The other parameters are the same as those in Figure 1.
The inclination angles are, from left to right, $\pi/16$, $\pi/8$ and $\pi/4$.
}
\label{fig:result3}
\end{figure*}

Figure 4 shows the results when the donor mass is varied from $8 M_{\odot}$ to $16 M_{\odot}$.
Again, the other parameters are the same as those in Figure 1.
The difference in the donor mass results in a remarkable difference in the duration of the super-Eddington bursts.
When the donor mass is low (left panel), the lifetime of the ULX increases significantly.
To a certain extent, this result is natural because the lifetime of the donor is significantly extended when its mass is small.
Nevertheless, if the donor mass is large (right panel), the donor radius is also large, such that the dense Be disc extends to the outer edge of the system, which would allow neutron stars to pass through the dense region of the disc even in long-period systems.
According to the computed results, however, this effect is limited.
Instead, the parameter range for giant bursts may become rather narrow because the conditions for the Be disc to warp are more difficult to attain for larger donor masses.

\begin{figure*}
\includegraphics[width=6.0cm]{modified_m08_b12_r511_d8_warp_line_new.eps}
\includegraphics[width=5.2cm]{modified_m12_b12_r511_d8_warp_line_new.eps}
\includegraphics[width=6.0cm]{modified_m16_b12_r511_d8_warp_line_new.eps}
\caption{
Figure 4: Possible super-Eddington luminosity bursting duration for different donor masses.
The other parameters are the same as those in Figure 1.
The donor masses are, from left to right, 8, 12 and $16 M_{\odot}$.
}
\label{fig:result4}
\end{figure*}

\section{discussion}

As seen in the previous section, BeHMXBs with relatively short orbital periods and large orbital eccentricities can produce X-ray bursts exceeding the Eddington luminosity during their binary evolution.
The duration of these bursts depends on the binary parameters but typically exceeds 1 Myr.
We compare the observed binary parameters of BeHMXBs with the computational results to investigate how many BeHMXBs can potentially evolve into ULXs.
In Figure 5, BeHMXBs observed in our galaxy and in nearby galaxies are plotted as filled circles in the orbital period versus orbital eccentricity plane (data are taken from \citealp{townsend2011} and \citealp{coe2015}).
Two observed ULXs with Be donors whose orbital eccentricities and orbital periods are known are plotted as triangles.
Given the colour contours shown in Figure 1, the edges of the parameter regime that satisfy the conditions for bursts above the Eddington luminosity are indicated by the two violet lines.
The systems in the region between the two lines can produce super-Eddington bursts for durations longer than 1 Myr.
The parameters of the binary are the same as those in Figure 1.
The vertical line shows the minimum orbital period of BeHMXBs, $P_{\rm{orb}} = 10 \, \rm{d}$, for reference.
The curve extending from the upper left to the lower right depicts the warp condition of the Be disc, as in Figure 1, at the terminal main sequence stage for a donor of $12 M_{\odot}$.
The region bounded by these four curves is thought to indicate the parameter regime in which BeHMXBs can be observed as transient ULXs.

If we restrict the sample to systems with known orbital periods and orbital eccentricities, more than half of the BeHMXBs are located in the region of potential ULXs.
Although there may be a bias in that BeHMXBs are more likely to be observed in systems with shorter orbital periods, this suggests that a substantial number of BeHMXBs will generate bursts that exceed the Eddington luminosity during their evolutionary process.
Among the observed ULXs with Be donors, Swift J0243.6+6124 \citep{doroshenko2018, wilson2018}, which has a short orbital period,
and another ULX, SMC X-3 \citep{tsygankov2017, weng2017}, which has a slightly larger eccentricity, are just outside of the edge of
the ULX region in Figure 5.
As will be shown later, however, they are included in the ULX region when the parameter is changed slightly.

\begin{figure}
\includegraphics[width=\columnwidth]{fig_mod_m12_r511_d8_newfig5.eps}
\caption{
Figure 5: Comparison of observed BeHMXBs with the parameter regimes where super-Eddington bursts can occur.
In the region bounded by the two violet curves, super-Eddington luminosity bursts occur for durations longer than 1 kyr.
The vertical line indicates an orbital period of 10 d, and the orange curve shows the warp conditions at the last stage of the main sequence for reference.
The filled circles indicate the BeHMXBs observed in our galaxy and in the Magellanic Clouds.
The triangles denote ULXs with Be donors (SMC X-3 and swift J0243.6+6124).
The shaded area indicates the region in which Roche lobe overflow starts while the donor is in the main sequence phase.
The binary and donor parameters are the same as those in Figure 1.
}
\label{fig:result5}
\end{figure}

Nevertheless, for a system with a short orbital period or a large orbital eccentricity, the donor could fill its Roche lobe during the evolutionary process.
Before the donor fills the Roche lobe and mass transfer starts, the orbital eccentricity rapidly decreases due to strong tide and the system approaches a circular orbit \citep{hurley2002}.
If the mass transfer is stable, the donor in such a system will lose its circumstellar disc and the system will evolve into a persistent ULX where a large mass transfer rate is achieved via RLOF.
The shaded region in the upper left part of Figure 5 represents the binary parameter regime in which the donor evolves to a larger radius and fills the Roche lobe before it leaves the main sequence.
Figure 5 shows that binary systems that proceed to RLOF during the main-sequence phase are limited to those with short periods or large eccentricities.
Among the observed BeHMXBs, there are few such systems.
Considering that BeHMXBs dominate a large fraction of HMXBs, there may actually be fewer neutron star ULXs accreting by the RLOF.
Rather, neutron stars that show transient bursts following accretion from the Be disc may account for the majority of neutron star ULXs.
Conversely, for systems with orbital periods shorter than 1000 d, RLOF will occur when the donor evolves beyond the main sequence.
Because mass transfer by RLOF beyond the main sequence is expected to be unstable, however, such systems are likely to quickly transition to common envelope evolution \citep{taam2010, ivanova2013}.
In that case, they will radiate brighter at redder wavelengths than at X-ray wavelengths.
The period during which these objects emit brightly is less than a few months.
This period is extremely short compared to the lifetime of the system, so it is rarely observed, although theoretical and observational studies have been carried out in recent years \citep{pastorello2019,blagorodnova2021,matsumoto2022}.

Observational evidence suggests that BeHMXBs with small orbital eccentricities are more likely to produce large bursts \citep{reig2011}.
This may be because there are few systems with large eccentricities, to begin with.
Otherwise, perhaps this is because only systems with small eccentricities can survive as BeHMXBs because systems with large eccentricities quickly enter the RLOF phase.
Because there should be a relationship between the orbital eccentricity and the kick at the time of neutron star formation, the identification of the central object of a ULX may also limit the process of supernova explosions.

As seen in the previous section, the observability as a ULX strongly depends on the parameters of the Be disc.
For example, even a small increase in the density of the Be disc will increase the number of BeHMXB systems that can satisfy the ULX condition at a certain epoch of the evolutionary process.
Figure 6 shows the results when the base density of the Be disc becomes higher: $\rho_{0} = 5 \times 10^{-10} \rm{g \, cm}^{-3}$.
The other parameters are the same as those in Figure 5.
When the disc density increases by one order of magnitude,
the area of possible ULXs becomes wider and two observed ULXs (Swift J0243.6+6124 and SMC X-3), which were not in this range in Figure 5, are included in the ULX region.
Conversely, as the disc density is reduced, fewer and fewer objects fall within the ULX region.
This suggests that relatively high-density circumstellar discs exist in ULXs with Be donors.

The inclination of the Be disc plane to the orbital plane ($\delta$) also affects the results.
That is, when the vectors of the rotational motion of the Be disc material and the orbital motion of the neutron star are close, the relative velocity of the disc material becomes smaller.
Thus, more material can be captured by the neutron star, the mass transfer rate becomes larger and a super-Eddington burst is more likely to occur.
Consequently, a smaller inclination has the same effect as increasing the density of the Be disc, in the sense that it increases the mass transfer rate.

\begin{figure}
\includegraphics[width=\columnwidth]{fig_mod_m12_r510_d8_newfig6.eps}
\caption{
Figure 6: As in Figure 5 but showing the results of increasing the base density of the Be disc by one order of magnitude.
}
\label{fig:result6}
\end{figure}

In the above results, however, the inclination of the disc plane is treated as being constant and independent of the warp condition.
The inclination of the disc plane and the orbital plane could affect the warp conditions.
It is also possible that the rotation axis of the donor and the Be disc plane evolve with time as a result of tidal effects arising from the neutron star.
A dynamic treatment of the warped disc requires numerical simulations and is not considered in depth here (see \citealp{martin2014, franchini2021}).

As well as the inclination of the disc plane, the base density of the disc could evolve as the donor evolves.
Be-type stars are generally thought to spin quickly \citep{rivinius2013}; nevertheless, it is not clear how long the disc is maintained when it spins down as a result of the angular momentum loss associated with the formation of the disc or the tidal forces of the neutron star.
There are many unknowns with respect to the formation of the discs and their densities.
As our understanding of Be-type stars increases, we will be able to learn more about the X-ray emission mechanism of BeHMXBs.

The typical duration for which BeHMXB can produce bursts of super-Eddington luminosity, as calculated here, is approximately 1 Myr.
Nevertheless, the typical lifetime of a BeHMXB is approximately 10 Myr.
Consequently, only approximately 10\% of binary systems in the above parameter regime will actually become ULXs.
Moreover, even if they are in this ULX region, such bursting objects are transient sources, which means that they are not always visible to observers.
Considering the frequency of Type II BeHMXB bursts (less than 10\% of cycles because of the disruption of the Be disc after a large burst \citep{reig2011}), a monitoring period of several years for each object may be necessary to detect a large burst.
Thus, less than approximately 1\% of BeHMXBs are likely to be recognised as active ULXs.

Based on this estimate, there are more than 200 nearby BeHMXBs (of these, 31 objects with known periods and eccentricities are shown in Figure 5).
Conversely, there are three ULX systems with Be donors in the vicinity; SMC X-3, Swift J2043.6+6124 and RX J0209.6-7427 \citep{chandra2020}.
This ratio is not bad compared with the previous estimate.
Because there is an observational bias, further observations are needed.
In principle, it is believed that the observed luminosity is amplified by beaming relative to the intrinsic luminosity \citep{king2017}.
Several recent studies have argued, however, that the effect of beaming is limited \citep{takahashi2017, erkut2020, mushtukov2021}.
The absence of a bias in the viewing angle of the object is convenient for comparisons of the sample with theoretical models.

The potential of a BeHMXB to be a ULX will help in estimates of the number of neutron star binary systems that could be gravitational wave sources and in understanding the ratio of neutron stars to black holes in ULXs \citep{kuranov2021}.
HMXBs containing neutron stars are believed to be the progenitors of double neutron star binary systems.
Particularly, the formation of double neutron star systems, which merge within the age of the universe, requires the passage of moderately spaced HMXBs during their evolution \citep{belczynski2002, tauris2015}.
BeHMXBs, with orbital periods ranging from a few to several hundreds of days, would represent an important fraction of progenitor candidates.
If BeHMXBs can be observed as ULXs, because they can be observed even outside of the Galaxy, their population may provide valuable information for understanding the merger rate of double neutron star binaries.

It is also interesting to determine how many of the ULXs have a black hole origin.
Particularly, the fraction of intermediate-mass black holes is an important question in terms of the formation of massive black holes \citep{kaaret2017}.
If we can classify the ULX sample into those with neutron star origins and those with black hole origins, we will be one step closer to estimating the number of intermediate-mass black holes.
The fact that BeHMXBs, which represent the largest fraction of neutron star binary systems, could be ULXs suggests that the fraction of neutron stars in the ULX sample is not small.
To promote these studies, the model calculations provided in this study can be further developed by adding a method for population synthesis.
Because many results have already been obtained in population syntheses considering ULXs supplied via RLOF \citep{fragos2015, shao2015, wiktorowicz2017, kuranov2020, kuranov2021}, it should be possible to estimate the ULX population more comprehensively by including major mass transport from the Be disc.
These findings are also expected to be useful in extracting information concerning X-ray sources from donor information that will be obtained by next-generation optical telescopes.

\section{Conclusions}

Recently, several ULXs have been discovered that are emitted from neutron stars.
Although the number of observations is still small, some are known to have Be-type donors.
In this paper, we considered the following conditions for BeHMXBs to be ULXs as considered in \citet{karino2021}.
\begin{itemize}
\item The donor does not fill the Roche lobe at the periastron passage.
\item The mass capture rate from the Be disc via the BHL mechanism exceeds the Eddington limit.
\item Even considering the outflow from a slim disc, mass accretion onto the neutron star occurs at a rate exceeding the Eddington limit.
\end{itemize}
Besides this, we considered the following highly indefinite condition.
\begin{itemize}
\item The Be disc warps inside the truncation radius.
\end{itemize}
We hypothesised that BeHMXBs could produce super-Eddington bursts if all these conditions were satisfied, and we searched for the corresponding parameter regime.

The ULX region changes with the evolution of the donor and the binary system.
Particularly, a system with a high-density Be disc aligned with the orbital plane tends to become a ULX.
Additionally, the ULX region depends on the radius of the donor; systems with an expanded donor are likely to produce super-Eddington bursts.
In this study, we also investigated how long a BeHMXB could be a transient ULX with bursts above the Eddington luminosity during the evolution of the donor and binary system.
By investigating the above conditions during the binary evolution, we found that BeHMXBs can typically be ULXs for a duration of approximately 1 Myr during their evolution.
A comparison with nearby BeHMXBs indicates that many binary systems can become ULXs during their evolution.
Because BeHMXBs are common objects and a significant number of them can become ULXs, we conclude that a reasonable fraction of the observed ULXs are BeHMXBs.
Additionally, a significant number of BeHMXBs evolve into stationary ULXs as a result of RLOF caused by the expansion of the donor.
If the ULX population can be deduced from the HMXB population, it can be used to estimate how many non-pulsating ULXs have a neutron star origin.
In the future, it will be necessary to construct a common evolutionary model for Be-type HMXBs and ULXs with Be donors in combination with binary population synthesis computations.

\section*{Data Availability}

The data underlying this article will be shared on reasonable request to the corresponding author.

\section*{Acknowledgements}
The author wishes to thank the referee for the helpful suggestions.
This work was supported by JSPS KAKENHI Grant Number 18K03706.

\bibliographystyle{mnras}


\section{#1}}%
\newcommand*\patchAmsMathEnvironmentForLineno[1]{%
  \expandafter\let\csname old#1\expandafter\endcsname\csname #1\endcsname
  \expandafter\let\csname oldend#1\expandafter\endcsname\csname end#1\endcsname
  \renewenvironment{#1}%
     {\linenomath\csname old#1\endcsname}%
     {\csname oldend#1\endcsname\endlinenomath}
\newcommand*\patchBothAmsMathEnvironmentsForLineno[1]{%
  \patchAmsMathEnvironmentForLineno{#1}%
  \patchAmsMathEnvironmentForLineno{#1*}}%
\AtBeginDocument{%
\patchBothAmsMathEnvironmentsForLineno{equation}%
\patchBothAmsMathEnvironmentsForLineno{align}%
\patchBothAmsMathEnvironmentsForLineno{flalign}%
\patchBothAmsMathEnvironmentsForLineno{alignat}%
\patchBothAmsMathEnvironmentsForLineno{gather}%
\patchBothAmsMathEnvironmentsForLineno{multline}%
}

\catcode`@=11
\long\def\@makecaption#1#2{
   \vskip 10pt
   \setbox\@tempboxa\hbox{{\small\bf #1.} \ {\small #2}}
   \ifdim \wd\@tempboxa >\hsize
   {\small\bf #1.} \ {\small #2}\par
   \else
        \hbox to\hsize{\hfil\box\@tempboxa\hfil}
   \fi}
\catcode`@=12

\catcode`@=11
\def\secteqno{\@addtoreset{equation}{section}%
\def\thesection.\arabic{equation}}{\thesection.\arabic{equation}}}
\def\endsecteqno{\def\thesection.\arabic{equation}}{\@ifundefined{chapter}%
{\arabic{equation}}{\thechapter.\arabic{equation}}}}
\newcounter{subequation}
\def\alph{subequation}{\alph{subequation}}
\def\sneqnarray{\stepcounter{equation}\let\@currentlabel=\thesection.\arabic{equation}}
\setcounter{subequation}{1}
\def\@eqnnum{{\rm (\thesection.\arabic{equation}}\alph{subequation})}}
\global\@eqcnt\z@\tabskip\@centering\let\\=\@eqncr\let\@@eqncr=\@@sneqncr
$$\halign to \displaywidth\bgroup\@eqnsel\hskip\@centering
 $\displaystyle\tabskip\z@{##}$&\global\@eqcnt\@ne
 \hskip 2\arraycolsep \hfil${##}$\hfil
 &\global\@eqcnt\tw@ \hskip 2\arraycolsep
$\displaystyle\tabskip\z@{##}$\hfil
tabskip\@centering&\llap{##}\tabskip\z@\cr}
\def\@@sneqncr\egroup $$\global\@ignoretrue{\@@sneqncr\egroup $$\global\@ignoretrue}
\def\@@sneqncr{\let\@tempa\relax
   \ifcase\@eqcnt \def\@tempa{& & &}\or \def\@tempa{& &}
   \else \def\@tempa{&}\fi
     \@tempa \if@eqnsw\@eqnnum\stepcounter{subequation}\fi
     \global\@eqnswtrue\global\@eqcnt\z@\cr}
\def\def\@lbibitem[##1]##2{\@bibitem{##2}}{\def\@lbibitem[##1]##2{\@bibitem{##2}}}
\catcode`@=12

\def\begin{equation}{\begin{equation}}
\def\end{equation}{\end{equation}}
\def\begin{eqnarray}{\begin{eqnarray}}
\def\end{eqnarray}{\end{eqnarray}}
\def\nonumber{\nonumber}

\def\alpha}  \def\g{\gamma} \def\G{\Gamma{\alpha}  \def\g{\gamma} \def\G{\Gamma}
\def\lambda} \def\lap{\lambda^{\prime}} \def\pa{\partial} \def\de{\delta} \def\De{\Delta} \def\dag{\dagger{\lambda} \def\lap{\lambda^{\prime}} \def\pa{\partial} \def\de{\delta} \def\De{\Delta} \def\dag{\dagger}
\def\epsilon} \def\nb{\bm{\nabla}}  \def\Oc{{\rm O}} \def\S{{\rm S}{\epsilon} \def\nb{\bm{\nabla}}  \def\Oc{{\rm O}} \def\S{{\rm S}}
\def{\bm \nabla}{{\bm \nabla}}
\def{\bm \sigma}{{\bm \sigma}}

\documentclass[preprint, aps, superscriptaddress, prd, showpacs,
amsmath, amssymb, nofootinbib, showkeys]{revtex4-1}

\usepackage[german,english]{babel}
\usepackage{hyperref}
\usepackage{bm
\usepackage{bbm}
\usepackage{braket}
\usepackage{slashed}
\usepackage{multirow}
\usepackage{epsfig}
\usepackage{epstopdf}
\usepackage{diagbox}
\usepackage[mathlines]{lineno
\usepackage[utf8]{inputenc}
\usepackage{subfigure}
\usepackage{diagbox}
\usepackage {xcolor}

\allowdisplaybreaks

\DeclareMathOperator\arctanh{arctanh}

\begin{document}

\title{Unveiling Nucleon 3D Chiral-Odd Structure with Jet Axes}
\author{Wai Kin Lai}
\email{wklai@m.scnu.edu.cn}
\affiliation{Guangdong Provincial Key Laboratory of Nuclear Science, Institute of Quantum Matter, South China Normal University, Guangzhou 510006, China}
\affiliation{Guangdong-Hong Kong Joint Laboratory of Quantum Matter,
Southern Nuclear Science Computing Center, South China Normal University, Guangzhou 510006, China}
\affiliation{Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA}
\author{Xiaohui Liu}
\email{xiliu@bnu.edu.cn}
\affiliation{Center of Advanced Quantum Studies, Department of Physics,
Beijing Normal University, Beijing 100875, China}
\affiliation{Center for High Energy Physics, Peking University, Beijing 100871, China}
\author{Manman Wang}
\affiliation{Center of Advanced Quantum Studies, Department of Physics,
Beijing Normal University, Beijing 100875, China}
\author{Hongxi Xing}
\email{hxing@m.scnu.edu.cn}
\affiliation{Guangdong Provincial Key Laboratory of Nuclear Science, Institute of Quantum Matter, South China Normal University, Guangzhou 510006, China}
\affiliation{Guangdong-Hong Kong Joint Laboratory of Quantum Matter,
Southern Nuclear Science Computing Center, South China Normal University, Guangzhou 510006, China}

\date{\today}

\begin{abstract}
We reinterpret jet clustering
as an axis-finding procedure which, along with the proton beam, defines the virtual-photon transverse momentum $q_T$ in deep inelastic scattering (DIS). In this way, we are able to probe the nucleon intrinsic structure using jet axes in a fully inclusive manner, similar to the Drell-Yan process. We present the complete list of azimuthal asymmetries and the associated factorization formulae at leading power for deep-inelastic scattering of a nucleon. The factorization formulae involve both the conventional time-reversal-even (T-even) jet function and the T-odd one, which have access to all transverse-momentum-dependent parton distribution functions (TMD PDFs) at leading twist.
Since the factorization holds as long as $q_T \ll Q$, where $Q$ is the photon virtuality, the jet-axis probe into the nucleon structure should be feasible for machines with relatively low energies such as the Electron-Ion Collider in China (EicC). We show that, within the winner-take-all (WTA) axis-finding scheme, the coupling between the T-odd jet function and the quark transversity or the Boer-Mulders function could induce sizable azimuthal asymmetries at the EicC, the EIC and HERA.
We also give predictions for the azimuthal asymmetry of back-to-back dijet production in $e^+e^-$ annihilation at Belle and other energies.

\end{abstract}

\maketitle

\section{Introduction}
Recently jets and jet substructure have been proposed as alternative probes for
portraying the full three-dimensional (3D) image of a nucleon and enriched the content of the transverse-momentum-dependent (TMD) spin physics~\cite{Kang:2017glf,Liu:2018trl,Gutierrez-Reyes:2019msa,Arratia:2020nxw,Liu:2020dct,Gutierrez-Reyes:2018qez,Gutierrez-Reyes:2019vbx,Kang:2022dpx,Kang:2020xyq}. The jet probe into the nucleon structure has been shown to be able to access the TMD parton distribution functions, including the Sivers function of a transversely polarized nucleon. Conventionally, we require the jets to acquire large transverse momenta, and therefore jets are regarded only feasible for high-energy colliders such as the LHC but practically challenging for machines with a relatively low center-of-mass energy such as the Electron-Ion Collider in China (EicC)~\cite{Anderle:2021wcy} or detectors more optimized for low energy scales, such as the EIC Comprehensive Chromodynamics Experiment (ECCE)~\cite{AbdulKhalek:2021gbh}. However, in this work, we will argue that this is not the case by reinterpreting jet clustering as an axis-finding procedure to measure the virtual photon $q_T$, which allows an inclusive probe of the TMD spin physics suitable also for low energy machines~\cite{liutalk}.

In order to maximize the full outreach of the jet probe into the complete list of the nucleon spin structure, the concept of the time-reversal-odd (T-odd) jet function was proposed recently~\cite{Liu:2021ewb}. The T-odd jet function couples directly to the chiral-odd nucleon parton distributions, such as the quark transversity and the Boer-Mulders function of the proton. It immediately opens up many unique opportunities for probing the nucleon intrinsic spin dynamics using jets, which were thought to be impossible.
Besides, the T-odd jet function is interesting by its own, since it could ``film" the QCD non-perturbative dynamics by continuously changing the jet axis from one to another.

In this work, we study the phenomenology of the T-odd jet function in deep-inelastic scattering (DIS) of a nucleon and $e^+e^-$ annihilation. In Section~\ref{sec:measure_qt}, we explain how the jet axis is used for measuring the photon $q_T$ in a fully inclusive way and argue why the jet-axis probe of the nucleon spin and the TMDs is feasible even for low-energy machines such as the Electron-Ion Collider in China (EicC) and Belle~\cite{Accardi:2022oog}.
In Section~\ref{sec:T_odd_jet}, we briefly review the notion of the T-odd jet function. In Section~\ref{sec:DIS}, we give the complete list of the azimuthal asymmetries in the jet-axis probe in deep-inelastic scattering of a nucleon. We give predictions on the azimuthal asymmetries associated with the couplings of the T-odd jet function with the quark transversity and the Boer-Mulders function at the EicC, the EIC~\cite{AbdulKhalek:2021gbh}, and HERA.
In Section~\ref{sec:ee}, we study the azimuthal asymmetry of back-to-back dijet production in $e^+e^-$ annihilation, which is induced by the T-odd jet function. In Section~\ref{sec:summary}, we give a summary and an outlook.

\section{Measuring photon $q_T$ in DIS}\label{sec:measure_qt}

All conventional probes of the TMDs and the spin structure are more or less equivalent to measuring the virtual-photon transverse momentum $q_T$ with respect to two pre-defined axes. For instance, in the Drell-Yan process, the incoming nucleon beams naturally set up the $\pm z$-axis and the photon transverse momentum is then straightforwardly determined. In DIS, since we only have one nucleon beam, we thus need another direction to define the photon $q_T$. Tagging a final-state hadron becomes a natural option for this purpose, and in this case, the photon $q_T$ is then measured with respect to the nucleon beam and the tagged hadron momentum $P_h$. This is nothing but the semi-inclusive deep inelastic scattering (SIDIS).

Finding an axis for measuring the photon $q_T$ in DIS is certainly not limited to tagging hadrons. Many other strategies could also help here, such as the final-state-particle clustering. The procedure follows exactly the jet clustering algorithms, but a with different emphasis. Here the jet clustering procedure is barely a recursive algorithm for us to determine the axes, which once being determined, we measure the photon $q_T$ with respect to one of them and the proton beam to probe the nucleon structure, while totally forget about the jet, as illustrated in Fig.~\ref{fg:qt-jet}. Therefore, the jet-axis probe is fully differential just like the SIDIS.

 \begin{figure}[htbp]
  \begin{center}
   \includegraphics[scale=0.35]{jet-axis-define.pdf}
\caption{Photon $q_T$ by the jet axis in DIS, where one of the jet axes and the proton beam determine the ${\bar n}$ and $n$ directions respectively, which decide the photon $q_T$. }
  \label{fg:qt-jet}
 \end{center}
\end{figure}

Based on what we have described, we can derive the factorization formula for the jet-axis probe. Formally, the factorization theorem reads
\begin{eqnarray}
d \sigma \propto f_i^{U/T}(x,p_T^2) \otimes {\cal J}_{1,i}(z,k_T^2)
+ g_i^{U/T}(x,p_T^2)\otimes{\cal J}_{T,i}(z,k_T^2)  \,, \label{eq:fact_1}
\end{eqnarray}
where
$f_i^{U/T}(x,p_T^2) $ and $g_i^{U/T}(x,p_T^2)$ are the proton partonic transverse-momentum distributions for parton flavor $i$, and
${\cal J}_{1,i}(z,k_T^2)$ and ${\cal J}_{T,i}(z,k_T^2)$ are both the jet-axis-finding functions (jet functions) which encode the perturbatively calculable jet clustering procedure.
The conventional jet function ${\cal J}_{1,i}(z, k_T^2)$ is induced by an unpolarized quark, while a transversely polarized quark  gives rise to the time-reversal odd (T-odd) jet function ${\cal J}_{T,i}(z,k_T^2)$. The detailed factorization form and the definition of the T-odd jet function, ${\cal J}_{T,i}(z,k_T^2)$, will be given in the following sections.
 Here we note that the factorization theorem holds as long as $Q \gg k_T $, which shares the same requirement for the SIDIS factorization to be valid. In this sense, just like the SIDIS,  the jet-axis probe will also be low-energy-machine friendly, and could likely be implemented at the EicC.

  To adapt the jet axis finding procedure to low-energy machines, instead of using the usual $k_T$-type jet algorithms that are widely used at the LHC, in DIS we default to the energy-type jet algorithms which is more feasible for clustering particles with low transverse momenta and populated in the forward/backward rapidities. For instance, we can adopt the spherically-invariant jet algorithm~\cite{Cacciari:2011ma}, defined by
\begin{align}
d_{ij}={\rm min}(E^{-2}_i,E^{-2}_j)\frac{1-\cos\theta_{ij}}{1-\cos R}\,,\quad d_{iB}=E_i^{-2}\,.\label{eq:spherical_jet_alg}
\end{align}
where $\theta_{ij}$ is the angle between particles $i$ and $j$, while $E_{i}$ and $E_j$ are the energy carried by them. For TMD studies, the radius parameter $R$ will be chosen such that $R \sim {\cal O}(1) \gg q_T/Q$.

\section{T-odd jet function} \label{sec:T_odd_jet}
The inclusive photon $q_T$ cross sections with respect to the proton beam and the jet axis can be written in terms of a factorization theorem Eq.~(\ref{eq:fact_1}) derived from the soft-collinear effective theory
(SCET)~\cite{Bauer:2000yr, Bauer:2001ct, Bauer:2001yt, Bauer:2002nz}. The factorization theorem involves the
transverse-momentum-dependent (TMD) correlator
\begin{align}
{\cal J}^{ij}(z,k_T)&=\frac{1}{2z}\sum_X\int \frac{dy^+d^2\bm{y}_T}{(2\pi)^3}e^{ik\cdot y}
\langle 0|{\chi}^i_{\bar{n}}(y)|JX\rangle\langle JX|\bar{\chi}^j_{\bar{n}}(0)|0\rangle |_{y^-=0}\,, \label{eq:jet_corr}
\end{align}
where $\bar{n}$ is a light-like vector along the direction of the jet, $\chi_{\bar{n}}=W^\dagger_{\bar{n}}\xi_{\bar{n}} $ is the product of the collinear quark field $\xi_{\bar{n}}$ and the collinear Wilson line  $W^\dagger_{\bar{n}}$. Here, $z$ is the momentum fraction of the jet with respect to the fragmenting parton which initiates the jet, i.e. $z=P_J^-/ k^-$, with $P_J$ being the jet momentum that defines the jet axis, and $k$ the momentum of the fragmenting quark.
 The jet algorithm dependence is implicit in Eq.~(\ref{eq:jet_corr}), which determines the $P_J$ and hence the jet axis, and can be calculated perturbatively.

Conventionally, only the chiral-even Dirac structure $\slashed{\bar{n}}$ in Eq.~(\ref{eq:jet_corr}) was considered. However, as noted in Ref.~\cite{Liu:2021ewb}, in the nonperturbative regime in which $k_T \sim \Lambda_{\rm QCD}$,
spontaneous chiral symmetry breaking
leads to a nonzero component of the jet which is both time-reversal-odd (T-odd) and chiral-odd, when the jet axis is different from the direction of the
fragmenting parton.
Therefore,
the correlator in Eq.~(\ref{eq:jet_corr}) in general is a sum of two structures:
\begin{align}
{\cal J}(z,k_T)&= {\cal J}_1(z,k_T^2)\frac{\slashed{\bar{n}}}{2}
+i{\cal J}_T(z,k_T^2)\frac{\slashed{k}_T\slashed{\bar{n}}}{2}\,,     \label{eq:jet_expd}
\end{align}
where ${\cal J}_1(z,k_T^2)$ is the traditional jet function, and ${\cal J}_T(z,k_T^2)$ is the T-odd jet function. Due to its chiral-odd nature, an immediate application of the T-odd jet function is to probe the chiral-odd TMD PDFs of the nucleons in DIS, such as the Boer-Mulder function and the transversity, which were thought to be impossible to access using jets.

The T-odd jet function has the following advantages:
\begin{itemize}
\item
{\textit{Universality}}\\
 Like the traditional jet function, the T-odd jet function is process independent.
\item
{\textit{Flexibility}}\\
The flexibility of choosing a jet recombination scheme and hence the jet axis allows us to adjust sensitivity of the jet function to different nonperturbative contributions. This provides an opportunity to “film” the QCD nonperturbative dynamics, if one continuously changes the axis from one to another.
\item
{\textit{Perturbative predictability}} \\
Since a jet contains many hadrons, the jet function has more perturbatively calculable degrees of freedom than the fragmentation function. For instance, in the winner-take-all (WTA) scheme, for $R\sim \mathcal{O}(1)\gg |\bm{q}_T|/E_J$, the $z$-dependence in the jet function is completely determined~\cite{Gutierrez-Reyes:2019vbx}:
\begin{align}
{\cal J}(z,k_T,R)=\delta(1-z) {\mathfrak J}(k_T)+ {\cal O}\left( \frac{k_T^2}{E_J^2R^2} \right)\,. \label{eq:WTA_J}
\end{align}
\item
{\textit{Nonperturbative predictability}}\\
Similar to the study in Ref.~\cite{Becher:2013iya}, the T-odd jet function can be factorized into a product of a perturbative coefficient
 and a nonperturbative factor. The nonperturbative factor has an operator definition~\cite{Vladimirov:2020umg},
and as a vacuum matrix element, it can be calculated on the lattice~\cite{Shanahan:2020zxr, Zhang:2020dbb}. This is unlike the TMD fragmentation function, which is an operator element with a final-state hadron tagged, making evaluation on the lattice impossible by known techniques.
\end{itemize}

The T-odd jet function will show up in various jet observables which are sensitive to nonperturbative physics. In the following,
we study the azimuthal asymmetries in the jet-axis probe in DIS and back-to-back dijet production in $e^+e^-$ annihilation.

\section{Photon $q_T$ with respect to the jet axis in deep-inelastic scattering} \label{sec:DIS}

Consider deep-inelastic scattering of an electron off a polarized nucleon $e^-(l)+N(P)\to  e^-(l')+J(P_J)+X$, ($N=p,n$), in which we tag a jet and specify the jet axis with some recombination scheme.
We define the ${\bm q}_T$ of the virtual photon by going to the so-called factorization frame,
in which the proton beam direction and the jet axis direction are exactly opposite to each other, as shown in Fig.~\ref{fg:frames} (a). Alternatively, one can go to the gamma-nucleon system (GNS), a frame in which the virtual photon momentum and the proton beam are head-to-head (including the case of proton being at rest), and define ${\bm P}_{J\perp}$ of the jet as in Fig.~\ref{fg:frames} (b). One can show that
${\bm q}_T=-{\bm P}_{J\perp}/z$ up to corrections of order $1/Q^2$. Therefore, measuring ${\bm P}_{J\perp}$ is
equivalent to measuring ${\bm q}_T$. In the following, we will describe the kinematics in the GNS system, which is a convention commonly used in SIDIS~\cite{Bacchetta:2006tn}.

 \begin{figure}[htbp]
  \begin{center}
 \subfigure[]{
   \includegraphics[scale=0.6]{frame_factor.pdf}
     }
   \subfigure[]{
   \includegraphics[scale=0.6]{frame_GNS.pdf}
  }
\caption{Axes in DIS in different frames: (a) the factorization frame, in which ${\bm q}_T$ is defined, and (b) the GNS system, in which ${\bm P}_{J\perp}$ is defined.}
  \label{fg:frames}
  \end{center}
\end{figure}

 \begin{figure}[htbp]
  \begin{center}
   \includegraphics[scale=0.3]{ep.pdf}
\caption{Kinematic configuration of DIS in the GNS system.}
  \label{fg:ep}
 \end{center}
\end{figure}

Let $M$ be the mass of the nucleon $N$ and $q=l-l'$ is the momentum carried by the virtual photon with virtuality $Q^2=-q^2$. We introduce the invariant variables
\begin{align}
x=\frac{Q^2}{2P\cdot q}\,,\quad
y=\frac{P\cdot q}{P\cdot l}\,,\quad
z=\frac{P\cdot P_J}{P\cdot q}\,,\quad
\gamma=\frac{2Mx}{Q}\,.
\end{align}
In the nucleon rest frame,
we can define the perpendicular component of any 3-vector as the component perpendicular to the virtual photon momentum,  ${\bm q}$. Equivalently, in Lorentz invariant notations, given any 4-vector $v^\mu$, we
define its perpendicular component by $v^\mu_\perp=g_\perp^{\mu\nu}v_\nu$, where
\begin{align}
g_\perp^{\mu\nu}&=g^{\mu\nu}-\frac{q^\mu P^\nu+P^\mu q^\nu}{P\cdot q(1+\gamma^2)}+\frac{\gamma^2}{1+\gamma^2}
\left(\frac{q^\mu q^\nu}{Q^2}-\frac{P^\mu P^\nu}{M^2}\right)\,.
\end{align}
The 3-momenta ${\bm l}$ and ${\bm l}'$ define a plane, with respect to which we can define the azimuthal angle of any 3-vector perpendicular to ${\bm q}$. In Lorentz invariant notations, this is equivalent to defining the azimuthal angle of any 4-vector $v^\mu$ by
\begin{align}
\cos \phi_v=-\frac{l_\mu v_\nu g_\perp^{\mu\nu}}{\sqrt{l_\perp^2 v^2_\perp}}\,,\quad
\sin\phi_v=-\frac{l_\mu v_\nu \epsilon^{\mu\nu}}{\sqrt{l_\perp^2 v^2_\perp}}\,,
\end{align}
where
\begin{align}
\epsilon_\perp^{\mu\nu}&=\epsilon^{\mu\nu\rho\sigma}\frac{P_\rho q_\sigma}{P\cdot q\sqrt{1+\gamma^2}}\,,
\end{align}
with $\epsilon^{0123}=1$. We denote the azimthal angles of the jet momentum $P_J$ and the nucleon spin $S$ by $\phi_J$ and $\phi_S$ respectively.
The definitions of $P_{J\perp}$, $\phi_J$, and $\phi_S$ are depicted pictorially in Fig.~\ref{fg:ep}. The nucleon spin is decomposed as sum of a longitudinal and a perpendicular component,
\begin{align}
S^\mu=S_\parallel\frac{P^\mu-\frac{M^2}{P\cdot q}q^\mu }{M\sqrt{1+\gamma^2}}+S^\mu_\perp\,,
\quad   S_\parallel=\frac{S\cdot q}{P\cdot q}\frac{M}{\sqrt{1+\gamma^2}}\,.
\end{align}
The helicity of the incoming electron is denoted by $\lambda_e$.
We define $\psi$ as the azimuthal angle of $\bm l'$ around $\bm l$.
The fully differential cross section has the most general form given by
\begin{align}
&\frac{d\sigma}{dxdyd\psi dz d\phi_J dP_{J\perp}^2}=
\frac{\alpha}{xyQ^2}\frac{y^2}{2(1-\epsilon)}\left(1+\frac{\gamma^2}{2x}\right)\left\{ F_{UU,T}+\epsilon F_{UU,L}
+\sqrt{2\epsilon(1+\epsilon)}\cos\phi_J F_{UU}^{\cos\phi_J}\right.\nonumber\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad
+\epsilon\cos(2\phi_J)F^{\cos 2\phi_J}_{UU}+\lambda_e\sqrt{2\epsilon(1-\epsilon)}\sin\phi_J F^{\sin\phi_J}_{LU}\nonumber\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad
+S_\parallel \left[\sqrt{2\epsilon(1+\epsilon)}\sin\phi_J F_{UL}^{\sin\phi_J}+\epsilon\sin(2\phi_J)F_{UL}^{\sin 2\phi_J}\right]\nonumber\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad
+S_\parallel \lambda_e\left[\sqrt{1-\epsilon^2}F_{LL}+\sqrt{2\epsilon(1-\epsilon)}\cos\phi_J F_{LL}^{\cos\phi_J}\right]\nonumber\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad +|{\bm S}_{\perp}|\left[\sin(\phi_J-\phi_S)\left(F_{UT,T}^{\sin(\phi_J-\phi_S)}+\epsilon F_{UT,L}^{\sin(\phi_J-\phi_S)}\right)\right.\nonumber\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad+\epsilon\sin(\phi_J+\phi_S)F_{UT}^{\sin(\phi_J+\phi_S)}+\epsilon\sin(3\phi_J-\phi_S)F_{UT}^{\sin(3\phi_J-\phi_S)}\nonumber\\
&\left.\quad\quad\quad\quad\quad\quad\quad\quad\quad +\sqrt{2\epsilon(1+\epsilon)}\sin\phi_S F_{UT}^{\sin\phi_S}+\sqrt{2\epsilon(1+\epsilon)}\sin(2\phi_J-\phi_S)F_{UT}^{\sin(2\phi_J-\phi_S)}\right]\nonumber\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad +|{\bm S}_{\perp}|\lambda_e\left[\sqrt{1-\epsilon^2}\cos(\phi_J-\phi_S)F_{LT}^{\cos(\phi_J-\phi_S)}
+\sqrt{2\epsilon(1-\epsilon)}\cos\phi_S F_{LT}^{\cos\phi_S}\right.\nonumber\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad \left.\left.+\sqrt{2\epsilon(1-\epsilon)}\cos(2\phi_J-\phi_S)F_{LT}^{\cos(2\phi_J-\phi_S)}\right]\right\}\,,\label{eq:diff_cross_section}
\end{align}
where $\alpha$ is the fine structure constant and
\begin{align}
\epsilon=\frac{1-y-\frac{1}{4}\gamma^2 y^2}{1-y+\frac{1}{2}y^2+\frac{1}{4}\gamma^2y^2}\,.
\end{align}
The structure functions $F$'s on the right-hand side of Eq.~(\ref{eq:diff_cross_section}) depend
on $x, Q^2, z$, and $P_{J\perp}^2$. The subscripts $A,B,C$ of $F_{AB,C}$ denote the polarizations of the incoming lepton, the incoming nucleon, and the virtual photon respectively, with $U$ standing for unpolarized, $T$ for transversely polarized, and $L$ for longitudinally polarized.
Up to corrections of $\mathcal{O}(M^2/Q^2)$, we have $\psi\approx \phi_S$, $\gamma^2\approx 0$, and the following approximations for the
coefficients of the $F$'s in Eq.~(\ref{eq:diff_cross_section}):
\begin{align}
&\frac{y^2}{2(1-\epsilon)}\approx 1-y+\frac{1}{2}y^2,\nonumber\\
&\frac{y^2}{2(1-\epsilon)}\epsilon \approx 1-y,\nonumber\\
&\frac{y^2}{2(1-\epsilon)}\sqrt{2\epsilon(1+\epsilon)} \approx (2-y)\sqrt{1-y},\nonumber\\
& \frac{y^2}{2(1-\epsilon)}\sqrt{2\epsilon(1-\epsilon)} \approx y\sqrt{1-y}\,,\nonumber\\
& \frac{y^2}{2(1-\epsilon)}\sqrt{1-\epsilon^2} \approx y\left(1-\frac{1}{2}y\right)\,.
\end{align}
The structure functions $F$'s are convolutions of the nucleon TMD PDFs and the jet functions. As noted in Eq.~(\ref{eq:jet_expd}), there are two
jet functions $\mathcal{J}_1(z,k_T^2)$ and $\mathcal{J}_T(z,k_T^2)$ at leading power in $\Lambda_{QCD}/Q$. At leading order in $M/Q$, there are eight
quark TMD PDFs for the nucleon~\cite{Angeles-Martinez:2015sea}, each encoding a specific correlation between the quark spin and the proton spin as depicted in Table~\ref{tb:TMD_PDF}. Among these eight TMD PDFs, $f_1$, $g_{1L}$, $f^\perp_{1T}$, and $g_{1T}$ are chiral-even, while
$h_1^\perp$, $h_{1L}^\perp$, $h_{1T}$, and $h_{1T}^\perp$ are chiral-odd.
For any functions
$w({\bm p}_T,{\bm k}_T)$, $f(x,p_T^2)$, and $\eta(z,k_T^2)$, we define
\begin{align}
{\cal C}[w f \eta]\equiv x\sum_a e_a^2\int d^2\bm{p}_T\int d^2\bm{k}_T \,
\delta^{(2)}\left(\bm{p}_T-\bm{k}_T-\bm{P}_{J\perp}/z\right)
w(\bm{p}_T,\bm{k}_T)f^a(x,p_T^2)\eta^a(z,k_T^2)\,,
\end{align}
where $a$ denotes a quark or antiquark flavor.
At leading order in $\alpha_s$ and $M/Q$, the nonvanishing $F$'s are given by
\begin{align}
& F_{UU,T}={\cal C}[f_1{\cal J}_1]\,,\label{eq:F_UU,T}\\
& F_{LL}={\cal C}[g_{1L}{\cal J}_1]\,,\\
& F^{\sin(\phi_J-\phi_S)}_{UT,T} ={\cal C}\left[-\frac{\hat{\bm{h}}\cdot \bm{p}_T}{M}
f^\perp_{1T}{\cal J}_1\right]\,,\\
& F^{\cos(\phi_J-\phi_S)}_{LT}={\cal C}\left[\frac{\hat{\bm{h}}\cdot \bm{p}_T}{M}
g_{1T}{\cal J}_1\right]\,,\\
& F_{UU}^{\cos(2\phi_J)} ={\cal C}\left[-\frac{(2(\hat{\bm{h}}\cdot \bm{k}_T)(\hat{\bm{h}}\cdot \bm{p}_T)-\bm{k}_T\cdot\bm{p}_T)}{M}
h_1^\perp{\cal J}_T\right]\,,\\
& F_{UL}^{\sin(2\phi_J)} ={\cal C}\left[-\frac{(2(\hat{\bm{h}}\cdot \bm{k}_T)(\hat{\bm{h}}\cdot \bm{p}_T)-\bm{k}_T\cdot\bm{p}_T)}{M}
h_{1L}^\perp{\cal J}_T\right]\,,\\
& F^{\sin(\phi_J+\phi_S)}_{UT} ={\cal C}\left[-\hat{\bm{h}}\cdot \bm{k}_T
h_1{\cal J}_T\right]\,,\\
& F_{UT}^{\sin(3\phi_J-\phi_S)} ={\cal C}\left[\frac{2(\hat{\bm{h}}\cdot \bm{p}_T)(\bm{p}_T\cdot \bm{k}_T)+
\bm{p}_T^2(\hat{\bm{h}}\cdot \bm{k}_T)-4(\hat{\bm{h}}\cdot \bm{p}_T)^2(\hat{\bm{h}}\cdot \bm{k}_T)}{2M^2}
h_{1T}^\perp{\cal J}_T\right]\,,\label{eq:F_UT_sin3}
\end{align}
where $\hat{\bm{h}}=\bm{P}_{J\perp}/|\bm{P}_{J\perp}|$ and $h_1= h_{1T}+\frac{\bm{p}_T^2}{2M^2}h_{1T}^\perp$\,.
From Eqs.~(\ref{eq:F_UU,T})-(\ref{eq:F_UT_sin3}), we see that the T-even jet function $\mathcal{J}_1$ couples to the chiral-even
TMD PDFs, while the T-odd jet function $\mathcal{J}_T$ couples to the chiral-odd TMD PDFs. With both the T-even and T-odd jet functions, one can thus access all eight TMD PDFs at leading twist.

\begin{table}[htbp]
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
\backslashbox{hadron}{quark} & unpolarized & chiral & transverse \\
\hline
$U$& $f_1$ &  & $h_1^\perp$ \\
\hline
$L$&  & $g_{1L}$ & $h^\perp_{1L}$ \\
\hline
$T$& $f^\perp_{1T}$ & $g_{1T}$ & $h_{1T},h^{\perp}_{1T}$\\
\hline
\end{tabular}
\caption{The eight TMD PDFs of a nucleon at leading twist.}
\end{center}
\label{tb:TMD_PDF}
\end{table}

With the known proton TMD PDFs and partial knowledge on the jet functions,
we can make preliminary predictions on the azimuthal asymmetries associated with the jet-axis probe in DIS machines such as the EIC, the EicC, and HERA.
As notes in Section~\ref{sec:T_odd_jet}, the $z$-dependence of the jet functions become trivial with the WTA jet-axis definition. In the following,
we will adopt the spherically-invariant jet algorithm Eq.~(\ref{eq:spherical_jet_alg}) with $R=1$ and the WTA scheme, so that
\begin{align}
{\cal J}(z,k_T,R)&=\delta(1-z) {\mathfrak J}(k_T)+ {\cal O}\left( \frac{k_T^2}{E_J^2R^2} \right)\,, \label{eq:WTA_EJ}
\end{align}
where
\begin{align}
{\mathfrak J}(k_T)&={J}(k_T^2)\frac{\slashed{\bar{n}}}{2}+i{J}_T(k_T^2)\frac{\slashed{k}_T\slashed{\bar{n}}}{2}\,.
\end{align}
We will study the azimuthal asymmetries associated with the terms $|{\bm S}_\perp|\epsilon\sin(\phi_J+\phi_S)F_{UT}^{\sin(\phi_J+\phi_S)}$
and $\epsilon\cos(2\phi_J)F_{UU}^{\cos 2\phi_J}$ in Eq.~(\ref{eq:diff_cross_section}). These terms probe the transversity $h_1$ and the Boer-Mulders function $h_1^\perp$. These terms
can be singled out by specific modulated cross sections. We define a $|{\bm P}_{J\perp}|$-distribution of asymmetry that probes the transversity by
\begin{align}
A^{\sin(\phi_J+\phi_S)}(|{\bm P}_{J\perp}|)&=\frac{2}{|\bm{S}_\perp|\int d\sigma \epsilon}\int d\sigma\sin(\phi_J+\phi_S)\label{eq:A_transversity_1} \\
&=\frac{\langle \epsilon F_{UT}^{\sin(\phi_J+\phi_S)}\rangle}{\bar{\epsilon} \langle F_{UU,T}\rangle}\,,\label{eq:A_transversity_2}
\end{align}
where by $\langle X \rangle$ we mean
\begin{align}
\langle X \rangle&=\int dx \int dy \int d\phi_S \int dz \int d\phi_J\, \frac{\alpha}{xyQ^2}\frac{y^2}{2(1-\epsilon)}X\,,
\end{align}
and
\begin{align}
\bar{\epsilon}=\frac{\int d\sigma \epsilon}{\int d\sigma}\,.
\end{align}
We can write $F_{UU,T}$ and $F_{UT}^{\sin(\phi_J+\phi_S)}$ as
\begin{align}
F_{UU,T}&=x\sum_a e_a^2\int \frac{d^2b}{(2\pi)^2}e^{-i{\bm P}_{J\perp}\cdot{\bm b}}\tilde{f}^a_1(x,b^2)\tilde{J}^a(b^2)\,,\label{eq:F_UU}\\
F_{UT}^{\sin(\phi_J+\phi_S)}&=-ix\sum_a e_a^2\int \frac{d^2b}{(2\pi)^2}\frac{P_{J\perp}^i}{|\bm{P}_{J\perp}|}
e^{-i{\bm P}_{J\perp}\cdot{\bm b}}\tilde{h}^a_1(x,b^2)\partial_{b^i}\tilde{J}^a_T(b^2)\,,\label{eq:F_UT}
\end{align}
where we have used the Fourier transforms $\tilde{f}_1(x,b^2)=\int d^2p_T\, e^{-i{\bm p}_T\cdot{\bm b}}f_1(x,p_T^2)$ and
$\tilde{J}(b^2)=\int d^2k_T\, e^{-i{\bm k}_T\cdot{\bm b}}J(k_T^2)$, and similarly for $\tilde{h}_1$ and $\tilde{J}_T$.
Similar to the treatment in Ref.~\cite{Kang:2014zza}, we include the effect of evolution by including a
Sudakov factor in $b$-space
\begin{align}
\tilde{f}_1(x,b^2,Q)=e^{-S_{\rm pert}-S_{\rm NP}}\tilde{f}_1(x,b^2,Q_0)\,,\\
\tilde{J}(b^2,Q)=e^{-S_{\rm pert}-S_{\rm NP}}\tilde{J}(b^2,Q_0)\,,
\end{align}
and similarly for $\tilde{h}_1$ and $\tilde{J}_T$, where
\begin{align}
S_{\rm pert}&=\int^Q_{\mu_b} \frac{d\mu}{\mu}\,\frac{\alpha_s(\mu)}{2\pi}C_F
\ln \frac{Q^2}{\mu^2}\,,\\
S_{\rm NP}&=\frac{g_2}{2}\ln\left(\frac{b}{b_*}\right)\ln\left(\frac{Q}{Q_0}\right)\,,
\end{align}
with $\mu_b=1.22/b_*$, $b_*=b/\sqrt{1+b^2/b_*^2}$, $Q_0=1.549$~GeV, and $g_2=0.84$. Here,
we have adopted the expression of $S_{\rm pert}$ at leading logarithm.

For the transversity, we use the fitted parametrized form Ref.~\cite{Martin:2014wua}. For the jet functions,
although not mandatory, we will apply the jet charge measurement in order to enhance flavor separation. This amounts to
replacing the overall normalizations of the jet functions by the charge bins $r_a$, whose values for ${J}^a$
and for jet charge $Q_J>0.25$ and $Q_J<-0.25$ have been obtained in Ref.~\cite{Kang:2020fka} and will be used in this work.
For the charge bins associated with $J^a_T$, we will take them as a product $N_a r_a$, where $N_a$ is the ratio
of the overall normalization of the pion Collins function $H_1^{\perp a}$ to the overall normalization of fragmentation function $D_1^a$ as
obtained in Ref.~\cite{deFlorian:2007aj}.
For the $p_T^2$-dependence of $\tilde{J}$ and $\tilde{J}_T$, we use that of the
fragmentation function and the Collins function of the pion as obtained in Ref.~\cite{deFlorian:2007aj}. Figure~\ref{fg:asym_transversity_EIC} (a)
shows the predictions of the asymmetry
$A^{\sin(\phi_J+\phi_S)}(|{\bm P}_{J\perp}|)$ at the EIC according to Eq.~(\ref{eq:A_transversity_2}) (solid lines) and
Eq.~(\ref{eq:A_transversity_1}) (data points)
from simulations using \textsc{Pythia} 8.2~\cite{Sjostrand:2014zea} with the package \textsc{StringSpinner}~\cite{Kerbizi:2021pzn}, which incorporates spin interactions in the event generator.
From Fig.~\ref{fg:asym_transversity_EIC} (a), we see that the theoretical predictions on the $A^{\sin(\phi_J+\phi_S)}(|{\bm P}_{J\perp}|)$ distribution from the factorization formula Eq.~(\ref{eq:diff_cross_section}) roughly agree with the event generator simulations.
In Figure~\ref{fg:asym_transversity_EIC} (b), we show the prediction of $A^{\sin(\phi_J+\phi_S)}(|{\bm P}_{J\perp}|)$ with the E-scheme for the jet-axis definition from \textsc{Pythia} 8.2~\cite{Sjostrand:2014zea} with \textsc{StringSpinner}, with the same kinematic setting as
Fig.~\ref{fg:asym_transversity_EIC} (a). We see that the asymmetry no longer exists in the E-scheme. This is because the asymmetry is
nonvanishing only when the direction of the fragmenting parton which initiates the jet differs with that of the jet axis, which
hardly the case in the E-scheme. In this sense, by choosing different jet axes we are able to ``film" the nonperturbative dynamics of QCD.

 \begin{figure}[htbp]
  \begin{center}
 \subfigure[]{
   \includegraphics[scale=0.3]{Asym_transversity_EIC.pdf}
     }
   \subfigure[]{
   \includegraphics[scale=0.3]{E_scheme_plot.pdf}
  }
\caption{Azimuthal asymmetry $A^{\sin(\phi_J+\phi_S)}$ in the jet-axis probe at the EIC in (a) the WTA scheme and (b) the E-scheme.
The data points with error bars are from event generator simulations. The solid lines in (a)
are from Eq.~(\ref{eq:A_transversity_2}). }
  \label{fg:asym_transversity_EIC}
 \end{center}
\end{figure}

Likewise, we can make predictions on the asymmetry that probes the Boer-Mulders function defined by
\begin{align}
A^{\cos(2\phi_J)}(|{\bm P}_{J\perp}|)&=\frac{2}{\int d\sigma \epsilon}\int d\sigma\cos(2\phi_J)\label{eq:A_Boer_1}\\
&=\frac{\langle \epsilon F_{UU}^{\cos(2\phi_J)}\rangle}{\bar{\epsilon} \langle F_{UU,T}\rangle}\,.\label{eq:A_Boer_2}
\end{align}
The structure function $F_{UU}^{\cos(2\phi_J)}$ can be written as
\begin{align}
F_{UU}^{\cos(2\phi_J)}&=-\frac{x}{M}\sum_a e_a^2\int \frac{d^2b}{(2\pi)^2}e^{-i{\bm P}_{J\perp}\cdot{\bm b}}
\left[\frac{2}{|{\bm P}_{J\perp}|^2}\left(P^i_{J\perp}\cdot\partial_i\tilde{h}^{\perp a}_1\right)\left(P^j_{J\perp}\cdot\partial_j\tilde{J}^{a}_T\right)-
\partial_i\tilde{h}^{\perp a}_1\partial_i\tilde{J}_T\right]
\,,\label{eq:F_UU_cos2}
\end{align}
where $\tilde{h}_1^\perp$ and $\tilde{J}_T$ are the Fourier transforms of $h_1^\perp$ and $J_T$ respectively. We adopt the Boer-Mulders
functions obtained from Ref.~\cite{Barone:2009hw}. The predictions on $A^{\cos(2\phi_J)}(|{\bm P}_{J\perp}|)$ at the EIC according to Eq.~(\ref{eq:A_Boer_2})
are shown in Fig.~(\ref{fg:asym_Boer_EIC}).

 \begin{figure}[htbp]
  \begin{center}
   \includegraphics[scale=0.3]{Asym_Boer_EIC.pdf}
   \caption{Azimuthal asymmetry $A^{\cos(2\phi_J)}$ in the jet-axis probe at the EIC as predicted from Eq.~(\ref{eq:A_Boer_2}). }
   \label{fg:asym_Boer_EIC}
 \end{center}
\end{figure}

The predictions of $A^{\sin(\phi_J+\phi_S)}$ and $A^{\cos(2\phi_J)}$ at the EicC are shown in Fig.~\ref{fg:asym_EICC}. As in Fig.~\ref{fg:asym_transversity_EIC} (a) and Fig.~\ref{fg:asym_Boer_EIC}, the data points with error bars are from event generator
simulations and the lines are from the factorization formulae.

 \begin{figure}[htbp]
  \begin{center}
 \subfigure[]{
   \includegraphics[scale=0.3]{Asym_transversity_EICC.pdf}
     }
   \subfigure[]{
   \includegraphics[scale=0.3]{Asym_Boer_EICC.pdf}
  }
\caption{Azimuthal asymmetries $A^{\sin(\phi_J+\phi_S)}$ (a) and $A^{\cos(2\phi_J)}$ (b) in the jet-axis probe at the EicC. The data points with error bars are from event generator
simulations and the lines are from the factorization formulae. }
  \label{fg:asym_EICC}
 \end{center}
\end{figure}

For the sake of comparison with data in SIDIS, in Fig.~\ref{fg:asym_Hermes} we show the asymmetries $A^{\sin(\phi_J+\phi_S)}(|{\bm P}_{J\perp}|)$
(a) and $A^{\cos(2\phi_J)}(|{\bm P}_{J\perp}|)$ (b) at HERA with predictions for jets from Eqs.~(\ref{eq:A_transversity_2}) and~(\ref{eq:A_Boer_2})
(dashed lines), predictions for pion production from the parallels of Eqs.~(\ref{eq:A_transversity_2}) and~(\ref{eq:A_Boer_2}) as
in~Refs.~\cite{Bacchetta:2006tn,Barone:2009hw} (solid lines), and data points for pion production from the HERMES experiment~\cite{HERMES:2010mmo,HERMES:2012kpt}
(data points with error bars). From Fig.~\ref{fg:asym_Hermes}, we see that the T-odd jet function
does give azimuthal asymmetries with sizes and shapes similar to those in SIDIS, and so should be observable
even at low-energy machines.

 \begin{figure}[htbp]
  \begin{center}
 \subfigure[]{
   \includegraphics[scale=0.3]{Asym_transversity_Hermes.pdf}
     }
   \subfigure[]{
   \includegraphics[scale=0.3]{Asym_Boer_Hermes.pdf}
  }
\caption{Azimuthal asymmetries $A^{\sin(\phi_J+\phi_S)}$ and $A^{\cos(2\phi_J)}$ at HERA. The data points with error bars are from the HERMES experiment for pion production. The solid lines are predictions from the factorization formulae for pion production. The dashed line are predictions from the factorization formulae for the jet-axis probe. }
  \label{fg:asym_Hermes}
 \end{center}
\end{figure}

\section{back-to-back dijet production in $e^+e^-$ annihilation} \label{sec:ee}

The T-odd jet function will give rise to novel jet phenomena in $e^+e^-$ annihilation, which are measurable at $e^+e^-$ machines.
For instance, consider back-to-back dijet production in $e^+e^-$ annihilation, as shown in Fig.~\ref{fg:ee}. We define
${\bm q}_T=-{\bm P}_{J_1\perp}$. The back-to-back limit corresponds to $|{\bm q}_{T}|\ll \sqrt{s} R$, where $R$ is the jet radius. The azimuthal asymmetry $A$~\cite{Kang:2014zza} is given by
\begin{align}
A &=  2 \int  d\cos\theta \, \frac{ d\phi_1}{\pi} \cos(2\phi_1)A^{J_1J_2}\,,\label{eq:A}
\end{align}
where
\begin{align}
 A^{J_1J_2} &=
 1+ \cos(2\phi_1) \frac{\sin^2\theta}{1+\cos^2\theta} \frac{F_T}{F_U}\,,   \label{eq:A_2}
\end{align}
with
\begin{align}
F_U &= |\bm{q}_T|\, \sum_q e_q^2\, \int \frac{{d}^2b}{(2\pi)^2} e^{i\bm{q}_T\cdot\bm{b}} \tilde{J}_1^q(b^2){ \tilde{J}}^{\bar q}_1(b^2) \,,
\label{eq:F_U}\\
F_T &=|\bm{q}_T| \, \sum_q \, e_q^2\, \int \frac{{d}^2b}{(2\pi)^2}e^{i\bm{q}_T\cdot\bm{b}}
\left( 2 \frac{q_T^i q_T^j}{|\bm{q}_T|^2}
- \delta^{ij}  \right) \partial_{b^i} \tilde{J}^q_{T}(b^2) \partial_{b^j}\tilde{J}^{\bar q}_{T}(b^2) \,.\label{eq:F_T}
\end{align}

 \begin{figure}[htbp]
  \begin{center}
   \includegraphics[scale=0.3]{ee.pdf}
\caption{Back-to-back dijet production in $e^+e^-$ annihilation.}
  \label{fg:ee}
 \end{center}
\end{figure}

In Fig.~\ref{fg:R_plot}, we plot the asymmetry $A$ as a function of $|\bm{q}_T|$ as predicted by
Eq.~(\ref{eq:A_2}) for four different values of $\sqrt{s}$. To enhance the sensitivity, we have demanded that $Q_J>0.25$ for one of the jets and $Q_J<-0.25$ for the other.
The value $\sqrt{s}=\sqrt{110}$~GeV corresponds to the Belle experiment. The values $\sqrt{s}=91.2$~GeV, $165$~GeV, and $240$~GeV correspond to the $Z$-threshold, the $W$-threshold, and the $Z$-Higgs-threshold respectively at LEP as well as the CEPC. One can see that the asymmetry is more significant at low-energy machines.

 \begin{figure}[htbp]
  \begin{center}
   \includegraphics[scale=0.4]{R_ee.pdf}
\caption{Azimuthal asymmetry $A$ for dijet production in $e^+e^-$ annihilation as predicted by Eq.~(\ref{eq:A_2}) at Belle, LEP and the CEPC, with $Q_J>0.25$ for one of the jets and $Q_J<-0.25$ for the other.}
 \label{fg:R_plot}
  \end{center}
\end{figure}

\section{Summary and outlook} \label{sec:summary}

In this work, we reinterpreted the jet clustering procedure as a way to define an axis, which together with the proton beam defines the transverse momentum of the vitual photon in DIS. In this way, one can use jet-axis measurements in DIS to probe the TMD PDFs of the nucleons, just like in the Drell Yan process. We provided the complete list of azimuthal asymmetries in the jet-axis probe in DIS at leading power. We showed that, by including the T-odd jet function in addition to the traditional one, all eight TMD PDFs of a nucleon at leading twist can be accessed by the jet-axis probe.  As concrete examples, within the WTA axis scheme, we demonstrated that with both event-generator simulations and predictions from the factorization formulae,
couplings of T-odd jet function with the quark transversity and the Boer-Mulders function give rise to sizable azimuthal asymmetries at DIS machines of various energy regimes, such as the EIC, the EicC, and HERA. We also demonstrated, with event-generator simulations, how the change of the jet-axis definition induces changes in the asymmetry distributions drastically.
We also gave predictions for the azimuthal asymmetry of back-to-back dijet production in $e^+e^-$ annihilation. The T-odd jet function has opened the door to a fully comprehensive study of nucleon 3D structure with jet probes. Further theoretical and phenomenological studies of the T-odd jet function, such as high-order corrections, evaluations of the soft function on the lattice, and fittings with experimental data, will empower the jet probe as a precision tool which is fully differential for the study of TMD physics.

\acknowledgments
W.~K.~L. and H.~X. are supported by
the National Natural Science Foundation of China under Grant No. 12022512, No. 12035007, and by the Guangdong Major Project of Basic and Applied Basic Research No. 2020B0301030008. X.~L. and M.~W. are supported by the National Natural Science Foundation of China under Grant No.~12175016.
W.~K.~L. acknowledges support by the UC Southern California Hub, with funding from the UC National Laboratories division of the University of California Office of the President.