# Rotating black strings in $f(R)$-Maxwell theory

A<EMAIL_ADDRESS>S. Salarpour 3 and Y. Bahrampour
<EMAIL_ADDRESS>1 Physics Department and Biruni Observatory, College of
Sciences, Shiraz University, Shiraz 71454, Iran
2 Research Institute for Astrophysics and Astronomy of Maragha (RIAAM), P.O.
Box 55134-441, Maragha, Iran
3 Department of Physics, Shahid Bahonar University, P.O. Box 76175, Kerman,
Iran
4 Department of Mathematics, Shahid Bahonar University, P.O. Box 76175,
Kerman, Iran

###### Abstract

In general, the field equations of $f(R)$ theory coupled to a matter field are
very complicated and hence it is not easy to find exact analytical solutions.
However, if one considers traceless energy-momentum tensor for the matter
source as well as constant scalar curvature, one can derive some exact
analytical solutions from $f(R)$ theory coupled to a matter field. In this
paper, by assuming constant curvature scalar, we construct a class of charged
rotating black string solutions in $f(R)$-Maxwell theory. We study the
physical properties and obtain the conserved quantities of the solutions. The
conserved and thermodynamic quantities computed here depend on function
$f^{\prime}(R_{0})$ and differ completely from those of Einstein theory in AdS
spaces. Besides, unlike Einstein gravity, the entropy does not obey the area
law. We also investigate the validity of the first law of thermodynamics as
well as the stability analysis in the canonical ensemble, and show that the
black string solutions are always thermodynamically stable in $f(R)$-Maxwell
theory with constant curvature scalar. Finally, we extend the study to the
case where the Ricci scalar is not a constant and in particular $R=R(r)$. In
this case, by using the Lagrangian multipliers method, we derive an analytical
black string solution from $f(R)$ gravity and reconstructed the function
$R(r)$. We find that this class of solutions has an additional logaritmic term
in the metric function which incorporates the effect of the $f(R)$ theory in
the solutions.

Keywords: modified gravity; string; thermodynamics.

##  I Introduction

There has been considerable attentions in the past years in modified gravity
theories, specially $f(R)$ theory which is one of the encouraging candidates
for explaining the current accelerating of the universe expansion Odin ; Capo
(see also Anto  for a comprehensive review on $f(R)$ theories). In fact $f(R)$
theories can be regarded as the simplest extension of general relativity. Many
$f(R)$ models have passed all the available experimental tests and fit the
cosmological data. To prevent a ghost state, $f^{\prime}(R)>0$ for $R\geq
R_{0}$ is required Nun ; Fara . $f^{\prime\prime}(R)>0$ for $R\geq R_{0}$, is
needed to avoid the negative mass squared of a scalar-field degree of freedom
(tachyon) Anto . $f(R)\rightarrow R-2\Lambda$ for $R\geq R_{0}$, is required
for the presence of the matter era and for consistency with local gravity
constraints Anto . It was shown that $f(R)$ theories can be considered as
general relativity with an additional scalar field that provide new insight in
the two cases of Brans-Dicke theory with ${\omega}_{0}=0$ and
${\omega}_{0}=-3/2$ Soti .

There have been a lot of works in the literature attempting to construct
static and stationary black hole solutions in $f(R)$ gravity theories. One may
expect that some signatures of black holes in $f(R)$ theories will be in
disagreement with the expected physical results of Einstein’s gravity. In Cruz
the authors studied general solutions in $f(R)$ theory using a perturbation
approach around the Einstein-Hilbert action. In Psal  black hole solutions
were found by adding dynamical vector and tensor degrees of freedom to the
Einstein-Hilbert action. Also, the transition from neutron stars to a strong
scalar- field state in $f(R)$ gravity has been studied in Nova . Physical
properties of the matter forming an accretion disk in the spherically
symmetric background in $f(R)$ theories were explored in Pun . In Ref. Lobo
the construction of traversable wormhole geometries was discussed in $f(R)$
gravity. The Schwarzschild-de Sitter black hole like solutions of $f(R)$
gravity were obtained for a positively constant and a non-constant curvature
scalar in Cogn  and Sebas , respectively. A black hole solution was obtained
from $f(R)$ theories by requiring the negative constant curvature scalar Cruz
. If $1+f^{\prime}(R_{0})>0$, this black hole is similar to the Schwarzschild-
AdS (SAdS) black hole. It was argued that $f(R)$ and SAdS black holes have no
big difference in thermodynamic quantities when using the Euclidean action
approach and replacing the Newtonian constant $G$ by $G_{\rm
eff}=G/(1+f^{\prime}(R_{0}))$ Cruz . It is also interesting to study black
hole solutions in $f(R)$ theory coupled to a matter field. In general, the
field equations of $f(R)$ theory coupled to the matter field are very
complicated and hence it is not easy to find exact analytical solutions. In
order to construct the constant curvature scalar black hole solutions from
$f(R)$ gravity coupled to the matter, the trace of its energy-momentum tensor
$T_{\mu\nu}$ should be zero Moon . Two examples for the traceless $T_{\mu\nu}$
are Maxwell and Yang-Mills fields which were studied in Moon ; Habib .
Thermodynamics and properties of these solutions were also studied in ample
details Moon . It was found that these solutions are similar to the Reissner-
Nordström–AdS (RNAdS) black hole when making appropriate replacements Moon .
The Kerr-Newman black hole solutions with non-zero constant scalar curvature
in $f(R)$-Maxwell theory, their thermodynamics, as well as their local and
global stability were also studied in Alex .

In this paper we would like to continue the investigation on the $f(R)$ black
holes, by constructing a new class of charged rotating black string solutions
in $R+f(R)$-Maxwell theory with constant curvature scalar. The traceless
property of the energy-momentum tensor of the Maxwell field plays a crucial
role in our derivation. With assumptions $R_{0}<0$ and $1+f^{\prime}(R_{0})>0$
our solution is similar to charged black string solution in AdS space with
suitable replacing the parameters. We will also suggest the suitable
counterterm which removes the divergences of the action. We calculate the
conserved and thermodynamic quantities of these black strings by using the
counterterm method. We obtain a Smarr-type formula for the mass of the black
string and check the validity of the first law of thermodynamics. We perform
the stability analysis in the canonical ensemble and show that the black
strings are always thermodynamically stable in $f(R)$-Maxwell theory with
constant curvature scalar. Finally, we extend the study to the case where the
Ricci scalar is not constant and in particular $R=R(r)$ and derive an
analytical black string solution.

##  II Field Equations and solutions

We start from the four-dimensional $R+f(R)$ theory coupled to the Maxwell
field

$\displaystyle I_{G}$ $\displaystyle=$
$\displaystyle-\frac{1}{16\pi}\int_{\mathcal{M}}d^{4}x\sqrt{-g}\left(R+f(R)-F_{\mu\nu}F^{\mu\nu}\right)-\frac{1}{8\pi}\int_{\partial\mathcal{M}}d^{3}x\sqrt{-h}\Theta(h),$
(1)

where ${R}$ is the Ricci scalar curvature,
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the
electromagnetic field tensor, and $A_{\mu}$ is the electromagnetic potential.
The last term in Eq. (1) is the Gibbons-Hawking boundary term. It is required
for the variational principle to be well-defined. The factor $\Theta$
represents the trace of the extrinsic curvature for the boundary
${\partial\mathcal{M}}$ and $h$ is the induced metric on the boundary. The
equations of motion can be obtained by varying the action (1) with respect to
the gravitational field $g_{\mu\nu}$ and the gauge field $A_{\mu}$ which
yields the following field equations

$\displaystyle{R}_{\mu\nu}\left(1+f^{\prime}(R)\right)-\frac{1}{2}g_{\mu\nu}(R+f(R))+\left(g_{\mu\nu}\nabla^{2}-\nabla_{\mu}\nabla_{\nu}\right)f^{\prime}(R)=8\pi
T_{\mu\nu},$ (2) $\displaystyle\nabla_{\mu}F^{\mu\nu}=0,$ (3)

with the energy-momentum tensor

$T_{\mu\nu}=\frac{1}{4\pi}\left(F_{\mu\eta}F_{\nu}^{\text{
}\eta}-\frac{1}{4}g_{\mu\nu}F_{\lambda\eta}F^{\lambda\eta}\right).$ (4)

The above energy-momentum tensor is traceless in four dimension, i. e.,
$T^{\mu}_{\text{ }\ \mu}=0$. As we mentioned already this property plays an
important role in our derivation. In Eq. (2) the “prime” denotes
differentiation with respect to curvature scalar $R$. Assuming the constant
curvature scalar $R=R_{0}$, the trace of Eq. (2) yields

$\displaystyle
R_{0}\left(1+f^{\prime}(R_{0})\right)-2\left(R_{0}+f(R_{0})\right)=0,$ (5)

Solving the above equation for negative $R_{0}$, gives

$\displaystyle R_{0}=\frac{2f(R_{0})}{f^{\prime}(R_{0})-1}\equiv
4{\Lambda_{\rm f}}<0.$ (6)

Substituting the above relation into Eq. (2), we obtain the following equation
for Ricci tensor

${R}_{\mu\nu}=\frac{1}{2}g_{\mu\nu}\left(\frac{f(R_{0})}{f^{\prime}(R_{0})-1}\right)+\frac{2}{1+f^{\prime}(R_{0})}T_{\mu\nu}.$
(7)

Now, we want to construct charged rotating black string solutions of the field
equations (2) and (3) and investigate their properties. We are looking for the
four-dimensional rotating solution with cylindrical or toroidal horizons. The
metric which describes such a spacetime can be written in the following form
Lem ; shey

$\displaystyle ds^{2}$ $\displaystyle=$ $\displaystyle-N(r)\left(\Xi dt-
ad\phi\right)^{2}+r^{2}\left(\frac{a}{l^{2}}dt-\Xi
d\phi\right)^{2}+\frac{dr^{2}}{N(r)}+\frac{r^{2}}{l^{2}}dz^{2},$
$\displaystyle\Xi^{2}$ $\displaystyle=$ $\displaystyle 1+\frac{a^{2}}{l^{2}},$
(8)

where $a$ is the rotation parameter. The function $N(r)$ should be determined
and $l$ has the dimension of length which is related to the constant
$\Lambda_{\rm f}$ by the relation $l^{2}=-3/\Lambda_{\rm f}$. The two
dimensional space, $t$=constant and $r$ =constant, can be (i) the flat torus
model $T^{2}$ with topology $S^{1}\times S^{1}$, and $0\leq\phi<2\pi$, $0\leq
z<2\pi l$, (ii) the standard cylindrical model with topology $R\times S^{1}$,
and $0\leq\phi<2\pi$, $-\infty<z<\infty$, and (iii) the infinite plane $R^{2}$
with $-\infty<\phi<\infty$ and $-\infty<z<\infty$. We will focus upon (i) and
(ii). The Maxwell equation (3) can be integrated immediately to give

$\displaystyle F_{tr}$ $\displaystyle=$ $\displaystyle\frac{q\Xi}{r^{2}},$
$\displaystyle F_{\phi r}$ $\displaystyle=$
$\displaystyle-\frac{a}{\Xi}F_{tr},$ (9)

where $q$ is the charge parameter of the black string. Substituting the
Maxwell fields (II) as well as the metric (II) in the field equation (2) with
constant curvature, the non-vanishing independent components of the field
equations for $a=0$ reduce to

$\displaystyle\left(1+{\it
f^{\prime}(R_{0})}\right)\left(2r^{4}\frac{d^{2}N(r)}{dr^{2}}+4r^{3}\frac{dN(r)}{dr}+R_{0}r^{4}\right)-4q^{2}=0,$
(10) $\displaystyle\left(1+{\it
f^{\prime}(R_{0})}\right)\left(4r^{3}\frac{dN(r)}{dr}+4r^{2}N(r)+R_{0}r^{4}\right)+4q^{2}=0.$
(11)

One can easily show that the above equations have the following solution

$N(r)=-\frac{2m}{r}+\frac{q^{2}}{(1+f^{\prime}(R_{0}))r^{2}}-\frac{R_{0}}{12}r^{2},$
(12)

where $m$ is an integration constant which is related to the mass of the
string. One can also check that these solutions satisfy equations (2)-(3) in
the rotating case where $a\neq 0$. It is apparent that this spacetime is
similar with asymptotically AdS black string. Indeed, with the following
replacement

$\displaystyle\frac{q^{2}}{\left(1+f^{\prime}(R_{0})\right)}\rightarrow Q^{2}$
(13) $\displaystyle\frac{R_{0}}{4}\rightarrow\Lambda$ (14)

the solution reduces to the asymptotically AdS charged black string for
$\Lambda=-3/l^{2}$ Lem . Next we study the physical properties of the
solutions. The Kretschmann scalar for this solution is given by

$\displaystyle
R_{\mu\nu\lambda\kappa}R^{\mu\nu\lambda\kappa}=\frac{8}{3r^{8}(1+f^{\prime}(R_{0}))^{2}}\left[r^{2}(\frac{1}{16}{R_{0}}^{2}r^{6}+18m^{2})(1+f^{\prime}(R_{0}))^{2}-36mrq^{2}(1+f^{\prime}(R_{0}))+21q^{4}\right].$
(15)

Figure 1: The function $N(r)$ versus $r$ for $m=2$, $f^{\prime}(R_{0})=2$ and
$q=1$. $R_{0}=12$ (bold line) and $R_{0}=-12$ (continuous line).

Figure 2: The function $N(r)$ versus $r$ for $m=2$, $q=1$ and $R_{0}=-12$.
$f^{\prime}(R_{0})=0.5$ (bold line), $f^{\prime}(R_{0})=0$ (continuous line)
and $f^{\prime}(R_{0})=-0.5$ (dashed line).

Figure 3: The function $N(r)$ versus $r$ for $m=2$, $f^{\prime}(R_{0})=-2$ and
$q=1$. $R_{0}=12$ (bold line) and $R_{0}=-12$ (continuous line).

When $r\rightarrow 0$, the dominant term in the Kretschmann scalar is
${56q^{4}}/[(1+f^{\prime}(R_{0}))^{2}r^{8}]$. Therefore we have an essential
singularity located at $r=0$. The Kretschmann scalar also approaches
${R_{0}}^{2}/6$ as $r\rightarrow\infty$. As one can see from Eq. (12), the
solution is ill-defined for $f^{\prime}(R_{0})=-1$. The cases with
$f^{\prime}(R_{0})>-1$ and $f^{\prime}(R_{0})<-1$ should be considered
separately. In the first case where $f^{\prime}(R_{0})>-1$, there exist a
cosmological horizon for $R_{0}>0$, while there is no cosmological horizons if
$R_{0}<0$ (see fig. 1). Indeed, for $1+f^{\prime}(R_{0})>0$ and $R_{0}<0$ the
black string can have two inner and outer horizons provided the parameters of
the solutions are chosen suitably (see fig. 2). In the latter case
($f^{\prime}(R_{0})<-1$), the signature of the spacetime changes and the
conserved quantities such as mass and angular momenta become negative, as we
will see in the next section, thus this is not a physical case and we rule it
out from our consideration (see fig. 3 ).

##  III Conserved and Thermodynamic quantities

Next, we calculate the conserved quantities of the solutions by using the
counterterm method inspired by (A)dS/CFT correspondence Mal . The spacetimes
under consideration in this paper has zero curvature boundary,
$R_{abcd}(h)=0$, and therefore the counterterm for the stress energy tensor
should be proportional to $h^{ab}$. We find the suitable counterterm which
removes the divergences of the action in the form

$I_{\rm
ct}=-\frac{1}{8\pi}\int_{\partial\mathcal{M}}d^{3}x\sqrt{-h}\sqrt{-\frac{R_{0}}{3}},$
(16)

where $R_{0}<0$. Having the total finite action $I=I_{G}+I_{\mathrm{ct}}$ at
hand, one can use the quasilocal definition to construct a divergence free
stress-energy tensor BY . Thus the finite stress-energy tensor can be written
as

$T^{ab}=\frac{1}{8\pi}\left[\Theta^{ab}-\Theta
h^{ab}-\sqrt{-\frac{R_{0}}{3}}h^{ab}\right].$ (17)

The first two terms in Eq. (17) are the variation of the action (1) with
respect to $h_{ab}$, and the last term is the variation of the boundary
counterterm (16) with respect to $h_{ab}$. To compute the conserved charges of
the spacetime, one should choose a spacelike surface $\mathcal{B}$ in
$\partial\mathcal{M}$ with metric $\sigma_{ij}$, and write the boundary metric
in ADM (Arnowitt-Deser-Misner) form:

$h_{ab}dx^{a}dx^{a}=-N^{2}dt^{2}+\sigma_{ij}\left(d\varphi^{i}+V^{i}dt\right)\left(d\varphi^{j}+V^{j}dt\right),$

where the coordinates $\varphi^{i}$ are the angular variables parameterizing
the hypersurface of constant $r$ around the origin, and $N$ and $V^{i}$ are
the lapse and shift functions respectively. When there is a Killing vector
field $\mathcal{\xi}$ on the boundary, then the quasilocal conserved
quantities associated with the stress tensors of Eq. (17) can be written as

$Q(\mathcal{\xi)}=\int_{\mathcal{B}}d^{2}x\sqrt{\sigma}T_{ab}n^{a}\mathcal{\xi}^{b},$
(18)

where $\sigma$ is the determinant of the metric $\sigma_{ij}$, $\mathcal{\xi}$
and $n^{a}$ are, respectively, the Killing vector field and the unit normal
vector on the boundary $\mathcal{B}$. The first Killing vector of the
spacetime is $\xi=\partial/\partial t$, and therefore its associated conserved
charge of the string is the mass per unit volume. A simple calculation gives

$\displaystyle M$ $\displaystyle=$
$\displaystyle\int_{\mathcal{B}}d^{2}x\sqrt{\sigma}T_{ab}n^{a}\xi^{b}=\frac{(3\Xi^{2}-1)m}{8\pi
l}\left[1+f^{\prime}(R_{0})\right].$ (19)

The second conserved quantity is the angular momentum per unit volume
associated with the rotational Killing vectors
$\varsigma=\partial/\partial\phi$ which can be calculated as

$\displaystyle
J=\int_{\mathcal{B}}d^{2}x\sqrt{\sigma}T_{ab}n^{a}\varsigma^{b}=\frac{3\Xi
m\sqrt{\Xi^{2}-1}}{8\pi}\left[1+f^{\prime}(R_{0})\right].$ (20)

For $a=0$ ($\Xi=1$), the angular momentum per unit volume vanishes, and
therefore $a$ is the rotational parameters of the spacetime. Next we calculate
the entropy of the black string. Let us first give a brief discussion
regarding the entropy of the black hole solutions in $f(R)$ gravity. To this
aim, we follow the arguments presented in Brevik . If one use the Noether
charge method for evaluating the entropy associated with black hole solutions
in $f(R)$ theory with constant curvature, one finds Cogn

${S}=\frac{A}{4G}f^{\prime}(R_{0}),$ (21)

where $A=4\pi r_{+}^{2}$ is the horizon area. As a result, in $f(R)$ gravity,
the entropy does not obey the area law and one obtains a modification of the “
area law”. Motivated by the above argument, for the rotating black string
solution in $R+f(R)$ gravity, we find the entropy per unit length of the
string as

${S}=\frac{r_{+}^{2}\Xi}{4l}\left[1+f^{\prime}(R_{0})\right].$ (22)

Then we obtain the temperature and angular velocity of the horizon by analytic
continuation of the metric. Although our solution is not static, the Killing
vector

$\chi=\partial_{t}+\Omega\partial_{\phi}$ (23)

is the null generator of the event horizon where $\Omega$ is the angular
velocity of the outer horizon. The analytical continuation of the Lorentzian
metric by $t\rightarrow i\tau$ and $a\rightarrow ia$ yields the Euclidean
section, whose regularity at $r=r_{+}$ requires that we should identify
$\tau\sim\tau+\beta_{+}$ and $\phi\sim\phi+i\beta_{+}\Omega_{+}$ where
$\beta_{+}$ and $\Omega_{+}$ are the inverse Hawking temperature and the
angular velocity of the horizon. We find

$\displaystyle T$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi\Xi}\left(\frac{dN(r)}{dr}\right)_{r=r_{+}}=-\frac{\left[R_{0}r_{+}^{4}(1+f^{\prime}(R_{0}))+4q^{2}\right]}{16\pi\Xi[1+f^{\prime}(R_{0})]r_{+}^{3}},$
(24) $\displaystyle\Omega$ $\displaystyle=$ $\displaystyle\frac{a}{\Xi
l^{2}},$ (25)

where we have used equation $N(r_{+})=0$ for omitting the mass parameter $m$
from temperature expression. Since $1+f^{\prime}(R_{0})>0$, therefore the
temperature is non negative provided

$\displaystyle R_{0}r_{+}^{4}(1+f^{\prime}(R_{0}))\leq-4q^{2}\rightarrow
R_{0}\leq-\frac{4q^{2}}{r_{+}^{4}(1+f^{\prime}(R_{0}))},$ (26)

where the equality holds for extremal black string with zero temperature. The
next quantity we are going to calculate is the electric charge of the string.
To determine the electric field we should consider the projections of the
electromagnetic field tensor on special hypersurface. The normal vectors to
such hypersurface are

$u^{0}=\frac{1}{N},\text{ \ }u^{r}=0,\text{ \ }u^{i}=-\frac{V^{i}}{N},$ (27)

where $N$ and $V^{i}$ are the lapse function and shift vector. Then the
electric field is $E^{\mu}=g^{\mu\rho}F_{\rho\nu}u^{\nu}$, and the electric
charge per unit length of the string can be found by calculating the flux of
the electric field at infinity,

${Q}=\frac{\Xi q}{4\pi l\sqrt{1+f^{\prime}(R_{0})}}.$ (28)

The electric potential $U$, measured at infinity with respect to the horizon,
is defined by Cal

$U=A_{\mu}\chi^{\mu}\left|{}_{r\rightarrow\infty}-A_{\mu}\chi^{\mu}\right|_{r=r_{+}},$
(29)

where $\chi$ is the null generator of the event horizon given in Eq. (23). One
can easily obtain the electric potential as

$U=\frac{q}{\Xi r_{+}}\sqrt{1+f^{\prime}(R_{0})}.$ (30)

Then, we consider the first law of thermodynamics for the black string. In
order to do this, we obtain the mass $M$ as a function of extensive quantities
$S$, ${J}$ and $Q$. Using the expression for the mass, the angular momenta,
the entropy and the charge given in Eqs. (19), (20), (22) and (28) and the
fact that $N(r_{+})=0$, one can obtain a Smarr-type formula as

$M(S,J,Q)=\frac{J(3Z-1)}{3l\sqrt{Z(Z-1)}},$ (31)

where $Z=\Xi^{2}$ is the positive real root of the following equation:

$\displaystyle\frac{3\sqrt{Z-1}\pi^{2}Q^{2}l^{2}[1+f^{\prime}(R_{0})]^{2}-2J\pi\sqrt{\sqrt{Z}(1+f^{\prime}(R_{0}))Sl}+3S^{2}\sqrt{Z-1}}{2\pi\sqrt{\sqrt{Z}(1+f^{\prime}(R_{0}))Sl}}=0.$
(32)

One may then regard the parameters $S$, ${J}$ and $Q$ as a complete set of
extensive parameters for the mass $M(S,{J},Q)$ and define the intensive
parameters conjugate to $S$, ${J}$ and $Q$. These quantities are the
temperature, the angular velocities and the electric potential

$\displaystyle T$ $\displaystyle=$ $\displaystyle\left(\frac{\partial
M}{\partial S}\right)_{J,Q}=-\frac{J\left\\{\left[\pi
Ql(1+f^{\prime}(R_{0}))\right]^{2}-3S^{2}\right\\}}{3Sl\sqrt{Z(Z-1)}\left\\{\left[\pi
Ql(1+f^{\prime}(R_{0}))\right]^{2}+S^{2}\right\\}},$ (33)
$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle\left(\frac{\partial
M}{\partial J}\right)_{S,Q}$ (34) $\displaystyle=$
$\displaystyle\frac{3\left(3Z-1\right)\left\\{\left[\pi
Ql(1+f^{\prime}(R_{0}))\right]^{2}+S^{2}\right\\}-4\pi
J\sqrt{\sqrt{Z}(Z-1)(1+f^{\prime}(R_{0}))Sl}}{9l\sqrt{Z(Z-1)}\left\\{\left[\pi
Ql(1+f^{\prime}(R_{0}))\right]^{2}+S^{2}\right\\}},$ $\displaystyle U$
$\displaystyle=$ $\displaystyle\left(\frac{\partial M}{\partial
Q}\right)_{S,J}=\frac{4\pi^{2}QlJ[1+f^{\prime}(R_{0})]^{2}}{3\sqrt{Z(Z-1)}\left\\{\left[\pi
Ql(1+f^{\prime}(R_{0}))\right]^{2}+S^{2}\right\\}}.$ (35)

Numerical calculations show that the intensive quantities calculated by Eqs.
(33)-(35) coincide with Eqs. (24), (25) and (30), respectively. Thus, these
thermodynamics quantities satisfy the first law of thermodynamics

$dM=TdS+\Omega d{J}+Ud{Q}.$ (36)

##  IV Thermal Stability of black string

Figure 4: The function $(\partial^{2}M/\partial S^{2})_{J,Q}$ versus $q$ for
$l=1$, $\Xi=1.25$, $r_{+}=0.7$ and $R_{0}=-12$. $f^{\prime}(R_{0})=0$ (bold
line), $f^{\prime}(R_{0})=1$ (continuous line) and $f^{\prime}(R_{0})=2$
(dashed line).

Figure 5: The function $(\partial^{2}M/\partial S^{2})_{J,Q}$ versus $r_{+}$
for $l=1$, $\Xi=1.25$, $R_{0}=-12$ and $f^{\prime}(R_{0})=1$. $q=0.5$ (bold
line), $q=1$ (continuous line) and $q=1.5$ (dashed line).

Figure 6: The function $(\partial^{2}M/\partial S^{2})_{J,Q}$ versus $q$ for
$l=1$, $f^{\prime}(R_{0})=1$ and $R_{0}=-12$. $\Xi=1.25$, (bold line),
$\Xi=1.75$, (continuous line) and $\Xi=2.25$, (dashed line).

Finally, we investigate the thermal stability of rotating black string
solutions in $f(R)$ gravity coupled to a matter field. The stability of a
thermodynamic system with respect to small variations of the thermodynamic
coordinates is usually performed by analyzing the behavior of the entropy
$S(M,J,Q)$ around the equilibrium. The local stability in any ensemble
requires that $S(M,{J},Q)$ be a convex function of the extensive variables or
its Legendre transformation must be a concave function of the intensive
variables. The stability can also be studied by the behavior of the energy
$M(S,J,Q)$ which should be a convex function of its extensive variables. Thus,
the local stability can in principle be carried out by finding the determinant
of the Hessian matrix of $M(S,{J},Q)$ with respect to its extensive variables
$X_{i}$, $\mathbf{H}_{X_{i}X_{j}}^{M}=[\partial^{2}M/\partial X_{i}\partial
X_{j}]$ Cal ; Gub . In our case the mass $M$ is a function of entropy, angular
momenta, and charge. The number of thermodynamic variables depends on the
ensemble that is used. In the canonical ensemble, the charge and the angular
momenta are fixed parameters, and therefore the positivity of the
$(\partial^{2}M/\partial S^{2})_{{J},Q}$ is sufficient to ensure local
stability. We find that the black string solutions are always thermally stable
independent of the value of the parameters $q$ and $\Xi$. We have shown the
behavior of $(\partial^{2}M/\partial S^{2})_{{J},Q}$ as a function $q$ and
$r_{+}$ for different value of $\Xi$ and $f^{\prime}(R_{0})$ in figures 4-6.
These figures show that the black string solutions in $f(R)$-Maxwell theory
with constant curvature scalar are always thermally stable.

##  V Solution with non constant Ricci scalar

In this section we would like to extend the study to the case where the Ricci
scalar is not a constant, instead we reconstruct it as $R=R(r)$ as a result of
our calculations. In this case, we find out that we can obtain solution only
in the absence of the matter field. As we mentioned in the introduction, in
general, the field equations of $f(R)$ theory coupled to the matter field are
very complicated and hence it is not easy to find exact analytical solutions.
Thus we only consider the uncharged black string solution. We start with the
following action

$S=\frac{1}{16\pi}\int d^{4}x\sqrt{-g}f(R)\,,$ (37)

where $f(R)$ is a generic function of the Ricci scalar $R$. We also modify our
metric (II) a bit as follow

$\displaystyle ds^{2}$ $\displaystyle=$
$\displaystyle-N(r)\mathrm{e}^{2\alpha(r)}\left(\Xi dt-
ad\phi\right)^{2}+r^{2}\left(\frac{a}{l^{2}}dt-\Xi
d\phi\right)^{2}+\frac{dr^{2}}{N(r)}+\frac{r^{2}}{l^{2}}dz^{2},$ (38)

where we have added an additional function $\alpha(r)$ in the metric
coefficients. For simplicity we only consider the non-rotating black string
with $a=0$, thus the above metric reduces to

$\displaystyle
ds^{2}=-N(r)\mathrm{e}^{2\alpha(r)}dt^{2}+\frac{dr^{2}}{N(r)}+r^{2}d\phi^{2}+\frac{r^{2}}{l^{2}}dz^{2}.$
(39)

The scalar curvature for metric (39) reads

$\displaystyle R$ $\displaystyle=$
$\displaystyle-3\,{\frac{dN\left(r\right)}{dr}}{\frac{d\alpha\left(r\right)}{dr}}-2\,N\left(r\right)\left[{\frac{d}{dr}}\alpha\left(r\right)\right]^{2}-{\frac{d^{2}N\left(r\right)}{d{r}^{2}}}-2\,N\left(r\right){\frac{d^{2}\alpha\left(r\right)}{d{r}^{2}}}$
(40)
$\displaystyle-\frac{4}{r}\frac{dN\left(r\right)}{dr}-\frac{4N\left(r\right)}{r}\frac{d\alpha\left(r\right)}{dr}-\frac{2N\left(r\right)}{{r}^{2}}\,.$

We use the Lagrangian multipliers method. In the framework of Friedmann-
Robertson-Walker universe this method was studied in vile ; Capozziello ;
Monica , while for static spherically symmetric black hole solutions it was
investigated in Sebas ; Capozziello2 . In this approach one may consider the
scalar curvature $R$ as independent Lagrangian coordinates in addition to the
functions $\alpha(r)$ and $N(r)$, which appear from the metric line element.

Introducing the Lagrangian multipliers $\lambda$, after using (40), the action
(37) can be written

$\displaystyle S$ $\displaystyle\equiv$ $\displaystyle\frac{1}{16\pi}\int
dt\int
d{r}\left(e^{\alpha(r)}r^{2}\right)\left\\{f(R)-\lambda\left\\{R+\left\\{3\,\left[{\frac{d}{dr}}N\left(r\right)\right]{\frac{d}{dr}}\alpha\left(r\right)\right.\right.\right.$
(41)
$\displaystyle+2\,N\left(r\right)\left[{\frac{d}{dr}}\alpha\left(r\right)\right]^{2}+{\frac{d^{2}N\left(r\right)}{d{r}^{2}}}+2\,N\left(r\right)\frac{d^{2}\alpha\left(r\right)}{d{r}^{2}}+\frac{4}{r}\frac{dN\left(r\right)}{dr}$
$\displaystyle+\frac{4N\left(r\right)}{r}\frac{d\alpha\left(r\right)}{dr}+{\frac{2N\left(r\right)}{{r}^{2}}}\left.\left.\left.\right\\}\right\\}\right\\}\,.$

Varying the above action with respect to $R$, one gets

$\lambda=f^{\prime}(R),$ (42)

where the prime denotes the derivative with respect to the scalar curvature
$R$. Substituting this value and integrating by part, the Lagrangian takes the
form

$\displaystyle L(\alpha,d\alpha/dr,N,dN/dr,R,dR/dr)$ $\displaystyle=$
$\displaystyle
e^{\alpha}\left\\{r^{2}\left[f(R)-Rf^{\prime}(R)\right]-2f^{\prime}(R)\left(r\frac{dN(r)}{dr}+N(r)\right)\right.$
(43)
$\displaystyle+\left.f^{\prime\prime}(R)\frac{dR}{dr}r^{2}\left(\frac{dN(r)}{dr}+2N(r)\frac{d\alpha(r)}{dr}\right)\right\\}\,.$

Making the variation with respect to $\alpha$, one gets the first equation of
motion

$\displaystyle\frac{Rf^{\prime}(R)-f(R)}{f^{\prime}(R)}+\frac{2}{r^{2}}\left[N(r)+r\frac{dN(r)}{dr}\right]$
(44)
$\displaystyle+\frac{2N(r)f^{\prime\prime}(R)}{f^{\prime}(R)}\left[\frac{d^{2}R}{dr^{2}}+\left(\frac{dN(r)/dr}{2N(r)}\right)\frac{dR}{dr}+\frac{f^{\prime\prime\prime}(R)}{f^{\prime\prime}(R)}\left(\frac{dR}{dr}\right)^{2}\right]=0\,.$

The variation with respect to $N(r)$ leads the second equation of motion

$\left[\frac{d\alpha(r)}{dr}\left(\frac{f^{\prime\prime}(R)}{f^{\prime}(R)}\frac{dR}{dr}\right)-\frac{f^{\prime\prime}(R)}{f^{\prime}(R)}\frac{d^{2}R}{dr^{2}}-\frac{f^{\prime\prime\prime}(R)}{f^{\prime}(R)}\left(\frac{dR}{dr}\right)^{2}\right]=0\,,$
(45)

while by making the variation with respect to $R$, we recover Eq. (40). Given
$f(R)$, together with equation (40), the above equations form a system of
three differential equations for the three unknown quantities $\alpha(r),N(r)$
and $R(r)$. We would like to note that one advantage of this approach is that
$\alpha$ does not appear in Eq.(44). In what follow, we will find exact
solutions of the above system of differential equations.

In the special case of constant curvature $R=R_{0}$ and $\alpha=\rm constant$,
it is easy to show that the only solution of Eqs. (40) and(44) is the
Schwarzshild de Sitter black string solution with flat horizon,

$N(r)=-\frac{2m}{r}-\frac{\Lambda}{3}r^{2},$ (46)

where $m$ is a constant of integration which can be interpreted as the mass
parameter of the black string and we have defined $2\Lambda\equiv
R_{0}-f(R_{0})/f^{\prime}(R_{0})$ and $R_{0}=4\Lambda$. Notice that in the
absence of the matter field, $q=0$, solution (12) coincides with the result
obtained in (46), as expected.

Next, we consider the case of non constant Ricci curvature, but still with
$\alpha=\rm constant$. From Eq. (45) we have

$f^{\prime\prime\prime}\left(\frac{dR}{dr}\right)^{2}+f^{\prime\prime}\left(\frac{d^{2}R}{dr^{2}}\right)=\frac{d^{2}}{dr^{2}}f^{\prime}(R)=0,$
(47)

which has the following solution,

$f^{\prime}(R)=mr+n,$ (48)

where $m$ and $n$ are two integration constants. Given the explicit form of
$R$, we may find $r$ as a function of Ricci scalar and reconstruct
$f^{\prime}(R)$ realizing such solution. From Eq. (40) with constant $\alpha$,
one gets

$R=-{\frac{d^{2}N\left(r\right)}{d{r}^{2}}}-\frac{4}{r}\frac{dN\left(r\right)}{dr}-\,{\frac{2N\left(r\right)}{{r}^{2}}}.$
(49)

Using the fact that $(f^{\prime\prime}(R))dR/dr=df^{\prime}(R)/dr=m$ and
$df(R)/dr=f^{\prime}(R)dR/dr$, and multiplying Eq. (44) by $f^{\prime}(R)$, we
arrive at

$-\frac{d^{2}N(r)}{dr^{2}}\left(m+\frac{n}{r}\right)+\frac{4mN(r)}{r^{2}}+\frac{2nN(r)}{r^{3}}-\frac{m}{r}\frac{dN(r)}{dr}=0\,.$
(50)

When $m=0$, the solution of the above equation is ones obtained in (46). For
$n=0$, the general solution is

$N(r)=-C_{1}r^{2}+\frac{C_{2}}{r^{2}}.$ (51)

Substituting in Eq. (49) we again arrive at constant Ricci scalar,
$R=12C_{1}$. Although in this case $f^{\prime}(R)=df(R)/dR=mr$ is not a
constant, but still we have $df(R)/dr=0$, which implies that $f(R)=\rm
constant$. Next we look for the most general solution of Eq.(50) with $n\neq
0$ and $m\neq 0$. Solving (50), we find

$\displaystyle
N\left(r\right)=-C_{1}{r}^{2}+\frac{C_{2}}{r}\left[2\,{n}^{3}-3mn^{2}r+6\,{m}^{2}{r}^{2}n-6{m}^{3}r^{3}\ln\left(m+\frac{n}{r}\right)\right],$
(52)

where $C_{1}$ and $C_{2}$ are two arbitrary constants. Given solution (52) one
can basically construct $f(R)$ by using Eqs. (48) and (49). In order to
simplify the above solution, we choose $n=1$ and $C_{2}=-1/m$,

$N\left(r\right)=3-\frac{2}{mr}-6\,{mr}-C_{1}{r}^{2}+6{m}^{2}r^{2}\ln\left(m+\frac{1}{r}\right).$
(53)

For this general case the Ricci scalar becomes

${R(r)}=12C_{1}-72m^{2}\ln\left(m+\frac{1}{r}\right)+\frac{6\left(12m^{3}r^{3}+18m^{2}r^{2}+4mr-1\right)}{r^{2}(mr+1)^{2}}.$
(54)

which is clearly not a constant. Now we want to reconstruct the corresponding
$f(R)$ theory. From Eq.(48), for $n=1$ one has

$f^{\prime}(R)=\frac{df(R)}{dR}=\frac{df(R)}{dr}\frac{dr}{dR}=mr+1.$ (55)

Integrating (55), by using (54), we get

$f[R(r)]=-36m^{2}\ln\left(m+\frac{1}{r}\right)+\frac{3\left(12m^{3}r^{3}+18m^{2}r^{2}+4mr-2\right)}{r^{2}(ar+1)^{2}}.$
(56)

Combining Eqs. (54), (55) and (56), one gets the following differential
equation for function $f(R)$,

$\frac{3m^{2}}{[f^{\prime}(R)]^{2}[f^{\prime}(R)-1]^{2}}+f(R)-\frac{R}{2}+6C_{1}=0.$
(57)

This equation has a simple solution as

$f(R)=\frac{R}{2}-6C_{1}-48m^{2},$ (58)

but it has also another complicated solution which we have not presented it
here. Thus we have found a black string solution in $f(R)$ gravity with non
constant Ricci scalar. This approach also leads to construct Ricci scalar as a
function of $r$, as given in Eq. (54). The obtained solutions in this section
differ from that presented in Capozziello2  for axially symmetric solutions in
$f(R)$ gravity. It is worth mentioning that metric (38) has a good property
for which its static and rotating solutions coincide and so in this section we
only study the static case. Following the approach of this section, one can
easily check that solution (52) can be deduced for rotating case where $a\neq
0$. Besides, in this section we only considered the case with $\alpha=\rm
constant$, and derived the metric function (52) as well as $R(r)$. The study
can also be generalized to the case where $\alpha=\alpha(r)$. We leave it and
also thermodynamic considerations of the obtained solution in this section,
for future investigations.

##  VI Conclusions

In order to obtain the constant curvature black hole solution in $f(R)$
gravity theory coupled to a matter field, the trace of the energy-momentum
tensor of the matter field should be zero Moon . Since the energy-momentum
tensor of Maxwell field is traceless in four dimensions, therefore spherically
symmetric black hole solutions from $f(R)$ theory coupled to Maxwell field was
derived in four dimensional spacetime Moon .

In this paper we continued the study by constructing a new class of charged
rotating solutions in $f(R)$-Maxwell theory with constant curvature scalar.
This class of solutions describe the four dimensional charged rotation black
string with cylindrical or toroidal horizons with zero curvature boundary.
These solutions are similar to asymptotically AdS black string of Einstein-
Maxwell gravity with suitably replacement of the parameters. However, the
solution presented in this paper has at least two differences from AdS black
string solutions of Einstein-Maxwell gravity. First, the conserved and
thermodynamic quantities computed here depend on function $f^{\prime}(R_{0})$
and differ completely from those of Einstein theory in AdS spaces. Clearly the
presence of the general function $f^{\prime}(R_{0})$ changes the physical
values of conserved and thermodynamic quantities. Second, unlike Einstein
gravity, the entropy does not obey the area law for black string solutions in
$f(R)$-Maxwell theory as one can see from Eq (22). We studied the physical
properties of the solutions and found a suitable conterterm which removes the
divergence of the action. We obtained mass and angular momenta of the string
through the use of conterterm method. We also derived the entropy of the black
string in $f(R)$ gravity which has a modification from the area law. We
obtained a Smarr-type formula for the mass namely $M(S,J,Q)$ and checked that
the obtained conserved and thermodynamic quantities satisfy the first law of
black hole thermodynamics. Finally, we explored the thermal stability of the
solutions in the canonical ensemble and showed that the black strings derived
from $f(R)$\- Maxwell theory are always thermally stable. This is commensurate
with the fact that there is no Hawking-Page phase transition for black objects
with zero curvature horizon Wit .

We also extend the study to the case where the Ricci scalar is not a constant.
For this purpose we used the Lagrangian multipliers method and found an exact
black string solution in $f(R)$ gravity. In this approach one may consider the
scalar curvature $R(r)$ as an independent Lagrangian coordinates in addition
to the metric functions and deduce $R(r)$ as a solution of the field
equations. We found the explicit form of $R(r)$ as well as the metric function
which has a logaritmic term. It is worth noting that since for the non
constant Ricci scalar, the field equations of $f(R)$ theory coupled to the
matter field become very complicated, in this case, we could only derived
analytical solution in the absence of the matter field.

###### Acknowledgements.

We thank the referee for constructive comments which helped us to improve the
paper significantly. We also grateful to S. H. Hendi for useful comments and
helpful discussions. The work of A. Sheykhi has been supported financially by
Research Institute for Astronomy & Astrophysics of Maragha (RIAAM), Iran. A.
Sheykhi also thank from the Research Council of Shiraz University.

## References

  * (1) S. Nojiri and S. D. Odintsov, Phys. Rev. D 74, 086005 (2006);
S. Nojiri and S. D. Odintsov, Int. J. Geom. Meth. Mod. Phys. 4, 115 (2007);
S. Nojiri and S. D. Odintsov, J. Phys. Conf. Ser. 66, 012005 (2007);
S. Nojiri and S. D. Odintsov, Phys. Rept.505, 59 (2011);
S. Capozziello, S. Nojiri, S. D. Odintsov and A. Troisi, Phys. Lett. B 639,
135 (2006);
S. Nojiri and S. D. Odintsov, Phys. Rev. D 78, 046006 (2008).

  * (2) S. Capozziello, V. F. Cardone, and A. Troisi, J. Cosmol. Astropart. Phys. 08, 001 (2006);
S. Capozziello, V. F. Cardone, and A. Troisi, Mon. Not. R. Astron. Soc. 375,
1423 (2007);
K. Atazadeh, M. Farhoudi and H. R. Sepangi, Phys. Lett. B 660, 275 (2008);
C. Corda, Astropart. Phys. 34, 587 (2011).

  * (3) A. De Felice, S. Tsujikawa, Living Rev. Rel. 13 (2010) 3.
  * (4) A. Nunez and S. Solganik, hep-th/0403159.
  * (5) V. Faraoni, Phys. Rev. D 74, 104017 (2006).
  * (6) T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82, 451 (2010).
  * (7) A. de la Cruz-Dombriz, A. Dobado and A. L. Maroto, Phys. Rev. D 80, 124011 (2009).
  * (8) D. Psaltis, D. Perrodin, K. R. Dienes and I. Mocioiu, Phys. Rev. Lett. 100, 091101 (2008).
  * (9) J. Novak, Phys. Rev. D 58, 064019 (1998).
  * (10) C. S. J. Pun, Z. Kovacs and T. Harko, Phys. Rev. D 78, 024043 (2008).
  * (11) F. S. N. Lobo and M. A. Oliveira, Phys. Rev. D 80, 104012 (2009).
  * (12) G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov and S. Zerbini, JCAP 0502, 010 (2005).
  * (13) L. Sebastiani and S. Zerbini, arXiv:1012.5230.
  * (14) T. Moon, Y. S. Myung and E. J. Son, Gen. Rel. Grav. 43, 3079 (2011), arXiv:1101.1153.
  * (15) S. Habib Mazharimousavi, M. Halilsoy, T. Tahamtan, Eur. Phys. J. C 72 (2012) 1958.
  * (16) A. Larranaga, Pramana Journal of Physics 78, 697 (2012) ;
J. A. R. Cembranos, et al., arXiv:1109.4519.

  * (17) J. P. S. Lemos, Class. Quantum Gravit. 12, 1081 (1995);
J. P. S. Lemos, Phys. Lett. B 353, 46 (1995).

  * (18) M. H. Dehghani, Phys. Rev. D 66, 044006 (2002);
M. H. Dehghani and A. Khodam-Mohammadi, Phys. Rev. D 67, 084006 (2003);
A. M. Awad, Class. Quant. Grav. 20, 2827 (2003);
M. H. Dehghani and N. Farhangkhah, Phys. Rev. D 71, 044008 (2005);
A. Sheykhi, M. H. Dehghani, N. Riazi, and J. Pakravan, Phys. Rev. D 74, 084016
(2006);
A. Sheykhi, Phys. Rev. D 78, 064055 (2008).

  * (19) J. Maldacena, Adv. Theor. Math. Phys. 2, 231 (1998);
E. Witten, Adv. Theor. Math. Phys. 2, 253 (1998);
O. Aharony, S. S. Gubser, J. Maldacena, H. Ooguri and Y. Oz, Phys. Rep. 323,
183 (2000).

  * (20) J. D. Brown and J. W. York, Phys. Rev. D 47, 1407 (1993).
  * (21) I. Brevik, S. Nojiri, S.D. Odintsov and L. Vanzo, Phys. Rev. D 70, 043520 (2004).
  * (22) M. Cvetic and S. S. Gubser, J. High Energy Phys. 04, 024 (1999);
M. M. Caldarelli, G. Cognola and D. Klemm, Class. Quant. Gravit. 17, 399
(2000).

  * (23) S. S. Gubser and I. Mitra, J. High Energy Phys., 08, 018 (2001).
  * (24) E. Witten, Adv. Theor. Math. Phys. 2, 505 (1998).
  * (25) A. Vilenkin, Phys. Rev. D 32 2511 (1985).
  * (26) S. Capozziello, Int. J. Mod. Phys. D 11, 483, (2002)
  * (27) G. Cognola, M. Gastaldi and S. Zerbini, Int. J. Theor. Phys. 47, 898 (2008).
  * (28) S. Capozziello and M. De Laurentis, Phys. Rept. 509, 167 (2011).


Lattice Parton Collaboration ($\rm{\bf LPC}$)

# Nonperturbative Determination of Collins-Soper Kernel from Quasi Transverse-
Momentum Dependent Wave Functions

Min-Huan Chu  INPAC, Key Laboratory for Particle Astrophysics and Cosmology
(MOE), Shanghai Key Laboratory for Particle Physics and Cosmology, School of
Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
Yang Yuanqing Scientific Computering Center, Tsung-Dao Lee Institute, Shanghai
Jiao Tong University, Shanghai 200240, China  Zhi-Fu Deng  INPAC, Key
Laboratory for Particle Astrophysics and Cosmology (MOE), Shanghai Key
Laboratory for Particle Physics and Cosmology, School of Physics and
Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China  Jun Hua
Guangdong Provincial Key Laboratory of Nuclear Science, Institute of Quantum
Matter, South China Normal University, Guangzhou 510006, China  Guangdong-Hong
Kong Joint Laboratory of Quantum Matter, Southern Nuclear Science Computing
Center, South China Normal University, Guangzhou 510006, China  INPAC, Key
Laboratory for Particle Astrophysics and Cosmology (MOE), Shanghai Key
Laboratory for Particle Physics and Cosmology, School of Physics and
Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China  Xiangdong Ji
Department of Physics, University of Maryland, College Park, MD 20742, USA
Center for Nuclear Femtography, 1201 New York Ave., NW, Washington DC, 20005,
USA  Key Laboratory for Particle Astrophysics and Cosmology (MOE), Shanghai
Key Laboratory for Particle Physics and Cosmology, Tsung-Dao Lee Institute,
Shanghai Jiao Tong University, Shanghai 200240, China  Andreas Schäfer
Institut für Theoretische Physik, Universität Regensburg, D-93040 Regensburg,
Germany  Yushan Su  Department of Physics, University of Maryland, College
Park, MD 20742, USA  Peng Sun  Nanjing Normal University, Nanjing, Jiangsu,
210023, China  Wei Wang  INPAC, Key Laboratory for Particle Astrophysics and
Cosmology (MOE), Shanghai Key Laboratory for Particle Physics and Cosmology,
School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China  Yi-Bo Yang  CAS Key Laboratory of Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China  School of Fundamental Physics and Mathematical Sciences, Hangzhou
Institute for Advanced Study, UCAS, Hangzhou 310024, China  International
Centre for Theoretical Physics Asia-Pacific, Beijing/Hangzhou, China  School
of Physical Sciences, University of Chinese Academy of Sciences, Beijing
100049, China  Jun Zeng  INPAC, Key Laboratory for Particle Astrophysics and
Cosmology (MOE), Shanghai Key Laboratory for Particle Physics and Cosmology,
School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China  Jialu Zhang  INPAC, Key Laboratory for Particle Astrophysics
and Cosmology (MOE), Shanghai Key Laboratory for Particle Physics and
Cosmology, School of Physics and Astronomy, Shanghai Jiao Tong University,
Shanghai 200240, China  Jian-Hui Zhang  Center of Advanced Quantum Studies,
Department of Physics, Beijing Normal University, Beijing 100875, China  Qi-An
Zhang  School of Physics, Beihang University, Beijing 102206, China  Key
Laboratory for Particle Astrophysics and Cosmology (MOE), Shanghai Key
Laboratory for Particle Physics and Cosmology, Tsung-Dao Lee Institute,
Shanghai Jiao Tong University, Shanghai 200240, China

###### Abstract

In the framework of large-momentum effective theory, we present a first-
principle determination of the Collins-Soper kernel which governs the rapidity
evolution of transverse-momentum dependent (TMD) distributions. We simulate
the quasi TMD wave functions on a Euclidean lattice, and remove the pertinent
linear divergences with Wilson loop renormalization. By implementing one-loop
matching of the light-cone wave functions, we determine the Collins-Soper
kernel with transverse separation up to 0.6 fm. We also discuss the systematic
uncertainties from operator mixing and scale dependence. The impact from
higher power corrections is also taken into account. Our results allow for a
more precise first-principle determination of the soft function and other
transverse-momentum dependent quantities from lattice QCD.

##  I Introduction

Understanding the internal, three-dimensional structure of hadrons is a
primary goal of nuclear and particle physics. In this regard, the transverse
momentum-dependent (TMD) parton distribution functions (TMDPDFs)
Collins:1981uk; Collins:1981va play an important role as they characterize the
intrinsic transverse partonic structure of protons. They are also essential
ingredients in the description of multi-scale and non-inclusive processes,
such as Drell-Yan production of electroweak gauge bosons or Higgs bosons or
semi-inclusive deep-inelastic scattering with small transverse momentum, in
the context of QCD factorization theorems. As a result, they have received
considerable attention in the past few decades (for a review, see Ref.
Angeles-Martinez:2015sea). More accurate experimental measurements are
expected in the coming decades from JLab 12 GeV and the Electron-Ion Colliders
in the US Accardi:2012qut; AbdulKhalek:2021gbh and in China Anderle:2021wcy.

In contrast to the TMDPDFs that characterize the probability density of parton
momenta in hadrons, the transverse momentum-dependent wave functions (TMDWFs)
offer a probability amplitude description of the partonic structure of
hadrons, from which one can extract all quark/gluon densities. In the TMD
factorization scheme, they are the most important ingredients to predict
physical observables in exclusive processes, for instance, weak decays of
heavy $B$ meson Keum:2000wi; Lu:2000em which are valuable to extract the CKM
matrix element and probe new physics beyond the standard model. However, due
to the lack of precise knowledge of TMDWFs, the one-dimensional lightcone
distribution amplitudes (LCDAs) are used instead in most analyses of $B$
decays Keum:2000wi; Lu:2000em; Nagashima:2002ia, resulting in uncontrollable
errors. The unprecedented precision of experimental measurements of $B$ decays
Cerri:2018ypt strongly requests a reliable theoretical knowledge of TMDWFs.

A common feature of TMDPDFs and TMDWFs is that they depend both on the
longitudinal momentum fraction $x$ and on the transverse spatial separation of
partons. Considerable theoretical efforts have been devoted in recent years to
determine these quantities by fitting the pertinent experimental data
Landry:1999an; Landry:2002ix; DAlesio:2014mrz; Sun:2014dqm; Konychev:2005iy;
Bacchetta:2017gcc; Scimemi:2017etj; Scimemi:2019cmh; Bacchetta:2019sam, which,
however, is limited by the imprecise knowledge of the nonperturbative
behaviour of TMDPDFs and TMDWFs. Thus, it is highly desirable to find a method
to calculate them from first-principle approaches such as lattice QCD.

This has been realized in the framework of large momentum effective
theoryJi:2013dva; Ji:2014gla, which offers a systematic way to calculate
light-cone correlations by simulating time-independent Euclidean correlations
on the lattice. Significant progress has been made in calculating various
parton quantities from LaMET. For recent reviews, see Ref. Cichy:2018mum;
Ji:2020ect.

A very important property of LaMET is that the TMDPDFs and TMDWFs can be
accessed through the so-called quasi-TMDPDFs and quasi-TMDWFs Ebert:2019okf;
Ji:2019sxk; Ji:2019ewn; Ji:2021znw. In Ref. Ji:2019sxk, it has been noticed
for the first time that a four-quark form factor can be directly calculated on
the lattice, and that TMD factorization ensures that this quantity is
factorized into TMDWFs, soft factor and the matching kernel. The matching
kernel is perturbative, but obtaining the full TMD distribution requires
detailed simulations of both four-quark form factors and quasi TMDWFs on the
lattice Ji:2019sxk; Ji:2019ewn. On the other hand, one can make use of the
multiplicative factorization to obtain the CS kernel from quasi TMDWFs. The
first results for the CS kernel based on these proposals have been published
recently Shanahan:2020zxr; LatticeParton:2020uhz; Schlemmer:2021aij;
Li:2021wvl; Shanahan:2021tst.

In this work, we present a state-of-the-art calculation of the CS kernel,
based on a lattice QCD analysis of quasi-TMDWFs with $N_{f}=2+1+1$ valence
clover fermions on a staggered quark sea ???. A single ensemble with lattice
spacing $a\simeq 0.12$fm, volume $n_{s}^{3}\times n_{t}=48^{3}\times 64$, and
physical sea-quark masses is used. In order to improve the signal-to-noise
ratio, we tune the light-valence quark masses such that $m_{\pi}=670$ MeV. The
CS kernel is then extracted through the ratios of the quasi-TMDWFs and the
perturbative matching kernels at different momenta,
$P^{z}=2\pi/n_{s}\times\\{8,10,12\\}=\\{1.72,~{}2.15,~{}2.58\\}$GeV. This
corresponds to the boost factor $\gamma=\\{2.57,~{}3.21,~{}3.85\\}$,
respectively. This analysis improves the previous ones LatticeParton:2020uhz;
Li:2021wvl by taking into account the one-loop perturbative contributions, and
analyzing systematic uncertainties from operator mixing, higher-order
corrections from the scale dependence, and higher power corrections in terms
of $1/P^{z}$.

The rest of this paper is organized as follows. In Sec. II, we give the
theoretical framework to extract the CS kernel from quasi-TMDWFs. Numerical
results for quasi-TMDWFs and CS kernel are presented in Sec. III. A brief
summary of this work is given in Sec. IV. More details about the analysis are
collected in the appendix.

##  II Theoretical framework

###  II.1 Collins-Soper Kernel and Rapidity Evolution

Unlike the lightcone PDFs and distribution amplitudes, the TMDPDFs and TMDWFs
depend on both the renormalization scale $\mu$ and an additional rapidity
renormalization scale. The latter arises because the matrix elements also
suffer from so-called rapidity divergences that require a dedicated regulator
Collins:1981uk; Becher:2010tm; Chiu:2011qc. In TMD factorization, the ???
contributions of hard, i.e. highly virtual, modes to the Born process are
usually calculated in the dimensional regularization scheme. ??? Collinear
modes in distinct directions and soft modes share the same virtuality ??? ,
and are only distinguishable by their rapidity. In calculations using
regularisation schemes such as dimensional regularisation, which only regulate
virtualities, one will encounter additional rapidity divergences that arise in
soft and collinear matrix elements when integrating over rapidity, and have to
be resolved using a dedicated regulator. After regularization, TMDPDFs and
TMDWFs acquire an additional rapidity scale dependence. This dependence should
cancel in theoretical predictions for physical observables.

The CS kernel $K\left(b_{\perp},\mu\right)$, known as the rapidity anomalous
dimension, encodes the rapidity dependence of the TMD distributions
Collins:1981va; Collins:1981uk:

$\displaystyle 2\zeta\frac{d}{d\zeta}\ln
f^{\mathrm{TMD}}\left(x,b_{\perp},\mu,\zeta\right)=K\left(b_{\perp},\mu\right),$
(1)

where $f^{\mathrm{TMD}}$ denotes any leading twist ??? TMDPDF or TMDWF. The
TMD distributions depend on the longitudinal momentum fraction $x$, transverse
spatial separation $b_{\perp}$, which is the Fourier-conjugate to the
transverse momentum $k_{\perp}$, as well as the renormalization scale $\mu$
and rapidity scale $\zeta$ which is related to the hadron momentum. The
$\mu$-dependence of CS kernel $K\left(b_{\perp},\mu\right)$ satisfies the
following renormalization group equation (RGE):

$\displaystyle\mu^{2}\frac{d}{d\mu^{2}}K\left(b_{\perp},\mu\right)=-\Gamma_{\text{cusp}}\left(\alpha_{s}\right).$
(2)

Here
$\Gamma_{\text{cusp}}\left(\alpha_{s}\right)=\alpha_{s}C_{F}/\pi+\mathcal{O}(\alpha_{s}^{2})$
is the cusp anomalous dimension, which has been calculated in perturbation
theory up to 2-loop in Ref. Li:2016ctv, and 3-loop in Ref. Moch:2017uml. The
solution to the RGE can be expressed as:

$\displaystyle
K\left(b_{\perp},\mu\right)=-2\int_{1/b_{\perp}}^{\mu}\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\Gamma_{\text{cusp}}\left(\alpha_{s}\left(\mu^{\prime}\right)\right)+K\left(\alpha_{s}\left(1/b_{\perp}\right)\right).$
(3)

For large $b_{\perp}$, where $b_{\perp}^{-1}\lesssim\Lambda_{\mathrm{QCD}}$,
the CS kernel becomes nonperturbative, which is represented by the non-cusp
anomalous dimension $K\left(\alpha_{s}\left(1/b_{\perp}\right)\right)$ in the
above equation.

In the past decades, the CS kernel has been widely studied in global fits of
TMD parton distributions Landry:1999an; Landry:2002ix; DAlesio:2014mrz;
Sun:2014dqm; Konychev:2005iy; Bacchetta:2017gcc; Scimemi:2017etj;
Scimemi:2019cmh; Bacchetta:2019sam. The explicit form in the nonperturbative
region can only be parametrized by extending the perturbative expressions at
small $b_{\perp}$, which inevitably introduces systematic uncertainties. A
direct calculation on the lattice of TMDPDFs and the relevant CS kernel was an
almost unsurmountable hurdle until the establishment of LaMET Ji:2013dva;
Ji:2014gla. A remarkable recent development in LaMET is that these quantities
can be accessed through the corresponding quasi observables Ebert:2019okf;
Ji:2019sxk; Ji:2019ewn; Ji:2021znw.

###  II.2 Quasi TMD Wave Functions

As stated above, one can define the quasi TMDWFs for a highly-boosted pion
along the $z$-direction with large momentum $P^{z}$ as:

$\displaystyle\tilde{\Psi}^{\pm}$
$\displaystyle\left(x,b_{\perp},\mu,\zeta_{z}\right)=\lim_{L\to\infty}\int\frac{dz}{2\pi}e^{ixzP^{z}}\frac{\tilde{\Phi}^{\pm
0}\left(z,b_{\perp},P^{z},a,L\right)}{\sqrt{Z_{E}(2L,b_{\perp},\mu,a)}},$ (4)

where the unsubtracted quasi TMDWF is defined as an equal-time correlator
containing nonlocal quark bilinear operator with staple-shaped gauge link:

$\displaystyle\tilde{\Phi}^{\pm 0}$
$\displaystyle\left(z,b_{\perp},P^{z},a,L\right)=\left\langle
0\right|\bar{\psi}\left(z\hat{n}_{z}/2+b_{\perp}\hat{n}_{\perp}\right)\Gamma$
$\displaystyle\times U_{\sqsupset,\pm
L}\left(z\hat{n}_{z}/2+b_{\perp}\hat{n}_{\perp},-z\hat{n}_{z}/2\right)\psi\left(-z\hat{n}_{z}/2\right)\left|P^{z}\right\rangle.$
(5)

For a pseudoscalar mesonic state, the Dirac structure $\Gamma$ can be chosen
as $\gamma^{z}\gamma_{5}$ or $\gamma^{t}\gamma_{5}$, which approaches the
leading-twist structure $\gamma^{+}\gamma_{5}$ in the light-cone limit. With a
large but finite $P^{z}$, the difference between $\gamma^{z}\gamma_{5}$ and
$\gamma^{t}\gamma_{5}$ is power suppressed by $P^{z}$. Technically, one can
also use a combination of them, such as
$\left(\gamma^{z}\gamma_{5}+\gamma^{t}\gamma_{5}\right)/2$ to diminish power
corrections, and more details can be found in Sec. III.4. Various combinations
were also explored in Ref. Li:2021wvl. The superscript “0” in
$\tilde{\Phi}^{\pm 0}$ indicates bare quantities. The linear divergences come
from the radiation of collinear gluons ??? by the gauge-link, and show up as
singularities in loop momentum integrals in the gauge-link’s self-energy.
Therefore, this divergence appears in the staple-shaped gauge link,

$\displaystyle U_{\sqsupset,\pm
L}\left(z\hat{n}_{z}/2+b_{\perp}\hat{n}_{\perp},-z\hat{n}_{z}/2\right)\equiv$
$\displaystyle
U_{z}^{\dagger}\left(z\hat{n}_{z}/2+b_{\perp}\hat{n}_{\perp};L\right)U_{\perp}\big{(}(L-z/2)\hat{n}_{z};{b}_{T}\big{)}U_{z}\left(-z\hat{n}_{z}/2;L\right),$
(6)

and cannot be regularized by the standard UV regulator in dimensional
regularization. The Euclidean gauge link in $U_{\sqsupset,\pm L}$ is defined
as

$\displaystyle U_{z}(\xi,\pm L)=\mathcal{P}\exp\left[-ig\int_{\xi^{z}}^{\pm
L}d\lambda n_{z}\cdot A\left(\vec{\xi}_{\perp}+n_{z}\lambda\right)\right],$
(7)

where $\xi^{z}=-\xi\cdot n_{z}$. The $\pm L$ corresponds to the farthest
difference that the gauge link can reach in positive or negative $n_{z}$
direction on a finite Euclidean lattice. This is depicted as the blue and red
lines in Fig. 1.

Figure 1: Illustration of the staple-shaped gauge-link included in
unsubtracted quasi-TMDWFs and related Wilson loop. The blue and red double
lines in the upper panel represent the $L$-shift direction on the Euclidean
lattice, and the lower panel shows the correspondingly Wilson loop, which will
subtract UV logarithmic and linear divergences in quasi-TMDWFs .

Since the linear divergence is associated with the gauge-link, it can be
removed by a similar gauge-link with the same total length. An optional choice
is to make use the Wilson loop, denoted as $Z_{E}$. The Wilson loop can be
chosen as the vacuum expectation of a flat rectangular Euclidean Wilson-loop
in the $z$-$\perp$ plane:

$\displaystyle
Z_{E}\left(2L,b_{\perp},\mu,a\right)=\frac{1}{N_{c}}\operatorname{Tr}\left\langle
0\left|U_{\perp}(0;b_{\perp})U_{z}\left(b_{\perp}\hat{n}_{\perp};2L\right)\right|0\right\rangle.$
(8)

Here the length of $Z_{E}$ is twice that of the staple-shaped gauge-link
$U_{z}^{(\dagger)}$ in the $z$ direction, and thus it is anticipated that the
square root of $Z_{E}\left(2L,b_{\perp},\mu\right)$ ??? cancels the linear
divergence and heavy quark potential in the gauge-link. There are residual
logarithmic divergences from the vertex of Wilson line and light quark, which
can be removed by forming an appropriate ratio with a by short distance quasi-
PDF matrix element Ji:2021uvr. As these logarithmic divergences are
independent of $z$, $b_{\perp}$, $P_{z}$ and $L$, they will explicitly cancel
out when the ratio of quasi TMDWFs is studied.

###  II.3 Factorization of Quasi TMDWFs

With the help of soft functions, the infrared contributions in the subtracted
quasi TMDWFs can be properly accounted for such that the same infrared
structures for the quasi TMDWFs and light-cone ones are guaranteed. This
implies a multiplicative factorization theorem in the framework of LaMET
Ebert:2019okf; Ji:2019sxk; Ji:2019ewn; Ji:2021znw:

$\displaystyle\tilde{\Psi}^{\pm}(x,b_{\perp},\mu,\zeta_{z})S_{r}^{1/2}(b_{\perp},\mu)$
$\displaystyle=H^{\pm}\left(\zeta_{z},\overline{\zeta}_{z},\mu^{2}\right)\exp\left[\frac{1}{2}K(b_{\perp},\mu)\ln\frac{\mp\zeta_{z}-i\epsilon}{\zeta}\right]$
$\displaystyle\qquad\times\Psi^{\pm}(x,b_{\perp},\mu,\zeta)+\mathcal{O}\left(\frac{\Lambda_{QCD}^{2}}{\zeta_{z}},M^{2}\zeta_{z},\frac{1}{b_{\perp}^{2}\zeta_{z}}\right),$
(9)

where the superscript $\pm$ in Eq. (9) corresponds to the direction in the
Wilson line, $\Psi^{\pm}$ is the TMDWFs defined in the inﬁnite momentum frame.
The reduced soft function $S^{1/2}_{r}(b_{\perp},\mu)$ emerges from the
different soft gluon radiation effects for $\tilde{\Psi}^{\pm}$ and
$\Psi^{\pm}$, and is defined as a combination of non-cancelling soft functions
??? Ji:2021znw. The mismatch of the rapidity scale $\zeta$ and $\zeta_{z}$ can
be compensated by the CS kernel $K(b_{\perp},\mu)$. Both $S$ and $K$ are
independent of the $\pm$ choice. $H^{\pm}$ is the 1-loop perturbative matching
kernel Ji:2021znw:

$\displaystyle H^{\pm}(\zeta_{z},\overline{\zeta}_{z},\mu)$
$\displaystyle=1+\frac{\alpha_{s}C_{F}}{4\pi}\bigg{(}-\frac{5\pi^{2}}{6}-4+{\ell}_{\pm}+\overline{\ell}_{\pm}-\frac{1}{2}({\ell}_{\pm}^{2}+\overline{\ell}_{\pm}^{2})\bigg{)}.$
(10)

With the abbreviations ${\ell}_{\pm}=\ln\left[(-\zeta_{z}\pm
i\epsilon)/\mu^{2}\right]$, and
$\overline{\ell}_{\pm}=\ln\left[(-\overline{\zeta}_{z}\pm
i\epsilon)/\mu^{2}\right]$, the rapidity scales $\zeta_{z}=(2xP^{z})^{2}$ and
$\overline{\zeta}_{z}=\left(2\bar{x}P^{z}\right)^{2}$, and $\bar{x}=1-x$. It
should be noticed that $H^{\pm}$ contains nonzero imaginary parts in
$\ell_{\pm}$ and $\bar{\ell}_{\pm}$. While the imaginary parts in $\ell_{\pm}$
are constants, the ones in the double logarithms $\bar{\ell}_{\pm}$ are
momentum-dependent.

A characteristic behavior of Eq.(9) is that this factorization is
multiplicative Note: A factor is by definiton multiplicative ???, which
indicates that hard gluon contributions are local. This is due to the fact
that hard gluon exchange between the quark and anti-quark sectors in quasi
TMDWFs is power suppressed: if there were such a hard gluon, the spatial
separation between its attachments is much smaller than $b_{\perp}$, resulting
in power suppression compared to the typical hard mode contributions. Thus, at
leading power the factorization of quasi TMDWFs are multiplicative Note: A
factor is by definiton multiplicative ???. This feature is illustrated in
Fig.2, in which the collinear, soft and hard sub-diagrams represent the
pertinent contributions. Further, $\zeta_{z}$ arising from Lorentz-invariant
combinations of collinear momentum modes, will provide the natural hard scale
of the hard sub-diagram. More detailed explanations for the factorization of
quasi TMDPDFs in LaMET are given in the recent review Ji:2019ewn.

Figure 2: Leading-power reduced graph for pion quasi TMDWFs. Here $C,~{}S$
denote the collinear and soft sectors of the infra-red structure, while $H$
denote the hard contributions. Since the hard-gluon exchange between the quark
and anti-quark is power suppressed, the hard parts are disconnected with each
others, and therefore the factorization of quasi-WF amplitude is
multiplicative.

###  II.4 Collins-Soper Kernel From Quasi TMDWFs

From Eq.(9), one can see that the momentum dependence in quasi TMDWFs provides
an option to determine the CS kernel. This can be written in a way similar to
Eq.(1) Ji:2021znw:

$\displaystyle
2\zeta_{z}\frac{d}{d\zeta_{z}}\ln\tilde{\Psi}^{\pm}\left(x,b_{\perp},\mu,\zeta_{z}\right)=$
$\displaystyle
K\left(b_{\perp},\mu\right)+\frac{1}{2}\mathcal{G}^{\pm}\left(x^{2}\zeta_{z},\mu\right)$
$\displaystyle+$
$\displaystyle\frac{1}{2}\mathcal{G}^{\pm}\left(\bar{x}^{2}\zeta_{z},\mu\right)+\mathcal{O}\left(\frac{1}{\zeta_{z}}\right),$
(11)

where $K\left(b_{\perp},\mu\right)$ denotes the same kernel as in Eq.(1), and
does not depend on the hard scale $\zeta_{z}$ for large $P^{z}$. Unlike
TMDWFs, the quasi distributions also contain hard contributions, whose
rapidity dependence is represented by the perturbative $\mathcal{G}^{\pm}$ as
a function of hard scale $\zeta_{z}$. From the $\zeta_{z}$ dependence of quasi
TMDWFs, we can see that when $P^{z}\to\infty$, the large logarithms in $P^{z}$
are partially absorbed into $K\left(b_{\perp},\mu\right)$ and the remanent is
incorporated in the perturbative matching kernel. Therefore, both the matching
kernel $H^{\pm}$ and an exponential of CS kernel $K\left(b_{\perp},\mu\right)$
are needed to describe the dependence on $\zeta_{z}$ of quasi TMDWFs.

In order to extract the CS kernel $K\left(b_{\perp},\mu\right)$ explicitly,
one can make use of Eq.(9) with two different large momenta $P_{1}^{z}\neq
P_{2}^{z}\gg 1/b_{\perp}$ but the same scale $\zeta_{z}$. Taking a ratio of
these two quantities gives

$\displaystyle\frac{\tilde{\Psi}^{\pm}(x,b_{\perp},\mu,P_{1}^{z})}{\tilde{\Psi}^{\pm}(x,b_{\perp},\mu,P_{2}^{z})}=\frac{H^{\pm}\left(xP_{1}^{z},\mu\right)}{H^{\pm}\left(xP_{2}^{z},\mu\right)}\exp\left[K(b_{\perp},\mu)\ln\frac{P_{1}^{z}}{P_{2}^{z}}\right],$
(12)

where the reduced soft function $S_{r}(b_{\perp},\mu)$ and TMDWFs
$\Psi^{\pm}(x,b_{\perp},\mu,\zeta)$ have been canceled in the ratio.
Therefore, the CS kernel $K\left(b_{\perp},\mu\right)$ can be extracted
through

$\displaystyle
K\left(b_{\perp},\mu\right)=\frac{1}{\ln(P_{1}^{z}/P_{2}^{z})}\ln\frac{H^{\pm}(xP_{2}^{z},\mu)\tilde{\Psi}^{\pm}(x,b_{\perp},\mu,P_{1}^{z})}{H^{\pm}(xP_{1}^{z},\mu)\tilde{\Psi}^{\pm}(x,b_{\perp},\mu,P_{2}^{z})}.$
(13)

Note that the extracted result is formally independent of $x$ and
$P_{1/2}^{z}$ at leading power, and both $\tilde{\Psi}^{+}$ and
$\tilde{\Psi}^{-}$ can be used to extract $K\left(b_{\perp},\mu\right)$. Such
behavior is derived at the leading power in the factorization scheme and might
be undermined by power corrections. Accordingly, in order to diminish the
systematic uncertainties, we take the average:

$\displaystyle K\left(b_{\perp},\mu\right)=$
$\displaystyle\frac{1}{2\ln(P_{1}^{z}/P_{2}^{z})}\left[\ln\frac{H^{+}(xP_{2}^{z},\mu)\tilde{\Psi}^{+}(x,b_{\perp},\mu,P_{1}^{z})}{H^{+}(xP_{1}^{z},\mu)\tilde{\Psi}^{+}(x,b_{\perp},\mu,P_{2}^{z})}\right.$
$\displaystyle+\left.\ln\frac{H^{-}(xP_{2}^{z},\mu)\tilde{\Psi}^{-}(x,b_{\perp},\mu,P_{1}^{z})}{H^{-}(xP_{1}^{z},\mu)\tilde{\Psi}^{-}(x,b_{\perp},\mu,P_{2}^{z})}\right].$
(14)

The details will be discussed in Sec. III.5.

##  III Numerical Simulations and Results

###  III.1 Lattice setup

The numerical simulation ??? use $N_{f}=2+1+1$ valence clover fermions on a
highly improved staggered quark (HISQ) sea ??? Follana:2006rc and a 1-loop
Symanzik improved gauge action Symanzik:1983dc, generated by the MILC
collaboration MILC:2012znn using periodic boundary conditions. In the
calculation, we use a single ensemble with the lattice spacing $a\simeq
0.12$fm and the volume $n_{s}^{3}\times n_{t}=48^{3}\times 64$ at physical
sea-quark masses. In order to increase the signal-to-noise ratio in the
numerical simulation, we tune the light-valence quark masses to the strange
one, namely $m_{\pi}^{\mathrm{sea}}=130$MeV and
$m_{\pi}^{\mathrm{val}}=670$MeV. (The Collins-Soper kernel should only depend
weakly on quark mass.) ??? To further improve the statistical signals, we
adopt hypercubic (HYP) smeared fat links Hasenfratz:2001hp for the gauge
ensembles. To access the infinite momentum limit for the CS kernel, we employ
three different hadron momenta as
$P^{z}=2\pi/n_{s}\times\\{8,10,12\\}=\\{1.72,~{}2.15,~{}2.58\\}$GeV
corresponding to the boost factor $\gamma=\\{2.57,~{}3.21,~{}3.85\\}$.

###  III.2 Quasi TMDWFs From Two Point Correlation Functions

In order to calculate the quasi TMDWFs defined in Eq.(4), we generate Coulomb
gauge fixed wall source propagators

$\displaystyle
S_{w}\left(x,t,t^{\prime};\vec{p}\right)=\sum_{\vec{y}}S\left(t,\vec{x};t^{\prime},\vec{y}\right)e^{i\vec{p}\cdot(\vec{y}-\vec{x})},$
(15)

where $(t^{\prime},\vec{y})$ and $(t,\vec{x})$ denote the space-time positions
of source and sink. Then one can construct the two-point function (2pt)
related to the quasi TMDWFs in Eq.(4):

$\displaystyle C_{2}^{\pm}$
$\displaystyle\left(z,b_{\perp},P^{z};p^{z},L,t\right)=\frac{1}{n_{s}^{3}}\sum_{\vec{x}}\mathrm{tr}e^{i\vec{P}\cdot\vec{x}}\left\langle
S_{w}^{\dagger}\left(\vec{x}_{1},t,0;-\vec{p}\right)\right.$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\times\left.\Gamma
U_{\sqsupset,\pm
L}\left(\vec{x}_{1},\vec{x}_{2}\right)S_{w}\left(\vec{x}_{2},t,0;\vec{p}\right)\right\rangle$
(16)

with $\vec{x}_{1}=\vec{x}+z\hat{n}_{z}/2+b_{\perp}\hat{n}_{\perp}$ and
$\vec{x}_{2}=\vec{x}-z\hat{n}_{z}/2$. The quark momentum
$\vec{p}=(\vec{0}_{\perp},p^{z})$ is along the $z$-direction, and each of two
quarks carries half of the hadron momentum. Thereby the hadron momentum
satisfies $\vec{P}=2\vec{p}$. The anti-quark propagator can be obtained from
Eq.(10) by applying $\gamma_{5}$-hermiticity
$S_{w}^{\dagger}(x,y)=\gamma_{5}S_{w}(y,x)\gamma_{5}$. As mentioned above, the
Dirac structures are chosen as $\Gamma=\gamma^{z}\gamma_{5}$ and
$\gamma^{t}\gamma_{5}$ that can be projected onto the leading twist light-cone
contributions in the large $P^{z}$ limit.

By generating the wall source propagators with quark momentum
$p^{z}=\pm\\{4,5,6\\}\times 2\pi/(n_{s}a)$, and three segments of gauge links
following Eq.(6), one can construct the 2pt function on the lattice. With the
help of reduction formulas, the $C_{2}^{\pm}$ can be parametrized as

$\displaystyle C_{2}^{\pm}$
$\displaystyle\left(z,b_{\perp},P^{z};p^{z},L,t\right)=\frac{A_{w}(p^{z})A_{p}}{2E}\tilde{\Phi}^{\pm
0}\left(z,b_{\perp},P^{z},L\right)$ $\displaystyle\times
e^{-Et}\left[1+c_{0}\left(z,b_{\perp},P^{z},L\right)e^{-\Delta Et}\right],$
(17)

where $A_{w}(p^{z})$ is the matrix element of the pion interpolating ﬁeld with
Coulomb gauge fixed wall source and $A_{p}$ is the one for a point source
(sink). These terms as well as the factor
$E^{-1}=1/\sqrt{m_{\pi}^{2}+\left(P^{z}\right)^{2}}$ are cancelled by the
local two-point function $C_{2}^{\pm}\left(0,0,P^{z};p^{z},0,t\right)$ at the
same time slice. Thus the remaining ground-state matrix element
$\tilde{\Phi}^{\pm 0}$ is normalized. The ratio of nonlocal and local 2pt
functions can be parametrized as

$\displaystyle
R^{\pm}\left(z,b_{\perp},P^{z},L,t\right)=\frac{C_{2}^{\pm}\left(z,b_{\perp},P^{z},L,t\right)}{C_{2}\left(0,0,P^{z},0,t\right)}$
$\displaystyle=$ $\displaystyle\tilde{\Phi}^{\pm
0}\left(z,b_{\perp},P^{z},L\right)\left[1+c_{0}\left(z,b_{\perp},P^{z},L\right)e^{-\Delta
Et}\right].$ (18)

Figure 3: Comparison of two-state fit and one-state fit to extract
$\tilde{\Phi}^{\pm 0}\left(z,b_{\perp},P^{z},L\right)$. Taking
$\\{z,b_{\perp},P^{z},L\\}=\\{0a,2a,24\pi/n_{s},6a\\}$ as example, we can see
that the two-state fit works for $t\in[2a,8a]$ while the one-state fit works
for $t\in[5a,8a]$. The fitted results are consistent with each other, while
the one-state fit is more conservative.

In the above parametrization, the excited-state contributions are collected
into the $c_{0}$ term, and $\Delta E$ denotes the mass gap between the ground
and first excited state. With increasing Euclidean time, contributions from
the excited state decay and the plateau obtained for
$R^{\pm}\left(z,b_{\perp},P^{z},L,t\right)$ at large times reflects the
ground-state contribution $\tilde{\Phi}^{\pm 0}$. We employ two methods to
extract $\tilde{\Phi}^{\pm 0}$, namely the two-state fit directly using Eq.
(18), and the one-state fit setting $c_{0}=0$. Fig. 3 compares the results for
both methods. From this figure, one can see that these results are consistent
within error. While the two-state fit works for $t\in[2a,8a]$, the one-state
fit works for $t\in[5a,8a]$. More results are given in Appendix A. In most
cases, we find that the one-state fit results are consistent with and more
conservative than the two-state fit ones. In the following analysis we chose
the more conservative results from the one-state fit as the central results.

###  III.3 Wilson Loop

The unsubtracted quasi TMDWF matrix elements $\tilde{\Phi}^{\pm
0}\left(z,b_{\perp},P^{z},a,L\right)$ extracted from the joint fit of the two-
point function contain the factor $e^{-\delta\bar{m}(2L+b_{\perp})}$ from the
linear divergence, the  heavy quark effective potential factor
$e^{-V(b_{\perp})L}$ and logarithmic divergences $Z_{O}$:

$\displaystyle\tilde{\Phi}^{\pm 0}\left(z,b_{\perp},P^{z},a,L\right)\propto
e^{-\delta\bar{m}(2L+b_{\perp})}e^{-V(b_{\perp})L}Z_{O}.$ (19)

The linear divergence in $e^{-\delta\bar{m}(2L+b_{\perp})}$ comes from the
self-energy of the Wilson line Ji:2017oey; Ishikawa:2017faj; Green:2017xeu;
Ji:2020brr, where $\delta\bar{m}$ contains a term proportional to $1/a$ and a
non-perturbative renormalon contribution $m_{0}$:

$\displaystyle\delta\bar{m}=\frac{m_{-1}(a)}{a}-m_{0}$ (20)

Note that the exponent of the linear divergence term is proportional to the
total length of the Wilson link, e.g. $2L+b_{\perp}$ for the staple link. Due
to this factor, the numerical value for a Wilson loop dramatically decreases
for small $a$ and large $L$.

The heavy quark effective potential term $e^{-V(b_{\perp})L}$ comes from
interactions between the two Wilson lines along the $z$ direction in the
staple link. The heavy quark effective potential $V(b_{\perp})$ is often used
to determine the lattice spacing of an ensemble.

The logarithmic divergence $Z_{O}$ comes from the vertices involving the
Wilson line and light quark. The logarithmic divergence up to leading order,
resumed by renormalization group equation, and matched to lattice is
Ji:1991pr; LatticePartonCollaborationLPC:2021xdx

$\displaystyle Z_{O}(1/a,\mu)=\left(\frac{\ln[1/(a\Lambda_{\rm
QCD})]}{\ln[\mu/\Lambda_{\rm QCD}]}\right)^{\frac{3C_{F}}{b_{0}}},$

where $\Lambda_{\rm QCD}$ is different from that of $\overline{\rm MS}$ scheme
and lattice perturbation theory. One can use it to effectively absorb higher
order contributions Lepage:1992xa.

In this work, the Wilson loop renormalization method Chen:2016fxx;
Zhang:2017bzy; Musch:2010ka; Green:2017xeu; Zhang:2017zfe is adopted, in which
the Wilson loop $Z_{E}$ defined in Eq.(8) contains linear divergence, heavy
quark effective potential and logarithimic divergence:

$\displaystyle Z_{E}\left(2L,b_{\perp},a\right)\propto
e^{-\delta\bar{m}(4L+2b_{\perp})}e^{-V(b_{\perp})2L}.$ (22)

According to Ref. LatticePartonCollaborationLPC:2021xdx, $\delta\bar{m}$ in
the Wilson loop is the same as that in hadron matrix elements, and thus it is
anticipated that the linear divergence is removed when dividing by
$\sqrt{Z_{E}\left(2L,b_{\perp},a\right)}$:

$\displaystyle\tilde{\Phi}^{\pm}(z,b_{\perp},P^{z},L)=\frac{\tilde{\Phi}^{\pm
0}\left(z,b_{\perp},P^{z},a,L\right)}{\sqrt{Z_{E}\left(2L,b_{\perp},a\right)}}.$
(23)

As shown in Fig. 4, the subtracted quasi TMDWFs tend to be a constant when
$L\geq 0.4$ fm. We then use the subtracted quasi TMDWFs defined as

$\displaystyle\tilde{\Phi}^{\pm}(z,b_{\perp},P^{z})=\lim_{L\to\infty}\frac{\tilde{\Phi}^{\pm
0}\left(z,b_{\perp},P^{z},a,L\right)}{\sqrt{Z_{E}\left(2L,b_{\perp},a\right)}}.$
(24)

However, it is anticipated that there is a residual logarithmic divergence
$Z_{O}$:

$\displaystyle\frac{\tilde{\Phi}^{\pm
0}\left(z,b_{\perp},P^{z},a,L\right)}{\sqrt{Z_{E}\left(2L,b_{\perp},a\right)}}\propto\frac{e^{-\delta\bar{m}(2L+b_{\perp})}e^{-V(b_{\perp})L}Z_{O}}{e^{-\delta\bar{m}(2L+b_{\perp})}e^{-V(b_{\perp})L}}=Z_{O}.$
(25)

In the extraction of the CS kernel, a ratio of quasi TMDWFs is adopted and
accordingly the residual logarithmic divergence $Z_{O}$ is canceled.

Figure 4: Results for the $L$-dependence of unsubtracted and subtracted quasi
TMDWFs: real part (upper panel) and imaginary part (lower panel), as well as
the square root of the Wilson loop. The case $\Gamma=\gamma^{z}\gamma_{5}$ and
$\\{P^{z},b_{\perp},z\\}=\\{16\pi/n_{s},2a,2a\\}$ is used for illustration.

###  III.4 Operator Mixing Effects

Figure 5: $\lambda$-dependence of quasi-WF matrix elements with different
Dirac structures. Here we take the case of
$\\{P^{z},b_{\perp}\\}=\\{24\pi/n_{s},3a\\}$ as an example, the deviation
between these two cases mainly comes from power corrections. Results for more
sets of $\\{P^{z},b_{\perp}\\}$ are shown in appendix C.

For a pseudoscalar meson on a Euclidean lattice, both
$\Gamma=\gamma^{t}\gamma_{5}$ and $\gamma^{z}\gamma_{5}$ project onto the
leading twist light-cone distribution amplitude, i.e. $\gamma^{+}\gamma_{5}$
in the large $P^{z}$ limit. The differences between them arises from power
corrections in terms of $M^{2}/\left(P^{z}\right)^{2}$.

Fig. 5 shows the comparison of the $\lambda=zP^{z}$-dependence of quasi TMDWFs
with $\Gamma=\gamma^{t}\gamma_{5}$ and $\gamma^{z}\gamma_{5}$ at
$P^{z}=24\pi/n_{s}\simeq 2.58$GeV. It can be seen from the plots that there
are some differences between the two sets of results for the real part in the
small $\lambda$ region. The differences are expected to decrease with
increasing $P^{z}$, and the correlators with $\gamma^{t}\gamma_{5}$ and
$\gamma^{z}\gamma_{5}$ will gradually converge to the light-cone from opposite
directions. Besides, in light-cone coordinate, $\gamma^{t}$ and $\gamma^{z}$
can be represented by $\gamma^{+}$ and $\gamma^{-}$,

$\displaystyle\gamma^{t}\gamma_{5}=\frac{1}{\sqrt{2}}(\gamma^{+}+\gamma^{-})\gamma_{5},$
$\displaystyle\gamma^{z}\gamma_{5}=\frac{1}{\sqrt{2}}(\gamma^{+}-\gamma^{-})\gamma_{5},$
(26)

Due to the momentum along the light-cone, operators with
$\gamma^{-}\gamma_{5}$ correspond to higher order terms of TMDWFs. Therefore,
power corrections arising from finite $P^{z}$ are likely to be eliminated in
the average of these two terms:

$\displaystyle\tilde{\Phi}^{\pm}=\frac{1}{2}\left[\tilde{\Phi}^{\pm}\left(\Gamma=\gamma^{t}\gamma_{5}\right)+\tilde{\Phi}^{\pm}\left(\Gamma=\gamma^{z}\gamma_{5}\right)\right].$
(27)

For a quantitative analysis see the appendix of LatticeParton:2020uhz. The
operator mixing effect reaches order 5%, which is much smaller than the
systematic uncertainties discussed in next subsection.

According to our numerical simulations, the subtracted quasi TMDWFs in
coordinate space $\tilde{\Psi}^{\pm}(x,b_{\perp},P^{z})$ as a function of
$\lambda=zP^{z}$ are complex, which is shown in Fig. 6. The examples are the
real and imaginary part of $\tilde{\Psi}^{\pm}(x,b_{\perp},P^{z})$ with
$P^{z}=24\pi/n_{s}$, $b_{\perp}=2a$ and $4a$. To determine quasi TMDWFs in
momentum space $\tilde{\Psi}^{\pm}(x,b_{\perp},P^{z})$, we use a brute-force
Fourier transformation (FT)

$\displaystyle\tilde{\Psi}^{\pm}(x,b_{\perp},P^{z})=\frac{1}{2\pi}\sum_{z_{min}}^{z_{max}}e^{ixzP^{z}}\tilde{\Phi}^{\pm}(z,b_{\perp},P^{z}).$
(28)

Due to the imaginary ??? part of $\tilde{\Phi}^{\pm}(z,b_{\perp},P^{z})$, also
$\tilde{\Psi}^{\pm}(x,b_{\perp},P^{z})$ has an ??? imaginary part. We obtain
the quasi TMDWFs in momentum space for both real and imaginary part of
$\tilde{\Psi}^{\pm}(x,b_{\perp},P^{z})$ shown in Fig.7, by taking
$P^{z}=24\pi/n_{s}$ and $b_{\perp}=2a$ and $4a$ as examples. ??? We truncate
the FT at $z_{min}$ and $z_{max}$. The deviation of
$\tilde{\Phi}^{\pm}(z_{min},b_{\perp},P^{z})$ and
$\tilde{\Phi}^{\pm}(z_{max},b_{\perp},P^{z})$ from zero is a measure of the
resulting truncation error. For the largest range of $z$ values we could
realize numerical, $z_{min}=-1.44$fm, $z_{max}=1.44$fm, this error is ???
still noticable. This brute-force truncation of the FT leads to an oscillatory
behavior of TMDWFs. ???

Figure 6: Examples for subtracted quasi TMDWFs in coordinate space. Here we
take the cases of $P^{z}=24\pi/n_{s}$. According to our results, the
subtracted quasi TMDWFs $\tilde{\Psi}^{+}(z,b_{\perp},P^{z})$ are complex,
thus the real and the imaginary part both need to be investigated. The upper
four figures are for subtracted quasi TMDWFs
$\tilde{\Psi}^{+}(z,b_{\perp},P^{z})$ as a function of $\lambda=zP^{z}$ after
the average over the two Dirac structures ($\gamma^{t}\gamma_{5}$ and
$\gamma^{z}\gamma_{5}$) is taken. The upper left figure shows the real part of
$\tilde{\Psi}^{+}(z,b_{\perp},P^{z})$ with $b_{\perp}=2a$, while the upper
right one is for $\tilde{\Psi}^{+}(z,b_{\perp},P^{z})$ with $b_{\perp}=4a$.
The lower two figures show the corresponding imaginary parts with
$b_{\perp}=2a$ and $4a$.

Figure 7: Examples for subtracted quasi TMDWFs in momentum space with hadron
momentum $P^{z}=24\pi/n_{s}$. The subtracted quasi TMDWFs in momentum space
$\tilde{\Psi}^{+}(x,b_{\perp},P^{z})$ are Fourier transformed from
$\tilde{\Phi}^{+}(z,b_{\perp},P^{z})$, which have real and imaginary part.
Both have to be investigated. Shown in Eq. (28), a brute-force FT is used to
determine $\tilde{\Psi}^{+}(x,b_{\perp},P^{z})$, where $z_{min}=-1.44$fm,
$z_{max}=1.44$fm. The upper left figure shows the real part for
$\tilde{\Psi}^{+}(x,b_{\perp},P^{z})$ with $b_{\perp}=2a$ case, while the
upper right one is for $\tilde{\Psi}^{+}(x,b_{\perp},P^{z})$ with
$b_{\perp}=4a$ case. Correspondingly, the lower two figures shows the
imaginary part for $\tilde{\Psi}^{+}(x,b_{\perp},P^{z})$ with $b_{\perp}=2a$
and $4a$.

###  III.5 Collins-Soper Kernel From Quasi TMDWFs

The CS kernel governs the rapidity evolution and thus is independent of the
momentum fraction of the involved parton. But as indicated in Eq.(9), the
factorization formula works only when $xP^{z}\gg\Lambda_{\rm QCD}$, and could
be invalid in the end-point regions $x\to 0,1$. Power corrections are likely
of the form $1/\left(xP^{z}\right)^{2}$ or $1/\left(\bar{x}P^{z}\right)^{2}$.
Therefore, the numerical CS kernel is fitted by a function of $x$, $P_{1}^{z}$
and $P_{2}^{z}$ and is written as $K(b_{\perp},\mu,x,P^{z}_{1},P^{z}_{2})$,

$\displaystyle K(b_{\perp},\mu,x,P^{z}_{1},P^{z}_{2})=$
$\displaystyle\frac{1}{2\ln(P_{1}^{z}/P_{2}^{z})}\left[\ln\frac{H^{+}(xP_{2}^{z},\mu)\tilde{\Psi}^{+}(x,b_{\perp},\mu,P_{1}^{z})}{H^{+}(xP_{1}^{z},\mu)\tilde{\Psi}^{+}(x,b_{\perp},\mu,P_{2}^{z})}\right.$
$\displaystyle+\left.\ln\frac{H^{-}(xP_{2}^{z},\mu)\tilde{\Psi}^{-}(x,b_{\perp},\mu,P_{1}^{z})}{H^{-}(xP_{1}^{z},\mu)\tilde{\Psi}^{-}(x,b_{\perp},\mu,P_{2}^{z})}\right].$
(29)

Here $K(b_{\perp},\mu,x,P^{z}_{1},P^{z}_{2})$ are extracted from the
perturbative matching kernels and quasi TMDWFs using 1-loop matching. They
will have power corrections of teh form
$\mathcal{O}\left(1/\left(xP^{z}\right)^{2}\right)$ and
$\mathcal{O}\left(1/\left(\bar{x}P^{z}\right)^{2}\right)$. In order to extract
the leading power contributions, we adopt the following parametrization

$\displaystyle K(b_{\perp},\mu,x,P^{z}_{1},P^{z}_{2})=K(b_{\perp},\mu)$
$\displaystyle\quad+A\left[\frac{1}{x^{2}(1-x)^{2}(P^{z}_{1})^{2}}-\frac{1}{x^{2}(1-x)^{2}(P^{z}_{2})^{2}}\right],$
(30)

where $A$ is the coefficient accounting for the magnitude of higher power
contributions, and can be determined through a joint fit of different lattice
data.

Fig. 8 presents the physical CS kernel $K(b_{\perp},\mu)$ at different
$b_{\perp}$. By employing three cases of quasi TMDWFs with
$P^{z}=\\{8,10,12\\}\times 2\pi/n_{s}$, one can extract
$K(b_{\perp},\mu,x,P^{z}_{1},P^{z}_{2})$ with $P_{1}^{z}/P_{2}^{z}=10/8$ and
$12/8$, shown as the different colored bands. Except in the end-point regions
($x<0.2$ or $x>0.8$), the lattice data is flat and reflects the leading power
contribution, which conforms with  expectations. Using the parametrization
formula Eq.(30), the physical CS kernel $K(b_{\perp},\mu)$ can be determined
by fitting the data, shown as the green band. As mentioned before, at large
$b_{\perp}$, the quasi TMDWFs show oscillations due to the truncation of the
Fourier transformation, which also affect the extracted
$K(b_{\perp},\mu,x,P^{z}_{1},P^{z}_{2})$, as shown in the lower panel of
Fig.8.

Figure 8: The fit results of $K(b_{\perp},\mu,x,P^{z}_{1},P^{z}_{2})$
extracted from the quasi-WFs $\tilde{\Psi}^{\pm}$. The chosen momentum pairs
$\\{P_{1}^{z},P^{z}_{2}\\}$ are denoted by $P^{z}_{1}/P_{2}^{z}$ in the
legend. The upper figure is for $b_{\perp}=2a$ and the lower one is for
$b_{\perp}=4a$. The horizontal shaded band shows the central value and
uncertainty of $K(b_{\perp},\mu)$, as well as the fit range of $x$, as
described in the text. The strong oscillation in the shaded area at both
edges, where LaMET approach is invalid, is caused by power corrections.

Theoretically the physical CS kernel is purely real, however, there still
exists a residual imaginary part at 1-loop matching. As discussed above, this
imaginary term comes from the perturbative matching kernel. It is easily to
prove that $H^{\pm}(zP_{2}^{z},\mu)/H^{\pm}(zP_{1}^{z},\mu)$ contains
imaginary part while the lattice calculated
$\tilde{\Psi}^{\pm}(x,b_{\perp},\mu,P_{1}^{z})/\tilde{\Psi}^{\pm}(x,b_{\perp},\mu,P_{2}^{z})$
is nearly real, shown as 9. Therefore, we consider this imaginary part as
systematic uncertainty from factorization theorem. This uncertainty can be
expressed as

$\displaystyle\sigma_{\mathrm{sys}}=\sqrt{K(b_{\perp},\mu)+\text{Im}^{2}\left[K^{+}(b_{\perp},\mu)\right]}-K(b_{\perp},\mu)$
(31)

where $\text{Im}\left[K^{+}(b_{\perp},\mu)\right]$ represents the numerical
imaginary part of extracting $K(b_{\perp},\mu)$ only by $\tilde{\Psi}^{+}$:

$\displaystyle K^{+}(b_{\perp},\mu)=$
$\displaystyle\frac{1}{\ln(P^{z}_{1}/P^{z}_{2})}\ln\frac{H^{+}(xP^{z}_{2},\mu)\tilde{\Psi}^{+}(x,b_{\perp},\mu,P^{z}_{1})}{H^{+}(xP^{z}_{1},\mu)\tilde{\Psi}^{+}(x,b_{\perp},\mu,P^{z}_{2})}.$
(32)

It should be noticed that the perturbative matching kernel $H^{+}(zP^{z},\mu)$
is the complex conjugate of $H^{-}(zP^{z},\mu)$, that is the imaginary parts
in these terms can be cancelled each other when we employ the average of
$H^{+}(zP_{2}^{z},\mu)/H^{+}(zP_{1}^{z},\mu)$ and
$H^{-}(zP_{2}^{z},\mu)/H^{-}(zP_{1}^{z},\mu)$. Therefore, as the final result,
we adopt
$K(b_{\perp},\mu)=\left[K^{+}(b_{\perp},\mu)+K^{-}(b_{\perp},\mu)\right]/2$ to
reserve the real part, and regard the imaginary contributions as our
systematic uncertainty.

Figure 9: Example of numerical results for the imaginary part of the ratio for
quasi TMDWFs at different $P^{z}$,
$\frac{\tilde{\Psi}^{+}(x,b_{\perp},\mu,P_{1}^{z})}{\tilde{\Psi}^{+}(x,b_{\perp},\mu,P_{2}^{z})}$,
with $b_{\perp}=3$a as a function of momentum fraction $x$. The imaginary
parts of both cases are close to zero.

###  III.6 Discussions

One should notice that the Wilson loop renormalized quasi TMDWF on the lattice
(Eq.(24)) has no dependence on scale $\mu$. A natural choice for the scale is
$1/a$, but it does not have the same meaning as the scale $\mu$ in the
$\overline{\rm MS}$ scheme. If one converts it to the $\overline{\rm MS}$
scheme through dividing it by $Z_{O}$ (Eq.(III.3)), the scale $\mu$ is
introduced. In principle, one should convert the Wilson loop renormalized
quasi TMDWF to the $\overline{\rm MS}$ scheme since our factorization formula
works there. However, since $Z_{O}$ has no dependence on momentum $P_{z}$, it
cancels in the ratio of quasi TMDWFs, so does the scale dependence. So, one
doesn’t need to do the scheme and scale conversion of the quasi TMDWF during
the extraction of CS kernel.

The extracted CS kernel from the combined fit of the ratios of quasi TMDWFs
with different momenta are shown as the red data points in Fig. 10, and the
two kinds of errors ???it is unclear which error bar corresponds to which
error??? for each result are the statistical and systematical uncertainties.
In the small-$b_{\perp}$ region, systematical uncertainties are dominant due
to the large power and the nonzero imaginary part.

As a comparison, we also give the tree-level matching result for the CS
kernel. With the leading order matching kernel
$H(xP^{z},\mu)=1+\mathcal{O}(\alpha_{s})$, Eq. (13) simplifies to the ratio of
quasi TMDWFs $\tilde{\Phi}$ at $z=0$ with momentum $P_{1}^{z}/P^{z}_{2}$. Blue
dots in Fig. 10 denote the results obtained for tree-level matching, for which
only statistical uncertainties are shown.

Figure 10: The upper panel shows the comparison of our results
$K(b_{\perp},\mu)$ and $K_{0}(b_{\perp},\mu)$ with the lattice calculations by
SWZ Shanahan:2021tst, LPC LatticeParton:2020uhz, ETMC/PKU Li:2021wvl and
SVZESSchlemmer:2021aij, as well as the perturbative calculations up to 3-loop.
$K(b_{\perp},\mu)$ denotes the CS kernel extracted through 1-loop matching,
and whose uncertainties correspond to the statistical errors and the
systematic ones from the non-zero imaginary part. $K_{0}(b_{\perp},\mu)$
denotes our tree-level results, only with statistical uncertainties. The lower
panel shows the comparison of our result with phenomenological extractions:
SV19 Scimemi:2019cmh, Pavia19 Bacchetta:2019sam and SIYY15 Sun:2014dqm give
phenomenological parameterizations of CS kernel fitted to data from high
energy collision processes like Drell-Yan.

Besides, we also show the CS kernel from perturbative calculations,
phenomenological extractions as well as the lattice results determined by
other collaborations. The black solid and dashed lines in the upper panel of
Fig. 10 indicate the perturbative results up to 3-loops, with coupling
constant $\alpha_{s}(\mu=1/b_{\perp})$. Similar as in this work, the results
of LPC LatticeParton:2020uhz and ETMC/PKU Li:2021wvl are also extracted from
quasi TMDWFs amplitudes through tree-level matching. While the SWZ
Shanahan:2021tst and SVZESSchlemmer:2021aij results are obtained using  the
factorization formula of quasi-TMDPDF. In the small $b_{\perp}$ region, our
results are consistent with the perturbative ones, which adds confidence in
this first-principle calculation of the CS kernel. In the large $b_{\perp}$
region, our results are consistent with other lattice calculations within
uncertainties. The lower panel shows the comparison with phenomenological
results. SV19 Scimemi:2019cmh and SIYY15 Sun:2014dqm use a parameterization
with perturbative and nonperturbative parts, however Pavia19 Bacchetta:2019sam
obtained their result with the factorization of TMDPDFs, obtaining the CS
kernel from the rapidity derivative. ??? In addition they fit parameters from
the Drell-Yan data to obtain their phenomenological CS kernel ???.

##  IV Summary and Outlook

In this work, we have calculated the CS kernel in the large momentum effective
theory framework. In particular the one-loop matching kernel has been adopted,
and several momenta were used to extract the CS kernel. We found that in the
small $b_{\perp}$ region, our results are consistent with the perturbative
ones, which adds confidence in this first-principle calculation of the CS
kernel In large $b_{\perp}$ region, our results are consistent with other
lattice calculations within uncertainties.

Our work shows that the methodology proposed in Ref. Ji:2019sxk; Ji:2019ewn is
indeed valid for the study of the CS kernel from first principles and that
LaMET provides important insights into transverse momentum-dependent hadron
structure.

## Acknowledgement

We thank Xu Feng, Yizhuang Liu, and Feng Yuan for useful discussions. This
work is supported in part by Natural Science Foundation of China under grant
No. 11735010, 11911530088, U2032102, 11653003, 11975127, 11975051, 12005130,
12147140. MC, JH, WW is also supported by Natural Science Foundation of
Shanghai under grant No. 15DZ2272100. PS is also supported by Jiangsu
Specially Appointed Professor Program. YBY is also supported by the Strategic
Priority Research Program of Chinese Academy of Sciences, Grant No.
XDB34030303. AS, PS, WW, YBY and JHZ are also supported by the NSFC-DFG joint
grant under grant No. 12061131006 and SCHA 458/22.

##  Appendix A Euclidean Time $t$ Dependence of Normalized $C_{2}$

In Sec. III.2 the ratio of nonlocal and local two-point functions is
parametrized in Eq.(18)

$\displaystyle
R^{\pm}\left(z,b_{\perp},P^{z},L,t\right)=\frac{C_{2}^{\pm}\left(z,b_{\perp},P^{z},L,t\right)}{C_{2}\left(0,0,P^{z},0,t\right)}$
$\displaystyle=$ $\displaystyle\tilde{\Phi}^{\pm
0}\left(z,b_{\perp},P^{z},L\right)\left[1+c_{0}\left(z,b_{\perp},P^{z},L\right)e^{-\Delta
Et}\right].$ (33)

From the above equation, one can see that
$R^{\pm}\left(z,b_{\perp},P^{z},L,t\right)$ decays exponentially with $t$. As
discussed in Sec. III.2, the one-state and two-state fits are both used to
extract $\tilde{\Phi}^{\pm 0}\left(z,b_{\perp},P^{z},L\right)$. Fig. 11 shows
four examples of fitted results with different
$\\{z,b_{\perp}\\}=\\{0a,1a\\},\\{0a,3a\\},\\{2a,2a\\},\\{2a,3a\\}$, and same
$\\{P^{z},L\\}=\\{24\pi/n_{s},6a\\}$, demonstrating that the one-state and
two-state fits are consistent with each other, but the one-state fit is more
conservative.

##  Appendix B Gauge Line Length $L$ Dependence of Quasi TMDWFs

In Sec. III.3, the Wilson loop is used to renormalize quasi TMDWFs, which
removes the linear divergence. Similar with the discussion in Sec. III.3, we
give results with the different $\\{P^{z},b_{\perp},z\\}$ in Fig. 12 to show
the Wilson line length $L$-dependence of Wilson loop, unsubtracted quasi
TMDWFs and subtracted quasi TMDWFs. At large $L$, $\tilde{\Phi}^{+0}$ decays
at the same speed of $\sqrt{Z_{E}}$, so Wilson loop cancels the linear
divergence in unsubtracted quasi TMDWFs.

##  Appendix C Operator Mixing Effects

As described in Sec. III.4, quasi TMDWFs for a pseudoscalar meson require the
projectors $\Gamma=\gamma^{t}\gamma_{5}$ or $\Gamma=\gamma^{z}\gamma_{5}$ onto
the leading twist light-cone distribution amplitude, i.e.
$\gamma^{+}\gamma_{5}$ in large $P^{z}$ limit. Fig. 13 shows examples with
different $\\{P^{z},b_{\perp}\\}$ for comparing quasi TMDWFs as functions of
$\lambda=zP^{z}$ of two Dirac matrices $\Gamma=\gamma^{t}\gamma_{5}$ and
$\Gamma=\gamma^{z}\gamma_{5}$. In small $\lambda$ area, the behaviors of quasi
TMDWFs for two Dirac matrices are a little different, which is expected to
decrease with the increase of $P^{z}$. So the average of these two cases is
likely to eliminate the power corrections.

##  Appendix D $x$-dependence in Collins-Soper kernel

As Sec. III.5 describes, there are power corrections to quasi TMDWFs with the
form of $1/(xP^{z})^{2}$ or $1/(\bar{x}P^{z})^{2}$. The $x$-dependence
analysis is needed for extracting CS kernel from quasi TMDWFs. Eq. 30 gives
parametrization form of CS kernel as a function of $x$, which is used for
$x$-dependence fit. Fig. 14 gives fit results with
$b_{\perp}=\\{0.12,0.36,0.6\\}$fm of CS kernel of $x$-dependence. When the
$b_{\perp}$ get larger, the oscillation becomes more severe due to the the
finite truncation in the Fourier transformation as described in Sec. III.4.

Figure 11: Four examples for comparing two-state fit and one-state fit to
extract the $\tilde{\Phi}^{\pm 0}(z,b_{\perp},P^{z},L)$ from
$R^{\pm}(z,b_{\perp},P^{z},L,t)$ as described in Sec. III.2 with
$\\{z,b_{\perp}\\}=\\{0a,1a\\},\\{0a,3a\\},\\{2a,2a\\},\\{2a,3a\\}$, and
$\\{P^{z},L\\}=\\{24\pi/n_{s},6a\\}$. The fit range for two-state fit is
$t\in[2a,8a]$, which for one-state fit is $t\in[5a,8a]$.

Figure 12: Results showing the dependence on the gauge line length $L$ of
unsubtracted and subtracted quasi TMDWFs as well as the Wilson loop with
$\\{P^{z},b_{\perp},z\\}$ shown in each figure. These results are for
$\Gamma=\gamma^{z}\gamma_{5}$.

Figure 13: Examples of comparisons for $\lambda$-dependence of quasi-WF matrix
elements with two Dirac matrices: $\Gamma=\gamma^{t}\gamma_{5}$ and
$\Gamma=\gamma^{z}\gamma_{5}$, with $\\{P^{z},b_{\perp}\\}$ shown in each
figure. Power corrections cause the deviation between both cases.

Figure 14: The fit results of $K(b_{\perp},\mu,x,P^{z}_{1},P^{z}_{2})$
extracted of quasi-WFs $\tilde{\Psi}^{\pm}$ for the cases
$b_{\perp}=\\{0.12,0.36,0.6\\}$fm.

## References

  * (1)
  * (2)