# Extraction of trans-helicity worm-gear distributions and opportunities at
the Electron-Ion Collider in China

Ke Yang<EMAIL_ADDRESS>School of Physics, Peking University,
Beijing 100871, 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  Southern Center for Nuclear-Science Theory (SCNT), Institute of Modern
Physics, Chinese Academy of Sciences, Huizhou 516000, China  Peng Sun
<EMAIL_ADDRESS>Institute of Modern Physics, Chinese Academy of Sciences,
Lanzhou, Gansu Province 730000, China  University of Chinese Academy of
Sciences, Beijing 100049, China  Yuxiang Zhao<EMAIL_ADDRESS>Southern
Center for Nuclear-Science Theory (SCNT), Institute of Modern Physics, Chinese
Academy of Sciences, Huizhou 516000, China  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  Bo-Qiang Ma<EMAIL_ADDRESS>School of
Physics, Peking University, Beijing 100871, China  Center for High Energy
Physics, Peking University, Beijing 100871, China  Collaborative Innovation
Center of Quantum Matter, Beijing, China

###### Abstract

We present a global analysis of the trans-helicity worm-gear distribution
function, $g_{1T}^{\perp}$, by fitting the longitudinal-transverse double spin
asymmetry data of the semi-inclusive deep inelastic scattering. The analysis
is performed within the framework of transverse momentum dependent
factorization and evolution. It is found that the $u$-quark favors a positive
distribution and the $d$-quark favors a negative distribution, which is
consistent with previous model calculations and phenomenological extractions.
Based on the fit to existing world data, we also study the impact of the
proposed electron-ion collider in China and conclude that it can significantly
improve the precision of the worm-gear distribution function and hence enhance
our understanding of nucleon spin structures.

##  I Introduction

Understanding the internal structure of nucleons is pivotal for comprehending
the strong force that binds quarks and gluons within nucleons, and for
shedding light on the fundamental properties of the matter. In recent years,
the pursuit of multi-dimensional tomography of the nucleon has emerged as a
cutting-edge approach to probe distributions of quarks and gluons within the
nucleon, offering a deeper understanding of its internal dynamics. Transverse
momentum dependent (TMD) parton distribution functions (PDFs) contain the
information of the parton transverse momentum with respect to the parent
nucleon, and hence provide three-dimensional imaging of the nucleon in the
momentum space.

At the leading twist, there are eight TMDs for quarks [1, 2, 3, 4]. Among
them, the worm-gear-T distribution $g_{1T}^{\perp}(x,k_{T}^{2})$, also known
as the trans-helicity distribution [5, 6] or the Kotzinian-Mulders function
[7, 8], describes the probability density of finding a longitudinally
polarized quark with longitudinal momentum fraction $x$ and transverse
momentum $k_{T}$ in a transversely polarized nucleon. As well as the worm-
gear-L, or longi-transversity, distribution $h_{1L}^{\perp}$ that describes
the probability density of finding a transversely polarized quark in a
longitudinally polarized nucleon, it can be expressed as the overlap between
wave functions differing by one unit of orbital angular momentum [9, 10, 11,
12, 13], and many efforts have been devoted to the worm-gear TMDs to
understand nucleon spin and flavor structures.

Although the two worm-gear distributions are defined as independent quantities
from the decomposition of the quark-quark correlator, some relation, such as
$g_{1T}^{\perp}=-h_{1L}^{\perp}$, is suggested based on quark model like
calculations [14, 11, 15, 16, 5, 6]. Following the SU(6) spin-flavor
structure, the $g_{1T}^{\perp}$ distribution of the up quark was predicted to
be positive and with a greater magnitude than the negative down quark
distribution, and explicit calculations have been done in the light-cone
constituent quark model [17, 18, 11, 19, 20], the spectator diquark model [20,
14, 21, 10], the MIT bag model [16], and the covariant parton model [15]. On
the other hand, the large-$N_{c}$ approximation [22] states that the worm-gear
distributions of up quark and down quark only differ by a sign and have the
same magnitude, i.e. $g_{1T}^{\perp u}=-g_{1T}^{\perp d}$. Besides, if taking
the Wandzura-Wilczek (WW)-type approximation [23, 24, 25, 26, 7, 8], which
neglects the contribution from quark-gluon-quark correlations, one may relate
the trans-helicity worm-gear distribution to the helicity distribution as

$g_{1T}^{\perp(1)}(x)\stackrel{{\scriptstyle\text{WW}}}{{\approx}}x\int_{x}^{1}\frac{dy}{y}g_{1}(y),$
(1)

where

$g_{1T}^{\perp(1)}(x)\equiv\pi\int
d^{2}{k}_{T}\frac{k_{T}^{2}}{2M^{2}}g_{1T}^{\perp}(x,{k}_{T}^{2}),$ (2)

is the first transverse moment that has also been studied in lattice QCD [27,
28, 29].

In experiment, the semi-inclusive deep inelastic scattering (SIDIS) is one of
the main processes to study TMDs. According to the TMD factorization, the
trans-helicity distribution $g_{1T}^{\perp}$ contributes to a double spin
asymmetry $A_{LT}$ with azimuthal modulation as $\cos(\phi_{h}-\phi_{S})$.
With the development of polarized beams and targets, this asymmetry has been
measured by HERMES [30], COMPASS [31, 32, 33] and Jefferson Lab (JLab) [34].
In some recent phenomenological analyses [35, 36], it was found that the
extracted worm-gear distributions supported the positive result for the up
quark and the negative result for the down quark as suggested by the model
calculations. However, due to the limited accuracy of existing world data, one
has to introduce some bias in the fit to obtain reasonable results and almost
no constraint is put on sea quarks.

The Electron-Ion Collider in China (EicC) is proposed as a future facility in
nuclear physics, and one of its main physics goals is to precisely measure
nucleon TMDs via the SIDIS process. It is designed to deliver a $3.5\,\rm GeV$
electron beam with 80% polarization colliding with various types of ion beams.
The designed energy of the proton beam is $20\,\rm GeV$ and correspondingly
the energy of the 3He beam is $40\,\rm GeV$. Both the proton and the 3He beams
can be longitudinally or transversely polarized with 70% polarization. The
instantaneous luminosity can reach about $2\times
10^{33}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$. The EicC kinematic coverage will
fill the gap between multi-hall SIDIS program at the $12\,\rm GeV$ upgraded
JLab, which covers relatively large-$x$ region dominated by valence quarks,
and the Electron-Ion Collider (EIC) to be built at the Brookhaven National
Laboratory (BNL), which can reach the small-$x$ region down to about $10^{-4}$
[37, 38]. Therefore, a combination of all these facilities is expected to
provide precise determination of TMDs in a full kinematic coverage [39],
towards a complete three-dimensional imaging of nucleon spin structures.

In this paper, we perform a global analysis of trans-helicity TMDs by fitting
the longitudinal-transverse double spin asymmetry data from HERMES, COMPASS
and JLab. Taking the world data fit result as the baseline, we further study
the impact of the EicC SIDIS program on the determination of the worm-gear
distribution. The rest of the paper is organized as follows. In Sec. II, we
briefly review the theoretical framework. In Sec. III, we present the
parametrization of the trans-helicity worm-gear distributions and the fit
results to world data. In Sec. IV, we study the EicC impact on the extraction
of the worm-gear distributions by adding simulated pseudodata in the fit. A
summary is drawn in Sec. V.

##  II Theoretical formalism

We consider the SIDIS process

$\ell(l)+N(P)\longrightarrow\ell(l^{\prime})+h(P_{h})+X,$ (3)

where $\ell$ represents the lepton, $N$ represents the nucleon, and $h$
represents the detected hadron. The four-momenta of corresponding particles
are given in parentheses. The commonly used kinematic variables for the SIDIS
process are defined as

$\displaystyle Q^{2}=-(l-l^{\prime})^{2}=-q^{2},$ (4) $\displaystyle
x=\frac{Q^{2}}{2P\cdot q},\quad y=\frac{P\cdot q}{P\cdot l},\quad
z=\frac{P\cdot P_{h}}{P\cdot q},$ (5)
$\displaystyle\gamma=\frac{2xM}{Q}=\frac{MQ}{P\cdot q},$ (6)

where $q=l-l^{\prime}$ is the transferred momentum and $M$ is the nucleon
mass.

For the SIDIS process with a transversely polarized target and a
longitudinally polarized lepton beam, one can write the differential cross
section within the one-photon-exchange approximation as [40]

$\displaystyle\frac{d\sigma}{dxdydzd\phi_{h}d\phi_{S}dP_{hT}^{2}}=\sigma_{0}\Big{\\{}F_{UU}$
(7)
$\displaystyle\quad+\lambda_{e}|S_{\perp}|\Big{[}\sqrt{1-\varepsilon^{2}}\cos(\phi_{h}-\phi_{S})F_{LT}^{\cos(\phi_{h}-\phi_{S})}$
$\displaystyle\quad+\sqrt{2\varepsilon(1-\varepsilon)}\cos(2\phi_{h}-\phi_{S})F_{LT}^{\cos(2\phi_{h}-\phi_{S})}$
$\displaystyle\quad+\sqrt{2\varepsilon(1-\varepsilon)}\cos(\phi_{S})F_{LT}^{\cos\phi_{S}}+\cdots\Big{]}\Big{\\}},$

where

$\sigma_{0}=\frac{\alpha^{2}}{xyQ^{2}}\frac{y^{2}}{2(1-\varepsilon)}\left(1+\frac{\gamma^{2}}{2x}\right),$
(8)

$\alpha$ is the electromagnetic fine structure constant, $|S_{\perp}|$
represents the transversal component of the nucleon spin vector, $\lambda_{e}$
represents the helicity of the lepton beam, and $\varepsilon$ is the ratio of
longitudinal and transverse photon flux,

$\varepsilon=\frac{1-y-\frac{1}{4}\gamma^{2}y^{2}}{1-y+\frac{1}{2}y^{2}+\frac{1}{4}\gamma^{2}y^{2}}.$
(9)

Figure 1: Trento conventions of the transverse momentum and azimuthal angles.

As shown in Fig. 1, we follow the Trento conventions [41], in which the
momenta of the virtual photon and the nucleon are chosen along the $\hat{z}$
direction. One can express the transverse momentum $P_{hT}$ of the hadron and
the azimuthal angles $\phi_{h}$ and $\phi_{S}$ in Lorentz invariant forms as

$\displaystyle P_{hT}=\sqrt{-g^{\mu\nu}_{\bot}P_{h\mu}P_{h\nu}},$ (10)
$\displaystyle\cos\phi_{h}=-\frac{l_{\mu}P_{h\nu}g^{\mu\nu}_{\bot}}{l_{\bot}P_{hT}},\quad\sin\phi_{h}=-\frac{l_{\mu}P_{h\nu}\epsilon^{\mu\nu}_{\bot}}{l_{\bot}P_{hT}},$
(11)
$\displaystyle\cos\phi_{S}=-\frac{l_{\mu}S_{\bot\nu}g^{\mu\nu}_{\bot}}{l_{\bot}S_{\bot}},\quad\sin\phi_{S}=-\frac{l_{\mu}S_{\bot\nu}\epsilon^{\mu\nu}_{\bot}}{l_{\bot}S_{\bot}},$
(12)

where $l_{\bot}=\sqrt{-g^{\mu\nu}_{\bot}l_{\mu}l_{\nu}}$ and
$S_{\bot}=\sqrt{-g^{\mu\nu}_{\bot}S_{\mu}S_{\nu}}$ with $S^{\mu}$ being the
spin vector of the nucleon. The transverse metric and the transverse
antisymmetry tensor are defined as

$\displaystyle
g^{\mu\nu}_{\bot}=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),$
(13)
$\displaystyle\epsilon^{\mu\nu}_{\bot}=\epsilon^{\mu\nu\rho\sigma}\frac{P_{\rho}q_{\sigma}}{P\cdot
q\sqrt{1+\gamma^{2}}},$ (14)

where $\epsilon^{\mu\nu\rho\sigma}$ is the totally antisymmetric tensor with
the convention $\epsilon^{0123}=1$.

The worm-gear distribution $g_{1T}^{\perp}$ can be extracted from the
longitudinal-transverse double spin asymmetry, which is given by the ratio
between the structure functions $F_{LT}^{\cos(\phi_{h}-\phi_{S})}$ and
$F_{UU}$. According to the TMD factorization [40], the structure functions at
low transverse momentum, i.e. small $\delta=|P_{hT}|/(zQ)$, can be
approximated in terms of TMD PDF and TMD fragmentation function (FF) as

$\displaystyle
F_{UU}=\big{|}C_{V}\big{(}Q^{2},\mu\big{)}\big{|}^{2}x\sum_{q}e_{q}^{2}\int_{0}^{\infty}\frac{b_{T}db_{T}}{2\pi}J_{0}\left(\frac{b_{T}P_{hT}}{z}\right)$
$\displaystyle\quad\times f_{1,q\leftarrow
H}(x,b_{T};\mu,\zeta)D_{1,q\rightarrow
h}(z,b_{T};\mu,\bar{\zeta})+\mathcal{O}\left(\frac{P_{hT}^{2}}{Q^{2}}\right),$
(15) $\displaystyle
F_{LT}^{\cos\left(\phi_{h}-\phi_{S}\right)}=\left|C_{V}\left(Q^{2},\mu\right)\right|^{2}x\sum_{q}e_{q}^{2}M\int_{0}^{\infty}\frac{b_{T}^{2}db_{T}}{2\pi}$
$\displaystyle\quad\times
J_{1}\left(\frac{b_{T}P_{hT}}{z}\right)g_{1T,q\leftarrow
H}^{\perp}(x,b_{T};\mu,\zeta)D_{1,q\rightarrow h}(z,b_{T};\mu,\bar{\zeta})$
$\displaystyle\quad+\mathcal{O}\left(\frac{P_{hT}^{2}}{Q^{2}}\right),$ (16)

where $e_{q}$ is the electric charge of the quark with flavor $q$, $C_{V}$ is
the hard factor that can be calculated via perturbative QCD, and $J_{0}$ and
$J_{1}$ are the first kind Bessel functions. Here the unpolarized TMD PDF
$f_{1}$, the worm-gear TMD PDF $g_{1T}^{\perp}$, and the unpolarized TMD FF
$D_{1}$ are given in $b_{T}$ space. They are related to corresponding
functions in the transverse momentum space through Fourier transforms,

$f_{1}(x,k_{T};\mu,\zeta)=\int_{0}^{\infty}\frac{b_{T}db_{T}}{2\pi}J_{0}(b_{T}k_{T})f_{1}(x,b_{T};\mu,\zeta),$
(17) $\displaystyle\frac{k_{T}}{M}g_{1T}^{\perp}$
$\displaystyle(x,k_{T};\mu,\zeta)$ (18)
$\displaystyle=\int_{0}^{\infty}\frac{b_{T}^{2}db_{T}}{2\pi}MJ_{1}(b_{T}k_{T})g_{1T}^{\perp}(x,b_{T};\mu,\zeta),$
$D_{1}(z,p_{T};\mu,\zeta)=\int_{0}^{\infty}\frac{b_{T}db_{T}}{2\pi}J_{0}(b_{T}p_{T})D_{1}(z,b_{T};\mu,\zeta),$
(19)

where $k_{T}$ represents the quark transverse momentum with respect to the
nucleon and $p_{T}$ represents the quark transverse momentum with respect to
the produced hadron. The details of the Fourier transformation are given in
Appendix A.

###  II.1 Evolution of TMD PDFs and FFs

The energy scale dependence on $\mu$ and $\zeta$ of the TMD functions are
given by the evolution equations,

$\displaystyle\mu^{2}\frac{d}{d\mu^{2}}F(x,b_{T};\mu,\zeta)=\frac{\gamma_{F}(\mu,\zeta)}{2}F(x,b_{T};\mu,\zeta),$
(20)
$\displaystyle\zeta\frac{d}{d\zeta}F(x,b_{T};\mu,\zeta)=-\mathcal{D}(\mu,b_{T})F(x,b_{T};\mu,\zeta),$
(21)

where $\gamma_{F}$ is the anomalous dimension, and $\mathcal{D}$ is the
rapidity anomalous dimension (RAD), also known as the Collins-Soper kernel.
The $F$ represents some TMD PDF or TMD FF, i.e. $f_{1}$, $D_{1}$, and
$g_{1T}^{\perp}$ in this study. One may have the formal solution,

$\displaystyle F$ $\displaystyle(x,b_{T};\mu,\zeta)$ (22)
$\displaystyle=R[\left(b_{T};\mu_{i},\zeta_{i}\right)\rightarrow(b_{T};\mu,\zeta)]F\left(x,b_{T};\mu_{i},\zeta_{i}\right),$

which relates the TMD PDF (or FF) at $(\mu,\zeta)$ to that at the initial
point $(\mu_{i},\zeta_{i})$. The evolution factor
$R\left[\left(b_{T};\mu_{i},\zeta_{i}\right)\rightarrow(b_{T};\mu,\zeta)\right]$
can be expressed as

$\displaystyle R[(b_{T};\mu_{i},\zeta_{i})\rightarrow(b_{T};\mu,\zeta)]$ (23)
$\displaystyle=\exp\left[\int_{\mathcal{P}}\left(\frac{\gamma_{F}(\mu,\zeta)}{\mu}d\mu-\frac{\mathcal{D}(\mu,b_{T})}{\zeta}d\zeta\right)\right],$

where $\mathcal{P}$ represents the path connecting the scales
$(\mu_{i},\zeta_{i})$ and $(\mu,\zeta)$. As a common choice, we set the energy
scales as $\mu^{2}=\zeta=Q^{2}$.

According to the integrability condition [43]

$\zeta\frac{d}{d\zeta}\gamma_{F}(\mu,\zeta)=-\mu\frac{d}{d\mu}\mathcal{D}(\mu,b_{T})=-\Gamma_{\text{cusp
}}(\mu),$ (24)

the evolution factor $R[(b_{T};\mu_{i},\zeta_{i})\rightarrow(b_{T};Q,Q^{2})]$
is in principle path independent. However, it differs from path to path when
truncating at some fixed order in perturbation theory. As suggested in Ref.
[42], the condition (24) allows one to construct a two-dimensional field
${\cal F}(\mu,\zeta)$, of which the gradient is given by
$\mathbf{E}=(\gamma_{F}/2,-\mathcal{D})$. Then $F(x,b_{T},\mu,\zeta)$ remains
unchanged if the path is along the equipotential line of ${\bf E}$, referred
to as a null-evolution line. In the $(\mu,\zeta)$ plane, there is a unique
saddle point $(\mu_{0},\zeta_{0})$ defined by

$\mathcal{D}(\mu_{0},b_{T})=0,\quad\gamma_{F}(\mu_{0},\zeta_{0})=0.$ (25)

Among the null-evolution lines, only the one passing through the saddle point
has finite $\zeta$ at all values of $\mu$. Hence, the $F(x,b_{T})\equiv
F(x,b_{T};\mu_{0},\zeta_{0})$ is referred to as the optimal TMD PDF or FF
[42]. Owing to the good properties of the null-evolution line and the saddle
point, we firstly evolve the $F$ from the saddle point along the null-
evolution line to the point with $\mu=Q$; secondly, we evolve the $F$ along
the straight line keeping $\mu=Q$ fixed until reaching the point with
$\zeta=Q^{2}$. The result for the evolution factor
$R[(b_{T};\mu_{i},\zeta_{i})\rightarrow(b_{T};Q,Q^{2})]$ along this path is
[42]

$R\left[\left(b_{T};\mu_{i},\zeta_{i}\right)\rightarrow\left(b_{T};Q,Q^{2}\right)\right]=\left(\frac{Q^{2}}{\zeta_{\mu}(Q,b_{T})}\right)^{-\mathcal{D}(Q,b_{T})}.$
(26)

The expressions for $\mathcal{D}(Q,b_{T})$ and $\zeta_{\mu}(Q,b_{T})$ can be
found in Appendix B. The precision for the perturbative calculation of various
factors in powers of $\alpha_{s}$ in this work is summarized in Table 1.

Table 1: Orders of perturbative calculations for anomalous dimensions and the $C(\mathbb{C})$ functions in the optimal TMD PDF and FF. | $\Gamma_{\rm cusp}$ | $\gamma_{V}$ | $\mathcal{D}_{\rm resum}$ | $\zeta_{\mu}^{\rm pert}$ | $\zeta_{\mu}^{\rm exact}$
---|---|---|---|---|---
$R$ | $\alpha_{s}^{3}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{1}$ | $\alpha_{s}^{1}$
| $f_{1}$ | $D_{1}$ | $g_{1T}^{\perp}$ |  |
$C(\mathbb{C})$ | $\alpha_{s}^{1}$ | $\alpha_{s}^{1}$ | $\alpha_{s}^{0}$ |  |

###  II.2 Unpolarized TMD PDF and FF

For unpolarized TMD PDFs and FFs, we adopt the SV19 parametrization [42]. The
optimal unpolarized TMD PDF and FF are expressed as

$\displaystyle f_{1,f\leftarrow
h}\left(x,b_{T}\right)=\sum_{f^{\prime}}\int_{x}^{1}\frac{dy}{y}C_{f\leftarrow
f^{\prime}}\left(y,b_{T},\mu_{\mathrm{OPE}}^{\mathrm{PDF}}\right)$
$\displaystyle\quad\quad\times f_{1,f^{\prime}\leftarrow
h}\left(\frac{x}{y},\mu_{\mathrm{OPE}}^{\mathrm{PDF}}\right)f_{\mathrm{NP}}(x,b_{T}),$
(27) $\displaystyle D_{1,f\rightarrow
h}\left(z,b_{T}\right)=\frac{1}{z^{2}}\sum_{f^{\prime}}\int_{z}\frac{dy}{y}y^{2}\mathbb{C}_{f\rightarrow
f^{\prime}}\left(y,b_{T},\mu_{\mathrm{OPE}}^{\mathrm{FF}}\right)$
$\displaystyle\quad\quad\times d_{1,f^{\prime}\rightarrow
h}\left(\frac{z}{y},\mu_{\mathrm{OPE}}^{\mathrm{FF}}\right)D_{\mathrm{NP}}(z,b_{T}),$
(28)

where $f_{1,f^{\prime}\leftarrow h}$ and $d_{1,f^{\prime}\rightarrow h}$ are
collinear PDFs and FFs. The scales $\mu_{\mathrm{OPE}}^{\mathrm{PDE}}$ and
$\mu_{\mathrm{OPE}}^{\mathrm{FF}}$ are chosen as

$\displaystyle\mu_{\mathrm{OPE}}^{\mathrm{PDE}}=\frac{2e^{-\gamma_{E}}}{b_{T}}+2\,\mathrm{GeV},$
(29)
$\displaystyle\mu_{\mathrm{OPE}}^{\mathrm{FF}}=\frac{2e^{-\gamma_{E}}z}{b_{T}}+2\,\mathrm{GeV},$
(30)

where $\gamma_{E}$ is the Euler-Mascheroni constant. The $2\,\mathrm{GeV}$
shift is introduced to keep the PDFs and FFs in perturbative region when
$b_{T}$ is large. The nonperturbative functions $f_{\rm NP}(x,b_{T})$ and
$D_{\rm NP}(z,b_{T})$ are to be parametrized.

For unpolarized TMD PDF, the coefficient function $C$ can be written as

$\displaystyle C_{f\leftarrow f^{\prime}}(x,b_{T},\mu)=$
$\displaystyle\delta(1-x)\delta_{ff^{\prime}}+$ (31) $\displaystyle
a_{s}(\mu)\left(-\mathbf{L}_{\mu}P_{f\leftarrow
f^{\prime}}^{(1)}+C_{f\leftarrow f^{\prime}}^{(1,0)}\right),$

up to NLO, where $a_{s}=\frac{g^{2}(\mu)}{(4\pi)^{2}}$ and $g(\mu)$ are the
QCD coupling constants. $\mathbf{L}_{\mu}$ is defined as

$\mathbf{L}_{\mu}=\ln\left(\frac{b_{T}^{2}\mu^{2}}{4e^{-2\gamma_{E}}}\right).$
(32)

$P_{f\leftarrow f^{\prime}}^{(1)}$ is the coefficient of the PDF evolution
kernel, which reads

$\displaystyle P_{q\leftarrow
q^{\prime}}^{(1)}(x)=2C_{F}\left(\frac{1+x^{2}}{1-x}\right)_{+}\delta_{qq^{\prime}},$
(33) $\displaystyle P_{q\leftarrow g}^{(1)}(x)=1-2x+2x^{2}.$ (34)

The “$+$” prescription is defined as

$\int_{x_{0}}^{1}dx[g(x)]_{+}f(x)=\int_{0}^{1}dxg(x)\left[f(x)\Theta\left(x-x_{0}\right)-f(1)\right],$
(35)

where $\Theta\left(x-x_{0}\right)$ is the Heaviside step function. $C_{F}=4/3$
is the quadratic Casimir eigenvalues of fundamental representation of $SU(3)$.
The expressions of $C_{f\leftarrow f^{\prime}}^{(n,0)}$ can be found in Ref.
[44], and their NLO terms are

$\displaystyle C_{q\leftarrow q^{\prime}}^{(1,0)}(x)$
$\displaystyle=C_{F}\left(2\bar{x}-\delta(\bar{x})\frac{\pi^{2}}{6}\right)\delta_{qq^{\prime}},$
(36) $\displaystyle C_{q\leftarrow g}^{(1,0)}(x)$ $\displaystyle=2x\bar{x},$
(37)

and $\bar{x}=1-x$.

For unpolarized TMD FF, one should replace $C_{f\leftarrow
f^{\prime}}(x,b_{T},\mu)$, $P_{f\leftarrow f^{\prime}}^{(1)}$ and
$C_{f\leftarrow f^{\prime}}^{(n,0)}(x)$ by $\mathbb{C}_{f\rightarrow
f^{\prime}}\left(z,b_{T},\mu\right)$, $\mathbb{P}_{f\rightarrow
f^{\prime}}^{(1)}$ and $\mathbb{C}_{f\rightarrow f^{\prime}}^{(n,0)}(z)$,
which at NLO are expressed as

$\displaystyle\mathbb{P}_{q\rightarrow q^{\prime}}^{(1)}(z)$
$\displaystyle=\frac{2C_{F}}{z^{2}}\bigg{(}\frac{1+z^{2}}{1-z}\bigg{)}_{+}\delta_{qq^{\prime}},$
(38) $\displaystyle\mathbb{P}_{q\rightarrow g}^{(1)}(z)$
$\displaystyle=\frac{2C_{F}}{z^{2}}\frac{1+(1-z)^{2}}{z},$ (39)
$\displaystyle\mathbb{C}_{q\rightarrow q^{\prime}}^{(1,0)}(z)$
$\displaystyle=\frac{C_{F}}{z^{2}}\bigg{[}2(1-z)+\frac{4(1+z^{2})\ln z}{1-z}$
(40)
$\displaystyle\quad-\delta(1-z)\frac{\pi^{2}}{6}\bigg{]}\delta_{qq^{\prime}},$
(41) $\displaystyle\mathbb{C}_{q\rightarrow g}^{(1,0)}(z)$
$\displaystyle=\frac{2C_{F}}{z^{2}}\bigg{[}z+2\big{(}1+(1-z)^{2}\big{)}\frac{\ln
z}{z}\bigg{]}.$ (42)

###  II.3 Worm-gear asymmetry

The longitudinal-transverse double spin asymmetry of the SIDIS process is
defined as

$\displaystyle A_{LT}=$
$\displaystyle\frac{1}{\left|S_{\perp}\right|\left|\lambda_{e}\right|}\frac{\left[\mathrm{d}\sigma_{LT}(+,\uparrow)-\mathrm{d}\sigma_{LT}(-,\uparrow)\right]-\left[\mathrm{d}\sigma_{LT}(+,\downarrow)-\mathrm{d}\sigma_{LT}(-,\downarrow)\right]}{\mathrm{d}\sigma_{LT}(+,\uparrow)+\mathrm{d}\sigma_{LT}(-,\uparrow)+\mathrm{d}\sigma_{LT}(+,\downarrow)+\mathrm{d}\sigma_{LT}(-,\downarrow)},$
(43)

where $+$ ($-$) represents the positive (negative) helicity state of the
electron beam and $\uparrow$ ($\downarrow$) represents the transverse spin
direction of nucleon $S_{\bot}$ to be parallel (anti-parallel) to the
designated positive transverse axis. The worm-gear asymmetry is defined as
$\cos\left(\phi_{h}-\phi_{S}\right)$ modulation of the double spin asymmetry,

$A_{LT}^{\cos\left(\phi_{h}-\phi_{S}\right)}=\frac{\left\langle
2\cos\left(\phi_{h}-\phi_{S}\right)\sigma_{\mathrm{LT}}\right\rangle}{\sqrt{1-\varepsilon^{2}}\left\langle\sigma_{UU}\right\rangle}=\frac{F_{LT}^{\cos\left(\phi_{h}-\phi_{S}\right)}}{F_{UU}}.$
(44)

and one can express it with functions we defined above as

$A_{LT}^{\cos\left(\phi_{h}-\phi_{S}\right)}=\frac{M\sum_{q}e_{q}^{2}\int_{0}^{\infty}\frac{b_{T}^{2}db_{T}}{2\pi}J_{1}\left(\frac{b_{T}P_{hT}}{z}\right)R^{2}\left[\left(b_{T};\mu_{i},\zeta_{i}\right)\rightarrow\left(b_{T};Q,Q^{2}\right)\right]g_{1T,q\leftarrow
N}^{\perp}(x,b_{T})D_{1,q\rightarrow
h}(z,b_{T})}{\sum_{q}e_{q}^{2}\int_{0}^{\infty}\frac{b_{T}db_{T}}{2\pi}J_{0}\left(\frac{b_{T}P_{hT}}{z}\right)R^{2}\left[\left(b_{T};\mu_{i},\zeta_{i}\right)\rightarrow\left(b_{T};Q,Q^{2}\right)\right]f_{1,q\leftarrow
N}(x,b_{T})D_{1,q\rightarrow h}(z,b_{T})},$ (45)

where $N$ is the target and $h$ is the detected hadron.

##  III Extraction of the worm-gear distributions

###  III.1 Fit to world SIDIS data

With the formalism above, we perform a global analysis of world SIDIS data to
extract the worm-gear distributions $g_{1T}^{\perp}$ of the nucleon. The
results will also serve as the baseline for the impact study of the EicC.

We parametrize the optimal worm-gear distributions of the proton at the
initial scale as

$\displaystyle g_{1T,q\leftarrow
p}^{\perp}(x,b_{T})=N_{q}\frac{(1-x)^{\alpha_{q}}x^{\beta_{q}}}{\mathrm{~{}B}(\alpha+1,\beta+1)}\exp\left(-r_{q}b_{T}^{2}\right)$
(46)

for $u$ and $d$ quarks and

$\displaystyle g_{1T,q\leftarrow
p}^{\perp}(x,b_{T})=N_{q}f_{1}(x,\mu_{0})\exp\left(-r_{q}b_{T}^{2}\right)$
(47)

for $\bar{u}$, $\bar{d}$, $s$, and $\bar{s}$ quarks with $\mu_{0}=2\,\rm GeV$.
Here $\mathrm{~{}B}(\alpha+1,\beta+1)$ is the Euler Beta function, introduced
to reduce the correlation among parameters. Assuming the isospin symmetry, we
can express corresponding distribution functions of the neutron as

$\displaystyle g_{1T,u\leftarrow n}^{\perp}(x,b_{T})$
$\displaystyle=g_{1T,d\leftarrow P}^{\perp}(x,b_{T}),$ (48) $\displaystyle
g_{1T,\bar{u}\leftarrow n}^{\perp}(x,b_{T})$
$\displaystyle=g_{1T,,\bar{d}\leftarrow P}^{\perp}(x,b_{T}),$ (49)
$\displaystyle g_{1T,d\leftarrow n}^{\perp}(x,b_{T})$
$\displaystyle=g_{1T,u\leftarrow P}^{\perp}(x,b_{T}),$ (50) $\displaystyle
g_{1T,\bar{d}\leftarrow n}^{\perp}(x,b_{T})$
$\displaystyle=g_{1T,\bar{u}\leftarrow P}^{\perp}(x,b_{T}),$ (51)
$\displaystyle g_{1T,s\leftarrow n}^{\perp}(x,b_{T})$
$\displaystyle=g_{1T,s\leftarrow P}^{\perp}(x,b_{T}),$ (52) $\displaystyle
g_{1T,\bar{s}\leftarrow n}^{\perp}(x,b_{T})$
$\displaystyle=g_{1T,\bar{s}\leftarrow P}^{\perp}(x,b_{T}).$ (53)

For unpolarized TMD PDFs and FFs, we adopt the SV19 fit [42], in which the
nonperturbative functions $f_{\mathrm{NP}}$ and $D_{\mathrm{NP}}$ are
parametrized as

$\displaystyle
f_{\mathrm{NP}}(x,b_{T})=\exp\left[-\frac{\lambda_{1}(1-x)+\lambda_{2}x+x(1-x)\lambda_{5}}{\sqrt{1+\lambda_{3}x^{\lambda_{4}}b_{T}^{2}}}b_{T}^{2}\right],$
(54) $\displaystyle
D_{\mathrm{NP}}(z,b_{T})=\exp\left[-\frac{\eta_{1}z+\eta_{2}(1-z)}{\sqrt{1+\eta_{3}(b_{T}/z)^{2}}}\frac{b_{T}^{2}}{z^{2}}\right]\left(1+\eta_{4}\frac{b_{T}^{2}}{z^{2}}\right).$
(55)

The values of the parameters $\lambda_{i}$ and $\eta_{i}$ are listed in Table
2, which can be also found in [42, 44]. For the FFs to charged hadrons, we
approximate them as

$\displaystyle D_{1,f\leftarrow h^{+}}$
$\displaystyle=D_{1,f\leftarrow\pi^{+}}+D_{1,f\leftarrow
K^{+}}+D_{1,f\leftarrow p},$ (56) $\displaystyle D_{1,f\leftarrow h^{-}}$
$\displaystyle=D_{1,f\leftarrow\pi^{-}}+D_{1,f\leftarrow
K^{-}}+D_{1,f\leftarrow\bar{p}}.$ (57)

Table 2: The parameters for nonperturbative functions of the optimal unpolarized TMD PDF and FF. The units are in $\mathrm{GeV}^{2}$ except that $\lambda_{4}$ is dimensionless.  $\lambda_{1}$ |  | $\lambda_{2}$ |  | $\lambda_{3}$ |  | $\lambda_{4}$ |  | $\lambda_{5}$
---|---|---|---|---|---|---|---|---
$0.198$ |  | $9.3$ |  | $431$ |  | $2.12$ |  | $-4.44$
$\eta_{1}$ |  | $\eta_{2}$ |  | $\eta_{3}$ |  | $\eta_{4}$ |  |
$0.260$ |  | $0.476$ |  | $0.478$ |  | $0.483$ |  |
Table 3: The SIDIS double spin asymmetry data by HERMES [30], COMPASS [31, 32, 33], and JLab [34]. The SFA refers to $\left\langle 2\cos\left(\phi_{h}-\phi_{S}\right)\sigma_{LT}\right\rangle/\left(\sqrt{1-\varepsilon^{2}}\left\langle\sigma_{UU}\right\rangle\right)$, and the CSA refers to $\left\langle 2\cos\left(\phi_{h}-\phi_{S}\right)\sigma_{LT}\right\rangle/\left\langle\sigma_{UU}\right\rangle$. Data set | Target | Beam | original data | data points | data points | Process | Measurement
---|---|---|---|---|---|---|---
|  |  | points | after cut | after cut |  |
|  |  |  | $\delta<0.5,$ | $\delta<0.3,$ |  |
|  |  |  | $Q>1\,\mathrm{GeV}$ | $Q>1\,\mathrm{GeV}$ |  |
HERMES [30] | H2 | $27.6\,\mathrm{GeV}\,e^{\pm}$ | 64 | 26 | 11 | $e^{\pm}p\rightarrow e^{\pm}\pi^{+}X$ | SFA
|  |  | 64 | 26 | 11 | $e^{\pm}p\rightarrow e^{\pm}\pi^{-}X$ | SFA
|  |  | 64 | 26 | 12 | $e^{\pm}p\rightarrow e^{\pm}K^{+}X$ | SFA
|  |  | 64 | 26 | 12 | $e^{\pm}p\rightarrow e^{\pm}K^{-}X$ | SFA
|  |  | 64 | 30 | 15 | $e^{\pm}p\rightarrow e^{\pm}PX$ | SFA
COMPASS [31] | NH3 | $160\,\mathrm{GeV}\,\mu^{+}$ | 66 | 28 | 9 | $\mu^{+}p\rightarrow\mu^{+}h^{+}X$ | SFA
|  |  | 66 | 26 | 8 | $\mu^{+}p\rightarrow\mu^{+}h^{-}X$ | SFA
JLab [34] |  3He | $5.9\,\mathrm{GeV}\,e^{-}$ | 4 | 2 | 1 | $e^{-}n\rightarrow e^{-}\pi^{+}X$ | CSA
|  |  | 4 | 2 | 1 | $e^{-}n\rightarrow e^{-}\pi^{-}X$ | CSA
Total |  |  | 460 | 192 | 80 |  |

In this analysis, we include the SIDIS longitudinal-transverse double spin
asymmetry data from HERMES [30], COMPASS [31, 32, 33] and JLab [34], as
summarized in Table 3. Since the TMD factorization is only valid at small
$\delta=|P_{hT}|/(zQ)$, only data with $\delta<0.5$ are included in the fit.

For the HERMES data, experimental results are provided in both one-dimensional
binning and three-dimensional binning. We only use the three-dimensional bins
in this study, because they are supposed to contain more information for the
study of TMDs, which are multidimensional functions. For the COMPASS data, the
experimental results are provided in one-dimensional binning but on $x$, $z$,
and $P_{hT}$ respectively. Since they tell the dependence on different
variable, we include all these bins in the fit. However, to avoid double
counting, we multiply a factor of $1/3$ when calculating the $\chi^{2}$ from
the COMPASS data. Then the total $\chi^{2}/N$ to be minimized in the fit is
defined as

$\chi^{2}/N=\frac{\frac{1}{3}\chi^{2}_{\mathrm{COMPASS}}+\chi^{2}_{\mathrm{HERMES}}+\chi^{2}_{\mathrm{JLab}}}{\frac{1}{3}N_{\mathrm{COMPASS}}+N_{\mathrm{HERMES}}+N_{\mathrm{JLab}}},$
(58)

where $N_{\rm data~{}set}$ represents the number of points for each data set.
For each data set, we have

$\chi^{2}_{\text{data set}}=\sum_{\begin{subarray}{c}i,j\in\\\ \text{data-
points
}\end{subarray}}\left(t_{i}-a_{i}\right)V_{ij}^{-1}\left(t_{j}-a_{j}\right),$
(59)

where $i$ and $j$ run over all points in each set, $t_{i}$ represent the
theoretical values, and $a_{i}$ represent experimental values. The $V$-matrix
is given by

$V_{ij}=\delta_{ij}\left(\sigma_{i}^{\rm uncor.}\right)^{2}+\sigma_{i}^{\rm
cor.}\sigma_{j}^{\rm cor.},$ (60)

where $\sigma_{i}^{\text{ uncor. }}$ and $\sigma_{i}^{\rm cor.}$ stand for
uncorrelated and correlated uncertainties respectively.

As the existing world data are not precise enough to constrain all parameters
introduced in Eq. (46), we practically reduce the number of parameters by
imposing the conditions,

$\alpha_{u}=\alpha_{d}=\alpha,\quad\beta_{u}=\beta_{d}=\beta,\quad
r_{u}=r_{d}=r,$ (61)

and assuming vanishing distributions for sea quarks, $\bar{u}$, $\bar{d}$,
$s$, and $\bar{s}$. In the end, we have five free parameters to be determined
as listed in Table 4.

To estimate the uncertainties, we create 1000 replicas of the data by smearing
the central values of each data point according to a Gaussian distribution
with data uncertainties being the widths. For each replica, we perform a fit.
Then the central values of all physical quantities are evaluated from the
average of the 1000 fits. More details of this approach are described in Ref.
[44].

In this study, we achieve total $\chi^{2}/N=0.84$ as listed in Table 5,
together with $\chi^{2}$ values for each data set. The expectation values and
uncertainties of the parameters are summarized in Table 4. In Figs. 2–6, we
show the comparison between the fit results and experimental data, in which
the filled points are included in fit while the open points are not. The
extracted worm-gear distribution functions $g_{1T}^{\perp}(x,k_{T})$ are shown
in Fig. 7 at several $x$-slices. As one can observe from the results, the $u$
quark distribution is positive, while the $d$ quark favors a negative
distribution though still consistent with zero. This finding qualitatively
agrees with the predictions from the quark model [11, 17, 18, 19, 20]. In
addition, we also evaluate the transverse moments of the worm-gear
distributions,

$\displaystyle g_{1T}^{\perp(0)}(x)$ $\displaystyle=2\pi\int_{0}^{k_{T}^{\rm
max}}k_{T}dk_{T}g_{1T}^{\perp}\left(x,k_{T}\right),$ (62) $\displaystyle
g_{1T}^{\perp(1)}(x)$ $\displaystyle=2\pi\int_{0}^{k_{T}^{\rm
max}}k_{T}dk_{T}\left(\frac{k_{T}^{2}}{2M^{2}}\right)g_{1T}^{\perp}\left(x,k_{T}\right),$
(63)

where the truncation is chosen as $k_{T}^{\rm max}=1\,\rm GeV$. The results
are shown in Figs. 8–11, which are compatible with Ref. [36].

Figure 2: Comparison between the fit results and the experimental data by
COMPASS [31] from the proton target with charged hadron $h^{\pm}$ measured in
the final state. The filled points are within the kinematic cuts, $\delta<0.5$
and $Q>1\,\rm GeV$, and included the fit, while the open points are not
included in the fit. The green lines and the bands are the mean values and the
standard deviations calculated from the fits to 1000 replicas.

Figure 3: Comparison between the fit results and the experimental data by
HERMES [30] from the proton target with $\pi^{\pm}$ measured in the final
state. The data points and the fit bands follow the same notations as Fig. 2.

Figure 4: Comparison between the fit results and experimental data by HERMES
[30] from the proton target with $K^{\pm}$ measured in the final state. The
data points and the fit bands follow the same notations as Fig. 2. Figure 5:
Comparison between the fit results and the experimental data by HERMES [30]
from the proton target with a fragmented proton measured in the final state.
The data points and the fit bands follow the same notation as Fig. 2. Figure
6: Comparison between the fit results and the experimental data by JLab [34]
from the effectively polarized neutron target with $\pi^{\pm}$ measured in the
final state. The data points and the fit bands follow the same notation as
Fig. 2.

Figure 7: The worm-gear distributions $g_{1T}^{\perp}(x,k_{T})$ at the scale
$Q=2\,\rm GeV$. The uncertainty bands correspond to 68% CL estimated from the
fits to 1000 replicas. The green bands are extracted distributions by fitting
the world SIDIS data, the red bands are EicC projections with only statistical
uncertainties, and the blue bands are EicC projections with both statistical
and systematic uncertainties.

Figure 8: The zeroth transverse moment of the worm-gear functions,
$g_{1T}^{\perp(0)}(x)$ as defined in Eq. (62), for $u$ and $d$ quarks at the
scale $Q=2\,\rm GeV$. The uncertainty bands correspond to 68% CL estimated
from the fits to 1000 replicas. The green bands are extracted distributions by
fitting the world SIDIS data, the red bands are EicC projections with only
statistical uncertainties, and the blue bands are EicC projections with both
statistical and systematic uncertainties.

Figure 9: EicC projections of the zeroth transverse moment of the worm-gear functions, $g_{1T}^{\perp(0)}(x)$ as defined in Eq. (62), for $\bar{u}$, $\bar{d}$, $s$, and $\bar{s}$ quarks at the scale $Q=2\,\rm GeV$. The uncertainty bands correspond to 68% CL estimated from the fits to 1000 replicas. The red bands only contain statistical uncertainties, and the blue bands contain both statistical and systematic uncertainties. Table 4: Results of parameters for world data fit. The central values are the average of the results from 1000 replicas, and the uncertainties correspond to $68\%$ CL. The value of $r$ is provided in unit of $\mathrm{GeV}^{2}$ and the others are dimensionless. Parameter | Value | Parameter | Value
---|---|---|---
$N_{u}$ | $0.0206^{+0.0058}_{-0.0050}$ | $\alpha$ | $16.59^{+65.88}_{-10.11}$
$N_{d}$ | $-0.0073^{+0.0079}_{-0.0082}$ | $\beta$ | $5.57^{+28.65}_{-3.87}$
$10^{9}r$ | ${1.52}^{+9.25}_{-1.49}$ |  |
Table 5: The list of results of $\chi^{2}/N$ of world data fit for each world data set. Data set | data points | $\chi^{2}/N$
---|---|---
HERMES $\pi^{+}$ [30] | 26 | $1.01$
HERMES $\pi^{-}$ [30] | 26 | $0.75$
HERMES $K^{+}$ [30] | 26 | $1.09$
HERMES $K^{-}$ [30] | 26 | $0.68$
HERMES $P$ [30] | 30 | $0.90$
COMPASS $h^{+}$ [31] | 28 | $0.44$
COMPASS $h^{-}$ [31] | 26 | $0.65$
JLab $\pi^{+}$ [34] | 2 | $0.61$
JLab $\pi^{-}$ [34] | 2 | $1.15$
Total | 192 | $0.84$

##  IV EicC projections

The EicC events are generated at the vertex level using the SIDIS Monte Carlo
generator, which has been used in previous studies [44, 45]. To select events
in the DIS region, we apply the cuts,

$\displaystyle Q^{2}>1\,\mathrm{GeV},\quad 0.3<z<0.7,$ (64) $\displaystyle
W>5\,\mathrm{GeV},\quad W^{\prime}>2\,\mathrm{GeV},$ (65)

where $W=\sqrt{(q+P)^{2}}$ is the invariant mass of the produced, and
$W^{\prime}=\sqrt{(q+P-P_{h})^{2}}$ is the missing mass. According to the
detection conditions of the designed EicC detector, we further require the
scattered electron momentum $P_{e}>0.35\,\mathrm{GeV}$ and the hadron momentum
$P_{h}>0.3\,\mathrm{GeV}$. In the simulation, we take the $3.5\,\mathrm{GeV}$
polarized electron beam with $80\%$ polarization, the $20\,\mathrm{GeV}$
transversely polarized proton beam with $70\%$ polarization, and the
$40\,\mathrm{GeV}$ transversely polarized ${}^{3}\mathrm{He}$ beam with $70\%$
polarization. Aiming at a complete separation of contributions from all light
flavor quarks, we take into account both $\pi^{\pm}$ and $K^{\pm}$ data.

To quantify the impact, we assume $50\,\mathrm{fb}^{-1}$ integrated
luminosities of $ep$ and $e^{3}\mathrm{He}$ collisions, which can be achieved
with about one-year run according to the proposed instantaneous luminosity.
For the systematic uncertainties, we assign 2% relative uncertainty to the
polarization of the electron beam, 3% relative uncertainty to the polarization
of the ion beam, and 5% relative uncertainty to the 3He nuclear effect. These
are expected the dominant sources of systematic uncertainties based on our
current knowledge from existing polarized SIDIS measurements. Because the
detailed design of the detectors are still unavailable, we leave more
realistic estimation of systematic uncertainties to future studies.

The central values of the worm-gear asymmetry for the EicC pseudodata are
evaluated from world data fit, which only include nonvanishing contributions
from $u$ and $d$ quarks. Owing to the considerable amount of the expected EicC
data, we can adopt stricter criteria to select data in the TMD region. Hence
we set the cut as $\delta<0.3$, and $5008$ pseudodata points are included.
Besides, we also free the parameters $N_{q}$ for sea quark distributions and
choose the same $r_{q}$ for $\bar{u}$, $\bar{d}$, $s$, and $\bar{s}$ as

$r_{s}=r_{\bar{u}}=r_{\bar{d}}=r_{\bar{s}}=(r_{u}+r_{d})/2.$ (66)

Then there are 12 free parameters in our fit summarized in Table 6. Following
the same procedure, we perform a simultaneous fit to the world data and the
EicC pseudodata. This analysis gives $\chi^{2}/N=1.09$ with corresponding
values and uncertainties of the parameters listed in Table 4, which are
evaluated from the fits to 1000 replicas. The EicC projections of the worm-
gear distributions $g_{1T}^{\perp}(x,k_{T})$ are shown in Fig. 7. The zeroth
transverse momentum moments $g_{1T}^{\perp(0)}(x)$ are shown in Figs. 8 and 9,
and the first transverse momentum moments $g_{1T}^{\perp(1)}(x)$ are shown in
Figs 10 and 11.

##  V SUMMARY

In this work, we perform a global fit to the worm-gear asymmetries from SIDIS
in a small transverse momentum region, including the TMD evolution effect at
the next-to-next-to-leading-logarithmic (NNLL) accuracy. Due to the fact that
the existing experimental uncertainties are too large to determine the worm-
gear distributions of sea quarks, only up and down quarks are considered in
the global fit. Then an impact study is performed by including the EicC pseudo
data in our global fit. For EicC pseudo data, the statistical uncertainties
and dominant systematic uncertainties are taken into account. The latter is
mainly due to the uncertainties from beam polarimetry and the uncertainties of
3He nuclear effects.

Once the precise data are available from EicC, the precision of the worm-gear
distributions for up and down quark will be significantly improved. Meanwhile,
it will also provide the opportunity to extract the worm-gear distributions of
sea quarks. With much more expected precise data from EicC, one can extract
the TMDs utilizing more flexible parametrizations and thus less biased
determination of the nucleon spin structures. Owing to the high precision and
a wide phase space coverage of EicC pseudo data, a more strict cut of
$\delta$, $W$ and $W^{\prime}$ will be feasible. It allows us to have clearer
selection of data required by the TMD factorization. On the other hand, the
events in the transition region are also valuable to test the matching between
TMD and collinear regions. The combination of polarized $ep$ and
$e^{3}\mathrm{He}$ data at similar kinematics are essential for a complete
flavor separation. It is important to remark that the kinematics coverage of
EicC will fill the gap between the on-going JLab-$12\,\mathrm{GeV}$ program
and the approved Electron-Ion Collider to be built at BNL. With all these
facilities, we will be able to have a complete physical picture of nucleon
three-dimensional structures, towards a profound understanding of strong
interactions.

Table 6: Results of parameters for EicC pseudo data fit. The central values are the average of the results from 1000 replicas, and the uncertainties correspond to $68\%$ CL. The values of $r_{u}$ and $r_{d}$ are provided in unit of $\mathrm{GeV}^{2}$ and the others are dimensionless. Parameter | Stat. | Stat. + Syst. | Parameter | Stat. | Stat. + Syst.
---|---|---|---|---|---
$N_{u}$ | $0.0209^{+0.0006}_{-0.0006}$ | $0.0209^{+0.0006}_{-0.0006}$ | $\alpha_{u}$ | $12.46^{+0.72}_{-0.67}$ | $12.47^{+0.67}_{-0.64}$
$N_{d}$ | $-0.0077^{+0.0008}_{-0.0009}$ | $-0.0077^{+0.0008}_{-0.0008}$ | $\alpha_{d}$ | $13.01^{+4.68}_{-2.70}$ | $12.94^{+6.22}_{-2.83}$
$N_{s}$ | $-0.00023^{+0.00044}_{-0.00046}$ | $-0.00026^{+0.00042}_{-0.00047}$ | $\beta_{u}$ | $4.46^{+0.26}_{-0.25}$ | $4.46^{+0.23}_{-0.23}$
$N_{\bar{u}}$ | $0.00019^{+0.00023}_{-0.00022}$ | $0.00020^{+0.00024}_{-0.00021}$ | $\beta_{d}$ | $4.31^{+1.47}_{-0.88}$ | $4.29^{+1.80}_{-0.90}$
$N_{\bar{d}}$ | $0.00021^{+0.00032}_{-0.00031}$ | $0.00022^{+0.00034}_{-0.00032}$ | $r_{u}$ | $0.0067^{+0.0050}_{-0.0048}$ | $0.0067^{+0.0053}_{-0.0050}$
$N_{\bar{s}}$ | $0.00038^{+0.00045}_{-0.00037}$ | $0.00038^{+0.00045}_{-0.00037}$ | $r_{d}$ | $0.016^{+0.025}_{-0.016}$ | $0.016^{+0.024}_{-0.016}$

Figure 10: The first transverse moment of the worm-gear functions,
$g_{1T}^{\perp(1)}(x)$ as defined in Eq. (63), for $u$ and $d$ quarks at the
scale $Q=2\,\rm GeV$. The uncertainty bands correspond to 68% CL estimated
from the fits to 1000 replicas. The green bands are extracted distributions by
fitting the world SIDIS data, the red bands are EicC projections with only
statistical uncertainties, and the blue bands are EicC projections with both
statistical and systematic uncertainties.

Figure 11: EicC projections of the first transverse moment of the worm-gear
functions, $g_{1T}^{\perp(1)}(x)$ as defined in Eq. (63), for $\bar{u}$,
$\bar{d}$, $s$, and $\bar{s}$ quarks at the scale $Q=2\,\rm GeV$. The
uncertainty bands correspond to 68% CL estimated from the fits to 1000
replicas. The red bands only contain statistical uncertainties, and the blue
bands contain both statistical and systematic uncertainties.

###### 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
Major Project of Basic and Applied Basic Research No. 2020B0301030008, the
Guangdong Provincial Key Laboratory of Nuclear Science with No.
2019B121203010, the National Natural Science Foundation of China under Grants
No. 12175117 and No. 12321005, and Shandong Province Natural Science
Foundation Grant No. ZFJH202303. P. Sun is supported by the Natural Science
Foundation of China under Grants No. 11975127 and No. 12061131006. B.-Q. Ma is
supported by National Natural Science Foundation of China under Grants No.
12075003 and No. 12335006. The authors also acknowledge the computing resource
available at the Southern Nuclear Science Computing Center.

##  Appendix A Fourier transforms for TMDs

The Fourier transforms for TMDs are

$\displaystyle f_{1}\left(x,k_{T}\right)$
$\displaystyle=\int\frac{d^{2}\mathbf{b}_{T}}{4\pi^{2}}e^{i\mathbf{b}_{T}\cdot\mathbf{k}_{T}}f_{1}(x,b_{T})$
(67)
$\displaystyle=\int_{0}^{+\infty}\frac{b_{T}db_{T}}{2\pi}J_{0}\left(b_{T}k_{T}\right)f_{1}(x,b_{T}),$
$\displaystyle f_{1}(x,b_{T})$ $\displaystyle=\int
d^{2}\mathbf{k}_{T}e^{-i\mathbf{b}_{T}\cdot\mathbf{k}_{T}}f_{1}\left(x,k_{T}\right)$
(68)
$\displaystyle=2\pi\int_{0}^{+\infty}k_{T}dk_{T}J_{0}\left(b_{T}k_{T}\right)f_{1}\left(x,k_{T}\right),$
$\displaystyle\frac{k_{T}}{M}g_{1T}^{\perp}\left(x,k_{T}\right)$
$\displaystyle=\int\frac{d^{2}\mathbf{b}_{T}}{4\pi^{2}}e^{i\mathbf{b}_{T}\cdot\mathbf{k}_{T}}(-i{b}_{T}M)g_{1T}^{\perp}(x,b_{T}),$
(69) $\displaystyle g_{1T}^{\perp}\left(x,k_{T}\right)$
$\displaystyle=\frac{M^{2}}{k_{T}}\int_{0}^{+\infty}\frac{b_{T}^{2}db_{T}}{2\pi}J_{1}\left(b_{T}k_{T}\right)g_{1T}^{\perp}(x,b_{T}),$
$\displaystyle(-iM{b}_{T})$ $\displaystyle g_{1T}^{\perp}(x,b_{T})=\int
d^{2}\mathbf{k}_{T}e^{-i\mathbf{b}_{T}\cdot\mathbf{k}_{T}}\frac{k_{T}}{M}g_{1T}^{\perp}\left(x,k_{T}\right),$
(70) $\displaystyle g_{1T}^{\perp}(x,b_{T})$
$\displaystyle=\frac{2\pi}{M^{2}b_{T}}\int_{0}^{+\infty}k_{T}^{2}dk_{T}J_{1}\left(b_{T}k_{T}\right)g_{1T}^{\perp}\left(x,k_{T}\right),$
$\displaystyle D_{1}\left(z,k_{T}\right)$
$\displaystyle=\int\frac{d^{2}\mathbf{b}_{T}}{4\pi^{2}}e^{-i\mathbf{b}_{T}\cdot\mathbf{k}_{T}}D_{1}(z,b_{T})$
(71)
$\displaystyle=\int_{0}^{+\infty}\frac{b_{T}db_{T}}{2\pi}J_{0}\left(b_{T}k_{T}\right)D_{1}(z,b_{T}),$
$\displaystyle D_{1}(z,b_{T})$ $\displaystyle=\int
d^{2}\mathbf{k}_{T}e^{i\mathbf{b}_{T}\cdot\mathbf{k}_{T}}D_{1}\left(z,k_{T}\right)$
(72)
$\displaystyle=2\pi\int_{0}^{+\infty}k_{T}dk_{T}J_{0}\left(b_{T}k_{T}\right)D_{1}\left(z,k_{T}\right).$

##  Appendix B Expressions for energy evolution factor

The $\mathcal{D}(Q,b_{T})$ is the rapidity anomalous dimension (RAD). At large
values of $b_{T}$, the $\mathcal{D}(Q,b_{T})$ behaves like a linear function
of $b_{T}$, which is suggested by some models such as [46, 47]. Therefore, we
parameterize the RAD as

$\mathcal{D}(\mu,b_{T})=\mathcal{D}_{\text{resum}}\left(\mu,b_{T}^{*}(b_{T})\right)+c_{0}b_{T}b_{T}^{*}(b_{T}),$
(73)

where the $\mathcal{D}_{\text{resum}}\left(\mu,b_{T}^{*}(b_{T})\right)$ is the
resummed perturbative expansion of RAD, and $b_{T}^{*}(b_{T})$ take the form

$b_{T}^{*}(b_{T})=\frac{b_{T}}{\sqrt{1+b_{T}^{2}/B_{\mathrm{NP}}^{2}}}.$ (74)

At small values of $b_{T}$, the term $c_{0}b_{T}b_{T}^{*}(b_{T})$ can be
ignored and the term $\mathcal{D}_{\text{resum}}$ is dominant, while at large
values of $b_{T}$, the $\mathcal{D}$ behave like $c_{0}B_{\text{NP}}b_{T}$. We
take $B_{\text{NP}}=1.93\,\mathrm{GeV}^{-1}$ and
$c_{0}=0.0391\,\mathrm{GeV}^{2}$ as determined in SV19 model [42]. The
$\mathcal{D}_{\text{resum}}$ can be represented as:

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

where $X=\beta_{0}a_{s}\mathbf{L}_{\mu}$, $\beta_{i}$ are coefficients of
anomaly dimension of strong coupling constant, which satisfies

$\mu^{2}\frac{da_{s}(\mu)}{d\mu^{2}}=-\beta\left(a_{s}\right)=-\sum_{i=0}^{\infty}a_{s}^{i+2}(\mu)\beta_{i}.$
(76)

$C_{A}=3$ and $T_{R}=1/2$ are color factors of the $SU(3)$. The $\Gamma_{i}$
are coefficients of expansion of CUSP anomaly dimension
$\Gamma_{\text{cusp}}(\mu)$, which is related with the integrability condition
(24) of the evolution equation. The $\Gamma_{i}$ are defined by

$\Gamma_{\text{cusp }}(\mu)=\sum_{i=0}^{\infty}a_{s}^{i+1}\Gamma_{i}.$ (77)

With the CUSP anomaly dimension, the anomaly dimension $\gamma_{V}$ can be
written as

$\gamma_{F}(\mu,\zeta)=\Gamma_{\text{cusp
}}(\mu)\ln\left(\frac{\mu^{2}}{\zeta}\right)-\gamma_{V}(\mu),$ (78)

and the $\gamma_{V}(\mu)$ can be expanded as

$\gamma_{V}(\mu)=\sum_{i=1}^{\infty}a_{s}^{n}\gamma_{i}.$ (79)

The $\Gamma_{i}$ and $\gamma_{i}$ can be determined perturbatively, and up to
two-loop level

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

where $N_{f}$ is the number of active quark flavors and have different values
at different energy scales (see Table 7), $\zeta_{3}\approx 1.202$ is the
Apéry’s constant.

Table 7: Values of $N_{f}$ at different values of energy scale. $\mu\leq 1.27\,\mathrm{GeV}$ | $N_{f}=3$
---|---
$1.27<\mu\leq 4.18\,\mathrm{GeV}$ | $N_{f}=4$
$\mu>4.18\,\mathrm{GeV}$ | $N_{f}=5$

Due to that the non-perturbative corrections to the RAD can not be ignored at
large-$b_{T}$, we need to use the exact solution of $\zeta_{\mu}$ at
large-$b_{T}$; while at very small-$b_{T}$, we use the perturbative solution.
In order to connect these two region, we introduce a
$e^{-b_{T}^{2}/B_{\mathrm{NP}}^{2}}$ factor, and the $\zeta_{\mu}$ is
expressed as [42]

$\zeta_{\mu}(b_{T})=\zeta_{\mu}^{\text{pert
}}(b_{T})e^{-\frac{b_{T}^{2}}{B_{\mathrm{NP}}^{2}}}+\zeta_{\mu}^{\operatorname{exact}}(b_{T})\left(1-e^{-\frac{b_{T}^{2}}{B_{\mathrm{NP}}^{2}}}\right).$
(81)

Therefore, at $b_{T}^{2}\ll B^{2}_{\mathrm{NP}}$, $\zeta_{\mu}$ is dominantly
given by perturbative solution, and at other regions, it will turn to exact
solution. We express the $\zeta_{\mu}^{\text{pert }}$ and
$\zeta_{\mu}^{\operatorname{exact}}$ here,

$\zeta_{\mu}^{\text{pert }}(\mu,b_{T})=\frac{2\mu
e^{-\gamma_{E}}}{b_{T}}e^{-v(\mu,b_{T})},$ (82)
$\zeta_{\mu}^{\operatorname{exact}}(\mu,b_{T})=\mu^{2}e^{-g(\mu,b_{T})/\mathcal{D}(\mu,b_{T})},$
(83)

where

$v(\mu,b_{T})=\frac{\gamma_{1}}{\Gamma_{0}}+a_{s}\left[\frac{\beta_{0}}{12}\mathbf{L}_{\mu}^{2}+\frac{\gamma_{2}+d_{2}(0)}{\Gamma_{0}}-\frac{\gamma_{1}\Gamma_{1}}{\Gamma_{0}^{2}}\right],$
(84)

and

$\displaystyle
g(\mu,b_{T})=\frac{1}{a_{s}}\frac{\Gamma_{0}}{2\beta_{0}^{2}}\left\\{e^{-p}-1+p+a_{s}\left[\frac{\beta_{1}}{\beta_{0}}\left(e^{-p}-1+p-\frac{p^{2}}{2}\right)\right.\left.-\frac{\Gamma_{1}}{\Gamma_{0}}\left(e^{-p}-1+p\right)+\frac{\beta_{0}\gamma_{1}}{\Gamma_{0}}p\right]\right.+$
(85) $\displaystyle
a_{s}^{2}\left[\left(\frac{\Gamma_{1}^{2}}{\Gamma_{0}^{2}}-\frac{\Gamma_{2}}{\Gamma_{0}}\right)(\cosh
p-1)\left.+\left(\frac{\beta_{1}\Gamma_{1}}{\beta_{0}\Gamma_{0}}-\frac{\beta_{2}}{\beta_{0}}\right)(\sinh
p-p)+\left(\frac{\beta_{0}\gamma_{2}}{\Gamma_{0}}-\frac{\beta_{0}\gamma_{1}\Gamma_{1}}{\Gamma_{0}^{2}}\right)\left(e^{p}-1\right)\right]\right\\}.$

In the $g(\mu,b_{T})$, the $p$ is

$p=\frac{2\beta_{0}\mathcal{D}(\mu,b_{T})}{\Gamma_{0}},$ (86)

and in the $v(\mu,b_{T})$, the $d_{2}(0)$ is

$d_{2}(0)=C_{F}C_{A}\left(\frac{404}{27}-14\zeta_{3}\right)-\frac{112}{27}T_{R}N_{f}C_{F}.$
(87)

## References

  * [1] R. D. Tangerman and P. J. Mulders, Intrinsic transverse momentum and the polarized Drell-Yan process, Phys. Rev. D 51, 3357 (1995) [arXiv:hep-ph/9403227].
  * [2] P. J. Mulders and R. D. Tangerman, The Complete tree level result up to order 1/Q for polarized deep inelastic leptoproduction, Nucl. Phys. B 461, 197 (1996); Erratum, Nucl. Phys. B 484, 538 (1997) [arXiv:hep-ph/9510301].
  * [3] D. Boer and P. J. Mulders, Time reversal odd distribution functions in leptoproduction, Phys. Rev. D 57, 5780-5786 (1998) [arXiv:hep-ph/9711485].
  * [4] X. d. Ji and F. Yuan, Parton distributions in light cone gauge: Where are the final state interactions?, Phys. Lett. B 543, 66-72 (2002) [arXiv:hep-ph/0206057].
  * [5] J. Zhu and B. Q. Ma, Three-dimensional parton distribution functions $g_{1T}$ and $h_{1L}^{\perp}$ in the polarized proton-antiproton Drell-Yan process, Eur. Phys. J. C 71, 1807 (2011) [arXiv:1104.5545].
  * [6] J. Zhu and B. Q. Ma, Proposal for measuring new transverse momentum dependent parton distributions $g_{1T}$ and $h_{1L}^{\perp}$ through semi-inclusive deep inelastic scattering, Phys. Lett. B 696, 246-251 (2011) [arXiv:1104.4564].
  * [7] Xuan Luo and H. Sun, Kotzinian-Mulders effect in semi-inclusive DIS within TMD factorization *, Chin. Phys. C 46, 023102 (2022) [arXiv:2004.13508].
  * [8] A. M. Kotzinian and P. J. Mulders, Longitudinal quark polarization in transversely polarized nucleons, Phys. Rev. D 54, 1229-1232 (1996) [arXiv:hep-ph/9511420].
  * [9] G. A. Miller, Densities, Parton Distributions, and Measuring the Non-Spherical Shape of the Nucleon, Phys. Rev. C 76, 065209 (2007) [arXiv:0708.2297].
  * [10] A. Bacchetta, F. Conti and M. Radici, Transverse-momentum distributions in a diquark spectator model, Phys. Rev. D 78, 074010 (2008) [arXiv:0807.0323].
  * [11] B. Pasquini, S. Cazzaniga and S. Boffi, Transverse momentum dependent parton distributions in a light-cone quark model, Phys. Rev. D 78, 034025 (2008) [arXiv:0806.2298].
  * [12] B. Pasquini, S. Boffi and P. Schweitzer, The Spin Structure of the Nucleon in Light-Cone Quark Models, Mod. Phys. Lett. A 24, 2903-2912 (2009) [arXiv:0910.1677].
  * [13] M. Burkardt, Spin-orbit correlations and single-spin asymmetries, [arXiv:0709.2966].
  * [14] R. Jakob, P. J. Mulders and J. Rodrigues, Modeling quark distribution and fragmentation functions, Nucl. Phys. A 626, 937-965 (1997) [arXiv:hep-ph/9704335].
  * [15] A. V. Efremov, P. Schweitzer, O. V. Teryaev and P. Zavada, Transverse momentum dependent distribution functions in a covariant parton model approach with quark orbital motion, Phys. Rev. D 80, 014021 (2009) [arXiv:0903.3490].
  * [16] H. Avakian, A. V. Efremov, P. Schweitzer and F. Yuan, The transverse momentum dependent distribution functions in the bag model, Phys. Rev. D 81, 074035 (2010) [arXiv:1001.5467].
  * [17] A. Bacchetta, M. Boglione, A. Henneman and P. J. Mulders, Bounds on transverse momentum dependent distribution and fragmentation functions, Phys. Rev. Lett. 85, 712-715 (2000) [arXiv:hep-ph/9912490].
  * [18] X. d. Ji, J. P. Ma and F. Yuan, Three quark light cone amplitudes of the proton and quark orbital motion dependent observables, Nucl. Phys. B 652, 383-404 (2003) [arXiv:hep-ph/0210430].
  * [19] S. Boffi, A. V. Efremov, B. Pasquini and P. Schweitzer, Azimuthal spin asymmetries in light-cone constituent quark models, Phys. Rev. D 79, 094012 (2009) [arXiv:0903.1271].
  * [20] A. Bacchetta, M. Radici, F. Conti and M. Guagnelli, Weighted azimuthal asymmetries in a diquark spectator model,’’ Eur. Phys. J. A 45, 373-388 (2010) [arXiv:1003.1328].
  * [21] L. P. Gamberg, G. R. Goldstein and M. Schlegel, Transverse Quark Spin Effects and the Flavor Dependence of the Boer-Mulders Function, Phys. Rev. D 77, 094016 (2008) [arXiv:0708.0324].
  * [22] P. V. Pobylitsa, Transverse momentum dependent parton distributions in large N(c) QCD, [arXiv:hep-ph/0301236].
  * [23] H. Avakian, A. V. Efremov, K. Goeke, A. Metz, P. Schweitzer and T. Teckentrup, Are there approximate relations among transverse momentum dependent distribution functions?, Phys. Rev. D 77, 014023 (2008) [arXiv:0709.3253].
  * [24] A. Accardi, A. Bacchetta, W. Melnitchouk and M. Schlegel, What can break the Wandzura-Wilczek relation?, J. High Energy Phys. 11 (2009) 093 [arXiv:0907.2942].
  * [25] K. Kanazawa, Y. Koike, A. Metz, D. Pitonyak and M. Schlegel, Operator Constraints for Twist-3 Functions and Lorentz Invariance Properties of Twist-3 Observables, Phys. Rev. D 93, 054024 (2016) [arXiv:1512.07233].
  * [26] I. Scimemi and A. Vladimirov, Matching of transverse momentum dependent distributions at twist-3, Eur. Phys. J. C 78, 802 (2018) [arXiv:1804.08148].
  * [27] P. Hagler, B. U. Musch, J. W. Negele and A. Schafer, Intrinsic quark transverse momentum in the nucleon from lattice QCD, EPL 88, 61001 (2009) [arXiv:0908.1283].
  * [28] B. U. Musch, P. Hagler, J. W. Negele and A. Schafer, Exploring quark transverse momentum distributions with lattice QCD, Phys. Rev. D 83, 094507 (2011) [arXiv:1011.1213].
  * [29] B. Yoon, M. Engelhardt, R. Gupta, T. Bhattacharya, J. R. Green, B. U. Musch, J. W. Negele, A. V. Pochinsky, A. Schäfer and S. N. Syritsyn, Nucleon Transverse Momentum-dependent Parton Distributions in Lattice QCD: Renormalization Patterns and Discretization Effects, Phys. Rev. D 96, 094508 (2017) [arXiv:1706.03406].
  * [30] 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 [arXiv:2007.07755].
  * [31] 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-145 (2017) [arXiv:1609.07374].
  * [32] B. Parsamyan, Measurement of target-polarization dependent azimuthal asymmetries in SIDIS and Drell-Yan processes at COMPASS experiment, PoS QCDEV2017, 042 (2018)
  * [33] H. Avakian, B. Parsamyan and A. Prokudin, Spin orbit correlations and the structure of the nucleon, Riv. Nuovo Cim. 42, 1-48 (2019) [arXiv:1909.13664].
  * [34] J. Huang et al. (Jefferson Lab Hall A Collaboration), Beam-Target Double Spin Asymmetry $A_{LT}$ in Charged Pion Production from Deep Inelastic Scattering on a Transversely Polarized 3He Target at 1.4 $<$ Q2 $<$ 2.7 GeV2, Phys. Rev. Lett. 108, 052001 (2012) [arXiv:1108.0489].
  * [35] S. Bhattacharya, Z. B. Kang, A. Metz, G. Penn and D. Pitonyak, First global QCD analysis of the TMD g1T from semi-inclusive DIS data, Phys. Rev. D 105, 034007 (2022) [arXiv:2110.10253].
  * [36] M. Horstmann, A. Schafer and A. Vladimirov, Study of the worm-gear-T function g1T with semi-inclusive DIS data, Phys. Rev. D 107, 034016 (2023) [arXiv:2210.07268].
  * [37] A. Accardi, J. L. Albacete, M. Anselmino, N. Armesto, E. C. Aschenauer, A. Bacchetta, D. Boer, W. K. Brooks, T. Burton and N. B. Chang, et al. Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all, Eur. Phys. J. A 52, 268 (2016) [arXiv:1212.1701].
  * [38] R. Abdul Khalek, A. Accardi, J. Adam, D. Adamiak, W. Akers, M. Albaladejo, A. Al-bataineh, M. G. Alexeev, F. Ameli and P. Antonioli, et al. Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report, Nucl. Phys. A 1026, 122447 (2022) [arXiv:2103.05419].
  * [39] 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 [arXiv:2103.10276].
  * [40] A. Bacchetta, M. Diehl, K. Goeke, A. Metz, P. J. Mulders and M. Schlegel, Semi-inclusive deep inelastic scattering at small transverse momentum, J. High Energy Phys. 02 (2007) 093 [arXiv:hep-ph/0611265].
  * [41] A. Bacchetta, U. D’Alesio, M. Diehl and C. A. Miller, Single-spin asymmetries: The Trento conventions, Phys. Rev. D 70, 117504 (2004) [hep-ph/0410050].
  * [42] 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 [arXiv:1912.06532].
  * [43] J. C. Collins and D. E. Soper, Back-To-Back Jets: Fourier Transform from B to K-Transverse, Nucl. Phys. B 197, 446-476 (1982)
  * [44] C. Zeng, T. Liu, P. Sun and Y. Zhao, Toward three-dimensional nucleon structures at the Electron-Ion Collider in China: A study of the Sivers function, Phys. Rev. D 106, 094039 (2022) [arXiv:2208.14620].
  * [45] C. Zeng, H. Dong, T. Liu, P. Sun and Y. Zhao, Role of sea quarks in the nucleon transverse spin, Phys. Rev. D 109, no.5, 056002 (2024) [arXiv:2310.15532].
  * [46] S. Tafat, Nonperturbative corrections to the Drell-Yan transverse momentum distribution, J. High Energy Phys. 05 (2001) 004 [arXiv:hep-ph/0102237].
  * [47] A. A. Vladimirov, Self-contained definition of the Collins-Soper kernel, Phys. Rev. Lett. 125, 192002 (2020) [arXiv:2003.02288].


October 17, 2016

# Spread of infectious diseases: Effects of the treatment of population

Kazue Kudo and Kanako Mizuno Department of Computer ScienceDepartment of
Computer Science Ochanomizu University Ochanomizu University 2-1-1 Ohtsuka
2-1-1 Ohtsuka Bunkyo-ku Bunkyo-ku Tokyo 112-8610 Tokyo 112-8610 Japan Japan
<EMAIL_ADDRESS>
###### Abstract

In a metapopulation network, infectious diseases spread widely because of the
travel of individuals. In the present study, we consider a modified
metapopulation Susceptible-Infected-Removed (SIR) model with a latent period,
which we call the SHIR model. In the SHIR model, an infectious period is
divided into two stages. In the first stage, which corresponds to the latent
period, infectious individuals can travel. However, in the second stage, the
same individuals cannot travel since they are seriously ill. Final size
distributions of the metapopulation SIR and SHIR models are simulated with two
different methods and compared. In Monte Carlo simulations, in which the
population is treated as an integer, the distributions show similar behavior.
However, in reaction-diffusion systems, in which the population is treated as
a real number, the final size distribution of the SHIR model has a
discontinuous jump, and that of the SIR model shows a continuous transition.
The discontinuous jump is found to be an artifact that occurs owing to an
inappropriate termination condition.

epidemic model, metapopulation network

##  1 Introduction

The spread of infectious diseases is often modeled on a complex network. If
the network represents a human network, each node corresponds to an
individual. Transmission occurs between connected nodes. The epidemic
condition depends on the structure of the network as well as the transmission
and recovery rates. The impact of an epidemic is characterized by the basic
reproduction number $R_{0}$, which is the expected number of infections caused
by a typical infectious individual in a completely susceptible population [1,
2]. $R_{0}=1$ represents the epidemic threshold. In the Susceptible-Infected-
Removed (SIR) model, $R_{0}$ is given simply by the ratio of the transmission
rate to the recovery rate when the network is a complete graph, or under a
mean-field approximation. However, when the network is a complex network,
e.g., a scale-free network, the epidemic threshold depends also on the average
connectivity between the nodes [3].

Human mobility is also an important factor behind the spread of infectious
diseases. Metapopulation models are often used to implement the movement of
individuals. In a metapopulation network, each node corresponds to a
subpopulation or patch, which contains individuals. In metapopulation SIR
models, while infection dynamics occur only within each patch, individuals can
move between neighboring patches [4, 5, 6, 7]. The epidemic threshold in each
patch is given by the basic reproduction number $R_{0}$. The epidemic in a
metapopulation system is characterized by not only the epidemic threshold but
also the global invasion threshold, which depends on the mobility rate and
network structure, as well as the basic reproduction number [5, 6, 7].

In the current study, we consider a modified SIR model with a latent period,
which we call the SHIR model, on a metapopulation network [8]. In the SHIR
model, an infectious period is divided into two stages: an infected stage
(denoted by “H”), which corresponds to the latent period, and a seriously-ill
stage (denoted by “I”), in which the individual is infected and cannot move to
another patch. In Ref. [8], unexpected results were reported; a discontinuous
jump appeared in the final size distribution, which expresses the number of
recovered individuals in the SHIR model, and the global invasion threshold was
almost independent of the mobility rate. In this paper, we explain the reason
behind those results, based on numerical simulations performed with two
different methods.

The rest of the paper is organized as follows. The metapopulation SIR and SHIR
models are introduced in Sec. 2. In Sec. 3, we demonstrate the final size
distributions simulated with two different methods. The time evolution of the
number of infected individuals is also shown to explain the discontinuous jump
in the final size distribution. Finally, conclusions are presented in Sec. 4.

##  2  Models

First, we introduce a metapopulation SIR model in which individuals can move
between connected patches. In the metapopulation SIR model, each patch can
contain many individuals, and each individual is in one of the three states:
Susceptible, Infected, or Recovered. Infection occurs in each patch, and a
disease is spread by the movement of individuals. Each individual can move
from one patch to another that is connected with it. If we assume that the
mobility rate $w$ is the same in the entire network, the time evolution of
$S_{i}$, $I_{i}$, and $R_{i}$ (the numbers of susceptible, infected, and
recovered individuals in the $i$-th patch, respectively) is described by the
following equations.

$\displaystyle\frac{dS_{i}}{dt}$
$\displaystyle=-\beta\frac{S_{i}I_{i}}{N_{i}}+w\sum_{j\in\partial
i}(S_{j}-S_{i}),$ (1a) $\displaystyle\frac{dI_{i}}{dt}$
$\displaystyle=\beta\frac{S_{i}I_{i}}{N_{i}}-\gamma I_{i}+w\sum_{j\in\partial
i}(I_{j}-I_{i}),$ (1b) $\displaystyle\frac{dR_{i}}{dt}$ $\displaystyle=\gamma
I_{i}+w\sum_{j\in\partial i}(R_{j}-R_{i}).$ (1c)

Here, $\beta$ and $\gamma$ are the transmission and recovery rates,
respectively, and $N_{i}=S_{i}+I_{i}+R_{i}$ is the total number of individuals
in the $i$-th patch. The notation $\partial i$ represents the set of neighbors
of the $i$-th patch. In other words, the summations are taken over all the
patches that are connected to the $i$-th patch.

Next, we modify the metapopulation SIR model to take into account a latent
period. If a person were sick, he/she probably would not go out. Then,
infected people who can travel are considered to be the ones who are in a
latent period. Under this assumption, an infectious period is divided into two
stages: infected (in a latent period) and seriously-ill stages. In the first
stage, an individual can infect other individuals and travel between
neighboring patches. However, in the second stage, an individual cannot
travel. In the modified model, the time evolution of $S_{i}$, $H_{i}$ ,
$I_{i}$ and $R_{i}$, (i.e., the number of susceptible, infected, seriously-
ill, and recovered individuals in the $i$-th patch) are described by the
following equations.

$\displaystyle\frac{dS_{i}}{dt}$
$\displaystyle=-\beta\frac{S_{i}(H_{i}+I_{i})}{N_{i}}+w\sum_{j\in\partial
i}(S_{j}-S_{i}),$ (2a) $\displaystyle\frac{dH_{i}}{dt}$
$\displaystyle=\beta\frac{S_{i}(H_{i}+I_{i})}{N_{i}}-\mu
H_{i}+w\sum_{j\in\partial i}(H_{j}-H_{i}),$ (2b)
$\displaystyle\frac{dI_{i}}{dt}$ $\displaystyle=\mu H_{i}-\gamma I_{i},$ (2c)
$\displaystyle\frac{dR_{i}}{dt}$ $\displaystyle=\gamma
I_{i}+w\sum_{j\in\partial i}(R_{j}-R_{i}).$ (2d)

Here, $\mu$ is the transition rate from the infected to the seriously-ill
stage, and $N_{i}=S_{i}+H_{i}+I_{i}+R_{i}$. We call this model the SHIR model.
This model is different from the SEIR model [9], which is an epidemic model in
which a latent period is expressed as an “Exposed” state. The exposed
individual is infected, but does not infect others. However, the SHIR model as
well as the SEIR model is one of the generalized SIR models that include
multiple infectious stages [2].

The basic reproduction number $R_{0}$ in the standard SIR model (without
mobility) is given by $\beta/\gamma$, and that in the generalized SIR model
that includes $n$ infectious stages is given by

$R_{0}=\sum_{i=1}^{n}\frac{\beta_{i}}{\gamma_{i}},$ (3)

where $\beta_{i}$ is the transmission rate in the $i$-th infectious stage, and
$1/\gamma_{i}$ is the mean duration an individual spends in the $i$-th
infectious stage [2, 10]. In Eq. (2), $\beta_{1}=\beta_{2}=\beta$,
$\gamma_{1}=\mu$ and $\gamma_{2}=\gamma$. Thus, we find
$R_{0}=\beta(\mu+\gamma)/(\mu\gamma)$ for the SHIR model (without mobility)
[8].

##  3  Final Size Distributions

In this section, we demonstrate the final size distributions of the
metapopulation SIR and SHIR models. The final size of an epidemic is the
cumulative number of infected individuals, which is equal to the total number
of recovered individuals at the end of the epidemic. The attack ratio, which
is the normalized final size, is given by the ratio of the number of recovered
individuals to the total number of individuals.

The metapopulation network we use is a scale-free network of $500$ nodes,
whose degree distribution is $P(k)\sim k^{-\Gamma}$ with $\Gamma=3.5$.
Initially, one individual is infected, and the others are susceptible. The
initial number of individuals at each node is taken as $\bar{N}k/\langle
k\rangle$, where $\langle k\rangle$ is the average degree, with $\bar{N}=100$,
unless otherwise specified. In the numerical simulations below, we fix
parameters as $\gamma=0.25$ for the SIR and $\gamma=\mu=0.5$ for the SHIR
models. Then, the epidemic threshold $R_{0}=1$ corresponds to $\beta=0.25$ in
both the models.

Figure 1: Final size distributions simulated with a Monte Carlo method for (a)
a metapopulation SIR model and (b) a metapopulation SHIR model. $w$ and
$\beta$ are the mobility and transmission rates, respectively.

We perform numerical simulations with two different methods. The first one is
the Monte Carlo method. One Monte Carlo step consists of transmission and
traveling processes. In the case of the SIR model, in the $i$-th patch, a
susceptible individual becomes infected with probability
$1-(1-\frac{\beta}{N_{i}}\tau)^{I_{i}}$, where $\tau$ is the time scale, and
each infected individual recovers with probability $\gamma\tau$ [6]. In the
SHIR model, a susceptible individual becomes infected with probability
$1-(1-\frac{\beta}{N_{i}}\tau)^{H_{i}+I_{i}}$, an infected individual turns
seriously ill with probability $\mu\tau$, and a seriously-ill individual
recovers with probability $\gamma\tau$. Each individual who can travel, moves
from one patch with degree $k$ to another neighboring patch with probability
$w\tau/k$. In Monte Carlo simulations, $S_{i}$, $H_{i}$, $I_{i}$, and $R_{i}$
are treated as integers. In one Monte Carlo step, in which the time scale is
taken as $\tau=1$, the transmission process is followed by the traveling
process. States of the individuals are updated after each process. Simulation
continues as long as $\sum_{i}I_{i}\geq 1$ for the SIR model and
$\sum_{i}(H_{i}+I_{i})\geq 1$ for the SHIR model. The final size distributions
for the metapopulation SIR and SHIR models are shown in Figs. 1(a) and 1(b),
respectively. In both the distributions, the attack ratio starts to increase
near $\beta=0.25$, which corresponds to the epidemic threshold of $R_{0}=1$.
When the mobility rate is small, the attack ratio is also small in both Figs.
1(a) and 1(b).

Figure 2: Final size distributions in reaction-diffusion systems for (a) a
metapopulation SIR model and (b) a metapopulation SHIR model. $w$ and $\beta$
are mobility and transmission rates, respectively.

The second numerical method involves solving the reaction-diffusion equations
(1) and (2) numerically. We employ the 4th Runge-Kutta method for time
evolution. In this method, $S_{i}$, $H_{i}$, $I_{i}$, and $R_{i}$ are treated
as real numbers. In Ref. [8], termination conditions, which are conditions to
stop a simulation, were $\sum_{i}I_{i}<1$ for the SIR model and
$\sum_{i}(H_{i}+I_{i})<1$ for the SHIR model. Figures 2(a) and 2(b) show the
final size distributions under these termination conditions for the
metapopulation SIR and SHIR models, respectively. In Fig. 2(a), we see a
continuous transition at $\beta=0.25$ ($R_{0}=1$). However, in Fig. 2(b), a
discontinuous jump appears at $\beta=0.5$. Another important point is that the
final size distributions are not related to the mobility rate $w$. These
results are almost the same as the ones reported in Ref. [8].

The reason why the final size distribution is independent of $w$ is owing to
the fact that the second method is a type of a mean-field approach in which
the network structure is irrelevant. In other words, the entire metapopulation
network can be regarded as a large single population in the mean-field
approach. In such a case, attack ratio does not depend on the mobility rate,
although the time scales for the epidemic to spread and cease depend on the
mobility rate. By contrast, in Monte Carlo simulations, the final size
distributions are strongly dependent on the mobility rate. Furthermore, in
this case, the network structure is relevant as well. When the mobility rate
is very small, the epidemic cannot spread to other patches even if the
transmission rate is large, although a local outbreak may occur. Since a large
fraction of patches stay uninfected for a small mobility rate, the attack
ratio is small even for a large transmission rate.

Figure 3: Final size distributions as functions of $w\bar{N}$ for (a) a
metapopulation SIR model and (b) a metapopulation SHIR model. The transmission
rate is $\beta=1$. Monte Carlo simulations for $\bar{N}=50$, $100$, $500$,
$1000$, $5000$, and $10000$ are compared with the result in a reaction-
diffusion system, namely, a mean-field (MF) approach.

Moreover, in Monte Carlo simulations, the attack ratio depends not only on the
mobility rate $w$ but also on the average population $\bar{N}$ of a patch. The
quantity that defines the global invasion threshold is proportional to
$w\bar{N}$, and depends on the network structure and other parameters [6]. In
Fig. 3, the attack ratio is plotted as a function of $w\bar{N}$ for $\beta=1$.
For different values of $\bar{N}$, the attack ratio simulated with the Monte
Carlo method collapses into a single curve. The attack ratio increases rapidly
at around $w\bar{N}\sim 1$ and saturates to a certain value. The saturated
value of attack ratio in the Monte Carlo simulations is approximately
$0.7$–$0.8$, while that in the mean-field approach is almost $1$. This is
because infected individuals disappear at some probability before the epidemic
grows. This type of behavior is often observed when the initial condition has
only one infected individual in the entire metapopulation.

Figure 4: (a) Final size distributions for the metapopulation SHIR model with
$w=1$. Termination conditions are $\sum_{i}(H_{i}+I_{i})<1$ for Condition 1
and $\sum_{i}(H_{i}+I_{i})=0$ for Condition 2. Interval $0<\beta<1$ is
separated into three regimes: (i) $0<\beta<0.25$, (ii) $0.25<\beta<0.5$, and
(iii) $0.5<\beta<1$. (b) Time evolution of the number of infected and
seriously-ill individuals $\sum_{i}(H_{i}+I_{i})$. $\beta=0.2,0.4$ and $0.6$
are representative points in Regimes (i), (ii), and (iii), respectively.

Here, we explain the reason for the discontinuous jump in the reaction-
diffusion system. The key is the termination condition. When the termination
threshold is almost zero, the discontinuous jump disappears. The final size
distributions under different termination conditions are shown in Fig. 4(a),
where the mobility rate is fixed at $w=1$. The green curve in Fig. 4(a) is the
same as the final size distribution in Fig. 2(b). The blue curve in Fig. 4(a)
is calculated under the termination condition $\sum_{i}(H_{i}+I_{i})=0$ (more
precisely, $\sum_{i}(H_{i}+I_{i})<\epsilon$ with $\epsilon=10^{-14}$). It has
a continuous transition at $\beta=0.25$, which corresponds to $R_{0}=1$. To
understand the effects of different termination conditions, we separate the
interval of $\beta$ into three regimes: (i) $0<\beta<0.25$, (ii)
$0.25<\beta<0.5$, and (iii) $0.5<\beta<1$. In Regime (i), the number of
infected individuals decreases monotonically as shown in Fig. 4(b). Since
$R_{0}$ is lower than the epidemic threshold, a disease cannot spread at all.
By contrast, in Regime (iii), initially, the number of infected individuals
increases rapidly. In these two regimes, termination conditions have no
influence on the attack ratio. However, in Regime (ii), the situation is a bit
complicated. The number of infected individuals decreases at first but
increases after a short interval of time. This initial decay in
$\sum_{i}(H_{i}+I_{i})$ is natural, since

$\frac{d}{dt}\sum_{i}(H_{i}+I_{i})=\sum_{i}\left[\beta\frac{S_{i}(H_{i}+I_{i})}{N_{i}}-\gamma
I_{i}\right]\simeq(\beta-\gamma)\sum_{i}(H_{i}+I_{i})$ (4)

with $H_{i}=0$ and $S_{i}\simeq N_{i}$ at $t=0$. Since $\beta<\gamma$ in
Regime (ii), the simulation stops immediately when the termination threshold
is $1$. When the threshold is closer to 0 (or $\epsilon$), the simulation
continues, and the number of infected individuals grows enough since
$R_{0}>1$. In the SIR model, $\beta=\gamma$ is equivalent to $R_{0}=1$. Thus,
Regime (ii) does not exist in the SIR model.

##  4  Conclusions

We introduced a modified metapopulation SIR model with a latent period, which
we call the SHIR model, and compared the final size distributions of the
metapopulation SIR and SHIR models. In Monte Carlo simulations, the
distributions of both the models appear similar and depend on the mobility
rate as well as the transmission rate. However, those in the reaction-
diffusion systems are independent of the mobility rate and have different
properties between the SIR and SHIR model. The difference between the two
methods arises from the treatment of population; population is treated as an
integer in Monte Carlo simulations and a real number in reaction-diffusion
systems. Although the final size distribution of the SIR model has a
continuous transition, that of the SHIR model shows a discontinuous jump,
depending on the termination condition of the simulation. The difference
caused by the termination condition indicates that a temporal decay of the
number of infected individuals can mislead us into underestimating the final
size of an epidemic.

## References

  * [1] R.M. Anderson, and R.M. May, Infectious diseases of humans: dynamics and control (Oxford University Press, Oxford, 1991).
  * [2] J. Ma and D.J.D. Earn, Bull. Math. Biol. 68 (2006) 679.
  * [3] R.M. May and A.L. Lloyd, Phys. Rev. E 64 (2001) 066112.
  * [4] M.J. Keeling and P. Rohani, Ecol. Lett. 5 (2002) 20.
  * [5] P.C. Cross, J.O. Lloyd-Smith, P.L.F. Johnson, and W.M. Getz, Ecol. Lett. 8 (2005) 587.
  * [6] V. Colizza and A. Vespignani, Phys. Rev. Lett. 99 (2007) 148701.
  * [7] V. Colizza and A. Vespignani, J. Theor. Biol. 251 (2008) 450.
  * [8] K. Mizuno and K. Kudo, in Proceedings of the International Conference on Social Modeling and Simulation, plus Econophysics Colloquium 2014, ed. H. Takayasu, N. Ito, I. Noda and M. Takayasu (Springer Proceedings in Complexity, Springer, 2015) p. 141.
  * [9] I. Schwartz and H. Smith, J. Math. Biol. 18 (1983) 233.
  * [10] J.M. Hyman, J. Li, and E.A. Stanley, Math. Biosci. 155 (1999) 77.