††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Active control of a diode laser with injection locking

Ziting Chen  Department of Physics, The Hong Kong University of Science and
Technology, Clear Water Bay, Kowloon, Hong Kong, China  Bojeong Seo
Department of Physics, The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Mingchen Huang  Department of
Physics, The Hong Kong University of Science and Technology, Clear Water Bay,
Kowloon, Hong Kong, China  Mithilesh K. Parit  Department of Physics, The Hong
Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong
Kong, China  Peng Chen<EMAIL_ADDRESS>Department of Physics, The Hong Kong
University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong,
China  Gyu-Boong Jo<EMAIL_ADDRESS>Department of Physics, The Hong Kong
University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong,
China

###### Abstract

We present a simple and effective method to implement an active stabilization
of a diode laser with injection locking, which requires minimal user
intervenes. The injection locked state of the diode laser is probed by a
photodetector, of which sensitivity is enhanced by a narrow laser-line filter.
Taking advantage of the characteristic response of laser power to spectral
modes from the narrow laser-line filter, we demonstrate that high spectral
purity and low intensity noise of the diode can be simultaneously maintained
by an active feedback to the injected laser. Our method is intrinsically cost-
effective, and does not require bulky devices, such as Fabry-Perot
interferometers or wavemeters, to actively stabilize the diode laser. Based on
successful implementation of this method in our quantum gas experiments, it is
conceivable that our active stabilization will greatly simplify potential
applications of injection locking of diode lasers in modularized or integrated
optical systems.

##  I Introduction

Injection locking of a diode laser is essential in a majority of fields, such
as laser spectroscopy Wieman and Hollberg (1991); Eng  _et al._ (1980), laser
cooling and trapping Schäfer  _et al._ (2015); Komori  _et al._ (2003); Hosoya
_et al._ (2015); Shimada  _et al._ (2013); Pagett  _et al._ (2016); Schkolnik
_et al._ (2020), optical communication Goldberg  _et al._ (1983); Pang  _et
al._ (2020); Paraïso  _et al._ (2019) and high-precision metrology Diddams
_et al._ (2004); Takamoto  _et al._ (2005); Hinkley  _et al._ (2013); Liu  _et
al._ (2015). It allows phase and frequency locking of a "slave" laser in
reference to a "master" laser by an optical link Hadley (1986). Among the
advantages of using injection locking of the diode are its cost-effectiveness,
flexibility and efficiency. For example, it offers a simple solution to
provide a clean laser source without a special design of laser diodes or
controllers. The spectral mode of injection locked slave laser is nearly
dependent on the master laser. Most of single mode diode lasers, including
some multimode diode lasers can be used for injection locking. This feature
allows us to harness injection locking for achieving a relatively high power
single-frequency diode laser source, without using an optical grating that
causes significant optical loss. Beside those advantages, the injection
locking scheme can also be exploited as active optical filters in quantum
optics Tistomo and Gee (2011).

In spite of its broad applications, it remains a challenge that an injection
locked laser requires mechanically and thermally isolated environment
otherwise necessitates an unpredictable re-locking by an user. A major
obstacle is the stability of injection locking which is hard to actively
maintain. Environmental variations and current drifts of laser controller can
frequently cause the slave laser out of lock. It is not convenient to manually
re-lock the slave laser, especially for injection locking of multiple laser
sources. In order to actively stabilize the injection locking, it is proposed
to use a Fabry-Perot interferometer (FPI) to monitor a locking state of the
slave laser Saxberg  _et al._ (2016). The FPI operates in a scanning mode to
detect the spectral mode of the slave laser and thereby determines its locking
state. As the FPI itself has frequency drifts, the scanning mode is necessary
for proper detection in their method. However, a FPI is sensitive to optical
beam alignment and the scanning mode significantly sacrifices the detection
bandwidth, making it unfavorable for general-purpose applications.

Figure 1: Measured laser power as a function of diode current. (A).
Transmitted laser power from a narrow laser-line filter shows a series of
spikes (black solid line). At those spikes, the injected slave laser operates
exactly on the same frequency of the master laser. (B). Reflected laser power
from the narrow laser-line filter exhibits a series of dips (red solid line),
which will be used for active locking.  Figure 2: Schematics of the active
stabilization scheme. A narrow laser-line filter is inserted into the beam
path of injection locking, from which the reflected beam is collected by a
photodetector (PD1) and sent to a PID controller. The PID controller generates
feedback signal to regulate the slave laser current. Optionally, the
transmitted power, which will be useful for applications, can be measured with
another photodetector (PD2).

In this work, we propose a simple way to effectively implement active control
of laser injection locking without requiring additional bulky devices. We
employ a photodetector and a narrow laser-line filter to probe the spectral
mode of slave laser, which can be easily integrated into the beam path. The
spectral mode of slave laser is monitored through the reflected signal from
the narrow-laser line filter, which is neither sensitive to beam alignment nor
bandwidth limited. With a simple feedback to the slave laser current, we
demonstrate that high spectral mode purity and low intensity noise of the
slave laser can be simultaneously obtained. We also show that the slave laser
can be locked even when the master laser frequency is dynamically scanned over
GHz range. We have successfully applied this method in laser cooling and
trapping of erbium atoms towards degenerate dipolar quantum gases Lu  _et al._
(2011); Aikawa  _et al._ (2012). Our method paves the way towards general
applications of injection locking. A promising scenario is coherent addition
of multiple laser diodes for light power scaling, where injection locking
plays a fundamental role Schimmel  _et al._ (2017).

##  II Active control scheme for injection locking of a diode laser

###  II.1 Working principle

A diode laser with light injection assumes power spikes versus laser current
above threshold. These power spikes happen when the internal cavity of diode
laser is near resonant to the frequency of injection light. The current
threshold for lasing decreases as coherent addition of light power from the
external injection beam. Therefore, the total output laser power increases in
the injection locking area. As the diode laser frequency is tunable with
current, the injection locking is usually done by adjusting the current.

Fig. 1(A) shows a typical output power-current curve from an injected laser
diode (transmitted from a laser-line filter), where output power is enhanced
only at certain currents. Those power spikes reflect that the internal cavity
mode of the diode is resonant to the injection light. The shape of spikes are
asymmetrical due to the nonlinear response of lasing to the external injection
light in contrast to the passive cavities, where the transmission follows a
Lorenz curve. As the power spikes are correlated to injection locking, we can
use them to monitor locking state of the diode laser. This can be done with a
simple photodetector to monitor the laser power. The detection sensitivity can
be further enhanced by a narrow laser-line filter Schkolnik  _et al._ (2020),
which typically has a transmission bandwidth of 2-3 nm and high optical
density (OD$>$5) for side bands rejection. The central transmission rate of
laser-line filters can be also very high ($>$90%). The photodetector measures
the reflection beam of the diode laser from the laser-line filter, where a
local minimum power indicates an injection locking point, conditioning that
injection light frequency is within the transmission bandwidth. With this
detection geometry, the DC offset of the photodetector can be largely removed,
resulting in an enhanced signal-noise ratio compared to monitoring the
transmitted laser power, as shown in Fig. 1(B).

###  II.2 Schematics of the active stabilization scheme

A schematic diagram of the active injection locking of the slave laser diode
is shown in Fig. 2. An external cavity diode laser (ECDL, Toptica DL Pro 401
nm), a master laser, provides 0.7 mW light for injecting a slave laser
(Nichia, NDV4316, 400-410 nm) through an optical isolator (Thorlabs
IO-3D-405-PBS). The optical isolator maintains bi-directional isolation for
master laser and slave laser beams to avoid optical crosstalk. A narrow laser-
line filter (Semrock LL01-405-25) is inserted in the beam path before the
optical isolator. This filter is slightly rotated in order to shift the
centeral frequency from 405 nm to 401 nm. The rotation angle is determined by
minimizing the reflective power of seeding light. A photodetector (Thorlabs
PDA36A2) is used to probe the light reflected from the laser-line filter. A
digital proportional-integral-derivative (PID) module (Toptica Digilock 110)
regulates the slave laser current based on the photodetector output.

As described above, a local minimum of the photodetector output corresponds to
an injection locking point of the slave laser. One may find an optimal
injection locking current by searching the desired output power of the slave
laser. However, we find it difficult to control the injected slave laser
through the peak-locking because of the asymmetric current-power response
including minimum plateaus. By noting that the locking state is less sensitive
to the current on the right side of the power dip of reflected laser beam, we
take advantage of the side-of-fringe locking for active control of the slave
laser current. (see Fig. 1 and Fig. 3 (A)). In contrast to the previous work
where a special algorithm is required to search for optimal current for
injection locking Saxberg  _et al._ (2016), our scheme does not require
constant searching mode set by an algorithm. Consequently, the active locking
scheme does not increase the intensity noise of the slave laser without
limiting the feedback bandwidth.

Figure 3: Emission spectrum of the injected slave laser under active
stabilization (A) We examine spectral purity of the slave laser at different
lock positions. Three positions $\alpha,\beta,\gamma$ are chosen in the dip of
the reflected laser power from the narrow laser-line filter. (B). FPI
transmission signals of the slave laser at locking positions
$\alpha,\beta,\gamma$ (C). Zoom-in profile of a single peak in curve $\alpha$
of (B).

###  II.3 Spectral purity with the active control

In Fig. 3, we examine how the side-of-fringe locking affects the spectral
purity for different locking position. With a slave laser operating at 25∘C
and 109 mA, we monitor the spectral mode of the injected laser with a Fabry-
Perot interferometry (FPI, Thorlabs SA210-3B, not shown in Fig. 2). In a
typical reflection signal (Fig. 3 (a)), we stabilize the injection locking on
the right side of the dip using the side-of-fringe locking scheme. At position
$\alpha$, closest to the minimum of the reflection signal, the FPI measurement
does not present any spurious modes, which suggests that the slave laser
frequency is faithfully locked to the master laser. At position $\beta$,
spurious modes slightly show up and become more severe at position $\gamma$,
which deviates most from the dip minimum among three positions. These results
demonstrate that injection locking can be maintained with a proper lock
position using the side-of-fringe locking. Our active locking scheme maintains
injection locking at the laser power slightly higher than the minimum of the
reflected laser power (e.g. the position $\alpha$), and stabilizes the power
fluctuation of the slave laser. To evaluate the stability of slave laser power
under the condition of injection locking with active control, we measure the
slave laser power transmitted from the laser-line filter. The ratio of the
root-mean-square (rms) value of laser power fluctuations to the mean value of
laser power is approximately 0.06% in the course of 6 minutes measurement.
Referring to the slave laser power stability without the active control, which
is about 0.08 to 0.1% (not shown in Fig. 4), we find that our active control
scheme can not only stabilize the injection locking but also reduce the laser
intensity noise.

Figure 4: Power stability of the injected slave laser. We monitor the emission
power of the slave laser for 6 minutes at the position of optional part (see
Fig. 2 ). Histogram of laser power is shown in the right panel.

###  II.4 Master laser frequency scan, locking bandwidth

Another advantage of this method is that the frequency of slave laser can be
dynamically tuned, following the master laser frequency. Fig. 5 shows the
wavelength (top) of the slave laser and the corresponding PID output voltage
(bottom) from the feedback controller when the master laser linearly scans
over 1 GHz. The change of the injection locking point induced by the master
laser is readily compensated by the active feedback in contrast to the passive
injection locking. Our results suggest that the demonstrated active feedback
method is useful for not only the frequency stabilization, but also
dynamically tracking the injection locking states. This feature would be
critical in applications where the slave laser frequency is also required to
vary, for example, coherent light detection and ranging (LIDAR) for range
detection Zhang  _et al._ (2019). The feedback bandwidth of our locking scheme
is mainly determined by the current controller (Thorlabs, LDC202C, 100 kHz) of
the slave laser. We have evaluated the feedback bandwidth by modulating the
current of the master laser.

Figure 5: Active feedback of injection locking with scanning of master laser
frequency. The wavelength of slave laser (top) is measured with scanning of
master laser frequency over 1.4 GHz. The corresponding PID controller output
voltage (bottom) shows that the active feedback compensates the shift of
injection locking point.

##  III Application in laser cooling of cold atoms

Recent developments in cooling non-alkali atoms for quantum simulation He  _et
al._ (2019) and computing Daley (2011) require stable high-power laser sources
at the wavelength of 400-500 nm. In early experiments, those wavelengths are
obtained by frequency doubling infrared lasers because of the limited power of
the blue laser diode. The active locking scheme demonstrated in this work will
provide a cost-effective way for cooling non-alkali cold atoms, such as
ytterbium Song  _et al._ (2016) and erbium Seo  _et al._ (2020) without using
such a frequency doubling setup.

The scheme is already being used in the daily production of Bose-Einstein
condensates of erbium atoms (168Er) in our lab Seo  _et al._ (2020). We employ
two slave lasers to provide sufficient power for slowing erbium atoms from an
effusive oven. One slave laser delivers its full power to the major Zeeman
slower, up to 30 mW after fiber coupling. The other slave laser provides
similar power, solely used for the second stage Zeeman slower Seo  _et al._
(2020). They are injection locked from the same master laser with the active
locking scheme demonstrated in this work. The master laser frequency is
externally locked to the erbium 401 nm transition ($4f^{12}6s^{2}\
{}^{3}H_{6}$ $\longrightarrow$ $4f^{12}(^{3}H_{6})6s6p(^{1}P_{1})$). We expect
our method to be applicable to general diode lasers for injection locking,
especially in ultraviolet wavelength, where high power light source is rare
and typically achieved by frequency doubling of infrared light Adams and
Ferguson (1992); Pizzocaro  _et al._ (2014).

##  IV Conclusion

In conclusion, we have demonstrated a novel scheme for active control of laser
injection locking. We have shown that the active stabilization of injection
locking can be easily implemented with simple optical elements but without a
special PID algorithm being required. The spectral purity and power noise of
the slave laser are evaluated. We demonstrate that the frequency of slave
laser can be dynamically tuned following the seeding light. All these features
ensure that our active control method offers a new possibility of harnessing
injection locking of diode lasers in modularized or integrated optical system.
Furthermore, we envision operating a diode laser with automatic injection
locking, which sets the stage for achieving an auto-lock laser system.

## Acknowledgments

P.C. acknowledges the support by Shanghai Natural Science Foundation (Grant
No. 18ZR1443800) and Innovation Promotion Association of the Chinese Academy
of Sciences. G.-B.J. acknowledges supports from the RGC, the Croucher
Foundation and the Hari Harilela foundation (project 16305317, 16304918,
16306119,16302420, C6005-17G and N-HKUST601/17).

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Wieman and Hollberg (1991) C. E. Wieman and L. Hollberg, Rev. Sci. Instrum. 62, 1 (1991).
  * Eng  _et al._ (1980) R. Eng, J. Butler, and K. Linden, Opt. Eng. 19, 196945 (1980).
  * Schäfer  _et al._ (2015) V. Schäfer, C. Ballance, C. Tock, and D. Lucas, Opt. Lett. 40, 4265 (2015).
  * Komori  _et al._ (2003) K. Komori, Y. Takasu, M. Kumakura, Y. Takahashi, and T. Yabuzaki, Jpn. journal applied physics 42, 5059 (2003).
  * Hosoya  _et al._ (2015) T. Hosoya, M. Miranda, R. Inoue, and M. Kozuma, Rev. Sci. Instrum. 86, 073110 (2015).
  * Shimada  _et al._ (2013) Y. Shimada, Y. Chida, N. Ohtsubo, T. Aoki, M. Takeuchi, T. Kuga, and Y. Torii, Rev. Sci. Instrum. 84, 063101 (2013).
  * Pagett  _et al._ (2016) C. Pagett, P. Moriya, R. Celistrino Teixeira, R. Shiozaki, M. Hemmerling, and P. W. Courteille, Rev. Sci. Instrum. 87, 053105 (2016).
  * Schkolnik  _et al._ (2020) V. Schkolnik, J. R. Williams, and N. Yu, arXiv preprint arXiv:2004.11732 (2020).
  * Goldberg  _et al._ (1983) L. Goldberg, H. Taylor, J. Weller, and D. Bloom, Electron. Lett. 19, 491 (1983).
  * Pang  _et al._ (2020) X.-L. Pang, A.-L. Yang, C.-N. Zhang, J.-P. Dou, H. Li, J. Gao, and X.-M. Jin, Phys. Rev. Appl. 13, 034008 (2020).
  * Paraïso  _et al._ (2019) T. K. Paraïso, I. De Marco, T. Roger, D. G. Marangon, J. F. Dynes, M. Lucamarini, Z. Yuan, and A. J. Shields, NPJ Quantum Inf. 5, 1 (2019).
  * Diddams  _et al._ (2004) S. A. Diddams, J. C. Bergquist, S. R. Jefferts, and C. W. Oates, Science 306, 1318 (2004).
  * Takamoto  _et al._ (2005) M. Takamoto, F.-L. Hong, R. Higashi, and H. Katori, Nature 435, 321 (2005).
  * Hinkley  _et al._ (2013) N. Hinkley, J. A. Sherman, N. B. Phillips, M. Schioppo, N. D. Lemke, K. Beloy, M. Pizzocaro, C. W. Oates, and A. D. Ludlow, Science 341, 1215 (2013).
  * Liu  _et al._ (2015) H. Liu, M. Yin, D. Kong, Q. Xu, S. Zhang, and H. Chang, Appl. Phys. Lett. 107, 151104 (2015).
  * Hadley (1986) G. Hadley, IEEE J. Quantum Electron. 22, 419 (1986).
  * Tistomo and Gee (2011) A. S. Tistomo and S. Gee, Opt. Express 19, 1081 (2011).
  * Saxberg  _et al._ (2016) B. Saxberg, B. Plotkin-Swing, and S. Gupta, Rev. Sci. Instrum. 87, 063109 (2016).
  * Lu  _et al._ (2011) M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Phys. Rev. Lett. 107, 190401 (2011).
  * Aikawa  _et al._ (2012) K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401 (2012).
  * Schimmel  _et al._ (2017) G. Schimmel, S. Janicot, M. Hanna, J. Decker, P. Crump, G. Erbert, U. Witte, M. Traub, P. Georges, and G. Lucas-Leclin, in  _High-Power Diode Laser Technology XV_ , Vol. 10086 (International Society for Optics and Photonics, 2017) p. 100860O.
  * Zhang  _et al._ (2019) X. Zhang, J. Pouls, and M. C. Wu, Opt. Express 27, 9965 (2019).
  * He  _et al._ (2019) C. He, E. Hajiyev, Z. Ren, B. Song, and G.-B. Jo, J. Phys. B: At. Mol. Opt. Phys. 52, 102001 (2019).
  * Daley (2011) A. J. Daley, Quantum Inf. Process. 10, 865 (2011).
  * Song  _et al._ (2016) B. Song, Y. Zou, S. Zhang, C.-w. Cho, and G.-B. Jo, Appl. Phys. B 122, 1 (2016).
  * Seo  _et al._ (2020) B. Seo, P. Chen, Z. Chen, W. Yuan, M. Huang, S. Du, and G.-B. Jo, Phys. Rev. A 102, 013319 (2020).
  * Adams and Ferguson (1992) C. Adams and A. Ferguson, Opt. Commun. 90, 89 (1992).
  * Pizzocaro  _et al._ (2014) M. Pizzocaro, D. Calonico, P. C. Pastor, J. Catani, G. A. Costanzo, F. Levi, and L. Lorini, Appl. optics 53, 3388 (2014).


# Towards three-dimensional nucleon structures at the Electron-Ion Collider in
China: A study of the Sivers function

Chunhua Zeng<EMAIL_ADDRESS>Institute of Modern Physics, Chinese
Academy of Sciences, Lanzhou, Gansu Province 730000, China  Lanzhou
University, Lanzhou, Gansu Province 730000, China  University of Chinese
Academy of Sciences, Beijing 100049, China  Tianbo Liu<EMAIL_ADDRESS>Key
Laboratory of Particle Physics and Particle Irradiation (MOE), Institute of
Frontier and Interdisciplinary Science, Shandong University, Qingdao, Shandong
266237, China  Peng Sun<EMAIL_ADDRESS>Institute of Modern Physics,
Chinese Academy of Sciences, Lanzhou, Gansu Province 730000, China  Department
of Physics and Institute of Theoretical Physics, Nanjing Normal University,
Nanjing, Jiangsu 210023, China  Yuxiang Zhao<EMAIL_ADDRESS>Institute of
Modern Physics, Chinese Academy of Sciences, Lanzhou, Gansu Province 730000,
China  University of Chinese Academy of Sciences, Beijing 100049, China  Key
Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle
Physics, Central China Normal University, Wuhan 430079, China

###### Abstract

Transverse momentum-dependent parton distribution functions (TMDs) provide
three-dimensional imaging of the nucleon in momentum space. With its
fundamental importance in understanding the spin structure of the nucleon, the
precise measurement of TMDs is considered one of the main physics topics of
the proposed Electron-Ion Collider in China (EicC). In this paper, we
investigate the impact of future semi-inclusive deep inelastic scattering
(SIDIS) data from EicC on the extraction of TMDs. Taking the Sivers function
as an example, we revisited the world SIDIS data, which serves as the input
for the simulated EicC data. By performing a combined analysis of the world
data and the EicC simulated data, we quantitatively demonstrate that the
Sivers functions can be precisely determined for various quark flavors,
especially for sea quarks, in the region, $x>0.005$, directly covered by the
EicC.

##  I Introduction

Nucleons are the basic building blocks of more than 98% of the visible matter
in the Universe. Understanding the nucleon structure in terms of the
underlying quarks and gluons degrees of freedom is a central issue in modern
nuclear and particle physics. Because of the nonperturbative property of
quantum chromodynamics (QCD) at hadronic energy scales, it is still a
challenging task to directly evaluate nucleon structures from first principles
[1], although much progress has been made in recent years [2, 3, 4]. The
nucleon spin structures have received great interest since the first
measurement of polarized lepton-nucleon deep inelastic scattering, which found
the quark spin only contributes a small fraction to the nucleon spin [5, 6]
and triggered the so-called proton spin puzzle [7, 8]. The remaining part of
the nucleon spin is nowadays attributed to the gluon angular momentum and the
quark orbital angular momentum, though one may have different decomposition
versions [9, 10, 11, 12]. After more than 30 years of efforts from both
experimental and theoretical aspects, the quark spin contribution is
relatively well determined [13, 14, 15] and the gluon spin part is also
started being known [16, 17, 18, 19]. However, we have very rare knowledge of
the orbital motion of quarks and gluons, which requires three-dimensional
imaging of the nucleon.

Since quarks and gluons cannot be directly observed in modern detectors, most
analyses of hadron involved high energy scatterings rely on the QCD
factorization [20, 21], which approximates the cross section by the
convolution of perturbatively calculable short-distance scattering off
partons, including quarks and gluons, and universal long-distance functions.
The measurement of partonic structures of the nucleon was pioneered by SLAC
[22] via the inclusive lepton-nucleon deep inelastic scattering (DIS) process,

$\displaystyle e(l)+N(P)\to e(l^{\prime})+X,$ (1)

where only the scattered lepton in the final state is detected and $X$ stands
for the undetected hadronic system. The labels $l$, $P$, and $l^{\prime}$ in
parentheses represent the four-momenta of corresponding particles. The large
momentum transfer $q=l-l^{\prime}$ with $Q^{2}=-q^{2}$ mediated by a virtual
photon (or some other gauge boson) provides a hard scale $Q\gg\Lambda_{\rm
QCD}$, which localizes the probe to “see” the partons in the nucleon. The
inclusive DIS cross section can be expressed as the convolution of the lepton-
parton hard scattering cross section and the collinear parton distribution
function (PDF) $f_{q}(x,\mu)$, which describes the density of finding a parton
of flavor $q$ carrying a fraction $x$ of the longitudinal momentum of the
parent nucleon with $\mu$ representing the factorization scale.

In semi-inclusive DIS (SIDIS) process,

$\displaystyle e(l)+N(P)\to e(l^{\prime})+h(P_{h})+X,$ (2)

a final-state hadron is identified in coincidence with the scattered lepton.
In addition to the hard scale $Q$, which localizes the probe as in the
inclusive DIS, the transverse momentum $P_{h\perp}$ of the hadron defined in
the virtual photon-nucleon frame provides an adjustable scale. When
$P_{h\perp}$ is integrated or compatible with $Q$, the SIDIS process
effectively becomes a single-scale problem and one can apply the collinear
factorization to approximate the cross section as the convolution of the hard
scattering cross section, the collinear PDF $f_{q}(x,\mu)$, and the collinear
fragmentation function $D_{q\to h}(z,\mu)$, which describes the probability
density of a parton of flavor $q$ fragmenting to the hadron carrying
longitudinal momentum fraction $z$. When $P_{h\perp}\ll Q$, the SIDIS is a
two-scale process, which is sensitive to the intrinsic transverse momenta of
partons, and one should apply the transverse momentum dependent (TMD)
factorization to express the cross section as the convolution of the hard
scattering cross section, the TMD PDF $f_{q}(x,k_{\perp},\mu)$, and the TMD FF
$D_{q\to h}(z,p_{\perp},\mu)$, where $k_{\perp}$ and $p_{\perp}$ represent the
transverse momenta of the parton with respect to the hadron. The SIDIS is one
of the main processes to measure TMD PDFs and hence to provide three-
dimensional partonic structures of the nucleon in the momentum space.

Taking the spin degrees of freedom into account, we can learn about much
richer nucleon structures, including the correlation between parton transverse
momentum and the spin of the parton or the nucleon. At the leading twist, one
can define eight TMD PDFs for the nucleon and only three of them have
collinear counterparts. The Sivers function $f_{1T}^{\perp}(x,k_{\perp})$, as
one of the eight leading-twist TMD PDFs, was originally introduced to explain
transverse single spin asymmetries in high energy scatterings [23, 24]. It
reflects the correlation between quark transverse momentum and the transverse
spin of the nucleon. The Sivers function, as well as the Boer-Mulders function
[25], is a naively time-reversal odd function, and thus , was believed to
vanish in QCD for a long time [26]. Motivated by some model calculation, it is
found that the nonvanishing Sivers function can arise from the final-state
interaction in the SIDIS process [27] or the initial-state interaction in the
Drell-Yan (DY) process [28]. This effect is formally from the nontrivial
Wilson line, which connects the quark field operator at different spacetime
points to ensure the color gauge invariance [29]. It is further proven in
perturbative QCD factorization that the quark Sivers function in the SIDIS
process and the one in DY process only differ by a sign,

$\displaystyle f_{1T}^{\perp}(x,k_{\perp})\Big{|}_{\rm
SIDIS}=-f_{1T}^{\perp}(x,k_{\perp})\Big{|}_{\rm DY},$ (3)

which receives great interest. Precise measurement of the Sivers function from
both processes will provide an important test of the perturbative QCD (pQCD).

During the last two decades, many SIDIS experiments have been carried out by
HERMES [30, 31, 32], COMPASS [33, 34, 35], and JLab [36, 37], for explorations
of nucleon TMD PDFs. The measured sizable target transverse single spin
asymmetry (SSA), usually referred to as the Sivers asymmetry, can be
interpreted at the leading twist via the convolution of the Sivers function
and the unpolarized TMD FF, and indicates a nonzero Sivers function. Although
it is the mostly studied polarized TMD PDF, the extractions of the Sivers
function from global analyses of the SSA data still have large uncertainties
[38, 39], especially for sea quarks. The sign change prediction was tested
through the $\pi N$ DY process [40] and the $pp$ collision with $W/Z$
production process [41]. The data do favor a sign flip of the Sivers function
[42], but the large uncertainties are not able to exclude nonflip cases and
other possibilities. All these pioneering experiments have provided valuable
data and revealed rich information about the Sivers function and more
generally the three-dimensional structures of the nucleon. However, our
current knowledge of the Sivers function is still far from a satisfactory
level, especially for sea quarks. Neither the shape of the function is well
determined, nor the sign is known for all flavors.

We are nowadays at the stage toward the precise measurement of TMD PDFs. The
multihall SIDIS program at the 12 GeV upgraded JLab aims at the measurement at
a relatively large-$x$ region, sometimes referred to as the valence quark
region. The Electron-Ion Collider (EIC) [43, 44] to be built at BNL has unique
advantages to reach the small-$x$ region, down to about $10^{-4}$. The
Electron-Ion Collider in China (EicC) [45] was recently proposed as one of the
next-generation colliders in the world. It is designed to deliver a $3.5\,\rm
GeV$ election beam colliding with a $20\,\rm GeV$ polarized proton beam or a
polarized ${}^{3}\rm He$ beam, as well as several other kinds of ion beams.
The ion beam can be longitudinally or transversely polarized with 70%
polarization, and the electron beam is polarized with 80% polarization. The
instantaneous luminosity of $ep$ collision can reach $2\times 10^{33}\,\rm
cm^{-2}\,s^{-1}$. The precise measurement of nucleon three-dimensional
structures, especially the sea quark distributions, is one of the main physics
goals at EicC. The combination of the polarized proton beam and the
effectively polarized neutron beam and the capability of particle
identifications of pions and kaons allows a complete separation of all light
quark flavors [46]. The kinematic coverage of EicC can fill the gap between
JLab and EIC. The combination of the three projects is expected to provide
complete three-dimensional imaging of the nucleon, varying from low scale to
high scale and from large $x$ to small $x$.

In this paper, we take the Sivers function as an example to investigate the
impact of the EicC three-dimensional nucleon spin structure program on the
determination of TMD PDFs. We perform a global analysis of existing world
SIDIS SSA data including the TMD evolution as the baseline. The improvement of
EicC is then estimated by adding simulated pseudodata of semi-inclusive
charged pion and charged kaon productions from both $ep$ and $e^{3}\rm He$
collisions. The paper is organized as follows. In Sec. II, we briefly
summarize the theoretical framework to extract the Sivers function from SIDIS
target transverse SSA data, leaving some detailed formulas in Appendix A. The
global analysis of world data and the EicC projection are presented in Sec.
III, followed by the summary in Sec. IV.

##  II Theoretical formalism

###  II.1 TMD factorization formula

We consider the SIDIS process (2) on a transversely polarized nucleon. With
the one-photon-exchange approximation, one can express the differential cross
section as

$\displaystyle\frac{d\sigma(S_{\perp})}{dx_{\mbox{\tiny\\!$B$}}dydz_{h}dP_{h\perp}^{2}d\phi_{h}d\phi_{S}}=\sigma_{0}\left[F_{UU}+\epsilon_{\perp\alpha\beta}S_{\perp}^{\alpha}F_{UT}^{\beta}+...\right],$
(4)

where

$\displaystyle\sigma_{0}=\frac{\alpha^{2}}{x_{\mbox{\tiny\\!$B$}}yQ^{2}}\frac{1-y+\frac{1}{2}y^{2}+\frac{1}{4}\gamma^{2}y^{2}}{1+\gamma^{2}}\left(1+\frac{\gamma^{2}}{2x_{\mbox{\tiny\\!$B$}}}\right),$
(5)

and $S_{\perp}^{\alpha}$ represents the transverse polarization of the
nucleon. As commonly used in the SIDIS process, we define the kinematic
variables

$\displaystyle Q^{2}$ $\displaystyle=-q^{2}=-(l-l^{\prime})^{2},$
$\displaystyle x_{\mbox{\tiny\\!$B$}}$ $\displaystyle=\frac{Q^{2}}{2P\cdot
q},\quad y=\frac{P\cdot q}{P\cdot l},\quad z_{h}=\frac{P\cdot P_{h}}{P\cdot
q},$ $\displaystyle\gamma$
$\displaystyle=\frac{2x_{\mbox{\tiny\\!$B$}}M}{Q}=\frac{MQ}{P\cdot q},$

where $l$ is the incoming lepton momentum, $l^{\prime}$ is the outgoing lepton
momentum, $P$ is the incoming nucleon momentum, $P_{h}$ is the detected
outgoing hadron momentum, and $M$ is the nucleon mass. The transverse
antisymmetric tensor is

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

with the convention $\epsilon^{0123}=1$. Following the Trento convention [47],
we define the hadron transverse momentum $P_{h\perp}$ and azimuthal angles in
the virtual photon-nucleon frame, as illustrated in Fig. 1. The $\phi_{h}$ is
the angle from the lepton plane to the hadron plane and the $\phi_{S}$ is the
angle from the lepton plane to the transverse spin $S_{\perp}$ of the nucleon.
These variables can also be expressed in Lorentz invariant forms as

$\displaystyle P_{h\perp}$
$\displaystyle=\sqrt{-g_{\perp}^{\mu\nu}P_{h\mu}P_{h\nu}},$ $\displaystyle
l_{\perp}$ $\displaystyle=\sqrt{-g_{\perp}^{\mu\nu}l_{\mu}l_{\nu}},$
$\displaystyle\cos\phi_{h}$
$\displaystyle=-\frac{l_{\mu}P_{h\nu}g_{\perp}^{\mu\nu}}{l_{\perp}P_{h\perp}},\quad\sin\phi_{h}=-\frac{l_{\mu}P_{h\nu}\epsilon_{\perp}^{\mu\nu}}{l_{\perp}P_{h\perp}},$
$\displaystyle\cos\phi_{S}$
$\displaystyle=-\frac{l_{\mu}S_{\perp\nu}g_{\perp}^{\mu\nu}}{l_{\perp}S_{\perp}},\quad\sin\phi_{S}=-\frac{l_{\mu}S_{\perp\nu}\epsilon_{\perp}^{\mu\nu}}{l_{\perp}S_{\perp}},$

where

$\displaystyle g_{\perp}^{\mu\nu}$
$\displaystyle=g^{\mu\nu}-\frac{q^{\mu}P^{\nu}+P^{\mu}q^{\nu}}{P\cdot
q\left(1+\gamma^{2}\right)}+\frac{\gamma^{2}}{1+\gamma^{2}}\left(\frac{q^{\mu}q^{\nu}}{Q^{2}}-\frac{P^{\mu}P^{\nu}}{M^{2}}\right).$
(7)

Figure 1: The Trento conventions of SIDIS kinematic variables.

The structure functions $F_{UU}$ and $F_{UT}^{\beta}$ are functions of
$x_{\mbox{\tiny\\!$B$}}$, $z_{h}$, $P_{h\perp}$, and $Q$. For the convenience
to apply the QCD factorization, one usually expresses the structure functions
in the “W+Y” formalism [48],

$\displaystyle F_{UU}(x_{\mbox{\tiny\\!$B$}},z_{h},P_{h\perp},Q)$
$\displaystyle=W_{UU}(x_{\mbox{\tiny\\!$B$}},z_{h},P_{h\perp},Q)$
$\displaystyle\quad+Y_{UU}(x_{\mbox{\tiny\\!$B$}},z_{h},P_{h\perp},Q),$ (8)
$\displaystyle F_{UT}^{\beta}(x_{\mbox{\tiny\\!$B$}},z_{h},P_{h\perp},Q)$
$\displaystyle=W_{UT}^{\beta}(x_{\mbox{\tiny\\!$B$}},z_{h},P_{h\perp},Q)$
$\displaystyle\quad+Y^{\beta}_{UT}(x_{\mbox{\tiny\\!$B$}},z_{h},P_{h\perp},Q).$
(9)

The “W” term dominates the cross section in the small $P_{hT}$ region, where
the TMD factorization works, and the “Y” term provides a matching to the
collinear factorization region at large $P_{hT}$. In this study, we restrict
to the small $P_{hT}$ region and neglect the “Y” term. The resummation formula
of the “W” term is derived in the impact parameter space via a Fourier
transform,

$\displaystyle W_{UU}(x_{\mbox{\tiny\\!$B$}},z_{h},P_{h\perp},Q)$
$\displaystyle=\int\frac{d^{2}\vec{b}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{b}}\widetilde{W}_{UU}(x_{\mbox{\tiny\\!$B$}},z_{h},b,Q),$
(10) $\displaystyle W_{UT}^{\beta}(x_{\mbox{\tiny\\!$B$}},z_{h},P_{h\perp},Q)$
$\displaystyle=\int\frac{d^{2}\vec{b}}{(2\pi)^{2}}e^{-i\vec{q}_{\perp}\cdot\vec{b}}\widetilde{W}_{UT}^{\beta}(x_{\mbox{\tiny\\!$B$}},z_{h},b,Q),$
(11)

where $q_{\perp}=P_{h\perp}/z_{h}$ can be understood as the transverse
momentum of the virtual photon in the nucleon-hadron back-to-back frame.

According to the TMD factorization, one has the factorized expressions as

$\displaystyle\widetilde{W}_{UU}(x_{\mbox{\tiny\\!$B$}},z_{h},b,Q)$
$\displaystyle=H(\mu,Q)\otimes f_{1}(x,b;\mu,\zeta_{1})\otimes
D_{1}(z,b;\mu,\zeta_{2}),$
$\displaystyle\widetilde{W}_{UT}^{\beta}(x_{\mbox{\tiny\\!$B$}},z_{h},b,Q)$
$\displaystyle=H(\mu,Q)\otimes
b^{\beta}Mf^{\bot}_{1T}(x,b;\mu,\zeta_{1})\otimes D_{1}(z,b;\mu,\zeta_{2}),$
(12)

with corrections suppressed by powers of $q_{\perp}/Q$. The $H(\mu,Q)$ is the
hard factor describing the short-distance scattering. In the calculation of
the nucleon transverse SSA, the hard factors cancel out between the
denominator and the numerator for fixed $\mu$ and $Q$, and hence we do not
consider its effect in this work. $f_{1}(x,b;\mu,\zeta)$,
$f^{\perp}_{1T}(x,b;\mu,\zeta)$ and $D_{1}(z,b;\mu,\zeta)$ are, respectively,
the unpolarized TMD PDF, the Sivers function, and the unpolarized TMD FF in
the $b$ space with $x$ representing the longitudinal momentum fraction carried
by the active quark inside the nucleon and $z$ representing the longitudinal
momentum fraction carried by the hadron from the fragmented quark. They are
related to the corresponding functions in transverse momentum space via a
Fourier transform,

$\displaystyle f_{1}(x,k_{\perp};\mu,\zeta)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{db\,b}{2\pi}J_{0}(bk_{\perp})f_{1}(x,b;\mu,\zeta)\
,$ $\displaystyle\frac{k_{\perp}}{M}f^{\perp}_{1T}(x,k_{\perp};\mu,\zeta)$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{db\,b^{2}M}{2\pi}J_{1}(bk_{\perp})f^{\perp}_{1T}(x,b;\mu,\zeta)\
,$ $\displaystyle D_{1}(z,p_{\perp};\mu,\zeta)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{db\,b}{2\pi}J_{0}(bp_{\perp})D_{1}(z,b;\mu,\zeta)\
,$ (13)

where $J_{0,1}$ are the Bessel functions. The $\mu$ is the factorization
scale, which one can choose arbitrarily. The scales $\zeta_{1}$ and
$\zeta_{2}$ in Eq. (12) are related to the invariant mass of virtual photon by
$\zeta_{1}\zeta_{2}=Q^{4}$, and the “W” term only depends on their product
$\zeta_{1}\zeta_{2}$. For convenience, we set the scales in the calculation as

$\displaystyle\mu^{2}=\zeta_{1}=\zeta_{2}=Q^{2}.$ (14)

The convolutions over $x$ and $z$ in (12) are indicated by “$\otimes$”. At the
leading order, i.e., the tree level, we have $x_{\mbox{\tiny\\!$B$}}=x$ and
$z_{h}=z$.

###  II.2 Energy evolution and TMD resummation

In perturbation theory, one expands the factorized formula (12) in powers of
the strong coupling constant $\alpha_{s}$. However, at each fixed order, e.g.,
the $n$th order, the TMD functions $f_{1}(x,b;Q,Q^{2})$,
$f_{1T}^{\perp}(x,b;Q,Q^{2})$, and $D_{1}(z,b;Q,Q^{2})$ contain logarithmic
enhanced terms of $[\alpha_{s}\ln^{2}(Qb)]^{n}$ and $[\alpha_{s}\ln(Qb)]^{n}$.
In order to obtain a reliable prediction from perturbation theory, one needs
to take into account the large logarithms of all orders. This can be achieved
by evolving the scales $Q$ and $Q^{2}$ to $\mu_{i}\sim 1/b$ and $\zeta_{i}\sim
1/b^{2}$ and in the meantime resuming the large logarithms into an evolution
factor $R[b;(\mu_{i},\zeta_{i})\to(Q,Q^{2})]$. Then one can formally relate
the TMD functions at $(Q,Q^{2})$ and $(\mu_{i},\zeta_{i})$ as

$\displaystyle f_{1}(x,b;Q,Q^{2})D_{1}(z,b;Q,Q^{2})$ $\displaystyle=$
$\displaystyle
R^{2}[b;(\mu_{i},\zeta_{i})\to(Q,Q^{2})]f_{1}(x,b;\mu_{i},\zeta_{i})D_{1}(z,b;\mu_{i},\zeta_{i}),$
$\displaystyle f^{\bot}_{1T}(x,b;Q,Q^{2})D_{1}(z,b;Q,Q^{2})$ $\displaystyle=$
$\displaystyle
R^{2}[b;(\mu_{i},\zeta_{i})\to(Q,Q^{2})]f^{\bot}_{1T}(x,b;\mu_{i},\zeta_{i})D_{1}(z,b;\mu_{i},\zeta_{i}).$
(15)

The TMD evolution equations are given by

$\displaystyle\frac{d}{d\ln\mu}f_{1}(x,b;\mu,\zeta)=\gamma_{F}(\mu,\zeta)f_{1}(x,b;\mu,\zeta),$
(16)
$\displaystyle\frac{d}{d\ln\zeta}f_{1}(x,b;\mu,\zeta)=-\mathcal{D}(\mu,b)f_{1}(x,b;\mu,\zeta),$
(17)

where $\gamma_{F}(\mu,\zeta)$ is the TMD anomalous dimension and
$\mathcal{D}(\mu,b)$ is the rapidity anomalous dimension, with similar
equations for $f_{1T}^{\perp}(x,b;\mu,\zeta)$ and $D_{1}(z,b;\mu,\zeta)$. By
solving the equations above, one can express the factor
$R[b;(\mu_{i},\zeta_{i})\to(Q,Q^{2})]$ as a path integral from
$(\mu_{i},\zeta_{i})$ to $(Q,Q^{2})$,

$\displaystyle R[b;(\mu_{i},\zeta_{i})\to(Q,Q^{2})]$
$\displaystyle=\exp\left[\int_{\cal
P}\Big{(}\frac{\gamma_{F}(\mu,\zeta)}{\mu}d\mu-\frac{\mathcal{D}(\mu,b)}{\zeta}d\zeta\Big{)}\right],$
(18)

where the path $\cal P$ can be chosen arbitrarily because of the integrability
condition

$\displaystyle\zeta\frac{d}{d\zeta}\gamma_{F}(\mu,\zeta)=-\mu\frac{d}{d\mu}{\cal
D}(\mu,b)=-\Gamma_{\rm cusp}(\mu),$ (19)

where $\Gamma_{\rm cusp}(\mu)$ is the cusp anomalous dimension. Then the
anomalous dimension $\gamma_{F}(\mu,\zeta)$ can be written as

$\displaystyle\gamma_{F}(\mu,\zeta)=\Gamma_{\rm
cusp}(\mu)\ln\Big{(}\frac{\mu^{2}}{\zeta}\Big{)}-\gamma_{V}(\mu),$ (20)

where $\gamma_{V}(\mu)$ is the finite part of the renormalization of the
vector form factor. Expanding these factors in powers of the strong coupling
constant,

$\displaystyle\Gamma_{\rm cusp}(\mu)$
$\displaystyle=\sum_{n=0}^{\infty}a_{s}^{n+1}\Gamma_{n},$ (21)
$\displaystyle\gamma_{V}(\mu)$
$\displaystyle=\sum_{n=1}^{\infty}a_{s}^{n}\gamma_{n},$ (22)

with $a_{s}=\alpha_{s}/(4\pi)$ introduced for convenience, one can calculate
the coefficients $\Gamma_{n}$ and $\gamma_{n}$ perturbatively when
$\mu\gg\Lambda_{\rm QCD}$, and up to two-loop level,

$\displaystyle\Gamma_{0}$ $\displaystyle=4C_{F},$ (23)
$\displaystyle\Gamma_{1}$
$\displaystyle=4C_{F}\Big{[}\big{(}\frac{67}{9}-\frac{\pi^{2}}{3}\big{)}C_{A}-\frac{20}{9}T_{R}N_{f}\Big{]},$
(24) $\displaystyle\gamma_{1}$ $\displaystyle=-6C_{F},$ (25)
$\displaystyle\gamma_{2}$ $\displaystyle=C_{F}^{2}(-3+4\pi^{2}-48\zeta_{3})$
$\displaystyle+C_{F}C_{A}\Big{(}-\frac{961}{27}-\frac{11\pi^{2}}{3}+52\zeta_{3}\Big{)}$
$\displaystyle+C_{F}T_{R}N_{f}\Big{(}\frac{260}{27}+\frac{4\pi^{2}}{3}\Big{)},$
(26)

where $C_{F}=4/3$, $C_{A}=3$, and $T_{R}=1/2$ are color factors of the
$SU(3)$, $N_{f}=4$ is the number of active quark flavors, and
$\zeta_{3}\approx 1.202$ is the Apéry’s constant. Similarly, one can expand
$\mathcal{D}(\mu,b)$ in powers of $a_{s}$ as,

$\displaystyle\mathcal{D}_{\rm
pert}(\mu,b)=\sum_{n=0}^{\infty}a_{s}^{n}d_{n}(\textbf{L}_{\mu}),$ (27)

which is only valid when $1/b\gg\Lambda_{\rm QCD}$, as indicated by the
subscript “pert”. The $\textbf{L}_{\mu}$ is defined as

$\displaystyle\textbf{L}_{\mu}=\ln(\frac{\mu^{2}b^{2}}{4e^{-2\gamma_{E}}}),$
(28)

with $\gamma_{E}$ the Euler-Mascheroni constant. Up to two-loop order, one has

$\displaystyle d_{0}(\textbf{L}_{\mu})$ $\displaystyle=0,$ (29) $\displaystyle
d_{1}(\textbf{L}_{\mu})$ $\displaystyle=\frac{\Gamma_{0}}{2}\textbf{L}_{\mu},$
(30) $\displaystyle d_{2}(\textbf{L}_{\mu})$
$\displaystyle=\frac{\Gamma_{0}}{4}\beta_{0}\textbf{L}_{\mu}^{2}+\frac{\Gamma_{1}}{2}\textbf{L}_{\mu}+d_{2}(0).$
(31)

where

$\displaystyle
d_{2}(0)=C_{F}C_{A}\Big{(}\frac{404}{27}-14\zeta_{3}\Big{)}-\frac{112}{27}T_{R}N_{f}C_{F}$
(32)

Following the treatment in [49], one can substitute $\mathcal{D}_{\rm
pert}(\mu,b)$ into (19) to get a set of equations of $d_{n}({\bf L}_{\mu})$.
By solving these equations, one can sum the series of logarithmic terms
$X^{n}$ for each $d_{n}(\textbf{L}_{\mu})$, with
$X=\beta_{0}a_{s}\textbf{L}_{\mu}$, and obtain the resummed expression,

$\displaystyle\mathcal{D}_{\rm
resum}(\mu,b)=-\frac{\Gamma_{0}}{2\beta_{0}}\ln(1-X)$
$\displaystyle+\frac{a_{s}}{2\beta_{0}(1-X)}\Big{[}-\frac{\beta_{1}\Gamma_{0}}{\beta_{0}}(\ln(1-X)+X)+\Gamma_{1}X\Big{]}$
$\displaystyle+\frac{a_{s}^{2}}{(1-X)^{2}}\Big{[}\frac{\Gamma_{0}\beta_{1}^{2}}{4\beta_{0}^{3}}(\ln^{2}(1-X)-X^{2})$
$\displaystyle+\frac{\beta_{1}\Gamma_{1}}{4\beta_{0}^{2}}\big{(}X^{2}-2X-2\ln(1-X)\big{)}$
$\displaystyle+\frac{\Gamma_{0}\beta_{2}}{4\beta_{0}^{2}}X^{2}-\frac{\Gamma_{2}}{4\beta_{0}}X(X-2)$
$\displaystyle+C_{F}C_{A}\Big{(}\frac{404}{27}-14\zeta_{3}\Big{)}-\frac{112}{27}T_{R}N_{f}C_{F}\Big{]}.$
(33)

However, the perturbative expansion is invalid when $b\gtrsim 1/\Lambda_{\rm
QCD}$, and hence one can only apply the above expression up to some $b$ value.
For a smooth connection, one usually introduces the modified $b$ variable as

$\displaystyle b^{*}=\frac{b}{\sqrt{1+b^{2}/B^{2}_{\rm NP}}},$ (34)

which serves as an effective truncation since $b^{*}<B_{\rm NP}$ holds for
arbitrary large $b$ and $b^{*}\approx b$ for small $b$. Meanwhile, the large
$b$ behavior of ${\cal D}(\mu,b)$ is modeled by a function $d_{\rm NP}$, which
we adopt the form in Ref. [50],

$\displaystyle d_{\rm NP}(b)=c_{0}bb^{*},$ (35)

which is linear function at large-$b$ region as suggested in Refs. [51, 52,
53, 54]. We set the free parameters $B_{\rm NP}=1.93\,\rm GeV^{-1}$ and
$c_{0}=0.0391\,\rm GeV^{2}$ as determined in Ref. [55] by fitting unpolarized
SIDIS and DY data. Finally, the $\mathcal{D}(\mu,b)$ is expressed as

$\displaystyle\mathcal{D}(\mu,b)=\mathcal{D}_{\rm resum}(\mu,b^{*})+d_{\rm
NP}(b)\ .$ (36)

According to the $\zeta$ prescription [55], by solving the equation,

$\displaystyle\frac{d\ln\zeta_{\mu}(\mu,b)}{d\ln\mu^{2}}=\frac{\gamma_{F}(\mu,\zeta_{\mu}(\mu,b))}{2\mathcal{D}(\mu,b)}\
,$ (37)

with the boundary conditions,

$\displaystyle\mathcal{D}(\mu_{0},b)=0\
,\quad\gamma_{F}(\mu_{0},\zeta_{\mu}(\mu_{0},b))=0\ ,$ (38)

one can find a special point $\zeta_{\mu}(Q,b)$, so that

$\displaystyle
R[b;(\mu_{i},\zeta_{i})\to(Q,Q^{2})]=\Big{(}\frac{Q^{2}}{\zeta_{\mu}(Q,b)}\Big{)}^{-\mathcal{D}(Q,b)}\
.$ (39)

Then the Eq. (37) can be rewritten as

$\displaystyle\frac{d\ln\zeta_{\mu}(\mathcal{D}(\mu,b),b)}{d\mathcal{D}(\mu,b)}\Gamma_{\rm
cusp}=\frac{\gamma_{F}(\mu,\zeta_{\mu}(\mu,b))}{2\mathcal{D}(\mu,b)}\ .$ (40)

Using Eq. (36) as an input, one can obtain a solution that is independent of
the form of $\mathcal{D}(\mu,b)$ as

$\displaystyle\zeta^{\rm exact}_{\mu}(\mu,b)$
$\displaystyle=\mu^{2}e^{-g(\mu,b)/\mathcal{D}(\mu,b)}\ .$ (41)

Up to two-loop order, $g(\mu,b)$ can be written as

$\displaystyle
g(\mu,b)=\frac{1}{a_{s}}\frac{\Gamma_{0}}{2\beta_{0}^{2}}\Bigg{\\{}e^{-p}-1+p+a_{s}\bigg{[}\frac{\beta_{1}}{\beta_{0}}\big{(}e^{-p}-1+p-\frac{p^{2}}{2}\big{)}$
$\displaystyle-\frac{\Gamma_{1}}{\Gamma_{0}}(e^{-p}-1+p)+\frac{\beta_{0}\gamma_{1}}{\Gamma_{0}}p\bigg{]}+a_{s}^{2}\bigg{[}\Big{(}\frac{\Gamma_{1}^{2}}{\Gamma_{0}^{2}}-\frac{\Gamma_{2}}{\Gamma_{0}}\Big{)}(\cosh
p-1)$
$\displaystyle+\Big{(}\frac{\beta_{1}\Gamma_{1}}{\beta_{0}\Gamma_{0}}-\frac{\beta_{2}}{\beta_{0}}\Big{)}(\sinh
p-p)+\Big{(}\frac{\beta_{0}\gamma_{2}}{\Gamma_{0}}-\frac{\beta_{0}\gamma_{1}\Gamma_{1}}{\Gamma_{0}^{2}}\Big{)}(e^{p}-1)\bigg{]}\Bigg{\\}}\
,$ (42)

where

$\displaystyle p=\frac{2\beta_{0}\mathcal{D}(\mu,b)}{\Gamma_{0}}\ .$ (43)

At extremely small $b$ region, there could be numerical difficulties for the
solution (41) to reduce to match the pQCD prediction [55]. On the other hand,
$\mathcal{D}(\mu,b)$ will reduce to $\mathcal{D}_{\rm pert}(Q,b)$ since it is
determined by pQCD in such a region. Instead of ${\cal D}(\mu,b)$, one use
$\mathcal{D}_{\rm pert}(Q,b)$ as an input to solve Eq. (37) and will obtain
the solution,

$\displaystyle\zeta^{\rm pert}_{\mu}(\mu,b)$ $\displaystyle=\frac{2\mu
e^{-\gamma_{E}}}{b}e^{-v(\mu,b)}\ ,$ (44)

which is consistent with the pQCD result by construction. Up to two-loop
order, $v(\mu,b)$ can be written as

$\displaystyle
v(\mu,b)=\frac{\gamma_{1}}{\Gamma_{0}}+a_{s}\Big{[}\frac{\beta_{0}}{12}\textbf{L}_{\mu}^{2}+\frac{\gamma_{2}+d_{2}(0)}{\Gamma_{0}}-\frac{\gamma_{1}\Gamma_{1}}{\Gamma_{0}^{2}}\Big{]}\
.$ (45)

To match the large-$b$ region solution $\zeta_{\mu}^{\rm exact}(\mu,b)$ and
the small-$b$ region solution $\zeta_{\mu}^{\rm pert}(\mu,b)$, we express
$\zeta_{\mu}(\mu,b)$ as [55]:

$\displaystyle\zeta_{\mu}(\mu,b)$ $\displaystyle=\zeta^{\rm
pert}_{\mu}(\mu,b)e^{-b^{2}/B_{\rm NP}^{2}}$ $\displaystyle\quad+\zeta^{\rm
exact}_{\mu}(\mu,b)\Big{(}1-e^{-b^{2}/B_{\rm NP}^{2}}\Big{)}\ .$ (46)

###  II.3 Unpolarized TMD PDF and FF

Taking the phenomenological ansatzes in Ref. [55], one can express the
unpolarized TMD PDF and FF as

$\displaystyle f_{1,f\leftarrow
h}(x,b;\mu_{i},\zeta_{i})=\sum_{f^{\prime}}\int_{x}^{1}\frac{dy}{y}C_{f\leftarrow
f^{\prime}}(y,b,\mu_{\rm OPE}^{\rm PDF})$
$\displaystyle\quad\quad\quad\quad\times f_{1,f^{\prime}\leftarrow
h}\Big{(}\frac{x}{y},\mu_{\rm OPE}^{\rm PDF}\Big{)}f_{\rm NP}(x,b),$ (47)
$\displaystyle D_{1,f\to
h}(z,b;\mu_{i},\zeta_{i})=\frac{1}{z^{2}}\sum_{f^{\prime}}\int_{z}^{1}\frac{dy}{y}y^{2}\mathbb{C}_{f\to
f^{\prime}}(y,b,\mu_{\rm OPE}^{\rm FF})$
$\displaystyle\quad\quad\quad\quad\times d_{1,f^{\prime}\to
h}\Big{(}\frac{z}{y},\mu_{\rm OPE}^{\rm FF}\Big{)}D_{\rm NP}(z,b),$ (48)

where $f_{\rm NP}(x,b)$ and $D_{\rm NP}(z,b)$ are nonperturbative functions,
$f_{1,f^{\prime}\leftarrow h}(x,\mu)$ and $d_{1,f^{\prime}\to h}(x,\mu)$ are
collinear PDF and FF, and $C_{f\leftarrow f^{\prime}}(y,b,\mu)$ and
$\mathbb{C}_{f\to f^{\prime}}(y,b,\mu)$ are matching coefficients that can be
calculated via pQCD. In this work, we consider the $C(\mathbb{C})$ functions
up to the one-loop order, and the expressions are provided in Appendix A. The
evolution of scales $\mu_{\rm OPE}^{\rm PDF}$ and $\mu_{\rm OPE}^{\rm FF}$ in
the $\zeta$ prescription are independent of external parameters. Here we adopt
the choices [55]

$\displaystyle\mu_{\rm OPE}^{\rm PDF}$
$\displaystyle=\frac{2e^{-\gamma_{E}}}{b}+2\,{\rm GeV}\ ,$ (49)
$\displaystyle\mu_{\rm OPE}^{\rm FF}$
$\displaystyle=\frac{2e^{-\gamma_{E}}z}{b}+2\,\rm GeV\ ,$ (50)

and $2\,\rm GeV$ is a typical reference scale for PDFs and FFs.
Correspondingly, we take the parametrization of the nonperturbative functions
$f_{\rm NP}(x,b)$ and $D_{\rm NP}(z,b)$ as [55]

$\displaystyle f_{\rm NP}(x,b)$
$\displaystyle=\exp\Big{[}-\frac{\lambda_{1}(1-x)+\lambda_{2}x+x(1-x)\lambda_{5}}{\sqrt{1+\lambda_{3}x^{\lambda_{4}}b^{2}}}b^{2}\Big{]},$
(51) $\displaystyle D_{\rm NP}(z,b)$
$\displaystyle=\exp\Big{[}-\frac{\eta_{1}z+\eta_{2}(1-z)}{\sqrt{1+\eta_{3}(b/z)^{2}}}\frac{b^{2}}{z^{2}}\Big{]}\big{(}1+\eta_{4}\frac{b^{2}}{z^{2}}\big{)},$
(52)

where the parameters $\lambda$’s and $\eta$’s parameters $\lambda$ and $\eta$
are extracted from unpolarized SIDIS and Drell-Yan data at small transverse
momentum. Their values in this study are listed in Table 1.

Table 1: The values of the parameters in the parametrization of unpolarized TMD PDF and FF. Their units are in GeV2 except for $\lambda_{4}$.  $\ \ \ \ \ \lambda_{1}$ |  $\lambda_{2}$ |  $\lambda_{3}$ |  $\lambda_{4}$ |  $\lambda_{5}$
---|---|---|---|---
0.198 | 9.30 | 431 | 2.12 | -4.44
$\eta_{1}$ | $\eta_{2}$ | $\eta_{3}$ | $\eta_{4}$ |
0.260 | 0.476 | 0.478 | 0.483 |

###  II.4  Sivers asymmetry

The Sivers asymmetry in SIDIS is a transverse single spin asymmetry. It can be
generally defined as

$\displaystyle A_{UT}^{\sin(\phi_{h}-\phi_{S})}=$ $\displaystyle\frac{2\int
d\phi_{S}d\phi_{h}[d\sigma^{\uparrow}-d\sigma^{\downarrow}]\sin(\phi_{h}-\phi_{S})}{\int
d\phi_{S}d\phi_{h}[d\sigma^{\uparrow}+d\sigma^{\downarrow}]},$ (53)

where $\sigma^{\uparrow}$ and $\sigma^{\downarrow}$ represent the cross
sections from transversely polarized nucleon. The azimuthal angles $\phi_{S}$
and $\phi_{h}$ are defined following the Trento convention as illustrated in
Fig. 1.

Table 2: The precision of the factors in powers of $\alpha_{s}$ in this work. |  $\Gamma_{cusp}$ |  $\gamma_{V}$ | $\mathcal{D}_{resum}$ |  $\zeta_{\mu}^{pert}$ |  $\zeta_{\mu}^{exact}$ |  $C(\mathbb{C}$)
---|---|---|---|---|---|---
$F_{UU}$ | $\alpha_{s}^{3}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{1}$ | $\alpha_{s}^{1}$ | $\alpha_{s}^{1}$
$F_{UT}$ | $\alpha_{s}^{3}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{1}$ | $\alpha_{s}^{1}$ | $\alpha_{s}^{0}$

Within the TMD formalism presented above, one can express the Sivers asymmetry
$A_{UT}^{\sin(\phi_{h}-\phi_{S})}$ as

$\displaystyle
A_{UT}^{\sin(\phi_{h}-\phi_{S})}=-\frac{M\sum_{q}e^{2}_{q}\int_{0}^{\infty}\frac{bdb}{2\pi}bJ_{1}(\frac{b|P_{hT}|}{z})R^{2}[b;(\mu_{i},\zeta_{i})\to(Q,Q^{2})]f^{\bot}_{1T,q\leftarrow
h_{1}}(x,b)D_{1,q\to
h_{2}}(z,b)}{\sum_{q}e^{2}_{q}\int_{0}^{\infty}\frac{bdb}{2\pi}J_{0}(\frac{b|P_{hT}|}{z})R^{2}[b;(\mu_{i},\zeta_{i})\to(Q,Q^{2})]f_{1,q\leftarrow
h_{1}}(x,b)D_{1,q\to h_{2}}(z,b)}\ .$ (54)

The precision for the perturbative calculation of the factors in powers of
$\alpha_{s}$ in this work is summarized in Table 2.

Table 3: World SIDIS data used in our analysis. The numbers in parentheses are the original number of data points before applying $\delta$ cut. Dataset |  Target |  Beam |  Data points |  Reaction
---|---|---|---|---
COMPASS [33] |  ${}^{6}\text{LiD}$ |  $160\,\rm GeV$ $\mu^{+}$ |  10(41) |  $\mu^{+}d\to\mu^{+}\pi^{+}X$
|  |  |  |  $\mu^{+}d\to\mu^{+}\pi^{-}X$
|  |  |  |  $\mu^{+}d\to\mu^{+}K^{+}X$
|  |  |  |  $\mu^{+}d\to\mu^{+}K^{-}X$
|  |  |  |  $\mu^{+}d\to\mu^{+}K^{0}X$
COMPASS [34] |  $\text{NH}_{3}$ |  $160\,\rm GeV$ $\mu^{+}$ |  10(36) |  $\mu^{+}p\to\mu^{+}\pi^{+}X$
|  |  |  |  $\mu^{+}p\to\mu^{+}\pi^{-}X$
|  |  |  |  $\mu^{+}p\to\mu^{+}K^{+}X$
|  |  |  |  $\mu^{+}p\to\mu^{+}K^{-}X$
HERMES [32] |  $\text{H}_{2}$ |  $27.6\,\rm GeV$ $e^{\pm}$ |  104(256) |  $e^{\pm}p\to e^{\pm}\pi^{+}X$
|  |  |  |  $e^{\pm}p\to e^{\pm}\pi^{-}X$
|  |  |  |  $e^{\pm}p\to e^{\pm}K^{+}X$
|  |  |  |  $e^{\pm}p\to e^{\pm}K^{-}X$
JLab [36] |  ${}^{3}\text{He}$ |  $5.9\,\rm GeV$ $e^{-}$ |  4(8) |  $e^{-}n\to e^{-}\pi^{+}X$
|  |  |  |  $e^{-}n\to e^{-}\pi^{-}X$
JLab [37] |  ${}^{3}\text{He}$ |  $5.9\,\rm GeV$ $e^{-}$ |  2(5) |  $e^{-}{}^{3}{\rm He}\to e^{-}K^{+}X$
|  |  |  |  $e^{-}{}^{3}{\rm He}\to e^{-}K^{-}X$

##  III Extract the Sivers function from the SIDIS data

We parametrize the Sivers function of each quark flavor as

$\displaystyle f_{1T;q\leftarrow
p}^{\perp}(x,b)=N_{q}\frac{(1-x)^{\alpha_{q}}x^{\beta_{q}}(1+\epsilon_{q}x)}{n(\beta_{q},\epsilon_{q},\alpha_{q})}\exp\big{(}-r_{q}b^{2}\big{)},$
(55)

where $N_{q}$, $\alpha_{q}$, $\beta_{q}$, $\epsilon_{q}$, and $r_{q}$ are free
parameters. The factor

$\displaystyle n(\beta,\epsilon,\alpha)$
$\displaystyle=\frac{\Gamma(\alpha+1)(2+\alpha+\beta+\epsilon+\epsilon\beta)\Gamma(\beta+1)}{\Gamma(\beta+\alpha+3)}$

is introduced to reduce the correlation between the parameters controlling the
shape and the parameter for the normalization.

To extract the Sivers function from the SIDIS SSA data, one also needs
unpolarized TMD PDF $f_{1}(x,b)$ and unpolarized TMD FF $D_{1}(z,b)$, which we
take from recent global analysis [55], as mentioned in the previous section.
For simplicity, we only consider nonvanishing Sivers functions for light
flavors, i.e., $u$, $d$, $s$, $\bar{u}$, $\bar{d}$, and $\bar{s}$, while heavy
quark Sivers functions are all set to zero. We also assume the isospin
symmetry to relate the Sivers function of the neutron and the Sivers function
of the proton as

$\displaystyle f^{\perp}_{1T,u\leftarrow n}(x,b)$
$\displaystyle=f^{\perp}_{1T,d\leftarrow p}(x,b),$ $\displaystyle
f^{\perp}_{1T,\bar{u}\leftarrow n}(x,b)$
$\displaystyle=f^{\perp}_{1T,\bar{d}\leftarrow p}(x,b),$ $\displaystyle
f^{\perp}_{1T,d\leftarrow n}(x,b)$ $\displaystyle=f^{\perp}_{1T,u\leftarrow
p}(x,b),$ $\displaystyle f^{\perp}_{1T,\bar{d}\leftarrow n}(x,b)$
$\displaystyle=f^{\perp}_{1T,\bar{u}\leftarrow p}(x,b),$ $\displaystyle
f^{\perp}_{1T,s\leftarrow n}(x,b)$ $\displaystyle=f^{\perp}_{1T,s\leftarrow
p}(x,b),$ $\displaystyle f^{\perp}_{1T,\bar{s}\leftarrow n}(x,b)$
$\displaystyle=f^{\perp}_{1T,\bar{s}\leftarrow p}(x,b).$ (56)

Since a free neutron target is not available for SIDIS experiments, the
polarized deuteron and polarized ${}^{3}\rm He$ are commonly used to obtain
parton distributions in the neutron. As an approximation, we set the Sivers
functions of the deuteron and the ${}^{3}\rm He$ via the weighted combination
of the proton Sivers function and the neutron Sivers function. For the
deuteron, the Sivers function is expressed as

$\displaystyle f^{\perp}_{1T,q\leftarrow d}(x,b)$
$\displaystyle=\frac{P_{d}^{n}f^{\perp}_{1T,q\leftarrow
n}(x,b)+P_{d}^{p}f^{\perp}_{1T,q\leftarrow p}(x,b)}{2},$ (57)

where $P_{d}^{n}=P_{d}^{p}=0.925$ are effective polarizations of the neutron
and the proton in a polarized deuteron [56]. Similarly, the Sivers function of
the ${}^{3}\rm He$ is

$\displaystyle f^{\perp}_{1T,q\leftarrow^{3}\rm He}(x,b)$
$\displaystyle=\frac{P_{\rm He}^{n}f^{\perp}_{1T,q\leftarrow n}(x,b)+2P_{\rm
He}^{p}f^{\perp}_{1T,q\leftarrow p}(x,b)}{3},$ (58)

where $P_{\rm He}^{n}=0.86$ and $P_{\rm He}^{p}=-0.028$ are effective
polarizations of the neutron and the proton in a polarized ${}^{3}\rm He$
[57]. Here we have followed the convention to normalize the distribution to
per nucleon.

This parametrization setup is applied in the following for both the fit to
world SIDIS data and the fit to EicC pseudodata.

###  III.1 Fit to world SIDIS data

We first perform the fit to existing SIDIS target transverse SSA data from
HERMES [32], COMPASS [33, 34], and JLab [36, 37]. For the COMPASS
measurements, we only include the analyzed data with particle identifications
of the final state hadron for the purpose of flavor separation, although the
results without hadron particle identifications for the same datasets have
also been published [35]. For the validity of the TMD factorization, only
small transverse momentum data with $\delta=|P_{h\perp}|/(zQ)<0.5$ are
selected. After applying this cut, there are 130 data points, as summarized in
Table 3. These data are all from fixed-target experiments with various beam
energies: $160\,\rm GeV$ muon beam for COMPASS, $27.6\,\rm GeV$
electron/positron beam for HERMES, and $5.9\,\rm GeV$ electron beam for JLab.

Since the existing world data are not precise enough to constrain so many
parameters introduced in (55) for all light flavors, we reduce the number of
parameters by setting $s$ and $\bar{s}$ quarks Sivers function to $0$ and
imposing the following assumptions for $u,d,\bar{u},\bar{d}$ quarks:

$\begin{split}&r_{\bar{u}}=r_{\bar{d}}=r_{\rm
sea},\quad\beta_{\bar{u}}=\beta_{\bar{d}}=\beta_{\rm sea},\\\
&\alpha_{u}=\alpha_{d}=3,\quad\alpha_{\bar{u}}=\alpha_{\bar{d}}=5,\quad\epsilon_{\bar{u}}=\epsilon_{\bar{d}}=0.\end{split}$
(59)

Then we have $12$ free parameters, as listed in Table 4. To estimate the
uncertainty, we randomly shift the central values of the data points by
Gaussian distributions with the Gaussian widths to be the experimental
uncertainties and then perform the fit to each smeared dataset. By repeating
this procedure, 100 replicas are created for calculations of central values
and uncertainties of physical quantities related to the Sivers function. The
total $\chi^{2}/N$ of the fit as well as its value for various experimental
datasets are listed in table 5. Here, $N$ is the number of experimental data
points. The details of the calculation are explained in Appendix B.

The central values of the parameters together with their uncertainties are
listed in Table 4. The $\beta$ parameters turn out to be negative for up and
down quarks while positive for sea quarks, which gives us some hints that in
small-$x$ region the Sivers effect for sea quarks is weaker than that of up
and down quarks. On the other hand, existing world data are not precise
enough, and therefore a decisive conclusion will rely on future data from
electron-ion colliders. The comparisons between experimental data and the
calculations by using replicas are shown in Figs. 2-5, where the filled data
points with $\delta<0.5$ are included in the fit while the open data points
with $\delta>0.5$ are not. The results of the Sivers functions are shown in
Fig. 7 via slices at various $x$ values. For better visualization of the $x$
dependence, we also present ${\bf k}_{\perp}$-integrated distributions of the
Sivers function via its zeroth and first transverse momentum moments,

$\displaystyle f_{1T}^{\perp(0)}(x)$ $\displaystyle=\pi\int d{\bf
k}_{\perp}^{2}\,f_{1T}^{\perp}(x,{\bf k}_{\perp}^{2}),$ (60) $\displaystyle
f_{1T}^{\perp(1)}(x)$ $\displaystyle=\pi\int d{\bf k}_{\perp}^{2}\,\frac{{\bf
k}_{\perp}^{2}}{2M^{2}}f_{1T}^{\perp}(x,{\bf k}_{\perp}^{2}).$ (61)

Since TMDs are well defined at small transverse momentum and the fit only
includes data in the small transverse momentum region, we truncate the
integrals at ${\bf k}_{\perp}^{\rm max}=0.6\,\rm GeV$. The truncated zeroth
and first transverse momentum moments are shown in Fig. 8 and Fig. 9
respectively. Our results are compared with BPV20 [39] in Fig. 10. Within
uncertainty, the results are consistent with each other. However, the
uncertainties of sea quark Sivers functions in Ref. [39] are much smaller than
our estimations because of strong assumptions in the parametrization. To be
more specific, in Ref. [39] the intrinsic $k_{\perp}$ dependence for all the
quarks is assumed to be the same, and additionally, the Sivers functions of
sea quarks ($\bar{u}$, $\bar{d}$, $\bar{s}$, ${s}$) are assumed to be
identical, except for a separate normalization parameter for the $s$ quark.

Table 4: The values of the parameters from the fit to world SIDIS SSA data. The central values are the average of the results from 100 fits, and the uncertainties are the standard deviations. The values of $r_{u}$, $r_{d}$, and $r_{\rm sea}$ are provided in unit of $\rm GeV^{2}$ and the others are dimensionless.  Parameter |  Value |  Parameter |  Value
---|---|---|---
$r_{u}$ |  $0.08^{+0.04}_{-0.03}$ |  $N_{u}$ |  $-0.08^{+0.02}_{-0.04}$
$r_{d}$ |  $0.2^{+0.9}_{-0.2}$ |  $N_{d}$ |  $1.0^{+2.6}_{-0.5}$
$r_{\rm sea}$ |  $0.2^{+1.5}_{-0.2}$ |  $N_{\bar{u}}$ |  $0.1^{+0.7}_{-0.1}$
$\beta_{u}$ |  $-0.5^{+0.2}_{-0.2}$ |  $N_{\bar{d}}$ |  $0.0^{+0.9}_{-0.2}$
$\beta_{d}$ |  $-0.97^{+0.12}_{-0.02}$ |  $\epsilon_{u}$ |  $10^{+7}_{-2}$
$\beta_{\rm sea}$ |  $0.4^{+2.2}_{-0.8}$ |  $\epsilon_{d}$ |  $113^{+215}_{-82}$
Table 5: The $\chi^{2}$ values for different datasets. $N$ is the number of data points for each experimental dataset. Dataset |  $N$ |  $\chi^{2}/N$
---|---|---
COMPASS [33] |  10 |  $1.08_{-0.41}^{+0.52}$
COMPASS [34] |  10 |  $1.29_{-0.44}^{+0.61}$
HERMES [32] $\pi^{+}$ |  26 |  $1.95_{-0.48}^{+0.48}$
HERMES [32] $\pi^{-}$ |  26 |  $1.83_{-0.37}^{+0.46}$
HERMES [32] $K^{+}$ |  26 |  $2.23_{-0.48}^{+0.63}$
HERMES [32] $K^{-}$ |  26 |  $2.35_{-0.48}^{+0.49}$
JLab [36][37] |  6 |  $1.03_{-0.46}^{+0.75}$
total |  130 |  $1.90_{-0.18}^{+0.18}$

Figure 2: Comparison with HERMES SSA data [32] from the proton target for
$\pi^{+}$ (upper left), $\pi^{-}$ (upper right), $K^{+}$ (lower left), and
$K^{-}$ (lower right) productions. The filled data points are included in the
fit, while the open data points are not. The green lines are the central value
calculated from the fit and the bands represent the 1 standard deviation of
the calculated asymmetries by using 100 replicas. Figure 3: Comparison with
COMPASS SSA data [34] from the proton target for $\pi^{+}$, $\pi^{-}$,
$K^{+}$, and $K^{-}$ productions. The markers and bands have the same meaning
as in Fig. 2. Figure 4: Comparison with COMPASS SSA data [33] from the
deuteron target for $\pi^{+}$, $\pi^{-}$, $K^{+}$, $K^{-}$, and $K_{S}^{0}$
productions. The markers and bands have the same meaning as in Fig. 2.

Figure 5: Comparison with JLab SSA data [36, 37] from the ${}^{3}\rm He$
target. The asymmetry values for $\pi^{+}$ and $\pi^{-}$ productions in the
left panel have been effectively converted to those of the neutron, and the
asymmetry values for $K^{+}$ and $K^{-}$ productions. The markers and bands
have the same meaning as in Fig. 2.

###  III.2 EicC projections of the Sivers function

The SIDIS events are generated at the vertex level using a Monte Carlo event
generator [58], which has been adopted in the simulation of several JLab12
SIDIS experiments. It allows us to run in collider mode. According to the
conceptual design of EicC, the electron beam energy is set at $3.5\,\rm GeV$,
the proton beam energy is set at $20\,\rm GeV$, and the 3He beam energy is set
at $40\,\rm GeV$. The unpolarized SIDIS differential cross section used in the
generator is based on global fit to HERMES and COMPASS multiplicity data. A
similar cross section parametrization is also adopted in SIDIS-RC EvGen [59],
a recently released generator for radiative correction studies in the SIDIS
process. For acceptance reasons, we require the electron momentum greater than
$0.35\,\rm GeV$ and the hadron momentum greater than $0.30\,\rm GeV$. The
pseudorapidity $\eta$ is restricted to $-3.5<\eta<3.5$ for both the electron
and the hadron. Full azimuthal angle coverage in the laboratory frame is
assumed. To have full flavor separation of all light quarks, we include both
charged pion, $\pi^{\pm}$, and charge kaon, $K^{\pm}$, productions in SIDIS.
For the particle identification (PID) of the final electron, pion, and kaon,
we apply a PID acceptance cut on the maximum particle momentum, $p_{\rm max}$,
in different pseudorapidity regions,

$\eta$ |  $[-3.5,-1]$ |  $(-1,1]$ |  $(1,3.5]$
---|---|---|---
$p_{\rm max}$ | $4\,\rm GeV$ | $6\,\rm GeV$ | $15\,\rm GeV$

where the initial ion beam is defined along the positive direction. We further
impose the physical cuts, $Q^{2}>1\,\rm GeV^{2}$, $0.3<z<0.7$, $W>5\,\rm GeV$,
and $W^{\prime}>2\,\rm GeV$, to select events in the deep inelastic region and
to exclude resonance regions.

We estimate the statistics by assuming $50\,\rm fb^{-1}$ for $ep$ collisions
and $50\,\rm fb^{-1}$ for $e{{}^{3}}\rm He$ collisions. Such accumulated
luminosities can be achieved with about one year running according to the
designed instantaneous luminosity $2\times 10^{33}\,\rm cm^{-2}s^{-1}$.
Keeping the statistical uncertainty at $10^{-3}$ level, we obtain 13545 data
points in four-dimensional bins in $x$, $Q^{2}$, $z$, and $P_{h\perp}$. Not
only the precision of the EicC pseudodata is much higher than existing world
data, but also the amount of data points is about 10 times more. It allows us
to apply more strict kinematic cuts for a cleaner selection of data in the TMD
region. In this study, we choose $\delta<0.3$, and 4983 EicC pseudodata points
are selected. The distributions of the EicC pseudodata are shown in Fig. 6,
where the colored points are selected in the fit while the gray ones are
excluded. Since the energies per nucleon are different for the proton beam and
the ${}^{3}\rm He$ ion beam, the kinematics of the proton data and the neutron
data are slightly different but still overlap in a wide range.

Figure 6: Kinematic distributions of the EicC pseudodata in $x-Q^{2}$ (left)
and $z-P_{h\perp}$ (right) planes. Each bin is plotted as a point at the bin
center kinematic values. Kinematics cuts, $Q^{2}>1\,\rm GeV^{2}$, $W>5\,\rm
GeV$, $W^{\prime}>2\,\rm GeV$, and $0.3<z<0.7$ have been applied. The blue
points are the proton data with $\delta<0.3$, the red points are the neutron
data with $\delta<0.3$, and the gray points are the data with $\delta>0.3$.

The SSA values of the pseudodata are calculated with the central value of the
result from the fit to the world data. Since a realistic estimation of
systematic uncertainties is only possible when the detailed designs of
detectors are available, we only consider some expected dominant sources of
systematic uncertainties. For the proton data, we assign $3\%$ relative
uncertainty to account for the polarization of the proton beam, and for the
neutron data, we assign $5\%$ relative uncertainty to account for the
polarization of the ${}^{3}\rm He$ ion beam and the nuclear effect. Total
uncertainties are evaluated via the quadrature combination of statistical
uncertainties and systematic uncertainties. The precise EicC data with wide
kinematics coverage allows us to adopt a more flexible parametrization of the
Sivers functions. Therefore, we remove the artificial assumptions in Eq. (59),
while still keep
$\epsilon_{\bar{u}}=\epsilon_{\bar{d}}=\epsilon_{s}=\epsilon_{\bar{s}}=0$, and
then we have 26 free parameters, as listed in Table 6. To estimate the impact
of the EicC on the extraction of the Sivers function, we perform a
simultaneous fit to the world data and the EicC pseudodata as described above.
Following the same procedure, 100 replicas are created by randomly shifting
the values according to the simulated statistical uncertainty. The fit reaches
$\chi^{2}/N=1.15$ for only statistical uncertainties and $\chi^{2}/N=1.13$ for
both statistical and systematic uncertainties. The average values of the
parameters and their uncertainties are provided in Table 6. The results of the
EicC projection of the Sivers functions are shown in Fig. 7 via slices at
various $x$ values, in Fig. 8 via the truncated zeroth transverse momentum
moment, and in Fig. 9 via the truncated first transverse momentum moment in
comparison with the results of the fit to existing world data.

Table 6: The parameters from the fit to world SIDIS data and EicC pseudodata. The central values are the average of the results from 100 fits, and the uncertainties are the standard deviations. The values of $r_{u}$, $r_{d}$, $r_{s}$, $r_{\bar{u}}$, $r_{\bar{d}}$, and $r_{\bar{s}}$ are provided in unit of $\rm GeV^{2}$ and the others are dimensionless. The “Stat.” column means that only statistical uncertainties of EicC pseudodata are considered in the fit, while “Stat. + Syst.” column means that both statistical and systematic uncertainties of EicC pseudodata are included in the fit.  Parameter |  Stat. |  Stat.+Syst. |  Parameter |  Stat. |  Stat.+Syst.
---|---|---|---|---|---
$r_{u}$ |  $0.068^{+0.002}_{-0.001}$ |  $0.067^{+0.002}_{-0.002}$ |  $N_{u}$ |  $-0.075^{+0.001}_{-0.001}$ |  $-0.075^{+0.001}_{-0.001}$
$r_{d}$ |  $0.092^{+0.003}_{-0.003}$ |  $0.091^{+0.003}_{-0.003}$ |  $N_{d}$ |  $0.72^{+0.02}_{-0.02}$ |  $0.72^{+0.02}_{-0.02}$
$r_{s}$ |  $0.005^{+0.044}_{-0.005}$ |  $0.009^{+0.054}_{-0.009}$ |  $N_{s}$ |  $-0.001^{+0.001}_{-0.001}$ |  $-0.001^{+0.001}_{-0.001}$
$r_{\bar{u}}$ |  $0.065^{+0.011}_{-0.008}$ |  $0.064^{+0.012}_{-0.009}$ |  $N_{\bar{u}}$ |  $0.012^{+0.001}_{-0.001}$ |  $0.012^{+0.001}_{-0.001}$
$r_{\bar{d}}$ |  $0.044^{+0.008}_{-0.006}$ |  $0.044^{+0.007}_{-0.007}$ |  $N_{\bar{d}}$ |  $-0.016^{+0.001}_{-0.001}$ |  $-0.016^{+0.001}_{-0.001}$
$r_{\bar{s}}$ |  $1.1^{+6.7}_{-1.0}\times 10^{-8}$ |  $1.5^{+12.2}_{-1.4}\times 10^{-8}$ |  $N_{\bar{s}}$ |  $0.002^{+0.001}_{-0.001}$ |  $0.001^{+0.001}_{-0.001}$
$\beta_{u}$ |  $-0.44^{+0.02}_{-0.03}$ |  $-0.44^{+0.02}_{-0.03}$ |  $\alpha_{u}$ |  $2.57^{+0.04}_{-0.04}$ |  $2.55^{+0.07}_{-0.05}$
$\beta_{d}$ |  $-0.9840^{+0.0003}_{-0.0003}$ |  $-0.9840^{+0.0003}_{-0.0004}$ |  $\alpha_{d}$ |  $2.52^{+0.06}_{-0.05}$ |  $2.56^{+0.08}_{-0.06}$
$\beta_{s}$ |  $0.2^{+3.8}_{-0.6}$ |  $0.2^{+4.1}_{-0.6}$ |  $\alpha_{s}$ |  $6^{+4}_{-2}$ |  $6^{+4}_{-2}$
$\beta_{\bar{u}}$ |  $-0.37^{+0.04}_{-0.04}$ |  $-0.38^{+0.05}_{-0.05}$ |  $\alpha_{\bar{u}}$ |  $4.60^{+0.36}_{-0.27}$ |  $4.53^{+0.39}_{-0.39}$
$\beta_{\bar{d}}$ |  $0.17^{+0.08}_{-0.07}$ |  $0.16^{+0.10}_{-0.08}$ |  $\alpha_{\bar{d}}$ |  $1.003^{+0.008}_{-0.002}$ |  $1.004^{+0.006}_{-0.003}$
$\beta_{\bar{s}}$ |  $-0.7^{+0.1}_{-0.1}$ |  $-0.7^{+0.2}_{-0.1}$ |  $\alpha_{\bar{s}}$ |  $7^{+3}_{-1}$ |  $7^{+3}_{-2}$
$\epsilon_{u}$ |  $5.5^{+0.9}_{-0.5}$ |  $5.5^{+1.2}_{-0.6}$ |  $\epsilon_{d}$ |  $24^{+11}_{-11}$ |  $25^{+13}_{-12}$

Figure 7: The transverse momentum distribution of the Sivers functions at
different $x$ values. The green bands represent the uncertainties of the fit
to world SIDIS data, the red bands represent the EicC projections with only
statistical uncertainties, and the blue bands represent the EicC projections
including systematic uncertainties as described in the text. Figure 8: The
zeroth transverse moment of the Sivers functions as defined in Eq. (60) with
the ${\bf k}_{\perp}$ integral truncated at $0.6\,\rm GeV$. The green bands
represent the uncertainties of the fit to world SIDIS data, the red bands
represent the EicC projections with only statistical uncertainties, and the
blue bands represent the EicC projections including systematic uncertainties
as described in the text. Figure 9: The first transverse moment of the Sivers
functions as defined in Eq. (61) with the ${\bf k}_{\perp}$ integral truncated
at $0.6\,\rm GeV$. The green bands represent the uncertainties of the fit to
world SIDIS data, the red bands represent the EicC projections with only
statistical uncertainties, and the blue bands represent the EicC projections
including systematic uncertainties as described in the text. Figure 10: The
comparisons between our result and the result from Ref. [39] for Sivers
functions of $u,d,\bar{u},\bar{d}$ quark in the momentum space at $x=0.1$ and
$Q=2$ GeV. Within uncertainty, the results are consistent with each other.

##  IV Summary

In this paper, we present a quantitative assessment of the impact of EicC
SIDIS program on the determination of TMDs. Taking the Sivers function as an
example, we perform a global fit of the Sivers asymmetry data in SIDIS at
small transverse momentum, including the TMD evolution at the NNLL accuracy.
The impact of EicC is studied by adding the EicC pseudodata. In this study,
both statistical uncertainties and dominant systematic uncertainties are taken
into account for the EicC pseudodata, while complete detailed systematic
uncertainty studies are left for the future when the detector design is ready.
It has been demonstrated that the Sivers functions can be precisely determined
for various quark flavors, and particularly the sea quark distributions,
including the strange and antistrange, can be extracted at high precision with
the future EicC SIDIS data.

Once the precise data are available from EicC, one will be able to have less
biased extractions of the Sivers functions by using much more flexible
parametrizations. Besides, in the EicC era, one can have a cleaner selection
of data for TMD studies, e.g., by applying a more strict requirement on
$\delta\equiv|P_{h\perp}|/(zQ)$ to restrict data in the low transverse
momentum region and higher $W$ and $W^{\prime}$ cuts to avoid the resonance
region. It is important to remark that both polarized electron-proton and
electron-3He data are necessary for a complete flavor separation. To fully
explore the potential of 3He as an effective neutron source, detailed nuclear
effect corrections should be further investigated both experimentally and
theoretically in the future, since it is the dominant source of systematic
uncertainty in reality by using 3He data.

In principle, the EicC enables us to measure all 18 TMDs-related structure
functions in SIDIS via the combination of different electron and ion beam
polarization configurations and the separation of different azimuthal
modulations. The study of the Sivers function as presented in this paper can
be extended to other TMDs. Multi dimensional binning on $x$, $Q^{2}$, $z$ and
$p_{h\perp}$ will be available for the spin asymmetry measurements, and the
coverage of $x$ by EicC can reach down to about $0.005$. Given the existing
fixed-target experiments covering the low-$Q^{2}$ and high-$x$ region and the
Electron-Ion Collider to be built at BNL in US (US-EIC) reaching much lower
$x$ values, EicC will fill the kinematics gap between the coverage between the
JLab-12GeV program and the US-EIC. Combining the measurements at all these
facilities, we will be able to have a complete physical picture of the three-
dimensional structures of the nucleon with systematically controllable
uncertainties. Therefore, EicC will play an important role in the
understanding of nucleon spin structures with its unique significance for sea
quarks.

###### Acknowledgements.

This work is supported in part by the Strategic Priority Research Program of
the Chinese Academy of Sciences under grant number XDB34000000, the Guangdong
Provincial Key Laboratory of Nuclear Science with No. 2019B121203010, and the
National Natural Science Foundation of China under Contract No. 12175117. P.
Sun is supported by the Natural Science Foundation of China under grant No.
11975127 and No. 12061131006. The authors acknowledge the computing resource
available at the Southern Nuclear Science Computing Center.

##  Appendix A EXPRESSION OF MATCHING FUNCTION

For TMD PDFs, the coefficient function $C$ up to NLO is [55]

$\displaystyle C_{f\leftarrow
f^{\prime}}(x,b,\mu)=\delta(1-x)\delta_{ff^{\prime}}$
$\displaystyle\quad\quad\quad\quad+a_{s}{(\mu)}\Big{(}-\textbf{L}_{\mu}P^{(1)}_{f\leftarrow
f^{\prime}}+C^{(1,0)}_{f\leftarrow f^{\prime}}\Big{)}\ ,$ (62)

where

$\displaystyle C^{(1,0)}_{q\leftarrow q^{\prime}}(x)$
$\displaystyle=C_{F}\Big{[}2(1-x)-\delta(1-x)\frac{\pi^{2}}{6}\Big{]}\delta_{qq^{\prime}},$
(63) $\displaystyle C^{(1,0)}_{q\leftarrow g}(x)$ $\displaystyle=2x(1-x),$
(64) $\displaystyle P^{(1)}_{q\leftarrow q^{\prime}}(x)$
$\displaystyle=2C_{F}\Big{[}\frac{2}{(1-x)_{+}}-1-x+\frac{3}{2}\delta(1-x)\Big{]}\delta_{qq^{\prime}},$
(65) $\displaystyle P^{(1)}_{q\leftarrow g}(x)$ $\displaystyle=1-2x+2x^{2}.$
(66)

For TMD FFs, the matching coefficient $\mathbb{C}$ up to NLO follows the same
pattern as in Eq. (A) with the replacement of the PDF DGLAP kernels
$P^{(1)}_{f\leftarrow f^{\prime}}(x)$ by the FF DGLAP kernels [60],

$\displaystyle{\mathbb{P}}^{(1)}_{q\to q^{\prime}}(z)=$
$\displaystyle\frac{2C_{F}}{z^{2}}\Big{(}\frac{1+z^{2}}{1-z}\Big{)}_{+}\delta_{qq^{\prime}},$
(67) $\displaystyle{\mathbb{P}}^{(1)}_{q\to g}(z)=$
$\displaystyle\frac{2C_{F}}{z^{2}}\frac{1+(1-z)^{2}}{z},$ (68)

and the replacement of $C^{(1,0)}_{f\leftarrow f^{\prime}}(x)$ by [55]

$\displaystyle\mathbb{C}^{(1,0)}_{q\to q^{\prime}}(z)=$
$\displaystyle\frac{C_{F}}{z^{2}}\Big{[}2(1-z)+\frac{4(1+z^{2})\ln{z}}{1-z}$
$\displaystyle-\delta(1-z)\frac{\pi^{2}}{6}\Big{]}\delta_{qq^{\prime}},$ (69)
$\displaystyle\mathbb{C}^{(1,0)}_{q\to g}(z)=$
$\displaystyle\frac{2C_{F}}{z^{2}}\Big{[}z+2\big{(}1+(1-z)^{2}\big{)}\frac{\ln{z}}{z}\Big{]}.$
(70)

The “+” prescription is defined as

$\displaystyle\int_{x_{0}}^{1}dx\,[g(x)]_{+}f(x)$
$\displaystyle=\int_{0}^{1}dx\,g(x)[f(x)\Theta(x-x_{0})-f(1)],$ (71)

where $\Theta(x-x_{0})$ is the Heaviside step function.

##  Appendix B ESTIMATION OF UNCERTAINTIES BY USING REPLICAS

We randomly shift the central values of the data points by Gaussian
distributions with the Gaussian widths to be the experimental uncertainties.
Repeating the process, we get 100 datasets
$y=\\{y_{1},y_{2},\cdots,y_{i},\cdots,y_{100}\\}$. Then each dataset is fitted
to yield the best estimation of various parameters, and hence a replica is
created. Afterward, one can calculate the set
$C=\\{\chi_{1}^{2}/N,\chi_{2}^{2}/N,\cdots,\chi_{i}^{2}/N,\cdots,\chi_{100}^{2}/N\\}$,
where $N$ is the total number of experimental data points. The central value,
$\chi^{2}/N$, is the mean value of the set $C$. In principle, $C$ could be any
physical quantity, for example, the Sivers function or the asymmetry
$A_{UT}^{\sin(\phi_{h}-\phi_{S})}$. In order to calculate the upper and lower
uncertainties of $\chi^{2}/N$, we define $C=\\{C_{up},C_{lo}\\}$, where
$C_{up}=\\{\chi_{1}^{2}/N,\chi_{2}^{2}/N,\cdots,\chi_{i}^{2}/N,\cdots,\chi_{k}^{2}/N\\}$
for $\chi_{i}^{2}/N>\chi^{2}/N$ and
$C_{lo}=\\{\chi_{1}^{2}/N,\chi_{2}^{2}/N,\cdots,\chi_{i}^{2}/N,\cdots,\chi_{100-k}^{2}/N\\}$
for $\chi_{i}^{2}/N<\chi^{2}/N$. Then we can calculate the upper and lower
uncertainties by

$\displaystyle\delta^{+}$
$\displaystyle=\sqrt{\sum_{1}^{k}\frac{(\chi^{2}_{i}/N-\chi^{2}/N)^{2}}{k}}$
(72)

for $\chi^{2}_{i}/N$ in $C_{up}$, and

$\displaystyle\delta^{-}$
$\displaystyle=\sqrt{\sum_{1}^{100-k}\frac{(\chi^{2}_{i}/N-\chi^{2}/N)^{2}}{100-k}}$
(73)

for $\chi^{2}_{i}/N$ in $C_{lo}$.

## References

  * [1] See e.g., M. Constantinou et al., Lattice QCD Calculations of Parton Physics, arXiv:2202.07193 [hep-lat].
  * [2] X. Ji, Parton Physics on a Euclidean Lattice, Phys. Rev. Lett. 110, 262002 (2013).
  * [3] Y.-Q. Ma and J.-W. Qiu, Extracting parton distribution functions from lattice QCD calculations, Phys. Rev. D 98, no.7, 074021 (2018).
  * [4] A. V. Radyushkin, Quasi-parton distribution functions, momentum distributions, and pseudo-parton distribution functions, Phys. Rev. D 96, 034025 (2017).
  * [5] J. Ashman et al. (European Muon Collaboration), A measurement of the spin asymmetry and determination of the structure function $g_{1}$ in deep inelastic muon-proton scattering, Phys. Lett. B 206, 364 (1988).
  * [6] J. Ashman et al. (European Muon Collaboration), An investigation of the spin structure of the proton in deep inelastic scattering of polarised muons on polarised protons, Nucl. Phys. B 328, 1 (1989).
  * [7] C. A. Aidala, S. D. Bass, D. Hasch, and G. K. Mallot, The spin structure of the nucleon, Rev. Mod. Phys. 85, 655 (2013).
  * [8] X. Ji, F. Yuan, and Y. Zhao, What we know and what we don’t know about the proton spin after 30 years, Nature Rev. Phys. 3, 27 (2021).
  * [9] R. L. Jaffe and A. Manohar, The $g_{1}$ problem: Deep inelastic electron scattering and the spin of the proton, Nucl. Phys. B 337, 509 (1990).
  * [10] X. D. Ji, Gauge-Invariant Decomposition of Nucleon Spin, Phys. Rev. Lett. 78, 610 (1997).
  * [11] E. Leader and C. Lorcé, The angular momentum controversy: What’s it all about and does it matter?, Phys. Rep. 541, 163 (2014) [Erratum: 802, 23 (2019)].
  * [12] M. Wakamatsu, Is gauge-invariant complete decomposition of the nucleon spin possible?, Int. J. Mod. Phys. A 29, 1430012 (2014).
  * [13] E. S. Ageev et al. (COMPASS Collaboration), Measurement of the spin structure of the deuteron in the DIS region, Phys. Lett. B 612, 154 (2005).
  * [14] V. Y. Alexakhin et al. (COMPASS Collaboration), The Deuteron Spin-dependent Structure Function g1(d) and its First Moment, Phys. Lett. B 647, 8 (2007).
  * [15] A. Airapetian et al. (HERMES Collaboration), Precise determination of the spin structure function g(1) of the proton, deuteron and neutron, Phys. Rev. D 75, 012007 (2007).
  * [16] E. R. Nocera et al. (NNPDF Collaboration), A first unbiased global determination of polarized PDFs and their uncertainties, Nucl. Phys. B 887, 276 (2014).
  * [17] D. de Florian, R. Sassot, M. Stratmann and W. Vogelsang, Evidence for polarization of gluons in the proton, Phys. Rev. Lett. 113, no.1, 012001 (2014).
  * [18] J. J. Ethier, N. Sato, and W. Melnitchouk, First simultaneous extraction of spin-dependent parton distributions and fragmentation functions from a global QCD analysis, Phys. Rev. Lett. 119, no.13, 132001 (2017).
  * [19] Y. Zhou, N. Sato, and W. Melnitchouk, How well do we know the gluon polarization in the proton?, Phys. Rev. D 105, no.7, 074022 (2022).
  * [20] J. C. Collins, D. E. Soper, and G. F. Sterman, Factorization of Hard Processes in QCD, Adv. Ser. Dir. High Energy Phys. 5, 1 (1989).
  * [21] J. Collins, Foundations of perturbative QCD, Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology Series Vol. 32 (Cambridge University Press, Cambridge, England, 2011).
  * [22] E. D. Bloom et al., High-Energy Inelastic e p Scattering at 6-Degrees and 10-Degrees, Phys. Rev. Lett. 23, 930(1969).
  * [23] D. W. Sivers, Single Spin Production Asymmetries from the Hard Scattering of Point-Like Constituents, Phys. Rev. D 41, 83 (1990).
  * [24] D. W. Sivers, Hard scattering scaling laws for single spin production asymmetries, Phys. Rev. D 43, 261 (1991).
  * [25] D. Boer and P. J. Mulders, Time reversal odd distribution functions in leptoproduction, Phys. Rev. D 57, 5780 (1998).
  * [26] J. C. Collins, Fragmentation of transversely polarized quarks probed in transverse momentum distributions, Nucl. Phys. B 396, 161 (1993).
  * [27] S. J. Brodsky, D. S. Hwang, and I. Schmidt, Final state interactions and single spin asymmetries in semiinclusive deep inelastic scattering, Phys. Lett. B 530, 99 (2002).
  * [28] S. J. Brodsky, D. S. Hwang, and I. Schmidt, Initial state interactions and single spin asymmetries in Drell-Yan processes, Nucl. Phys. B 642, 344 (2002).
  * [29] J. C. Collins, Leading twist single transverse-spin asymmetries: Drell-Yan and deep inelastic scattering, Phys. Lett. B 536, 43 (2002).
  * [30] A. Airapetian et al. (HERMES Collaboration), Observation of the Naive-T-odd Sivers Effect in Deep-Inelastic Scattering, Phys. Rev. Lett. 103, 152002 (2009).
  * [31] A. Airapetian et al. (HERMES Collaboration), Effects of transversity in deep-inelastic scattering by polarized protons, Phys. Lett. B 693, 11 (2010).
  * [32] A. Airapetian et al. (HERMES Collaboration), Azimuthal single- and double-spin asymmetries in semi-inclusive deep-inelastic lepton scattering by transversely polarized protons, J. High Energy Phys. 12 (2020) 010.
  * [33] M. Alekseev et al. (COMPASS Collaboration), Collins and Sivers asymmetries for pions and kaons in muon-deuteron DIS, Phys. Lett. B 673, 127 (2009).
  * [34] C. Adolph et al. (COMPASS Collaboration), Collins and Sivers asymmetries in muonproduction of pions and kaons off transversely polarised protons, Phys. Lett. B 744, 250 (2015).
  * [35] C. Adolph et al. (COMPASS Collaboration), Sivers asymmetry extracted in SIDIS at the hard scales of the Drell–Yan process at COMPASS, Phys. Lett. B 770, 138 (2017).
  * [36] X. Qian et al. (Jefferson Lab Hall A Collaboration), Single Spin Asymmetries in Charged Pion Production from Semi-Inclusive Deep Inelastic Scattering on a Transversely Polarized 3He Target, Phys. Rev. Lett. 107, 072003 (2011).
  * [37] Y. X. Zhao et al. (Jefferson Lab Hall A Collaboration), Single spin asymmetries in charged kaon production from semi-inclusive deep inelastic scattering on a transversely polarized ${}^{3}He$ target, Phys. Rev. C 90, no.5, 055201 (2014).
  * [38] M. G. Echevarria, Z.-B. Kang, and J. Terry, Global analysis of the Sivers functions at NLO+NNLL in QCD, J. High Energy Phys. 01 (2021) 126.
  * [39] M. Bury, A. Prokudin and A. Vladimirov, Extraction of the Sivers function from SIDIS, Drell-Yan, and $W^{\pm}/Z$ boson production data with TMD evolution, J. High Energy Phys. 05 (2021) 151.
  * [40] M. Aghasyan et al. (COMPASS Collaboration), First measurement of transverse-spin-dependent azimuthal asymmetries in the Drell-Yan process, Phys. Rev. Lett. 119, no.11, 112002 (2017).
  * [41] L. Adamczyk et al. (STAR Collaboration), Measurement of the transverse single-spin asymmetry in $p^{\uparrow}+p\to W^{\pm}/Z^{0}$ at RHIC, Phys. Rev. Lett. 116, no.13, 132301 (2016).
  * [42] M. Anselmino, M. Boglione, U. D’Alesio, F. Murgia, and A. Prokudin, Study of the sign change of the Sivers function from STAR Collaboration W/Z production data, J. High Energy Phys. 04 (2017) 046.
  * [43] A. Accardi et al., Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all, Eur. Phys. J. A 52, no.9, 268 (2016).
  * [44] R. Abdul Khalek et al., Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report, Nucl. Phys. A 1026, 122447 (2022).
  * [45] D. P. Anderle et al., Electron-ion collider in China, Front. Phys. 16, no.6, 64701 (2021).
  * [46] D. P. Anderle, T.-J. Hou, H. Xing, M. Yan, C. P. Yuan, and Y. Zhao, Determining the helicity structure of the nucleon at the Electron Ion Collider in China, J. High Energy Phys. 08 (2021) 034.
  * [47] A. Bacchetta, U. D’Alesio, M. Diehl, and C. A. Miller, Single-spin asymmetries: The Trento conventions, Phys. Rev. D 70, 117504 (2004).
  * [48] J. C. Collins, D. E. Soper, and G. F. Sterman, Transverse Momentum Distribution in Drell-Yan Pair and W and Z Boson Production, Nucl. Phys. B 250, 199 (1985).
  * [49] M. G. Echevarria, A. Idilbi, A. Schäfer, and I. Scimemi, Model-Independent Evolution of Transverse Momentum Dependent Distribution Functions (TMDs) at NNLL, Eur. Phys. J. C 73, no.12, 2636 (2013).
  * [50] V. Bertone, I. Scimemi, and A. Vladimirov, Extraction of unpolarized quark transverse momentum dependent parton distributions from Drell-Yan/Z-boson production, J. High Energy Phys. 06 (2019) 028.
  * [51] S. Tafat, Nonperturbative corrections to the Drell-Yan transverse momentum distribution, J. High Energy Phys. 05 (2001) 004.
  * [52] A. A. Vladimirov, Self-contained definition of the Collins-Soper kernel, Phys. Rev. Lett. 125, no.19, 192002 (2020).
  * [53] F. Hautmann, I. Scimemi, and A. Vladimirov, Non-perturbative contributions to vector-boson transverse momentum spectra in hadronic collisions, Phys. Lett. B 806, 135478 (2020).
  * [54] J. Collins and T. Rogers, Understanding the large-distance behavior of transverse-momentum-dependent parton densities and the Collins-Soper evolution kernel, Phys. Rev. D 91, no.7, 074020 (2015).
  * [55] I. Scimemi and A. Vladimirov, Non-perturbative structure of semi-inclusive deep-inelastic and Drell-Yan scattering at small transverse momentum, J. High Energy Phys. 06 (2020) 137.
  * [56] R. B. Wiringa, V. G. J. Stoks and R. Schiavilla, Accurate nucleon-nucleon potential with charge-independence breaking, Phys. Rev. C 51, 38 (1995).
  * [57] J. L. Friar, B. F. Gibson, G. L. Payne, A. M. Bernstein, and T. E. Chupp, Neutron polarization in polarized He-3 targets, Phys. Rev. C 42, 2310 (1990).
  * [58] https://github.com/TianboLiu/LiuSIDIS.
  * [59] https://github.com/duanebyer/sidis.
  * [60] M. Stratmann and W. Vogelsang, Next-to-leading order evolution of polarized and unpolarized fragmentation functions, Nucl. Phys. B 496, 41 (1997).