# MHD turbulence formation in solar flares: 3D simulation and synthetic
observations

W. Ruan Centre for mathematical Plasma Astrophysics, Department of
Mathematics, KU Leuven,
Celestijnenlaan 200B, B-3001 Leuven, Belgium<EMAIL_ADDRESS>L. Yan
Key Laboratory of Earth and Planetary Physics, Institute of Geology and
Geophysics, Chinese Academy of Sciences  R. Keppens  Centre for mathematical
Plasma Astrophysics, Department of Mathematics, KU Leuven,
Celestijnenlaan 200B, B-3001 Leuven, Belgium

(Accepted Oct 17, 2022)

###### Abstract

Turbulent plasma motion is common in the universe, and invoked in solar flares
to drive effective acceleration leading to high energy electrons. Unresolved
mass motions are frequently detected in flares from extreme ultraviolet (EUV)
observations, which are often regarded as turbulence. However, how this plasma
turbulence forms during the flare is still largely a mystery. Here we
successfully reproduce observed turbulence in our 3D magnetohydrodynamic
simulation where the magnetic reconnection process is included. The turbulence
forms as a result of an intricate non-linear interaction between the
reconnection outflows and the magnetic arcades below the reconnection site, in
which the shear-flow driven Kelvin-Helmholtz Instability (KHI) plays a key
role for generating turbulent vortices. The turbulence is produced above high
density flare loops, and then propagates to chromospheric footpoints along the
magnetic field as Alfvénic perturbations. High turbulent velocities above 200
km s-1 can be found around the termination shock, while the low atmosphere
reaches turbulent velocities of 10 km s-1 at a layer where the number density
is about 1011 cm-3. The turbulent region with maximum non-thermal velocity
coincides with the region where the observed high-energy electrons are
concentrated, demonstrating the potential role of turbulence in acceleration.
Synthetic views in EUV and fitted Hinode-EIS spectra show excellent agreement
with observational results. An energy analysis demonstrates that more than 10%
of the reconnection downflow kinetic energy can be converted to turbulent
energy via KHI.

††journal: ApJ

##  1 Introduction

Solar flares actively drive space weather, which affects interplanetary space
and the Earth environment and atmosphere through energetic particles, strong
EUV/X-ray emissions and associated coronal mass ejections (CMEs). In a solar
flare event, up to 1033 erg of free energy stored in the solar magnetic field
can be released via the magnetic reconnection process, and up to half of this
energy is used to produce energetic electrons (Aulanier et al., 2013;
Aschwanden et al., 2017). Magnetic reconnection can develop at current sheets,
marking locations where a component of the magnetic field reverses (Parker,
1963). This changes the configuration of magnetic field lines, converts
magnetic energy to thermal and kinetic plasma energy and produces fast flows
that leave the reconnection site along the current sheet (a recent review see
Yamada et al. (2010)). In a typical flare event, reconnection is generally
believed to occur at coronal heights, and the downward reconnection outflows
lead to – and interact with – magnetic arcades in lower regions of the corona.
These arcades are filled with hot and dense ionized gas and they become bright
in the images at soft X-rays (SXR) and specific EUV passbands (Shibata et al.,
1995). Information on energetic electrons produced in solar flares is often
derived from hard X-rays (HXR, photons with energy $>20$ keV), as energetic
photons are believed to be produced by electrons via the bremsstrahlung
mechanism (e.g. Kontar et al., 2011). Strong and isolated HXR sources are
often observed (1) at the footpoints of the hot and dense arcade flare loops
at chromospheric heights (0.5 - 3 Mm above the photosphere) and (2) near the
top of the flare loops (e.g. Masuda et al., 1994; Tomczak & Ciborski, 2007; Su
et al., 2013). Hence, a large amount of energetic electrons of energy $>20$
keV must be produced during these explosive events (Hudson & Ryan, 1995;
Tomczak, 2001; Krucker et al., 2008).

In the study of solar flares, promising mechanisms to produce energetic
electrons are turbulence acceleration, as well as shock acceleration and
direct current (DC) electric field acceleration (Aschwanden, 2005). For the
former mechanism, unresolved mass motions with velocities exceeding 100 km s-1
have been discovered in flares for decades (e.g. Doschek et al., 1980; Gabriel
et al., 1981; Antonucci et al., 1982). The unresolved mass motions are often
thought of as turbulence, although it is impossible to know whether typical
energy cascades are happening from remote observations. Recent observations
with higher resolution data demonstrate that this plasma turbulence is not
really localized, and covers the entire flare region including flare loop top,
legs, footpoints and the region above the looptop (e.g. Doschek et al., 2014;
Kontar et al., 2017; Jeffrey et al., 2018; Stores et al., 2021). Peaks in
turbulence velocity values tend to show up above the high density flare loops,
reaching 100-200 km s-1, while footpoints have lower turbulent velocities of a
few tens km s-1. Turbulence has been frequently reproduced in 2D simulations,
but the obtained turbulence is more localized (e.g. Fang et al., 2016; Ruan et
al., 2018; Ye et al., 2019; Wang et al., 2022). Wang et al. (2022) invoke the
Kelvin-Helmholtz instability in their 2D settings, to explain the observed
wiggling of the current sheet above. In contributions to solar flare research,
we investigated (Ruan et al., 2020) the interplay between 2D
magnetohydrodynamic flare evolutions as coupled to analytic prescriptions of
the accelerated electron beams, where we could reproduce both HXR source
regions; and evolve the flare far into the postflare regime to reproduce
flare-driven coronal rain (Ruan et al., 2021). Recently, a fully 3D solar
flare simulation shows that turbulence can be produced in flares, with a clear
role played by a mixture of the Rayleigh-Taylor (RTI) and the Richtmyer-
Meshkov instability (RMI) at the interface between the reconnection
termination shock and the flare arcade, leading to finger-like supra-arcade
downflows (Shen et al., 2022). Our present 3D simulation will augment this
result, by looking especially at the spatial distribution and the values of
the non-thermal velocities obtained, and by making detailed comparisons to
observations. Our present 3D simulation goes one step further to investigate
the cause and consequence of the plasma turbulence in the entire flare region
by successfully reproducing the observational features of the plasma
turbulence, including the spatial distribution and characteristic velocities.

##  2 Method

The simulation is performed with the open-source MPI-AMRVAC code (Xia et al.,
2018; Keppens et al., 2021). The square simulation box has a domain of -50 Mm
$\leq x\leq$ 50 Mm, -50 Mm $\leq y\leq$ 50 Mm and 0 $\leq z\leq$ 100 Mm. The
minimal resolution is $64\times 64\times 64$, but an equivalent high
resolution of $1024\times 1024\times 1024$ is achieved via employing 5 levels
in our adaptive mesh refinement strategy, making our smallest cells less than
100 km across (the Shen et al. (2022) grid cell size was at 260 km). Gravity,
thermal conduction and optically thin radiative losses are included, where the
cooling curve comes from Colgan et al. (2008). A spatially varying, but
temporally invariant background heating is employed to balance the radiative
loss and maintain a corona. The governing equations are identical to Ruan et
al. (2020). A magnetic-field-line-based transition region adaptive conduction
(TRAC) method is adopted to ensure that underresolving the sharp transition
region variations – which would need better than 30 km resolution (Johnston et
al., 2020; Zhou et al., 2021) – does not lead to erroneous coronal temperature
and density evolutions.

The initial conditions for number density, temperature and background heating
are similar to our 2.5D flare simulation presented in Ruan et al. (2021),
where all initial profiles are functions of height $z$ only. The initial
coronal region has an electron density of order $\sim 10^{9}$ cm-3 and
temperature of $2-3$ MK. The number density and temperature profiles are
obtained from a relaxation in which the model C7 temperature profile in Avrett
& Loeser (2008) and a density profile calculated based on hydrostatic
equilibrium are employed as initial conditions, where the number density at
the bottom boundary is 3.7 $\times$ 1014 cm-3. The initial conditions for
magnetic field are modified from those in Ruan et al. (2021), which is given
by

$\displaystyle B_{x}$ $\displaystyle=$
$\displaystyle\sqrt{B_{0}^{2}-B_{z}^{2}},$ (1) $\displaystyle B_{y}$
$\displaystyle=$ $\displaystyle 0,$ (2) $\displaystyle B_{z}$ $\displaystyle=$
$\displaystyle\begin{cases}-B_{0},&\quad y<-\lambda\\\ B_{0},&\quad
y>\lambda\\\ B_{0}\sin[\pi y/(2\lambda)],&\quad else\\\ \end{cases}$ (3)

where $B_{0}=30$ G is the initial magnetic field strength and $\lambda=10$ Mm.
We dissipate the thick current sheet into a thin current sheet in the pre-
flare phase by adopting a resistivity given by

$\eta(y,z)=\eta_{1}\exp(-y^{2}/w_{\eta y}^{2})\exp[-(z-h_{\eta 1})^{2}/w_{\eta
z}^{2}]\,,$ (4)

where $\eta_{1}=10^{-2}$, $w_{\eta y}=10\ \rm Mm$, $h_{\eta 1}=30\ \rm Mm$ and
$w_{\eta z}=15\ \rm Mm$. Magnetic arcades are produced in the lower atmosphere
during the slow dissipation of the thick current sheet, without generating
strong chromospheric evaporation. Hot and dense flare loops will only be
generated above the magnetic arcades in the upcoming impulsive phase, that
leads to formation of clear bright loops in the high temperature EUV images
(e.g. 131 Å). The anomalous resistivity prescription is changed when a thin
current sheet is formed at $t=450$ s. The new resistivity adopted then is
given by

$\eta(x,y,z)=\\\ \begin{cases}\eta_{2}+\eta_{3}\exp\\{-[x^{2}+y^{2}+(z-h_{\eta
2})^{2}]/r_{\eta}^{2}\\},&480\ \mathrm{s}\geq t>450\ \mathrm{s}\\\
\eta_{2},&\quad t>480\ \mathrm{s}\end{cases}$ (5)

where $\eta_{2}=10^{-3}$, $\eta_{3}=10^{-1}$, $h_{\eta 2}=50\ \rm Mm$ and
$r_{\eta}=5\ \rm Mm$. Fast reconnection is triggered by this localized strong
resistivity at the center of our thinned current sheet. The magnetic field
inside the initial current sheet contains an $x$ component, therefore the
earliest formed magnetic arcades (mainly located inside/below the dense flare
loop) also contain an $x$ component, and hence have a shear angle to the
magnetic neutral line. The magnetic arcades generated at the impulsive phase
are almost in $y$-$z$ planes, as the $x$ component of the magnetic field
inside the reconnection current sheet has been transported to the lower
atmosphere and the region above the simulation box.

Periodic boundaries are employed at the $x$ boundaries. Symmetric boundary
conditions are adopted for density, pressure and $x$/$z$ components of
velocity/magnetic field at the $y$ boundaries, while anti-symmetric conditions
are used for $y$-components of velocity and magnetic field. Fixed boundary
conditions are employed at the bottom boundary. We employ zero gradient
extrapolation at the (two-layer) ghost cells of the upper boundary for
density, velocity and magnetic field. The temperature at the top ghost cells
is forced to decrease at a rate $dT/dz=-T_{b}/$(50 Mm) to avoid runaway high
temperatures at the boundary, where $T_{b}$ is the instantaneous local
temperature at the boundary. A combination of the ‘HLL’ (initials of authors
Harten, Lax and van Leer in Harten et al. (1983)) approximate Riemann solver
and ‘Cada3’ (first author from Čada & Torrilhon (2009)) flux limiter is
employed at the low level ($\leq 3$) grids, as the third-order limiter ‘Cada3’
achieves high accuracy at low spatial resolution. Another combination, ‘HLL’
with ‘Vanleer’ (from van Leer (1974)) limiter, is adopted at high level grids
located at the low atmosphere and flaring regions, as the second-order
‘Vanleer’ limiter has better performance in high gradient regions. A strong
stability preserving, three-step Runge-Kutta method is employed in time
integration (Ruuth, 2006).

Contribution functions from the CHIANTI atomic database have been used in
synthesizing EUV emissions (Del Zanna et al., 2015). The spatial resolution of
the synthesized EUV images is the same as given by the observation (using a
pixel size of 435 km for image at 131 Å passband). The synthesized 255 Å
spectra are assumed to have a slit width of 1450 km and a pixel resolution of
725 km $\times$ 22 mÅ. The 131 Å line is mainly related to emission by Fe VIII
and Fe XXI ions with peak formation temperatures of 105.6 K and 107.0 K, and
the 255 Å line is mainly associated with Fe XXIV ions with a peak formation
temperature of 107.2 K (Culhane et al., 2007; Lemen et al., 2012). Scattering
effects by the instrument have been included via multiplying by a Point Spread
Function (PSF), where we assume the PSF is a Gaussian function of standard
deviation of 1 pixel (Grigis et al., 2013). Note that only the contribution of
hot plasma ($>3$ MK) is considered in synthesizing 131 Å emission, to avoid
unrealistic strong emission at the low atmosphere due to underresolved
transition region physics. The GOES SXR flux is calculated with the method
given in Pinto et al. (2015).

##  3 Results

###  3.1 Bright EUV arcades evolution in flaring region

Our simulation starts with a single vertical ($x-z$ oriented) current sheet
all along $y=0$, which runs through the chromosphere and the corona. Magnetic
reconnection happens inside this current sheet at coronal heights, leading to
formation of an extended coronal magnetic arcade system below the reconnection
site. Some magnetic arcades are filled with hot ($\sim 10$ MK) and dense
plasma ($\sim 10^{10}$ cm-3) by downward flows from the reconnection site and
(thermal conduction driven) upward flows from the chromosphere, respectively.
This then forms high density flare loops, which are bright in images at the
EUV 131 Å passband. Photons at this passband are mainly released by Fe XXI
ions at temperatures $\sim$10 MK in flare events. Figure 1 gives synthesized
images at 131 Å passband of the simulated flare loops, where for comparison,
an actual flare event observed by the Atmospheric Imaging Assembly (AIA)
onboard spacecraft Solar Dynamics Observatory (SDO) is also given (Lemen et
al., 2012). The axes orientation is shown in both top panels, where $Z$ is the
vertical direction in our simulation, the bottom of the chromosphere being
located at $Z=0$. We indicated a reference bar of 10 arcsec length, to better
compare with observations, noting that 1 arcsec $=$ 725 km. Here we use
capital letters ($X$, $Y$, $Z$) to indicate directions when the length unit is
arcsec and use lower case letters ($x$, $y$, $z$) when the length unit is Mm.
The hot flare loops start to form at $t\approx 7$ min due to the triggering of
fast reconnection, and then the flare enters an impulsive phase. Here we focus
on the generation and development of turbulence in the impulsive phase, at the
time after the formation of the flare loops. Most of the figures in our paper
use the data at $t\approx 9$ min (shown with a black dashed line in Fig. 1d),
when finger-like structures have not yet appeared. The wide initial current
sheet is slowly dissipated by resistivity and then becomes a thin current
sheet in the period $t<7$ min. Magnetic arcades are also formed in this
period, but without the generation of strong evaporations. Consequently, there
are no bright 131 Å loops or strong SXR emission during this period.

The finger-like structures, which are caused by RTI/RMI according to Shen et
al. (2022), also appear in our simulation, but clearly at a phase later than
the one we analyze in detail here (see Fig. 1d). RMI is an instability that
often happens at the surface between two fluids of different densities, which
can be regarded as a special type of RTI (Richtmyer, 1960; Meshkov, 1969; Zhou
et al., 2021). The ‘traditional’ RTI is caused by a constant acceleration
(e.g. due to external gravity), and the initial perturbation grows
exponentially, whereas the RMI is the result of an impulsive acceleration
(e.g. caused by a shock wave), and the initial perturbation grows linearly.
Both RTI/RMI instabilities may end up in turbulent dynamics, in which
additionally KHI is involved for generation of vortices.

Figure 1: (a) & (b): Synthesized 131 Å views of the flare loop systems
obtained in our simulation. The two views are synthesized with different line
of sight (LOS) directions, where the corresponding orientation axes are given.
The red curves in (b) are projections of several magnetic field lines. The
approximate locations of the magnetic reconnection X-point (white cross) and
the reconnection outflows (white arrows) are also given in (b). Panel (c)
shows a 131 Å image of a flare event from Dec 26, 2011, where the LOS
direction is similar to that in (a), as reported by Cheng & Qiu (2016). Panel
(d) gives the time development of synthesized GOES SXR flux at 1-8 Å passband,
where the black dashed line gives the corresponding time of the synthesized
views. A cross-section in the $X-Z$ plane, in a zoomed-in view on the top of
the arcade, shown in panel (d), confirms the RTI/RMI process from Shen et al.
(2022), happening at a later time referring to the blue dotted line in this
panel. An animation of this figure is available. Panel (a)-(c) of the
animation give synthesized 131 Å views obtained with different LOS, where the
corresponding orientation axes are given. Panel (d) of the animation shows the
synthesized GOES SXR flux at 1-8 Å, the same as panel (d) of this figure. The
animation covers $\sim 10$ minutes of physical time starting at $t=0$ (real-
time duration 14 s).

###  3.2 Where is the turbulence?

Figure 2: (a) Non-thermal velocity distribution obtained from Gaussian fitting
of synthesized EIS/255 spectra, where the fitting method refers to Stores et
al. (2021). (b) Non-thermal velocity distribution obtained from plasma
density/velocity distribution. Cyan contours in (b) give the locations of the
bright loop in Fig. 1b, where the contour lines show the intensity levels 25%,
50% and 80% of peak 131 Å flux. The regions inside the blue contour have an
average magnetic field strength lower than 25 G. The green line gives the
approximate location of termination shock (also see Fig. 5c), but note that
the shock locations are different at different X-slices due to interactions
between reconnection outflows and magnetic arcades. Panel (c) gives the
EIS/255 spectrum inside the blue box in (a) and (d) gives the spectrum inside
the green box in (a) (marked with ‘+’). Solid lines show Gaussian fitting
results of the spectra. Panel (e) shows a non-thermal velocity map from the
EIS/255 spectra of an observed flare, where the cyan contours give the
location of hot/dense flare loop and white contours give the location of
footpoints at the solar surface, as taken from Stores et al. (2021). The
spatial resolution of the synthesized EIS/255 data is reduced before fitting,
as in Stores et al. (2021). Panel (a) has a pixel size of 4 arcsec $\times$ 4
arcsec, which is close to that in panel in (e) (6 arcsec $\times$ 4 arcsec).
Panel (b) has a smaller pixel size of 0.138 arcsec $\times$ 0.138 arcsec. An
animation of this figure is available, showing the evolution of spatial
average electron-number density in $x$-direction (panel a) and non-thermal
velocity distribution obtained from plasma density/velocity distribution
(panel b). The contours in the panels show the intensity levels 10%, 25%, 50%
and 80% of peak 131 Å flux. The animation covers 2.37 minutes of physical time
starting at t = 7.51 minutes ( real-time duration 3 s).

We obtain the turbulent, non-thermal velocities from synthesized spectral
profiles of the emission line at 255.1136 Å, corresponding to a peak
temperature 107.2 K, as in Stores et al. (2021). The non-thermal velocity
distribution derived from the spectral profiles is given in Fig. 2a, where the
observational result from Stores et al. (2021) is also replicated in Fig. 2e.
The observational data were obtained with the EUV Imaging Spectrometer (EIS)
onboard Hinode (Culhane et al., 2007). The synthesized spectral data adopted
the same spatial and wavelength resolutions as the observational data, and the
instrument effect has also been included. In order to obtain a clear view on
the distribution of turbulence relative to the high density flare loop, we
employ a LOS parallel to our $x$ direction (i.e. perpendicular to our flare
loops) when synthesizing these EIS/255 spectra. Figure 2a gives the non-
thermal velocity distribution across the $y-z$ plane obtained from the EIS
spectra, while a higher resolution equivalent calculated from the available
magnetohydrodynamic (MHD) plasma parameters (electron number density
$N_{\mathrm{e}}$ and $x$-component of velocity $v_{x}$) by means of

$v_{\mathrm{nth}}(y,z)=\sqrt{\frac{\int
N_{\mathrm{e}}(x,y,z)v_{x}^{2}(x,y,z)\mathrm{d}x-[\int
N_{\mathrm{e}}(x,y,z)v_{x}(x,y,z)\mathrm{d}x]^{2}/\int
N_{\mathrm{e}}(x,y,z)\mathrm{d}x}{\int N_{\mathrm{e}}(x,y,z)\mathrm{d}x}}\,,$
(6)

is shown in panel b. The 131 Å contours indicate the location of the flare EUV
loops.

The non-thermal velocity distribution from panels a/b shows very similar
features to the observational distribution in panel e, where the turbulence
has a wide spatial distribution above the AIA/131 flare loops and the maximum
non-thermal velocity is located above the looptop. The maximum non-thermal
velocity is located at a weak magnetic field region (see panel b), which is a
promising electron acceleration site suggested by Chen et al. (2020) and
observationally demonstrated in Fleishman et al. (2022). The spectral profiles
of the high non-thermal velocity regions have perfect Gaussian distributions
as demonstrated in Fig. 2 panels c and d, which is a distinctive feature of
flare turbulence. We confirm from the high resolution non-thermal velocity map
(panel b) that there are two isolated high $v_{\mathrm{nth}}$ regions above
the high density flare loop, which also show up in panel a. The high non-
thermal velocity below the loops in panel a is mainly caused by chromospheric
evaporations related to shocks rather than by turbulence.

###  3.3 How is the turbulence produced?

In order to investigate how the turbulence is produced, we analyze the
dynamics in a slice (parallel to the $x$-$z$ plane at $y=3$ Mm or $Y\sim 4$
arcsec) that runs across a region with maximal non-thermal velocities (as
indicated by the white dashed vertical line in Fig. 2b). For that vertical
plane, Figure 3a (top left panel) shows the localized AIA/131 emission flux
distribution as well as the ‘velocity field’ (through streamlines of
$v_{x}\vec{e}_{x}+v_{z}\vec{e}_{z}$) on the slice. On its right hand side, the
same panel Figure 3a gives the non-thermal velocity distribution as function
of the height as obtained from Fig. 2b. It is obvious in the velocity
streamline view that there are a lot of vortices around the height with
maximal non-thermal velocity (i.e. at $Z\approx 35$ arcsec).

Figure 3b-e (four top right panels) show in colorscales the spatial
distributions of number density $N_{\mathrm{e}}$, $v_{x}$, $v_{z}$, and
thermal pressure $P$ (the latter in dimensionless unit) of a selected
turbulent region: namely in the red box shown in Fig. 3a. We each time overlay
the streamlines as well. Vortices appear prominently near $Z\sim 35$ arcsec,
where we find large shear velocities (e.g. note the velocity jump $>1000$ km
s-1 in $v_{z}$ from $X\sim 21$ arcsec to $X\sim 23$ arcsec), which indicates
the turbulence is produced via KHI. KHI may happen at the interface between
high speed shear flows to produce vortices at their interface. This
instability is often inhibited by magnetic tension in magnetized plasmas. When
the flow velocity vector is parallel to the magnetic field lines, triggering
KHI requires that (half) the velocity jump is larger than the local Alfvén
speed (Keppens et al., 1999). Here, in the slice analysed, we find conditions
that are quite different from this standard 2D situation. We have here a
guiding magnetic field that is locally almost perpendicular to the $x$-$z$
plane shown here (shape of the guiding field refers to the white solid line in
Fig. 4a), while there are large shears in both $v_{x}$ and $v_{z}$. As a
result, KHI vortices can be produced on this $x$-$z$ plane without distorting
the magnetic field too much, and hence without producing a stabilizing
magnetic tension. Such a condition is unconditionally unstable for KHI. The
relative orientation of magnetic field and flow shear is similar in the case
of so-called TWIKH rolls from transverse-wave-induced KHI in coronal loops
(e.g. Antolin et al., 2014; Guo et al., 2019; Shi et al., 2021). A simplified
growth rate estimate of KHI in such a condition $\omega=kv_{0}$ (using
equation 13.38 in Goedbloed et al., 2019) gives a growth timescale
($1/\omega$) of seconds, the same order of magnitude as the vortex generation
time scale seen in our simulation. Here, we used $k=1/(\mathrm{several\ Mm})$
since our vortices have a length scale of Mm and a half shear speed $v_{0}$ of
hundreds up to one thousand km s-1, as found in the flow distribution.

Figure 3: Color map in top left panel (a): AIA/131 emission flux distribution
at slice $\mathrm{Y}\sim 4$ arcsec, along the white dashed vertical line in
Fig. 2b. At the RHS of panel(a) we show non-thermal velocity. Top right panels
(b)-(e): $N_{\mathrm{e}}$, $v_{{x}}$, $v_{{z}}$ and thermal pressure inside
the red box in panel (a). Bottom left panels (f)&(g): $N_{\mathrm{e}}$ and
$v_{{z}}$ at a slice marked with a green dashed line in (b)-(e). Bottom right
panels (h)&(i): $N_{\mathrm{e}}$ and $v_{{z}}$ at a slice marked with the cyan
dashed line in (b)-(e). The red vertical dashed lines in all bottom panels
(f)-(i) give the location of the slice shown in top right panels (b)-(e).
Streamlines in (a)-(e) show the velocity field given by
$v_{x}\vec{e}_{x}+v_{z}\vec{e}_{z}$, while that in (f)-(i) show the field
given by $v_{y}\vec{e}_{y}+v_{z}\vec{e}_{z}$. The width of stream lines in
panel (a) is in proportion to $\sqrt{v_{x}^{2}+v_{z}^{2}}$ to highlight the
turbulent regions, while that in other panels is constant.

The initial shear motions come from non-linear interactions observed in the
downward reconnection flows with the magnetic arcades above the AIA/131 loop.
Collisions between the downflows and the arcades leads to the formation of a
termination shock and complicated reflected upflows. Their spatial
distribution on the vertical $y$-$z$ plane is different for different $x$ (as
shown in the bottom panels Fig. 3g & i), which enhances shear motion between
downward flows and upward reflection flows. Combining all views collected in
Fig. 3 we see that vortices prefer to form on $x$-$z$ planes, which is
perpendicular to the local magnetic field direction. A related, detailed study
on the interaction between reconnection flows and their arcades in 2D
conditions is found in Ye et al. (2021). The magnetic field strength inside
the flaring region is spatial-dependent, where the legs of the guiding line
shown in Fig. 4a have an average strength of $\sim 40$ G, while the apex of
the line has a weaker strength of less than 20 G due to the impact of
turbulent motion. In real flares, the magnetic field can be much stronger,
with field strength higher than 100 G above the 131 Å loops. However,
turbulence can still be produced with the mechanism we mentioned above, as it
is difficult for the magnetic field to inhibit the growth of KHI when the
flows velocity vectors are perpendicular to the field lines.

Note further that the turbulent upper border of the AIA/131 loops in Fig. 3a
is attributed by Shen et al. (2022) to RTI/RMI effects. These authors focused
on the dark fingers as also seen in an as yet early stage in our Fig. 3a,
where the connection with Supra-Arcade Downflows was made by Shen et al.
(2022). The RTI/RMI region does not lead to large turbulent velocities during
our simulation and is clearly situated below the region of maximal non-thermal
velocities. As demonstrated in Fig. 3a, the non-thermal velocity around the
upper border of the AIA/131 loops is smaller than 50 km s-1. Our (higher
resolution) simulation demonstrates that the turbulent zone is foremost
located above the EUV loops (see also Fig. 2) and most likely driven by KHI
mechanisms.

###  3.4 Turbulence at the lower atmosphere

Turbulence can also be found at lower atmospheric regions in our simulation.
The bottom slice in the 3D view at left in Fig. 4a demonstrates the $v_{x}$
distribution at a layer at height $z=2$ Mm. Turbulent motions up to $\sim$ 10
km s-1 appear at the footpoints of the high density arcades. These are shown
in the vertical slice in the 3D view of Fig. 4a, in their number density
distribution at $x=30$ Mm ($X\sim$ 41 arcsec). The turbulent region at the
$z=2$ Mm layer has an average number density of $\sim$1011 cm-3 and an average
temperature of $\sim$1 MK. The temperature actually drops rapidly from several
MK to a typical chromospheric temperature of $\sim$0.01 MK around this height.
Turbulent motions can still be found at a lower $z=1$ Mm layer where the
average number density is $\sim$1013 cm-3, but with reduced speeds by one
order of magnitude compared to those at $z=2$ Mm.

Figure 4: Horizontal slice in the left 3D view panel (a): $v_{x}$ distribution
at $z=2$ Mm, while the vertical slice shows the $N_{\mathrm{e}}$ distribution
at $x=30$ Mm. Right panels (b)-(e): Time-space plots of $B_{x}$, $v_{x}$,
Alfvén speed $v_{\mathrm{A}}$ and $N_{\mathrm{e}}$ along the white solid
(field) line shown in (a), where the minimum $s$-value is located at the left
end of this line. The horizontal dashed lines in these panels (b)-(e) and the
vertical dotted line in Fig. 1d give the corresponding time for the 3D view
(a) at left.

We find that the turbulence at this lower atmosphere region gets propagated
from higher regions downwards, rather than being generated locally. We select
a (time-invariant) curve which gives a general direction of the magnetic field
(the white solid line in Fig. 4a) and study the time development of typical
MHD quantities along that curve. The time-space plots of $B_{x}$ and $v_{x}$
along the curve (top right panels) demonstrate that there are structures
propagating from the middle of this curve to its ends. The $B_{x}$ and $v_{x}$
are seen to show an anti-phase when structures are propagating toward $s>0$,
while they are in phase when the structures are propagating toward $s<0$.
Since the magnetic field vector is pointing from $s<0$ to $s>0$, such a phase
relationship indicates that the structures are Alfvénic perturbations. In
fact, the propagation speeds (fitted speeds are indicated in Fig. 4b) of the
structures are close to the local Alfvén speeds which are quantified in panel
(d). Sudden decreases in propagating speed around $s=30$ Mm are a result of a
change in the local Alfvén speed, where this change in Alfvén speed is caused
by plasma density variation that is due to chromospheric evaporation
processes, as evident from the time-space plot of number density in panel (e).
The downward propagating Alfvénic perturbations and the turbulence in the
chromosphere may contribute to the generation of energetic electrons (e.g.
Fletcher & Hudson, 2008). The Alfvénic perturbations may bring a lot of energy
to the lower atmosphere, and contribute to the heating of the chromosphere and
the generation of evaporations (e.g. Russell & Fletcher, 2013; Reep & Russell,
2016). Those effects are worth investigating in future studies.

The bigger picture obtained from all previous sections is as follows:
turbulence gets produced above the looptops due to KHI, and this turbulence
gets propagated along magnetic fields to lower regions, leading to a spatially
distributed turbulent plasma at all heights. The arcade-shape spatial
distribution of the non-thermal velocity regions seen in Fig. 2a and b also
supports this bigger picture, in line with the arcade-shape magnetic field
that forms during the flare process.

###  3.5 Turbulent Kinetic Energy

When connecting MHD turbulence with actual non-thermal particle acceleration,
we should still verify whether the released magnetic energy from reconnection
can be efficiently converted to turbulence energy. For KHI driven turbulence,
the energy in the turbulence comes from the kinetic energy of bulk flows, the
reconnection downflows in particular. Therefore, here we compare the time-
evolving turbulence energy with the time-integrated reconnection downflow
kinetic energy.

Turbulence energy consist of kinetic energy of the turbulent velocity field
and magnetic energy of the turbulent magnetic field. The spatial distribution
of the kinetic energy density and the magnetic energy density at $t\approx 9$
min is shown in Fig. 5a and b, respectively. The turbulent kinetic energy
density is calculated from

$Ev_{\mathrm{tur}}(y,z)=\int\rho(x,y,z)v_{x}^{2}(x,y,z)\mathrm{d}x-[\int\rho(x,y,z)v_{x}(x,y,z)\mathrm{d}x]^{2}/\int\rho(x,y,z)\mathrm{d}x,$
(7)

and the turbulent magnetic energy density is calculated from

$Eb_{\mathrm{tur}}(y,z)=\int B_{x}^{2}(x,y,z)\mathrm{d}x-[\int
B_{x}(x,y,z)\mathrm{d}x]^{2}/\int\mathrm{d}x,$ (8)

where $\rho$ is density. An assumption applied in the calculation is that the
turbulence is anisotropic (Alfvénic), such that the turbulent motion is (only)
freely developed in directions perpendicular to the (average) magnetic field,
which is why we only incorporate $x$-components in both flow and magnetic
field turbulent quantifications for this $y-z$ view. Fig. 5a and b demonstrate
that $Ev_{\mathrm{tur}}$ and $Eb_{\mathrm{tur}}$ have similar density and
spatial distribution, again supporting our assumption that the turbulence is
indeed Alfvénic. The time-integrated kinetic energy of the reconnection
downflow is calculated via integrating the kinetic energy flux that goes
across a slice at height $z=45$ Mm ($Z\sim 62$ arcsec, the black dashed line
in Fig. 5c), where we do the integration from $t=450$ s, the time where the
fast reconnection regime starts. The color map shows the spatial average
vertical velocity distribution given by

$\bar{v}_{z}(y,z)=\int N_{\mathrm{e}}(x,y,z)v_{z}(x,y,z)\mathrm{d}x/\int
N_{\mathrm{e}}(x,y,z)\mathrm{d}x,$ (9)

which gives locations of the reconnection downflow.

Figure 5: (a) Spatial distribution of turbulent kinetic energy. (b) Spatial
distribution of turbulent magnetic energy. (c) X-averaged vertical speed
distribution. Taken along the black solid line from this panel (c), we show in
(d) the time-integrated reconnection downflow kinetic energy (black solid
line) that goes through the surface. Panel (d) further shows in blue solid
line the instantaneous total turbulence energy, and its division over kinetic
energy and turbulent magnetic energy. In this panel (d), the vertical red
solid line corresponds to the time for panels (a), (b) and (c). As in Fig. 2,
the contours in panels (a),(b) and (c) show the AIA/131 loop location.

Figure 5d gives a comparison of the (time) integrated reconnection downflow
kinetic energy $Ek$, the instantaneous total turbulence energy
$E_{\mathrm{tur}}$, and the kinetic and magnetic components of this turbulence
energy ($E_{\mathrm{tur}}=$$Ev_{\mathrm{tur}}$ $+$ $Eb_{\mathrm{tur}}$).
$Ev_{\mathrm{tur}}$ closely follows $Eb_{\mathrm{tur}}$ during the entire time
evolution as expected, a feature of Alfvénic motions. $E_{\mathrm{tur}}$ shows
a similar tendency as $Ek$, where we find $E_{\mathrm{tur}}\approx 0.1Ek$, so
roughly a 10 % fraction at each time. The total, time-integrated amount that
is transferred from bulk kinetic energy of the reconnection downward flows to
turbulence energy is more than the instantaneous turbulence energy $Ek$,
considering that turbulence energy is continuously (transported to the
chromosphere and) diffused.

##  4 Conclusions

The spatial distribution and the characteristic velocity of the plasma
turbulence in flare observation are successfully reproduced in our 3D MHD
flare simulation. The maximum turbulent velocity reached is higher than 200 km
s-1, and is located in close connection to the reconnection termination shock.
Reconnection downflows collide with magnetic arcades near the termination
shock and then produce complicated reflection upflows. Non-linear interactions
between the downflows and upflows lead to spatially extended regions where
high speed turbulent motion forms naturally, in which KHI dominates. More than
10% of the downflow kinetic energy is converted into turbulence energy via
this mechanism. The turbulence is anisotropic due to magnetic tension, as the
turbulent vortices only appear most clearly inside the plane that is
perpendicular to the guiding magnetic field. The locally generated turbulence
can propagate along the guiding magnetic field as Alfvénic perturbations, and
then leads to a wider spatial distribution of turbulent motions. The lower
atmosphere, where the number density reaches $10^{11}$ cm-3, still has
turbulent motions of order 10 km s-1.

Our model works for any flare that involves magnetic reconnection, no matter
how big the flare is or how long it takes to release energy, since it is
governed by the scale-invariant MHD description. Downwards reconnection
outflows invariably encounter lower-lying magnetic fields, by which the
outflows are stopped and diverted in the direction perpendicular to the
magnetic field, making the flow shear orthogonal to the field. Such a
condition is a perfect environment for KHI in full 3D (and missed in 2D
settings), regardless of the magnetic field strength. KHI can be easily
triggered for various reasons in such a condition, where the interaction
between the downflows and reflected upflows is (only) one of them.

Our simulations used various algorithmic improvements that render it possible
to resolve details down to few 100 km, thanks to adaptive mesh refinement, and
the TRAC treatment to properly evolve multi-dimensional MHD from chromosphere
to corona. Application of this TRAC method is not crucial for this KHI
turbulence at the studied impulsive phase, but will be important for properly
generating coronal rain in the gradual phase (Ruan et al., 2021). A 3D
numerical study on the formation of post-flare coronal rain can serve as
interesting follow-up to this work. Then, the role of other instabilities
(like RTI/RMI and thermal instabilities) in various evolutionary stages of a
3D flare loop can be clarified in detail. The maximum non-thermal velocity is
located at a weak magnetic field region, where the high-energy electrons are
accelerated. It indicates that plasma turbulence probably plays an essential
role in accelerating flare electrons.

The next step of this study should address the consequences of the fluid
turbulence for charged particle dynamics more seriously, either through test
particle assessment of the attainable particle acceleration efficiency, or
even evolving to hybrid (fluid/particle) or fully kinetic models. There is no
doubt that turbulence can accelerate electrons/ions to high energies, but how
much of the energy released by reconnection can go to energetic electrons/ions
in this way is still a question. According to Aschwanden et al. (2017), about
half of the energy will go to the energetic electrons and about 17% of the
energy will go to the energetic ions on average. Emslie et al. (2012) quotes
lower acceleration efficiency, but still about 20% of the released energy
would go to the energetic electrons and ions. It might be difficult to achieve
this in a multi-step process such as turbulence or shock acceleration
(Cargill, 1996; Miller et al., 1997). Most likely, the tens of percent energy
transfer efficiency is the result of multiple acceleration mechanisms
cooperating (e.g., turbulence acceleration and shock acceleration). Note that
other scenarios can achieve high acceleration efficiencies, such as the
acceleration at fragmented current sheets (Cargill et al., 2012).

We thank the referee Peter Cargill for very constructive comments. WR was
supported by a postdoctoral mandate (PDMT1/21/027) by KU Leuven. LY was
supported by the Youth Innovation Promotion Association of CAS (2021064). RK
is supported by Internal funds KU Leuven through the project C14/19/089
TRACESpace, and an FWO project G0B4521N. RK also received funding from the
European Research Council (ERC) under the European Union Horizon 2020 research
and innovation programme (grant agreement No. 833251 PROMINENT ERC-ADG 2018).

## References

  * Antolin et al. (2014) Antolin, P., Yokoyama, T., & Van Doorsselaere, T. 2014, ApJ, 787, L22
  * Antonucci et al. (1982) Antonucci, E., Gabriel, A. H., Acton, L. W., et al. 1982, Sol. Phys., 78, 107
  * Aschwanden (2005) Aschwanden, M. J. 2005, Physics of the Solar Corona. An Introduction with Problems and Solutions (2nd edition)
  * Aschwanden et al. (2017) Aschwanden, M. J., Caspi, A., Cohen, C. M. S., et al. 2017, ApJ, 836, 17
  * Aulanier et al. (2013) Aulanier, G., Démoulin, P., Schrijver, C. J., et al. 2013, A&A, 549, A66
  * Avrett & Loeser (2008) Avrett, E. H., & Loeser, R. 2008, ApJS, 175, 229
  * Cargill (1996) Cargill, P. 1996, EOS Transactions, 77, 353
  * Cargill et al. (2012) Cargill, P. J., Vlahos, L., Baumann, G., Drake, J. F., & Nordlund, Å. 2012, Space Sci. Rev., 173, 223
  * Chen et al. (2020) Chen, B., Shen, C., Gary, D. E., et al. 2020, Nature Astronomy, 4, 1140
  * Cheng & Qiu (2016) Cheng, J. X., & Qiu, J. 2016, ApJ, 825, 37
  * Colgan et al. (2008) Colgan, J., Abdallah, J., J., Sherrill, M. E., et al. 2008, ApJ, 689, 585
  * Culhane et al. (2007) Culhane, J. L., Harra, L. K., James, A. M., et al. 2007, Sol. Phys., 243, 19
  * Del Zanna et al. (2015) Del Zanna, G., Dere, K. P., Young, P. R., Landi, E., & Mason, H. E. 2015, A&A, 582, A56
  * Doschek et al. (1980) Doschek, G. A., Feldman, U., Kreplin, R. W., & Cohen, L. 1980, ApJ, 239, 725
  * Doschek et al. (2014) Doschek, G. A., McKenzie, D. E., & Warren, H. P. 2014, ApJ, 788, 26
  * Emslie et al. (2012) Emslie, A. G., Dennis, B. R., Shih, A. Y., et al. 2012, ApJ, 759, 71
  * Fang et al. (2016) Fang, X., Yuan, D., Xia, C., Van Doorsselaere, T., & Keppens, R. 2016, ApJ, 833, 36
  * Fleishman et al. (2022) Fleishman, G. D., Nita, G. M., Chen, B., Yu, S., & Gary, D. E. 2022, Nature, doi:10.1038/s41586-022-04728-8
  * Fletcher & Hudson (2008) Fletcher, L., & Hudson, H. S. 2008, ApJ, 675, 1645
  * Gabriel et al. (1981) Gabriel, A. H., Phillips, K. J. H., Acton, L. W., et al. 1981, ApJ, 244, L147
  * Goedbloed et al. (2019) Goedbloed, H., Keppens, R., & Poedts, S. 2019, Shear flow and rotation (Cambridge University Press), 473–524
  * Grigis et al. (2013) Grigis, P., Yingna, S., & Weber, M. 2013, AIA PSF characterization and image deconvolution, Tech. rep., Tech. Rep., AIA team
  * Guo et al. (2019) Guo, M., Van Doorsselaere, T., Karampelas, K., & Li, B. 2019, ApJ, 883, 20
  * Harten et al. (1983) Harten, A., Lax, P. D., & Leer, B. v. 1983, SIAM Review, 25, 35
  * Hudson & Ryan (1995) Hudson, H., & Ryan, J. 1995, ARA&A, 33, 239
  * Jeffrey et al. (2018) Jeffrey, N. L. S., Fletcher, L., Labrosse, N., & Simões, P. J. A. 2018, Science Advances, 4, 2794
  * Johnston et al. (2020) Johnston, C. D., Cargill, P. J., Hood, A. W., et al. 2020, A&A, 635, A168
  * Keppens et al. (2021) Keppens, R., Teunissen, J., Xia, C., & Porth, O. 2021, Computers & Mathematics With Applications, 81, 316
  * Keppens et al. (1999) Keppens, R., Tóth, G., Westermann, R. H. J., & Goedbloed, J. P. 1999, Journal of Plasma Physics, 61, 1
  * Kontar et al. (2017) Kontar, E. P., Perez, J. E., Harra, L. K., et al. 2017, Phys. Rev. Lett., 118, 155101
  * Kontar et al. (2011) Kontar, E. P., Brown, J. C., Emslie, A. G., et al. 2011, Space Sci. Rev., 159, 301
  * Krucker et al. (2008) Krucker, S., Battaglia, M., Cargill, P. J., et al. 2008, A&A Rev., 16, 155
  * Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17
  * Masuda et al. (1994) Masuda, S., Kosugi, T., Hara, H., Tsuneta, S., & Ogawara, Y. 1994, Nature, 371, 495
  * Meshkov (1969) Meshkov, E. 1969, Fluid Dynamics, 4, 101
  * Miller et al. (1997) Miller, J. A., Cargill, P. J., Emslie, A. G., et al. 1997, J. Geophys. Res., 102, 14631
  * Parker (1963) Parker, E. N. 1963, ApJS, 8, 177
  * Pinto et al. (2015) Pinto, R. F., Vilmer, N., & Brun, A. S. 2015, A&A, 576, A37
  * Reep & Russell (2016) Reep, J. W., & Russell, A. J. B. 2016, ApJ, 818, L20
  * Richtmyer (1960) Richtmyer, R. D. 1960, Communications on Pure and Applied Mathematics, 13, 297
  * Ruan et al. (2018) Ruan, W., Xia, C., & Keppens, R. 2018, A&A, 618, A135
  * Ruan et al. (2020) —. 2020, ApJ, 896, 97
  * Ruan et al. (2021) Ruan, W., Zhou, Y., & Keppens, R. 2021, ApJ, 920, L15
  * Russell & Fletcher (2013) Russell, A. J. B., & Fletcher, L. 2013, ApJ, 765, 81
  * Ruuth (2006) Ruuth, S. 2006, Mathematics of Computation, 75, 183
  * Shen et al. (2022) Shen, C., Chen, B., Reeves, K. K., et al. 2022, Nature Astronomy, 6, 317
  * Shi et al. (2021) Shi, M., Van Doorsselaere, T., Antolin, P., & Li, B. 2021, ApJ, 922, 60
  * Shibata et al. (1995) Shibata, K., Masuda, S., Shimojo, M., et al. 1995, ApJ, 451, L83
  * Stores et al. (2021) Stores, M., Jeffrey, N. L. S., & Kontar, E. P. 2021, ApJ, 923, 40
  * Su et al. (2013) Su, Y., Veronig, A. M., Holman, G. D., et al. 2013, Nature Physics, 9, 489
  * Tomczak (2001) Tomczak, M. 2001, A&A, 366, 294
  * Tomczak & Ciborski (2007) Tomczak, M., & Ciborski, T. 2007, A&A, 461, 315
  * van Leer (1974) van Leer, B. 1974, Journal of Computational Physics, 14, 361
  * Čada & Torrilhon (2009) Čada, M., & Torrilhon, M. 2009, Journal of Computational Physics, 228, 4118
  * Wang et al. (2022) Wang, Y., Cheng, X., Ren, Z., & Ding, M. 2022, ApJ, 931, L32
  * Xia et al. (2018) Xia, C., Teunissen, J., El Mellah, I., Chané, E., & Keppens, R. 2018, ApJS, 234, 30
  * Yamada et al. (2010) Yamada, M., Kulsrud, R., & Ji, H. 2010, Reviews of modern physics, 82, 603
  * Ye et al. (2021) Ye, J., Cai, Q., Shen, C., et al. 2021, ApJ, 909, 45
  * Ye et al. (2019) Ye, J., Shen, C., Raymond, J. C., Lin, J., & Ziegler, U. 2019, MNRAS, 482, 588
  * Zhou et al. (2021) Zhou, Y., Williams, R. J., Ramaprabhu, P., et al. 2021, Physica D: Nonlinear Phenomena, 423, 132838
  * Zhou et al. (2021) Zhou, Y.-H., Ruan, W.-Z., Xia, C., & Keppens, R. 2021, A&A, 648, A29


# Unveiling Nucleon 3D Chiral-Odd Structure with Jet Axes

Wai Kin Lai<EMAIL_ADDRESS>Guangdong Provincial Key Laboratory of
Nuclear Science, Institute of Quantum Matter, South China Normal University,
Guangzhou 510006, China  Guangdong-Hong Kong Joint Laboratory of Quantum
Matter, Southern Nuclear Science Computing Center, South China Normal
University, Guangzhou 510006, China  Department of Physics and Astronomy,
University of California, Los Angeles, CA 90095, USA  Xiaohui Liu
<EMAIL_ADDRESS>Center of Advanced Quantum Studies, Department of Physics,
Beijing Normal University, Beijing 100875, China  Center for High Energy
Physics, Peking University, Beijing 100871, China  Manman Wang  Center of
Advanced Quantum Studies, Department of Physics, Beijing Normal University,
Beijing 100875, China  Hongxi Xing<EMAIL_ADDRESS>Guangdong Provincial
Key Laboratory of Nuclear Science, Institute of Quantum Matter, South China
Normal University, Guangzhou 510006, China  Guangdong-Hong Kong Joint
Laboratory of Quantum Matter, Southern Nuclear Science Computing Center, South
China Normal University, Guangzhou 510006, China

(August 27, 2024)

###### 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.

##  I 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 Kang  _et
al._ (2017); Liu  _et al._ (2019); Gutierrez-Reyes  _et al._ (2019a); Arratia
_et al._ (2020); Liu  _et al._ (2020); Gutierrez-Reyes  _et al._ (2018,
2019b); Kang  _et al._ (2022, 2020a). 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) Anderle  _et
al._ (2021) or detectors more optimized for low energy scales, such as the EIC
Comprehensive Chromodynamics Experiment (ECCE) Abdul Khalek  _et al._ (2021).
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 Liu (2021).

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 Liu and Xing (2021). 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 II, 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 Accardi  _et al._ (2022). In
Section III, we briefly review the notion of the T-odd jet function. In
Section IV, 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 Abdul Khalek  _et al._ (2021), and HERA. In Section V, 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 VI, we
give a summary and an outlook.

##  II Measuring photon $q_{T}$ in DIS

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. 1. Therefore, the jet-axis probe is
fully differential just like the SIDIS.

Figure 1: 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}$.

Based on what we have described, we can derive the factorization formula for
the jet-axis probe. Formally, the factorization theorem reads

$\displaystyle 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})\,,$ (1)

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 Cacciari  _et al._ (2012), defined by

$\displaystyle d_{ij}={\rm
min}(E^{-2}_{i},E^{-2}_{j})\frac{1-\cos\theta_{ij}}{1-\cos R}\,,\quad
d_{iB}=E_{i}^{-2}\,.$ (2)

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$.

##  III T-odd jet function

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. (1)
derived from the soft-collinear effective theory (SCET) Bauer  _et al._
(2001); Bauer and Stewart (2001); Bauer  _et al._ (2002a, b). The
factorization theorem involves the transverse-momentum-dependent (TMD)
correlator

$\displaystyle{\cal J}^{ij}(z,k_{T})$
$\displaystyle=\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}\,,$ (3)

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. (3), which determines the $P_{J}$ and hence the
jet axis, and can be calculated perturbatively.

Conventionally, only the chiral-even Dirac structure $\not{\bar{n}}$ in Eq.
(3) was considered. However, as noted in Ref. Liu and Xing (2021), 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. (3)
in general is a sum of two structures:

$\displaystyle{\cal J}(z,k_{T})$ $\displaystyle={\cal
J}_{1}(z,k_{T}^{2})\frac{\not{\bar{n}}}{2}+i{\cal
J}_{T}(z,k_{T}^{2})\frac{\not{k}_{T}\not{\bar{n}}}{2}\,,$ (4)

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:

  * •

Universality
Like the traditional jet function, the T-odd jet function is process
independent.

  * •

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.

  * •

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 Gutierrez-Reyes  _et al._ (2019b):

$\displaystyle{\cal J}(z,k_{T},R)=\delta(1-z){\mathfrak{J}}(k_{T})+{\cal
O}\left(\frac{k_{T}^{2}}{E_{J}^{2}R^{2}}\right)\,.$ (5)

  * •

Nonperturbative predictability
Similar to the study in Ref. Becher and Bell (2014), 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
Vladimirov (2020), and as a vacuum matrix element, it can be calculated on the
lattice Shanahan  _et al._ (2020); Zhang  _et al._ (2020). 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.

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.

##  IV Photon $q_{T}$ with respect to the jet axis in deep-inelastic
scattering

Consider deep-inelastic scattering of an electron off a polarized nucleon
$e^{-}(l)+N(P)\to e^{-}(l^{\prime})+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. 2 (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. 2 (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 Bacchetta  _et
al._ (2007).

Figure 2: 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. Figure 3: Kinematic configuration of DIS in
the GNS system.

Let $M$ be the mass of the nucleon $N$ and $q=l-l^{\prime}$ is the momentum
carried by the virtual photon with virtuality $Q^{2}=-q^{2}$. We introduce the
invariant variables

$\displaystyle 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}\,.$ (6)

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

$\displaystyle g_{\perp}^{\mu\nu}$
$\displaystyle=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)\,.$
(7)

The 3-momenta ${\bm{l}}$ and ${\bm{l}}^{\prime}$ 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

$\displaystyle\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}}}\,,$
(8)

where

$\displaystyle\epsilon_{\perp}^{\mu\nu}$
$\displaystyle=\epsilon^{\mu\nu\rho\sigma}\frac{P_{\rho}q_{\sigma}}{P\cdot
q\sqrt{1+\gamma^{2}}}\,,$ (9)

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. 3. The nucleon spin is decomposed as sum of a longitudinal
and a perpendicular component,

$\displaystyle 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}}}\,.$ (10)

The helicity of the incoming electron is denoted by $\lambda_{e}$. We define
$\psi$ as the azimuthal angle of $\bm{l}^{\prime}$ around $\bm{l}$. The fully
differential cross section has the most general form given by

$\displaystyle\frac{d\sigma}{dxdyd\psi
dzd\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.$
$\displaystyle\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}$
$\displaystyle\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]$
$\displaystyle\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]$
$\displaystyle\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.$
$\displaystyle\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})}$
$\displaystyle\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]$
$\displaystyle\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.$
$\displaystyle\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\\}\,,$
(11)

where $\alpha$ is the fine structure constant and

$\displaystyle\epsilon=\frac{1-y-\frac{1}{4}\gamma^{2}y^{2}}{1-y+\frac{1}{2}y^{2}+\frac{1}{4}\gamma^{2}y^{2}}\,.$
(12)

The structure functions $F$’s on the right-hand side of Eq. (11) 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. (11):

$\displaystyle\frac{y^{2}}{2(1-\epsilon)}\approx 1-y+\frac{1}{2}y^{2},$
$\displaystyle\frac{y^{2}}{2(1-\epsilon)}\epsilon\approx 1-y,$
$\displaystyle\frac{y^{2}}{2(1-\epsilon)}\sqrt{2\epsilon(1+\epsilon)}\approx(2-y)\sqrt{1-y},$
$\displaystyle\frac{y^{2}}{2(1-\epsilon)}\sqrt{2\epsilon(1-\epsilon)}\approx
y\sqrt{1-y}\,,$
$\displaystyle\frac{y^{2}}{2(1-\epsilon)}\sqrt{1-\epsilon^{2}}\approx
y\left(1-\frac{1}{2}y\right)\,.$ (13)

The structure functions $F$’s are convolutions of the nucleon TMD PDFs and the
jet functions. As noted in Eq. (4), 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 Angeles-Martinez  _et al._ (2015), each encoding a
specific correlation between the quark spin and the proton spin as depicted in
Table 1. 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

$\displaystyle{\cal C}[wf\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})\,,$
(14)

where $a$ denotes a quark or antiquark flavor. At leading order in
$\alpha_{s}$ and $M/Q$, the nonvanishing $F$’s are given by

$\displaystyle F_{UU,T}={\cal C}[f_{1}{\cal J}_{1}]\,,$ (15) $\displaystyle
F_{LL}={\cal C}[g_{1L}{\cal J}_{1}]\,,$ (16) $\displaystyle
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]\,,$ (17) $\displaystyle 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]\,,$
(18) $\displaystyle 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]\,,$ (19) $\displaystyle 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]\,,$ (20) $\displaystyle F^{\sin(\phi_{J}+\phi_{S})}_{UT}={\cal
C}\left[-\hat{\bm{h}}\cdot\bm{k}_{T}h_{1}{\cal J}_{T}\right]\,,$ (21)
$\displaystyle 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]\,,$ (22)

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.
(15)-(22), 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.

hadron quark | unpolarized | chiral | transverse
---|---|---|---
$U$ | $f_{1}$ |  | $h_{1}^{\perp}$
$L$ |  | $g_{1L}$ | $h^{\perp}_{1L}$
$T$ | $f^{\perp}_{1T}$ | $g_{1T}$ | $h_{1T},h^{\perp}_{1T}$
Table 1: The eight TMD PDFs of a nucleon at leading twist.

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 III, 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. (2) with $R=1$ and the WTA scheme, so
that

$\displaystyle{\cal J}(z,k_{T},R)$
$\displaystyle=\delta(1-z){\mathfrak{J}}(k_{T})+{\cal
O}\left(\frac{k_{T}^{2}}{E_{J}^{2}R^{2}}\right)\,,$ (23)

where

$\displaystyle{\mathfrak{J}}(k_{T})$
$\displaystyle={J}(k_{T}^{2})\frac{\not{\bar{n}}}{2}+i{J}_{T}(k_{T}^{2})\frac{\not{k}_{T}\not{\bar{n}}}{2}\,.$
(24)

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. (11). 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

$\displaystyle A^{\sin(\phi_{J}+\phi_{S})}(|{\bm{P}}_{J\perp}|)$
$\displaystyle=\frac{2}{|\bm{S}_{\perp}|\int d\sigma\epsilon}\int
d\sigma\sin(\phi_{J}+\phi_{S})$ (25) $\displaystyle=\frac{\langle\epsilon
F_{UT}^{\sin(\phi_{J}+\phi_{S})}\rangle}{\bar{\epsilon}\langle
F_{UU,T}\rangle}\,,$ (26)

where by $\langle X\rangle$ we mean

$\displaystyle\langle X\rangle$ $\displaystyle=\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\,,$ (27)

and

$\displaystyle\bar{\epsilon}=\frac{\int d\sigma\epsilon}{\int d\sigma}\,.$
(28)

We can write $F_{UU,T}$ and $F_{UT}^{\sin(\phi_{J}+\phi_{S})}$ as

$\displaystyle F_{UU,T}$
$\displaystyle=x\sum_{a}e_{a}^{2}\int\frac{d^{2}b}{(2\pi)^{2}}e^{-i{\bm{P}}_{J\perp}\cdot{\bm{b}}}\tilde{f}^{a}_{1}(x,b^{2})\tilde{J}^{a}(b^{2})\,,$
(29) $\displaystyle F_{UT}^{\sin(\phi_{J}+\phi_{S})}$
$\displaystyle=-ix\sum_{a}e_{a}^{2}\int\frac{d^{2}b}{(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})\,,$
(30)

where we have used the Fourier transforms $\tilde{f}_{1}(x,b^{2})=\int
d^{2}p_{T}\,e^{-i{\bm{p}}_{T}\cdot{\bm{b}}}f_{1}(x,p_{T}^{2})$ and
$\tilde{J}(b^{2})=\int
d^{2}k_{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. Kang
_et al._ (2015), we include the effect of evolution by including a Sudakov
factor in $b$-space

$\displaystyle\tilde{f}_{1}(x,b^{2},Q)=e^{-S_{\rm pert}-S_{\rm
NP}}\tilde{f}_{1}(x,b^{2},Q_{0})\,,$ (31)
$\displaystyle\tilde{J}(b^{2},Q)=e^{-S_{\rm pert}-S_{\rm
NP}}\tilde{J}(b^{2},Q_{0})\,,$ (32)

and similarly for $\tilde{h}_{1}$ and $\tilde{J}_{T}$, where

$\displaystyle S_{\rm pert}$
$\displaystyle=\int^{Q}_{\mu_{b}}\frac{d\mu}{\mu}\,\frac{\alpha_{s}(\mu)}{2\pi}C_{F}\ln\frac{Q^{2}}{\mu^{2}}\,,$
(33) $\displaystyle S_{\rm NP}$
$\displaystyle=\frac{g_{2}}{2}\ln\left(\frac{b}{b_{*}}\right)\ln\left(\frac{Q}{Q_{0}}\right)\,,$
(34)

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. Martin  _et
al._ (2015). 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. Kang  _et al._ (2020b) 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. de
Florian  _et al._ (2007). 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. de Florian  _et al._ (2007). Figure 4
(a) shows the predictions of the asymmetry
$A^{\sin(\phi_{J}+\phi_{S})}(|{\bm{P}}_{J\perp}|)$ at the EIC according to Eq.
(26) (solid lines) and Eq. (25) (data points) from simulations using Pythia
8.2 Sjöstrand  _et al._ (2015) with the package StringSpinner Kerbizi and
Lönnblad (2021), which incorporates spin interactions in the event generator.
From Fig. 4 (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. (11) roughly agree with the event generator
simulations. In Figure 4 (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 Pythia 8.2 Sjöstrand  _et al._ (2015) with
StringSpinner, with the same kinematic setting as Fig. 4 (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.

Figure 4: 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. (26).

Likewise, we can make predictions on the asymmetry that probes the Boer-
Mulders function defined by

$\displaystyle A^{\cos(2\phi_{J})}(|{\bm{P}}_{J\perp}|)$
$\displaystyle=\frac{2}{\int d\sigma\epsilon}\int d\sigma\cos(2\phi_{J})$ (35)
$\displaystyle=\frac{\langle\epsilon
F_{UU}^{\cos(2\phi_{J})}\rangle}{\bar{\epsilon}\langle F_{UU,T}\rangle}\,.$
(36)

The structure function $F_{UU}^{\cos(2\phi_{J})}$ can be written as

$\displaystyle F_{UU}^{\cos(2\phi_{J})}$
$\displaystyle=-\frac{x}{M}\sum_{a}e_{a}^{2}\int\frac{d^{2}b}{(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]\,,$ (37)

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. Barone  _et al._ (2010). The predictions on
$A^{\cos(2\phi_{J})}(|{\bm{P}}_{J\perp}|)$ at the EIC according to Eq. (36)
are shown in Fig. (5).

Figure 5: Azimuthal asymmetry $A^{\cos(2\phi_{J})}$ in the jet-axis probe at
the EIC as predicted from Eq. (36).

The predictions of $A^{\sin(\phi_{J}+\phi_{S})}$ and $A^{\cos(2\phi_{J})}$ at
the EicC are shown in Fig. 6. As in Fig. 4 (a) and Fig. 5, the data points
with error bars are from event generator simulations and the lines are from
the factorization formulae.

Figure 6: 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.

For the sake of comparison with data in SIDIS, in Fig. 7 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. (26) and (36) (dashed lines), predictions for pion production
from the parallels of Eqs. (26) and (36) as in Refs. Bacchetta  _et al._
(2007); Barone  _et al._ (2010) (solid lines), and data points for pion
production from the HERMES experiment Airapetian  _et al._ (2010, 2013) (data
points with error bars). From Fig. 7, 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.

Figure 7: 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.

##  V back-to-back dijet production in $e^{+}e^{-}$ annihilation

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. 8. 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$ Kang  _et al._ (2015) is given by

$\displaystyle A$ $\displaystyle=2\int
d\cos\theta\,\frac{d\phi_{1}}{\pi}\cos(2\phi_{1})A^{J_{1}J_{2}}\,,$ (38)

where

$\displaystyle A^{J_{1}J_{2}}$
$\displaystyle=1+\cos(2\phi_{1})\frac{\sin^{2}\theta}{1+\cos^{2}\theta}\frac{F_{T}}{F_{U}}\,,$
(39)

with

$\displaystyle F_{U}$
$\displaystyle=|\bm{q}_{T}|\,\sum_{q}e_{q}^{2}\,\int\frac{{d}^{2}b}{(2\pi)^{2}}e^{i\bm{q}_{T}\cdot\bm{b}}\tilde{J}_{1}^{q}(b^{2}){\tilde{J}}^{\bar{q}}_{1}(b^{2})\,,$
(40) $\displaystyle F_{T}$
$\displaystyle=|\bm{q}_{T}|\,\sum_{q}\,e_{q}^{2}\,\int\frac{{d}^{2}b}{(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})\,.$
(41)

Figure 8: Back-to-back dijet production in $e^{+}e^{-}$ annihilation.

In Fig. 9, we plot the asymmetry $A$ as a function of $|\bm{q}_{T}|$ as
predicted by Eq. (39) 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.

Figure 9: Azimuthal asymmetry $A$ for dijet production in $e^{+}e^{-}$
annihilation as predicted by Eq. (39) at Belle, LEP and the CEPC, with
$Q_{J}>0.25$ for one of the jets and $Q_{J}<-0.25$ for the other.

##  VI Summary and outlook

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.

###### Acknowledgements.

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.

## References

  * Kang  _et al._ (2017) Z.-B. Kang, X. Liu, F. Ringer, and H. Xing, JHEP 11, 068 (2017), arXiv:1705.08443 [hep-ph] .
  * Liu  _et al._ (2019) X. Liu, F. Ringer, W. Vogelsang, and F. Yuan, Phys. Rev. Lett. 122, 192003 (2019), arXiv:1812.08077 [hep-ph] .
  * Gutierrez-Reyes  _et al._ (2019a) D. Gutierrez-Reyes, Y. Makris, V. Vaidya, I. Scimemi, and L. Zoppi, JHEP 08, 161 (2019a), arXiv:1907.05896 [hep-ph] .
  * Arratia  _et al._ (2020) M. Arratia, Z.-B. Kang, A. Prokudin, and F. Ringer, Phys. Rev. D 102, 074015 (2020), arXiv:2007.07281 [hep-ph] .
  * Liu  _et al._ (2020) X. Liu, F. Ringer, W. Vogelsang, and F. Yuan, Phys. Rev. D 102, 094022 (2020), arXiv:2007.12866 [hep-ph] .
  * Gutierrez-Reyes  _et al._ (2018) D. Gutierrez-Reyes, I. Scimemi, W. J. Waalewijn, and L. Zoppi, Phys. Rev. Lett. 121, 162001 (2018), arXiv:1807.07573 [hep-ph] .
  * Gutierrez-Reyes  _et al._ (2019b) D. Gutierrez-Reyes, I. Scimemi, W. J. Waalewijn, and L. Zoppi, JHEP 10, 031 (2019b), arXiv:1904.04259 [hep-ph] .
  * Kang  _et al._ (2022) Z.-B. Kang, K. Lee, D. Y. Shao, and F. Zhao, (2022), arXiv:2201.04582 [hep-ph] .
  * Kang  _et al._ (2020a) Z.-B. Kang, K. Lee, and F. Zhao, Phys. Lett. B 809, 135756 (2020a), arXiv:2005.02398 [hep-ph] .
  * Anderle  _et al._ (2021) D. P. Anderle _et al._ , Front. Phys. (Beijing) 16, 64701 (2021), arXiv:2102.09222 [nucl-ex] .
  * Abdul Khalek  _et al._ (2021) R. Abdul Khalek _et al._ , (2021), arXiv:2103.05419 [physics.ins-det] .
  * Liu (2021) X. Liu (Jets at EicC: jet axes for TMDs presented at EicC CDR meeting 2021, (online), 20-21 November, 2021).
  * Liu and Xing (2021) X. Liu and H. Xing, (2021), arXiv:2104.03328 [hep-ph] .
  * Accardi  _et al._ (2022) A. Accardi _et al._ , in  _2022 Snowmass Summer Study_ (2022) arXiv:2204.02280 [hep-ex] .
  * Cacciari  _et al._ (2012) M. Cacciari, G. P. Salam, and G. Soyez, Eur. Phys. J. C 72, 1896 (2012), arXiv:1111.6097 [hep-ph] .
  * Bauer  _et al._ (2001) C. W. Bauer, S. Fleming, D. Pirjol, and I. W. Stewart, Phys. Rev. D63, 114020 (2001), arXiv:hep-ph/0011336 [hep-ph] .
  * Bauer and Stewart (2001) C. W. Bauer and I. W. Stewart, Phys. Lett. B516, 134 (2001), arXiv:hep-ph/0107001 [hep-ph] .
  * Bauer  _et al._ (2002a) C. W. Bauer, D. Pirjol, and I. W. Stewart, Phys. Rev. D65, 054022 (2002a), arXiv:hep-ph/0109045 [hep-ph] .
  * Bauer  _et al._ (2002b) C. W. Bauer, S. Fleming, D. Pirjol, I. Z. Rothstein, and I. W. Stewart, Phys. Rev. D66, 014017 (2002b), arXiv:hep-ph/0202088 [hep-ph] .
  * Becher and Bell (2014) T. Becher and G. Bell, Phys. Rev. Lett. 112, 182002 (2014), arXiv:1312.5327 [hep-ph] .
  * Vladimirov (2020) A. A. Vladimirov, Phys. Rev. Lett. 125, 192002 (2020), arXiv:2003.02288 [hep-ph] .
  * Shanahan  _et al._ (2020) P. Shanahan, M. Wagman, and Y. Zhao, Phys. Rev. D 102, 014511 (2020), arXiv:2003.06063 [hep-lat] .
  * Zhang  _et al._ (2020) Q.-A. Zhang _et al._ (Lattice Parton), Phys. Rev. Lett. 125, 192001 (2020), arXiv:2005.14572 [hep-lat] .
  * Bacchetta  _et al._ (2007) A. Bacchetta, M. Diehl, K. Goeke, A. Metz, P. J. Mulders, and M. Schlegel, JHEP 02, 093 (2007), arXiv:hep-ph/0611265 .
  * Angeles-Martinez  _et al._ (2015) R. Angeles-Martinez _et al._ , Acta Phys. Polon. B 46, 2501 (2015), arXiv:1507.05267 [hep-ph] .
  * Kang  _et al._ (2015) Z.-B. Kang, A. Prokudin, P. Sun, and F. Yuan, Phys. Rev. D 91, 071501 (2015), arXiv:1410.4877 [hep-ph] .
  * Martin  _et al._ (2015) A. Martin, F. Bradamante, and V. Barone, Phys. Rev. D 91, 014034 (2015), arXiv:1412.5946 [hep-ph] .
  * Kang  _et al._ (2020b) Z.-B. Kang, X. Liu, S. Mantry, and D. Y. Shao, Phys. Rev. Lett. 125, 242003 (2020b), arXiv:2008.00655 [hep-ph] .
  * de Florian  _et al._ (2007) D. de Florian, R. Sassot, and M. Stratmann, Phys. Rev. D 75, 114010 (2007), arXiv:hep-ph/0703242 .
  * Sjöstrand  _et al._ (2015) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, Comput. Phys. Commun. 191, 159 (2015), arXiv:1410.3012 [hep-ph] .
  * Kerbizi and Lönnblad (2021) A. Kerbizi and L. Lönnblad, (2021), arXiv:2105.09730 [hep-ph] .
  * Barone  _et al._ (2010) V. Barone, S. Melis, and A. Prokudin, Phys. Rev. D 81, 114026 (2010), arXiv:0912.5194 [hep-ph] .
  * Airapetian  _et al._ (2010) A. Airapetian _et al._ (HERMES), Phys. Lett. B 693, 11 (2010), arXiv:1006.4221 [hep-ex] .
  * Airapetian  _et al._ (2013) A. Airapetian _et al._ (HERMES), Phys. Rev. D 87, 012010 (2013), arXiv:1204.4161 [hep-ex] .