# Non-perturbative renormalization of asymmetric staple-shaped
operators in twisted mass lattice QCD

Constantia Alexandrou<EMAIL_ADDRESS>Department of Physics, University of
Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus  Computation-based Science and
Technology Research Center, The Cyprus Institute, 20 Kavafi Str., Nicosia
2121, Cyprus  Simone Bacchio  Computation-based Science and Technology
Research Center, The Cyprus Institute, 20 Kavafi Str., Nicosia 2121, Cyprus
Krzysztof Cichy  Faculty of Physics, Adam Mickiewicz University, ul.
Uniwersytetu Poznańskiego 2, 61-614 Poznań, Poland  Martha Constantinou
Temple University,1925 N. 12th Street, Philadelphia, PA 19122-1801, USA  Xu
Feng  School of Physics and State Key Laboratory of Nuclear Physics and
Technology, Peking University, Beijing 100871, China  Collaborative Innovation
Center of Quantum Matter, Beijing 100871, China  Center for High Energy
Physics, Peking University, Beijing 100871, China  Karl Jansen  NIC, DESY,
Platanenallee 6, D-15738 Zeuthen, Germany  Chuan Liu  School of Physics and
State Key Laboratory of Nuclear Physics and Technology, Peking University,
Beijing 100871, China  Collaborative Innovation Center of Quantum Matter,
Beijing 100871, China  Center for High Energy Physics, Peking University,
Beijing 100871, China  Aniket Sen<EMAIL_ADDRESS>Helmholtz Institut
für Strahlen- und Kernphysik, Rheinische Friedrich-Wilhelms-Universität Bonn,
Nussallee 14-16, 53115 Bonn  Bethe Center for Theoretical Physics, Rheinische
Friedrich-Wilhelms-Universität Bonn, Nussallee 12, 53115 Bonn, Germany
Gregoris Spanoudes  Department of Physics, University of Cyprus, P.O. Box
20537, 1678 Nicosia, Cyprus  Computation-based Science and Technology Research
Center, The Cyprus Institute, 20 Kavafi Str., Nicosia 2121, Cyprus  Fernanda
Steffens  Helmholtz Institut für Strahlen- und Kernphysik, Rheinische
Friedrich-Wilhelms-Universität Bonn, Nussallee 14-16, 53115 Bonn  Bethe Center
for Theoretical Physics, Rheinische Friedrich-Wilhelms-Universität Bonn,
Nussallee 12, 53115 Bonn, Germany  Jacopo Tarello  Department of Physics,
University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus  Computation-based
Science and Technology Research Center, The Cyprus Institute, 20 Kavafi Str.,
Nicosia 2121, Cyprus

###### Abstract

Staple-shaped Wilson line operators are necessary for the study of transverse
momentum-dependent parton distribution functions (TMDPDFs) in lattice QCD and
beyond. In this work, we study the renormalization of such operators in the
general case of an asymmetric staple. We analyze the mixing pattern of these
operators using their symmetry properties, where we find that the possible
mixing is restricted within groups of four operators. We then present
numerical results using the regularization independent scheme, RI/MOM, to
study the importance of mixing using one operator in particular, the
$\gamma_{0}$ operator. Based on these results, we consider the short distance
ratio (SDR) scheme, which is desirable in the absence of mixing. Finally, we
investigate a variant of the RI/MOM scheme, where the renormalization factors
are computed at short distances.

##  I Introduction

Collinear parton distribution functions (PDFs) probe the hadron structure from
the perspective of the spin and longitudinal momentum distributions of the
quarks and gluons that make up the hadron. The vast amount of work to
determine these distributions over the last five decades, both theoretically
and experimentally, has enormously expanded our view of the hadron structure,
the proton in particular Amoroso et al. (2022). And yet, these developments
are mostly limited to probing the one-dimensional structure of the proton. In
order to have a wider understanding of the proton structure, we also need to
understand how the momentum and the spin are distributed in the transverse
plane. For that, we need to measure and compute the generalized parton
distributions and the transverse-momentum-dependent PDFs (TMDPDFs), the latter
being the main focus of this present work. Although there has been an effort
to obtain TMDPDFs from phenomenological fits to experimental data Bacchetta et
al. (2017); Bertone et al. (2019); Scimemi and Vladimirov (2020); Bury et al.
(2022); Bacchetta et al. (2022); Barry et al. (2023), their accuracy is still
far from being at the same level as the collinear PDFs. This status will
change in the coming years with new data coming from Jefferson Lab Dudek et
al. (2012) and from the future Electron-Ion-Collider to be built at Brookhaven
National Lab Abdul Khalek et al. (2021). It is, thus, of great importance to
extract TMDPDFs from first principles calculations, namely, lattice QCD.

In the last 10 years, there has been an enormous advance in computing the
collinear PDFs using lattice QCD Ji (2013); Radyushkin (2017); Cichy and
Constantinou (2019); Radyushkin (2020); Constantinou (2021); Ji et al.
(2021a); Cichy (2022a, b); Constantinou et al. (2022). By comparison, the
computation of TMDPDFs is still in its infancy Ji et al. (2015, 2019, 2020a);
Ebert et al. (2022), although progress has been made in the computation of key
elements required to build the TMDPDFS, like the soft functions Zhang et al.
(2020); Li et al. (2022) and the Collins-Soper kernel Ebert et al. (2019a,
2020); Shanahan et al. (2020). Such computations have been restricted to
ratios of TMDPDFs Hagler et al. (2009); Musch et al. (2012); Yoon et al.
(2015, 2017), and only recently, a first full lattice QCD calculation of the
TMDPDFs themselves has been presented He et al. (2022). One fundamental
difficulty in these calculations is to have control over the renormalization
procedure, which is more involved than in the case of collinear PDFs. In
particular, the staple-shaped gauge link that enters in the evaluation of
TMDPDFs has three types of divergences: i) linear divergence coming from the
Wilson line, which connects the quark fields and which depends on the length
of the staple-shaped link; ii) logarithmic divergences coming form the
endpoints of the staple link, similarly to the case of straight gauge links;
and iii) logarithmic divergences coming from the presence of cusps in the
staple. In addition, in the limit of infinite-length staple, $L$, (which is
the case of interest), pinch-pole singularities arise as positive powers of
$L$, coming from the gluon exchange between the transverse segments of the
staple. Moreover, as in the case of straight gauge links, staple-shaped
operators of different Dirac structures can and will mix on the lattice among
certain groups (when chiral symmetry is broken), as dictated by discrete
symmetries. However, in this case, the mixing pattern of the operators
employed in lattice regularization can be significantly more involved than in
the case of the straight Wilson line.

A first study within lattice perturbation theory to one-loop for the case of
the symmetric staple Constantinou et al. (2019) showed mixing between certain
pairs of Dirac structures. The mixing depended solely on the direction of the
staple link entering the end points of the staple regardless of the shape of
the staple. This implies that the same mixing pattern occurs also for
asymmetric staple shapes. In Ref. Ji et al. (2021b), the mixing pattern of
these staple-shaped operators111For conciseness, we will refer to staple-
shaped Wilson quark bilinear operators as staple-shaped operators. has been
studied using symmetry considerations. It was found that more mixing is
present than observed in Ref. Constantinou et al. (2019). This demonstrates
that one-loop perturbation theory cannot fully reveal the mixing of the staple
operators, unlike the case of straight Wilson-line operators Constantinou and
Panagopoulos (2017), and a higher-loop computation is needed to confirm the
additional mixing patterns found by symmetry arguments. The authors of Refs.
Shanahan et al. (2020); Zhang et al. (2020, 2022) examined 16 independent
projectors in the RI-type prescription for each Dirac structure. While Ref.
Shanahan et al. (2020) identified several finite contributions in the
asymmetric staple-shaped operator, Ref. Zhang et al. (2020, 2022) found
negligible contributions, at least at small transverse separations, by setting
specific momentum components to be zero. However, we emphasize that this does
not imply nonzero mixing, as these finite contributions may appear in the
Green’s functions defined in the continuum and, thus, vanish in the
renormalization prescription that requires the difference between lattice and
continuum expressions. Additional finite contributions can also arise due to
mixing with higher dimensional operators of different Dirac structures which
vanish in the continuum limit. Moreover, these contributions are of higher
twist and are disregarded in the present calculation. In this work, we use
symmetry arguments to restrict the set of operators that are allowed to mix.

The common conclusion of these studies is that there are remaining non-
perturbative divergences in the RI/MOM scheme that make it hard to renormalize
at large transverse separations. An alternative approach that one can adopt is
the so-called ratio scheme as proposed for the quasi-PDF case Radyushkin
(2017); Orginos et al. (2017); Braun et al. (2019). In this approach, one
subtracts the Ultra-Violet (UV) divergences by taking the ratio with a
suitable object at a fixed short distance, where perturbation theory applies.
Many choices of suitable objects for the ratio have been proposed in Ref. Ji
et al. (2021c). The authors of Refs. Zhang et al. (2022); He et al. (2022)
choose the same matrix element but at zero longitudinal momentum for taking
the ratio and name this scheme as short distance ratio (SDR). As for the
remaining divergences associated with the asymmetric staple-shaped link, one
can cancel them by taking an appropriate ratio with the vacuum expectation
value of a rectangular Wilson loop Ji et al. (2020b). However, we would like
to stress that the SDR scheme is valid when operator mixing is absent or
negligible. We check that in the case under study the mixing is indeed
negligible and thus, one can employ the ratio scheme. Another option, as
proposed in Ref. Ji et al. (2021b), is to employ RI/MOM in the spirit of the
SDR scheme by fixing the dimensions of the staple at a short perturbative
range (we will call this scheme RI-short).

The paper is organized as follows. Section II discusses how TMDs can be
accessed on a Euclidean lattice using the quasi-distribution approach. In
Section III, we study operator mixing using symmetry arguments. Section IV
presents our lattice setup. The size of the mixing is examined in Section V.
In Sections V.1, V.2 and V.3 we employ RI/MOM, SDR and RI-short schemes,
respectively, by following the procedures described in Ref. Zhang et al.
(2022); He et al. (2022); Ji et al. (2021b). In Section VI, we present the
renormalized beam functions in the three renormalization schemes discussed
here, and then show the corresponding results in the $\overline{\text{MS}}$
scheme. The conclusions are presented in Section VII.

##  II Quasi-TMDPDFs

A first-principle computation of TMDPDFs in the context of LaMET is more
involved than collinear PDFs. The main obstacle is to correctly subtract the
so-called rapidity divergences, which are associated with the soft gluon
radiation of the infinitely long Wilson lines present in the staple-shaped
operator. This subtraction is made through the use of a soft function Collins
(2011); Echevarria et al. (2016). Different ways how to perform this
subtraction in the context of LaMET can be found in Ref. Ebert et al. (2019b).
In general, the TMD soft function involves two opposite light-like directions,
and this makes a lattice calculation significantly challenging. In Refs. Ji et
al. (2020a, b), the authors define a rapidity-independent reduced soft
function $S_{r}$ and show that it can be extracted from a form factor and the
quasi-TMD wave function of a light meson. Using $S_{r}$, the rapidity scheme-
independent TMDPDF $f^{TMD}(x,b,\mu,\zeta)$ can be written as Ji et al.
(2020b, 2021a)

$\begin{split}&f^{TMD}(x,b,\mu,\zeta)=H\left(\frac{\zeta_{z}}{\mu^{2}}\right)\,e^{-\ln{\left(\frac{\zeta_{z}}{\zeta}\right)}K(b,\mu)}\tilde{f}(x,b,\mu,\zeta_{z})\,S_{r}^{\frac{1}{2}}(b,\mu)+\dots\end{split}$
(1)

where $\tilde{f}(x,b,\mu,\zeta_{z})$ is the so-called quasi-TMDPDF, $x$ is the
longitudinal momentum fraction, $b$ the transverse separation. The scale $\mu$
defines the renormalization scale and $\zeta_{z}=(2xP^{z})^{2}$ is the
Collins-Soper scale of the quasi-TMDPDF, with $\zeta$ being the scale for the
light-cone correlation. The factor $H\left(\frac{\zeta_{z}}{\mu^{2}}\right)$
is the perturbative matching kernel that perturbatively connects the TMDPDFs
to the quasi-TMDPDFs, and $K(b,\mu)$ is the Collins-Soper kernel. The
renormalization scheme in Eq. (1) is left unspecified; usually, it is computed
in the $\overline{\text{MS}}$ scheme. The quasi-TMDPDF on the lattice is
defined as Ebert et al. (2019b); Ji et al. (2019)

$\tilde{f}(x,b,\mu,\zeta_{z})=\lim_{L\to\infty}\int\frac{dz}{2\pi}e^{-iz\zeta_{z}}\frac{P^{z}}{E_{P}}B(z,b,L,P^{z},\mu).$
(2)

Here, $B(z,b,P^{z})$ is known as the quasi-beam function, and it is defined in
terms of a bare matrix element with a non-local quark bilinear operator as

$\displaystyle B_{0,\Gamma}(z,b,L,P^{z})$ $\displaystyle=$
$\displaystyle\langle N(P^{z})|\mathcal{O}^{\Gamma}(z,b,L)|N(P^{z})\rangle$
(3) $\displaystyle=$ $\displaystyle\langle
N(P^{z})|\bar{\psi}(b+z)\Gamma\mathcal{W}(b+z;L)\psi(0)|N(P^{z})\rangle,$ (4)

where $N(P^{z})$ is a nucleon state with momentum boost of $(0,0,P^{z})$ and
$\psi(x)$ is the standard up and down quark doublet. For the unpolarized
TMDPDF, $\Gamma$ can be either $\gamma_{0}$ or $\gamma_{3}$ or in general a
combination of the two. In this work, we show results for $\Gamma=\gamma_{0}$,
since we observe a better signal for this operator compared to
$\Gamma=\gamma_{3}$. $\gamma_{0}$ was also used for the quasi-PDF case,
strongly due to this operator being free from any mixing, while $\gamma_{3}$
mixed with $\mathbf{1}$. In the case of the staple, however, no such advantage
exists, as will be shown in Section III. In general, the asymmetric staple-
shaped Wilson, $\mathcal{W}(b+z;L)$, is given by:

$\begin{split}\mathcal{W}(b+z;L)&=W_{z}(\vec{x};L)W_{\perp}(\vec{x}+L\hat{z};b)W_{z}(\vec{x}+L\hat{z}+b\hat{y};-L+z)\,,\\\
W_{z}(\vec{x};L)&={\cal P}\exp\left[-ig\int_{0}^{L}d\lambda\hat{z}\cdot
A(\vec{x}+\hat{z}\lambda)\right]\,.\end{split}$ (5)

with $L$ the length of the symmetric part, and $x$ an arbitrary point, while
$W_{z}$ defines a Wilson line along the boost direction and $W_{\perp}$
defines one along the transverse direction. The shape of the asymmetric staple
is shown in Fig. 1.

$\vec{x}$$\vec{x}+L\hat{z}$$\vec{x}+L\hat{z}+b\hat{y}$$\vec{x}+b\hat{y}+z\hat{z}$
Figure 1: The shape of the asymmetric staple defined in the operator of the
quasi-beam function

##  III Operator mixing through symmetry

In order to use symmetry arguments to determine which operators mix, we first
organize the staple-shaped operator using all possible ways to connect the
fermion fields. We show in Fig. 2 two of the possible ways to connect the
quark and antiquark fields with coordinates $(\hat{y},\hat{z})$. In this
notation, $b$ is always in the $\hat{y}$ direction, while $L$ and $z$ are in
the $\hat{z}$ direction. We define the operator shown on the left of Fig. 2 as
$O^{++}(\Gamma)$, where the two plus signs define the positive direction of
the Wilson line along transverse ($\hat{y}$) and longitudinal ($\hat{z}$)
axes. Under this notation, the operator shown on the right of Fig. 2 is
$O^{+-}(\Gamma)$. We have two other possibilities, $O^{-+}(\Gamma)$ and
$O^{--}(\Gamma)$, along with the corresponding charge-conjugated operators,
leading to a total of eight operators. The charge conjugated operators $O_{C}$
are given by the same geometry as that of the original staple, only with the
path going in the opposite direction. Thus, $O_{C}^{++}(\Gamma)$ is also given
by the left staple of Fig. 2, only the arrows are in the opposite direction,
and similarly for the remaining three operators. In the case of the symmetric
staple ($z=0$), the charge-conjugated operators are redundant.

$q(0,0)$$(0,L+z)$$(b,L+z)$$\overline{q}(b,z)$$q(0,0)$$(0,-L-z)$$(b,-L-z)$$\overline{q}(b,-z)$
Figure 2: The shape of the asymmetric staple defined in the operator of the
quasi-beam function

We analyze the symmetry properties using generalized parity and time reversal
with discrete flavor rotation, and charge conjugation for the fermion fields
in the twisted mass basis, $\chi(x_{\alpha},\textbf{x})$, and the gauge link,
$U(x_{\alpha},\textbf{x};\alpha)$, in some direction $\alpha$. The generalized
parity in the $\alpha$-direction combined with discrete flavor rotation
Shindler (2008) is then given by

$\mathcal{P}_{F\alpha}^{1,2}:\begin{cases}U(x_{\alpha},\textbf{x};\alpha)\to
U(x_{\alpha},-\textbf{x};\alpha),\\\
U(x_{\alpha},\textbf{x};\beta\neq\alpha)\to
U^{-1}(x_{\alpha},-\textbf{x}-a\hat{\beta};\beta),\\\
\chi(x_{\alpha},\textbf{x})\to
i\gamma_{\alpha}\tau_{1,2}\chi_{l}(x_{\alpha},-\textbf{x}),\\\
\bar{\chi}(x_{\alpha},\textbf{x})\to-i\bar{\chi}(x_{\alpha},-\textbf{x})\tau_{1,2}\gamma_{\alpha},\end{cases}$
(6)

where the ”standard” parity is the one with $\alpha=0$, and the $3$-vector x
is what remains from the $4$-vector $x$ after removing $x_{\alpha}$.
$\tau_{1,2}$ are the Pauli spin matrices in flavor space. Similarly, the
generalized time reversal combined with discrete flavor rotation is given by

$\mathcal{T}_{F\alpha}^{1,2}:\begin{cases}U(x_{\alpha},\textbf{x};\alpha)\to
U^{-1}(-x_{\alpha}-a,\textbf{x};\alpha),\\\
U(x_{\alpha},\textbf{x};\beta\neq\alpha)\to
U(-x_{\alpha},\textbf{x};\beta),\\\ \chi(x_{\alpha},\textbf{x}\to
i\gamma_{\alpha}\gamma_{5}\tau_{1,2}\chi_{l}(-x_{\alpha},\textbf{x}),\\\
\bar{\chi}(x_{\alpha},\textbf{x})\to-i\bar{\chi}(-x_{\alpha},\textbf{x})\tau_{1,2}\gamma_{5}\gamma_{\alpha}.\end{cases}$
(7)

Finally, the charge conjugation transformation is given by

$\mathcal{C}:\begin{cases}U(x)\to(U^{\dagger}(x))^{T},\\\ \chi(x)\to
C^{-1}\bar{\chi}(x)^{T},\\\ \bar{\chi}(x)\to-\chi(x)^{T}C.\end{cases}$ (8)

To study the symmetry properties, we need to consider all eight different
staple-shaped operators. To summarize the effect of the transformations, we
use the notation $(ijkl)_{c}$. Here $i,j,k,l$ denote the signs of
$O^{++},O^{--},O^{+-},O^{-+}$ and $c=\pm$ denotes the sign of the charge
conjugated operators ($c=+$ if it stays the same, $c=-$ if it is opposite).
The relevant combinations that have definite symmetry properties are:

$\displaystyle(++++)_{c}$ $\displaystyle\equiv(----)_{c},$
$\displaystyle(+-+-)_{c}$ $\displaystyle\equiv(-+-+)_{c},$
$\displaystyle(++--)_{c}$ $\displaystyle\equiv(--++)_{c},$
$\displaystyle(+--+)_{c}$ $\displaystyle\equiv(-++-)_{c}.$

As an example, we look at the symmetry properties of $\gamma_{0}$ and
$\gamma_{0}\gamma_{3}$ given in Tables 1 and 2, respectively. The Dirac
structure $\Gamma^{{}^{\prime}}$ is the transformation applied to the
insertion Dirac structure $\Gamma$, and we restrict ourselves here to the
flavor non-singlet case, $\tau_{3}$. The combination of operators that mix
are:

$\displaystyle(++++)_{c}\;\;$ $\displaystyle\text{with}\;\;(+--+)_{c},$
$\displaystyle(+-+-)_{c}\;\;$ $\displaystyle\text{with}\;\;(++--)_{c},$
$\displaystyle(++--)_{c}\;\;$ $\displaystyle\text{with}\;\;(+-+-)_{c},$
$\displaystyle(+--+)_{c}\;\;$ $\displaystyle\text{with}\;\;(++++)_{c}.$

Although we show the resulting mixing pattern only for the $\gamma_{0}$ and
$\gamma_{0}\gamma_{3}$ cases, this mixing holds true for all pairs
$(\Gamma,\Gamma\gamma_{3})$. The symmetry properties for all $\Gamma$’s are
summarized in Table 4 in Appendix A. Based on this, we conclude that the
possible mixing is between $\Gamma$ and
$\left\\{\Gamma\gamma_{2},\Gamma\gamma_{3},\Gamma\gamma_{2}\gamma_{3}\right\\}$.
Since, in this work, we only consider $\Gamma=\gamma_{0}$, the relevant mixing
for us is between $\gamma_{0}$ and
$\left\\{\gamma_{0}\gamma_{2},\gamma_{0}\gamma_{3},\gamma_{5}\gamma_{1}\right\\}$.

We point out that different mixing patterns have been considered in recent
studies of other groups, including mixing among all 16 operators of different
Dirac structures Shanahan et al. (2020), or at least mixing in pairs
($\Gamma,\Gamma\gamma_{3}$) Ji et al. (2021b). In our work, we choose to
consider the minimal set of staple-shaped operators (of the same dimension)
that are allowed to mix by the above-mentioned C, P, and T symmetries. We have
not considered mixing with higher-dimensional operators allowed by Lorentz
symmetry (see, e.g., Ref. Chen et al. (2019) for the straight Wilson-line
case), since it is power suppressed and not relevant when one takes the
continuum limit $a\rightarrow 0$. This can explain the nonzero (but small)
mixings among certain Dirac structures observed in Ref. Shanahan et al.
(2020). Also, in contrast to Ref. Ji et al. (2021b), in our analysis, the
staple line has been chosen to be in a 2-dimensional (D), and not 3D, plane in
Euclidean space formed by the transverse ($\hat{y}$) and longitudinal
($\hat{z}$) directions. In this respect, we end up with a basis of 8 (instead
of 16) operators, which are eigenstates of C, P and T transformations. We also
note that calculations in one-loop lattice perturbation theory Constantinou et
al. (2019) show a smaller mixing pattern compared to the symmetries; however,
this cannot guarantee a reduced mixing in higher loops.

| $\Gamma^{\prime}$ |  $(++++)_{c}$ $(----)_{c}$ |  $(+-+-)_{c}$ $(-+-+)_{c}$ |  $(++--)_{c}$ $(--++)_{c}$ |  $(+--+)_{c}$ $(-++-)_{c}$
---|---|---|---|---|---
$\mathcal{P}_{F0}^{1,2}$ | $\gamma_{0}$ | - | + | - | +
$\mathcal{P}_{F1}^{1,2}$ | $-\gamma_{0}$ | + | - | + | -
$\mathcal{P}_{F2}^{1,2}$ | $-\gamma_{0}$ | + | + | - | -
$\mathcal{P}_{F3}^{1,2}$ | $-\gamma_{0}$ | + | - | - | +
$\mathcal{T}_{F0}^{1,2}$ | $-\gamma_{0}$ | + | + | + | +
$\mathcal{T}_{F1}^{1,2}$ | $\gamma_{0}$ | - | - | - | -
$\mathcal{T}_{F2}^{1,2}$ | $\gamma_{0}$ | - | + | + | -
$\mathcal{T}_{F3}^{1,2}$ | $\gamma_{0}$ | - | - | + | +
$C$ | $\gamma_{0}$ | $c$ | $c$ | $c$ | $c$
Table 1: Symmetry properties of $\gamma_{0}$. | $\Gamma^{\prime}$ |  $(++++)_{c}$ $(----)_{c}$ |  $(+-+-)_{c}$ $(-+-+)_{c}$ |  $(++--)_{c}$ $(--++)_{c}$ |  $(+--+)_{c}$ $(-++-)_{c}$
---|---|---|---|---|---
$\mathcal{P}_{F0}^{1,2}$ | $-\gamma_{0}\gamma_{3}$ | + | - | + | -
$\mathcal{P}_{F1}^{1,2}$ | $\gamma_{0}\gamma_{3}$ | - | + | - | +
$\mathcal{P}_{F2}^{1,2}$ | $\gamma_{0}\gamma_{3}$ | - | - | + | +
$\mathcal{P}_{F3}^{1,2}$ | $-\gamma_{0}\gamma_{3}$ | + | - | - | +
$\mathcal{T}_{F0}^{1,2}$ | $-\gamma_{0}\gamma_{3}$ | + | + | + | +
$\mathcal{T}_{F1}^{1,2}$ | $\gamma_{0}\gamma_{3}$ | - | - | - | -
$\mathcal{T}_{F2}^{1,2}$ | $\gamma_{0}\gamma_{3}$ | - | + | + | -
$\mathcal{T}_{F3}^{1,2}$ | $-\gamma_{0}\gamma_{3}$ | + | + | - | -
$C$ | $-\gamma_{0}\gamma_{3}$ | $c$ | $c$ | $c$ | $c$
Table 2: Symmetry properties of $\gamma_{0}\gamma_{3}$.

##  IV Lattice Setup

For the lattice simulation, we use an $N_{f}=2+1+1$ clover-improved twisted
mass fermion ensemble of size $24^{3}\times 48$, produced by the Extended
Twisted Mass Collaboration (ETMC) Alexandrou et al. (2018). Our study is done
using $N_{conf}$ configurations with $N_{src}$ source positions for each
configuration. To increase statistics, the boost is taken in all three
directions, and both positive and negative. For each such direction of boost,
the staple is then constructed in both the remaining transverse directions.
This gives us 12 measurements ($6$ boost directions $\times$ $2$ transverse
directions) for each source position. The details of the lattice simulation
are summarized in Table 3.

Lattice size |  $a$ [fm] | $a\mu_{l}$ |  $m_{\pi}$ [MeV] | $N_{conf}$ | $N_{src}$ | $N_{meas}$
---|---|---|---|---|---|---
$24^{3}\times 48$ | $0.093$ | $0.00530$ | $350$ | $100$ | 8 | 9600
Table 3: We give the parameters of the lattice ensemble and measurements used
in the calculation. $a$ is the lattice spacing, $\mu_{l}$ is the bare twisted
light quark mass, $m_{\pi}$ the pion mass, $N_{conf}$ the number of
configurations, $N_{src}$ the number of source positions and $N_{meas}$ the
total number of measurements.

The bare matrix element for the quasi-beam function is calculated through a
ratio of a 3-point to a 2-point function,

$\begin{split}B^{\Gamma}(z,b,L,P^{z})&=\frac{\langle
C^{3pt}_{\Gamma}(z,b,L,P^{z};t_{s},\tau)\rangle}{\langle
C^{2pt}(P^{z};t_{s})\rangle}\\\ &=\frac{{\cal
P}\sum_{\mathbf{x}}e^{-i\mathbf{P}\cdot\mathbf{x}}\langle
0|N(\mathbf{x},t_{s})\mathcal{O}^{\Gamma}(z,b,L;\tau)\bar{N}(\mathbf{0},0)|0\rangle}{{\cal
P}\sum_{\mathbf{x}}e^{-i\mathbf{P}\cdot\mathbf{x}}\langle
0|N(\mathbf{x},t_{s})\bar{N}(\mathbf{0},0)|0\rangle},\end{split}$ (9)

where $\tau$ is the insertion time of the operator $\mathcal{O}^{\Gamma}$ and
$t_{s}$ defines the source-sink time separation. ${\cal P}$ is the parity
projector. In this work, we show results for a single source-sink separation
of $t_{s}=10a$ and a longitudinal momentum of $P^{z}=6\pi/24a\sim
1.7\,\mathrm{GeV}$.
The 3-point function is constructed for the isovector combination $u-d$, by
inserting a $\tau_{3}$ in flavor space. This choice ensures the elimination of
the disconnected contributions and only connected diagrams need to be
calculated.
Momentum smearing Bali et al. (2016) is applied to improve the signal for
large boosts. It is observed that applying stout smearing to the gauge links
used in construction of the staple also reduces the statistical errors. Here,
we have applied 5 steps of stout smearing to the staple-shaped Wilson line.
In Fig. 3 we show the bare matrix elements as a function of $z$ at a fixed
$b/a=1$ for different values of $L$. We notice that the real part of the
matrix elements are symmetric in the simultaneous change $z\rightarrow-z$,
$b\rightarrow-b$, while the imaginary part is anti-symmetric, for fixed
$P^{z}$. For $P^{z}>0$, we use the operator $O^{++}$ for $z>0$, and the
operator $O^{--}_{c}$ for $z<0$. Using these properties, we symmetrize our
results between positive and negative values of $z$. Our results show a clear
divergence in $L$, as expected Ji et al. (2021a).

Figure 3: Real and imaginary parts of the bare matrix elements for a
transverse separation of $b/a=1$.

In Fig. 4 we show the dependence of the bare beam function on $b$. As expected
Ji et al. (2021a), it decays exponentially with $b$.

Figure 4: Dependence of the bare matrix element on the transverse separation
$b$ for $L/a=10$ and $z/a=2$.

##  V Non-perturbative renormalization

As in the case of the straight Wilson line, there are logarithmic and linear
divergences associated with the length $z$ of the staple-shaped link. However,
here there are two extra divergences: one, which is associated with the cusps
of the staple, and one associated with the gluon exchange between the gauge
links of length $L\rightarrow\infty$, the so-called pinch-pole singularities.
On the other hand, the finite transverse separation $b$ mitigates possible
discretization effects associated with small $z$ separation, mainly in the
nonperturbative region (large $b$), which is the main interest of a lattice
calculation. In the following, we will analyze three different ways to carry
out the renormalization, namely the RI/MOM scheme, the short distance ratio
scheme, and a modified version of the RI/MOM scheme, where the renormalization
constants are computed at short distances.

###  V.1 RI/MOM

The RI/MOM scheme Martinelli et al. (1995) was first adapted for non-local
operators employed in the quasi-distribution approach in Refs. Constantinou
and Panagopoulos (2017); Alexandrou et al. (2017). The RI/MOM renormalization
constants $Z^{RI}_{\cal O}$ are defined by the condition

$\frac{Z^{\text{RI}}_{\cal
O}(z,b,L,\mu_{0};1/a)}{Z^{\text{RI}}_{q}(\mu_{0};1/a)}\frac{1}{12}{\rm
Tr}\left[\frac{\Lambda^{\Gamma}_{0}(z,b,L,p;1/a)\Gamma^{{}^{\prime}}}{e^{ip^{z}z+ip_{\perp}b}}\right]\Bigg{|}_{p^{2}{=}\mu_{0}^{2}}{=}1\,,$
(10)

where $\Lambda^{\Gamma}_{0}$ is defined in terms of the amputated Green’s
function

$\Lambda^{\Gamma}_{0}(z,b,L,p;1/a)=S_{q}^{-1}G^{\Gamma}(z,b,L,p;1/a)S_{q}^{-1}\,,$
(11)

with $S_{q}$ the off-shell quark propagator. The Green’s function is
calculated as

$G^{\Gamma}(z,b,p,L;1/a)=\langle q(p)|{\cal
O}^{\Gamma}(z,b,p,L;1/a)|q(p)\rangle\,.$ (12)

Because $G^{\Gamma}$, and thus $\Lambda^{\Gamma}$, have the same divergences
as $B_{0,\Gamma}$, all the divergences, in principle, cancel in the
renormalization procedure. $Z^{\text{RI}}_{q}$ is the quark wave function
renormalization defined as

$Z^{\text{RI}}_{q}(\mu_{0};1/a)=\frac{1}{12}{\rm
Tr}\left[(S_{q}(p;1/a))^{-1}\,S_{q}^{\rm
Born}(p)\right]\Bigg{|}_{p^{2}=\mu_{0}^{2}}\,,$ (13)

and the corresponding renormalized beam function is then

$B^{\text{RI}}_{\Gamma}(z,b,\mu_{0},P^{z})=\sum_{\Gamma^{{}^{\prime}}}\left[Z^{\text{RI}}_{\cal
O}(z,b,\mu_{0};1/a)\right]_{\Gamma\Gamma^{{}^{\prime}}}B_{0,\Gamma^{{}^{\prime}}}(z,b,P^{z};1/a)\,.$
(14)

The case of interest in this work is $\Gamma=\gamma_{0}$, computed using the
lattice ensemble of Table 3. As pointed out in Section III, $\gamma^{0}$ mixes
with $\gamma_{0}\gamma_{2}$, $\gamma_{0}\gamma_{3}$, and
$\gamma_{5}\gamma_{1}$. The renormalization factors for the diagonal and off-
diagonal cases are shown in Figs. 5 and 6, respectively. As can be seen, the
contributions from the off-diagonal mixing terms are $\lesssim 4\%$ of the
diagonal contribution for all $z$ when $b=1a,\,2a$.

Figure 5: Divergence of the RI/MOM renormalization factor as a function of $z$
for $b/a=1,2$ and at fixed $L/a=10$. Figure 6: Contribution of the off-
diagonal mixing terms compared to the diagonal one, in the RI/MOM
renormalization matrix at fixed $L/a=10$.

It is interesting to examine how the relative size of the off-diagonal to the
diagonal renormalization factors changes as $b$ increases, as the primary
focus of lattice computation of TMDPDFs is the the non-perturbative region in
the transverse separation $b$. We show in Fig. 7 the renormalization factors
as a function of $b$ for the diagonal and off-diagonal contributions. We see a
contribution of $\lesssim 7-8\%$ up to $b=6a$ for the off-diagonal terms.
Although their size shows a possible tendency of growth with increasing $b$,
this growth seems to be only moderate. We notice that the authors of Ref.
Shanahan et al. (2020) considered the entire set of Dirac structures for the
operator mixing, and observe similar results to ours, where only the relevant
operators allowed by symmetry were used. The contribution from
$\gamma_{5}\gamma_{1}$ is much larger than that of $\gamma_{0}\gamma_{2}$ and
$\gamma_{0}\gamma_{3}$, and there is a steady increase going to larger $b$
values. The magnitude of mixing for these operators found in Ref. Shanahan et
al. (2020) is also comparable to our findings. However, considering the
dominant diagonal renormalization factors, it seems that there are remaining
divergences related to large valus of $z$ and $b$, as also observed in Ref.
Zhang et al. (2022). Combined with the fact that the bare matrix element
decays exponentially with $b$ Ji et al. (2021a), it is significantly hard to
renormalize for large transverse separations using this scheme. However, in
the small $b$ region where we can extract information using RI/MOM, it can be
confirmed that the mixing is negligible. As an example, we show in Fig. 8 the
renormalized beam function for the $b=1a$ case, with circles denoting the
RI/MOM procedure using the full mixing matrix and with the triangles denoting
the case of only taking into account the diagonal contribution. Considering
that including the mixing has a negligible effect, one can safely ignore
operator mixing for these values of $b$.

Figure 7: The diagonal and off-diagonal contributions to the RI/MOM
renormalization factors as a function of $b$ at fixed $L/a=10$ and $z/a=2$.

Figure 8: Real and imaginary parts of the RI/MOM renormalized matrix elements
for a transverse separation of $b/a=1$.

###  V.2 The short distance ratio scheme

According to Figs. 6 and 7, mixing with different operators can be neglected
at least up to $b/a\lesssim 6$, and we can assume that the renormalization of
the staple-shaped link is multiplicative. This justifies the approach taken in
Refs. Zhang et al. (2022); He et al. (2022).

We first note that the vacuum expectation value of a rectangular Wilson loop
$Z_{E}$ with sides $2L+z$ and $b$,

$Z_{E}(b,2L+z;1/a)=\frac{1}{3}\operatorname{Tr}\langle 0|{\cal
W}(b;2L+z)W_{\perp}(x+b;b)|0\rangle,$ (15)

has, by construction, the same divergences as that of the staple-shaped gauge
link as defined in Eq. (5). Thus, it should cancel the pinch-pole singularity
associated with the length $L$ of the staple, as well as the divergences
associated with the cusps. Furthermore, as the sides of $Z_{E}$ are also
dependent on the longitudinal displacement $z$, the exponential divergence
associated with $z$ present in the staple-shaped link must also be cancelled
if an appropriated ratio is taken, namely the one proposed in Refs. Ji et al.
(2020a); Zhang et al. (2022):

$B_{\Gamma}(z,b,P^{z};1/a)=\lim_{L\rightarrow\infty}\frac{B_{0,\Gamma}(z,b,L,P^{z};1/a)}{\sqrt{Z_{E}(b,2L+z;1/a)}}$
(16)

To illustrate this point, we show in Fig. 9 the beam function as a function of
$L$, for a fixed $z=2a$ and two values of $b$, before and after dividing by
$\sqrt{Z_{E}}$. The cancellation of the divergences related to the length of
the Wilson line is explicitly shown. This ratio, thus, takes care of the
divergence associated with the length $L$ and width $b$ of the staple-shaped
operator.

Figure 9: The effect of taking the ratio of the bare matrix elements with
$\sqrt{Z_{E}}$ at a fixed $z/a=2$.

After dividing the beam function by the square root of $Z_{E}$, the only
remaining divergences are the UV divergences associated with the quark field
and its end-points connecting to the gauge links. As such, they can be
cancelled by taking ratios, as they have a multiplicative nature Ebert et al.
(2020); Green et al. (2020). Hence, one can define

$B^{SDR}_{\Gamma}(z,z_{0},b,b_{0},P^{z})=Z^{SDR}(z_{0},b_{0};1/a)\,B_{\Gamma}(z,b,P^{z};1/a)\,,$
(17)

where,

$Z^{SDR}(z_{0},b_{0};1/a)=\frac{1}{B_{\Gamma}(z=z_{0},b=b_{0},P^{z}=0;1/a)}\,.$
(18)

Because the remaining divergences are independent of the length of the Wilson
line, one is free to choose $z_{0}$ and $b_{0}$. However, because we want to
connect these quantities to the $\overline{\text{MS}}$ scheme via a
perturbative scheme conversion, $z_{0}$ and $b_{0}$ should be small enough for
perturbation theory to be valid.

The renormalized matrix elements are converted to the $\overline{\text{MS}}$
scheme using perturbation theory. We have computed the vertex, sail, and
tadpole one-loop diagrams for external quark states with a general momentum
$p^{z}$. For $p^{z}\rightarrow 0$, we obtain

$Z^{\overline{\text{MS}},SDR}(z_{0},b_{0},\mu_{0})=1+\frac{\alpha_{s}C_{F}}{2\pi}\left(\frac{1}{2}+\frac{3}{2}\ln\left(\frac{b_{0}^{2}+z_{0}^{2}}{4e^{-2\gamma_{E}}}\mu_{0}^{2}\right)-2\frac{z_{0}}{b_{0}}\textrm{arctan}\frac{z_{0}}{b_{0}}\right)\,,$
(19)

which agrees with Eq. (6) of Ref. Zhang et al. (2022). Note that this factor
equals $B_{\Gamma}^{\overline{\text{MS}}}(z=z_{0},b=b_{0},P^{z}=0,\mu_{0})$.
Details of this calculation for a general external momentum will be presented
in Ref. Constantinou et al. (2023). The renormalized beam function in the
$\overline{\text{MS}}$ scheme is then given by

$B^{\overline{\text{MS}}}_{\Gamma}(z,b,P^{z})=Z^{\overline{\text{MS}}}(z_{0},b_{0},\mu_{0})B^{SDR}_{\Gamma}(z,z_{0},b,b_{0},P^{z}).$
(20)

Results for $B^{\overline{\text{MS}}}(z,b,P^{z})$ are presented in Section VI,
where we also discuss the cancellation of the $z_{0},b_{0}$ dependence in Eq.
(20).

###  V.3 Short distance RI/MOM

As discussed in Section V.1, the usual RI scheme may be problematic at large
$z$ and $b$ as the magnitude of the $Z^{\text{RI}}$ factors grows
exponentially. Also, as shown in Refs. Zhang et al. (2021, 2022), the usual RI
scheme may still contain a residual linear divergence, which may not be
properly canceled. On the other hand, as discussed in Section V.2, the
$\sqrt{Z_{E}}$ factor cancels all divergences in $z$ and $b$ present in the
staple. Hence, we can define a vertex function that is free of such
divergences,

$\Lambda^{\Gamma}(z,b,p;1/a)=\frac{\Lambda^{\Gamma}_{0}(z,b,p;1/a)}{\sqrt{Z_{E}(b,2L+z;1/a)}}.$
(21)

Because the divergences related to the lengths of the Wilson line have been
removed, we can compute the renormalization factors as in Section V.2 at some
fixed $z_{0},b_{0}$ Ji et al. (2021b)

$\displaystyle\frac{Z^{\text{RI}-\text{short}}_{\cal
O}(z_{0},b_{0},\mu_{0};1/a)}{Z^{\text{RI}}_{q}(\mu_{0};1/a)}\frac{1}{12}{\rm
Tr}\left[\frac{\Lambda^{\Gamma}(z,b,p;1/a)\Gamma^{{}^{\prime}}}{e^{ip^{z}z+ibp_{\perp}}}\right]\Bigg{|}_{p^{2}{=}\mu_{0}^{2},z=z_{0},b=b_{0}}{=}1\,.$
(22) $\displaystyle\,.$ (23)

In principle, one can do this subtraction from the vertex in the standard
RI/MOM. This would reduce the growth of the $Z$-factors with increasing $b$.
However, in combining $Z^{\text{RI}}$ with the bare $B_{\gamma_{0}}$, the
$\sqrt{Z_{E}}$ factor cancels, and therefore this subtraction has no effect on
the renormalized matrix elements. But if we fix the RI/MOM factors at a short
perturbative scale, then this subtraction must be performed in order to cancel
the cusp and linear divergences.

As in Section V.2, we choose the pair $z_{0},b_{0}$ to be in the perturbative
region in order to make the perturbative conversion to the
$\overline{\text{MS}}$ scheme reliable. We label this procedure RI-short.

Figure 10: The $Z^{\text{RI}-\text{short}}$ renormalization factor within the
perturbative range $b_{0}/a=1$ as a function of the length $L$ of the staple.

The corresponding renormalized beam function in this scheme is then given by

$B^{\text{RI}-\text{short}}_{\Gamma}(z,z_{0},b,b_{0},\mu_{0},P^{z})=\sum_{\Gamma^{{}^{\prime}}}\left[Z^{\text{RI}-\text{short}}_{\cal
O}(z_{0},b_{0},\mu_{0};1/a)\right]_{\Gamma\Gamma^{{}^{\prime}}}B_{\Gamma^{{}^{\prime}}}(z,b,P^{z};1/a).$
(24)

For illustration, we show in Fig 10 the renormalization factors
$Z^{\text{RI}-\text{short}}_{\cal O}(z_{0},b_{0},\mu_{0};1/a)$ for
$z_{0}=0,\,1a$ and $b_{0}=1a$. We also show the off-diagonal renormalization
factors, which, in principle, mix with $\gamma^{0}$. As expected, they are
also independent of $L$ and can be omitted in the full calculation as their
contribution is negligible as compared to the diagonal contribution. Finally,
we convert $B^{\text{RI}-\text{short}}_{\Gamma}$ to the $\overline{\text{MS}}$
scheme using one-loop perturbation theory,

$B^{\overline{\text{MS}}}_{\Gamma}(z,b,\mu_{0},P^{z})=\sum_{\Gamma^{{}^{\prime}}}\left[Z^{\overline{\text{MS}},\text{RI}-\text{short}}_{\cal
O}(z_{0},b_{0},\mu_{0})\right]_{\Gamma\Gamma^{{}^{\prime}}}B^{\text{RI}-\text{short}}_{\Gamma^{{}^{\prime}}}(z,z_{0},b,b_{0},\mu_{0},P^{z})\,.$
(25)

The conversion matrix $Z_{O}^{\overline{\text{MS}},\text{RI}-\text{short}}$ is
calculated in dimensional regularization for arbitrary values of the momentum
scale $\mu_{0}$. Explicit expressions for all $\Gamma,\Gamma^{\prime}$ are
given in Ref. Constantinou et al. (2023). Similar perturbative studies can be
found in Refs. Constantinou et al. (2019); Ebert et al. (2020). A comparison
between the conversion factor in RI-short and the SDR scheme is shown in Fig.
11 for the case $\Gamma=\Gamma^{\prime}=\gamma_{0}$. We observe that the real
parts are compatible between the two schemes, while the RI-short scheme shows
a non-zero imaginary part, in contrast to the ratio scheme. This is a
consequence of using a non-zero momentum scale.

Figure 11: Conversion factor to the $\overline{\text{MS}}$ scheme (top: real
part, bottom: imaginary part) for the RI-short and SDR schemes at $b/a=1$
(left) and $b/a=2$ (right).

##  VI Results

In this Section, we present the renormalized beam functions in the 3
renormalization schemes discussed in Section V, as well as the corresponding
results in the $\overline{\text{MS}}$ scheme.

In Fig. 12, we show the real and imaginary parts of
$B_{\gamma_{0}}(z,b,\mu_{0},P^{z})$ as a function of $z$ for two values of the
transverse separation, $b/a=1,2$. For the real parts, the schemes agree within
errors for $|z/a|\lessapprox 6$. For $|z/a|\gtrapprox 6$,
$B^{SDR}_{\gamma_{0}}$ and $B^{\text{RI}-\text{short}}_{\gamma_{0}}$ are
consistent, while $B^{\text{RI}}_{\gamma_{0}}$ becomes increasingly more
negative with increasingly larger errors which may indicate the presence of a
residual divergence in the standard RI/MOM. For the imaginary parts, we
observe a similar behavior. Namely, the SDR and the RI-short schemes are
almost identical for any $|z/a|$ for the two values of $b/a$ considered. For
the usual RI/MOM, we observe deviations not only for $|z/a|\gtrapprox 6$ but
also for most of the $|z/a|$ region with an increase in errors with increasing
$|z/a|$.

Figure 12: Renormalized matrix elements for $z_{0}/a=b_{0}/a=1$ for the SDR
and the $\text{RI}-\text{short}$ schemes, and $P^{z}\simeq 1.7$ GeV.

The conversion to the $\overline{\text{MS}}$ scheme is done using one-loop
perturbation theory for all three schemes considered. For the two cases where
the renormalization factors are fixed at a short perturbative scale, we use
$z_{0}/a=b_{0}/a=1$ when computing
$B^{\text{RI}-\text{short}}_{\Gamma^{{}^{\prime}}}(z,z_{0},b,b_{0},\mu_{0},P^{z})$
and $B^{SDR}_{\Gamma}(z,z_{0},b,b_{0},P^{z})$. Since the final results in the
$\overline{\text{MS}}$ scheme should be independent of the values we choose
for the $z_{0},b_{0}$ pair, we examine the stability of our results on the
choice of the values for this pair. To this end, we show in Fig. 13 the
product of the renormalization factors for each scheme with the corresponding
conversion factors to the $\overline{\text{MS}}$ scheme,
$Z^{\overline{\text{MS}},SDR}_{\gamma_{0}}Z^{SDR}_{\gamma_{0}}$ and
$Z^{\overline{\text{MS}},\text{RI}}_{\gamma_{0}}Z^{\text{RI}-\text{short}}_{\gamma_{0}}$.
These products should be independent on the choice of values of $z_{0},b_{0}$
as long as the one-loop perturbative calculation is a reliable approximation.
We notice that these products are equal to the renormalization factors in the
$\overline{\text{MS}}$ scheme in the absence of mixing. For the ratio scheme,
there is a $10\%-20\%$ correction when going from $b_{0}/a=1$ to $b_{0}/a=2$
for the three values of $z_{0}$ considered here, with the correction
increasing for larger values of $z_{0}$. This indicates a limitation on the
use of one-loop perturbation theory to perform the conversion to the
$\overline{\text{MS}}$ scheme already at a transverse separation as low as
$b/a=2$. As a consequence, $\mathcal{O}(\alpha_{s}^{2})$ corrections cannot be
disregarded. For the $z_{0}$ dependence, the choice of $z_{0}/a=0$ and
$z_{0}/a=1$ is nearly equivalent within one sigma error. For $z/a=2$, the
discrepancy grows, especially for $b/a=2$. For the RI-short scheme, the
situation is not significantly changed. From this discussion, we conclude that
our choice for the pair $z_{0},b_{0}$ is a likely safe region to fix the short
perturbative scale.

Figure 13: Dependence of the product between the $\overline{\text{MS}}$
conversion factors and the SDR and $\text{RI}-\text{short}$ renormalization
factors on the choice of the short perturbative scale.

We show in Fig. 14 the renormalized $\overline{\text{MS}}$ matrix elements
computed in the ratio and the RI-short schemes at $P^{z}\simeq 1.7$ GeV. The
imaginary part of the matrix elements are in full agreement, the same
happening in the the real part, for large $z$ values. In the small $z$ region
of the real part there is a tendency to discrepancy, as already observed in
the intermediate schemes shown in Fig. 12, although they agree within one
sigma error.

Figure 14: $\overline{\text{MS}}$ renormalized matrix elements converted from
the SDR and the $\text{RI}-\text{short}$ schemes, using $z_{0}/a=b_{0}/a=1$
and $P^{z}\simeq 1.7$ GeV

##  VII Conclusion

In this work, we have addressed the non-perturbative renormalization of the
asymmetric staple-shaped quark bilinear operators. We used symmetry arguments
to restrict the possible mixing between the operators, and we found that
mixing is allowed within sets of four operators. We then employed the RI/MOM
scheme to estimate the importance of mixing coming from the non-diagonal
terms, and we found that they can be neglected up to transverse separations of
at least $6a$. This result justifies the use of multiplicative renormalization
for the asymmetric staple-shaped operator. Based on this conclusion, we
computed the renormalization factors in the SDR scheme, and in the RI-short
scheme, where the renormalization factors are computed at short distances. We
found that both schemes are nearly equivalent over all the $|z|$ region
considered. Finally, we converted both schemes to the $\overline{\text{MS}}$
scheme and presented the corresponding results for the renormalized beam
function. Given our previous computation of the Collins-Soper kernel and of
the soft function Li et al. (2022), the next natural step will be to compute
the TMDPDFs themselves over a wide range of $b$, and we plan to present
results for them in the near future.

###### Acknowledgements.

A.S and F.S are funded by the NSFC and the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) through the funds provided to the Sino-
German Collaborative Research Center TRR110 “Symmetries and the Emergence of
Structure in QCD” (NSFC Grant No. 12070131001, DFG Project-ID 196253076 - TRR
110). K.C. is supported by the National Science Centre (Poland) grant SONATA
BIS No.2016/22/E/ST2/00013 and OPUS grant No. 2021/43/13/ST2/00497. M.C.
acknowledges financial support by the U.S. Department of Energy, Office of
Nuclear Physics, Early Career Award under Grant No. DE-SC0020405. C.A. and
J.T. acknowledge partial funding from the project ”Unraveling the 3D Parton
structure of the nucleon with lattice QCD” (contract id EXCELLENCE/0421/0043)
co-financed by the European Regional Development Fund and the Republic of
Cyprus through the Research and Innovation Foundation. G.S. acknowledges
funding from the European High-Performance Computing Joint Undertaking (JU)
under grant agreement No. 951732 and from the European Regional Development
Fund and the Republic of Cyprus through the Research and Innovation Foundation
under contract No. EXCELLENCE/0421/0025. X.F. is supported in part by NSFC of
China under Grant No. 11775002 and National Key Research and Development
Program of China under Contracts No.2020YFA0406400. X.F and C.L. are supported
in part by NSFC of China under Grant No. 12070131001. C.L. is also supported
in part by NSFC of China under Grant No. 11935017, No. 12293060, No. 12293062
and No. 12293063. This work used computational resources from the John von
Neumann-Institute for Computing on the Juwels booster system at the research
center in Juelich, under the project with id TMDPDF1.

## References

  * Amoroso et al. (2022) S. Amoroso et al., Acta Phys. Polon. B 53, A1 (2022), eprint 2203.13923.
  * Bacchetta et al. (2017) A. Bacchetta, F. Delcarro, C. Pisano, M. Radici, and A. Signori, JHEP 06, 081 (2017), [Erratum: JHEP 06, 051 (2019)], eprint 1703.10157.
  * Bertone et al. (2019) V. Bertone, I. Scimemi, and A. Vladimirov, JHEP 06, 028 (2019), eprint 1902.08474.
  * Scimemi and Vladimirov (2020) I. Scimemi and A. Vladimirov, JHEP 06, 137 (2020), eprint 1912.06532.
  * Bury et al. (2022) M. Bury, F. Hautmann, S. Leal-Gomez, I. Scimemi, A. Vladimirov, and P. Zurita, JHEP 10, 118 (2022), eprint 2201.07114.
  * Bacchetta et al. (2022) A. Bacchetta, V. Bertone, C. Bissolotti, G. Bozzi, M. Cerutti, F. Piacenza, M. Radici, and A. Signori (MAP), JHEP 10, 127 (2022), eprint 2206.07598.
  * Barry et al. (2023) P. C. Barry, L. Gamberg, W. Melnitchouk, E. Moffat, D. Pitonyak, A. Prokudin, and N. Sato (2023), eprint 2302.01192.
  * Dudek et al. (2012) J. Dudek et al., Eur. Phys. J. A 48, 187 (2012), eprint 1208.1244.
  * Abdul Khalek et al. (2021) R. Abdul Khalek et al. (2021), eprint 2103.05419.
  * Ji (2013) X. Ji, Phys. Rev. Lett. 110, 262002 (2013), eprint 1305.1539.
  * Radyushkin (2017) A. V. Radyushkin, Phys. Rev. D 96, 034025 (2017), eprint 1705.01488.
  * Cichy and Constantinou (2019) K. Cichy and M. Constantinou, Adv. High Energy Phys. 2019, 3036904 (2019), eprint 1811.07248.
  * Radyushkin (2020) A. Radyushkin, Int. J. Mod. Phys. A 35, 2030002 (2020), eprint 1912.04244.
  * Constantinou (2021) M. Constantinou, Eur. Phys. J. A 57, 77 (2021), eprint 2010.02445.
  * Ji et al. (2021a) X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang, and Y. Zhao, Rev. Mod. Phys. 93, 035005 (2021a), eprint 2004.03543.
  * Cichy (2022a) K. Cichy, PoS LATTICE2021, 017 (2022a), eprint 2110.07440.
  * Cichy (2022b) K. Cichy, EPJ Web Conf. 258, 01005 (2022b), eprint 2111.04552.
  * Constantinou et al. (2022) M. Constantinou et al. (2022), eprint 2202.07193.
  * Ji et al. (2015) X. Ji, P. Sun, X. Xiong, and F. Yuan, Phys. Rev. D 91, 074009 (2015), eprint 1405.7640.
  * Ji et al. (2019) X. Ji, L.-C. Jin, F. Yuan, J.-H. Zhang, and Y. Zhao, Phys. Rev. D 99, 114006 (2019), eprint 1801.05930.
  * Ji et al. (2020a) X. Ji, Y. Liu, and Y.-S. Liu, Nucl. Phys. B 955, 115054 (2020a), eprint 1910.11415.
  * Ebert et al. (2022) M. A. Ebert, S. T. Schindler, I. W. Stewart, and Y. Zhao, JHEP 04, 178 (2022), eprint 2201.08401.
  * Zhang et al. (2020) Q.-A. Zhang et al. (Lattice Parton), Phys. Rev. Lett. 125, 192001 (2020), eprint 2005.14572.
  * Li et al. (2022) Y. Li et al., Phys. Rev. Lett. 128, 062002 (2022), eprint 2106.13027.
  * Ebert et al. (2019a) M. A. Ebert, I. W. Stewart, and Y. Zhao, Phys. Rev. D 99, 034505 (2019a), eprint 1811.00026.
  * Ebert et al. (2020) M. A. Ebert, I. W. Stewart, and Y. Zhao, JHEP 03, 099 (2020), eprint 1910.08569.
  * Shanahan et al. (2020) P. Shanahan, M. L. Wagman, and Y. Zhao, Phys. Rev. D 101, 074505 (2020), eprint 1911.00800.
  * Hagler et al. (2009) P. Hagler, B. U. Musch, J. W. Negele, and A. Schafer, EPL 88, 61001 (2009), eprint 0908.1283.
  * Musch et al. (2012) B. U. Musch, P. Hagler, M. Engelhardt, J. W. Negele, and A. Schafer, Phys. Rev. D 85, 094510 (2012), eprint 1111.4249.
  * 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 _33rd International Symposium on Lattice Field Theory_ (SISSA, 2015), eprint 1601.05717.
  * 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. D 96, 094508 (2017), eprint 1706.03406.
  * He et al. (2022) J.-C. He, M.-H. Chu, J. Hua, X. Ji, A. Schäfer, Y. Su, W. Wang, Y. Yang, J.-H. Zhang, and Q.-A. Zhang (LPC) (2022), eprint 2211.02340.
  * Constantinou et al. (2019) M. Constantinou, H. Panagopoulos, and G. Spanoudes, Phys. Rev. D 99, 074508 (2019), eprint 1901.03862.
  * Ji et al. (2021b) Y. Ji, J.-H. Zhang, S. Zhao, and R. Zhu, Phys. Rev. D 104, 094510 (2021b), eprint 2104.13345.
  * Constantinou and Panagopoulos (2017) M. Constantinou and H. Panagopoulos, Phys. Rev. D 96, 054506 (2017), eprint 1705.11193.
  * Zhang et al. (2022) K. Zhang, X. Ji, Y.-B. Yang, F. Yao, and J.-H. Zhang ([Lattice Parton Collaboration (LPC)]), Phys. Rev. Lett. 129, 082002 (2022), eprint 2205.13402.
  * Orginos et al. (2017) K. Orginos, A. Radyushkin, J. Karpie, and S. Zafeiropoulos, Phys. Rev. D 96, 094503 (2017), eprint 1706.05373.
  * Braun et al. (2019) V. M. Braun, A. Vladimirov, and J.-H. Zhang, Phys. Rev. D 99, 014013 (2019), eprint 1810.00048.
  * Ji et al. (2021c) X. Ji, Y. Liu, A. Schäfer, W. Wang, Y.-B. Yang, J.-H. Zhang, and Y. Zhao, Nucl. Phys. B 964, 115311 (2021c), eprint 2008.03886.
  * Ji et al. (2020b) X. Ji, Y. Liu, and Y.-S. Liu, Phys. Lett. B 811, 135946 (2020b), eprint 1911.03840.
  * Collins (2011) J. Collins, Int. J. Mod. Phys. Conf. Ser. 4, 85 (2011), eprint 1107.4123.
  * Echevarria et al. (2016) M. G. Echevarria, I. Scimemi, and A. Vladimirov, Phys. Rev. D 93, 054004 (2016), eprint 1511.05590.
  * Ebert et al. (2019b) M. A. Ebert, I. W. Stewart, and Y. Zhao, JHEP 09, 037 (2019b), eprint 1901.03685.
  * Shindler (2008) A. Shindler, Phys. Rept. 461, 37 (2008), eprint 0707.4093.
  * Chen et al. (2019) J.-W. Chen, T. Ishikawa, L. Jin, H.-W. Lin, J.-H. Zhang, and Y. Zhao (LP3), Chin. Phys. C 43, 103101 (2019), eprint 1710.01089.
  * Alexandrou et al. (2018) C. Alexandrou et al., Phys. Rev. D 98, 054518 (2018), eprint 1807.00495.
  * Bali et al. (2016) G. S. Bali, B. Lang, B. U. Musch, and A. Schäfer, Phys. Rev. D 93, 094515 (2016), eprint 1602.05525.
  * Martinelli et al. (1995) G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa, and A. Vladikas, Nucl. Phys. B 445, 81 (1995), eprint hep-lat/9411010.
  * Alexandrou et al. (2017) C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, H. Panagopoulos, and F. Steffens, Nucl. Phys. B 923, 394 (2017), eprint 1706.00265.
  * Green et al. (2020) J. R. Green, K. Jansen, and F. Steffens, Phys. Rev. D 101, 074509 (2020), eprint 2002.09408.
  * Constantinou et al. (2023) M. Constantinou, H. Panagopoulos, and G. Spanoudes, In preparation (2023).
  * Zhang et al. (2021) K. Zhang, Y.-Y. Li, Y.-K. Huo, A. Schäfer, P. Sun, and Y.-B. Yang ($\chi$QCD), Phys. Rev. D 104, 074501 (2021), eprint 2012.05448.

##  Appendix A Symmetry properties

As explained in Section III, we look at the mixing properties of the staple-
shaped operator using symmetry arguments. We observe the transformation of the
object $(ijkl)_{c}$, where $i,j,k,l$ denote the signs of
$O^{++},O^{--},O^{+-}$ and $O^{-+}$ respectively, and $c=\pm$ denoted the sign
of the charge conjugated operators. This gives rise to 4 relevant
combinations: $(++++)_{c},(+-+-)_{c},(++--)_{c}$, and $(+--+)_{c}$. In Table 4
we summarize the symmetries of the 16 Dirac structures, where the 4 signs
under any transformation denote the symmetry properties of the 4 relevant
combinations, respectively.

| $\gamma_{0}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$ | $\mathrm{1}$ | $\gamma_{5}$ | $\gamma_{5}\gamma_{0}$ | $\gamma_{5}\gamma_{1}$
---|---|---|---|---|---|---|---|---
${\cal P}_{F0}^{1,2}$ | $-+-+$ | $+-+-$ | $+-+-$ | $+-+-$ | $-+-+$ | $+-+-$ | $+-+-$ | $-+-+$
${\cal P}_{F1}^{1,2}$ | $+-+-$ | $-+-+$ | $+-+-$ | $+-+-$ | $-+-+$ | $+-+-$ | $-+-+$ | $+-+-$
${\cal P}_{F2}^{1,2}$ | $++--$ | $++--$ | $--++$ | $++--$ | $--++$ | $++--$ | $--++$ | $--++$
${\cal P}_{F3}^{1,2}$ | $+--+$ | $+--+$ | $+--+$ | $-++-$ | $-++-$ | $+--+$ | $-++-$ | $-++-$
${\cal T}_{F0}^{1,2}$ | $++++$ | $----$ | $----$ | $----$ | $----$ | $++++$ | $----$ | $++++$
${\cal T}_{F1}^{1,2}$ | $----$ | $++++$ | $----$ | $----$ | $----$ | $++++$ | $++++$ | $----$
${\cal T}_{F2}^{1,2}$ | $-++-$ | $-++-$ | $+--+$ | $-++-$ | $-++-$ | $+--+$ | $+--+$ | $+--+$
${\cal T}_{F3}^{1,2}$ | $--++$ | $--++$ | $--++$ | $++--$ | $--++$ | $++--$ | $++--$ | $++--$
$C$ | $c$ | $c$ | $c$ | $c$ | $-c$ | $-c$ | $-c$ | $-c$

| $\gamma_{5}\gamma_{2}$ | $\gamma_{5}\gamma_{3}$ | $\gamma_{0}\gamma_{1}$ | $\gamma_{0}\gamma_{2}$ | $\gamma_{0}\gamma_{3}$ | $\gamma_{1}\gamma_{2}$ | $\gamma_{1}\gamma_{3}$ | $\gamma_{2}\gamma_{3}$
---|---|---|---|---|---|---|---|---
${\cal P}_{F0}^{1,2}$ | $-+-+$ | $-+-+$ | $+-+-$ | $+-+-$ | $+-+-$ | $-+-+$ | $-+-+$ | $-+-+$
${\cal P}_{F1}^{1,2}$ | $-+-+$ | $-+-+$ | $+-+-$ | $-+-+$ | $-+-+$ | $+-+-$ | $+-+-$ | $-+-+$
${\cal P}_{F2}^{1,2}$ | $++--$ | $--++$ | $--++$ | $++--$ | $--++$ | $++--$ | $--++$ | $++--$
${\cal P}_{F3}^{1,2}$ | $-++-$ | $+--+$ | $-++-$ | $-++-$ | $+--+$ | $-++-$ | $+--+$ | $+--+$
${\cal T}_{F0}^{1,2}$ | $++++$ | $++++$ | $++++$ | $++++$ | $++++$ | $----$ | $----$ | $----$
${\cal T}_{F1}^{1,2}$ | $++++$ | $++++$ | $++++$ | $----$ | $----$ | $++++$ | $++++$ | $----$
${\cal T}_{F2}^{1,2}$ | $-++-$ | $+--+$ | $-++-$ | $+--+$ | $-++-$ | $+--+$ | $-++-$ | $+--+$
${\cal T}_{F3}^{1,2}$ | $++--$ | $--++$ | $--++$ | $--++$ | $++--$ | $--++$ | $++--$ | $++--$
$C$ | $-c$ | $-c$ | $c$ | $c$ | $c$ | $c$ | $c$ | $c$

Table 4: General symmetry properties of all possible Dirac structures. The 4
signs correspond to symmetry properties of combinations $(++++)_{c}$ (or its
equivalent $(----)_{c}$), $(+-+-)_{c}$ (or $(-+-+)_{c}$), $(++--)_{c}$ (or
$(--++)_{c}$) and $(+--+)_{c}$ (or $(-++-)_{c}$). The last row indicates the
symmetry properties with respect to charge conjugation ($c$ \- symmetric, $-c$
\- antisymmetric).

##  Appendix B Mixing sets

In this appendix we provide all the mixing sets of asymmetric staple-shaped
operators with different Dirac structures, as dictated by the generalized C,
P, T symmetries of Eqs. (6 – 8), in explicit form:

  * •

$\\{1,\gamma_{2},\gamma_{3},\gamma_{2}\gamma_{3}\\}$,

  * •

$\\{\gamma_{5},\gamma_{5}\gamma_{2},\gamma_{5}\gamma_{3},\gamma_{0}\gamma_{1}\\}$,

  * •

$\\{\gamma_{0},\gamma_{0}\gamma_{2},\gamma_{0}\gamma_{3},\gamma_{5}\gamma_{1}\\}$,

  * •

$\\{\gamma_{1},\gamma_{1}\gamma_{2},\gamma_{1}\gamma_{3},\gamma_{5}\gamma_{0}\\}$.

In the case of symmetric staple-shaped operators (z=0), the mixing sets are
reduced to:

  * •

$\\{\gamma_{2},\gamma_{3},\gamma_{2}\gamma_{3}\\}$,

  * •

$\\{\gamma_{5},\gamma_{5}\gamma_{2},\gamma_{5}\gamma_{3}\\}$,

  * •

$\\{\gamma_{0},\gamma_{0}\gamma_{2},\gamma_{0}\gamma_{3}\\}$,

  * •

$\\{\gamma_{1},\gamma_{1}\gamma_{2},\gamma_{1}\gamma_{3}\\}$,

while the remaining operators
$1,\gamma_{0}\gamma_{1},\gamma_{5}\gamma_{1},\gamma_{5}\gamma_{0}$ are
multiplicative renormalizable.