# The Sivers asymmetry in charged Kaon and $\Lambda$ hyperon produced SIDIS
process at electron ion colliders

Shuailiang Yang  School of Physics and Microelectronics, Zhengzhou University,
Zhengzhou, Henan 450001, China  Jianxi Song  School of Physics and
Microelectronics, Zhengzhou University, Zhengzhou, Henan 450001, China  Xiaoyu
Wang<EMAIL_ADDRESS>School of Physics and Microelectronics, Zhengzhou
University, Zhengzhou, Henan 450001, China  De-Min Li<EMAIL_ADDRESS>School
of Physics and Microelectronics, Zhengzhou University, Zhengzhou, Henan
450001, China  Zhun Lu<EMAIL_ADDRESS>School of Physics, Southeast
University, Nanjing, Jiangsu 211189, China

###### Abstract

We investigate the single transverse-spin asymmetry with a
$\sin(\phi_{h}-\phi_{S})$ modulation in the charged Kaon produced and in
$\Lambda$ hyperon produced SIDIS process within the theoretical framework of
transverse momentum dependent (TMD) factorization at the next-to-leading-
logarithmic order. The asymmetry is contributed by the convolution of Sivers
function and the unpolarized fragmentation function $D_{1}$ for the produced
hadron. The parametrization for the proton Qiu-Sterman function, which is
closely related to the Sivers function, is adopted to numerically estimate the
Sivers asymmetry at the kinematical region of Electron Ion Collider (EIC) and
Electron Ion Collider in China (EicC). The TMD evolution of the TMD parton
distribution functions are considered by employing the nonperturbative Sudakov
form factor. It is found that the predicted Sivers asymmetries
$A_{UT}^{\sin(\phi_{h}-\phi_{S})}$ as functions of $x$, $z$ and $P_{hT}$ are
sizable at the kinematical configurations of both EIC and EicC. The strange
constituent of the produced charged Kaon and $\Lambda$ hyperon in the final
state can be a promising probe of the sea quark Sivers function as well as the
flavor dependence in the proton target. Therefore, it is important to utilize
the future EIC facilities to constrain the sea quark distribution functions as
well as the validity of the generalized universality of the Sivers function.

##  I INTRODUCTION

Since the measurement by the European Muon Collaboration Ashman et al. (1988,
1989) showed that the spin fraction carried by the internal quarks is much
smaller than the spin of the proton, which is contradict to the conventional
theoretical prediction that the constituent quark spin contributes the total
proton spin, numerous studies have been carried out to explore the nucleon
spin structure from both theoretical and experimental aspects. Among the spin-
related observables, the transverse single spin asymmetries (TSSAs) can be the
key access to the information of transverse momentum structure of nucleon,
which is encoded in the transverse momentum dependent parton distribution
functions (TMD PDFs). In leading twist there are eight TMD PDFs, each one
describes a distribution of three-dimensional motion of partons with specified
polarization inside the nucleon. Particularly, the time-reversal odd (T-odd)
Sivers function $f_{1T}^{\perp}\left(x,p_{T}\right)$ Sivers (1990, 1991)
denotes the asymmetric distribution of unpolarized quarks inside a
transversely polarized nucleon, which arises from the correlation between the
quark transverse momentum and the nucleon transverse spin. Due to its T-odd
property, Sivers function as well as its chiral-odd partner the Boer-Mulders
function has been assumed to be forbidden by the naive time-reversal
invariance of QCD Collins (1993), the very existence of the two T-odd
distribution functions was not so obvious. However, the situation has changed
since the calculations in Refs. Brodsky et al. (2002a, b); Boer et al. (2003),
which showed that the T-odd distributions can actually survive using spectator
model calculations incorporating gluon exchange between the struck quark and
the spectator. In Ref. Collins (2002), the time-reversal-invariant argument
was reexamined and showed the gauge-link in the operator definition of the
correlators guarantee the T-odd distribution functions to be nonzero. More
importantly, the presence of the gauge-link indicates that Sivers function or
the Boer-Mulders function has opposite sign between semi-inclusive deeply
inelastic scattering (SIDIS) and Drell-Yan processes Brodsky et al. (2002b,
a); Collins (2002)

$\displaystyle
f_{1T}^{\perp}\left(x,p_{T}\right)_{[\mathrm{SIDIS}]}=-f_{1T}^{\perp}\left(x,p_{T}\right)_{[\mathrm{DY}]},$
(1)

which is a significant prediction by QCD. The verification of the sign change
is one of the most fundamental tests of QCD prediction, and it is also the
main pursue of the existing and future Drell-Yan facilities.

The transverse single spin asymmetry can be utilized to extract the
information of the Sivers function, and has been intensively investigated in
the past two decades from both experimental and theoretical aspects. The first
non-zero Sivers asymmetry was measured by the HERMES Collaboration at DESY in
electroproduction of charge pions off the transversely polarized hydrogen
target Airapetian et al. (2005). Updated measurements on the Sivers asymmetry
in pion produced as well as those in Kaon and $p/\bar{p}$ produced in three-
dimensional kinematic bin and enlarged phase space were reported in Refs.
Airapetian et al. (2009, 2020). The COMPASS Collaboration at CERN also
measured the Sivers asymmetries in charged hadrons produced SIDIS process
through muon beam scattering off the transversely polarized proton and
deuteron targets Alexakhin et al. (2005); Ageev et al. (2007); Alekseev et al.
(2009, 2010); Adolph et al. (2012, 2017). In addition, the data on the
weighted Sivers asymmetry are also released in Ref. Alexeev et al. (2019),
allowing for the extractions of the Sivers function and its first transverse
moment. The Hall A Collaboration at Jefferson Lab presented the measurement of
TSSA in charged pion produced SIDIS process with a transversely polarized
${}^{3}\mathrm{He}$ target Qian et al. (2011); Zhao et al. (2014). Besides the
SIDIS process, COMPASS also measured the Sivers asymmetry in Drell-Yan process
via $\pi\,N$ collision Aghasyan et al. (2017). Measurement of TSSAs sensitive
to Sivers function in the $W^{\pm}$ boson produced in proton-proton collisions
has also been performed by the STAR experiments at RHIC Adamczyk et al.
(2016). The data from these measurements have been applied to extract the
Sivers function using parametrizations and phenomenological approaches
Anselmino et al. (2005); Efremov et al. (2005); Collins et al. (2006);
Vogelsang and Yuan (2005); Anselmino et al. (2009, 2012); Bacchetta and Radici
(2011); Echevarria et al. (2014a); Anselmino et al. (2017); Martin et al.
(2017); Boglione et al. (2018); Bury et al. (2021).

From the theoretical aspect, intensive studies on the quark Sivers function
were performed using the QCD-inspired models, such as the spectator model
Brodsky et al. (2002a); Boer et al. (2003); Bacchetta et al. (2004a), the
light-cone quark model Lu and Ma (2004); Pasquini and Yuan (2010), the light-
front quark-diquark model Maji et al. (2018, 2017), the non-relativistic
constituent quark model Courtoy et al. (2008), the MIT bag model Yuan (2003);
Courtoy et al. (2009), and the Holographic QCD Lyubovitskij et al. (2022). The
sea quark Sivers function has been estimated from the light-cone wave function
in Refs. Dong et al. (2019); He and Wang (2019); Luan and Lu (2022). However,
so far only the valence quark Sivers functions are constrained in the valence
region with relatively large uncertainties in the transverse momentum space.
One of the reasons is that it is difficult to describe the corresponding
physical observables since there are complicated TMD effects. The TMD
evolution effects are encoded in the Sudakov-like form factor which also
includes details of the non-perturbative QCD dynamics. Therefore, this part of
the Sudakov-like form factor can not be calculated from the perturbative QCD
and is mostly unknown. Another reason is that the knowledge of the scale
dependence of the TSSAs is very limited since the measurements are mostly
performed in the fixed-target experiments with similar hard scales. With the
expected high energy and high precision of the Electron Ion Collider (EIC)
Accardi et al. (2016); Abdul Khalek et al. (2022) and the Electron Ion
Collider in China (EicC) Zeng et al. (2022), the precise knowledge on the TMD
distribution functions may be gained, not only for the valence quarks, but
also for sea quarks and gluons. Concerning the sea quark TMDs, the charged
Kaon produced or the Lambda hyperon produced in SIDIS can be recognized as an
ideal probe to the sea quark distribution of the proton due to the strange
constituent quark inside the kaon and the Lambda hyperon. Therefore, through
the Sivers asymmetry in $K^{\pm}$ produced and in Lambda produced off
transversely polarized nucleon at EIC and EicC, there might be an opportunity
to obtain the information of the Sivers distribution function of the sea quark
as well as its flavor dependence.

The purpose of this work is to evaluate the Sivers asymmetry in
$ep^{\uparrow}\rightarrow eK^{\pm}X$ and in $ep^{\uparrow}\rightarrow e\Lambda
X$ at the kinematical region of EIC and EicC. The theoretical tool adopted in
this study is the TMD factorization formalism Collins and Soper (1981);
Collins et al. (1985); Collins (2013); Ji et al. (2004), which has been widely
applied to various high energy processes, such as SIDIS Collins and Soper
(1981); Collins (2013); Ji et al. (2005); Aybat and Rogers (2011); Collins and
Rogers (2013); Echevarria et al. (2013), $e^{+}e^{-}$ annihilation Collins
(2013); Pitonyak et al. (2014); Boer (2009), Drell-Yan Collins (2013); Arnold
et al. (2009), and $\mathrm{W}/\mathrm{Z}$ produced in hadron collision
Collins et al. (1985); Collins (2013); Lambertsen and Vogelsang (2016). In
this framework, the differential cross section can be written as the
convolution of the well-defined TMD PDFs and/or fragmentation functions (FFs)
at the small transverse momentum region $P_{hT}\ll Q$ as an approximation. The
energy dependence of the TMD PDFs and FFs is encoded in the TMD evolution
equations, their solutions are usually given in $b$ space, which is conjugate
to the transverse momentum space Collins et al. (1985); Collins (2013) through
Fourier transformation. After solving the TMD evolution equations, the scale
dependence of the TMDs may be included in the exponential form of the so-
called Sudakov-like form factor Collins et al. (1985); Collins (2013); Aybat
and Rogers (2011); Collins and Hautmann (2000). The Sudakov-like form factor
can be further separated into perturbatively calculable part and the
nonperturbative part, the latter one can not be calculated through
perturbative theory and may be obtained by fitting experimental data. We will
consider the TMD evolution effects of the corresponding TMDs to obtain the
numerical results for the Sivers asymmetry in charged Kaon $K^{\pm}$ and
$\Lambda$ hyperon produced SIDIS process.

The rest of the paper is organized as follows. In Sec. II, we provide the
theoretical framework of Sivers asymmetry $A_{UT}^{\sin(\phi_{h}-\phi_{S})}$
in the charged Kaon produced and $\Lambda$ hyperon produced in SIDIS process
within the TMD factorization formalism. In Sec. III, we perform the numerical
estimation of the Sivers asymmetry at the kinematical region of EIC and EicC.
In Sec. IV, we summarize the work and discuss the results.

##  II FORMALISM OF THE SIVERS ASYMMETRY IN SIDIS PROCESS

In this section, we will set up the detailed formalism of the Sivers asymmetry
with a modulation of $\sin(\phi_{h}-\phi_{S})$ in the charged Kaon or Lambda
produced SIDIS process with an unpolarized electron beam scattered off a
transversely polarized proton target

$\displaystyle e(\ell)+p(P)^{\uparrow}\rightarrow
e\left(\ell^{\prime}\right)+K^{\pm}/\Lambda\left(P_{h}\right)+X,$ (2)

where $\ell$ and $\ell^{\prime}$ represent the four-momenta of the incoming
and outgoing electrons, $P$ is the four-momentum of the target proton, the up-
arrow represents the proton is transversely polarized, $P_{h}$ stands for the
four-momentum of the final-state hadron, which can be charged Kaon $K^{\pm}$
or Lambda hyperon. The four-momentum of the exchanged virtual photon is
$q=\ell-\ell^{\prime}$ and the usual defined energy scale is $Q^{2}=-q^{2}$.
We denote the masses of the proton target and the final-state hadron by $M$
and $M_{h}$. To express the differential cross section as well as the physical
observables, we adopt the following Lorentz invariants

$\displaystyle S=(P+\ell)^{2}\,,\quad x=\frac{Q^{2}}{2P\cdot q}\,,\quad
y=\frac{P\cdot q}{P\cdot\ell}\,,\quad z=\frac{P\cdot P_{h}}{P\cdot q}\,,$ (3)

where $S$ is the squared center of mass energy, $x$ represents the Bjorken
variable, $y$ represents the lepton (quark) energy momentum transferring
fraction, and $z$ represents the longitudinal momentum fraction of the final
fragmented hadron to the parent quark.

The reference frame applied in our study is shown in Fig. 1. According to the
Trento convention Bacchetta et al. (2004b), the $z$-axis is defined by the
direction of the virtual photon momentum. The azimuthal angle $\phi_{h}$ of
the outgoing hadron (charged Kaon or Lambda) is defined by

$\cos\phi_{h}=-\frac{\ell_{\mu}P_{h\nu}g_{\perp}^{\mu\nu}}{\sqrt{\ell_{T}^{2}P_{hT}^{2}}},$
(4)

with $\ell^{\mu}_{T}=g_{\perp}^{\mu\nu}\ell_{\nu}$ and
$P_{hT}^{\mu}=g_{\perp}^{\mu\nu}P_{h\nu}$ being the transverse components of
$\ell$ and $P_{h}$ respect to $z$-axis. The tensor $g_{\perp}^{\mu\nu}$ is

$g_{\perp}^{\mu\nu}=g^{\mu\nu}-\frac{q^{\mu}P^{\nu}+P^{\mu}q^{\nu}}{P(1+\gamma^{2})}+\frac{\gamma^{2}}{1+\gamma^{2}}(\frac{q^{\mu}q^{\nu}}{Q^{2}}-\frac{P^{\mu}P^{\nu}}{M^{2}}),$
(5)

with $\gamma=\frac{2Mx}{Q}$. The azimuthal angle $\phi_{S}$ of the proton spin
vector $S$ is defined by replacing $P_{h}$ by $S$ in Eq. (4), and the
transverse component of $S$ is $S_{T}^{\mu}=g_{\perp}^{\mu\nu}S_{\nu}$ similar
to the definition of $P_{hT}$.

Figure 1:  The reference frame in SIDIS process.

Assuming one photon exchange, the model-independent differential cross section
can be written as a set of structure functions with the general form as
Bacchetta et al. (2007)

$\displaystyle\frac{d\sigma}{dxdydzd\phi_{S}d\phi_{h}dP_{hT}^{2}}=\frac{\alpha_{em}^{2}}{xyQ^{2}}\frac{y}{2(1-\epsilon)}(1+{\gamma^{2}\over
2x})\left\\{F_{UU,T}+\left|S_{\perp}\right|\sin\left(\phi_{h}-\phi_{S}\right)F_{UT,T}^{\sin\left(\phi_{h}-\phi_{S}\right)}+\ldots\right\\},$
(6)

where $F_{UU,T}$ stands for the unpolarized structure function,
$F_{UT,T}^{\sin(\phi_{h}-\phi_{s})}$ is the transverse spin-dependent
structure function contributed by the Sivers function, $\epsilon$ is the ratio
of the longitudinal flux and the transverse flux of the photon which has the
definition $\epsilon=\frac{1-y-{1\over 4}\gamma^{2}y^{2}}{1-y+{1\over
2}y^{2}+{1\over 4}\gamma^{2}y^{2}}$ and the ellipsis denotes other structure
functions, which will not be considered in this work. The three subscripts in
the structure functions $F_{XY,Z}$ stand for the polarization of the lepton
beam ($X$), the target proton ($Y$) and the virtual photon ($Z$) with U being
unpolarized, T being transversely polarized. The Sivers asymmetry is defined
as the ratio of the difference between the spin-dependent differential cross
sections and the unpolarized differential cross section

$\displaystyle
A_{UT}^{\sin\left(\phi_{h}-\phi_{s}\right)}\equiv\frac{d\sigma^{\uparrow}-d\sigma^{\downarrow}}{d\sigma^{\uparrow}+d\sigma^{\downarrow}}=\frac{\sigma_{0}\left(x,y,Q^{2}\right)}{\sigma_{0}\left(x,y,Q^{2}\right)}\frac{F_{UT}^{\sin\left(\phi_{h}-\phi_{i}\right)}}{F_{UU}},$
(7)

where $\sigma_{0}=2\pi\alpha_{\mathrm{em}}^{2}\left[1+(1-y)^{2}\right]$. The
structure functions in Eq. (7) can be expressed as the convolution of the
corresponding TMD PDFs and FFs as Bacchetta et al. (2007)

$\displaystyle F_{UU,T}$ $\displaystyle=\mathcal{C}\left[f_{1}D_{1}\right],$
(8) $\displaystyle F_{UT,T}^{\sin\left(\phi_{h}-\phi_{s}\right)}$
$\displaystyle=\mathcal{C}\left[-\frac{\bm{\hat{h}}\cdot\bm{p}_{T}}{M}f_{1T}^{\perp}D_{1}\right].$
(9)

Here, $f_{1}(x,\bm{{p}_{T}})$ is the unpolarized TMD PDF, and
$f_{1T}^{\perp}(x,\bm{p}_{T})$ is the Sivers function. $D_{1}(z,\bm{k}_{T})$
is the unpolarized TMD FF, which depends on the longitudinal momentum fraction
$z$ and the transverse momentum $\bm{k}_{T}$ of the final-state quark.
$\hat{h}=\frac{\bm{P_{hT}}}{|P_{hT}|}$ is the unit vector along $P_{hT}$. The
notation $\mathcal{C}$ represents the convolution among the transverse momenta

$\mathcal{C}\bigl{[}\omega fD\bigr{]}=x\sum_{q}e_{q}^{2}\int
d^{2}\bm{p}_{T}d^{2}\bm{k}_{T}\delta^{(2)}\bigl{(}\bm{p}_{T}-\bm{k}_{T}-\bm{P}_{hT}/z\bigr{)}\omega(\bm{p}_{T},\bm{k}_{T})f^{q}(x,p_{T}^{2})\,D^{q}(z,k_{T}^{2}).$
(10)

Substituting Eq. (10) into Eq. (8), we can expand the unpolarized structure
function $F_{UU}$ as

$\displaystyle F_{UU}\left(Q;P_{hT}\right)$
$\displaystyle=\mathcal{C}\left[f_{1}D_{1}\right]$
$\displaystyle=x\sum_{q}e_{q}^{2}\int
d^{2}\bm{p}_{T}d^{2}\bm{k}_{T}\delta^{(2)}\left(\bm{p}_{T}-\bm{k}_{T}-\bm{P}_{hT}/z\right)f_{1}^{q}\left(x,p_{T}^{2}\right)D_{1}^{q}\left(z,k_{T}^{2}\right)$
$\displaystyle=\frac{x}{z^{2}}\sum_{q}e_{q}^{2}\int
d^{2}\bm{p}_{T}d^{2}\bm{K}_{T}\delta^{(2)}\left(\bm{p}_{T}+\bm{K}_{T}/z-\bm{P}_{hT}/z\right)f_{1}^{q}\left(x,p_{T}^{2}\right)D_{1}^{q}\left(z,\frac{K_{T}^{2}}{z^{2}}\right)$
$\displaystyle=\frac{x}{z^{2}}\sum_{q}e_{q}^{2}\int
d^{2}\bm{p}_{T}d^{2}\bm{K}_{T}\int\frac{d^{2}b}{(2\pi)^{2}}e^{-i\left(p_{T}+\bm{K}_{T}/z-\bm{P}_{hT}/z\right)\cdot\bm{b}}f_{1}^{q}\left(x,p_{T}^{2}\right)D_{1}^{q}\left(z,\frac{K_{T}^{2}}{z^{2}}\right)$
$\displaystyle=\frac{x}{z^{2}}\sum_{q}e_{q}^{2}\int\frac{d^{2}b}{(2\pi)^{2}}e^{i\bm{P}_{hT}\cdot\bm{b}/z}\tilde{f}_{1}^{q/p}\left(x,b\right)\tilde{D}_{1}^{h/q}\left(z,b\right),$
(11)

where $\bm{K}_{T}$ represents the transverse momentum of the final state
hadron with respect to the fragmentation quark, which has the relation
$\bm{K}_{T}=-z\bm{k}_{T}$ with $\bm{k_{T}}$ being the final-state quark
transverse momentum respect to $z$ axis. The $\delta$-function Fourier
Transformation was performed in the fourth line. The TMD distribution function
$\tilde{f_{1}}(x,b)$ and TMD fragmentation function $\tilde{D_{1}}(z,b)$ in
the $b$ space can be obtained by performing a Fourier Transformation from
momentum space to $b$ space

$\displaystyle\int
d^{2}\bm{p}_{T}e^{-i\bm{p}_{T}\cdot\bm{b}}f_{1}^{q}(x,p_{T}^{2})$
$\displaystyle=\tilde{f}_{1}^{q/p}(x,b),$ (12) $\displaystyle\int
d^{2}\bm{K}_{T}e^{-i\bm{K}_{T}/z\cdot\bm{b}}D_{1}^{q}(z,K_{T}^{2})$
$\displaystyle=\tilde{D}_{1}^{h/q}\left(z,b\right),$ (13)

hereafter, the term with a tilde denotes it is in the $b$ space. Similarly,
substituting Eq. (10) into Eq. (9), we can obtain expansions for the spin-
dependent structure function $F_{UT}^{\sin\left(\phi_{h}-\phi_{s}\right)}$ as

$\displaystyle
F_{UT,T}^{\sin\left(\phi_{h}-\phi_{S}\right)}\left(Q;P_{hT}\right)$
$\displaystyle=\mathcal{C}\left[-\frac{\hat{\bm{h}}\cdot\bm{p}_{T}}{M}f_{1T}^{\perp}D_{1}\right]$
$\displaystyle=x\sum_{q}e_{q}^{2}\int
d^{2}\bm{p}_{T}d^{2}\bm{k}_{T}\delta^{(2)}\left(\bm{p}_{T}-\bm{k}_{T}-\bm{P}_{hT}/z\right)\left[-\frac{\hat{\bm{h}}\cdot\bm{p}_{T}}{M}f_{1T}^{\perp}\left(x,p_{T}^{2}\right)D_{1}^{q}\left(z,k_{T}^{2}\right)\right]$
$\displaystyle=\frac{x}{z^{2}}\sum_{a}e_{q}^{2}\int
d^{2}\bm{p}_{T}d^{2}\bm{K}_{T}\delta^{(2)}\left(\bm{p}_{T}+\bm{K}_{T}/z-\bm{P}_{hT}/z\right)\left[-\frac{\hat{\bm{h}}\cdot\bm{p}_{T}}{M}f_{1T}^{\perp}\left(x,p_{T}^{2}\right)D_{1}^{q}\left(z,\frac{K_{T}^{2}}{z^{2}}\right)\right]$
$\displaystyle=\frac{x}{z^{2}}\sum_{q}e_{q}^{2}\int
d^{2}\bm{p}_{T}d^{2}\bm{K}_{T}\int\frac{d^{2}b}{(2\pi)^{2}}e^{-i\left(p_{T}+\bm{K}_{T}/z-\bm{P}_{hT}/z\right)\cdot
b}\left[-\frac{\hat{\bm{h}}\cdot\bm{p}_{T}}{M}f_{1T}^{\perp}\left(x,p_{T}^{2}\right)D_{1}^{q}\left(z,\frac{K_{T}^{2}}{z^{2}}\right)\right]$
$\displaystyle=-\frac{x}{2z^{2}}\sum_{q}e_{q}^{2}\int\frac{d^{2}b}{(2\pi)^{2}}e^{i\bm{P}_{h}\cdot\bm{b}/z}i\hat{h}_{\alpha}b^{\alpha}T_{q,F}\left(x,x\right)\tilde{D}_{1}^{h/q}\left(z,b\right).$
(14)

The Sivers function in $b$ space can also be obtained from momentum space
$f_{1}^{\perp}(x,\bm{p}_{T})$ to $b$ space by Fourier Transformation as

$\displaystyle f_{1T}^{\perp q(\alpha)}(x,b)$ $\displaystyle=\frac{1}{M}\int
d^{2}p_{\perp}e^{-ip_{\perp}\cdot b}p_{\perp}^{\alpha}f_{1T}^{\perp
q}\left(x,p_{\perp}^{2}\right)=\frac{ib^{\alpha}}{2}T_{q,F}\left(x,x\right),$
(15)

where $T_{q,F}\left(x,x\right)$ is the Qiu-Sterman (QS) function. Therefore,
the Sivers asymmetry can be rewritten as

$\displaystyle
A_{UT}^{\sin\left(\phi_{h}-\phi_{s}\right)}=\frac{\sigma_{0}\left(x,y,Q^{2}\right)}{\sigma_{0}\left(x,y,Q^{2}\right)}\frac{-\frac{x}{2z^{2}}\sum_{q}e_{q}^{2}\int\frac{d^{2}b}{(2\pi)^{2}}e^{i\bm{P}_{h}\cdot\bm{b}/z}i\hat{\bm{h}}\cdot
b^{\alpha}T_{q,F}\left(x,x\right)\tilde{D}_{1}^{h/q}\left(z,b\right)}{\frac{x}{z^{2}}\sum_{q}e_{q}^{2}\int\frac{d^{2}b}{(2\pi)^{2}}e^{i\bm{P}_{hT}\cdot\bm{b}/z}\tilde{f}_{1}^{q/p}\left(x,b\right)\tilde{D}_{1}^{h/q}\left(z,b\right)}.$
(16)

One should note that the energy dependence of the TMD structure functions were
not encoded in the above formalism, which will be studied in details in the
following subsections.

###  II.1 TMD evolution effects

First we set up the basic formalism of the TMD evolution effects for TMD PDFs
and FFs, which mainly serves to solve the energy dependence of the TMD PDFs
$f_{1}(x,\bm{{p}_{T}})$, $f_{1T}^{\perp}(x,\bm{p}_{T})$ and the TMD FF
$D_{1}(z,\bm{k}_{T})$. Since the complicated convolution among the transverse
momenta can be transformed into a simple product after performing the Fourier
Transformation, it is convenient to solve the energy dependence in the $b$
space.

Particularly, there are two different energy dependencies $\mu$ and
$\zeta_{F}~{}(\zeta_{D})$ of the TMD PDF $\tilde{F}(x,b)$ and the TMD FF
$\tilde{D}(z,b)$ in $b$ space according to TMD factorization. $\mu$ is the
renormalization scale related to the corresponding collinear PDFs/FFs, and
$\zeta_{F}~{}(\zeta_{D})$ is the energy scale serving as a cutoff to
regularize the light-cone singularity in the operator definition of the TMDs.
The $\mu$ and $\zeta_{F}~{}(\zeta_{D})$ dependencies are encoded in different
TMD evolution equations. The energy evolution for the
$\zeta_{F}~{}(\zeta_{D})$ dependence is encoded in the Collins-Soper (CS)
equationCollins et al. (1985); Collins (2013); Idilbi et al. (2004)

$\displaystyle\frac{\partial\ln\tilde{F}\left(x,b;\mu,\zeta_{F}\right)}{\partial\ln\sqrt{\zeta_{F}}}=\frac{\partial\ln\tilde{D}\left(z,b;\mu,\zeta_{D}\right)}{\partial\ln\sqrt{\zeta_{D}}}=\tilde{K}(b;\mu),$
(17)

while the $\mu$ dependence is given by the renormalization group equation

$\displaystyle\frac{d\tilde{K}}{d\ln\mu}=-\gamma_{K}\left(\alpha_{s}(\mu)\right),$
(18)
$\displaystyle\frac{d\ln\tilde{F}\left(x,b;\mu,\zeta_{F}\right)}{d\ln\mu}=\gamma_{F}\left(\alpha_{s}(\mu);\frac{\zeta_{F}^{2}}{\mu^{2}}\right),$
(19)
$\displaystyle\frac{d\ln\tilde{D}\left(z,b;\mu,\zeta_{D}\right)}{d\ln\mu}=\gamma_{D}\left(\alpha_{s}(\mu);\frac{\zeta_{D}^{2}}{\mu^{2}}\right),$
(20)

with $\alpha_{s}$ being the running strong coupling at the energy scale $\mu$,
$\tilde{K}$ being the CS evolution kernel, and $\gamma_{K},\gamma_{F}$ and
$\gamma_{D}$ being the anomalous dimensions. Hereafter, we will assume
$\mu=\sqrt{\zeta_{F}}=\sqrt{\zeta_{D}}=Q$, then the TMD PDFs and FFs can be
written as $\tilde{F}(x,b;Q)$ and $\tilde{D}(z,b;Q)$ for simplicity.

Solving these TMD evolution equations, one can obtain the solution of the
energy dependence for TMD parton distribution functions and the fragmentation
functions, of which the solution has the general form as

$\displaystyle\tilde{F}_{q/p}(x,b;Q)=\mathcal{F}\times
e^{-S}\times\tilde{F}_{q/p}(x,b;\mu_{B}),$ (21)
$\displaystyle\tilde{D}_{h/q}(z,b;Q)=\mathcal{D}\times
e^{-S}\times\tilde{D}_{h/q}(z,b;\mu_{B}),$ (22)

where $\mathcal{F}$ and $\mathcal{D}$ is the factor related to the hard
scattering and depend on the factorization schemes, $S$ is the Sudakov-like
form factor. Eq. (21) and Eq. (22) show that the energy evolution of TMD PDFs
and TMD fragmentation functions from an initial energy $\mu_{B}$ to another
energy $Q$ is encoded in the Sudakov-like form factor $S$ by the exponential
form $\mathrm{exp}(-S)$.

By performing the reverse Fourier Transformation of the TMDs in $b$ space, the
TMDs in momentum space can be obtained, thus it is of great importance to
study the $b$ space behavior of the TMDs. In the small $b$ region ($b\ll
1/\Lambda_{\textrm{QCD}}$), the $b$ dependence of TMDs is perturbative and can
be calculated by perturbative QCD. However, the dependence in large $b$ region
turns to be nonperturbative, since the operators are separated by a large
distance. To include the evolution effect in this region, a nonperturbative
Sudakov-like form factor $S_{\rm NP}$ is introduced and is usually given in a
parameterized form. The parameters of $S_{\rm NP}$ can be determined by
analyzing experimental data, given the lack of non-perturbative calculations.
In order to combine the information from both the small $b$ region and the
large $b$ region, a matching procedure is applied with a parameter
$b_{\mathrm{max}}$ serving as the boundary between the two regions.
Furthermore, one can define a $b$-dependent function $b_{\ast}$, which has the
property $b_{\ast}\approx b$ in small $b$ region and $b_{\ast}\approx
b_{\mathrm{max}}$ in large $b$ region

$\displaystyle b_{\ast}=\frac{b}{\sqrt{1+b^{2}/b_{\rm
max}^{2}}}\,\,,\,\,b_{\rm max}<1/\Lambda_{\textrm{QCD}}\,,$ (23)

as given in the original CSS prescription Collins et al. (1985). The
prescription also allows for a smooth transition from perturbative to
nonperturbative regions and avoids the Landau pole singularity in
$\alpha_{s}(\mu_{B})$. The typical value of $b_{\mathrm{max}}$ is chosen
around 1.5 GeV-1 to guarantee that $b_{\ast}$ is always in the perturbative
region. With the constraint of $b_{\ast}$, we can calculate TMDs within a
small $b$ region.

In the small $b$ region, the TMDs can be expressed as the convolutions of the
perturbatively calculable hard coefficients and the corresponding collinear
counterparts at fixed energy $\mu_{B}$, which could be the collinear PDFs/FFs
or the multiparton correlation functions Collins and Soper (1981); Bacchetta
and Prokudin (2013)

$\displaystyle\tilde{F}_{q/p}(x,b;\mu_{B})=C_{q\leftarrow i}\otimes
F_{i/p}(x,\mu_{B}),$ (24)
$\displaystyle\tilde{D}_{h/q}\left(z,b;\mu_{B}\right)=\hat{C}_{j\leftarrow
q}\otimes D_{h/j}\left(z,\mu_{B}\right),$ (25)

where $\mu_{B}=c_{0}/b_{*}$ and $c_{0}=2e^{-\gamma_{E}}$ and the Euler
constant $\gamma_{E}\approx 0.577$ Collins and Soper (1981), the $\otimes$
stands for the convolution in the momentum fraction $x$

$\displaystyle C_{q\leftarrow i}\otimes F_{i/p}(x,\mu_{B})$
$\displaystyle\equiv\sum_{i}\int_{x}^{1}\frac{d\xi}{\xi}C_{q\leftarrow
i}\left(\frac{x}{\xi},\mu_{B}\right)F_{i/p}(\xi,\mu_{B}),$ (26)
$\displaystyle\hat{C}_{j\leftarrow q}\otimes D_{h/j}(z,\mu_{B})$
$\displaystyle\equiv\sum_{j}\int_{z}^{1}\frac{d\xi}{\xi}\hat{C}_{j\leftarrow
q}\left(\frac{z}{\xi},\mu_{B}\right)D_{h/j}\left(\xi,\mu_{B}\right),$ (27)

where $C$ coefficients in the formula has different values in different
processes, and its specific value will be given in the subsequent calculation.
In addition, the sum $\sum_{i}$ runs over all parton flavors. Now we can
combine all the above information to get the expression for TMD distribution
function and the fragmentation function in $b$ space as

$\displaystyle\tilde{F}_{q/p}(x,b;Q)$ $\displaystyle=\mathcal{F}\times
e^{-S}\times C_{q\leftarrow i}\otimes F_{i/p}(x,\mu_{B})$
$\displaystyle=\mathcal{F}\times
e^{-S}\times\sum_{i}\int_{x}^{1}\frac{d\xi}{\xi}C_{q\leftarrow
i}\left(\frac{x}{\xi},\mu_{B}\right)F_{i/p}(\xi,\mu_{B}),$ (28)
$\displaystyle\tilde{D}_{h/q}(z,b;Q)$ $\displaystyle=\mathcal{D}\times
e^{-S}\times\hat{C}_{q\leftarrow i}\otimes D_{h/j}(z,\mu_{B})$
$\displaystyle=\mathcal{D}\times
e^{-S}\times\sum_{j}\int_{z}^{1}\frac{d\xi}{\xi}\hat{C}_{j\leftarrow
q}\left(\frac{z}{\xi},\mu_{B}\right)D_{h/j}\left(\xi,\mu_{B}\right).$ (29)

The Sudakov-like form factor $S$ can be separated into the perturbatively
calculable part ${S}_{\rm pert}(Q;b_{*})$ and the nonperturbative part $S_{\rm
NP}(Q;b)$

${S}(Q;b)={S}_{\rm pert}(Q;b_{*})+S_{\rm NP}(Q;b).$ (30)

The perturbative part ${S}_{\rm pert}(Q;b_{*})$ has a general form and can be
expanded as the series of $({\alpha_{s}\over\pi})$ Echevarria et al. (2013);
Echevarria et al. (2014a); Kang et al. (2011); Aybat et al. (2012); Echevarria
et al. (2014b).

${S}_{\rm
pert}(Q;b_{*})=\int^{Q^{2}}_{\mu_{b}^{2}}\frac{d\bar{\mu}^{2}}{\bar{\mu}^{2}}\left[A(\alpha_{s}(\bar{\mu}))\ln\left(\frac{Q^{2}}{\bar{\mu}^{2}}\right)+B(\alpha_{s}(\bar{\mu}))\right].$
(31)

In Eq. (31) the coefficients $A$ and $B$ can be expanded as following

$\displaystyle
A=\sum_{n=1}^{\infty}A^{(n)}\left(\frac{\alpha_{s}}{\pi}\right)^{n},$ (32)
$\displaystyle
B=\sum_{n=1}^{\infty}B^{(n)}\left(\frac{\alpha_{s}}{\pi}\right)^{n}.$ (33)

In our calculation we take $A^{(n)}$ up to $A^{(2)}$ and $B^{(n)}$ up to
$B^{(1)}$Collins et al. (1985); Aybat and Rogers (2011); Echevarria et al.
(2013); Kang et al. (2011); Landry et al. (2003); Qiu and Zhang (2001),

$\displaystyle A^{(1)}$ $\displaystyle=C_{F},$ (34) $\displaystyle A^{(2)}$
$\displaystyle=\frac{C_{F}}{2}\left[C_{A}\left(\frac{67}{18}-\frac{\pi^{2}}{6}\right)-\frac{10}{9}T_{R}n_{f}\right],$
(35) $\displaystyle B^{(1)}$ $\displaystyle=-\frac{3}{2}C_{F}.$ (36)

For the non-perturbative Sudakov-like form factor $S_{\rm NP}(Q;b)$, it cannot
be obtained from perturbation calculation, and it is usually extracted from
the experimental data. Inspired by Refs. Landry et al. (2003); Konychev and
Nadolsky (2006), a widely used parametrization (the EIKV parametrization) of
$S_{\rm NP}$ for TMD PDFs or fragmentation functions was proposed Bacchetta
and Prokudin (2013); Echevarria et al. (2014a); Landry et al. (2003); Konychev
and Nadolsky (2006); Davies et al. (1984); Ellis et al. (1997)

$\displaystyle
S_{\mathrm{NP}}^{\mathrm{pdf}/\mathrm{ff}}=b^{2}\left(g_{1}^{\mathrm{pdf}/\mathrm{ff}}+\frac{g_{2}}{2}\ln\frac{Q}{Q_{0}}\right),$
(37)

where the initial energy is $Q^{2}_{0}=2.4~{}\mathrm{GeV}^{2}$, and the factor
$1/2$ in front of $g_{2}$ comes from the fact that only one hadron is involved
for the parametrization of $S_{\rm NP}^{\rm pdf/ff}$. The parameter
$g_{1}^{\rm pdf/ff}$ in Eq. (37) depends on the type of TMDs, which can be
regarded as the width of the intrinsic transverse momentum for the relevant
TMDs at the initial energy scale $Q_{0}$ Aybat and Rogers (2011); Qiu and
Zhang (2001); Anselmino et al. (2012). Assuming a Gaussian form for the
dependence for the transverse momentum, one can obtain

$\displaystyle g_{1}^{\rm pdf}=\frac{\langle
p_{\perp}^{2}\rangle_{Q_{0}}}{4},\qquad g_{1}^{\rm ff}=\frac{\langle
k_{\perp}^{2}\rangle_{Q_{0}}}{4z_{h}^{2}},$ (38)

where $\langle p_{\perp}^{2}\rangle_{Q_{0}}$ and $\langle
k_{\perp}^{2}\rangle_{Q_{0}}$ represent the averaged intrinsic transverse
momenta squared for TMD PDFs and FFs at the initial scale $Q_{0}$,
respectively. $g_{2}=0.16$ Echevarria et al. (2014a) has the same form in all
forms of TMDs.

As the information on the Sudakov-like form factor for the Kaon fragmentation
functions can not be determined, we assume that the evolution of the TMD
distribution function and the fragmentation function for producing the $K$
meson from the initial energy scale $\mu$ to another energy scale $Q$ follows
the Gaussian form of $g_{2}(b)$ in Eq. (37), so we can obtain the
nonperturbative Sudakov-like form factor for the PDF and FF in the produced of
the $K$ mesons as

$\displaystyle S_{\mathrm{NP}}^{\mathrm{pdf}}(Q;b)$
$\displaystyle=\frac{g_{2}}{2}\ln\left(\frac{Q}{Q_{0}}\right)b^{2}+g_{1}^{\mathrm{pdf}}b^{2},$
(39) $\displaystyle S_{\mathrm{NP}}^{\mathrm{ff}}(Q;b)$
$\displaystyle=\frac{g_{2}}{2}\ln\left(\frac{Q}{Q_{0}}\right)b^{2}+g_{1}^{\mathrm{ff}}b^{2}.$
(40)

Combining all the steps mentioned above together, the scale-dependent TMD PDFs
and FFs in $b$ space as functions of $x$ (or $z$), $b$, and $Q$ can be
rewritten as

$\displaystyle\tilde{F}_{q/p}(x,b;Q)$
$\displaystyle=e^{-\frac{1}{2}S_{\mathrm{Pert}}(Q;b_{\ast})-S^{\rm
pdf}_{\mathrm{NP}}(Q;b)}\mathcal{F}(\alpha_{s}(Q))\sum_{i}\int_{x}^{1}\frac{d\xi}{\xi}C_{q\leftarrow
i}\left(\frac{x}{\xi},\mu_{B}\right)F_{i/p}(\xi,\mu_{B}),$ (41)
$\displaystyle\tilde{D}_{h/q}(z,b;Q)$
$\displaystyle=e^{-\frac{1}{2}S_{\mathrm{Pert}}(Q;b_{\ast})-S^{\rm
ff}_{\mathrm{NP}}(Q;b)}\mathcal{D}(\alpha_{s}(Q))\sum_{j}\int_{z}^{1}\frac{d\xi}{\xi}\hat{C}_{j\leftarrow
q}\left(\frac{z}{\xi},\mu_{B}\right)D_{h/j}\left(\xi,\mu_{B}\right).$ (42)

###  II.2 The solution of the unpolarized structure function

In the following we solve the denominator of the Sivers asymmetry in details,
which is the unpolarized structure function $F_{UU}$ in Eq. (8). Since we have
expanded $F_{UU}$ in Eq. (11), we will directly follow Eq. (11) to give the
complete expression for $F_{UU}$. We can divide Eq. (11) into three parts to
solve the problem: the weighted Bessel function, the distribution function
$\tilde{f}_{1}^{q/p}(x,b)$ and the fragmentation function
$\tilde{D}_{1}^{h/q}(z,b)$. The weighted Bessel function can be obtained as

$\displaystyle\frac{x}{z^{2}}\sum_{q}e_{q}^{2}\int\frac{d^{2}b}{(2\pi)^{2}}e^{i\bm{P_{hT}}\cdot\bm{b}/z}$
$\displaystyle=\frac{x}{z^{2}}\sum_{q}e_{q}^{2}\int_{0}^{\infty}\frac{dbb}{2\pi}J_{0}\left(\frac{P_{hT}b}{z}\right).$
(43)

According to Eq. (41) and Eq. (42), $\tilde{f}_{1}^{q}(x,b)$ and
$\tilde{D}_{1}^{h/q}(z,b)$ can be expanded as

$\displaystyle\tilde{f}_{1}^{q/p}(x,b;Q)$
$\displaystyle=e^{-\frac{1}{2}S_{\mathrm{Pert}}(Q;b_{\ast})-S^{f_{1}}_{\mathrm{NP}}(Q;b)}\mathcal{F}(\alpha_{s}(Q))\sum_{i}\int_{x}^{1}\frac{d\xi}{\xi}C_{q\leftarrow
i}\left(\frac{x}{\xi},\mu_{B}\right)f_{1}^{i/p}(\xi,\mu_{B}),$ (44)
$\displaystyle\tilde{D}_{1}^{h/q}(z,b;Q)$
$\displaystyle=e^{-\frac{1}{2}S_{\mathrm{Pert}}(Q;b_{\ast})-S^{D_{1}}_{\mathrm{NP}}(Q;b)}\mathcal{D}(\alpha_{s}(Q))\sum_{j}\int_{z}^{1}\frac{d\xi}{\xi}\hat{C}_{j\leftarrow
q}\left(\frac{z}{\xi},\mu_{B}\right)D_{1}^{h/j}\left(\xi,\mu_{B}\right),$ (45)

where the Sudakov-like form factor is expressed in Eq. (31), Eq. (39) and Eq.
(40). The hard scattering coefficients $\mathcal{F}(\alpha_{s}(Q))$ and
$\mathcal{D}(\alpha_{s}(Q))$ can be set equal to 1. In addition, the $C$
coefficient in Eq. (44) and Eq. (45) can be expressed as Kang et al. (2016)

$\displaystyle C_{q\leftarrow q^{\prime}}^{\rm(SIDIS)}(x,\mu_{B})$
$\displaystyle=\delta_{q^{\prime}q}[\delta(1-x)+\frac{\alpha_{s}}{\pi}(\frac{C_{F}}{2}(1-x)-2C_{F}\delta(1-x))]\;,$
(46) $\displaystyle C_{q\leftarrow g}^{\rm(SIDIS)}(x,\mu_{B})$
$\displaystyle=\frac{\alpha_{s}}{\pi}{T_{R}}\,x(1-x)\;,$ (47)
$\displaystyle\hat{C}_{q^{\prime}\leftarrow q}^{\rm(SIDIS)}(z,\mu_{B})$
$\displaystyle=\delta_{q^{\prime}q}[\delta(1-z)+\frac{\alpha_{s}}{\pi}(\frac{C_{F}}{2}(1-z)-2C_{F}\delta(1-z)+P_{q\leftarrow
q}(z)\,\ln z)]\;,$ (48) $\displaystyle\hat{C}_{g\leftarrow
q}^{\rm(SIDIS)}(z,\mu_{B})$
$\displaystyle=\frac{\alpha_{s}}{\pi}\left(\frac{C_{F}}{2}z\;+P_{g\leftarrow
q}(z)\,\ln z\right),$ (49)

where $\alpha_{s}$ represents the strong coupling, and the expansion at the
next-to-leading order can be written as follows

$\displaystyle\alpha_{s}\left(Q^{2}\right)=\frac{12\pi}{\left(33-2n_{f}\right)\ln\left(Q^{2}/\Lambda_{QCD}^{2}\right)}\left\\{1-\frac{6\left(153-19n_{f}\right)}{\left(33-2n_{f}\right)^{2}}\frac{\ln\ln\left(Q^{2}/\Lambda_{QCD}^{2}\right)}{\ln\left(Q^{2}/\Lambda_{QCD}^{2}\right)}\right\\}.$
(50)

In Eq. (50), $Q^{2}$ is the running energy scale and $n_{f}=5$,
$\Lambda_{QCD}=0.225$ GeV. The splitting functions $P_{q\leftarrow q}$ and
$P_{g\leftarrow q}$ in Eq. (48) and Eq. (49) have the general form

$\displaystyle P_{q\leftarrow q}(z)$
$\displaystyle=C_{F}\left[\frac{1+z^{2}}{(1-z)_{+}}+\frac{3}{2}\delta(1-z)\right],$
(51) $\displaystyle P_{g\leftarrow q}(z)$
$\displaystyle=C_{F}\frac{1+(1-z)^{2}}{z},$ (52)

where $C_{F}=4/3$, $T_{R}=1/2$, and the subscript symbol “$+$” denotes the
following prescription

$\displaystyle\int_{0}^{1}dz{f(z)\over(1-z)_{+}}=\int_{0}^{1}dz{f(z)-f(1)\over(1-z)}$
(53)

Combining the Bessel function part, the PDF part $\tilde{f}_{1}^{q/p}(x,b)$
and the FF part $\tilde{D}_{1}^{q}(z,b)$, we can expand the denominator of
Sivers asymmetry as

$\displaystyle F_{UU}(Q;P_{hT})$
$\displaystyle=\frac{x}{z^{2}}\sum_{q}e_{q}^{2}\int_{0}^{\infty}\frac{bdb}{(2\pi)}J_{0}\left(\frac{P_{hT}b}{z}\right)e^{-{S}_{\rm
pert}(Q;b_{*})-S_{\rm NP}^{\rm SIDIS}(Q;b)}$
$\displaystyle~{}~{}\left(\sum_{i}C\int_{x}^{1}\frac{d\xi}{\xi}C_{q\leftarrow
i}^{\rm(SIDIS)}\left(\frac{x}{\xi},\mu_{B}\right)f_{1}^{i/p}(\xi,\mu_{B})\right)\times\left(\sum_{j}\int_{z}^{1}\frac{d\xi}{\xi}\hat{C}_{j\leftarrow
q}^{\rm(SIDIS)}\left(\frac{z}{\xi},\mu_{B}\right)D_{1}^{h/j}\left(\xi,\mu_{B}\right)\right),$
(54)

where the nonperturbative Sudakov-like form factor
$S_{\mathrm{NP}}^{\operatorname{SIDIS}}(Q;b)$ receives the contribution from
the unpolarized TMD PDF and FF:

$\displaystyle S_{\mathrm{NP}}^{\operatorname{SIDIS}}(Q;b)$
$\displaystyle=S_{\mathrm{NP}}^{\mathrm{f}_{1}}(Q;b)+S_{\mathrm{NP}}^{\mathrm{D}_{1}}(Q;b)$
$\displaystyle=g_{2}\ln\left(\frac{Q}{Q_{0}}\right)b^{2}+g_{1}^{\mathrm{f}_{1}}b^{2}+g_{1}^{\mathrm{D}_{1}}b^{2}.$
(55)

Here, $g_{1}^{\mathrm{f}_{1}}$ and $g_{1}^{\mathrm{D}_{1}}$ are obtained from
Eq. (38), so $g_{1}^{\mathrm{f}_{1}}=g_{1}^{\rm pdf}=\frac{\langle
p_{\perp}^{2}\rangle_{Q_{0}}}{4}$, $g_{1}^{\mathrm{D}_{1}}=g_{1}^{\rm
ff}=\frac{\langle k_{\perp}^{2}\rangle_{Q_{0}}}{4z_{h}^{2}}$.

###  II.3 The solution of the transverse spin-dependent structure function

Similar to the previous subsection, we further obtain the expression for the
transverse spin-dependent structure function
$F_{UT}^{\sin(\phi_{h}-\phi_{s})}$. Eq. (14) can be separated to three parts
to solve it, namely, the weighted Bessel function, the Qiu-Sterman (QS)
function $T_{q,F}\left(x^{\prime},x^{\prime\prime}\right)$ and the
fragmentation function $\tilde{D}_{1}^{h/q}(z,b)$. The solution for the
weighted Bessel function is similar to Eq. (43)

$\displaystyle-\frac{x}{2z^{2}}\sum_{q}e_{q}^{2}\int\frac{d^{2}b}{(2\pi)^{2}}i\hat{h}_{\alpha}b^{\alpha}e^{i\bm{P_{hT}}\cdot\bm{b}/z}$
$\displaystyle=\frac{x}{2z^{2}}\sum_{q}e_{q}^{2}\int_{0}^{\infty}\frac{b^{2}db}{2\pi}J_{1}\left(\frac{P_{hT}b}{z}\right).$
(56)

The QS function part $T_{q,F}\left(x^{\prime},x^{\prime\prime}\right)$ in the
$b$ space has a similar solution form after solving the TMD evolution
equations

$\displaystyle T_{q,F}\left(x,x;Q\right)=\mathcal{T}\times e^{-S}\times
T_{q,F}\left(x,x;\mu_{B}\right).$ (57)

In the small $b$ region, the QS function can be calculated from the
convolution of the hard scattering coefficients and the collinear counterpart
utilizing the perturbative QCD

$\displaystyle T_{q,F}\left(x,x;\mu_{B}\right)=\Delta C_{q\leftarrow
i}^{T}\otimes T_{i,F}\left(x,x;\mu_{B}\right),$ (58)

where $\otimes$ denotes convolution and can be expended as

$\displaystyle\Delta C_{q\leftarrow i}^{T}\otimes
T_{i,F}\left(x,x;\mu_{B}\right)=\sum_{i}\int_{x}^{1}\frac{d\xi}{\xi}\Delta
C_{q\leftarrow
i}^{T}\left(\frac{x}{\xi},\mu_{B}\right)T_{i,F}\left(\xi,\xi;\mu_{B}\right),$
(59)

thus the QS function $T_{q,F}\left(x^{\prime},x^{\prime\prime}\right)$ and the
fragmentation function $\tilde{D}_{1}^{h/q}(z,b)$ can be expanded as

$\displaystyle
T_{q,F}\left(x,x;Q\right)=e^{-\frac{1}{2}S_{\mathrm{Pert}}(Q;b_{\ast})-S^{f_{1}}_{\mathrm{NP}}(Q;b)}\mathcal{T}(\alpha_{s}(Q))\sum_{i}\int_{x}^{1}\frac{d\xi}{\xi}\Delta
C_{q\leftarrow
i}^{T}\left(\frac{x}{\xi},\mu_{B}\right)T_{i,F}\left(\xi,\xi;\mu_{B}\right),$
(60)
$\displaystyle\tilde{D}_{1}^{h/q}(z,b;Q)=e^{-\frac{1}{2}S_{\mathrm{Pert}}(Q;b_{\ast})-S^{D_{1}}_{\mathrm{NP}}(Q;b)}\mathcal{D}(\alpha_{s}(Q))\sum_{j}\int_{z}^{1}\frac{d\xi}{\xi}\hat{C}_{j\leftarrow
q}\left(\frac{z}{\xi},\mu_{B}\right)D_{1}^{h/j}\left(\xi,\mu_{B}\right),$ (61)

where the Sudakov-like form factor is expressed in Eq. (31), Eq. (39) and Eq.
(40). And the hard scattering coefficient $\mathcal{T}(\alpha_{s}(Q))$ and
$\mathcal{D}(\alpha_{s}(Q))$ equal to 1. The $C$ coefficient o
$\tilde{D}_{1}^{q}(z,b)$ is expressed in Eqs. (48), (49), (51) and (52). The
expression of $\Delta C$ in Eq. (60) is Sun and Yuan (2013)

$\displaystyle\Delta
C_{qq}^{T(DIS)}(x,\mu_{B})=\delta(1-x)+\frac{\alpha_{s}(\mu_{B})}{\pi}\left(-\frac{1}{4N_{c}}(1-x)-2C_{F}\delta(1-x)\right),$
(62)

where $N_{c}=3$. Combining all the ingredients together, the numerator part of
the Sivers asymmetry can be rewritten as

$\displaystyle F_{UT,T}^{\sin\left(\phi_{h}-\phi_{S}\right)}=$
$\displaystyle-\frac{x}{2z^{2}}\sum_{q}e_{q}^{2}\int\frac{d^{2}b}{(2\pi)^{2}}e^{i\bm{P}_{hT}/z\cdot\bm{b}}\hat{h}_{\alpha}\cdot
b^{\alpha}T_{q,F}\left(x,x\right)\tilde{D}_{1}^{h/q}(z,b)$ $\displaystyle=$
$\displaystyle\frac{x}{2z^{2}}\sum_{q}e_{q}^{2}\int_{0}^{\infty}\frac{b^{2}db}{2\pi}J_{1}\left(\frac{P_{hT}b}{z}\right)e^{-S_{\text{pert
}}\left(Q;b_{*}\right)-S_{\mathrm{NP\ Sivers}}^{\text{SIDIS }}(Q;b)}$
$\displaystyle\left(\sum_{i}\int_{x}^{1}\frac{d\xi}{\xi}\Delta C_{q\leftarrow
i}^{T(DIS)}\left(\frac{x}{\xi},\mu_{B}\right)T_{i,F}\left(\xi,\xi;\mu_{B}\right)\right)\times\left(\sum_{j}\int_{z}^{1}\frac{d\xi}{\xi}\hat{C}_{j\leftarrow
q}^{(SIDIS)}\left(\frac{z}{\xi},\mu_{B}\right)D_{1}^{h/j}\left(\xi,\mu_{B}\right)\right),$
(63)

where the nonperturbative Sudakov-like form factor $S_{\mathrm{NP\
Sivers}}^{\operatorname{SIDIS}}(Q;b)$ is the combination of and the one for
the Sivers function the one for the unpolarized TMD FF

$\displaystyle S_{\mathrm{NP\ Sivers}}^{\operatorname{SIDIS}}(Q;b)$
$\displaystyle=S_{\mathrm{NP}}^{Sivers}(Q;b)+S_{\mathrm{NP}}^{\mathrm{D}_{1}}(Q;b)$
$\displaystyle=g_{2}\ln\left(\frac{Q}{Q_{0}}\right)b^{2}+g_{1}^{Sivers}b^{2}+g_{1}^{\mathrm{D}_{1}}b^{2},$
(64)

where $g_{1}^{Sivers}=\frac{\langle k_{s\perp}^{2}\rangle}{4}$,
$g_{1}^{\mathrm{D}_{1}}=\frac{\langle
k_{\perp}^{2}\rangle_{Q_{0}}}{4z_{h}^{2}}$.

##  III NUMERICAL ESTIMATE

In this section, we present the numerical estimate for the Sivers asymmetry in
charged Kaon and $\Lambda$ hyperon produced SIDIS process at the kinematical
configurations of EIC and EicC.

In order to obtain the results for the Sivers asymmetry, we shall have the
collinear distribution functions and fragmentation functions as input. For the
unpolarized proton collinear distribution function $f_{1}(x,\mu_{B})$, we
apply the parametrization from MSTW2008 Martin et al. (2009). For the
collinear unpolarized fragmentation function $D_{1}(z)$ of charged Kaon, the
DSS parameterization results at NLO accuracy is used de Florian et al. (2007).
For the collinear unpolarized fragmentation function $D_{1}^{\Lambda}(z)$ of
$\Lambda$ hyperon, we adopt the model results from the diquark spectator model
Yang et al. (2017)

$\displaystyle D_{1}^{\Lambda}(z)=$
$\displaystyle\frac{g_{s}^{2}}{4(2\pi)^{2}}\frac{e^{-\frac{2m_{q}^{2}}{\Lambda^{2}}}}{z^{4}L^{2}}\left\\{z(1-z)\left(\left(m_{q}+M_{\Lambda}\right)^{2}-m_{D}^{2}\right)\right.\times\exp\left(\frac{-2zL^{2}}{(1-z)\Lambda^{2}}\right)$
$\displaystyle+\left((1-z)\Lambda^{2}-2\left(\left(m_{q}+M_{\Lambda}\right)^{2}-m_{D}^{2}\right)\right)\left.\times\frac{z^{2}L^{2}}{\Lambda^{2}}\Gamma\left(0,\frac{2zL^{2}}{(1-z)\Lambda^{2}}\right)\right\\}.$
(65)

The values of the free parameters in Eq. (65) are taken from Ref. Yang et al.
(2017). Since the model result is obtained at the initial energy of 0.23 GeV2,
to make it applicable to a more general energy range, we use the QCDNUM
evolution package Botje (2011) to evolve the unpolarized fragmentation
function $D_{1}(z)$ from the initial energy of 0.23 GeV2 to another energy
scale.

In our calculation we adopt the Qiu-Sterman (QS) function
$T_{q,F}\left(x,x;\mu_{B}\right)$ from the parameterization in Ref. Echevarria
et al. (2014a)

$\displaystyle
T_{q,F}(x,x,\mu_{B})=N_{q}\frac{\left(\alpha_{q}+\beta_{q}\right)^{\left(\alpha_{q}+\beta_{q}\right)}}{\alpha_{q}^{\alpha_{q}}\beta_{q}^{\beta^{q}}}x^{\alpha_{q}}(1-x)^{\beta_{q}}f_{1}(x,\mu_{B}).$
(66)

The parameters obtained by fitting to the data in SIDIS processes from HERMES,
COMPASS and JLab Airapetian et al. (2009); Alekseev et al. (2009); Adolph et
al. (2012); Qian et al. (2011), are listed in Table. 1.

$\alpha_{u}$ | $\alpha_{d}$ | $\alpha_{sea}$ | $\beta$ | $N_{u}$ | $N_{d}$
---|---|---|---|---|---
1.051 | 1.552 | 0.851 | 4.857 | 0.106 | -0.163
$N_{\bar{u}}$ | $N_{\bar{d}}$ | $N_{s}$ | $N_{\bar{s}}$ | $\langle k_{s\perp}^{2}\rangle$ |
-0.012 | -0.105 | 0.103 | -1.000 | 0.282 $\mathrm{GeV}^{2}$ |
Table 1: Values of the parameters in the parameterization of the Qiu-Sterman
function in Ref. Echevarria et al. (2014a)

We should notice that in Ref. Echevarria et al. (2014a), the QS function is
assumed to be proportional to the collinear unpolarized distribution function
$f_{1}(x,\mu_{B})$ in which the DGLAP evolution effect is not included. To
investigate the impact of the QS function DGLAP evolution effect on the Sivers
asymmetry, we take $Q_{0}^{2}=2.4~{}\mathrm{GeV}^{2}$ as the initial energy
and evolve the QS function in Eq. (66) into another energy. We adopt two
different approaches to evolve the QS function: one is to assume the QS
function follows the same evolution effect as that for unpolarized
distribution function, the other one is to change the evolution kernel in the
QCDNUM evolution package to include the QS evolution kernel by considering the
homogenous terms [the terms containing $T_{q,F}(x,x,\mu_{B})$] in the
evolution kernel as an approximation Kang and Qiu (2009):

$\displaystyle
P_{qq}^{\mathrm{Sivers}}\approx\frac{4}{3}\left(\frac{1+z^{2}}{(1-z)_{+}}+\frac{3}{2}\delta(1-z)\right)-\frac{3}{2}\frac{1+z^{2}}{1-z}-3\delta(1-z).$
(67)

In Eq. (37) and Eq. (38), the free parameter $g_{1}$ and the universal
parameter $g_{2}$ contain information about the evolution of TMDs and are the
key parameters that determine the evolution of TMDs from one initial energy
$\mu$ to another $Q$. Here we adopt the results given in Ref. Echevarria et
al. (2014a) for the mean transverse momentum squared

$\left\langle p_{\perp}^{2}\right\rangle=0.38~{}\mathrm{GeV}^{2},\left\langle
k_{\perp}^{2}\right\rangle=0.19~{}\mathrm{GeV}^{2}.$

For the universal parameter $g_{2}$ in the nonperturbative Sudakov-like form
factor, the specific value $g_{2}=0.16$ is also given in Ref. Echevarria et
al. (2014a).

The kinematical region that is available at EIC is chosen as follows Accardi
et al. (2016)

$\displaystyle 0.001<x<0.4,\quad 0.07<y$ $\displaystyle<0.9,\quad
0.2<z<0.8\text{, }$ $\displaystyle 1\mathrm{GeV}^{2}<Q^{2},\quad
W>5\mathrm{GeV},\quad\sqrt{s}$ $\displaystyle=100\mathrm{GeV},\quad
P_{hT}<0.5\mathrm{GeV}.$ (68)

As for the EicC, the following kinematical cuts are adopted

$\displaystyle 0.005<x<0.5,\quad 0.07<$ $\displaystyle y<0.9,\quad 0.2<z<0.7,$
$\displaystyle 1\mathrm{GeV}^{2}<Q^{2}<200\mathrm{GeV}^{2},\quad
W>2\mathrm{GeV}$ $\displaystyle,\quad\sqrt{s}=16.7\mathrm{GeV},\quad
P_{hT}<0.5\mathrm{GeV},$ (69)

where $W^{2}=(P+q)^{2}\approx\frac{1-x}{x}Q^{2}$ is invariant mass of the
virtual photon-nucleon system. Since TMD factorization is proved to be valid
to describe the physical observables in the region $P_{hT}\ll Q$,
$P_{hT}<0.5~{}\mathrm{GeV}$ is chosen to guarantee the validity of TMD
factorization. Combining Eq. (7), Eq. (54) and Eq. (63) and the kinematical
regions of EIC and EicC, we can calculate the single-spin dependent Sivers
asymmetry of charged Kaon and $\Lambda$ hyperon produced SIDIS process within
the EIC and EicC kinematical region.

The results are shown in Fig. 2. The six rows in the figure depict the results
of Sivers asymmetry of $K^{+}$ production, $K^{-}$ production, and $\Lambda$
hyperon production in the EIC and EicC kinematical region, respectively. The
left, central and right panels denote the Sivers asymmetry as the function of
$x$, $z$ and $P_{hT}$, respectively. In the figure, the black solid line
represents the Sivers asymmetry that is obtained when the QS evolution kernel
is included to evolve Eq. (66). The red dashed line represents the Sivers
asymmetry obtained by using the evolution kernel for the unpolarized
distribution function from QCDNUM evolution package to evolve Eq. (66).

Figure 2: The Sivers asymmetry in semi-inclusive charged Kaon and $\Lambda$
hyperon produced SIDIS process at the kinematics of EIC and EicC as functions
of $x$ (left panels), $z$ (middle panels), and $P_{hT}$ (right panels).

As can be seen from Fig. 2, the Sivers asymmetries in the charged Kaon
produced and in $\Lambda$ hyperon produced SIDIS process are sizable in the
EIC and EicC kinematical region. Thus, the measurements of the Sivers
asymmetry in these facilities provide an ideal tool to obtain the sea quark
Sivers function as well as the flavor dependence of the Sivers function. Also,
it seems that the magnitude of the asymmetry in EicC is larger than that in
EIC. In addition, Fig. 2 shows that the value of the solid line is larger than
the dashed line, which indicates that the Sivers asymmetry obtained using the
evolution kernel of QS function is greater than that obtained using the
unpolarized distribution function evolution kernel. One should note that the
DGLAP evolution in the TMD effects may play some role in the future
phenomenological analysis of the Sivers asymmetry.

We can also find from Fig. 2 that both the $z$-dependent and
$P_{hT}$-dependent Sivers asymmetries in $K^{+}$ produced SIDIS process at the
kinematical region of EIC and EicC are negative, and the magnitude of Sivers
asymmetry increases with increasing $z$ or increasing $P_{hT}$. While for $x$
dependent Sivers asymmetry of $K^{+}$ production in the EIC and EicC
kinematical region, there is a node around $x=0.15$. On the other hand, for
the Sivers asymmetry of $K^{-}$ production, it is always positive in all
cases. The Sivers asymmetry increases with $x$ at small $x$ region and
decreases along $x$ at large $x$ region, and gradually increases with $z$ and
$P_{hT}$. In addition, by comparing the magnitudes of the Sivers asymmetry of
$K^{+}$ produced with that of $K^{-}$ produced, we can see that the former one
is significantly larger than the latter one, since for $K^{+}$ meson the
constituent quarks are $u$ and $\bar{s}$, and for $K^{-}$ meson are $\bar{u}$
and $s$ quarks. Therefore there is relatively large contribution of valence
$u$ quark Sivers function for $K^{+}$, while the contributions for $K^{-}$ are
both from sea quark Sivers function. It is known from Ref. Echevarria et al.
(2014a) that the QS function of the valence quark is significantly larger than
the that of the sea quarks, so the Sivers asymmetry in $K^{+}$ produced
process is greater than that in $K^{-}$ produced process. Finally, the Sivers
asymmetry of $\Lambda$ hyperon production in both the EIC and EicC kinematics
is positive, and the tendency of $x,z$ and $P_{hT}$ dependent asymmetry is
generally consistent with that of $K^{-}$ meson. In addition, the Sivers
asymmetry in $\Lambda$ hyperon produced process is larger than that in charged
Kaon $K^{\pm}$ produced process.

Therefore, the future higher precision EICs can provide a unique opportunity
to extract the proton Sivers function of valence quark and sea quark from the
charged Kaon $K^{+}$ production, to analyze the flavor dependence of the
Sivers distribution function among sea quarks from $K^{-}$ produced SIDIS
process, to investigate the flavor dependence from $\Lambda$ hyperon produced
SIDIS process. Combining the experimental data from Drell-Yan process, it also
provides an ideal tool to study the sign change of the Sivers function and the
generalized universality of the T-odd distribution functions.

##  IV CONCLUSION

In this work, we apply the TMD factorization formalism to study the Sivers
asymmetry with $\sin\left(\phi_{h}-\phi_{s}\right)$ modulation in charged
$K^{\pm}$ produced and $\Lambda$ hyperon produced in SIDIS process at the
kinematical configurations of EIC and EicC. We take into account the TMD
evolution effects of distribution functions as well as fragmentation
functions. In our calculation we adopt the EIKV parametrization for the non-
perturbative Sudakov-like form factor, the accuracy of the perturbative
Sudakov-like form factor as well as the hard coefficients is kept at the NLL
order. Two different ways to deal with the energy dependence of Qiu-Sterman
function associated with the Sivers function are applied. The first approach
is to assume the Qiu-Sterman function evolves as the same way as the
unpolarized distribution functions $f_{1}(x,Q^{2})$. The second one is to
evolve Qiu-Sterman function considering an approximate evolution kernel for
the Qiu-Sterman function containing the homogenous terms by customising the
DGLAP kernel. The Sivers asymmetries are calculated as the functions of $x,z$,
and $P_{hT}$. Our numerical results demonstrate that the Sivers asymmetries of
charged Kaon production and $\Lambda$ hyperon production are measurable at the
kinematics of EIC and EicC, with the magnitudes of around several percents.
The results show that the Sivers asymmetry in $K^{\pm}$ and $\Lambda$ hyperon
production in SIDIS process can serve as an ideal tool to extract the
information of sea quark Sivers function as well as to constrain the flavor
dependence of Sivers function by utilizing future high energy and high
luminosity EICs. In addition, the difference between the Sivers asymmetries
from considering two different evolution kernels suggests that the DGLAP
evolution of the Qiu-Sterman function in the TMD evolution schemes will play a
role in the phenomenological calculation, which should be considered in the
future interpretation of experimental data as well as theoretical studies.

##  V ACKNOWLEDGMENTS

This work is partially supported by the NSFC (China) grants 11905187,11847217
and 12150013. X. Wang is supported by the China Postdoctoral Science
Foundation under Grant No. 2018M640680.

## References

  * Ashman et al. (1988) J. Ashman et al. (European Muon), Phys. Lett. B 206, 364 (1988).
  * Ashman et al. (1989) J. Ashman et al. (European Muon), Nucl. Phys. B 328, 1 (1989).
  * Sivers (1990) D. W. Sivers, Phys. Rev. D 41, 83 (1990).
  * Sivers (1991) D. W. Sivers, Phys. Rev. D 43, 261 (1991).
  * Collins (1993) J. C. Collins, Nucl. Phys. B 396, 161 (1993).
  * Brodsky et al. (2002a) S. J. Brodsky, D. S. Hwang, and I. Schmidt, Phys. Lett. B 530, 99 (2002a).
  * Brodsky et al. (2002b) S. J. Brodsky, D. S. Hwang, and I. Schmidt, Nucl. Phys. B 642, 344 (2002b).
  * Boer et al. (2003) D. Boer, S. J. Brodsky, and D. S. Hwang, Phys. Rev. D 67, 054003 (2003).
  * Collins (2002) J. C. Collins, Phys. Lett. B 536, 43 (2002).
  * Airapetian et al. (2005) A. Airapetian et al. (HERMES), Phys. Rev. Lett. 94, 012002 (2005).
  * Airapetian et al. (2009) A. Airapetian et al. (HERMES), Phys. Rev. Lett. 103, 152002 (2009).
  * Airapetian et al. (2020) A. Airapetian et al. (HERMES), J. High Energy Phys. 12, 010 (2020).
  * Alexakhin et al. (2005) V. Y. Alexakhin et al. (COMPASS), Phys. Rev. Lett. 94, 202002 (2005).
  * Ageev et al. (2007) E. S. Ageev et al. (COMPASS), Nucl. Phys. B 765, 31 (2007).
  * Alekseev et al. (2009) M. Alekseev et al. (COMPASS), Phys. Lett. B 673, 127 (2009).
  * Alekseev et al. (2010) M. G. Alekseev et al. (COMPASS), Phys. Lett. B 692, 240 (2010).
  * Adolph et al. (2012) C. Adolph et al. (COMPASS), Phys. Lett. B 717, 383 (2012).
  * Adolph et al. (2017) C. Adolph et al. (COMPASS), Phys. Lett. B 770, 138 (2017).
  * Alexeev et al. (2019) M. G. Alexeev et al. (COMPASS), Nucl. Phys. B 940, 34 (2019).
  * Qian et al. (2011) X. Qian et al. (Jefferson Lab Hall A), Phys. Rev. Lett. 107, 072003 (2011).
  * Zhao et al. (2014) Y. X. Zhao et al. (Jefferson Lab Hall A), Phys. Rev. C 90, 055201 (2014).
  * Aghasyan et al. (2017) M. Aghasyan et al. (COMPASS), Phys. Rev. Lett. 119, 112002 (2017).
  * Adamczyk et al. (2016) L. Adamczyk et al. (STAR), Phys. Rev. Lett. 116, 132301 (2016).
  * Anselmino et al. (2005) M. Anselmino, M. Boglione, U. D’Alesio, A. Kotzinian, F. Murgia, and A. Prokudin, Phys. Rev. D 72, 094007 (2005), [Erratum: Phys.Rev.D 72, 099903 (2005)].
  * Efremov et al. (2005) A. V. Efremov, K. Goeke, S. Menzel, A. Metz, and P. Schweitzer, Phys. Lett. B 612, 233 (2005).
  * Collins et al. (2006) J. C. Collins, A. V. Efremov, K. Goeke, S. Menzel, A. Metz, and P. Schweitzer, Phys. Rev. D 73, 014021 (2006).
  * Vogelsang and Yuan (2005) W. Vogelsang and F. Yuan, Phys. Rev. D 72, 054028 (2005).
  * Anselmino et al. (2009) M. Anselmino, M. Boglione, U. D’Alesio, A. Kotzinian, S. Melis, F. Murgia, A. Prokudin, and C. Turk, Eur. Phys. J. A 39, 89 (2009).
  * Anselmino et al. (2012) M. Anselmino, M. Boglione, and S. Melis, Phys. Rev. D 86, 014028 (2012).
  * Bacchetta and Radici (2011) A. Bacchetta and M. Radici, Phys. Rev. Lett. 107, 212001 (2011).
  * Echevarria et al. (2014a) M. G. Echevarria, A. Idilbi, Z.-B. Kang, and I. Vitev, Phys. Rev. D 89, 074013 (2014a).
  * Anselmino et al. (2017) M. Anselmino, M. Boglione, U. D’Alesio, F. Murgia, and A. Prokudin, J. High Energy Phys. 04, 046 (2017).
  * Martin et al. (2017) A. Martin, F. Bradamante, and V. Barone, Phys. Rev. D 95, 094024 (2017).
  * Boglione et al. (2018) M. Boglione, U. D’Alesio, C. Flore, and J. O. Gonzalez-Hernandez, J. High Energy Phys. 07, 148 (2018).
  * Bury et al. (2021) M. Bury, A. Prokudin, and A. Vladimirov, J. High Energy Phys. 05, 151 (2021).
  * Bacchetta et al. (2004a) A. Bacchetta, A. Schaefer, and J.-J. Yang, Phys. Lett. B 578, 109 (2004a).
  * Lu and Ma (2004) Z. Lu and B.-Q. Ma, Phys. Rev. D 70, 094044 (2004).
  * Pasquini and Yuan (2010) B. Pasquini and F. Yuan, Phys. Rev. D 81, 114013 (2010).
  * Maji et al. (2018) T. Maji, D. Chakrabarti, and A. Mukherjee, Phys. Rev. D 97, 014016 (2018).
  * Maji et al. (2017) T. Maji, D. Chakrabarti, and O. V. Teryaev, Phys. Rev. D 96, 114023 (2017).
  * Courtoy et al. (2008) A. Courtoy, F. Fratini, S. Scopetta, and V. Vento, Phys. Rev. D 78, 034002 (2008).
  * Yuan (2003) F. Yuan, Phys. Lett. B 575, 45 (2003).
  * Courtoy et al. (2009) A. Courtoy, S. Scopetta, and V. Vento, Phys. Rev. D 79, 074001 (2009).
  * Lyubovitskij et al. (2022) V. E. Lyubovitskij, I. Schmidt, and S. J. Brodsky, Phys. Rev. D 105, 114032 (2022).
  * Dong et al. (2019) H. Dong, D.-X. Zheng, and J. Zhou, Phys. Lett. B 788, 401 (2019).
  * He and Wang (2019) F. He and P. Wang, Phys. Rev. D 100, 074032 (2019).
  * Luan and Lu (2022) X. Luan and Z. Lu, Phys. Lett. B 833, 137299 (2022).
  * Accardi et al. (2016) A. Accardi et al., Eur. Phys. J. A 52, 268 (2016).
  * Abdul Khalek et al. (2022) R. Abdul Khalek et al., Nucl. Phys. A 1026, 122447 (2022).
  * Zeng et al. (2022) C. Zeng, T. Liu, P. Sun, and Y. Zhao, Phys. Rev. D 106, 094039 (2022).
  * Collins and Soper (1981) J. C. Collins and D. E. Soper, Nucl. Phys. B 193, 381 (1981), [Erratum: Nucl.Phys.B 213, 545 (1983)].
  * Collins et al. (1985) J. C. Collins, D. E. Soper, and G. F. Sterman, Nucl. Phys. B 250, 199 (1985).
  * Collins (2013) J. Collins, _Foundations of perturbative QCD_ , vol. 32 (Cambridge University Press, 2013), ISBN 978-1-107-64525-7, 978-1-107-64525-7, 978-0-521-85533-4, 978-1-139-09782-6.
  * Ji et al. (2004) X.-D. Ji, J.-P. Ma, and F. Yuan, Phys. Lett. B 597, 299 (2004).
  * Ji et al. (2005) X.-D. Ji, J.-P. Ma, and F. Yuan, Phys. Rev. D 71, 034005 (2005).
  * Aybat and Rogers (2011) S. M. Aybat and T. C. Rogers, Phys. Rev. D 83, 114042 (2011).
  * Collins and Rogers (2013) J. C. Collins and T. C. Rogers, Phys. Rev. D 87, 034018 (2013).
  * Echevarria et al. (2013) M. G. Echevarria, A. Idilbi, A. Schäfer, and I. Scimemi, Eur. Phys. J. C 73, 2636 (2013).
  * Pitonyak et al. (2014) D. Pitonyak, M. Schlegel, and A. Metz, Phys. Rev. D 89, 054032 (2014).
  * Boer (2009) D. Boer, Nucl. Phys. B 806, 23 (2009).
  * Arnold et al. (2009) S. Arnold, A. Metz, and M. Schlegel, Phys. Rev. D 79, 034005 (2009).
  * Lambertsen and Vogelsang (2016) M. Lambertsen and W. Vogelsang, Phys. Rev. D 93, 114013 (2016).
  * Collins and Hautmann (2000) J. C. Collins and F. Hautmann, Phys. Lett. B 472, 129 (2000).
  * Bacchetta et al. (2004b) A. Bacchetta, U. D’Alesio, M. Diehl, and C. A. Miller, Phys. Rev. D 70, 117504 (2004b).
  * Bacchetta et al. (2007) A. Bacchetta, M. Diehl, K. Goeke, A. Metz, P. J. Mulders, and M. Schlegel, JHEP 02, 093 (2007).
  * Idilbi et al. (2004) A. Idilbi, X.-D. Ji, J.-P. Ma, and F. Yuan, Phys. Rev. D 70, 074021 (2004).
  * Bacchetta and Prokudin (2013) A. Bacchetta and A. Prokudin, Nucl. Phys. B 875, 536 (2013).
  * Kang et al. (2011) Z.-B. Kang, B.-W. Xiao, and F. Yuan, Phys. Rev. Lett. 107, 152002 (2011).
  * Aybat et al. (2012) S. M. Aybat, J. C. Collins, J.-W. Qiu, and T. C. Rogers, Phys. Rev. D 85, 034043 (2012).
  * Echevarria et al. (2014b) M. G. Echevarria, A. Idilbi, and I. Scimemi, Phys. Rev. D 90, 014003 (2014b).
  * Landry et al. (2003) F. Landry, R. Brock, P. M. Nadolsky, and C. P. Yuan, Phys. Rev. D 67, 073016 (2003).
  * Qiu and Zhang (2001) J.-W. Qiu and X.-F. Zhang, Phys. Rev. Lett. 86, 2724 (2001).
  * Konychev and Nadolsky (2006) A. V. Konychev and P. M. Nadolsky, Phys. Lett. B 633, 710 (2006).
  * Davies et al. (1984) C. T. H. Davies, B. R. Webber, and W. J. Stirling, 1, I.95 (1984).
  * Ellis et al. (1997) R. K. Ellis, D. A. Ross, and S. Veseli, Nucl. Phys. B 503, 309 (1997).
  * Kang et al. (2016) Z.-B. Kang, A. Prokudin, P. Sun, and F. Yuan, Phys. Rev. D 93, 014009 (2016).
  * Sun and Yuan (2013) P. Sun and F. Yuan, Phys. Rev. D 88, 114012 (2013).
  * Martin et al. (2009) A. D. Martin, W. J. Stirling, R. S. Thorne, and G. Watt, Eur. Phys. J. C 63, 189 (2009).
  * de Florian et al. (2007) D. de Florian, R. Sassot, and M. Stratmann, Phys. Rev. D 75, 114010 (2007).
  * Yang et al. (2017) Y. Yang, Z. Lu, and I. Schmidt, Phys. Rev. D 96, 034010 (2017).
  * Botje (2011) M. Botje, Comput. Phys. Commun. 182, 490 (2011).
  * Kang and Qiu (2009) Z.-B. Kang and J.-W. Qiu, Phys. Rev. D 79, 016003 (2009).


* van Dishoeck & Black (1988) van Dishoeck, E. F. & Black, J. H. 1988, ApJ, 334, 771
  * Watanabe et al. (2017) Watanabe, Y., Nishimura, Y., Harada, N., et al. 2017, ApJ, 845, 116
  * Wilson et al. (2005) Wilson, B. A., Dame, T. M., Masheder, M. R. W., & Thaddeus, P. 2005, A&A, 430, 523
  * Wilson & Rood (1994) Wilson, T. L. & Rood, R. 1994, ARA&A, 32, 191
  * Wolfire et al. (2010) Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191
  * Yun et al. (2021) Yun, H.-S., Lee, J.-E., Choi, Y., et al. 2021, ApJS, 256, 16
  * Zari et al. (2016) Zari, E., Lombardi, M., Alves, J., Lada, C. J., & Bouy, H. 2016, A&A, 587, A106
  * Zhang et al. (2018) Zhang, G.-Y., Xu, J.-L., Vasyunin, A. I., et al. 2018, A&A, 620, A163
  * Zucker et al. (2019) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2019, ApJ, 879, 125

##  Appendix A Maps of sample points for California and Orion A

Figures 12 and 13 represent in color scale the column density maps of
California and Orion A derived by Lada et al. (2017) and Lombardi et al.
(2014). The numbers mark the positions selected using the stratified random
sampling technique, and numerical values indicate the column density bin to
which the position belongs.

Figure 12: Sampling of the California cloud. The color images represent the
distribution of column density as derived by Lada et al. (2017) from Herschel
Space Observatory data, and the numbers indicate the positions chosen using
stratified random sampling. The right panel is an expanded view of the L1482
region, which contains some of the highest column density positions. The
numbers refer to the column density bin of the position (1 is lowest and 8 is
highest), the 10 pc scale bar assumes a distance of 470 pc (Zucker et al.
2019), and the wedge scale is in units of cm-2. Figure 13: Sampling of the
Orion A cloud. The color images represent the distribution of column density
as derived by Lombardi et al. (2014) from Herschel Space Observatory data, and
the numbers indicate the positions chosen using stratified random sampling.
The right panel is an expanded view of the ISF, which contains some of the
highest column density positions. The numbers refer to the column density bin
of the position (1 is lowest and 12 is highest), the 10 pc scale bar assumes a
distance of 432 pc (Zucker et al. 2019), and the wedge scale is in units of
cm-2.

##  Appendix B Sample positions and 3 mm line intensities

Tables 8 and 9 present the coordinates, H2 column density estimates, and 3 mm
line intensities for all the target positions in the California and Orion
clouds. Full versions of the tables are available at the CDS.

Table 8: Sample positions and line intensities for California.

Positionaa$a$The first number in the position name indicates the column density bin and the second one indicates the order in our sampling sequence. Numbers in parentheses indicate 1$\sigma$ errors. | RA(J2000) | Dec(J2000) |  $N$(H2) |  12CO(1–0) |  13CO(1–0) | C18O(1–0) | HCN(1–0) | CS(2–1)
---|---|---|---|---|---|---|---|---
| (h m s) | (° ′ ″) | (cm-2) | (K km s-1) | (K km s-1) | (K km s-1) | (K km s-1) | (K km s-1)
CAL-08_10 | 04 19 28.8 | +37 59 04 | $3.81(0.04)\times 10^{22}$ | 21 (2) | 4.6 (0.5) | 2.1 (0.2) | 0.69 (0.07) | 0.65 (0.07)
CAL-08_09 | 04 25 36.1 | +37 06 05 | $4.46(0.05)\times 10^{22}$ | 19 (2) | 6.6 (0.7) | 2.3 (0.2) | 1.1 (0.1) | 0.93 (0.09)
CAL-08_08 | 04 25 40.2 | +37 07 01 | $3.83(0.03)\times 10^{22}$ | 18 (2) | 5.6 (0.6) | 2.4 (0.2) | 0.72 (0.07) | 0.69 (0.07)
CAL-08_07 | 04 30 48.5 | +34 58 03 | $4.31(0.14)\times 10^{22}$ | 30 (3) | 9.9 (1.0) | 1.3 (0.1) | 4.8 (0.5) | 2.8 (0.3)
CAL-08_06 | 04 30 39.4 | +35 29 02 | $4.08(0.05)\times 10^{22}$ | 43 (4) | 17 (2) | 4.3 (0.4) | 1.6 (0.2) | 1.4 (0.1)
CAL-08_05 | 04 25 32.0 | +37 07 01 | $3.92(0.02)\times 10^{22}$ | 14 (1) | 4.3 (0.4) | 1.5 (0.2) | 1.0 (0.1) | 0.50 (0.05)
CAL-08_04 | 04 30 36.8 | +35 54 03 | $4.88(0.12)\times 10^{22}$ | 56 (6) | 17 (2) | 3.6 (0.4) | 7.2 (0.7) | 6.0 (0.6)
CAL-08_03 | 04 21 41.0 | +37 33 05 | $4.07(0.07)\times 10^{22}$ | 31 (3) | 8.5 (0.9) | 2.5 (0.3) | 2.6 (0.3) | 2.6 (0.3)
CAL-08_02 | 04 30 40.4 | +35 29 01 | $4.37(0.04)\times 10^{22}$ | 50 (5) | 19 (2) | 4.7 (0.5) | 1.4 (0.1) | 1.4 (0.1)
CAL-08_01 | 04 21 42.0 | +37 33 04 | $3.84(0.07)\times 10^{22}$ | 30 (3) | 7.1 (0.7) | 2.3 (0.2) | 1.8 (0.2) | 2.3 (0.2)

101010A full version of this table is available at the CDS.

Table 9: Sample positions and line intensities for Orion A.

Positionaa$a$The first number in the position name indicates the column density bin and the second one indicates the order in our sampling sequence. Numbers in parentheses indicate 1$\sigma$ errors. | RA(J2000) | Dec(J2000) |  $N$(H2) |  12CO(1–0) |  13CO(1–0) | C18O(1–0) | HCN(1–0) | CS(2–1)
---|---|---|---|---|---|---|---|---
| (h m s) | (° ′ ″) | (cm-2) | (K km s-1) | (K km s-1) | (K km s-1) | (K km s-1) | (K km s-1)
ORIA-12_10 | 05 35 14.9 | -06 37 50 | $1.8(0.1)\times 10^{23}$ | 1380 (140) | 124 (12) | 14.8 (2) | 611 (62) | 110 (11)
ORIA-12_09 | 05 35 12.7 | -06 37 50 | $1.8(0.1)\times 10^{23}$ | 916 (92) | 73 (7) | 7.4E (0.7) | 285 (29) | 47 (5)
ORIA-12_08 | 05 35 14.5 | -06 39 10 | $2.3(0.4)\times 10^{23}$ | 333 (33) | 56 (6) | 5.2 (0.5) | 68 (7) | 30 (3)
ORIA-12_07 | 05 35 14.5 | -06 38 03 | $1.8(0.1)\times 10^{23}$ | 1000 (100) | 84 (8) | 8.5 (8) | 365 (37) | 59 (6)
ORIA-12_06 | 05 35 10.9 | -06 37 03 | $1.8(0.4)\times 10^{23}$ | 428 (43) | 52 (5) | 4.0 (0.4) | 63 (6) | 25 (2)
ORIA-12_05 | 05 35 14.0 | -06 37 43 | $1.8(0.1)\times 10^{23}$ | 258 (26) | 172 (17) | 14 (1) | 1050 (110) | 220 (22)
ORIA-12_04 | 05 35 15.4 | -06 37 36 | $1.8(0.1)\times 10^{23}$ | 2100 (220) | 136E (14) | 15 (2) | 816 (82) | 133 (13)
ORIA-12_03 | 05 35 18.5 | -06 40 31 | $2.6(0.6)\times 10^{23}$ | 337 (34) | 57 (6) | 4.8 (0.5) | 58 (6) | 21 (2)
ORIA-12_02 | 05 35 12.3 | -06 37 30 | $1.8(0.1)\times 10^{23}$ | 642 (64) | 80 (8) | 7.7 (0.8) | 253 (26) | 53 (5)
ORIA-12_01 | 05 35 23.9 | -06 58 30 | $1.9(0.3)\times 10^{23}$ | 98 (10) | 31 (3) | 5.6 (0.6) | 25 (3) | 11 (1)

111111A full version of this table is available at the CDS.

##  Appendix C Comparison with the Orion A data from Yun et al. (2021)

Figure 14: Comparison between the Orion A sampling observations and the
mapping data of Yun et al. (2021). The figure presents a scanned copy of Fig.
15 from Yun et al. (2021), with red circles representing the intensities
derived from the sampling observations. The solid blue line represents the
approximate $3\sigma$ limit of the mapping data calculated using the velocity
range of integration given by Yun et al. (2021) in their Table 2 and the rms
noise in the spectra given in their Fig. 8. The dashed line represents the
equivalent limit for the sampling data. See the main text for further details.

Paper I presented a comparison between the Perseus 12CO(1–0) and 13CO(1–0)
sampling data and the mapping observations of the same lines by Ridge et al.
(2006) using the FCRAO 14 m telescope. This comparison showed that sampling
observations can reproduce at the 25% level both the mean intensity and the
intensity dispersion as a function of $N$(H2) derived from the mapping
observations (Paper I, Sect. 4.1). For the Orion A cloud, Yun et al. (2021)
have recently presented maps of several 3 mm emission lines carried out with
the twin TRAO 14 m telescope, and this new data set provides an excellent
opportunity to test again the ability of the sampling method to reproduce the
global emission of a molecular cloud. The observations of Yun et al. (2021)
consist of maps in the lines of 13CO(1–0), C18O(1–0), HCN(1–0), CS(2–1),
HCO+(1–0), and N2H+(1–0) that cover almost the full extent of the Orion A
cloud emission (see the coverage in their Fig. 1). In contrast with the data
from Ridge et al. (2006), the Yun et al. (2021) data are not publicly
available, so our comparison with the Orion A data has been made by
superposing our intensity distributions over those presented by Yun et al.
(2021) in their Fig. 15. Since Yun et al. (2021) presented their intensities
in $T_{\mathrm{A}}^{*}$ units (Neal Evans, private communication), which
correspond to about half the intensities in units of $T_{\mathrm{MB}}$ (see
the beam efficiencies in their Table 1), to correctly simulate a line
comparison in $T_{\mathrm{MB}}$ units (as done with the data from Ridge et al.
2006), we divided our main-beam intensities by a factor of 2. The results from
the comparison are shown in Fig. 14, which presents a superposition of our
scaled-down intensities (red circles) over the intensities presented by Yun et
al. (2021) in their Fig. 15.

As can be seen in Fig. 14, our sampling data match well the distribution of
intensity as a function of $N$(H2) derived from the mapping data for all
observed lines. The best agreement is seen in the 13CO(1–0) data since this
line presents the highest signal to noise, and its distribution can be
reliably followed over two orders of magnitude in $N$(H2). The intensity of
the other lines cannot be followed over such a large range of column
densities, but in all cases, the sampling data overlap well with the
distribution of mapping intensities and presents a similar level of dispersion
over the range of column densities for which the mapping data can be
considered complete (over the blue solid line). Given the results in Fig. 14,
we conclude that the sampling observations provide a good approximation to the
distribution of intensities with $N$(H2) found by mapping, and that they do
not seem to miss any particular regime of intensities that is present in the
maps.

A more quantitative comparison between our sampling data and the mapping
observations of Yun et al. (2021) could be carried out using line
luminosities. As discussed in our Sect. 3.7, line luminosities can be easily
estimated from the sampling data by multiplying the observed line intensities
by the surface area of the different column density bins, and luminosity
estimates for the main transitions observed in our three target clouds are
presented in Table. 6. Yun et al. (2021) have also calculated luminosities for
the lines observed in their survey, but unfortunately the values given by
these authors in their Sect. 5.4 suffer from several inaccuracies that make it
necessary a revision (Neal Evans, private communication). A preliminary
comparison with the revised luminosities of 13CO(1–0) and C18O(1–0) suggests a
level of agreement with our sampling results at better than the 50% level,
which is similar to the agreement found in Sect. 3.7 when comparing our
luminosities with other estimates in the literature. A definitive comparison
between our values and those of Yun et al. (2021) awaits the final estimate of
the revised luminosities by these authors.

##  Appendix D Temperature effects

###  D.1 Gas temperature estimates

Figure 15: C18O(2–1)/C18O(1–0) intensity ratio as a function of the gas
temperature determined using three different methods: peak of the 12CO(1–0)
line (top), HCN(1–0)/HNC(1–0) line ratio method of Hacar et al. (2020)
(middle), and the empirical fit based on the dust temperature described in the
text (bottom). The data are color-coded as in Fig. 1, and the purple points in
the bottom-right panel indicate the Orion A positions where the NH3-based
temperature estimate of Friesen et al. (2017) has been used. The dashed lines
represent the expected line ratio under optically thin LTE conditions, and the
ability of the data to match this curve reflects the quality of the
temperature determination.  Figure 16: Distribution of dust temperature as a
function of H2 column density for the sampled positions of California,
Perseus, and the regions of Orion A for which no NH3-based gas temperature is
available. The gray line represents the simple analytic fit to the equilibrium
dust temperature described in the text. Dust temperature data are from Lada et
al. (2017, California), Zari et al. (2016, Perseus), and Lombardi et al.
(2014, Orion A).

To investigate the effect of temperature in the emission of the different
lines, we first needed to determine its mean value along each observed line of
sight. The ideal tool for this is a tracer that is bright enough to be
detected at all column densities, has an excitation close to the gas kinetic
temperature, and is optically thin so it samples the entire column density of
gas at each position. The tracer in our data set that best approaches these
characteristics is C18O, and in Paper I we used the C18O(2–1)/C18O(1–0) ratio
to determine that the temperature in the Perseus positions is approximately
constant at about 11 K, in agreement with previous ammonia results (Rosolowsky
et al. 2008). While a similar analysis could be done with the Orion A data,
our C18O(2–1) observations of California were limited to 13 positions (out of
80) due to observing-time constraints. As a result, it is not possible to
derive a C18O line ratio over most of this cloud, and an alternative
temperature determination is needed. In this appendix we explore three
possible methods and judge their merit testing them against the C18O line
ratio, which we consider the best available temperature indicator.

Our first (and simplest) temperature estimator is the peak intensity of the
12CO(1–0) line. 12CO is easily thermalized, and the $J$=1–0 line is optically
thick. Under these conditions, the gas kinetic temperature $T_{\mathrm{k}}$
can be assumed to equal the excitation temperature, which is given by

$T_{\mathrm{ex}}=\frac{h\nu/k}{\ln[1+\frac{h\nu/k}{T_{\mathrm{B}}-J_{\nu}(T_{\mathrm{bg}})}]},$
(8)

where $\nu$ is the line frequency, $J_{\nu}$ the radiation temperature, and
$T_{\mathrm{b}g}$ the cosmic background temperature (=2.73 K). To determine
the 12CO(1–0) peak brightness temperature at each position, we fitted each
spectrum using a Gaussian profile. This type of profile does not always
represent the best fit to the complex 12CO(1–0) spectrum, but a visual
inspection of the fits shows that the method provides a good approximation to
the line peak. Using this estimate and Eq. D.1, we determined a
12CO(1–0)-based kinetic temperature at each cloud position in out sample. The
top panel of Fig. 15 presents the C18O(2–1)/C18O(1–0) line ratio derived from
our available data as a function of the gas temperature estimated using the
12CO(1–0) peak temperature. As can be seen, the line ratios for Perseus and
Orion A systematically lie to the right of the dashed line, which corresponds
to the prediction for the optically thin local thermodynamic equilibrium (LTE)
conditions expected to apply. This systematic shift indicates that the
12CO(1–0) peak method overestimates the value of the C18O-inferred gas kinetic
temperature by a significant margin. A similar overestimate of the gas
temperature by the 12CO data was first noticed by Castets et al. (1990) in
their multiline CO study of Orion A, and recent discussions of this effect
have been presented by Hacar et al. (2020) and Roueff et al. (2021). The most
likely cause of this overestimate is the high optical depth of the 12CO(1–0)
line, which makes the emission originate from the cloud outer layers. These
layers are exposed to the external UV radiation field and therefore warmer
than the attenuated interior (e.g., Castets et al. 1990). A comparison of the
CO-derived temperatures in Orion A with the NH3-based temperature estimates of
Friesen et al. (2017) supports this interpretation, and suggests that the 12CO
data can overestimate the internal gas temperature by up to a factor of 2.

An alternative method for estimating the gas kinetic temperature is to use the
HCN(1–0)/HNC(1–0) line ratio, which behaves as a chemical thermometer due to
the high sensitivity to temperature of the relative abundance of the two HCN
isomers (Hacar et al. 2020). Using NH3 data from the Orion A cloud, Hacar et
al. (2020) have calibrated the relation between gas temperature and HCN/HNC
ratio, and proposed a simple prescription to derive the gas kinetic
temperature (their Eqs. 3 and 4). We applied this prescription to our survey
data, and determined a gas temperature at each observed position. In the
middle row of panels of Fig. 15 we present the C18O(2–1)/C18O(1–0) line ratio
as a function of the gas temperature estimated using the HCN/HNC line method.
As can be seen, the Orion A data follow more closely the optically thin LTE
limit line than when using the 12CO(1–0) line method, although this may simply
reflect the fact that the HCN/HNC ratio was calibrated using NH3 data from
this cloud. In contrast, the HCN/HNC ratio method does not align well the
Perseus line ratios with the analytic curve, indicating that the method does
not provide a reliable estimate of the gas kinetic temperature in this cloud.
This should have been expected since the HCN/HNC ratio method loses accuracy
for temperatures below 15 K (Hacar et al. 2020).

As a third alternative to estimate gas temperatures, we explored the use of
the dust temperature as a reference. Dust temperatures are derived together
with H2 column densities in the analysis of the Herschel data (Lombardi et al.
2014; Zari et al. 2016; Lada et al. 2017), so they are available at each cloud
position. While the gas and dust components are only thermally coupled at
densities in excess of $\approx 10^{4}$ cm-3 (Goldsmith 2001), the dust
temperature can still serve as an indicator of the presence (or absence) of
energetic perturbations that may affect the gas temperature. In our clouds,
the largest temperature perturbation occurs in Orion A, where the gas and dust
in the vicinity of the ONC are heated by the newly formed stars. Since this
temperature perturbation can be accurately characterized using the NH3-derived
temperatures of (Friesen et al. 2017), we used the NH3 temperatures for the
ONC vicinity, and used the dust temperature method only in the remaining parts
of the clouds. In these parts, the temperature excursions are likely very
small, as suggested by the dust temperature profiles shown in Fig. 16. These
profiles indicate that the distribution of dust temperature in the three
clouds has little dispersion, and that most values cluster along a well-
defined band that runs almost horizontally. This band likely traces the
equilibrium temperature of the dust for unperturbed cloud conditions since is
very similar in the three clouds despite their significant differences in
star-formation activity. It presents a slight decrease as $N$(H2) increases
that likely results from the inward attenuation of the interstellar radiation
field, which is the main heating agent of the dust (Evans et al. 2001; Planck
Collaboration et al. 2011). Following work by Hacar et al. (2017b) and Hocuk
et al. (2017), we fitted the trend in the three clouds using the expression
$T_{\mathrm{dust}}=19.0-4.0\log_{10}(N(\mathrm{H}_{2})/10^{21}\,\mathrm{cm}^{-2})$
(gray line in the figure).

To use the dust temperature as a proxy for the gas temperature, we followed a
simple empirical approach. We first subtracted from the dust temperature the
observed gradient with $N$(H2) since this gradient is caused by the
attenuation of the interstellar radiation field, and is not expected to affect
the gas component because is heated by cosmic rays and not thermally coupled
to the dust in the outer cloud (Goldsmith 2001; Galli et al. 2002). To the
result, we simply added a global offset of 8 K, which has been derived from an
empirical match to the C18O(2–1)/C18O(1–0) line ratio. As seen in the bottom
row of panels in Fig. 15, the empirical dust fit, when complemented with the
NH3-derived temperatures of Friesen et al. (2017), provides a better match to
the C18O(2–1)/C18O(1–0) line ratio than the fits obtained using either the
CO(1–0) peak intensity or the HCN/HNC line ratio. For this reason, we used
this fit as the temperature indicator in the three clouds of our sample.
Further work is needed to understand the correct level of coupling between the
gas and dust temperatures through the clouds.

###  D.2 Temperature-correction factors

To determine (and possibly correct for) the effect of temperature variations
in the emission profile of the different clouds, we needed to determine how
the intensity of any line depends on the gas kinetic temperature. This
determination requires a certain degree of approximation because the intensity
of a line depends on a complex combination of the physical and chemical
conditions of the emitting cloud, which are in principle not well known. The
observed systematic behavior of the line intensities, however, suggests that
the main characteristics of the emission are dominated by relatively few cloud
properties that could be approximated using a simplified model. If this is the
case, we could define a series of correction factors that relate the intensity
emitted by gas at an arbitrary temperature $T_{\mathrm{k}}$ with the emission
that the same gas would emit if it were a reference temperature of 10 K. These
factors are defined as

$f_{10\mathrm{K}}(T_{\mathrm{k}})\equiv\frac{I\;(T_{\mathrm{k}})}{I\;(10\mathrm{\leavevmode\nobreak\
K})},$ (9)

were $I(T_{\mathrm{k}})$ and $I(10\leavevmode\nobreak\ K)$ are the intensities
at the two temperatures. To first approximation, we could use these factors to
convert any observed line intensity into the expected intensity that the same
line of sight would emit if it were at 10 K.

As a first step to understand the behavior of the correction factors, we
considered two simple limiting cases. If the line is optically thick and
thermalized ($T_{\mathrm{e}x}=T_{\mathrm{k}}$), which may be the case for a
low $J$ 12CO transition, the equation of radiative transfer states that the
integrated intensity will be equal to
$(J_{\nu}(T_{\mathrm{k}})-J_{\nu}(T_{\mathrm{b}g}))\;\Delta\mathrm{V}$, where
$J_{\nu}$ is the radiation temperature, $T_{\mathrm{b}g}$ is the cosmic
background temperature, and $\Delta\mathrm{V}$ the linewidth in units of
velocity. Under these conditions, the correction factor defined as the line
intensity at an arbitrary kinetic temperature divided by the intensity at 10 K
will be

$f^{\mathrm{t}hick}_{10\mathrm{K}}(T_{\mathrm{k}})=\frac{J_{\nu}(T_{\mathrm{k}})-J_{\nu}(T_{\mathrm{b}g})}{J_{\nu}(10\leavevmode\nobreak\
\mathrm{K})-J_{\nu}(T_{\mathrm{b}g})}.$ (10)

Figure 17: Intensity predictions of the observed lines as a function of gas
temperature normalized to their value at 10 K. The intensities have been
calculated by running the model of the Perseus cloud presented in Paper I for
12 different gas temperatures. For each temperature, six different H2 column
densities in the range $2\times 10^{21}$ to $2\times 10^{23}$ cm-2 have been
assumed, and each circle represents the intensity average with an error bar
that indicates the 1$\sigma$ dispersion (and is often smaller than the marker
size). The upper block of panels contains the CO isotopologs, the middle block
contains the traditional dense-gas tracers and their isotopologs, and the
bottom panel contains N2H+(1–0). The solid red lines represent the empirical
fits described in the text and used in the analysis to correct for gas
temperature variations. The dashed black lines in the CO and C18O panels
represent the optically thick and thin LTE solutions discussed in the text.

Another simple illustrative case is the optically thin line of a species in
LTE, in which the excitation temperature of all transitions equals the kinetic
temperature. For this case, which may be applicable to the C18O lines, the
integrated intensity will be
$[J_{\nu}(T_{\mathrm{k}})-J_{\nu}(T_{\mathrm{b}g})]\;\tau_{\nu}\;\Delta\mathrm{V}$,
where the optical depth $\tau_{\nu}$ is given by

$\tau_{\nu}=\frac{c^{3}}{{8\pi\nu^{3}}}\;\frac{g_{\mathrm{u}}A_{\mathrm{ul}}N}{{\Delta
V}}\;\frac{(1-e^{{-h\nu}/kT_{\mathrm{k}}})\;e^{-E_{\mathrm{l}}/kT_{\mathrm{k}}}}{Q(T_{\mathrm{k}})},$
(11)

and where $g_{\mathrm{u}}$ is the statistical weight of the upper level,
$A_{\mathrm{ul}}$ is the Einstein coefficient of spontaneous emission, $N$ the
molecular column density, and $Q$ the partition function calculated as a sum
over all relevant energy levels. As in the thick case, we can derive a
correction factor by dividing the expected line intensity at any given
temperature by the intensity at 10 K.

In real clouds, most molecular tracers are not thick and thermalized or thin
and in LTE, so the conversion factors between temperatures will likely deviate
from the previous two simple limits. To calculate these factors, we therefore
had to solve the equation of radiative transfer under physical and chemical
conditions that match those of real clouds. In Paper I, we presented a
radiative transfer model of the Perseus cloud that was able to reproduce the
intensity distribution of multiple lines (see Tables 3 and D.1 in Paper I ).
This model used a simple parameterization of the gas physical and chemical
conditions, and in particular, assumed a constant gas temperature of 11 K as
indicated by previous ammonia estimates (Rosolowsky et al. 2008).

To model how the intensity of the different lines is expected to vary with the
gas kinetic temperature, we reran the Perseus model for 12 different
temperatures that cover the range 8-100 K. For each gas temperature, we
estimated the emergent intensity of the different lines as a function of H2
column density using six logarithmically spaced H2 column densities that cover
the range $2\times 10^{21}$-$2\times 10^{23}$ cm-2. For each line, we
normalized the intensity at any given temperature by the value at a reference
temperature of 10 K as described in Equ. 9.

Figure 17 (top panels) shows the results from our radiative transfer model for
the observed transitions of the CO isotopologs. The circles represent the
intensity ratio averaged over the range of column densities covered by each
model, and the error bars represent the corresponding dispersion. The red
lines represent a set of analytic fits described below, and the black dashed
lines represent the two previously described limits: the optically thick
thermalized limit for 12CO and the optically thin LTE case for C18O. As can be
seen, the model CO intensities approximately follow the optically thick
prediction at low temperatures. As this parameter increases over about 50 K,
however, the higher $J$ levels of CO become more populated and the optical
depth of the lower transitions starts to drop, decreasing the match with the
thick approximation. As also shown in the figure, the model C18O intensities
follow the optically thin LTE limit at low temperatures, but as the
temperature increases, the LTE approximation breaks down and the intensities
deviate from the dashed line.

Model results for the intensity ratio of the traditional dense-gas tracers and
their isotopologs are presented in the middle group of panels of Fig. 17, and
the special case of N2H+(1–0) is shown in the bottom panel. As can be seen,
the optically thick lines from the main isotopologs present a higher
sensitivity to the temperature increase, but the general behavior of both
thick and thin lines is similar. It is characterized by an increase of the
line ratio at low temperatures and a gradual flattening as the temperature
exceeds 30-50 K.

We can summarize the results of our modeling by saying that the intensity of
all the observed lines is expected to vary smoothly and systematically with
the gas kinetic temperature in the range of interest (8-100 K). Our model
results provide a simple recipe to quantify and, if necessary, correct for
changes in the gas temperature across each cloud. Since our model results have
been calculated for a discrete number of temperatures, it may be necessary to
interpolate the corrections to the exact temperature value of each position.
This could be done by interpolating the model results, but it seems more
convenient to find analytic expressions that fit the model results and can be
used to quickly compute correction factors for any gas temperature.

Inspired by the similar behavior of all the different dense-gas tracers, and
by some of the factors that appear in the LTE equations, we explored the
possibility of fitting the dense-gas tracer data using the simple expression

$f_{10\mathrm{K}}(T_{\mathrm{k}})=\frac{\exp(-T_{0}/T_{\mathrm{k}})}{\exp(-T_{0}/10\mathrm{\leavevmode\nobreak\
K})},$ (12)

where $T_{0}$ is a reference temperature that has been determined by fitting
the data in Fig. 17. The fit results are represented with red lines in the
figure, and the $T_{0}$ values for each dense gas transition are given in
Table 10. As can be seen, a simple exponential factor provides a reasonable
fit to all dense-gas tracer lines and requires $T_{0}$ values in the narrow
range 7-13 K.

Fitting the intensity ratios of the CO isotopologs requires more complex
expressions because the model results present a diversity of behaviors that
range from the quasi-linear increase of the thick main isotopologs to the
systematic decrease of the C18O lines. Using again the LTE formulas as an
inspiration, and favoring the use of similar factors for similar lines, we
have found that it is possible to fit the model results using an expression of
the form

$f_{10\mathrm{K}}(T_{\mathrm{k}})=\frac{g_{1}(T_{\mathrm{k}})}{g_{1}(10\mathrm{\leavevmode\nobreak\
K})}\times\frac{g_{2}(T_{\mathrm{k}})}{g_{2}(10\mathrm{\leavevmode\nobreak\
K})},$ (13)

where $g_{1}(T_{\mathrm{k}})$ and $g_{2}(T_{\mathrm{k}})$ are the empirical
functions given in Table 11. As can be seen in Fig. 17, this factorization
provides a reasonable fit to intensity correction factors of all CO
isotopologs in the 8-100 K range. While convenient, these expressions only
represent empirical fits to the model results, and do not have any underlying
theoretical basis.

Table 10: Reference T0 temperature for dense-gas tracers.  Line |  $T_{0}$ [K] | Line |  $T_{0}$ [K]
---|---|---|---
HCN(1–0) | 13 | H13CN(1–0) | 12
CS(2–1) | 12 | C34S(2–1) | 9
HNC(1–0) | 10 | HN13C(1–0) | 7
HCO+(1–0) | 13 | H13CO+(1–0) | 8
N2H+(1–0) | 7 |  |
Table 11: Temperature correction factors for CO isotopologs.  Line | $g_{1}(T_{\mathrm{k}})$ | $g_{2}(T_{\mathrm{k}})$ |  | $g_{1}(T_{\mathrm{k}})$ | $g_{2}(T_{\mathrm{k}})$
---|---|---|---|---|---
CO(1-0) | $T_{\mathrm{k}}^{1.3}$ | $(1-e^{-50/T_{\mathrm{k}}})$ | CO(2-1) | $T_{\mathrm{k}}^{1.4}$ | $(1-e^{-80/T_{\mathrm{k}}})$
13CO(1-0) | $e^{-4/T_{\mathrm{k}}}$ | $(1-e^{-50/T_{\mathrm{k}}})$ |  13CO(2-1) | $e^{-14/T_{\mathrm{k}}}$ | $(1-e^{-80/T_{\mathrm{k}}})$
C18O(1-0) | 1 | $(1-e^{-50/T_{\mathrm{k}}})$ | C18O(2-1) | $e^{-10/T_{\mathrm{k}}}$ | $(1-e^{-80/T_{\mathrm{k}}})$

##  Appendix E Plots of isotopic ratios

Figs. 18 and 19 present intensity ratios as a function of H2 column density
for the rare CO isotopologs and the traditional dense-gas tracers. A
comparison with the standard isotopic value (represented with a dashed line)
indicates that the emission from the rare CO isotopologs is optically thin,
while the emission from the dense-gas tracers is optically thick.

Figure 18: C18O(1–0)/C17O(1–0) intensity ratio as a function of the H2 column
density used to determine the optical depth of the rare isotopologs. The
dashed horizontal line represents the standard isotopic ratio of 3.2 (Wilson &
Rood 1994), and the good match with the observations indicates that both
transitions are optically thin in all three clouds. Note the restricted range
of the x-axis due to the low S/N of the line ratio for very low values of
$N$(H2). Figure 19: Isotopic intensity ratios for the traditional dense-gas
tracers as a function of H2 column density for California, Perseus, and Orion
A. The $N$(H2) range has been reduced compared to other plots due to the low
signal-to-noise ratio of the data with low $N$(H2). The dashed horizontal
lines represent the standard interstellar medium isotopic ratios (77 for
12C/13C and 22 for 32S/34S), and the low value of the observed ratios
indicates that the emission of the main isotopolog is optically thick. Note
how the outflow points from Orion A (near $2\times 10^{23}$ cm-2) have lower
optical depths in some tracers but are still optically thick. No temperature
correction has been applied to the data.

##  Appendix F Luminosities from mapping and sampling methods

Table 12 presents the set of luminosity values used in the comparison between
sampling and mapping luminosities shown in Fig. 6. The luminosities from
mapping observations were chosen from the literature requiring that they were
well documented and reliable calibrated, or that the mapping data were
publicly available so that the luminosities could be independently estimated.
If necessary, luminosities have been re-scaled to the Gaia-based distances of
Zucker et al. (2019) given in Sect. 1. To compare with the Perseus estimates
from Ridge et al. (2006) and the Orion A estimates from Nishimura et al.
(2015), we restricted our sampling estimates to the same cloud area not masked
by these authors, for which we used their publicly available data. A full
discussion of the luminosity comparison is given in Sect. 3.7.

Table 12: Line luminosities from mapping and sampling measurements.

Line | $L_{\mathrm{sampling}}$ | $L_{\mathrm{mapping}}$ | Map ref.aa$a$Reference code for mapping luminosities: Dam23 = Dame & Lada (2023), Lew22 = Lewis et al. (2022), Nis15 = Nishimura et al. (2015), Rid06 = Ridge et al. (2006)
---|---|---|---
| (K km s-1 pc2) | (K km s-1 pc2) |
Orion A
12CO(1–0) | 28,200 | 21,300 | Lew22
12CO(2–1) | 19,100 | 16,000 | Nis15
13CO(2–1) | 3,350 | 2,030 | Nis15
C18O(2–1) | 91 | 53 | Nis15
Perseus
12CO(1–0) | 3,890 | 3,920 | Rid06
12CO(1–0) | 6,900 | 8,530 | Lew22
13CO(1–0) | 700 | 810 | Rid06
HCN(1-0) | 77 | 67 | Dam23
California
12CO(1–0) | 35,200 | 30,500 | Lew22
121212