# Lattice QCD calculation of the Collins-Soper kernel from quasi TMDPDFs

Phiala Shanahan<EMAIL_ADDRESS>Center for Theoretical Physics, Massachusetts
Institute of Technology, Cambridge, MA, USA 02139  Michael Wagman
<EMAIL_ADDRESS>Fermi National Accelerator Laboratory, Batavia, IL 60510, USA
Yong Zhao<EMAIL_ADDRESS>Physics Division, Argonne National Laboratory,
Lemont, IL 60439, USA  Physics Department, Brookhaven National Laboratory,
Bldg. 510A, Upton, NY 11973, USA

###### Abstract

This work presents a lattice quantum chromodynamics (QCD) calculation of the
nonperturbative Collins-Soper kernel, which describes the rapidity evolution
of quark transverse-momentum-dependent parton distribution functions. The
kernel is extracted at transverse momentum scales in the range 400 MeV
$<q_{T}<1.7$ GeV in a calculation with dynamical fermions and quark masses
corresponding to a larger-than-physical pion mass, $m_{\pi}=538(1)$ MeV. It is
found that different approaches to extract the Collins-Soper kernel from the
same underlying lattice QCD matrix elements yield significantly different
results and uncertainty estimates, revealing that power corrections, such as
those associated with higher-twist effects, and perturbative matching between
quasi and light-cone beam functions, cannot be neglected.

††preprint: FERMILAB-PUB-21-326-T,MIT-CTP/5316

##  I Introduction

Transverse-momentum-dependent parton distribution functions (TMDPDFs) describe
the intrinsic transverse momentum $q_{T}$ of the partonic constituents of
hadrons Collins and Soper (1981, 1982); Collins  _et al._ (1985). These non-
perturbative functions can be accessed directly in high-energy scattering
processes such as Drell-Yan production and semi-inclusive deep-inelastic
scattering with small transverse hadron momentum $q_{T}$ Scimemi and
Vladimirov (2019); Bacchetta  _et al._ (2019), and indirectly through other
processes such as studies of hadrons in jets Buffing  _et al._ (2018);
Gutierrez-Reyes  _et al._ (2019). Significant efforts are underway to improve
constraints on TMDPDFs both from current and planned experiments Gautheron
_et al._ (2010); Dudek  _et al._ (2012); Aschenauer  _et al._ (2015); Accardi
_et al._ (2016); Abdul Khalek  _et al._ (2021) and through theory calculations
in the framework of lattice quantum chromodynamics (QCD) Musch  _et al._
(2011, 2012); Engelhardt  _et al._ (2016); Yoon  _et al._ (2015, 2017);
Shanahan  _et al._ (2019, 2020); Zhang  _et al._ (2020); Schlemmer  _et al._
(2021); Li  _et al._ (2021) using approaches such as large-momentum effective
theory (LaMET) Ji (2013, 2014); Ji  _et al._ (2020) or the Lorentz-invariant
method based on ratios of TMDPDFs Musch  _et al._ (2011).

TMDPDFs $f_{i}^{\text{TMD}}(x,b_{T},\mu,\zeta)$, defined for a parton of
flavor $i$ in a given hadron state, are functions of the longitudinal momentum
fraction $x$ of the parton, the Fourier conjugate $b_{T}$ of $q_{T}$, the
virtuality scale $\mu$, and the rapidity scale $\zeta$ which is related to the
hadron momentum. While the $\mu$-evolution of TMDPDFs is perturbative for
perturbative scales $\mu$ and $\zeta$, the $\zeta$-evolution is governed by
the Collins-Soper evolution kernel (or anomalous dimension)
$\gamma_{\zeta}^{i}(\mu,b_{T})$, which is nonperturbative for scales
$b_{T}\sim q_{T}^{-1}\sim\Lambda^{-1}_{\mathrm{QCD}}$, even if both $\mu$ and
$\zeta$ are perturbative. Constraints on the kernel
$\gamma_{\zeta}^{i}(\mu,b_{T})$ in the nonperturbative region are necessary in
order to relate TMDPDFs determined from experiment or lattice QCD at different
scales.

Recently, it was shown in Refs. Ebert  _et al._ (2019a, b, 2020) that the
Collins-Soper kernel can be calculated from ratios of quasi TMDPDFs at
different hadron momenta, quantities which are both calculable in lattice QCD
and which can be related to TMDPDFs Ji  _et al._ (2015, 2019a); Ebert  _et
al._ (2019b); Ji  _et al._ (2019b, c). This provides a pathway to first-
principles QCD calculations of the kernel in the nonperturbative region, which
will provide valuable complementary information to constraints from global
analyses of experimental data. This prospect has motivated a series of proof-
of-principle lattice QCD investigations of the Collins-Soper kernel both
directly Shanahan  _et al._ (2020, 2019); Zhang  _et al._ (2020); Schlemmer
_et al._ (2021); Li  _et al._ (2021) through the approach of Refs. Ebert  _et
al._ (2019a, b, 2020) and via related prescriptions Vladimirov and Schäfer
(2020).

In this work, a direct calculation of the Collins-Soper kernel is presented,
based on a lattice QCD study with dynamical fermions and quark masses
corresponding to a larger-than-physical pion mass $m_{\pi}=538(1)$ MeV, and a
single value of the lattice spacing and volume. The kernel is extracted at
transverse momentum scales in the range 400 MeV $<q_{T}<1.7$ GeV and compared
with phenomenological parametrizations and existing lattice QCD calculations.
This analysis includes several advances over previous lattice QCD studies of
the Collins-Soper kernel via the same approach. In particular, matching of
quasi TMDPDFs and TMDPDFs is performed to one-loop order, the mixing of
different TMDPDFs under renormalization is fully accounted for, and the
analysis includes improved treatments of power corrections and systematic
effects arising from the finite lattice volume and various statistical
limitations of the calculation. It is found that different approaches to
extract the Collins-Soper kernel from the same underlying lattice QCD matrix
elements yield significantly different results and uncertainty estimates,
revealing that power corrections, such as those associated with higher-twist
effects, and perturbative matching between quasi and light-cone beam
functions, cannot be neglected.

The method by which the Collins-Soper kernel can be computed following Refs.
Ebert  _et al._ (2019a, b) is detailed in Section II. The lattice QCD
calculation is reported in Section III, while a summary and outlook is
provided in Section IV.

##  II Quasi TMDPDFs and the Collins-Soper kernel

The quark Collins-Soper kernel $\gamma^{q}_{\zeta}(\mu,b_{T})$ can be computed
in lattice QCD from a ratio of nonsinglet quasi TMDPDFs
$\tilde{f}_{{\text{ns}}}^{\mathrm{TMD}}$ at different hadron momenta (taken in
the $z$-direction) $P^{z}_{i}\gg\Lambda_{\text{QCD}}$ Ebert  _et al._ (2019a,
b); Ji  _et al._ (2019c):

$\displaystyle\gamma_{\zeta}^{q}(\mu,b_{T})$
$\displaystyle=\frac{1}{\ln(P^{z}_{1}/P^{z}_{2})}$
$\displaystyle\quad\times\ln\frac{C^{\mathrm{TMD}}_{\text{ns}}(\mu,xP_{2}^{z})\,\tilde{f}_{{\text{ns}}}^{\mathrm{TMD}}(x,\vec{b}_{T},\mu,P_{1}^{z})}{C^{\mathrm{TMD}}_{\text{ns}}(\mu,xP_{1}^{z})\,\tilde{f}_{{\text{ns}}}^{\mathrm{TMD}}(x,\vec{b}_{T},\mu,P_{2}^{z})}$
$\displaystyle\quad+{\cal O}\Big{(}{1\over(xP^{z}b_{T})^{2}}\,,{\Lambda_{\rm
QCD}^{2}\over(xP^{z})^{2}}\Big{)}\,.$ (1)

The perturbative matching coefficient $C^{\mathrm{TMD}}_{\text{ns}}$ relates
the quasi TMDPDFs, which are defined in terms of Euclidean-space matrix
elements as detailed below, to the corresponding light-cone TMDPDFs through a
factorization theorem based on an expansion in powers of the nucleon momentum
Ebert  _et al._ (2019a, b); Ji  _et al._ (2019b, c). Additional
nonperturbative factors related to the soft sector Ji  _et al._ (2019a); Ebert
_et al._ (2019b) cancel in the ratio; recently exploratory lattice QCD studies
of these factors have been performed Zhang  _et al._ (2020); Li  _et al._
(2021) following the approach proposed in Refs. Ji  _et al._ (2019b, c). The
flavor nonsinglet unpolarized quark quasi TMDPDF is defined as
$\tilde{f}^{\mathrm{TMD}}_{\text{ns}}=\tilde{f}^{\mathrm{TMD}}_{u}-\tilde{f}^{\mathrm{TMD}}_{d}$,
where

$\displaystyle\tilde{f}_{i}^{\mathrm{TMD}}\big{(}x,\vec{b}_{T},\mu,P^{z}\big{)}\equiv\lim_{\begin{subarray}{c}a\to
0\\\
\eta\to\infty\end{subarray}}\int\frac{\mathrm{d}b^{z}}{2\pi}e^{-\mathrm{i}b^{z}\left(xP^{z}\right)}\mathcal{Z}^{\overline{\mathrm{MS}}}_{\gamma^{4}\Gamma}(\mu,b^{z}\\!,a)$
$\displaystyle\qquad\qquad\times{P^{z}\over
E_{\vec{P}}}\tilde{B}^{\Gamma}_{i}\big{(}b^{z},\vec{b}_{T},a,\eta,P^{z}\big{)}\tilde{\Delta}_{S}\left(b_{T},a,\eta\right).$
(2)

Here $a$ denotes the lattice spacing, and
$E_{\vec{P}}=\sqrt{\vec{P}^{2}+m_{h}^{2}}$ where $\vec{P}=P^{z}\vec{e}_{z}$ is
the hadron three-momentum and $m_{h}$ is the hadron mass. The factor
$\mathcal{Z}^{\overline{\mathrm{MS}}}_{\gamma^{4}\Gamma}(\mu,b^{z},a)$, where
$\Gamma$ is a Dirac matrix label, renormalizes the quasi TMDPDF and matches it
to the $\overline{\mathrm{MS}}$ scheme at scale $\mu$ Constantinou  _et al._
(2019); Ebert  _et al._ (2020); Shanahan  _et al._ (2019), and the quasi soft
factor $\tilde{\Delta}_{S}$ Ji  _et al._ (2015, 2019a); Ebert  _et al._
(2019a, b) and quasi beam function $\tilde{B}^{\Gamma}_{i}$ are both
calculable in lattice QCD. Summation over $\Gamma$ is implied, accounting for
operator mixing between quasi TMDPDFs with different Dirac structures
resulting from the breaking of rotational and chiral symmetries in lattice QCD
calculations Constantinou  _et al._ (2019); Shanahan  _et al._ (2019); Green
_et al._ (2020); Ji  _et al._ (2021). Mixing with gluon operators is neglected
in Eq. (II), but cancels in the nonsinglet combination of quasi TMDPDFs. It
should be noted that the choice of the Dirac structure $\gamma^{4}$ in Eq.
(II) is not unique; the quasi TMDPDF with Dirac structure $\gamma^{3}$ can
also be boosted onto $\gamma^{+}$ and thus be matched to the spin-independent
TMDPDF in the infinite-momentum limit (in that case, the factor of
$P^{z}/E_{\vec{P}}$ in Eq. (II) is replaced by 1). While the notation is
specialized to $\gamma^{4}$ for clarity throughout this exposition, numerical
results are presented for both choices of Dirac structure in Sec. III.

The quasi beam function $\tilde{B}^{\Gamma}_{i}$ is defined as the matrix
element of a nonlocal quark bilinear operator with a staple-shaped Wilson line
in a boosted hadron state:

$\displaystyle\tilde{B}^{\Gamma}_{i}(b^{z},\vec{b}_{T},a,\eta,P^{z})=$
$\displaystyle\Bigl{\langle}h(P^{z})\big{|}\mathcal{O}_{\Gamma}^{i}(b^{\mu},0,\eta)\big{|}h(P^{z})\Bigr{\rangle}\,,$
(3)

where $h(P^{z})$ denotes the state of hadron $h$ with four-momentum
$P^{\mu}=(0,0,P^{z},E_{\vec{P}})$. The operator
$\mathcal{O}_{\Gamma}^{i}(b^{\mu},0,\eta)$, depicted in Fig. 1, is defined as

$\displaystyle\mathcal{O}^{i}_{\Gamma}(b^{\mu},z^{\mu},\eta)\equiv$
$\displaystyle\
\bar{q}_{i}(z^{\mu}+b^{\mu})\frac{\Gamma}{2}W_{\hat{z}}(z^{\mu}+b^{\mu};\eta-b^{z})$
$\displaystyle\ \times
W^{\dagger}_{T}(z^{\mu}+\eta\hat{z};b_{T})W^{\dagger}_{\hat{z}}(z^{\mu};\eta)q_{i}(z^{\mu})$
$\displaystyle\equiv$ $\displaystyle\
\bar{q}_{i}(z^{\mu}+b^{\mu})\frac{\Gamma}{2}\widetilde{W}(\eta;b^{\mu};z^{\mu})q_{i}(z^{\mu}),$
(4)

where $b^{\mu}=(\vec{b}_{T},b^{z},0)$, and $W_{\hat{\alpha}}(x^{\mu};\eta)$
denotes a Wilson line beginning at $x^{\mu}$ with length $\eta$ in the
direction of ${\hat{\alpha}}$. The subscript $T$ denotes that the associated
Wilson line is in a direction transverse to $\hat{z}^{\mu}=(0,0,1,0)$.

Figure 1: Diagrammatic representation of the Wilson line included in the
operators $\mathcal{O}^{i}_{\Gamma}(b^{\mu},z^{\mu},\eta)$, see Eq. (4).

In practice, it is useful to define a dimensionless ‘bare’ nonsinglet quasi
beam function:

$\displaystyle B^{\text{bare}}_{\Gamma}(b^{z},\vec{b}_{T},a,\eta,P^{z})\equiv$
$\displaystyle\frac{1}{2E_{\vec{P}}}\left(\tilde{B}^{\Gamma}_{u}(b^{z},\vec{b}_{T},a,\eta,P^{z})\right.$
$\displaystyle\left.\hskip
17.07164pt-\tilde{B}^{\Gamma}_{d}(b^{z},\vec{b}_{T},a,\eta,P^{z})\right),$ (5)

as well as a modified $\overline{\mathrm{MS}}$-renormalized quasi beam
function $B^{\overline{\mathrm{MS}}}_{\Gamma}$ Shanahan  _et al._ (2020):

$\displaystyle
B^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,b^{z},\vec{b}_{T},a,\eta,P^{z})$
$\displaystyle\equiv
Z_{\mathcal{O}_{\gamma^{4}\Gamma}}^{\overline{\mathrm{MS}}}(\mu,b^{z},b_{T}^{R},a,\eta)$
$\displaystyle\times\tilde{R}(b_{T},b_{T}^{R},a,\eta)B^{\text{bare}}_{\Gamma}(b^{z},\vec{b}_{T},a,\eta,P^{z}).$
(6)

Compared with the standard $\overline{\mathrm{MS}}$-renormalized quasi beam
function, this definition includes the additional factor $\tilde{R}$,
described further below. The renormalization factor
$Z_{\mathcal{O}_{\gamma^{4}\Gamma}}^{\overline{\mathrm{MS}}}$ is defined as
the product of a regularization-independent momentum subtraction scheme
($\mathrm{RI}^{\prime}\mathrm{/MOM}$) factor
$Z^{\mathrm{RI}^{\prime}\mathrm{/MOM}}_{\mathcal{O}_{\gamma^{4}\Gamma}}$ and a
perturbative matching factor to the $\overline{\text{MS}}$ scheme,
$\mathcal{R}^{\overline{\mathrm{MS}}}_{\mathcal{O}_{\gamma^{4}\Gamma}}$, which
has been calculated at next-to-leading order in continuum perturbation theory
with dimensional regularization ($D=4-2\epsilon$) Constantinou  _et al._
(2019); Ebert  _et al._ (2020):

$\displaystyle
Z_{\mathcal{O}_{\gamma^{4}\Gamma}}^{\overline{\mathrm{MS}}}(\mu,b^{z}\\!,b_{T},a,\eta)=$
$\displaystyle\mathcal{R}^{\overline{\mathrm{MS}}}_{\mathcal{O}_{\gamma^{4}\Gamma}}(\mu,p_{R},b^{z},\vec{b}_{T},\eta)$
$\displaystyle\times
Z^{\mathrm{RI}^{\prime}\mathrm{/MOM}}_{\mathcal{O}_{\gamma^{4}\Gamma}}(p_{R},b^{z}\\!,\vec{b}_{T},a,\eta).$
(7)

In this expression, the dependence on the $\mathrm{RI}^{\prime}\mathrm{/MOM}$
scale $p_{R}$ and on the direction of $\vec{b}_{T}$ cancels between
$Z^{\mathrm{RI}^{\prime}\mathrm{/MOM}}_{\mathcal{O}_{\gamma^{4}\Gamma}}$ and
$\mathcal{R}^{\overline{\mathrm{MS}}}_{\mathcal{O}_{\gamma^{4}\Gamma}}$ at all
orders in perturbation theory (up to discretization artifacts). The quasi beam
function renormalization factor is related to the TMDPDF renormalization
factor $\mathcal{Z}^{\overline{\mathrm{MS}}}_{\gamma^{4}\Gamma}$ in Eq. (II)
as

$\displaystyle\mathcal{Z}_{\gamma^{4}\Gamma}^{\overline{\mathrm{MS}}}(\mu,b^{z}\\!,a)=Z_{\mathcal{O}_{\gamma^{4}\Gamma}}^{\overline{\mathrm{MS}}}(\mu,b^{z}\\!,b_{T},a,\eta)Z_{S}^{\overline{\mathrm{MS}}}(\mu,b_{T},a,\eta),$
(8)

where $Z_{S}^{\overline{\mathrm{MS}}}$ renormalizes the quasi soft factor
$\tilde{\Delta}_{S}$. The $\eta$ and $b_{T}$-dependence on the right-hand side
of Eq. (8) describes linear power divergences proportional to $\eta/a$ and
$b_{T}/a$ which cancel between the two terms.

To ensure that the matching factor
$\mathcal{R}^{\overline{\mathrm{MS}}}_{\mathcal{O}_{\gamma^{4}\Gamma}}$ is in
the perturbative region, both
$Z^{\mathrm{RI}^{\prime}\mathrm{/MOM}}_{\mathcal{O}_{\gamma^{4}\Gamma}}$ and
$\mathcal{R}^{\overline{\mathrm{MS}}}_{\mathcal{O}_{\gamma^{4}\Gamma}}$ should
be computed at a scale $b_{T}^{R}\ll\Lambda_{\rm QCD}^{-1}$. In Eq. (6), this
scale is taken to be distinct from $\vec{b}_{T}$ which is associated with the
staple geometry of the operator defining the bare quasi beam function. As a
result, $Z_{\mathcal{O}_{\gamma^{4}\Gamma}}^{\overline{\mathrm{MS}}}$ cannot
completely cancel the ultraviolet (UV) divergence in the bare quasi beam
function, and remnant linear divergences $\sim|b_{T}-b_{T}^{R}|/a$ appear in
Eq. (6). The factor $\tilde{R}$ is included to cancel such divergences. One
possible choice of $\tilde{R}$ is Shanahan  _et al._ (2019)

$\displaystyle\tilde{R}(b_{T},b_{T}^{R},a,\eta)=\frac{{Z}_{\mathcal{O}_{\gamma^{4}\gamma^{4}}}^{\mathrm{RI}^{\prime}\mathrm{/MOM}}(p_{R}=\tilde{p}_{R},b^{z}=0,\vec{b}_{T},a,\eta)}{{Z}_{\mathcal{O}_{\gamma^{4}\gamma^{4}}}^{\mathrm{RI}^{\prime}\mathrm{/MOM}}(p_{R}=\tilde{p}_{R},b^{z}=0,\vec{b}_{T}^{R},a,\eta)}\,,$
(9)

defined for a fixed choice of $\tilde{p}_{R}$, and of the directions of
$\vec{b}_{T}$ and $\vec{b}^{R}_{T}$. An alternative choice of $\tilde{R}$ is
defined and used in Ref. Ebert  _et al._ (2020).

In terms of the modified $\overline{\mathrm{MS}}$-renormalized quasi beam
functions, the Collins-Soper kernel may be computed as

$\displaystyle\gamma^{q}_{\zeta}(\mu,b_{T})=\frac{1}{\ln(P^{z}_{1}/P^{z}_{2})}\ln\Biggr{[}\frac{C^{\mathrm{TMD}}_{\text{ns}}(\mu,xP_{2}^{z})}{C^{\mathrm{TMD}}_{\text{ns}}(\mu,xP_{1}^{z})}$
$\displaystyle\\!\times\\!\frac{\int\\!\mathrm{d}b^{z}e^{-ib^{z}\\!xP_{1}^{z}}P_{1}^{z}\lim_{\begin{subarray}{c}a\to
0\\\
\eta\to\infty\end{subarray}}B^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,b^{z},\vec{b}_{T},a,\eta,P_{1}^{z})}{\int\\!\mathrm{d}b^{z}e^{-ib^{z}\\!xP_{2}^{z}}\\!P_{2}^{z}\lim_{\begin{subarray}{c}a\to
0\\\
\eta\to\infty\end{subarray}}B^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,b^{z},\vec{b}_{T},a,\eta,P_{2}^{z})}\Biggr{]}$
$\displaystyle\qquad\qquad+{\cal
O}\Big{(}{1\over(xP^{z}b_{T})^{2}}\,,{\Lambda_{\rm
QCD}^{2}\over(xP^{z})^{2}}\Big{)}\,.$ (10)

Since $\tilde{\Delta}_{S}$ and its renormalization factor
$Z_{S}^{\overline{\mathrm{MS}}}$, as well as the factor $\tilde{R}$ included
in the definition of $B^{\overline{\mathrm{MS}}}_{\gamma^{4}}$, are
independent of $b^{z}$, this expression is equivalent to that in Eq. (II). The
specific choice of $\tilde{R}$ (including the choices of $\tilde{p}_{R}$ and
$\vec{b}_{T}$ and $\vec{b}^{R}_{T}$ orientations) does not affect the value of
$\gamma^{q}_{\zeta}$.

##  III Numerical study

The Collins-Soper kernel is computed via Eq. (II) in a lattice QCD calculation
with four dynamical quark flavors, using an ensemble of gauge field
configurations generated by the MILC collaboration with the one-loop Symanzik
improved gauge action and the highly improved staggered quark action. A single
ensemble is studied, with a lattice volume of $L^{3}\times T=48^{3}\times 64$,
a lattice spacing corresponding to $a=0.12$ fm, and sea quark masses tuned to
approximately match the physical quark masses in nature; see Ref. Bazavov  _et
al._ (2013) for further details of the ensemble generation. Calculations are
performed in a partially-quenched mixed-action setup, with the tree-level
$\mathcal{O}(a)$-improved Wilson clover fermion action used for the valence
quarks, with $\kappa=0.1241$ tuned such that the pion mass is $m_{\pi}=538(1)$
MeV. The gauge fields used in the calculation have been subjected to Wilson
flow to flow-time $\mathfrak{t}=1.0$ Lüscher (2010), to enhance the signal-to-
noise ratio in the numerical results111Note that the flowed gauge fields were
also used for constructing $\not{D}$.. The following sections detail the
computation of each element of $\gamma^{q}_{\zeta}$.

###  III.1 Bare quasi beam functions

Quasi beam functions in a pion external state are computed from ratios of
three-point and two-point correlation functions:

$\displaystyle\mathcal{R}_{\Gamma}(t,\tau,b^{\mu},a,\eta,P^{z})$
$\displaystyle\qquad=\frac{C_{\text{3pt}}^{\Gamma,u}(t,\tau,b^{\mu},a,\eta,P^{z}\vec{e}_{z})-C_{\text{3pt}}^{\Gamma,d}(t,\tau,b^{\mu},a,\eta,P^{z}\vec{e}_{z})}{C_{\text{2pt}}(t,P^{z}\vec{e}_{z})}$
$\displaystyle\qquad\xrightarrow{t\gg\tau\gg
0}B^{\text{bare}}_{\Gamma}(b^{z},\vec{b}_{T},a,\eta,P^{z})+\ldots,$ (11)

where

$\displaystyle
C_{\text{3pt}}^{\Gamma,i}(t,\tau,b^{\mu},a,\eta,\vec{P}=P^{z}\vec{e}_{z})$
$\displaystyle\ =\sum_{\vec{x},\vec{z}}e^{i\vec{P}\cdot\vec{x}}\langle
0|\pi_{\vec{P},S}(\vec{x},t)\mathcal{O}^{i}_{\Gamma}(b^{\mu},(\vec{z},\tau),\eta)\pi_{\vec{P},S}^{\dagger}(0)|0\rangle$
$\displaystyle\xrightarrow{t\gg\tau\gg
0}\frac{Z_{\vec{P}}}{4aE^{2}_{\vec{P}}}e^{-E_{\vec{P}}t}\tilde{B}^{\Gamma}_{i}(b^{z},\vec{b}_{T},a,\eta,P^{z})+\ldots$
(12)

and

$\displaystyle C_{\text{2pt}}(t,\vec{P})$
$\displaystyle=\sum_{\vec{x}}e^{i\vec{P}\cdot\vec{x}}\langle
0|\pi_{\vec{P},S}(\vec{x},t)\pi_{\vec{P},S}^{\dagger}(0)|0\rangle$
$\displaystyle\overset{t\gg
0}{\longrightarrow}\frac{Z_{\vec{P}}}{2aE_{\vec{P}}}e^{-E_{\vec{P}}t}+\ldots.$
(13)

In the construction of the correlation functions, momentum-smeared
interpolating operators
$\pi_{\vec{P},S}(\vec{x},t)=\overline{u}_{S(\vec{P}/2)}(\vec{x},t)\gamma_{5}d_{S(\vec{P}/2)}(\vec{x},t)$
are built from quasi local smeared quark fields $q_{S(\vec{P})}(\vec{x},t)$
constructed with 50 steps of iterative Gaussian momentum-smearing Bali  _et
al._ (2016) with smearing radius $\varepsilon=0.2$. $Z_{\vec{P}}$ denotes the
combination of overlap factors for the source and sink interpolating
operators.

Two and three-point functions are constructed for three choices of pion boost
$\vec{P}=P^{z}\vec{e}_{z}$, with $P^{z}=n^{z}2\pi/L$ for
$n_{z}\in\\{3,5,7\\}$. An effective energy function that asymptotes to
$E_{\vec{P}}$ can be defined as

$\displaystyle E^{\text{eff}}_{\vec{P}}(t)$
$\displaystyle=\frac{1}{a}\text{arccosh}{\left({\frac{C_{\text{2pt}}(t+a,\vec{P})+C_{\text{2pt}}(t-a,\vec{P})}{2C_{\text{2pt}}(t,\vec{P})}}\right)}$
$\displaystyle\overset{t\gg 0}{\longrightarrow}E_{\vec{P}}+\ldots,$ (14)

where the ellipses denote exponentially-suppressed corrections from excited
states. This function is shown for the ensemble investigated here in Fig. 2,
including the results of fits to the two-point correlation functions. The fits
are performed as described in Appendix A of Ref. Shanahan  _et al._ (2020),
with the number of states in each fit chosen by maximizing an information
criterion, and different choices of fit range sampled and combined in a
weighted average. The relative deviations in the extracted energies from the
continuum dispersion relation $E_{\vec{P}}=\sqrt{m_{\pi}^{2}+|\vec{P}|^{2}}$
are at most 5% for all momenta studied, increasing with increasing
$|\vec{P}|$, but consistent with the expected size of lattice artifacts.

Figure 2:  Effective energy function defined in Eq. (III.1) for pion states boosted in the $\hat{z}$-direction with $P^{z}=n^{z}2\pi/L$. Shaded bands display the result of single-exponential fits to the two-point correlation functions with the largest two momenta, and two-exponential fits to those with the smallest two momenta, obtained as described in the text.  $n_{z}$ |  $P^{z}$ [GeV] | $\eta/a$ | $n_{\text{src}}$ | $n_{\text{cfg}}$
---|---|---|---|---
3 | 0.65 | {12,14} | 4 | 96
3 | 0.65 | 23 | 16 | 100
5 | 1.1 | {12,14} | 4 | 449
7 | 1.5 | {12,14} | 16 | 596
Table 1:  Quasi beam functions are computed for $n_{\text{src}}$ source
locations on each of $n_{\text{cfg}}$ configurations for each pion boost
$P^{z}$. Matrix elements of operators with staple widths $\vec{b}_{T}$ in the
positive $\hat{x}$ direction with $0\leq|b_{T}|\leq\eta$ and asymmetries
$-\eta\leq b^{z}\leq\eta$ are computed.

The ratio $\mathcal{R}_{\Gamma}$ defined in Eq. (11) is constructed for all
Dirac structures $\Gamma$ and a range of operators with different staple
widths and asymmetries, detailed in Table 1. As indicated, the number of
measurements is varied for the different boost momenta, to partially
compensate for the differences in statistical noise. All ratios are computed
for sink times $t\in\\{7,9,11,13\\}$ and with all operator insertion times
$\tau$ such that $0<\tau<t$. The fitting procedure used to determine
$B^{\text{bare}}_{\Gamma}$ at each set of parameters is precisely the same
procedure as detailed in Appendix A of Ref. Shanahan  _et al._ (2020). Several
examples of the fits in $t$ and $\tau$ used to extract the quasi beam
functions are shown in Appendix A. An example of the resulting bare quasi beam
functions, for particular choices of $b_{T}$ and $\eta$, is shown in Fig. 3.
Additional examples are provided in Appendix B.

Figure 3:  An example of a bare quasi beam function, computed as described in
the text, for $b_{T}/a=1$ and $\eta/a=14$. Further examples are included in
Appendix B.

###  III.2 Non-perturbative renormalization

Computing the modified $\overline{\mathrm{MS}}$-renormalized quasi beam
function $B^{\overline{\mathrm{MS}}}_{\Gamma}$ by Eq. (6) requires, in
addition to the bare quasi beam functions, the
$\mathrm{RI}^{\prime}\mathrm{/MOM}$ renormalization factor
$Z^{\mathrm{RI}^{\prime}\mathrm{/MOM}}_{\mathcal{O}_{\gamma^{4}\Gamma}}$. This
factor enters the calculation of the renormalized quasi beam function both
through the renormalization itself (Eq. (7)) and in the computation of the
factor $\tilde{R}$ (Eq. (9)). The calculation of the nonperturbative
renormalization undertaken here follows closely the presentation of Ref.
Shanahan  _et al._ (2019).

The matrix
$Z^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\Gamma\Gamma^{\prime}}}$
is defined by the renormalization condition

$Z_{q}^{-1}(p_{R})Z^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\Gamma\Gamma^{\prime}}}(p_{R})\Lambda^{\mathcal{O}_{\Gamma^{\prime}}}_{\alpha\beta}(p)\big{|}_{p^{\mu}=p^{\mu}_{R}}=\Lambda_{\alpha\beta}^{\mathcal{O}_{\Gamma};\text{tree}}(p)\,,$
(15)

where dependence on the lattice spacing is left implicit and
$\Lambda^{\mathcal{O}_{\Gamma};(\text{tree})}$ denotes the bare (tree-level)
amputated Green’s function of the operator $\mathcal{O}_{\Gamma}$ defined in
Eq. (4) in an off-shell quark state in the Landau gauge:

$\Lambda^{\mathcal{O}_{\Gamma}}(p)=S^{-1}(p)G^{\mathcal{O}_{\Gamma}}(p)S^{-1}(p)\,.$
(16)

Here $S(p)$ is the quark propagator projected to momentum $p$:

$\displaystyle S_{\alpha\beta}(p,x)$ $\displaystyle=\sum_{y}e^{-ip\cdot
y}\langle q_{\alpha}(x)\bar{q}_{\beta}(y)\rangle,$ (17) $\displaystyle
S_{\alpha\beta}(p)$ $\displaystyle=\frac{1}{V}\sum_{x}e^{ip\cdot
x}S_{\alpha\beta}(p,x),$ (18)

and $G^{\mathcal{O}_{\Gamma}}$ denotes the non-amputated quark-quark Green’s
function with one insertion of $\mathcal{O}_{\Gamma}$, which implicitly
depends on the geometry of the staple-shaped Wilson line defining the
operator:

$\displaystyle G^{\mathcal{O}_{\Gamma}}_{\alpha\beta}(p)$ $\displaystyle=\
\frac{1}{V}\sum_{x,y,z}{\rm e}^{{\mathrm{i}}p\cdot(x-y)}\langle
q_{\alpha}(x)\mathcal{O}_{\Gamma}(z+b,z)\bar{q}_{\beta}(y)\rangle,$
$\displaystyle=\
\frac{1}{V}\sum_{z}\langle\\!\gamma_{5}S^{\dagger}(p,\\!b+z)\gamma_{5}\widetilde{W}(\eta;b+z,\\!z)\frac{\Gamma}{2}S(p,\\!z)\\!\rangle_{\alpha\beta},$
(19)

where $V=L^{3}\times T$ is the lattice volume. The quark wavefunction
renormalization $Z_{q}$ is defined as

$\displaystyle Z_{q}(p_{R})S(p)\big{|}_{p^{2}=p_{R}^{2}}=S^{\text{tree}}(p)$
$\displaystyle\implies$ $\displaystyle Z_{q}(p_{R})=\
\frac{1}{12}\text{Tr}\left[S^{-1}(p)S^{\text{tree}}(p)\right]\bigg{|}_{p^{2}=p_{R}^{2}}$
(20) $\displaystyle\hphantom{Z_{q}(p_{R})}=\ \frac{{\rm Tr}\left[{\rm
i}\sum_{\lambda}\gamma_{\lambda}\sin(ap_{\lambda})S^{-1}(p)\right]}{12\sum_{\lambda}\sin^{2}(ap_{\lambda})}\bigg{|}_{p^{2}=p_{R}^{2}}.$

From Eq. (15), the matrix of $\mathrm{RI}^{\prime}\mathrm{/MOM}$
renormalization factors may be computed as

$\left(Z^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\Gamma\Gamma^{\prime}}}(p_{R})\right)^{-1}=\frac{\mathcal{V}^{\mathcal{O}_{\Gamma\Gamma^{\prime}}}(p)}{6e^{ip_{R}\cdot
b}Z_{q}(p_{R})}\bigg{|}_{p^{\mu}=p^{\mu}_{R}},$ (21)

where $b^{\mu}$ denotes the separation of the endpoints of the nonlocal
operator $\mathcal{O}_{\Gamma}$, the projected vertex function is defined as

$\mathcal{V}^{\mathcal{O}_{\Gamma\Gamma^{\prime}}}(p)\equiv\text{Tr}\left[\Lambda^{\mathcal{O}_{\Gamma}}(p)\Gamma^{\prime}\right],$
(22)

and the replacement

$\displaystyle\text{Tr}\left[\Lambda^{\mathcal{O}_{\Gamma};\text{tree}}(p_{R})\Gamma^{\prime}\right]=6e^{ip_{R}\cdot
b}\delta^{\Gamma\Gamma^{\prime}}$ (23)

has been made. It should be noted that since the operator
$\mathcal{O}^{q}_{\Gamma}$ is nonlocal and frame dependent, different
directions in $p_{R}^{\mu}$ constitute different renormalization schemes,
related by finite renormalization factors. As such,
$Z^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\Gamma\Gamma^{\prime}}}$
depends on $p_{R}^{\mu}$ itself rather than only on its magnitude, which acts
as the renormalization scale.

The complete $16\times 16$ matrix of $\mathrm{RI}^{\prime}\mathrm{/MOM}$
renormalization factors
$Z^{\mathrm{RI}^{\prime}\mathrm{/MOM}}_{\mathcal{O}_{\Gamma\Gamma^{\prime}}}$
is computed for the same set of operator geometries defined in Table 1, on
$n_{\text{cfg}}=50$ gauge field configurations. For all operator geometries
with $\eta/a\in\\{12,14\\}$, the renormalization factors are computed for a
range of four-momenta tabulated in Tab. 2, to enable an investigation of
residual dependence on the renormalization scale (which would be canceled by
an all-orders matching to the $\overline{\mathrm{MS}}$ scheme) and
discretization artifacts. For operator geometries with $\eta/a=23$, only the
four-momentum with $n^{\mu}=(12,12,12,12)$ is used. An example of the
resulting $\mathrm{RI}^{\prime}\mathrm{/MOM}$ renormalization matrices is
shown in Fig. 4, which displays a subset of the off-diagonal renormalization
factors normalized relative to the diagonal components:

$\displaystyle\mathcal{M}^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\Gamma\mathcal{P}}}\equiv$
$\displaystyle\frac{\text{Abs}[Z^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\Gamma\mathcal{P}}}]}{\frac{1}{16}\sum_{i}\text{Abs}[Z^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\Gamma_{i}\Gamma_{i}}}]}\,,$
(24)

computed for quark bilinear operators with straight Wilson lines and a
particular choice of momentum.

$n^{\mu}$ |  $\sqrt{p^{2}}$ [GeV] |  $p^{z}$ [GeV] | $p^{[4]}/(p^{2})^{2}$
---|---|---|---
(6,6,6,6) | 2.5 | 1.3 | 0.26
(6,6,6,9) | 2.7 | 1.3 | 0.26
(6,6,6,12) | 3.0 | 1.3 | 0.30
(8,8,8,8) | 3.3 | 1.7 | 0.26
(8,8,8,12) | 3.6 | 1.7 | 0.26
(8,8,8,16) | 4.0 | 1.7 | 0.30
(12,12,12,12) | 4.9 | 2.6 | 0.26
(12,12,12,18) | 5.4 | 2.6 | 0.25
Table 2: Four-momenta used in the calculation of nonperturbative
renormalization factors as described in the text, where $p^{\mu}$ is the four-
momentum in physical units corresponding to $n^{\mu}$ in lattice units. The
H(4) invariant $p^{[4]}$ is defined as $p^{[4]}=\sum_{\mu=1}^{4}p_{\mu}^{4}$.
Figure 4: Submatrix of the $\mathrm{RI}^{\prime}\mathrm{/MOM}$ mixing matrix
$\mathcal{M}^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\Gamma\mathcal{P}}}$,
defined in Eq. (24), computed for quark bilinear operators with straight
Wilson lines ($b_{T}=0$) with various extents $b^{z}$, for momentum
$n^{\nu}=(12,12,12,12)$ in lattice units. Points representing Dirac structures
in the upper triangle of the mixing matrix are slightly offset on the
horizontal axis for clarity.

The $\overline{\mathrm{MS}}$ renormalization factors are computed via Eq. (7)
from the $\mathrm{RI}^{\prime}\mathrm{/MOM}$ renormalization matrices and the
perturbative matching factor
$\mathcal{R}^{\overline{\mathrm{MS}}}_{\mathcal{O}_{\gamma^{4}\Gamma}}$,
calculated at next-to-leading order in continuum perturbation theory with
dimensional regularization ($D=4-2\epsilon$) Constantinou  _et al._ (2019);
Ebert  _et al._ (2020). In this work, the scale $\tilde{p}_{R}$ is set to 4.9
GeV. Examples of the resulting $\overline{\mathrm{MS}}$ renormalization
factors are presented in Figs. 5 and 6. This renormalization is independent of
$p_{R}$ up to discretization artifacts, two-loop perturbative matching
corrections which are neglected here, and nonperturbative effects that vanish
at asymptotically large $p^{2}_{R}$. While in principle one might fit to
results generated at different $p_{R}$ values to constrain discretization
effects, this is not feasible with the data generated here, and the covariance
matrices for the nonlocal operators cannot be reliably estimated. The
renormalization constants computed with $n^{\mu}=(12,12,12,12)$ are thus taken
as best-estimates and used in the further analysis of the Collins-Soper
kernel. For those operator structures where the
$\mathrm{RI}^{\prime}\mathrm{/MOM}$ renormalization factors have been computed
at other choices of $p_{R}$, this yields indistinguishable results for the
Collins-Soper kernel compared with results obtained using the weighted
averaging procedure over momenta which is employed in Ref. Shanahan  _et al._
(2019). The components of
$\mathcal{M}^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\Gamma\mathcal{P}}}$
are of similar magnitude for both $\Gamma\in\\{\gamma^{3},\gamma^{4}\\}$,
which is consistent with the conclusion of Ref. Ji  _et al._ (2021) but is in
contrast with the results of Refs. Constantinou  _et al._ (2019); Green  _et
al._ (2020) which predict no mixing effects for $\Gamma=\gamma^{3}$ at ${\cal
O}(a^{0})$. The quasi beam functions with both Dirac structures are thus
treated on equal footing in this work.

Figure 5: Example of numerical results for
$Z^{\text{$\mathrm{RI}^{\prime}\mathrm{/MOM}$}}_{\mathcal{O}_{\gamma^{4}\gamma^{4}}}$
(orange circles) and
$Z^{\overline{\mathrm{MS}}}_{\mathcal{O}_{\gamma^{4}\gamma^{4}}}$ (blue
squares), computed using the different values of $p_{R}$ given in Table 2, for
operator geometry $\eta/a=14$, $b^{z}/a=3$, $b_{T}/a=3$, at renormalization
scale $\mu=2$ GeV. The blue shaded band shows the value used in further
analysis, obtained as described in the text.

(a)

(b)

Figure 6:  Diagonal $\overline{\text{MS}}$ renormalization factors
$Z^{\overline{\text{MS}}}_{\mathcal{O}_{\gamma^{4}\gamma^{4}}}$, computed for
operator geometry $\eta/a=14$ and various values of $b_{T}$ and $b^{z}$, at
renormalization scale $\mu=2$ GeV.

###  III.3 Renormalized quasi beam functions

Modified $\overline{\mathrm{MS}}$-renormalized quasi beam functions are
computed via Eq. (6) from the bare quasi beam functions and
$\overline{\mathrm{MS}}$ renormalization factors calculated as described in
Secs. III.1 and III.2; the uncertainties of the two components are combined in
quadrature. While for clarity all equations in this section are expressed for
renormalized quasi beam functions defined with the Dirac structure
$\gamma^{4}$, both $B^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ and
$B^{\overline{\mathrm{MS}}}_{\gamma^{3}}$ are computed and analyzed
independently.

(a) Example of the $b^{z}$-dependence of the asymmetry in the renormalized
quasi beam functions for a fixed choice of operator geometry with $\eta/a=14$,
$b_{T}/a=1$.  (b) Example of the $b_{T}$-dependence of the asymmetry parameter
$\Delta=-V(b_{T})+V(b_{T}^{R})$ for operators with $\eta/a=14$, fit to
renormalized quasi beam functions $B^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ as
described in the text. The approximate linear dependence of $\Delta$ on
$b_{T}$ can be explained by a linear term in $V(b_{T})$; the approximate
independence on $n^{z}$ is expected, as $V(b_{T})$ should be independent of
the external state.

Figure 7: Illustration of the asymmetry in the renormalized quasi beam
functions computed in this work; additional examples are provided in Appendix
B.

The real and imaginary parts of the renormalized quasi beam functions should
be symmetric and antisymmetric functions of $b^{z}$ respectively in the
$\eta\rightarrow\infty$ limit. The numerical results obtained in this work,
however, show significant departures from these expectations at finite $\eta$
as was also observed in the quenched calculation of Ref. Shanahan  _et al._
(2019). The $b^{z}$-dependence of the asymmetry in the absolute value of the
renormalized quasi beam functions is illustrated in Fig. 7a. This asymmetry
can be understood as an incomplete cancellation of the Wilson-line self-energy
correction proportional to $e^{-V(b_{T})\,b^{z}}$ between the
$\overline{\mathrm{MS}}$ renormalization factor and the bare quasi beam
function, where $V(b_{T})$ is the static quark-antiquark potential at distance
$b_{T}$. Such an effect yields a $b^{z}$-dependent asymmetry proportional to
$e^{-\Delta b^{z}}$, depending on an asymmetry parameter
$\Delta=-V(b_{T})+V(b_{T}^{R})$. This is in fact a good model of the asymmetry
observed in the numerical calculations of this work; fits of
$|B^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,-b^{z},\vec{b}_{T},a,\eta,P^{z})|/|B^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,b^{z},\vec{b}_{T},a,\eta,P^{z})|$
(and the analogous expression constructed from
$B^{\overline{\mathrm{MS}}}_{\gamma^{3}}$) to this functional form, fit
separately for fixed values of $b_{T}$, $\eta$, and $P^{z}$, achieve
acceptable goodness-of-fit with an average $\chi^{2}/\text{d.o.f.}=0.56$.
These fits, and the resulting values of the asymmetry parameter $\Delta$, are
illustrated in Fig. 7 and in Appendix B. Asymmetry-corrected modified
$\overline{\mathrm{MS}}$-renormalized quasi beam functions are thus defined as

$\displaystyle B^{\overline{\mathrm{MS}};\text{corr}}_{\gamma^{4}}(\mu,b^{z},$
$\displaystyle\vec{b}_{T},a,\eta,P^{z})=\
e^{\Delta(\vec{b}_{T},a,\eta,P^{z})|b^{z}|}$
$\displaystyle\times\left(\text{Re}\vphantom{B^{\overline{\mathrm{MS}}}_{\gamma^{4}}}\right.$
$\displaystyle\\!\\!\left[B^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,|b^{z}|,\vec{b}_{T},a,\eta,P^{z})\right]$
$\displaystyle+$
$\displaystyle\left.\text{sign}(b^{z})\,\text{Im}\\!\\!\left[B^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,|b^{z}|,\vec{b}_{T},a,\eta,P^{z})\right]\right).$
(25)

Here, the uncertainties in the factor $e^{\Delta|b^{z}|}$ and the quasi beam
functions are added in quadrature and, after asymmetry correction, the more
precise results for beam functions with $b^{z}>0$ (which involve shorter
Wilson lines) are mirrored to $b^{z}<0$. An example of the modified
$\overline{\mathrm{MS}}$-renormalized quasi beam functions before and after
asymmetry correction is given in Fig. 8.

(a)

(b)

Figure 8:  Comparison of an example of the modified
$\overline{\mathrm{MS}}$-renormalized quasi beam functions before (gray) and
after (color) asymmetry correction via Eq. (25), for beams functions computed
with operator geometry $\eta/a=14$, $b_{T}/a=1$.

(a)  (b)

Figure 9:  Example of the asymmetry-corrected modified
$\overline{\mathrm{MS}}$-renormalized quasi beam function
$B^{\overline{\mathrm{MS}};\text{corr}}_{\gamma^{4}}$ defined in Eq. (6). The
horizontal shaded bands show the results of weighted averages of this quantity
over choices of $b_{T}^{R}$ and $\eta$ (as a function of $b^{z}$, $b_{T}$, and
$P^{z}$), as described in the text.

It is expected that the renormalized beam functions should be independent of
the matching scale $b_{T}^{R}$. This expectation is satisfied within
uncertainties for the numerical investigations of this work as illustrated in
Fig. 9; for use in the calculation of the Collins-Soper kernel, a weighted
average Aoki  _et al._ (2020) over possible choices of $b_{T}^{R}$ in the
window $a\ll b_{T}^{R}\ll\Lambda_{\text{QCD}}^{-1}$ is thus implemented,
performed precisely as detailed in Appendix C of Ref. Shanahan  _et al._
(2020). Similarly, the asymmetry-corrected renormalized quasi beam functions
do not depend on $\eta$ within uncertainties, and the formal extrapolation to
$\eta\rightarrow\infty$ is implemented with an analogous weighted average. The
resulting averaged values of the asymmetry-corrected modified
$\overline{\mathrm{MS}}$-renormalized quasi beam function are denoted by
$\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,b^{z},b_{T},P^{z})$,
and have no dependence on $\eta$ or $b^{R}_{T}$; the dependence on $a$ is also
dropped in this notation, as a single lattice spacing is used in this
numerical study. Figures showing
$\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ and
$\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$ are provided in Appendix
B. The contributions to these quantities from mixing between operators with
different Dirac structures is generally found to be less than $10\%$; see Fig.
10 for an illustration. Additional examples are provided in Appendix B.

(a)  (b)

Figure 10: Percentage contribution to the renormalized quasi beam functions
from mixing of operators with different Dirac structures. Note that the ratios
shown are outside of the plot range near the nodes of the beam functions; in
this example the maximum mixing that is resolved from zero at greater than
$2\sigma$ is $0.32(5)$, and occurs in the real component of the beam function
for $n_{z}=5$, $b_{z}/a=5$.

###  III.4 Collins-Soper kernel

To determine the Collins-Soper kernel via Eq. (II) from the averaged
asymmetry-corrected modified $\overline{\mathrm{MS}}$-renormalized quasi beam
functions $\overline{B}_{\gamma^{4}}^{\overline{\mathrm{MS}}}$ and
$\overline{B}_{\gamma^{3}}^{\overline{\mathrm{MS}}}$, defined in the previous
section, requires a Fourier transform of the quasi beam functions in $b^{z}$.
As is clear from Fig. 23, however, the $b^{z}$-range of the data is not
sufficient for the tails of the quasi beam functions at large $|b^{z}|$ to
decay to plateaus consistent with zero, particularly at the largest $b_{T}$
and smallest $P^{z}$ values used in this numerical study. As such, it is to be
expected that a discrete Fourier transform (DFT) will have significant
systematic uncertainties from the truncation of the data in $P^{z}b^{z}$;
details of a DFT analysis of the quasi beam functions are presented in
Appendix C. Instead, Fourier transforms are taken after fitting the quasi beam
functions to functional forms that allow extrapolation in $b^{z}$. This
approach trades the systematic uncertainties of a DFT for model-dependence in
the fit form used to extrapolate the quasi beam functions. While this model-
dependence is difficult to quantify rigorously, this approach allows the
physical expectation that the quasi beam functions should decay smoothly at
large $|b^{z}|$ to be incorporated (a DFT effectively models the beam
functions as zero outside the $b^{z}$-range in which lattice QCD results are
computed).

In particular, the quasi beam functions are modeled as a sum of polynomials in
$b^{z}$ times Gaussian functions, which provide an appropriate basis, since it
is expected that the quasi beam functions should be analytic functions of
$b^{z}$ (the $b^{z}\to 0$ limit at fixed $b_{T}$ does not introduce additional
divergences), and yield high-quality fits with few free parameters.
Specifically, for each choice of $b_{T}$ and $P^{z}$, the real and imaginary
parts of the quasi beam function
$\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ (and independently
$\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$) are fit with even and
odd functions of $b^{z}$ respectively, defined as

$\displaystyle f_{\text{Re}}(\sigma,\\{r_{n}\\};b^{z})$
$\displaystyle=\text{exp}[-(b^{z})^{2}/(2\sigma^{2})]\sum_{n=0}^{n_{\text{max}}}r_{n}(b^{z})^{2n}$
(26) $\displaystyle f_{\text{Im}}(\sigma,\\{r_{n}\\};b^{z})$
$\displaystyle=\text{exp}[-(b^{z})^{2}/(2\sigma^{2})]\sum_{n=0}^{n_{\text{max}}}i_{n}(b^{z})^{2n+1}.$
(27)

The value of $n_{\text{max}}$ is chosen to minimize the Akaike information
criterion (AIC) Akaike (1974) and corresponds to
$n_{\text{max}}\in\\{2,3,4\\}$ for all cases. The fits with these optimal
values of $n_{\text{max}}$ are of high quality in all cases, with an average
$\chi^{2}/\text{d.o.f.}=0.41$. The resulting models of the quasi beam
functions are denoted $\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ (and,
correspondingly $\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$), and are
illustrated in Fig. 11 (with further examples provided in Appendix B).

(a)

(b)

Figure 11:  Example of fits by Eqs. (26) and (27) (shaded bands) to the real
and imaginary parts of the quasi beam functions, for $b_{T}/a=1$.

Finally, in terms of the models
$\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$, the relation defining the
Collins-Soper kernel in Eq. (II) is realized as

$\displaystyle\hat{\gamma}^{q}_{\zeta}(\mu,b_{T};P_{1}^{z},P_{2}^{z},x)$
$\displaystyle\equiv\frac{1}{\ln(P^{z}_{1}/P^{z}_{2})}\ln\Biggr{[}\frac{C^{\mathrm{TMD}}_{\text{ns}}(\mu,xP_{2}^{z})}{C^{\mathrm{TMD}}_{\text{ns}}(\mu,xP_{1}^{z})}$
$\displaystyle\\!\times\\!\frac{\int\\!\mathrm{d}b^{z}e^{-ib^{z}\\!xP_{1}^{z}}P_{1}^{z}\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,b^{z},b_{T},P_{1}^{z})}{\int\\!\mathrm{d}b^{z}e^{-ib^{z}\\!xP_{2}^{z}}\\!P_{2}^{z}\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,b^{z},b_{T},P_{2}^{z})}\Biggr{]},$
(28)

which coincides with $\gamma^{q}_{\zeta}(\mu,b_{T})$ up to power corrections
such as higher-twist corrections in the factorization formula for the quasi
TMDPDF, and discretization artifacts, which introduce the dependence on
$P_{1}^{z}$, $P_{2}^{z}$, and $x$. One approach to determine
$\gamma^{q}_{\zeta}(\mu,b_{T})$ from
$\hat{\gamma}^{q}_{\zeta}(\mu,b_{T};P_{1}^{z},P_{2}^{z},x)$ is to model, fit,
and subtract, these various artifacts. However, the most straightforward
models of these effects do not provide good fits to the numerical data of this
study, as detailed in Appendix D. Instead, since the contamination in
$\hat{\gamma}^{q}_{\zeta}$ will be different at each choice of $P_{1}^{z}$,
$P_{2}^{z}$, and $x$, and the effects can be expected to be larger222The
matching coefficient includes large logarithms of $xP^{z}_{i}$ at small $x$,
while the quasi beam functions at $x\to 0$ and $x\to 1$ are sensitive to the
long-range correlations in $b^{z}$ and are thus affected by the truncation of
the data in $P^{z}b^{z}$. In addition, the power corrections are expected to
be enhanced at small $x$. at large and small values of $x$ and at small values
of $P^{z}$, the variation of $\hat{\gamma}^{q}_{\zeta}$ over these choices is
used to define an estimate of the systematic uncertainty.

Precisely, a best value for the Collins-Soper kernel is determined from
$\hat{\gamma}^{q}_{\zeta}$ via a multi-step procedure. First, the largest
window of $x$ is determined for which the data for all choices of the pair
$\\{P_{1}^{z},P_{2}^{z}\\}$ are consistent with a common constant value. In
practice this region is defined as the largest window in which a constant fit
to the data at a set of discrete $x$ points has a $\chi^{2}/\text{d.o.f.}\leq
1$. The central value and uncertainty are defined as the median and the 68%
confidence interval of the union of the bootstrap data in that $x$ window,
including all $\\{P_{1}^{z},P_{2}^{z}\\}$ pairs. The result of this procedure
is robust to changes in the discretization of $x$, for sufficiently fine
discretization scales (100 points spanning $0<x<1$ uniformly are used in the
analysis presented). This procedure is performed separately for
$\hat{\gamma}^{q}_{\zeta}$ determined from beam functions calculated with
Dirac structures $\gamma^{4}$ and $\gamma^{3}$; examples of the resulting
values are shown with $\hat{\gamma}^{q}_{\zeta}$ in Fig. 12, with the
remainder presented in Appendix B. The central values of the independent
calculations with Dirac structures $\gamma^{4}$ and $\gamma^{3}$ are averaged,
and the average uncertainty is added in quadrature with half the difference
between the central values obtained using each Dirac structure, to yield the
final results of this work which are shown as a function of $b_{T}$ in Fig.
13, and are tabulated in Table. 3.

$b_{T}$ [fm] | 0.12 | 0.24 | 0.36 | 0.48
---|---|---|---|---
$\gamma^{q,\overline{\mathrm{MS}}}_{\zeta}$ | -0.419(53)(50) | -0.49(5)(12) | -0.76(9)(8) | -0.82(15)
Table 3: Collins-Soper kernel with $\mu=2\text{ GeV}$ as a function of
$b_{T}$. The first uncertainty is the average of that determined from
calculations using $\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ and
$\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$ as described in the text,
while the second is a systematic uncertainty computed as half the difference
of the central values of the results obtained using quasi beam functions
defined with the two Dirac structures.

(a) $\hat{\gamma}^{q}_{\zeta}$ computed from quasi beam functions
$\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$.

(b) $\hat{\gamma}^{q}_{\zeta}$ computed from quasi beam functions
$\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$.

Figure 12: $\hat{\gamma}^{q}_{\zeta}$, computed as defined in Eq. (28) for all
momentum pairs $\\{P_{1}^{z},P_{2}^{z}\\}$, denoted by $P_{1}^{z}/P_{2}^{z}$
in the legend. The horizontal shaded band shows the fit window in $x$, as well
as the total uncertainty of the best result, determined as described in the
text. Figure 13: Collins-Soper kernel as a function of $b_{T}$, determined
from $\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ (purple circles)
$\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$ (red triangles), and the final
combined results of this work (green squares), computed as described in the
text. For the latter points, the inner (outer) error bars show the first
(quadrature-combined) uncertainties given in Table. 3. No result computed from
$\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$ is shown at the largest
$b_{T}$ value because in this case $\hat{\gamma}^{q}_{\zeta}$ cannot be fit to
a constant with the method described in the text, as shown in Fig. 25.  Figure
14:  $b_{T}$-dependence of the Collins-Soper kernel computed from the same
quasi beam functions via the different approaches defined in Sec. III.4. All
points other than the primary results of this work (“NLO”) are offset on the
horizontal axis for clarity. For the “NLO” and “LO” approaches, results
computed based on quasi beam functions with Dirac structures $\gamma^{4}$ and
$\gamma^{3}$ are combined as described in the text; the outer error bars
include half the difference between the results with $\gamma^{4}$ and
$\gamma^{3}$ combined in quadrature with the average uncertainty, shown by the
inner error bars. For the other approaches the empty (filled) points show
results obtained with Dirac structure $\gamma^{4}$ ($\gamma^{3}$).
“Hermite/Bernstein” points with Dirac structure $\gamma_{3}$ are not shown at
$b_{T}/a=4$ because the corresponding fits of the $P^{z}b^{z}$-dependence of
the relevant quasi beam functions were of poor quality, as described in the
text.

In addition to the approach followed here, there are a number of alternative
methods of extracting the Collins-Soper kernel that have been proposed or
employed in other studies, for example:

  * •

“LO”: The perturbative matching coefficient $C^{\mathrm{TMD}}_{\text{ns}}$
computed to leading-order (LO), instead of NLO, can be used in an analysis
otherwise mirroring that presented here;

  * •

“Hermite/Bernstein”: As proposed in Ref. Shanahan  _et al._ (2020), the
$P^{z}b^{z}$-dependence of the quasi beam functions can be fit to models based
on Hermite and Bernstein polynomial bases constructed to yield $x$-independent
Collins-Soper kernels via Eq. (28), taking the LO value of the perturbative
matching coefficient $C^{\mathrm{TMD}}_{\text{ns}}$;

  * •

“$b^{z}=0$”: An approximation of the Collins-Soper kernel can be computed with
LO matching using only the quasi beam functions evaluated at $b^{z}=0$ (this
approach does not require a Fourier transform in $b^{z}$):

$\hskip
25.60747pt[{\gamma}^{q}_{\zeta}(\mu,b_{T})]^{b^{z}=0}\equiv\frac{1}{\ln(P^{z}_{1}/P^{z}_{2})}\ln\Biggr{[}\frac{\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,0,b_{T},P_{1}^{z})}{\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,0,b_{T},P_{2}^{z})}\Biggr{]};$
(29)

  * •

“$b^{z}=0$, bare”: As proposed in Ref. Zhang  _et al._ (2020), the same
approach described for “$b^{z}=0$” can be followed, using bare quasi beam
functions ${B}^{\text{bare}}_{\gamma^{4}}$ rather than renormalized quasi beam
functions (i.e., neglecting operator mixing between different Dirac
structures);

  * •

“$b^{z}=0/b_{T}=0$, bare”: As proposed in Ref. Li  _et al._ (2021), a
variation of the ‘$b^{z}=0$” approach can be used, approximating the Collins
Soper kernel as

$\displaystyle\hskip 11.38109pt[{\gamma}^{q}_{\zeta}$
$\displaystyle(\mu,b_{T})]^{b^{z}=0/b_{T}=0}\equiv\frac{1}{\ln(P^{z}_{1}/P^{z}_{2})}$
$\displaystyle\times\ln\Biggr{[}\frac{{B}^{\text{bare}}_{\gamma^{4}}(0,b_{T},a,\eta,P_{1}^{z}){B}^{\text{bare}}_{\gamma^{4}}(0,0,a,\eta,P_{2}^{z})}{{B}^{\text{bare}}_{\gamma^{4}}(0,b_{T},a,\eta,P_{2}^{z}){B}^{\text{bare}}_{\gamma^{4}}(0,0,a,\eta,P_{1}^{z})}\Biggr{]}.$
(30)

Each of these methods can be followed using the quasi beam functions computed
in this work; a comparison of the results is provided in Fig. 14. For the “LO”
approach, the same procedure is followed to combine the results obtained using
quasi beam functions with Dirac structures $\gamma^{3}$ and $\gamma^{4}$ as
for the “NLO” method which yields the main results of this work. For the other
approaches the results obtained with the two Dirac structures are not always
consistent at one standard deviation, and are shown separately; for
$b_{T}/a=4$ no results for the “Hermite/Bernstein” approach are shown with
Dirac structure $\gamma^{3}$ as the model fits were of poor quality with
$\chi^{2}/\text{d.o.f.}>2$. In the case of the “$b^{z}=0/b_{T}=0$, bare”
approach, bare quasi beam functions with $\eta/a=14$, which is the largest
extent studied in this work for all $P^{z}$, are used in the analysis.

(a) Comparison with the SV19 Scimemi and Vladimirov (2019) and Pavia19
Bacchetta  _et al._ (2019) phenomenological parameterizations and the next-to-
next-to-next-to-leading order (N3LO) perturbative result Li and Zhu (2017);
Henn  _et al._ (2020). (b)  Comparison with quenched results of Ref. Shanahan
_et al._ (2020) (SWZ), as well as results from the LPC Zhang  _et al._ (2020),
Regensburg/NMSU Schlemmer  _et al._ (2021), and ETMC/PKU Li  _et al._ (2021)
collaborations. Different sets of points with the same color show different
sets of results from the same collaboration.

Figure 15:  $b_{T}$-dependence of the Collins-Soper kernel as determined in
this work (green squares in both panels) compared with (a) phenomenological
results, and (b) the results of other lattice QCD calculations of this
quantity.

Clearly, although the same quasi beam functions are used, the Collins-Soper
kernel determined via each of these approaches is very different, and many of
the results are inconsistent with the best results of this study at several
standard deviations, with uncertainties as much as an order of magnitude
smaller. This is to be expected if the effects of higher-order matching,
renormalization and mixing, and power corrections, are significant, as each of
the approaches listed above treats one or more of these systematic effects
differently than in the primary analysis presented here.

##  IV Outlook

This work presents a determination of the Collins-Soper kernel from a
dynamical lattice QCD calculation following the approach of Refs. Ebert  _et
al._ (2019a, b). Several systematic uncertainties remain to be addressed; in
particular, the quark masses used correspond to an unphysically-large pion
mass of $m_{\pi}=538(1)$ MeV, and the results are obtained using a single
ensemble of gauge field configurations such that effects from the
discretization and finite lattice volume cannot be fully quantified. A fully
model-independent calculation will require these systematics to be addressed,
lattice QCD calculations to be performed over a larger range of $P^{z}b^{z}$
to eliminate the need to extrapolate the quasi beam functions to large
$|b^{z}|$ and enable the DFT approach to be used, and larger values of $P^{z}$
to be included to reduce the contributions from power corrections and higher-
twist effects which dominate the uncertainties of this calculation. With these
caveats in mind, the results of this work may be compared with
phenomenological extractions of the Collins-Soper kernel, as shown in Fig.
15a. The lattice QCD and phenomenological determinations are broadly
consistent at large $b_{T}$, with clear deviations at the smallest $b_{T}$
values studied; discretization effects are expected to be largest at small
$b_{T}$ and might be relevant for understanding this effect. It is clear that,
while challenging to achieve computationally, future fully-controlled
calculations by this approach have the potential to differentiate different
models of the Collins-Soper kernel and will provide important input for the
analysis of low-energy SIDIS data and the determinations of the TMDPDFs.

In considering the prospects for such future controlled determinations of the
Collins-Soper kernel from lattice QCD, it is informative to contrast the
results of this study with those of other lattice QCD investigations; a
comparison of existing calculations Shanahan  _et al._ (2020); Zhang  _et al._
(2020); Schlemmer  _et al._ (2021); Li  _et al._ (2021) is provided in Fig.
15b. All dynamical calculations use quark masses resulting in similar values
of the pion mass to that of the calculation presented here (ranging from the
lightest ensemble with $m_{\pi}=350$ MeV in Ref. Li  _et al._ (2021) to
$m_{\pi}=547$ MeV in Ref. Zhang  _et al._ (2020)), while the quenched
calculation of Ref. Shanahan  _et al._ (2020), in which the kernel should not
depend on the valence quark masses since it is independent of the external
state, is performed at $m_{\pi}=1.207$ GeV. Each calculation uses a slightly
different approach to constrain the Collins-Soper kernel from quasi beam
functions. In particular, the “Hermite/Bernstein” approach is followed in Ref.
Shanahan  _et al._ (2020) (“SWZ”), the calculation of Ref. Zhang  _et al._
(2020) (“LPC”) uses the “$b^{z}=0$, bare” approach, that of Ref. Schlemmer
_et al._ (2021) (“Regensburg/NMSU”) uses an approach similar to the
“$b^{z}=0$, bare” approach but with NLO matching, and Ref. Li  _et al._ (2021)
(“ETMC/PKU) applies the “$b^{z}=0/b_{T}=0$, bare” approach. While the various
calculations exhibit similar dependence on $b_{T}$, there are some significant
discrepancies between the numerical results, and a wide range of uncertainty
estimates. Given the analysis of Sec. III.4, this is to be expected; even when
the same quasi beam function data is used, following the various “$b^{z}=0$”
approaches and the approach presented here result in significant systematic
differences, and significantly different uncertainty estimates. Since Refs.
Zhang  _et al._ (2020); Schlemmer  _et al._ (2021); Li  _et al._ (2021) all
use somewhat larger maximum $P^{z}$ values than that of the present study, the
effects of power corrections and higher-twist contamination can be expected to
be smaller than those found in Sec. III.4, but these effects, together with
the difference between NLO and LO matching illustrated in Appendix E, could
nevertheless be responsible for the discrepancies visible in Fig. 15b.
Clearly, a fully-controlled determination of the Collins-Soper kernel from
lattice QCD will require NLO or even higher-order matching or resummation and
a treatment of power corrections that is more robust than that achieved in any
of the studies performed to date. The analysis presented here constitutes an
important step towards that goal, revealing clearly that these important
sources of systematic uncertainty cannot be neglected.

###### Acknowledgements.

The authors thank Will Detmold, Iain Stewart and Alexey Vladimirov for useful
discussions, Xu Feng for providing numerical data for inclusion in Fig. 15,
and Steven Gottlieb and the MILC collaboration for the use of the gauge
configurations used in this project, which were generated at Indiana
University on Big Red 2+ and Big Red 3. This research was supported in part by
Lilly Endowment, Inc., through its support for the Indiana University
Pervasive Technology Institute that provides the Big Red supercomputers. This
work is supported in part by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, under grant Contract Numbers DE-SC0011090, DE-
AC02-06CH11357, DE-SC0012704, and within the framework of the TMD Topical
Collaboration. PES is additionally supported by the U.S. DOE Early Career
Award DE-SC0021006, by a NEC research award, by the Carl G and Shirley
Sontheimer Research Fund, and by the U.S. National Science Foundation under
Cooperative Agreement PHY-2019786 (The NSF AI Institute for Artificial
Intelligence and Fundamental Interactions, http://iaifi.org/). This research
used resources of the National Energy Research Scientific Computing Center
(NERSC), a U.S. Department of Energy Office of Science User Facility operated
under Contract No. DE-AC02-05CH11231, the Extreme Science and Engineering
Discovery Environment (XSEDE) Bridges-2 at the Pittsburgh Supercomputing
Center (PSC) through allocation TG-PHY200036, which is supported by National
Science Foundation grant number ACI-1548562, and facilities of the USQCD
Collaboration, which are funded by the Office of Science of the U.S.
Department of Energy. This manuscript has been authored by Fermi Research
Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of
Energy, Office of Science, Office of High Energy Physics. The Chroma Edwards
and Joo (2005), QLua Pochinsky , QUDA Clark  _et al._ (2010); Babich  _et al._
(2011); Clark  _et al._ (2016), QDP-JIT Winter  _et al._ (2014), and QPhiX Joó
_et al._ (2016) software libraries were used in this work. Data analysis used
NumPy Harris  _et al._ (2020) and Julia Bezanson  _et al._ (2017); Mogensen
and Riseth (2018), and figures were produced using Mathematica Inc. .

##  Appendix A Fits to ratios of three and two point functions

As detailed in Sec. III.1, ratios
$\mathcal{R}_{\Gamma}(t,\tau,b^{\mu},a,\eta,P^{z})$ (Eq. (11)) of three-point
and two-point correlation functions asymptote to the bare quasi beam functions
in the limit $\\{\tau,t-\tau\\}\rightarrow\infty$, with contamination from
matrix elements in excited states present at finite $t$ and $\tau$. The $t$
and $\tau$-dependence of the ratios is fit precisely as defined in Appendix A
of Ref. Shanahan  _et al._ (2020) to extract the bare quasi beam functions:
fits are performed for all different possible cuts on source/operator/sink
separations, with the AIC Akaike (1974) used to select the number of states
included in the spectral representation for each fit. The results are combined
via a weighted averaging procedure. Some examples of the results of this
fitting procedure are given in Fig. 16.

Figure 16: Examples of fits to ratios of three and two-point functions
$\mathcal{R}_{\Gamma}(t,\tau,b^{\mu},a,\eta,P^{z})$ defined in Eq. (11),
performed as described in the text. The shaded bands matching the color of
each set of points show 68% bootstrap confidence intervals of the fits from
the fit range that has the highest weight in the weighted average that
determines the final result, which is indicated by the gray horizontal bands.
These bands show the total uncertainty on the bare quasi beam functions,
including both the statistical uncertainty and the systematic uncertainty
which arises from variation of the results between different fit range
choices.

##  Appendix B Additional beam function results

This section collates additional examples of results at intermediate states of
the numerical calculation and analysis presented in Sec. III.4:

  * •

Bare quasi beam functions
$B_{\Gamma}^{\text{bare}}(b^{z},\vec{b}_{T},a,\eta,P^{z})$: supplementing Fig.
3 of the main text, additional examples of the bare quasi beam functions are
provided in Fig. 17 for $B_{\gamma^{4}}^{\text{bare}}$ and Fig. 18 for
$B_{\gamma^{3}}^{\text{bare}}$.

  * •

Renormalized quasi beam functions
$B^{\overline{\mathrm{MS}}}_{\Gamma}(\mu,b^{z},\vec{b}_{T},a,\eta,P^{z})$:
examples of the modified $\overline{\mathrm{MS}}$-renormalized quasi beam
functions $B^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ and
${B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$, defined in Sec. III.3, are given
in Figs. 19 and 20 respectively.

  * •

Renormalized quasi beam function asymmetry fits: supplementing Fig. 7 of the
main text, fits to the $b^{z}$-dependence of the asymmetry in the renormalized
quasi beam functions, performed as detailed in Sec. III.3, are illustrated in
Fig. 21.

  * •

Asymmetry-corrected modified $\overline{\mathrm{MS}}$-renormalized quasi beam
functions
$B^{\overline{\mathrm{MS}};\text{corr}}_{\Gamma}(\mu,b^{z},\vec{b}_{T},a,\eta,P^{z})$:
an example of the dependence of
$B^{\overline{\mathrm{MS}};\text{corr}}_{\gamma^{3}}$ on $b_{T}^{R}$ and
$\eta$ is provided in Fig. 22, supplementing the analogous figure for
$B^{\overline{\mathrm{MS}};\text{corr}}_{\gamma^{4}}$ which is presented in
Fig. 9 in the main text.

  * •

Averaged asymmetry-corrected modified $\overline{\mathrm{MS}}$-renormalized
quasi beam functions
$\overline{B}^{\overline{\mathrm{MS}}}_{\Gamma}(\mu,b^{z},b_{T},P^{z})$: in
addition to the example provided in Fig. 11 of the main text, Figs. 23 and 24
give examples of $\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ and
$\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$, including fits to the
$P^{z}b^{z}$-dependence of these functions and extrapolations to larger
$P^{z}b^{z}$, performed as described in Sec. III.4.

  * •

Estimator $\hat{\gamma}^{q}_{\zeta}(\mu,b_{T};P_{1}^{z},P_{2}^{z},x)$ for the
Collins-Soper kernel (Eq. (28)): supplementing Fig. 12 of the main text, the
remaining numerical results for $\hat{\gamma}^{q}_{\zeta}$ as a function of
$x$ at different values of $b_{T}$ are displayed in Fig. 25.

Figure 17: Examples of the bare quasi beam functions
$B^{\text{bare}}_{\gamma^{4}}$ determined as described in Sec. III.1 (note
that $B^{\text{bare}}_{\gamma^{4}}$ for $b_{T}=0.12$ fm and $\eta=1.68$ fm is
shown in Fig. 3 in the main text).

##  Appendix C Discrete Fourier transform analysis

As discussed in Sec. III.4, a model-independent lattice QCD extraction of the
Collins-Soper kernel by the approach followed here would require that Eq. (28)
is evaluated with a DFT of $\overline{B}^{\overline{\mathrm{MS}}}_{\Gamma}$
replacing the Fourier transform of analytic fits
$\hat{B}^{\overline{\mathrm{MS}}}_{\Gamma}$ to the $P^{z}b^{z}$-dependence of
the quasi beam functions, and that the results are stable under truncation of
the data in $P^{z}b^{z}$. The $P^{z}b^{z}$-range of the data presented here
is, however, not sufficient for the tails of the quasi beam functions at large
$|P^{z}b^{z}|$ to decay to plateaus consistent with zero, particularly at the
largest $b_{T}$ and smallest $P^{z}$ values used in this study. It is thus to
be expected that a DFT-based analysis has qualitative and quantitative
differences from the analytic model approach. These differences can be seen by
comparison of Fig. 26, which displays the results of a DFT-based analysis,
with Figs. 12 and 25, which show the results of the analysis of Sec. III.4. As
anticipated, the DFT approach yields numerical values which are significantly
different from those achieved by modeling rather than truncating the tails of
the quasi beam functions in $P^{z}b^{z}$, particularly for evaluations
including quasi beam functions computed with the smallest boost corresponding
to $n^{z}=3$. As a result, the values obtained with different choices of
$P^{z}$ are not consistent within uncertainties at intermediate values of $x$.
The differences are, naturally, less significant for results computed with the
largest choices of $P^{z}$, supporting the expectation that future studies
constraining larger values of $P^{z}b^{z}$ will achieve model-independent
results via the DFT approach.

##  Appendix D Power corrections and higher-twist effects

The estimator $\hat{\gamma}^{q}_{\zeta}(\mu,b_{T};P_{1}^{z},P_{2}^{z},x)$ (Eq.
(28)) coincides with the Collins-Soper kernel up to power corrections, such as
higher-twist corrections in the factorization formula for the quasi TMDPDF,
and discretization artifacts; in the absence of contamination from these
effects, $\hat{\gamma}^{q}_{\zeta}$ should be independent of $x$, $P_{1}^{z}$
and $P_{2}^{z}$. Clearly, the results shown in Figs. 12 and 25 indicate that
this contamination is not negligible relative to the uncertainties of this
calculation. As discussed in Sec. III.4, it is natural to attempt to model,
fit, and subtract this contamination in order to determine a best value for
the Collins-Soper kernel.

One possible approach is to consider a simple model of corrections to the
factorization formula Ebert  _et al._ (2019b); Ji  _et al._ (2019c) for the
quasi TMDPDF, for example through the inclusion of free parameters
$\lambda_{1}$ and $\lambda_{2}$ parameterizing power corrections as:

$\displaystyle\tilde{f}_{{\text{ns}}}^{\mathrm{TMD}}(x,\vec{b}_{T},$
$\displaystyle\mu,P^{z})=\left[C^{\mathrm{TMD}}_{\text{ns}}(\mu,xP^{z})+{\lambda_{1}\Lambda_{\text{QCD}}^{2}\over(xP^{z})^{2}}\right]$
$\displaystyle\times\exp\left[{1\over
2}\gamma^{q}_{\zeta}(\mu,b_{T})\ln{\frac{(2xP^{z})^{2}}{\zeta_{0}}}\right]g_{S}^{q}(b_{T},\mu)$
$\displaystyle\times
f_{{\text{ns}}}^{\mathrm{TMD}}(x,\vec{b}_{T},\mu,\zeta_{0})\left[1+{\lambda_{2}\over(xP^{z}b_{T})^{2}}\right]\,$
(31)

(here $g_{S}^{q}(b_{T},\mu)$ is defined as in Ref. Ebert  _et al._ (2019a) as
the mismatch between the lightlike and quasi soft factors). This form is
chosen since in the $\overline{\mathrm{MS}}$ scheme the higher-twist
corrections must appear in even powers, with a suppression through
$k_{T}^{2}/k_{z}^{2}$, which in Fourier space becomes $1/(xP^{z}b_{T})^{2}$.
At tree-level, the factor $\lambda_{2}$ is a constant. The $1/(xP^{z})^{2}$
power correction comes from the renormalon ambiguity in the perturbative
series for the matching coefficient.

The relationship between the quasi TMDPDF and the quasi beam functions (Eq.
(II)) then suggests a model of the Fourier transform of the quasi beam
functions, defined as

$\widetilde{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,x,b_{T},P^{z})\equiv\int
db^{z}\,e^{ib^{z}xP^{z}}P^{z}\hat{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}(\mu,b^{z},b_{T},P^{z}),$
(32)

of the form

$\displaystyle\widetilde{B}^{\text{fit}}(\mu,x,b_{T},P^{z})$
$\displaystyle=\left[C^{\mathrm{TMD}}_{\text{ns}}(\mu,xP^{z})+{\lambda_{1}\Lambda_{\rm
QCD}^{2}\over(xP^{z})^{2}}\right]\left({(2xP^{z})^{2}\over\zeta_{0}}\right)^{{1\over
2}\gamma_{\zeta}(\mu,b_{T})}$ $\displaystyle\quad\times
F(x,b_{T},\mu,\zeta_{0})\left[1+{\lambda_{2}\over(xP^{z}b_{T})^{2}}\right]\,,$
(33)

where $\gamma_{\zeta}$, $F$, $\lambda_{1}$ and $\lambda_{2}$ are free
parameters and $\zeta_{0}$ can be chosen freely. For each $x$, the model can
be fit to $\widetilde{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ (and
separately to $\widetilde{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$) for all
choices of $P^{z}$ and $b_{T}$ simultaneously, and the Collins-Soper kernel
extracted as the fit parameter $\gamma_{\zeta}$.

The results of this procedure yield results for the Collins-Soper kernel which
are not more consistent with a constant in $x$ than the results without the
correction applied. A global fit at discretized $x$ values is of poor quality,
with $\chi^{2}/\text{d.o.f.}\gtrsim 2$. There are a range of reasons that the
model form above may not a good description of the numerical data; for
example, the assumption that the $1/(xP^{z}b_{T})^{2}$ corrections are
proportional to the leading power contribution in Eq. (D), may not be a good
approximation. This approach is thus not taken in the main analysis presented
here, but may be worthwhile to consider for future studies with larger values
of $P^{z}$ where the power corrections will be suppressed relative to those in
the current work.

##  Appendix E NLO matching effects

A key difference between the approach followed to obtain the primary results
of this work and a number of the alternative approaches explored in Sec. III.4
is the inclusion of the perturbative matching coefficient
$C^{\mathrm{TMD}}_{\text{ns}}$ computed to NLO, instead of to LO, in the
calculation of the estimator
$\hat{\gamma}^{q}_{\zeta}(\mu,b_{T};P_{1}^{z},P_{2}^{z},x)$ via Eq. (28). To
illustrate the importance of this effect, Fig. 27 displays the relevant
contribution from the NLO matching coefficient to $\hat{\gamma}^{q}_{\zeta}$,
computed in Refs. Ebert  _et al._ (2019a, b), for each of the momentum
combinations used in the numerical study of this work. Clearly this
contribution, which is of the order of the Collins-Soper kernel itself for
$x<0.5$, is significant, and its inclusion affects not only the value but also
the $x$-dependence of the estimator $\hat{\gamma}^{q}_{\zeta}$.

Figure 18:  Examples of the bare quasi beam functions
$B^{\text{bare}}_{\gamma^{3}}$ determined as described in Sec. III.1.

Figure 19:  Examples of the modified $\overline{\mathrm{MS}}$-renormalized
quasi beam functions $B^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ determined as
described in Sec. III.3.

Figure 20:  Examples of the modified $\overline{\mathrm{MS}}$-renormalized
quasi beam functions $B^{\overline{\mathrm{MS}}}_{\gamma^{3}}$ determined as
described in Sec. III.3.

Figure 21: Fits to the $b^{z}$-dependent asymmetry in the modified
$\overline{\mathrm{MS}}$-renormalized quasi beam functions, as detailed in
Sec. III.3 (the asymmetry in the ratio of beam functions
$B^{\overline{\mathrm{MS}}}_{\gamma^{4}}$ with $b_{T}$=0.12 fm, $\eta=1.68$ fm
is provided in Fig. 7 of the main text).

Figure 22:  Example of the asymmetry-corrected modified
$\overline{\mathrm{MS}}$-renormalized quasi beam function
$B^{\overline{\mathrm{MS}};\text{corr}}_{\gamma^{3}}$, including the results
of weighted averages of this quantity over choices of $b_{T}^{R}$ and $\eta$
(as a function of $b^{z}$, $b_{T}$, and $P^{z}$), as described in the main
text. Fig. 9 displays the analogous figure for
$B^{\overline{\mathrm{MS}};\text{corr}}_{\gamma^{4}}$.

Figure 23: Examples of the averaged asymmetry-corrected modified
$\overline{\mathrm{MS}}$-renormalized quasi beam functions
$\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{4}}$, including fits by Eqs.
(26) and (27) to the real and imaginary parts, shown as shaded bands. Fig. 11
of the main text shows the example for $b_{T}=0.12$ fm.

Figure 24: As in Fig. 23, for
$\overline{B}^{\overline{\mathrm{MS}}}_{\gamma^{3}}$.

Figure 25:  $\hat{\gamma}^{q}_{\zeta}$ computed as defined in Eq. (28) for
momentum pairs $\\{P_{1}^{z},P_{2}^{z}\\}$, denoted by $P_{1}^{z}/P_{2}^{z}$
in the legend. The horizontal shaded band shows the total uncertainty of the
best result, and the corresponding $x$-window, determined as described in the
text. Fig. 12 of the main text shows the analogous results for $b_{T}=0.12$
fm.

Figure 26: $\hat{\gamma}^{q}_{\zeta}$ computed as defined in Eq. (28) with a
DFT of $\overline{B}^{\overline{\mathrm{MS}}}_{\Gamma}$ replacing the Fourier
transform of analytic fits $\hat{B}^{\overline{\mathrm{MS}}}_{\Gamma}$ to the
$P^{z}b^{z}$-dependence of the quasi beam functions, for momentum pairs
$\\{P_{1}^{z},P_{2}^{z}\\}$ denoted by $P_{1}^{z}/P_{2}^{z}$ in the legend.
Figure 27: Contribution of the NLO matching coefficient to the estimator
$\hat{\gamma}^{q}_{\zeta}$ for momentum pairs $\\{P_{1}^{z},P_{2}^{z}\\}$
denoted by $P_{1}^{z}/P_{2}^{z}$ in the legend. The contribution of this term
with LO matching, i.e., $C_{\text{ns}}(\mu,xP^{z})=1$, would be zero.

## References

  * Collins and Soper (1981) J. C. Collins and D. E. Soper, Nucl. Phys. B193, 381 (1981), [Erratum: Nucl. Phys.B213,545(1983)].
  * Collins and Soper (1982) J. C. Collins and D. E. Soper, Nucl. Phys. B197, 446 (1982).
  * Collins  _et al._ (1985) J. C. Collins, D. E. Soper, and G. F. Sterman, Nucl. Phys. B250, 199 (1985).
  * Scimemi and Vladimirov (2019) I. Scimemi and A. Vladimirov, (2019), arXiv:1912.06532 [hep-ph] .
  * Bacchetta  _et al._ (2019) A. Bacchetta, V. Bertone, C. Bissolotti, G. Bozzi, F. Delcarro, F. Piacenza, and M. Radici, (2019), arXiv:1912.07550 [hep-ph] .
  * Buffing  _et al._ (2018) M. G. A. Buffing, Z.-B. Kang, K. Lee, and X. Liu, (2018), arXiv:1812.07549 [hep-ph] .
  * Gutierrez-Reyes  _et al._ (2019) D. Gutierrez-Reyes, I. Scimemi, W. J. Waalewijn, and L. Zoppi, JHEP 10, 031 (2019), arXiv:1904.04259 [hep-ph] .
  * Gautheron  _et al._ (2010) F. Gautheron _et al._ (COMPASS), (2010).
  * Dudek  _et al._ (2012) J. Dudek _et al._ , Eur. Phys. J. A48, 187 (2012), arXiv:1208.1244 [hep-ex] .
  * Aschenauer  _et al._ (2015) E.-C. Aschenauer _et al._ , (2015), arXiv:1501.01220 [nucl-ex] .
  * Accardi  _et al._ (2016) A. Accardi _et al._ , Eur. Phys. J. A52, 268 (2016), arXiv:1212.1701 [nucl-ex] .
  * Abdul Khalek  _et al._ (2021) R. Abdul Khalek _et al._ , (2021), arXiv:2103.05419 [physics.ins-det] .
  * Musch  _et al._ (2011) B. U. Musch, P. Hagler, J. W. Negele, and A. Schafer, Phys. Rev. D83, 094507 (2011), arXiv:1011.1213 [hep-lat] .
  * Musch  _et al._ (2012) B. U. Musch, P. Hagler, M. Engelhardt, J. W. Negele, and A. Schafer, Phys. Rev. D85, 094510 (2012), arXiv:1111.4249 [hep-lat] .
  * Engelhardt  _et al._ (2016) M. Engelhardt, P. Hägler, B. Musch, J. Negele, and A. Schäfer, Phys. Rev. D93, 054501 (2016), arXiv:1506.07826 [hep-lat] .
  * Yoon  _et al._ (2015) B. Yoon, T. Bhattacharya, M. Engelhardt, J. Green, R. Gupta, P. Hägler, B. Musch, J. Negele, A. Pochinsky, and S. Syritsyn, in  _Proceedings, 33rd International Symposium on Lattice Field Theory (Lattice 2015): Kobe, Japan, July 14-18, 2015_ , SISSA (SISSA, 2015) arXiv:1601.05717 [hep-lat] .
  * Yoon  _et al._ (2017) B. Yoon, M. Engelhardt, R. Gupta, T. Bhattacharya, J. R. Green, B. U. Musch, J. W. Negele, A. V. Pochinsky, A. Schäfer, and S. N. Syritsyn, Phys. Rev. D96, 094508 (2017), arXiv:1706.03406 [hep-lat] .
  * Shanahan  _et al._ (2019) P. Shanahan, M. Wagman, and Y. Zhao, (2019), arXiv:1911.00800 [hep-lat] .
  * Shanahan  _et al._ (2020) P. Shanahan, M. Wagman, and Y. Zhao, Phys. Rev. D 102, 014511 (2020), arXiv:2003.06063 [hep-lat] .
  * Zhang  _et al._ (2020) Q.-A. Zhang _et al._ (Lattice Parton), Phys. Rev. Lett. 125, 192001 (2020), arXiv:2005.14572 [hep-lat] .
  * Schlemmer  _et al._ (2021) M. Schlemmer, A. Vladimirov, C. Zimmermann, M. Engelhardt, and A. Schäfer, (2021), arXiv:2103.16991 [hep-lat] .
  * Li  _et al._ (2021) Y. Li _et al._ , (2021), arXiv:2106.13027 [hep-lat] .
  * Ji (2013) X. Ji, Phys. Rev. Lett. 110, 262002 (2013), arXiv:1305.1539 [hep-ph] .
  * Ji (2014) X. Ji, Sci. China Phys. Mech. Astron. 57, 1407 (2014), arXiv:1404.6680 [hep-ph] .
  * Ji  _et al._ (2020) X. Ji, Y. Liu, Y.-S. Liu, J.-H. Zhang, and Y. Zhao, (2020), arXiv:2004.03543 [hep-ph] .
  * Ebert  _et al._ (2019a) M. A. Ebert, I. W. Stewart, and Y. Zhao, Phys. Rev. D99, 034505 (2019a), arXiv:1811.00026 [hep-ph] .
  * Ebert  _et al._ (2019b) M. A. Ebert, I. W. Stewart, and Y. Zhao, (2019b), arXiv:1901.03685 [hep-ph] .
  * Ebert  _et al._ (2020) M. A. Ebert, I. W. Stewart, and Y. Zhao, JHEP 03, 099 (2020), arXiv:1910.08569 [hep-ph] .
  * Ji  _et al._ (2015) X. Ji, P. Sun, X. Xiong, and F. Yuan, Phys. Rev. D91, 074009 (2015), arXiv:1405.7640 [hep-ph] .
  * Ji  _et al._ (2019a) X. Ji, L.-C. Jin, F. Yuan, J.-H. Zhang, and Y. Zhao, Phys. Rev. D99, 114006 (2019a), arXiv:1801.05930 [hep-ph] .
  * Ji  _et al._ (2019b) X. Ji, Y. Liu, and Y.-S. Liu, (2019b), arXiv:1910.11415 [hep-ph] .
  * Ji  _et al._ (2019c) X. Ji, Y. Liu, and Y.-S. Liu, (2019c), arXiv:1911.03840 [hep-ph] .
  * Vladimirov and Schäfer (2020) A. A. Vladimirov and A. Schäfer, (2020), arXiv:2002.07527 [hep-ph] .
  * Constantinou  _et al._ (2019) M. Constantinou, H. Panagopoulos, and G. Spanoudes, Phys. Rev. D99, 074508 (2019), arXiv:1901.03862 [hep-lat] .
  * Green  _et al._ (2020) J. R. Green, K. Jansen, and F. Steffens, (2020), arXiv:2002.09408 [hep-lat] .
  * Ji  _et al._ (2021) Y. Ji, J.-H. Zhang, S. Zhao, and R. Zhu, (2021), arXiv:2104.13345 [hep-ph] .
  * Bazavov  _et al._ (2013) A. Bazavov _et al._ (MILC), Phys. Rev. D 87, 054505 (2013), arXiv:1212.4768 [hep-lat] .
  * Lüscher (2010) M. Lüscher, JHEP 08, 071 (2010), [Erratum: JHEP03,092(2014)], arXiv:1006.4518 [hep-lat] .
  * Bali  _et al._ (2016) G. S. Bali, B. Lang, B. U. Musch, and A. Schäfer, Phys. Rev. D93, 094515 (2016), arXiv:1602.05525 [hep-lat] .
  * Aoki  _et al._ (2020) S. Aoki _et al._ (Flavour Lattice Averaging Group), Eur. Phys. J. C 80, 113 (2020), arXiv:1902.08191 [hep-lat] .
  * Akaike (1974) H. Akaike, IEEE Transactions on Automatic Control 19, 716 (1974).
  * Li and Zhu (2017) Y. Li and H. X. Zhu, Phys. Rev. Lett. 118, 022004 (2017), arXiv:1604.01404 [hep-ph] .
  * Henn  _et al._ (2020) J. M. Henn, G. P. Korchemsky, and B. Mistlberger, JHEP 04, 018 (2020), arXiv:1911.10174 [hep-th] .
  * Edwards and Joo (2005) R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD),  _Lattice field theory. Proceedings, 22nd International Symposium, Lattice 2004, Batavia, USA, June 21-26, 2004_ , Nucl. Phys. Proc. Suppl. 140, 832 (2005), [,832(2004)], arXiv:hep-lat/0409003 [hep-lat] .
  * (45) A. Pochinsky, Qlua. https://usqcd.lns.mit.edu/qlua.
  * Clark  _et al._ (2010) M. Clark, R. Babich, K. Barros, R. Brower, and C. Rebbi, Comput. Phys. Commun. 181, 1517 (2010), arXiv:0911.3191 [hep-lat] .
  * Babich  _et al._ (2011) R. Babich, M. Clark, B. Joo, G. Shi, R. Brower, and S. Gottlieb, in  _SC11 International Conference for High Performance Computing, Networking, Storage and Analysis_ (2011) arXiv:1109.2935 [hep-lat] .
  * Clark  _et al._ (2016) M. A. Clark, B. Joo, A. Strelchenko, M. Cheng, A. Gambhir, and R. Brower, (2016), arXiv:1612.07873 [hep-lat] .
  * Winter  _et al._ (2014) F. T. Winter, M. A. Clark, R. G. Edwards, and B. Joó, in  _2014 IEEE 28th International Parallel and Distributed Processing Symposium_ (2014) pp. 1073–1082.
  * Joó  _et al._ (2016) B. Joó, D. D. Kalamkar, T. Kurth, K. Vaidyanathan, and A. Walden, in  _High Performance Computing_ , edited by M. Taufer, B. Mohr, and J. M. Kunkel (Springer International Publishing, Cham, 2016) pp. 415–427.
  * Harris  _et al._ (2020) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature 585, 357 (2020).
  * Bezanson  _et al._ (2017) J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM Review 59, 65 (2017).
  * Mogensen and Riseth (2018) P. K. Mogensen and A. N. Riseth, Journal of Open Source Software 3, 615 (2018).
  * (54) W. R. Inc., “Mathematica, Version 12.2,” Champaign, IL, 2020\.