# Double parton distributions in the pion from lattice QCD

Gunnar S. Bali a Luca Castagnini c,d Markus Diehl e Jonathan R. Gaunt f
Benjamin Gläßle a Andreas Schäfer a and Christian Zimmermann Institute for
Theoretical Physics, University of Regensburg, 93040 Regensburg,
GermanyDepartment of Theoretical Physics, Tata Institute of Fundamental
Research, Homi Bhabha Road, Mumbai 400005, IndiaFachbereich Physik, University
of Hamburg, 22761 Hamburg, GermanyDeutsches Elektronen-Synchroton DESY, 22607
Hamburg, GermanyCERN Theory Division, 1211 Geneva 23, SwitzerlandZentrum für
Datenverarbeitung, Universität Tübingen, Wächterstr. 76, 72074 Tübingen,
Germany

###### Abstract

We perform a lattice study of double parton distributions in the pion, using
the relationship between their Mellin moments and pion matrix elements of two
local currents. A good statistical signal is obtained for almost all relevant
Wick contractions. We investigate correlations in the spatial distribution of
two partons in the pion, as well as correlations involving the parton
polarisation. The patterns we observe depend significantly on the quark mass.
We investigate the assumption that double parton distributions approximately
factorise into a convolution of single parton distributions.

††preprint:  DESY 20-098, CERN-TH-2020-086

##  1 Introduction

Matrix elements of currents in a hadron offer a variety of ways to quantify
and study hadron structure. In particular, information about _correlations_
inside the hadron can be obtained from the matrix elements of two currents
that are separated by a space-like distance. Such matrix elements can be
calculated in lattice QCD, and there has been considerable activity in this
area over the years Barad:1984px ; Barad:1985qd ; Wilcox:1986ge ;
Wilcox:1986dk ; Wilcox:1990zc ; Chu:1990ps ; Lissia:1991gv ; Burkardt:1994pw ;
Alexandrou:2002nn ; Alexandrou:2003qt ; Alexandrou:2008ru . These studies
address a broad range of physics questions, such as confinement Barad:1984px ;
Barad:1985qd , the size of hadrons Wilcox:1986ge ; Wilcox:1986dk ;
Wilcox:1990zc ; Burkardt:1994pw , density correlations Chu:1990ps , comparison
with quark models Lissia:1991gv , or the non-spherical shape of hadrons with
spin $1$ or larger Alexandrou:2002nn ; Alexandrou:2003qt ; Alexandrou:2008ru .

We continued this line of investigation in a recent paper Bali:2018nde . We
performed a lattice computation of the matrix elements of two scalar,
pseudoscalar, vector, or axial vector currents in the pion and compared our
results with predictions of chiral perturbation theory. For the first time, we
computed all Wick contractions that contribute to these matrix elements,
whilst earlier work had focused on the case in which the two currents are
inserted on different quark lines between the hadron source and sink operators
(see graph $C_{1}$ in figure 4). We obtained signals with a good statistical
accuracy for almost all contractions and were thus able to study their
relative importance. Our results were compared with different models in
Rinaldi:2020ybv ; Courtoy:2020tkd .

Extending our work in Bali:2018nde , we will in the present paper use two-
current matrix elements from the lattice to obtain information about double
parton distributions (DPDs). DPDs describe the correlated distribution of two
partons inside a hadron and appear in the cross sections for double parton
scattering, which occurs when there are two separate hard-scattering processes
in a single hadron-hadron collision. The study of this mechanism has a long
history in collider physics, from early theoretical papers such as
Landshoff:1978fq ; Kirschner:1979im ; Politzer:1980me ; Paver:1982yp ;
Shelest:1982dg ; Mekhfi:1983az ; Sjostrand:1986ep  to the detailed
investigation of QCD dynamics and factorisation that started about ten years
ago Blok:2010ge ; Diehl:2011tt ; Gaunt:2011xd ; Ryskin:2011kk ; Blok:2011bu ;
Diehl:2011yj ; Manohar:2012jr ; Manohar:2012pe ; Ryskin:2012qx ; Gaunt:2012dd
; Blok:2013bpa ; Diehl:2017kgu . After early experimental studies
Akesson:1986iv ; Alitti:1991rd , a multitude of double parton scattering
processes has been measured at the Tevatron and the LHC, see Abe:1997xk ;
Abazov:2015nnn ; Aaij:2016bqq ; Aaboud:2018tiq ; Sirunyan:2019zox  and
references therein. Some final states produced by double parton scattering are
of particular interest because they are a background to search channels for
new physics. A prominent example are like-sign gauge boson pairs $W^{+}W^{+}$
and $W^{-}W^{-}$ Kulesza:1999zh ; Gaunt:2010pi ; Sirunyan:2019zox ;
Ceccopieri:2017oqe ; Cotogno:2018mfv ; Cotogno:2020iio , the decay of which
can yield like-sign lepton pairs. A wealth of further information about double
parton scattering can be found in the monograph Bartalini:2017jkk .

Double parton distributions remain poorly known, and their extraction from
experimental data is considerably more difficult than the extraction of single
parton distributions (PDFs). It is therefore important to have as much
theoretical guidance as possible about the properties and behaviour of DPDs.
Apart from approaches that focus on fulfilling theoretical constraints
Gaunt:2009re ; Golec-Biernat:2014bva ; Golec-Biernat:2015aza ; Diehl:2020xyg ,
there exists a large number of model calculations for the DPDs of the nucleon
Chang:2012nw ; Rinaldi:2013vpa ; Broniowski:2013xba ; Rinaldi:2014ddl ;
Broniowski:2016trx ; Kasemets:2016nio ; Rinaldi:2016jvu ; Rinaldi:2016mlk  and
a smaller number for those of the pion Rinaldi:2018zng ; Courtoy:2019cxq ;
Broniowski:2019rmu ; Broniowski:2020jwk .

A relation between the Mellin moments of DPDs and two-current matrix elements
that can be computed on the lattice was written down in Diehl:2011tt ;
Diehl:2011yj . This generalises the relation between matrix elements of one
current and the Mellin moments of PDFs, which has been extensively exploited
in lattice studies, as reviewed for instance in Hagler:2009ni ; Lin:2017snn ;
Lin:2020rut . Whilst knowledge of a few Mellin moments is insufficient for
reconstructing the full DPDs, it allows one to investigate crucial features of
these functions, such as their dependence on the distance between the two
partons and on the parton polarisation. In the present paper, we pursue this
idea for the DPDs of the pion, focusing on their lowest Mellin moments. We use
the same lattice data as in our study Bali:2018nde . Corresponding work on the
DPDs of the nucleon is in progress, and preliminary results have been
presented in Zimmermann:2019quf .

This paper is organised as follows. In section 2, we recapitulate some basics
about DPDs and then elaborate on the relation between their Mellin moments and
the two-current matrix elements we compute on the lattice. This will in
particular lead us to introduce the concept of skewed DPDs. In section 3, we
describe the main elements of our lattice simulations (a full account is given
in Bali:2018nde ) and investigate several lattice artefacts that are present
in our data. Our results for zero pion momentum are presented and discussed in
section 4. In section 5, we develop a parametrisation of the data for both
zero and nonzero pion momenta, which will allow us to reconstruct the Mellin
moments of pion DPDs, albeit in a model-dependent fashion. Our main findings
are summarised in section 6.

##  2 Theory

###  2.1 Double parton distributions

To begin with, we recall some basics about double parton distributions. An
extended introduction to the subject can be found in Diehl:2017wew .

Factorisation for a double parton scattering process means that its cross
section is given in terms of hard-scattering cross sections at parton level
and double parton distributions for each of the colliding hadrons. For pair
production of colourless particles, such as $Z$, $W$ or Higgs bosons, this
factorisation can be proven rigorously. A DPD gives the joint probability for
finding in a hadron two partons with longitudinal momentum fractions $x_{1}$
and $x_{2}$ at a transverse distance $\boldsymbol{y}$ from each other. The
distributions for quarks and antiquarks are defined by operator matrix
elements as

$\displaystyle F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})=2p^{+}\int
dy^{-}\int\frac{dz^{-}_{1}}{2\pi}\,\frac{dz^{-}_{2}}{2\pi}\,$ $\displaystyle
e^{i\mskip 1.5mu(x_{1}z_{1}^{-}+x_{2}z_{2}^{-})\mskip 1.5mup^{+}}$
$\displaystyle\times\bra{h(p)}\mathcal{O}_{a_{1}}(y,z_{1})\,\mathcal{O}_{a_{2}}(0,z_{2})\ket{h(p)}\,.$
(1)

We use light-cone coordinates $v^{\pm}=(v^{0}\pm v^{3})/\sqrt{2}$ and boldface
letters for the transverse part $\boldsymbol{v}=(v^{1},v^{2})$ for any four-
vector $v^{\mu}$. The definition (2.1) refers to a reference frame in which
the transverse hadron momentum is zero, $\boldsymbol{p}=\boldsymbol{0}$. In a
frame where the hadron moves fast in the positive $z$ direction, $x_{1}$ and
$x_{2}$ can be interpreted as longitudinal momentum fractions. The vectors
$z_{1}^{\mu}$ and $z_{2}^{\mu}$ are lightlike with only $z_{1}^{-}$ and
$z_{2}^{-}$ nonzero, whereas $y^{\mu}$ is spacelike with $y^{+}=0$. The hadron
state is denoted by ${h(p)}$, and it is understood that an average over its
polarisation is taken on the r.h.s. of (2.1) if the hadron has nonzero spin.
Unless specified otherwise, the expressions of the present section hold both
for a pion and for the nucleon (and in fact for any unpolarised hadron or
nucleus).

The matrix element in (2.1) involves the same twist-two operators that appear
in the definition of ordinary PDFs. For quarks, one has

$\displaystyle\mathcal{O}_{a}(y,z)$
$\displaystyle=\bar{q}\bigl{(}y-{\textstyle\frac{1}{2}}z\bigr{)}\,\Gamma_{a}\,q\bigl{(}y+{\textstyle\frac{1}{2}}z\bigr{)}\Big{|}_{z^{+}=y^{+}=0,\,\boldsymbol{z}=\boldsymbol{0}}$
(2)

with spin projections

$\displaystyle\Gamma_{q}$ $\displaystyle={\textstyle\frac{1}{2}}\gamma^{+}\,,$
$\displaystyle\Gamma_{\Delta q}$
$\displaystyle={\textstyle\frac{1}{2}}\gamma^{+}\gamma_{5}\,,$
$\displaystyle\Gamma_{\delta
q}^{j}={\textstyle\frac{1}{2}}i\sigma^{j+}\gamma_{5}\quad(j=1,2)\,.$ (3)

The analogous expressions for antiquarks can e.g. be found in (Diehl:2011yj, ,
section 2.2). The form (2) holds in light-cone gauge $A^{+}=0$, whereas in
other gauges a Wilson line is to be inserted between the fields. Since the two
fields in (2) have light-like separation from each other, their product
requires renormalisation. This results in a scale dependence of the two
operators and of the DPD in (2.1), which we do not indicate for the sake of
brevity.

Lorentz invariance implies that one can write

$\displaystyle F_{q_{1}q_{2}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle=f_{q_{1}q_{2}}(x_{1},x_{2},y^{2})\,,$ $\displaystyle F_{\Delta
q_{1}\Delta q_{2}}(x_{1},x_{2},\boldsymbol{y})$ $\displaystyle=f_{\Delta
q_{1}\Delta q_{2}}(x_{1},x_{2},y^{2})\,,$ $\displaystyle F_{\delta
q_{1}q_{2}}^{j_{1}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle=\epsilon^{j_{1}k}\boldsymbol{y}^{k}\,mf_{\delta
q_{1}q_{2}}(x_{1},x_{2},y^{2})\,,$ $\displaystyle F_{q_{1}\delta
q_{2}}^{j_{2}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle=\epsilon^{j_{2}k}\boldsymbol{y}^{k}\,mf_{q_{1}\delta
q_{2}}(x_{1},x_{2},y^{2})\,,$ $\displaystyle F_{\delta q_{1}\delta
q_{2}}^{j_{1}j_{2}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle=\delta^{j_{1}j_{2}}f_{\delta q_{1}\delta
q_{2}}(x_{1},x_{2},y^{2})+\bigl{(}2\boldsymbol{y}^{j_{1}}\boldsymbol{y}^{j_{2}}-\delta^{j_{1}j_{2}}\boldsymbol{y}^{2}\bigr{)}\mskip
1.5mum^{2}f^{\mskip 1.5mut}_{\delta q_{1}\delta q_{2}}(x_{1},x_{2},y^{2})$ (4)

with $y^{2}=y^{\mu}y_{\mu}=-\boldsymbol{y}^{2}$. Due to parity invariance, one
has $F_{q_{1}\Delta q_{2}}=F_{\Delta q_{1}q_{2}}=0$, and time reversal
invariance implies $F_{\delta q_{1}\Delta q_{2}}=F_{\Delta q_{1}\delta
q_{2}}=0$. The hadron mass $m$ has been introduced on the r.h.s. of (2.1) so
that all scalar functions $f$ have the dimension of an inverse area. The
operator $\mathcal{O}^{j}_{\delta q}$ is a vector whose direction gives the
transverse quark spin direction, and $\epsilon^{jk}$ is the two-dimensional
antisymmetric tensor with $\epsilon^{12}=+1$. The density interpretation of
the different distributions in (2.1) is then as follows:

  * •

The unpolarised distribution $f_{q_{1}q_{2}}$ gives the probability density to
find two quarks with momentum fractions $x_{1}$ and $x_{2}$ at a transverse
distance $\boldsymbol{y}$, regardless of their polarisation.

  * •

$f_{\Delta q_{1}\Delta q_{2}}$ is the density for finding two quarks with
their longitudinal polarisations aligned minus the density for finding them
with their longitudinal polarisations anti-aligned.

  * •

$f_{\delta q_{1}\delta q_{2}}$ is the analogue of $f_{\Delta q_{1}\Delta
q_{2}}$ for transverse quark polarisations.

  * •

$f_{\delta q_{1}\mskip 1.5muq_{2}}$ describes a correlation between the
transverse polarisation of the quark $q_{1}$ and the distance $\boldsymbol{y}$
of that quark from the unpolarised quark $q_{2}$. In $f_{q_{1}\delta q_{2}}$,
the first quark is unpolarised and the second quark has transverse
polarisation.

  * •

$f^{t}_{\delta q_{1}\delta q_{2}}$ describes a correlation between the
transverse polarisations of the two quarks and their transverse distance
$\boldsymbol{y}$.

Decompositions of the same form as (2.1) can be given for the cases where one
replaces one or both of the quarks by an antiquark, with the same physical
interpretation as given above for two quarks.

Note that the polarisation dependence of DPDs is not only interesting from the
point of view of hadron structure, but can have measurable implications on
double parton scattering, as was for instance shown in Diehl:2011yj ;
Kasemets:2012pr ; Cotogno:2018mfv ; Cotogno:2020iio . Lattice calculations can
give information about the strength of the different spin correlations we just
discussed.

We note that cross sections for double parton scattering involve the product
of two DPDs integrated over the interparton distance,

$\int
d^{2}\boldsymbol{y}\;F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})\,F_{b_{1}b_{2}}(x_{1}^{\prime},x_{2}^{\prime},\boldsymbol{y})\,.$
(5)

The dependence of DPDs on $\boldsymbol{y}$ can hence not be directly inferred
from experimental observables. If $\boldsymbol{y}$ is small, one can use
perturbation theory to compute $F_{a_{1}a_{2}}$ in terms of PDFs and splitting
functions Diehl:2011yj ; Diehl:2019rdh . By contrast, for large distances the
$\boldsymbol{y}$ dependence is fully non-perturbative. Lattice studies can
give information about this dependence, whose knowledge is crucial for
computing double parton scattering cross sections.

Both unpolarised and polarised DPDs can exhibit correlations in their
dependence on $x_{1}$, $x_{2}$ and $\boldsymbol{y}$. We cannot address this
aspect in our present study, because the matrix elements we compute are
related to the lowest Mellin moments of DPDs, i.e. their integrals over both
$x_{1}$ and $x_{2}$. In principle, one could investigate higher Mellin
moments, i.e. integrals weighted with powers of $x_{1}$ and $x_{2}$. This
would require extending the set of currents in (9) to currents that involve
covariant derivatives and is beyond the scope of the present work.

Phenomenological analyses often make the assumption that in unpolarised DPDs
the two partons are independent of each other. This gives the relation

$\displaystyle F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle\overset{?}{=}\int
d^{2}\boldsymbol{b}\;f_{a_{1}}(x_{1},\boldsymbol{b}+\boldsymbol{y})\,f_{a_{2}}(x_{2},\boldsymbol{b})\,,$
(6)

where $f_{a}(x,\boldsymbol{b})$ is an unpolarised impact parameter dependent
single parton distribution. The question mark above the equal sign in (6)
indicates that this is a hypothesis. Our lattice study allows us to test this
indirectly in two different ways, as discussed in sections 2.4, 4.4 and 5.4.

A related but different simplifying assumption is that unpolarised DPDs can be
written as

$\displaystyle F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle\overset{?}{=}f_{a_{1}}(x_{1})\,f_{a_{2}}(x_{2})\,G(\boldsymbol{y})\,,$
(7)

where $f_{a}(x)$ denotes a standard PDF and $G(\boldsymbol{y})$ is a factor
describing the dependence on the transverse parton distance. This assumption
leads to the so-called “pocket formula”, which expresses double parton
scattering cross sections in terms of the cross sections for each single
scattering and a universal factor $\sigma_{\text{eff}}^{-1}=\int
d^{2}\boldsymbol{y}\;[\mskip 1.5muG(\boldsymbol{y})\mskip 1.5mu]^{2}$. Whilst
our study cannot address the factorisation between the $x_{1}$, $x_{2}$ and
$\boldsymbol{y}$ dependence assumed in (7), we can investigate the assumption
that the $\boldsymbol{y}$ dependence is the same for all parton combinations
$(a_{1},a_{2})$ in a given hadron. We will do this in section 4.5.

###  2.2 Matrix elements of local currents

The matrix element (2.1) involves fields at light-like distances and is hence
not suitable for direct evaluation on a Euclidean lattice. What we can study
in Euclidean space-time are the matrix elements

$\displaystyle M^{\mu_{1}\cdots\mu_{2}\cdots}_{q_{1}q_{2},i_{1}i_{2}}(p,y)$
$\displaystyle=\bra{h(p)}J^{\mu_{1}\cdots}_{q_{1},i_{1}}(y)\,J^{\mu_{2}\cdots}_{q_{2},i_{2}}(0)\ket{h(p)}\,,$
(8)

where as in (2.1) a polarisation average is understood if the hadron $h$
carries spin. The local currents $J_{q,i}^{\mu\cdots}$ we will consider here
are

$\displaystyle J_{q,V}^{\mu}(y)$ $\displaystyle=\bar{q}(y)\mskip
1.5mu\gamma^{\mu}\mskip 1.5muq(y)\,,$ $\displaystyle J_{q,A}^{\mu}(y)$
$\displaystyle=\bar{q}(y)\mskip 1.5mu\gamma^{\mu}\gamma_{5}\,q(y)\,,$
$\displaystyle J_{q,T}^{\mu\nu}(y)$ $\displaystyle=\bar{q}(y)\mskip
1.5mu\sigma^{\mu\nu}\mskip 1.5muq(y)\,.$ (9)

For spacelike distances $y$, which we assume throughout this work, the two
currents in (8) commute, so that one has

$\displaystyle M^{\mu_{1}\cdots\mu_{2}\cdots}_{q_{1}q_{2},i_{1}i_{2}}(p,y)$
$\displaystyle=M^{\mu_{2}\cdots\mu_{1}\cdots}_{q_{2}q_{1},i_{2}i_{1}}(p,-y)\,.$
(10)

Together with the fact that the currents in (9) are Hermitian, it follows that
the matrix elements (8) are real valued.

The currents transform under charge conjugation ($C$) and under the
combination of parity and time reversal ($PT$) as

$\displaystyle J_{q,i}^{\mu\cdots}(y)$
$\displaystyle\underset{C}{\to}\eta_{\textit{\tiny
C}}^{i}\,J_{q,i}^{\mu\cdots}(y)\,,$ $\displaystyle J_{q,i}^{\mu\cdots}(y)$
$\displaystyle\underset{PT}{\to}\eta_{\textit{\tiny
PT}}^{i}\,J_{q,i}^{\mu\cdots}(-y)$ (11)

with sign factors

$\displaystyle\eta_{\textit{\tiny C}}^{i}$
$\displaystyle=+1\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ i=A\,,$
$\displaystyle\eta_{\textit{\tiny C}}^{i}$
$\displaystyle=-1\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ i=V,T$ (12) and
$\displaystyle\eta_{\textit{\tiny PT}}^{i}$
$\displaystyle=+1\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ i=V\,,$
$\displaystyle\eta_{\textit{\tiny PT}}^{i}$
$\displaystyle=-1\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ i=A,T\,.$ (13)

The combination of a parity and time reversal transformation gives

$\displaystyle M^{\mu_{1}\cdots\mu_{2}\cdots}_{q_{1}q_{2},i_{1}i_{2}}(p,y)$
$\displaystyle=\eta_{\textit{\tiny PT}}^{i_{1}}\,\eta_{\textit{\tiny
PT}}^{i_{2}}\,M^{\mu_{1}\cdots\mu_{2}\cdots}_{q_{1}q_{2},i_{1}i_{2}}(p,-y)$
(14)

and thus relates the matrix elements for $y$ and $-y$.

#### Symmetry relations for pion matrix elements.

For pion matrix elements, one has additional relations due to charge
conjugation and isospin invariance. For $\eta_{\textit{\tiny
C}}^{i_{1}}\,\eta_{\textit{\tiny C}}^{i_{2}}=1$, which is the case for all
current combinations considered in our work, one has

$\displaystyle M_{q_{1}q_{2}}(y,p)\mskip 1.5mu\big{|}_{\pi^{+}}$
$\displaystyle=M_{q_{1}q_{2}}(y,p)\mskip 1.5mu\big{|}_{\pi^{-}}\,,$ (15)

where we indicated for which hadron the matrix element is taken but for
brevity omitted the Lorentz indices and the labels $i_{1},i_{2}$ specifying
the currents. Still for $\eta_{\textit{\tiny C}}^{i_{1}}\,\eta_{\textit{\tiny
C}}^{i_{2}}=1$, one finds

$\displaystyle M_{uu}(y,p)\mskip 1.5mu\big{|}_{\pi^{c}}$
$\displaystyle=M_{dd}(y,p)\mskip 1.5mu\big{|}_{\pi^{c}}\,,$ $\displaystyle
M_{ud}(y,p)\mskip 1.5mu\big{|}_{\pi^{c}}$ $\displaystyle=M_{du}(y,p)\mskip
1.5mu\big{|}_{\pi^{c}}$ (16)

for $c=+,-,0$, as well as

$\displaystyle M_{ud}(y,p)\mskip 1.5mu\big{|}_{\pi^{+}}+M_{uu}(y,p)\mskip
1.5mu\big{|}_{\pi^{+}}$ $\displaystyle=M_{ud}(y,p)\mskip
1.5mu\big{|}_{\pi^{0}}+M_{uu}(y,p)\mskip 1.5mu\big{|}_{\pi^{0}}\,.$ (17)

A derivation of these relations can be found in (Bali:2018nde, , section 2.1).

#### Tensor decomposition and extraction of twist-two functions.

The matrix elements in (8) are related to the lowest Mellin moments of DPDs as

$\displaystyle\int_{-\infty}^{\infty}dy^{-}\,M_{q_{1}q_{2},VV}^{++}(p,y)\,\Big{|}_{y^{+}=0,\,\boldsymbol{p}=\boldsymbol{0}}$
$\displaystyle=2p^{+}\mskip 1.5muI_{q_{1}q_{2}}(y^{2})\,,$
$\displaystyle\int_{-\infty}^{\infty}dy^{-}\mskip
1.5muM_{q_{1}q_{2},AA}^{++}(p,y)\,\Big{|}_{y^{+}=0,\,\boldsymbol{p}=\boldsymbol{0}}$
$\displaystyle=2p^{+}\mskip 1.5muI_{\Delta q_{1}\Delta q_{2}}(y^{2})\,,$
$\displaystyle\int_{-\infty}^{\infty}dy^{-}\mskip
1.5muM_{q_{1}q_{2},TV}^{k_{1}++}(p,y)\,\Big{|}_{y^{+}=0,\,\boldsymbol{p}=\boldsymbol{0}}$
$\displaystyle=2p^{+}\mskip 1.5mu\boldsymbol{y}^{k_{1}}mI_{\delta
q_{1}q_{2}}(y^{2})\,,$ $\displaystyle\int_{-\infty}^{\infty}dy^{-}\mskip
1.5muM_{q_{1}q_{2},VT}^{+k_{2}+}(p,y)\,\Big{|}_{y^{+}=0,\,\boldsymbol{p}=\boldsymbol{0}}$
$\displaystyle=2p^{+}\mskip 1.5mu\boldsymbol{y}^{k_{2}}mI_{q_{1}\delta
q_{2}}(y^{2})\,,$ $\displaystyle\int_{-\infty}^{\infty}dy^{-}\mskip
1.5muM_{q_{1}q_{2},TT}^{k_{1}+k_{2}+}(p,y)\,\Big{|}_{y^{+}=0,\,\boldsymbol{p}=\boldsymbol{0}}$
$\displaystyle=2p^{+}\mskip 1.5mu\bigl{[}\mskip
1.5mu\delta^{k_{1}k_{2}}\,I_{\delta q_{1}\delta q_{2}}(y^{2})$
$\displaystyle\qquad\quad-\bigl{(}2\boldsymbol{y}^{k_{1}}\boldsymbol{y}^{k_{2}}-\delta^{k_{1}k_{2}}\boldsymbol{y}^{2}\bigr{)}m^{2}\mskip
1.5muI^{t}_{\delta q_{1}\delta q_{2}}(y^{2})\bigr{]}$ (18)

with the Mellin moments given by

$\displaystyle I_{a_{1}a_{2}}(y^{2})$
$\displaystyle=\int_{-1}^{1}dx_{1}\int_{-1}^{1}dx_{2}\;f_{a_{1}a_{2}}(x_{1},x_{2},y^{2})$
$\displaystyle=\int_{0}^{1}dx_{1}\int_{0}^{1}dx_{2}\,\Bigl{[}f_{a_{1}a_{2}}(x_{1},x_{2},y^{2})+\eta_{\textit{\tiny
C}}^{i_{1}}f_{\bar{a}_{1}a_{2}}(x_{1},x_{2},y^{2})$
$\displaystyle\qquad\qquad\qquad\qquad+\eta_{\textit{\tiny
C}}^{i_{2}}\,f_{a_{1}\bar{a}_{2}}(x_{1},x_{2},y^{2})+\eta_{\textit{\tiny
C}}^{i_{1}}\,\eta_{\textit{\tiny
C}}^{i_{2}}\,f_{\bar{a}_{1}\bar{a}_{2}}(x_{1},x_{2},y^{2})\Bigr{]}\,.$ (19)

Here $i_{1}$ and $i_{2}$ refer to the currents in the matrix elements on the
l.h.s. of (2.2). An analogous relation holds between $I^{t}$ and the lowest
moment of $f^{t}$. The relations (2.2) extend the well-known connection
between the Mellin moments of PDFs and the matrix elements of a single local
current to the case of two partons.

In analogy to the case of PDFs, the matrix element (2.1) defining a DPD has
support for both positive and negative $x_{1}$ and $x_{2}$, with positive
$x_{i}$ corresponding to a parton $a_{i}$ and negative $x_{i}$ to its
antiparton $\bar{a}_{i}$. On the r.h.s. of (2.2), we have limited the
integration region to positive momentum fractions. Note that if $a_{1}$ and
$a_{2}$ are quarks and if $i=V$ or $T$ (but not $A$), then the quark-antiquark
distributions on the r.h.s. enter with a minus sign. This is of special
importance for distributions in a pion, whose valence Fock state consists of a
quark and an antiquark. Relations analogous to (2.2) exist for higher Mellin
moments in $x_{1}$ and $x_{2}$ and involve local currents with covariant
derivatives Diehl:2011yj , as is the case for PDFs.

Contrary to $\Gamma_{\delta q}^{j}$ in (3), the tensor current
$J_{q,T}^{\mu\nu}$ in (9) is defined without $\gamma_{5}$. As a consequence,
the vector indices $k_{1}$ and $k_{2}$ in (2.2) do _not_ give the transverse
quark spin direction but the transverse quark spin direction rotated by
$+90^{\circ}$ in the $x-y$ plane. This follows from the relation
$i\sigma^{j+}\gamma_{5}=\epsilon^{jk}\sigma^{k+}$.

The relations (2.2) still refer to Minkowski space, because they involve plus-
components. To make contact with matrix elements evaluated in Euclidean space,
we decompose the matrix elements (8) in terms of basis tensors and of Lorentz
invariant functions $A$, $B$, $C$, $D$, $E$ that depend on
$y^{2}=y^{\mu}y_{\mu}$ and $py=p^{\mu}y_{\mu}$. We write

$\displaystyle\tfrac{1}{2}\bigl{[}M^{\mu\nu}_{q_{1}q_{2},VV}(p,y)+M^{\nu\mu}_{q_{1}q_{2},VV}(p,y)\bigr{]}$
$\displaystyle=t_{VV,A}^{\mu\nu}\,A_{q_{1}q_{2}}+t_{VV,B}^{\mu\nu}\,m^{2}\mskip
1.5muB_{q_{1}q_{2}}+t_{VV,C}^{\mu\nu}\,m^{4}\mskip 1.5muC_{q_{1}q_{2}}$
$\displaystyle\quad+t_{VV,D}^{\mu\nu}\,m^{2}\mskip 1.5muD_{q_{1}q_{2}}\,,$
$\displaystyle\operatorname{T}M^{\mu\nu\rho}_{q_{1}q_{2},TV}(p,y)$
$\displaystyle=u_{TV,A}^{\mu\nu\rho}\,m\mskip 1.5muA_{\delta
q_{1}q_{2}}+u_{TV,B}^{\mu\nu\rho}\,m^{3}\mskip 1.5muB_{\delta q_{1}q_{2}}\,,$
$\displaystyle\tfrac{1}{2}\,\bigl{[}M^{\mu\nu\rho\sigma}_{q_{1}q_{2},TT}(p,y)+M^{\rho\sigma\mu\nu}_{q_{1}q_{2},TT}(p,y)\bigr{]}$
$\displaystyle=u_{TT,A}^{\mu\nu\rho\sigma}\,A_{\delta q_{1}\delta
q_{2}}+u_{TT,B}^{\mu\nu\rho\sigma}\,m^{2}\mskip 1.5muB_{\delta q_{1}\delta
q_{2}}+u_{TT,C}^{\mu\nu\rho\sigma}\,m^{2}\mskip 1.5muC_{\delta q_{1}\delta
q_{2}}$ $\displaystyle\quad+u_{TT,D}^{\mu\nu\rho\sigma}\,m^{4}\mskip
1.5muD_{\delta q_{1}\delta q_{2}}+u_{TT,E}^{\mu\nu\rho\sigma}\,m^{2}\mskip
1.5muE_{\delta q_{1}\delta q_{2}}\,.$ (20)

For the operator combination $TV$, we subtract trace terms according to

$\displaystyle\operatorname{T}u^{\mu\nu\rho}$
$\displaystyle=u^{\mu\nu\rho}+\tfrac{1}{3}\mskip
1.5mu\bigl{(}g^{\mu\rho}u^{\nu\alpha}{}_{\\!\alpha}-g^{\nu\rho}u^{\mu\alpha}{}_{\\!\alpha}\bigr{)}\,,$
(21)

where it is understood that $u^{\mu\nu\rho}$ is antisymmetric in $\mu$ and
$\nu$. The decomposition for $M_{q_{1}q_{2},AA}$ has the same form as the one
for $M_{q_{1}q_{2},VV}$, involving the same basis tensors but different
invariant functions $A_{\Delta q_{1}\Delta q_{2}}$, …, $D_{\Delta q_{1}\Delta
q_{2}}$. The decomposition for $M_{q_{1}q_{2},VT}$ is like the one for
$M_{q_{1}q_{2},TV}$ with an appropriate change in the role of the Lorentz
indices. In the following, we will not discuss the combination $VT$ any
further, because it can be traded for $TV$ using the relation (10). The basis
tensors are chosen as

$\displaystyle t_{VV,A}^{\mu\nu}$
$\displaystyle=2p^{\mu}p^{\nu}-{\textstyle\frac{1}{2}}\mskip
1.5mug^{\mu\nu}p^{2}\,,$ $\displaystyle t_{VV,B}^{\mu\nu}$
$\displaystyle=p^{\mu}y^{\nu}+p^{\nu}y^{\mu}-{\textstyle\frac{1}{2}}\mskip
1.5mug^{\mu\nu}py\,,$ $\displaystyle t_{VV,C}^{\mu\nu}$
$\displaystyle=2y^{\mu}y^{\nu}-{\textstyle\frac{1}{2}}\mskip
1.5mug^{\mu\nu}y^{2}\,,$ $\displaystyle t_{VV,D}^{\mu\nu}$
$\displaystyle=g^{\mu\nu}\,,$ $\displaystyle u_{TV,A}^{\mu\nu\rho}$
$\displaystyle=2(y^{\mu}p^{\nu}-p^{\mu}y^{\nu})\mskip
1.5mup^{\rho}+\tfrac{2}{3}\mskip
1.5mu(g^{\mu\rho}y^{\nu}-g^{\nu\rho}y^{\mu})\mskip
1.5mup^{2}-\tfrac{2}{3}\mskip
1.5mu(g^{\mu\rho}p^{\nu}-g^{\nu\rho}p^{\mu})\mskip 1.5mupy\,,$ $\displaystyle
u_{TV,B}^{\mu\nu\rho}$ $\displaystyle=2(y^{\mu}p^{\nu}-p^{\mu}y^{\nu})\mskip
1.5muy^{\rho}+\tfrac{2}{3}\mskip
1.5mu(g^{\mu\rho}y^{\nu}-g^{\nu\rho}y^{\mu})\mskip 1.5mupy-\tfrac{2}{3}\mskip
1.5mu(g^{\mu\rho}p^{\nu}-g^{\nu\rho}p^{\mu})\mskip 1.5muy^{2}\,,$
$\displaystyle u_{TT,A}^{\mu\nu\rho\sigma}$ $\displaystyle=-2\mskip
1.5mu\bigl{(}g^{\mu\rho}p^{\nu}p^{\sigma}-g^{\mu\sigma}p^{\nu}p^{\rho})+\tfrac{1}{2}\mskip
1.5mu(g^{\mu\rho}g^{\nu\sigma}-g^{\mu\sigma}g^{\nu\rho})\mskip
1.5mup^{2}-\\{\mu\leftrightarrow\nu\\}\,,$ $\displaystyle
u_{TT,B}^{\mu\nu\rho\sigma}$
$\displaystyle=-y^{2}\,u_{TT,A}^{\mu\nu\rho\sigma}-4\mskip
1.5mu(y^{\mu}p^{\nu}-p^{\mu}y^{\nu})\mskip
1.5mu(y^{\rho}p^{\sigma}-p^{\rho}y^{\sigma})+\tfrac{2}{3}\mskip
1.5mu(g^{\mu\rho}g^{\nu\sigma}-g^{\mu\sigma}g^{\nu\rho})\mskip
1.5mu\bigl{[}p^{2}y^{2}-(py)^{2}\bigr{]}\,,$ $\displaystyle
u_{TT,C}^{\mu\nu\rho\sigma}$
$\displaystyle=-(g^{\mu\rho}p^{\nu}y^{\sigma}-g^{\mu\sigma}p^{\nu}y^{\rho}+g^{\mu\rho}y^{\nu}p^{\sigma}-g^{\mu\sigma}y^{\nu}p^{\rho})+\tfrac{1}{2}\mskip
1.5mu(g^{\mu\rho}g^{\nu\sigma}-g^{\mu\sigma}g^{\nu\rho})\mskip
1.5mupy-\\{\mu\leftrightarrow\nu\\}\,,$ $\displaystyle
u_{TT,D}^{\mu\nu\rho\sigma}$ $\displaystyle=-2\mskip
1.5mu\bigl{(}g^{\mu\rho}y^{\nu}y^{\sigma}-g^{\mu\sigma}y^{\nu}y^{\rho})+\tfrac{1}{2}\mskip
1.5mu(g^{\mu\rho}g^{\nu\sigma}-g^{\mu\sigma}g^{\nu\rho})\mskip
1.5muy^{2}-\\{\mu\leftrightarrow\nu\\}\,,$ $\displaystyle
u_{TT,E}^{\mu\nu\rho\sigma}$
$\displaystyle=g^{\mu\rho}g^{\nu\sigma}-g^{\mu\sigma}g^{\nu\rho}\,.$ (22)

The tensor components related to twist-two matrix elements can be identified
from the l.h.s. of (2.2), taking into account that $y^{+}=0$ and
$\boldsymbol{p}=\boldsymbol{0}$ in that equation. For the basis tensors, a
nonzero plus-component requires the vector $p$ on the r.h.s. of (2.2), whilst
a nonzero transverse component requires the vector $y$ or the metric tensor.
One thus finds that the invariant functions corresponding to operators of
twist two are $A_{q_{1}q_{2}}$, $A_{\Delta q_{1}\Delta q_{2}}$, $A_{\delta
q_{1}q_{2}}$, $A_{\delta q_{1}\delta q_{2}}$ and $B_{\delta q_{1}\delta
q_{2}}$. We will call them “twist-two functions” in the remainder of this
work. All of them are even functions of $py$ due to the symmetry relation
(14).

One can project out the invariant functions by multiplying the matrix elements
with suitable linear combinations of basis tensors. For the twist-two
functions, the relevant projections read

$\displaystyle A_{q_{1}q_{2}}$
$\displaystyle=\frac{1}{8N^{2}}\,\Bigr{\\{}3(y^{2})^{2}\,t_{VV,A}^{\mu\nu}-6y^{2}py\;t_{VV,B}^{\mu\nu}+\bigl{[}p^{2}y^{2}+2(py)^{2}\bigr{]}\mskip
1.5mut_{VV,C}^{\mu\nu}\Bigr{\\}}\bigl{[}M_{q_{1}q_{2},VV}\bigr{]}_{\mu\nu}\,,$
$\displaystyle m\mskip 1.5muA_{\delta q_{1}q_{2}}$
$\displaystyle=\frac{3}{16N^{2}}\,\Bigr{\\{}y^{2}\,u_{TV,A}^{\mu\nu\rho}-py\,u_{TV,B}^{\mu\nu\rho}\Bigr{\\}}\operatorname{T}\bigl{[}M_{q_{1}q_{2},TV}\bigr{]}_{\mu\nu\rho}\,,$
$\displaystyle A_{\delta q_{1}\delta q_{2}}$
$\displaystyle=\frac{1}{64N^{2}}\,\Bigr{\\{}3(y^{2})^{2}\,u_{TT,A}^{\mu\nu\rho\sigma}-6y^{2}py\;u_{TT,C}^{\mu\nu\rho\sigma}+\bigl{[}p^{2}y^{2}+2(py)^{2}\bigr{]}\mskip
1.5muu_{TT,D}^{\mu\nu\rho\sigma}\Bigr{\\}}\,\bigl{[}M_{q_{1}q_{2},TT}\bigr{]}_{\mu\nu\rho\sigma}\,,$
$\displaystyle m^{2}B_{\delta q_{1}\delta q_{2}}$
$\displaystyle=\frac{1}{64N^{2}}\Bigr{\\{}3\mskip
1.5muu_{TT,B}^{\mu\nu\rho\sigma}+6\mskip
1.5mupy\;u_{TT,C}^{\mu\nu\rho\sigma}-3p^{2}\,u_{TT,D}^{\mu\nu\rho\sigma}\Bigr{\\}}\,\bigl{[}M_{q_{1}q_{2},TT}\bigr{]}_{\mu\nu\rho\sigma}$
(23)

with a normalisation factor

$\displaystyle N=p^{2}y^{2}-(py)^{2}\,.$ (24)

For spacelike $y^{\mu}$, which we are interested in, one has $N<0$, so that
the projections are always well defined.

Using (2.2) and (2.2), one can derive the relation between Mellin moments of
DPDs and integrals of twist-two functions over $py$:

$\displaystyle I_{a_{1}a_{2}}(y^{2})$
$\displaystyle=\int_{-\infty}^{\infty}d(py)\,A_{a_{1}a_{2}}(py,y^{2})\,,$
$\displaystyle I^{t}_{\delta q_{1}\delta q_{2}}(y^{2})$
$\displaystyle=\int_{-\infty}^{\infty}d(py)\,B_{\delta q_{1}\delta
q_{2}}(py,y^{2})\,,$ (25)

where in the first line we have all combinations of $(a_{1},a_{2})$ that
appear on the r.h.s. of (2.2).

The matrix elements (8) can be evaluated in Euclidean space-time at $y^{0}$,
i.e. with the two current operators taken at equal Euclidean time. This
entails the important restriction

$\displaystyle(py)^{2}=(\vec{\mskip 0.5mup}\mskip 1.5mu\vec{\mskip
0.5muy}\mskip 1.5mu)^{2}\leq\vec{\mskip 0.5mup}\mskip 1.5mu^{2}\,\vec{\mskip
0.5muy}\mskip 1.5mu^{2}\,,$ (26)

where $\vec{\mskip 0.5muv}\mskip 1.5mu=(v^{1},v^{2},v^{3})$ denotes the
spatial components of a four-vector $v^{\mu}$. Since the range of accessible
hadron momenta $\vec{p}$ in a lattice calculation is finite, the range of the
variable $py$ is limited, and one cannot directly evaluate the integrals in
(2.2). In addition, one needs data for nonzero hadron momentum $\vec{\mskip
0.5mup}\mskip 1.5mu$ to access even a finite range in $py$.

We note that the restriction (26) also applies if one computes the Mellin
moments of transverse-momentum dependent single parton distributions (TMDs) on
the lattice Hagler:2009mb ; Musch:2010ka ; Yoon:2017qzo . In that case,
$y^{\mu}$ is the distance between the quark and the antiquark field in the
matrix elements that define the distributions. The same holds for lattice
studies of single parton distributions in $x$ space. There has been an
enormous amount of activity in this area in recent years; we can only cite a
few papers here Ji:2013dva ; Ji:2014gla ; Radyushkin:2017cyf ; Orginos:2017kos
; Ji:2018hvs ; Braun:2018brg ; Alexandrou:2019lfo ; Ebert:2019okf ; Ji:2019ewn
and refer to the recent reviews Cichy:2018mum ; Ji:2020ect  for an extended
bibliography.

###  2.3 Skewed double parton distributions

Together with the restriction (26), the necessity to perform an integral over
all $py$ in (2.2) presents a significant complication for relating matrix
elements calculated on a Euclidean lattice with the Mellin moments of DPDs.
This prompts us to extend the theoretical framework in such a way that we can
discuss the physical meaning of the twist-two functions $A_{a_{1}a_{2}}$ and
$B_{\delta q_{1}\delta q_{2}}$ at a given value of $py$.

To this end, we introduce _skewed_ double parton distributions111The term
“skewed” refers to the parton momenta here, whilst the hadron momentum is the
same in the bra and ket vector of (2.3). This is different from “skewed parton
distributions”, now commonly called “generalised parton distributions”, which
involve two instead of four parton fields, such that there is an asymmetry
both in the parton and in the hadron momenta.

$\displaystyle F_{a_{1}a_{2}}(x_{1},x_{2},\zeta,\boldsymbol{y})=2p^{+}\int
dy^{-}e^{-i\zeta y^{-}p^{+}}$
$\displaystyle\int\frac{dz^{-}_{1}}{2\pi}\,\frac{dz^{-}_{2}}{2\pi}\,e^{i\mskip
1.5mu(x_{1}z_{1}^{-}+x_{2}z_{2}^{-})\mskip 1.5mup^{+}}$
$\displaystyle\quad\times\bra{h(p)}\mathcal{O}_{a_{1}}(y,z_{1})\,\mathcal{O}_{a_{2}}(0,z_{2})\ket{h(p)}\,.$
(27)

Compared with the definition (2.1) of ordinary DPDs, we have an additional
exponential $e^{-i\zeta p^{+}y^{-}}$ here. As a consequence, the partons
created or annihilated by the fields $\bar{q}$ and $q$ in
$\mathcal{O}_{a_{1}}$ and $\mathcal{O}_{a_{2}}$ have different longitudinal
momentum fractions. A sketch is given in figure 1 for $(a_{1},a_{2})=(u,d)$
and the case where $x_{1}-{\textstyle\frac{1}{2}}\zeta$,
$x_{1}+{\textstyle\frac{1}{2}}\zeta$, $x_{2}-{\textstyle\frac{1}{2}}\zeta$ and
$x_{2}+{\textstyle\frac{1}{2}}\zeta$ are all positive. If
$x_{1}-{\textstyle\frac{1}{2}}\zeta$ becomes negative, the $u$ quark in the
wave function of $\ket{h}$ becomes an antiquark $\bar{u}$ with momentum
fraction $-x_{1}+{\textstyle\frac{1}{2}}\zeta$ in the wave function of
$\bra{h}$. Corresponding statements hold for
$x_{1}+{\textstyle\frac{1}{2}}\zeta$, $x_{2}-{\textstyle\frac{1}{2}}\zeta$ and
$x_{2}+{\textstyle\frac{1}{2}}\zeta$.

Figure 1:  Graphical representation of a skewed DPD for quark flavours $u$ and
$d$ in the hadron $h$. The configuration shown is for the case where all
momentum fractions given at the top of the graph are positive.

For nonzero $\zeta$, the distributions (2.3) do not appear in cross sections
for double parton scattering, but they may be regarded as a rather
straightforward extension of the DPD concept. Let us take a closer look at
some of their properties. The support region of the matrix element (2.3) in
the momentum fraction arguments is the same as if all four parton fields were
at the same transverse position. In that case, we would have a collinear
twist-four distribution. The support properties of these distributions were
derived in Jaffe:1983hp , and the argument given there does not depend on the
transverse position arguments of the parton fields. The result given in
Jaffe:1983hp  is equivalent to the interpretation of
$x_{1}-{\textstyle\frac{1}{2}}\zeta$, $x_{1}+{\textstyle\frac{1}{2}}\zeta$,
$x_{2}-{\textstyle\frac{1}{2}}\zeta$ and $x_{2}+{\textstyle\frac{1}{2}}\zeta$
as positive or negative momentum fractions, as described in the previous
paragraph. For nonzero $\zeta$ there are hence different regions, in which one
has either 1, 2 or 3 partons in the wave function of $\ket{h}$. With the
constraints that the partons in the wave function of $\ket{h}$ must carry the
same total longitudinal momentum as those in the wave function of $\bra{h}$,
and that this cannot be larger than the longitudinal hadron momentum, one
obtains the constraint

$-1\leq\zeta\leq 1$ (28)

and the support region for $(x_{1},x_{2})$ shown in figure 2. For $\zeta=0$
this region becomes a square with corners $(0,\pm 1)$ and $(\pm 1,0)$, whereas
for $\zeta=\pm 1$ it becomes a square with corners
$(\pm{\textstyle\frac{1}{2}},\pm{\textstyle\frac{1}{2}})$.

Figure 2:  Support region of the distribution
$F_{ud}(x_{1},x_{2},\zeta,\boldsymbol{y})$ in the momentum fraction arguments.
The notation $d\mskip 1.5mu|du\bar{u}$ means that one has one $d$ quark in the
wave function of $\ket{h}$ and $du\bar{u}$ in the wave function of $\bra{h}$.
In both panels, the triangle for the region $ud|du$ has the corners
$\bigl{(}\mskip
1.5mu{\textstyle\frac{1}{2}}|\zeta|,{\textstyle\frac{1}{2}}|\zeta|\mskip
1.5mu\bigr{)}$, $\bigl{(}\mskip
1.5mu{\textstyle\frac{1}{2}}|\zeta|,1-{\textstyle\frac{1}{2}}|\zeta|\mskip
1.5mu\bigr{)}$ and $\bigl{(}\mskip
1.5mu1-{\textstyle\frac{1}{2}}|\zeta|,{\textstyle\frac{1}{2}}|\zeta|\mskip
1.5mu\bigr{)}$. Notice that the parton configuration in each of the four
triangles is the same for positive and negative $\zeta$, whereas the
configuration in each of the squares is different.

Using $PT$ symmetry, one finds that

$\displaystyle F_{a_{1}a_{2}}(x_{1},x_{2},\zeta,\boldsymbol{y})$
$\displaystyle=\eta_{\textit{\tiny PT}}^{i_{1}}\,\eta_{\textit{\tiny
PT}}^{i_{2}}\,F_{a_{1}a_{2}}(x_{1},x_{2},-\zeta,-\boldsymbol{y})\,,$ (29)

where $\eta_{\textit{\tiny PT}}^{i}=+1$ for an unpolarised parton and
$\eta_{\textit{\tiny PT}}^{i}=-1$ for a polarised one. The skewed DPDs can be
decomposed in terms of scalar distributions as in (2.1), with the
distributions on both sides depending additionally on $\zeta$. The symmetry
property (29) then implies

$\displaystyle f_{a_{1}a_{2}}(x_{1},x_{2},\zeta,y^{2})$
$\displaystyle=f_{a_{1}a_{2}}(x_{1},x_{2},-\zeta,y^{2})$ (30)

and an analogous relation for $f^{t}$.

#### Mellin moments.

We define the lowest Mellin moments of skewed DPDs as

$\displaystyle I_{a_{1}a_{2}}(y^{2},\zeta)$
$\displaystyle=\int_{-1}^{1}dx_{1}\int_{-1}^{1}dx_{2}\;f_{a_{1}a_{2}}(x_{1},x_{2},\zeta,y^{2})$
(31)

and likewise for $f^{t}$, where the integration region in $x_{1},x_{2}$
follows from figure 2. The moments are nonzero for $\zeta$ in the interval
$[-1,1]$. The generalisation of (2.2) to nonzero $\zeta$ reads

$\displaystyle I_{a_{1}a_{2}}(y^{2},\zeta)$
$\displaystyle=\int_{-\infty}^{\infty}d(py)\,e^{-i\zeta
py}\,A_{a_{1}a_{2}}(y^{2},py)\,,$ (32)

which can readily be inverted for the function $A_{a_{1}a_{2}}(y^{2},py)$. In
particular, one finds

$\displaystyle
A_{a_{1}a_{2}}(y^{2},py=0)=\frac{1}{\pi}\int_{0}^{1}d\zeta\,I_{a_{1}a_{2}}(y^{2},\zeta)\,,$
(33)

where we have used the symmetry relation (30) to reduce the integration region
to positive $\zeta$. Rather than the Mellin moment of a DPD, a twist-two
function at $py=0$ is thus the average of the Mellin moment of a skewed DPD
over the skewness parameter $\zeta$.

Quantities that characterise the $\zeta$ dependence of
$I_{a_{1}a_{2}}(y^{2},\zeta)$ are the even moments in $\zeta$,

$\displaystyle\langle\zeta^{2m}\rangle_{a_{1}a_{2}}(y^{2})$
$\displaystyle=\frac{\int_{-1}^{1}d\zeta\;\zeta^{2m}\,I_{a_{1}a_{2}}(y^{2},\zeta)}{\int_{-1}^{1}d\zeta\,I_{a_{1}a_{2}}(y^{2},\zeta)}=\Biggl{[}\mskip
1.5mu\frac{(-1)^{m}}{A_{a_{1}a_{2}}(y^{2},py)}\frac{\partial^{2m}A_{a_{1}a_{2}}(y^{2},py)}{(\partial\mskip
1.5mupy)^{2m}}\mskip 1.5mu\Biggr{]}_{py=0}\,.$ (34)

Odd moments $\langle\zeta^{2m+1}\rangle$ are zero because of the symmetry
(30). To compute the moments $\langle\zeta^{2m}\rangle$, one needs
$A_{a_{1}a_{2}}(y^{2},py)$ in the vicinity of $py=0$. According to (26), this
can be evaluated from Euclidean data with nonzero hadron momentum $\vec{\mskip
0.5mup}\mskip 1.5mu$.

Relations analogous to (32) to (34) can be written down for $I^{t}_{\delta
q_{1}\delta q_{2}}$ and $B_{\delta q_{1}\delta q_{2}}$ in the place of
$I_{a_{1}a_{2}}$ and $A_{a_{1}a_{2}}$.

###  2.4 Factorisation hypotheses

We now discuss how the factorisation hypothesis (6) for DPDs can be formulated
at the level of Mellin moments and twist-two functions. At this point, we
specialise to the case where the hadron $h$ is a $\pi^{+}$. This avoids
complications due to the proton spin, which are discussed in (Diehl:2011yj, ,
section 4.3.1).

Let us take the lowest Mellin moment in $x_{1}$ and $x_{2}$ of (6). The Mellin
moment of an unpolarised impact parameter dependent parton distribution is

$\displaystyle\int_{-1}^{1}dx\,f_{q}(x,\boldsymbol{b})$
$\displaystyle=\int\frac{d^{2}\boldsymbol{\Delta}}{(2\pi)^{2}}\;e^{-i\boldsymbol{b}\boldsymbol{\Delta}}\,F_{q,V}(-\boldsymbol{\Delta}^{2})\,,$
(35)

where $F_{q,V}(t)$ is the form factor of the vector current

$\displaystyle\bra{\pi^{+}(p^{\prime})}J_{q,V}^{\mu}(0)\ket{\pi^{+}(p)}$
$\displaystyle=(p+p^{\prime})^{\mu}\,F_{q,V}(t)$ $\displaystyle\text{with }t$
$\displaystyle=(p-p^{\prime})^{2}\,.$ (36)

We then obtain from (6)

$\displaystyle I_{ud}(-\boldsymbol{y}^{2})$
$\displaystyle\overset{?}{=}\int\frac{d^{2}\boldsymbol{\Delta}}{(2\pi)^{2}}\;e^{-i\boldsymbol{y}\boldsymbol{\Delta}}\,F_{u,V}(-\boldsymbol{\Delta}^{2})\,F_{d,V}(-\boldsymbol{\Delta}^{2})\,.$
(37)

We note that thanks to isospin invariance, one has $F_{u,V}=-F_{d,V}$. As this
is not essential in the present context, we will not use it here.

Since one cannot directly determine $I_{ud}(-\boldsymbol{y}^{2})$ from
Euclidean correlation functions, one cannot directly test (37) with lattice
data. We therefore derive an analogous relation for the twist-two function
$A_{ud}(y^{2},py)$ at $py=0$.

We recall from Diehl:2011yj  that (6) can be obtained by inserting a complete
set of intermediate states between the operators
$\mathcal{O}_{a_{1}}(y,z_{1})$ and $\mathcal{O}_{a_{2}}(0,z_{2})$ in the DPD
definition (2.1) and then _assuming_ that the dominant term in this sum is the
ground state. Following exactly the same steps for the skewed DPD (2.3), one
obtains

$\displaystyle
F_{ud}(x_{1},x_{2},\zeta,\boldsymbol{y})\overset{?}{=}\int\frac{d^{2}\boldsymbol{\Delta}}{(2\pi)^{2}}\;e^{-i\boldsymbol{y}\boldsymbol{\Delta}}\;\frac{1}{1-\zeta}\;$
$\displaystyle H_{u}\biggl{[}\frac{2\mskip
1.5mux_{1}}{2-\zeta},\frac{\zeta}{2-\zeta},t(\boldsymbol{\Delta}^{2},\zeta)\mskip
1.5mu\biggr{]}\,$ $\displaystyle\times$
$\displaystyle\,H_{d}\biggl{[}\frac{2\mskip
1.5mux_{2}}{2-\zeta},\frac{\zeta}{2-\zeta},t(\boldsymbol{\Delta}^{2},\zeta)\mskip
1.5mu\biggr{]}$ (38)

with

$\displaystyle t(\boldsymbol{\Delta}^{2},\zeta)$
$\displaystyle=-\frac{\zeta^{2}m^{2}+\boldsymbol{\Delta}^{2}}{1-\zeta}\,.$
(39)

Here $H_{q}(x,\xi,t)$ is the generalised parton distribution (GPD) for
unpolarised quarks in a pion; its definition can be found e.g. in
(Diehl:2003ny, , section 3.2). The momentum fraction arguments $x$ and $\xi$
of $H_{q}$ are defined in a symmetric way between the incoming and outgoing
hadron and parton momenta, with $x$ referring to the sum of parton momenta and
$\xi$ to their difference, and with momentum fractions normalised to the sum
of hadron momenta in the bra and the ket state. Both $x$ and $\xi$ are limited
to the interval $[-1,1]$. A pictorial representation of the GPDs that appear
on the r.h.s. of (2.4) is given in figure 3(a).

(a) (b)

Figure 3:  (a): Pictorial representation of the r.h.s. of the factorisation
hypothesis (2.4). This is obtained by inserting a full set of intermediate
states between the operators in the matrix element
$\bra{h}\mathcal{O}_{u}\mskip 1.5mu\mathcal{O}_{d}\ket{h}$ and then retaining
only the ground state. (b): The representation obtained when inserting the
full set of states after reordering the operators to
$\bra{h}\mathcal{O}_{d}\mskip 1.5mu\mathcal{O}_{u}\ket{h}$. All momentum
fractions refer to the hadron $h$ in the matrix element. We use form (a) for
$\zeta\geq 0$ and form (b) for $\zeta<0$.

At this point, we must critically examine the support properties of the two
sides of (2.4) in $x_{1}$ and $x_{2}$. The support of the l.h.s. is shown in
figure 2, whereas the one of the r.h.s. is the square delineated by
$-1+{\textstyle\frac{1}{2}}\zeta\leq x_{1,2}\leq
1-{\textstyle\frac{1}{2}}\zeta$ in the $(x_{1},x_{2})$ plane. For $\zeta\geq
0$, this misses the kinematic constraint $|x_{1}|+|x_{2}|\leq 1$ in $F_{ud}$,
whereas for $\zeta<0$ it is even larger.

In the matrix element (2.3), the order of the two operators can be
interchanged, because the respective fields are separate by spacelike
distances. In a schematic notation, we thus have $\bra{h}\mathcal{O}_{u}\mskip
1.5mu\mathcal{O}_{d}\ket{h}=\bra{h}\mathcal{O}_{d}\mskip
1.5mu\mathcal{O}_{u}\ket{h}$. If we insert a set of intermediate states in the
latter matrix element, we obtain (2.4) with $\zeta$ replaced by $-\zeta$ on
the r.h.s. This is represented in figure 3(b). In that case, the mismatch
between the support regions of the two sides is less bad for $\zeta\leq 0$
than for $\zeta>0$. We therefore retain (2.4) for $\zeta\geq 0$ and its
analogue with $\zeta\to-\zeta$ on the r.h.s. for $\zeta<0$. This also
satisfies the symmetry in $\zeta$ required by $\mathit{PT}$ invariance and
stated in (29), which is violated if one uses (2.4) for positive and negative
$\zeta$.

The mismatch of support properties just discussed also affects the case
$\zeta=0$ and is hence not special to the skewed kinematics we are considering
here. In fact, it is well known that the factorisation hypothesis (6) for DPDs
violates the momentum conservation constraint $x_{1}+x_{2}\leq 1$. From a
theoretical point of view, inserting a full set of intermediate states between
the operators in the DPD definition (2.1) or its skewed analogue (2.3) is of
course a legitimate manipulation, but we see that the restriction of this set
to the ground state leads to theoretical inconsistencies such as an incorrect
support region or the loss of a symmetry required by $\mathit{PT}$ invariance.
How the sum over all states manages to restore the correct properties is
difficult to understand in an intuitive manner. We note that a similar
observation was made in Jaffe:1983hp  when discussing the support properties
of PDFs and of higher-twist distributions.

Integrating both sides of (2.4) over their respective support regions in
$x_{1}$ and $x_{2}$ and using the sum rule $\int
dx\,H_{q}(x,\xi,t)=F_{q,V}(t)$, one obtains

$\displaystyle I_{ud}(-\boldsymbol{y}^{2},\zeta)$
$\displaystyle\overset{?}{=}\frac{(1-{\textstyle\frac{1}{2}}\zeta)^{2}}{1-\zeta}\int\frac{d^{2}\boldsymbol{\Delta}}{(2\pi)^{2}}\;e^{-i\boldsymbol{y}\boldsymbol{\Delta}}\,F_{u,V}\bigl{(}t(\boldsymbol{\Delta}^{2},\zeta)\bigr{)}\,F_{d,V}\bigl{(}t(\boldsymbol{\Delta}^{2},\zeta)\bigr{)}\,.$
(40)

Using this for $\zeta\geq 0$ and inserting it into (33), we obtain

$\displaystyle A_{ud}(-\boldsymbol{y}^{2},py=0)$
$\displaystyle\overset{?}{=}\frac{1}{\pi}\int_{0}^{1}d\zeta\,\frac{(1-{\textstyle\frac{1}{2}}\zeta)^{2}}{1-\zeta}\int\frac{d^{2}\boldsymbol{\Delta}}{(2\pi)^{2}}\;e^{-i\boldsymbol{y}\boldsymbol{\Delta}}\,F_{u,V}\bigl{(}t(\boldsymbol{\Delta}^{2},\zeta)\bigr{)}\,F_{d,V}\bigl{(}t(\boldsymbol{\Delta}^{2},\zeta)\bigr{)}\,.$
(41)

We note that (41) is expressed in terms of a two-dimensional vector
$\boldsymbol{y}$. This is different from the factorisation hypothesis we
derived in (Bali:2018nde, , section 5.3), which involved the zero-components
of currents and a three-dimensional vector $\vec{\mskip 0.5muy}\mskip 1.5mu$.

Note that the two hypotheses (41) and (37) are based on the same assumption
but are not equivalent to each other. Both are special cases of (40), obtained
by either setting $\zeta=0$ or by integrating over $\zeta$ from $0$ to $1$.
The assumption that the ground state dominates the sum over intermediate
states could be a better approximation in one or the other case. Using our
lattice results, we will investigate (41) in section 4.4 and (37) in section
5.4.

##  3 Lattice computation and lattice artefacts

We performed lattice simulations for the matrix elements (8) in a pion with
the currents given in (9). We set $y^{0}=0$, so that on the lattice the two
currents are inserted at the same Euclidean time, but with a spatial
separation $\vec{\mskip 0.5muy}\mskip 1.5mu$. We generated data both for zero
and nonzero pion momentum $\vec{\mskip 0.5mup}\mskip 1.5mu$. The lattice
techniques we employed are explained in detail in (Bali:2018nde, , section 3).
In the following, we recall only the basic steps described in that work and
then proceed to the specifics of our present analysis.

###  3.1 Lattice graphs

To evaluate the two-current matrix elements (8), we compute the four-point
correlation function of a pion source operator at Euclidean time $0$, a pion
sink operator at Euclidean time $t$, and the two currents $J_{i}$ and $J_{j}$
at Euclidean time $\tau$. The correlation function receives contributions from
a large number of Wick contractions, which are shown in figure 4. We will also
refer to these contractions as “lattice graphs” or simply as “graphs”.

Figure 4:  Lattice graphs for the correlation functions used to extract the
matrix elements (8) in a pion. The dependence on the pion momentum
$\vec{\mskip 0.5mup}\mskip 1.5mu$ is not indicated for brevity.
$\eta_{\textit{\tiny C}}^{i}$ denotes the charge conjugation parity of the
current $J_{i}$ and is defined in (12).

The relation between pion matrix elements and lattice graphs depends on the
product of $C$ parities of the currents. Omitting Lorentz indices and the
dependence on the pion momentum $p$, and using the shorthand notation

$\displaystyle C_{1}$ $\displaystyle=C_{1}^{ij}(y)\,,$ $\displaystyle C_{2}$
$\displaystyle={\textstyle\frac{1}{2}}\bigl{[}C_{2}^{ij}(y)+C_{2}^{ji}(-y)\bigr{]}\,,$
$\displaystyle A$
$\displaystyle={\textstyle\frac{1}{2}}\bigl{[}A^{ij}(y)+A^{ji}(-y)\bigr{]}\,,$
$\displaystyle S_{1}$
$\displaystyle={\textstyle\frac{1}{2}}\bigl{[}S_{1}^{ij}(y)+S_{1}^{ji}(-y)\bigr{]}\,,$
$\displaystyle S_{2}$ $\displaystyle=S_{2}^{ij}(y)\,,$ $\displaystyle D$
$\displaystyle=D^{ij}(y)$ (42)

for the graphs or their symmetrised combinations, we have

$\displaystyle M_{ud,ij}(y)\mskip 1.5mu\big{|}_{\pi^{+}}$
$\displaystyle=C_{1}+\bigl{[}2S_{1}+D\bigr{]}\,,$ $\displaystyle
M_{uu,ij}(y)\mskip 1.5mu\big{|}_{\pi^{+}}$
$\displaystyle=\bigl{[}2C_{2}+S_{2}\bigr{]}+\bigl{[}2S_{1}+D\bigr{]}\,,$
$\displaystyle M_{ud,ij}(y)\mskip 1.5mu\big{|}_{\pi^{0}}$
$\displaystyle=\bigl{[}2S_{1}+D\bigr{]}-A\,,$ $\displaystyle
M_{uu,ij}(y)\mskip 1.5mu\big{|}_{\pi^{0}}$
$\displaystyle=C_{1}+\bigl{[}2S_{1}+D\bigr{]}+\bigl{[}2C_{2}+S_{2}\bigr{]}+A$
(43)

for $\eta_{\textit{\tiny C}}^{i\phantom{j}}\eta_{\textit{\tiny C}}^{j}=+1$,
which is satisfied for all combinations of currents considered in the present
study. We note that this is no longer the case if one includes operators with
covariant derivatives (corresponding to higher Mellin moments). One readily
checks that (3.1) satisfies the general symmetry relation (17), as it must.

To compute the different graphs on the lattice, we use a variety of techniques
as detailed in (Bali:2018nde, , section 3.3). We make extensive use of
stochastic sources, and for graph $C_{2}$ we use a hopping parameter expansion
to reduce statistical noise for the propagation between the two currents.

For the disconnected graphs $S_{2}$ and $D$ we need to subtract vacuum
contributions, namely the product of a two-point correlation function of the
pion source and sink with a two-point correlation function of the two
currents. The latter corresponds to the vacuum expectation value
$\bra{0}J_{i}(y)J_{j}(0)\ket{0}$. The vacuum subtraction for the disconnected
graph $S_{1}$ involves $\bra{0}J_{i}(y)\ket{0}$ or $\bra{0}J_{j}(0)\ket{0}$,
which is zero because our currents carry Lorentz indices.

We anticipate that the doubly disconnected graph $D$ in general gives a good
signal for the four-point correlation function, but that there is a near-
perfect cancellation between this correlator and its vacuum subtraction term.
The result after subtraction is consistent with zero and has huge statistical
uncertainties compared with those of any other graph. We will hence not be
able to report useful results for graph $D$. Fortunately, we encounter no such
problem for graph $S_{2}$.

###  3.2 Lattice simulation and extraction of twist-two functions

We perform our simulations using the Wilson gauge action and $n_{F}=2$ mass
degenerate flavours of non-perturbatively improved Sheikholeslami-Wohlert (NPI
Wilson-clover) fermions. The gauge configurations were generated by the RQCD
and QCDSF collaborations. We use two gauge ensembles, whose parameters are
given in table 1. They have different spatial sizes, $L=32$ and $L=40$, which
allows us to study finite volume effects in section 3.5. Despite having data
for only a single lattice spacing, $a=0.071\operatorname{fm}$, we are also
able to investigate discretisation effects, as discussed in section 3.3.

ensemble | $\beta$ | $a\,[\operatorname{fm}]$ | $\kappa$ | $L^{3}\times T$ | $m_{\pi}\,[\operatorname{MeV}]$ | $Lm_{\pi}$ | $N_{\text{full}}$ | $N_{\text{used}}$
---|---|---|---|---|---|---|---|---
IV | 5.29 | 0.071 | 0.13632 | $32^{3}\times 64$ | $294.6\pm 1.4$ | 3.42 | 2023 | 960
V | 5.29 | 0.071 | 0.13632 | $40^{3}\times 64$ | $288.8\pm 1.1$ | 4.19 | 2025 | 984
Table 1:  Details of the gauge ensembles used in this analysis.
$N_{\text{full}}$ is the total number of available gauge configurations, and
$N_{\text{used}}$ is the number of configurations used in our simulations.
More detail can be found in Bali:2014gha ; Bali:2014nma .

For the ensemble with $L=40$, we performed simulations with different $\kappa$
values in the quark propagator:

light quarks: $\displaystyle\kappa$ $\displaystyle=0.13632\,,$ $\displaystyle
m_{\pi}$ $\displaystyle=293\operatorname{MeV}\,,$ strange:
$\displaystyle\kappa$ $\displaystyle=0.135616\,,$ $\displaystyle m_{\pi}$
$\displaystyle=691\operatorname{MeV}\,,$ charm: $\displaystyle\kappa$
$\displaystyle=0.125638\,,$ $\displaystyle m_{\pi}$
$\displaystyle=3018\operatorname{MeV}\,.$ (44)

Here “light quarks” refers to the $\kappa$ value used for simulating the sea
quarks, whereas the other two values correspond to the physical strange and
charm quark masses, as determined in Bali:2016lvx  and Bali:2017pdv  by tuning
the pseudoscalar ground state mass to $685.8\operatorname{MeV}$ in the first
case and the spin-averaged $S$-wave charmonium mass to
$3068.5\operatorname{MeV}$ in the second case. Since our simulations are
performed with an $n_{F}=2$ fermion action, the strange and charm quarks are
partially quenched.

The values of $m_{\pi}$ in (3.2) are obtained from exponential fits of the
pion two-point function. We quote them only for orientation and do not attempt
to quantify their errors. These masses are in reasonable agreement with the
value in table 1 for light quarks, and with the mass of the pseudoscalar
ground state quoted below (3.2) for strange quarks.

#### Pion matrix elements.

For all lattice graphs, we compute the correlation functions with zero three-
momentum $\vec{\mskip 0.5mup}\mskip 1.5mu$ of the pion. For the connected
graphs $C_{1}$ and $C_{2}$, we additionally have data with finite pion
momenta. These data are restricted to the $L=40$ lattice and to light quarks.
The pion momenta that can be realised on the lattice are given by

$\displaystyle\vec{\mskip 0.5mup}\mskip 1.5mu$
$\displaystyle=\frac{2\pi}{La}\,\vec{\mskip 0.5muP}\mskip 1.5mu\,,$ (45)

where the components of $\vec{\mskip 0.5muP}\mskip 1.5mu$ are integers and
$2\pi/(La)\approx 437\operatorname{MeV}$ in our case. For simplicity we write
$P=|\vec{\mskip 0.5muP}\mskip 1.5mu|$. Graph $C_{1}$ is computed for all 18
nonzero momenta with $P^{2}\leq 2$, for 6 momenta with $P^{2}=3$, and for one
momentum with $P^{2}=4$. For $C_{2}$, we have results for all 6 momenta with
$P=1$.

The distance between the pion source and sink in the correlation functions is
fixed to $t=15a\approx 1.07\operatorname{fm}$ as a default. To investigate the
influence of excited states, we also calculate graphs $C_{1}$, $C_{2}$ and
$S_{1}$ with $t=32a$. The matrix element (8) is extracted from the ratio
between the four-point correlation function around $\tau=t/2$ and the pion
two-point function. For graphs $C_{1}$ and $A$, we measure the $\tau$
dependence of the four-point function and fit to a plateau in the $\tau$
ranges specified in (Bali:2018nde, , equation (4.1)). The quality of the
corresponding plateaus is good for matrix elements that have a nonzero value
within statistical uncertainties. For the remaining graphs, we extract the
matrix element from data at $\tau/a=7$ and $8$ if $t/a=15$. For the $C_{2}$
and $S_{1}$ data with $t/a=32$, we use $\tau/a=16$. A comparison of data with
$t=15a$ and $t=32a$ is shown in section 3.4.

All lattice currents are converted to the $\overline{\text{MS}}$ scheme at the
renormalisation scale

$\displaystyle\mu$ $\displaystyle=2\operatorname{GeV}\,.$ (46)

As described in (Bali:2018nde, , section 3.4), this is done using a
combination of non-perturbative and perturbative renormalisation and includes
an estimate of the quark mass dependent order $a$ improvement term.

#### Invariant functions.

From the matrix elements (8), we determine the invariant functions for each
individual value of $\vec{\mskip 0.5muy}\mskip 1.5mu$ and $\vec{\mskip
0.5mup}\mskip 1.5mu$. This is done using a minimum $\chi^{2}$ fit of the data
for all tensor components to the decomposition (2.2). For invariant functions
of twist two, we also use the projector method (2.2). In both cases, the
statistical error of an invariant function at given $\vec{\mskip 0.5muy}\mskip
1.5mu$ and $\vec{\mskip 0.5mup}\mskip 1.5mu$ is computed using the jackknife
method. To eliminate autocorrelations, we take the number of jackknife samples
as $1/8$ times the number $N_{\text{used}}$ of gauge configurations given in
table 1.

For $P=0$, the twist-two functions extracted with one or the other method show
excellent agreement with each other and have statistical uncertainties of
almost the same size. For $P>0$, the values obtained with the projection
method have much larger statistical errors than those obtained with a fit and
provide only a very weak cross check. All data shown in the following are
obtained by the fit method, both for $P=0$ and $P>0$.

In the remainder of this section, we investigate the extent to which our data
are affected by lattice artefacts, largely following the corresponding studies
in (Bali:2018nde, , section 4). We only consider data with $py=0$ here,
because they have much smaller statistical errors than the data for $py\neq
0$. We will return to the case of nonzero $py$ in section 5.

When discussing twist-two functions extracted from the correlation functions
for particular lattice graphs, we will generically write $A_{qq}$, $A_{\Delta
q\Delta q}$, …, $B_{\delta q\delta q}$, without reference to specific quark
flavors $q_{1}$ and $q_{2}$. This is because the distinction between $u$ and
$d$ quarks in a pion only appears when lattice graphs are combined as
specified in (3.1).

###  3.3 Isotropy and boost invariance

The decomposition (2.2) of matrix elements in terms of basis tensors and
functions of $y^{2}$ and $py$ assumes Lorentz invariance and thus requires
both the continuum and the infinite-volume limit. If our lattice simulations
are sufficiently close to these limits, then the values of twist-two functions
extracted for individual points $\vec{\mskip 0.5muy}\mskip 1.5mu$ and
$\vec{\mskip 0.5mup}\mskip 1.5mu$ with $\vec{\mskip 0.5mup}\mskip
1.5mu\vec{\mskip 0.5muy}\mskip 1.5mu=0$ must not depend on the directions of
$\vec{\mskip 0.5muy}\mskip 1.5mu$ or $\vec{\mskip 0.5mup}\mskip 1.5mu$ or on
the size of $\vec{\mskip 0.5mup}\mskip 1.5mu$.

Let us test whether this is the case in our simulations for light quarks on
the lattice with $L=40$. We restrict our attention to graphs $C_{1}$ and
$C_{2}$, for which statistical errors are small enough to reveal the effects
of interest. For the sake of legibility, we henceforth write $y=|\vec{\mskip
0.5muy}\mskip 1.5mu|$ for the length of the spatial distance $\vec{\mskip
0.5muy}\mskip 1.5mu$ between the two currents. We continue to use $y^{2}$ and
$py$ to denote the products $y^{\mu}y_{\mu}$ and $p^{\mu}y_{\mu}$ of four-
vectors in Minkowski space. Since $y^{\mu}$ is always spacelike in our
context, this implies that $y^{2}<0$.

(a) $A_{qq}$, graph $C_{1}$ (b) $A_{\delta q\mskip 1.5muq}$, graph $C_{1}$

Figure 5:  Examples for the anisotropy in graph $C_{1}$ at large
$y=|\vec{\mskip 0.5muy}\mskip 1.5mu|$. The data shown are for $L=40$, zero
pion momentum, and light quarks. The results in this and all the following
figures are given in the $\overline{\mathrm{MS}}$ scheme at the scale
$\mu=2\operatorname{GeV}$. The error bars shown in the plots are statistical
and obtained with the jackknife method.

At large $y$ of order $La/2$, we see a clear anisotropy with a saw-tooth
pattern in all twist-two functions that have sufficiently small errors.
Examples are shown in figure 5. This pattern is expected on a lattice with
periodic boundary conditions and can be understood in terms of “mirror
images”. The same effect has been seen and discussed in previous lattice
studies of two-current correlators Burkardt:1994pw ; Alexandrou:2008ru ,
including our study in Bali:2018nde  that employed the same lattice data as
the present work. As shown in Burkardt:1994pw , the effect of mirror charges
at a given distance $y$ is smallest for points $\vec{\mskip 0.5muy}\mskip
1.5mu$ close to one of the space diagonals, i.e. the lines given by
$\vec{\mskip 0.5muz}\mskip 1.5mu=(z^{1},z^{2},z^{3})$ with
$|z^{1}|=|z^{2}|=|z^{3}|$. To quantify this, we define $\theta(\vec{\mskip
0.5muy}\mskip 1.5mu)$ as the angle between $\vec{\mskip 0.5muy}\mskip 1.5mu$
and the space diagonal in the same octant as $\vec{\mskip 0.5muy}\mskip
1.5mu$. In (Bali:2018nde, , section 4.2), we found that a cut

$\displaystyle\cos\theta(\vec{\mskip 0.5muy}\mskip 1.5mu)$ $\displaystyle\geq
0.9$ (47)

on the data efficiently removes the effect of mirror charges at large $y$,
whilst keeping sufficient statistics.

A different type of anisotropy in the $C_{1}$ data is observed at small $y$
and shown in figure 6. For $A_{\Delta q\Delta q}$, $A_{\delta q\mskip 1.5muq}$
and $B_{\delta q\delta q}$, the data with zero pion momentum exhibit a clear
discrepancy between points $\vec{\mskip 0.5muy}\mskip 1.5mu$ on a coordinate
axis (i.e. with two components being zero) and all other points. This
discrepancy is very strong for $y$ below $5a$ and ceases to be visible above
$7a$. The data for $A_{\delta q\delta q}$ (not shown in the figure) have
larger errors and show only a weak anisotropy for $y<4a$. Only the function
$A_{qq}$ is not affected by this phenomenon, for which we have no explanation.

(a) $A_{qq}$, graph $C_{1}$ (b) $A_{\Delta q\Delta q}$, graph $C_{1}$ (c)
$A_{\delta q\mskip 1.5muq}$, graph $C_{1}$ (d) $B_{\delta q\delta q}$, graph
$C_{1}$

Figure 6:  Twist-two functions at small $y$ for graph $C_{1}$, with scaled
pion momenta $P=0$ and $P=1$ as defined below (45). All points have $py=0$ and
are for $L=40$ and light quarks. The data for $P=\sqrt{2}$, $P=\sqrt{3}$ and
$P=2$ agree with those for $P=1$ within errors but are not shown for the sake
of clarity. Data with $\cos\theta=\sqrt{1/3}$ correspond to $\vec{\mskip
0.5muy}\mskip 1.5mu$ on a coordinate axis.

By contrast, we find that the $C_{1}$ data with nonzero pion momentum and
$py=0$ are isotropic in $\vec{\mskip 0.5muy}\mskip 1.5mu$. For nonzero
momenta, we can hence average all data with the same values of $y$ and $P$,
which greatly decreases statistical errors. We find good agreement between the
$P>0$ data and the $P=0$ data with $\vec{\mskip 0.5muy}\mskip 1.5mu$ on a
coordinate axis for all twist-two functions except $A_{\Delta q\Delta q}$,
where the agreement is only approximate.

(a) $|A_{qq}|$, graph $C_{2}$ (b) $|A_{\Delta q\Delta q}|$, graph $C_{2}$ (c)
$|A_{\delta q\mskip 1.5muq}|$, graph $C_{2}$ (d) $|B_{\delta q\delta q}|$,
graph $C_{2}$

Figure 7:  Twist-two functions at small $y$ for graph $C_{2}$. All points are
for $L=40$, zero pion momentum, and light quarks. For $A_{\delta q\delta q}$
(not shown), one finds a clear anisotropy at $y<4a$, whilst at larger $y$ the
statistical errors are too large for drawing firm conclusions.

We now turn our attention to graph $C_{2}$ at small $y$. Here, we find a very
strong anisotropy in the $P=0$ data. This is shown in figure 7, where we
distinguish points $\vec{\mskip 0.5muy}\mskip 1.5mu$ on the coordinate axes,
which have $\cos\theta=\sqrt{1/3}$, points with
$\sqrt{1/3}<\cos\theta\leq\sqrt{2/3}$, and points with
$\sqrt{2/3}<\cos\theta$. We note that points in a coordinate plane, i.e. with
at least one component of $\vec{\mskip 0.5muy}\mskip 1.5mu$ equal to zero,
have $\sqrt{1/3}\leq\cos\theta\leq\sqrt{2/3}$. In all channels, we see a clear
discrepancy between the points $\vec{\mskip 0.5muy}\mskip 1.5mu$ on a
coordinate axis and all other points. In addition, there is a significant
mismatch between points with $\cos\theta$ above or below $\sqrt{2/3}$ in
several channels, most strongly so in $A_{\delta q\mskip 1.5muq}$. We recall
that a strong anisotropy for $C_{2}$ at small $y$ was also seen for the
correlation functions in our study Bali:2018nde . In section 4.2 of that work,
we argued that this reflects an anisotropy in the lattice propagator between
the two currents, and that points selected by the cut (47) should give the
most reliable results according to the analysis in Cichy:2012is .

We also have $P=1$ data for $C_{2}$, which we can compare with those for
$P=0$. As seen in figure 8, for $A_{qq}$ and $A_{\delta q\mskip 1.5muq}$ the
data at $P=1$ are inconsistent with those at $P=0$, regardless of the value of
$\cos\theta$ in the latter. Since for $P=1$ the condition $py=0$ requires
$\vec{\mskip 0.5muy}\mskip 1.5mu$ to lie in a coordinate plane, we can in fact
not select points satisfying the cut (47) in this case. We therefore discard
our data with nonzero $P$ for $C_{2}$. Testing boost invariance of twist-two
functions at $py=0$ in the presence of the cut (47) would require data with at
least $P=\sqrt{2}$, which we do not have for $C_{2}$.

(a) $A_{qq}$, graph $C_{2}$ (b) $A_{\delta q\mskip 1.5muq}$, graph $C_{2}$

Figure 8:  Twist-two functions at intermediate $y$ for graph $C_{2}$. All
points are for $py=0$, $L=40$ and light quarks.

At $P=0$ and small $y$, we are now in a difficult situation. Points with large
$\cos\theta$ are preferred for $C_{2}$, while for $C_{1}$ the points with the
smallest possible value $\cos\theta=\sqrt{1/3}$ seem to be more reliable,
given that they agree with the $P>0$ data. Applying different cuts in
$\cos\theta$ to the data for $C_{1}$ and $C_{2}$ would prevent us from taking
linear combinations of those graphs at a given $\vec{\mskip 0.5muy}\mskip
1.5mu$. However, we regard combining data point by point in $\vec{\mskip
0.5muy}\mskip 1.5mu$ as highly desirable for a transparent and consistent
treatment of statistical correlations in the jackknife analysis.

To avoid this problem, we choose to discard points with $y<5a$ in our further
analysis, and to apply the cut in (47) to the $P=0$ data for all lattice
graphs. After this cut, data points with equal values of $y$ are averaged also
for $P=0$. We thus avoid the regions where the anisotropy for $C_{1}$ and
$C_{2}$ seen at $P=0$ is most severe. For $C_{1}$, a small discrepancy between
the data with $P=0$ and $P>0$ is still visible up to about $y\sim 8a$, but we
consider this to be at an acceptable level. The result of this procedure is
shown for graph $C_{1}$ in figure 9. The agreement between the data for
different pion momenta is quite good, except for the function $A_{\Delta
q\Delta q}$.

(a) $A_{qq}$, graph $C_{1}$ (b) $A_{\Delta q\Delta q}$, graph $C_{1}$ (c)
$A_{\delta q\mskip 1.5muq}$, graph $C_{1}$ (d) $A_{\delta q\delta q}$, graph
$C_{1}$

Figure 9:  Comparison of twist-two functions for graph $C_{1}$ at different
values of the scaled pion momentum $P$. All points are for $py=0$, $L=40$ and
light quarks. For $B_{\delta q\delta q}$ (not shown), one finds good agreement
between all points, with a similar quality as for $A_{qq}$.

As an exception to the selection just described, we will in section 4.4 use
the $C_{1}$ data for $A_{qq}$ down to $y=3a$, given that in this particular
channel there are no indications of anisotropy or a pion momentum dependence,
as can be seen in figure 6(a).

###  3.4 Excited state contributions

As specified in section 3.2, we have a limited set of data with a separation
of $t=32a$ between pion source and sink. Comparing this with our results for
$t=15a$ allows us to assess the relevance of excited state contributions in
our extraction of the pion matrix elements (8).

On our lattice with size $L=32$, we have $t=32a$ data for graphs $S_{1}$ and
$C_{2}$. Unfortunately, these graphs give a statistical signal consistent with
zero for all twist-two functions and for both source-sink separations. We
hence limit the following discussion to graph $C_{1}$ on our $L=40$ lattice.

In general, we find that the results for the two source-sink separation agree
reasonably well for light quarks, as illustrated in the upper panels of figure
10. For strange quarks, the data have smaller statistical errors and we can
clearly see discrepancies between $t=15a$ and $32a$, as shown in the lower
panels of the figure. Except for the case of $A_{\Delta q\Delta q}$, these
discrepancies are, however, small when compared with the size of the twist-two
functions.

In our data for charm quarks, the statistical signal and the agreement between
the two source-sink separations is excellent for all twist-two functions, and
even better than the one in figure 10(a). With the exception mentioned above,
we thus find no indication for a sizeable contamination from excited states in
our results.

(a) $A_{qq}$, graph $C_{1}$, light quarks (b) $B_{\delta q\delta q}$, graph
$C_{1}$, light quarks (c) $A_{\Delta q\Delta q}$, graph $C_{1}$, strange
quarks (d) $B_{\delta q\delta q}$, graph $C_{1}$, strange quarks

Figure 10:  Comparison of twist-two functions extracted for graph $C_{1}$ with
source-sink separations $t/a=15$ or $32$. All data are for $L=40$, zero pion
momentum, and subject to the cut (47).

###  3.5 Volume dependence

(a) $A_{qq}$, graph $C_{1}$ (b) $A_{\delta q\mskip 1.5muq}$, graph $C_{1}$ (c)
$B_{\delta q\delta q}$, graph $C_{1}$ (d) $A_{\delta q\mskip 1.5muq}$, graph
$C_{2}$

(e) $B_{\delta q\delta q}$, graph $A$ (f) $A_{\delta q\mskip 1.5muq}$, graph
$S_{2}$

Figure 11:  Comparison of data for the two different lattice sizes in our
study. All points are for zero pion momentum, light quarks, and subject to the
cut (47).

Let us finally compare our simulations for light quarks on the lattices with
$L=40$ and $L=32$. In general, the data for the smaller lattice have larger
jackknife errors. This is to be expected from the parameters that determine
the statistical averaging in our simulations. Details for these are given in
table 2 of Bali:2018nde .

For twist-two functions with a small relative error, we typically find a weak
volume dependence compared with the size of the functions themselves, as shown
in panels (a) to (c) of figure 11. In the case of panel (b), this dependence
is, however, statistically significant. For functions that have large relative
errors, the volume dependence appears to be more pronounced in some cases,
especially at low $y$. An example is figure 11(e). One may take this as a
general warning against over-interpreting statistically weak signals in our
simulations.

##  4 Results for zero pion momentum

In this section, we present our results for the twist-two functions at $py=0$.
All data shown in the following are for zero pion momentum and have been
extracted from the lattice with $L=40$ with our standard source-sink
separation $t=15a$. The data selection described at the end of section 3.3
removes regions in which we see strong lattice artefacts in the form of broken
rotational or boost symmetry.

As we explained in section 2.3, twist-two functions at $py=0$ are not directly
related to the Mellin moments of DPDs. Instead, they are Mellin moments of
skewed DPDs, integrated over the skewness parameter $\zeta$. As seen in figure
2, these moments receive contributions from parton configurations that are
different from those in a DPD at $\zeta=0$. When interpreting the results of
the present section, we will _assume_ that these configurations are not
dominant, and that the qualitative features of invariant functions at $py=0$
are the same as for Mellin moments of DPDs at $\zeta=0$. The results presented
in section 5.3 will lend support to this assumption.

Figure 12:  Examples for the partonic regimes of graph $C_{1}$ in a $\pi^{+}$.
The notation $u\bar{u}\bar{d}\mskip 1.5mu|\bar{d}$ is the same as in figure 2,
i.e. the partons to the left of the vertical bar belong to the pion on the
left of the graph.

Figure 13:  As figure 12, but for graph $C_{2}$.

Notice that each of the lattice graphs in figure 4 can contribute to each of
the partonic regimes shown in figure 2. Examples for different regimes of the
connected graphs $C_{1}$ and $C_{2}$ are shown in figures 12 and 13.

###  4.1 Comparison of graphs

In figures 14 and 15, we compare the contributions from different lattice
graphs to the twist-two functions for light quarks. The contributions from
graphs $S_{1}$ and $C_{2}$ are multiplied with a factor $2$ in the figures,
since they always appear with this weight in physical matrix elements
according to (3.1).

For all twist-two functions except $A_{\Delta q\Delta q}$, graph $C_{1}$ gives
a very clear signal, which is positive for $A_{\delta q\delta q}$ and negative
for the other functions. By comparison, the signal for the annihilation graph
$A$ is smaller than the one for $C_{1}$ by an order of magnitude or more,
except for $y>20a$, where the statistical errors prevent us from making a
clear statement. The function $A_{\Delta q\Delta q}$ shows a different
behaviour, with $C_{1}$ and $A$ being of similar size and much smaller than
$C_{1}$ for all other twist-two functions. We recall from section 3.3 that
$A_{\Delta q\Delta q}$ is more strongly affected by lattice artefacts than the
other channels, see figure 9(b). We nevertheless discuss this function here
and in the following, because even with the large systematic uncertainties we
have seen, its qualitative behaviour and overall size compared with the other
polarisation channels are significant results of our simulations.

A clear signal for the connected graph $C_{2}$ is only seen for $A_{qq}$ and
$A_{\delta q\mskip 1.5muq}$, with a sign opposite to the one for graph
$C_{1}$. This signal is most important at small $y$. For the graphs $S_{1}$
and $S_{2}$ with one disconnected fermion loop, the signal we obtain is rather
noisy in all channels. For graph $D$ with two disconnected fermion loops, the
signal after vacuum subtraction is even more noisy and not shown.

(a) $A_{qq}$, light quarks

(b) $A_{\Delta q\Delta q}$, light quarks

(c) $A_{\delta q\mskip 1.5muq}$, light quarks

Figure 14:  Contributions of the different lattice graphs to twist-two
functions for light quarks.

(a) $A_{\delta q\delta q}$, light quarks

(b) $B_{\delta q\delta q}$, light quarks

Figure 15:  Continuation of figure 14.

From our simulations with the strange quark mass, we have only data for graphs
$C_{1}$ and $A$. In all channels, we obtain an excellent signal for $C_{1}$,
whereas for $A$ the statistical significance is typically not much larger than
one standard deviation. In the region $5a\leq y\leq 15a$, we find that $A$ is
smaller than $C_{1}$ by one to two orders of magnitude, except for $A_{\Delta
q\Delta q}$. For $A_{qq}$ and $A_{\delta q\delta q}$, we see in figure 16 that
the behaviour of $A$ is quite flat, unlike the one of $C_{1}$, so that at
large $y$ the two graphs become more comparable in size. As in the case of
light quarks, the function $A_{\Delta q\Delta q}$ behaves differently, with
graph $A$ being smaller than $C_{1}$ at $y\sim 5a$ and the data for both
graphs having a zero crossing a bit below $y=9a$. Recall, however, that also
for strange quarks we see stronger lattice artefacts in $A_{\Delta q\Delta q}$
than in other channels, as seen in figure 10(c).

(a) $|A_{qq}|$, graphs $C_{1}$ and $A$ (b) $|A_{\delta q\delta q}|$, graphs
$C_{1}$ and $A$

Figure 16:  Comparison of graphs $A$ and $C_{1}$ for strange quarks. Here and
in subsequent figures with a logarithmic scale, we stop showing data for
individual quark masses at a value of $y$ beyond which the error bars become
so large that they would obscure the plot.

From our simulations with the charm quark mass, we have data for all graphs
except $S_{2}$. A clear nonzero signal is seen for $C_{1}$ and $C_{2}$ up to
$y\sim 10a$ to $15a$, with $2C_{2}$ being smaller than $C_{1}$ by at least one
order of magnitude. The signal for $A$ and $S_{1}$ is in general consistent
with zero. The only exception to this is $A_{\delta q\mskip 1.5muq}$. For this
function, we see a clear signal for $2S_{1}$ at $y$ around $5a$, which is
about 50 times smaller than the one for $C_{1}$. We also see a weak $1\sigma$
signal for $A$, which we do not wish to over-interpret.

By and large, we find that for all quark masses the only graphs that give
signals of appreciable size are $C_{1}$ and, in several cases, $C_{2}$. We
therefore take a closer look at these graphs in the next subsection. The
annihilation graph is negligible, except in the case of $A_{\Delta q\Delta q}$
for light or strange quarks, where the signal from graph $C_{1}$ is small by
itself. Disconnected graphs either have a negligibly small signal or large
statistical errors.

###  4.2 Results for connected graphs

The contribution of graph $C_{1}$ to the twist-two function $A_{qq}$ for
unpolarised partons is negative for all three quark masses in our study. We
recall from (2.2) that the regime with a quark and an antiquark in the pion
contributes with a negative sign to the lowest Mellin moment of a DPD. The
same holds for the Mellin moment of a skewed DPD, and hence for $A_{qq}$ at
$py=0$. A negative sign of $A_{qq}$ is easily understood by the dominance of
the valence $q\bar{q}$ Fock state, which is probed by graph $C_{1}$ as shown
in the first panel of figure 12.

The situation is different for graph $C_{2}$, whose partonic representation
always involves a higher Fock state of the pion. The $Z$-graphs in figure 13
probe the $qq$, $\bar{q}\bar{q}$ and $q\bar{q}$ regimes in a similar manner.
We find that for all quark masses, the contribution of $C_{2}$ to $A_{qq}$ is
positive, which means that for a given distance $y$ this graph gives a larger
probability for finding a $qq$ or $\bar{q}\bar{q}$ rather than a $q\bar{q}$
pair.

(a) $-A_{qq}$, graph $C_{1}$ (b) $-mA_{\delta q\mskip 1.5muq}$, graph $C_{1}$
(c) $A_{\delta q\delta q}$, graph $C_{1}$ (d) $-m^{2}B_{\delta q\delta q}$,
graph $C_{1}$

Figure 17:  Mass dependence of twist-two functions for graph $C_{1}$.

(a) $A_{\Delta q\Delta q}$, graph $C_{1}$

Figure 18:  As figure 18, for $A_{\Delta q\Delta q}$ and with a linear instead
of a logarithmic scale.

Let us now take a closer look at the mass dependence of our results for graph
$C_{1}$. We multiply $A_{\delta q\mskip 1.5muq}$ and $B_{\delta q\delta q}$
with the power of the meson mass $m$ with which they appear in the
decomposition (2.2) of two-current matrix elements. We see in figure 18 that
for all twist-two functions except $A_{\Delta q\Delta q}$, the decrease with
$y$ becomes stronger with increasing quark mass, which simply reflects the
decreasing size of the meson. At $y\sim 5a$, the functions $A_{qq}$,
$mA_{\delta q\mskip 1.5muq}$ and $m^{2}B_{\delta q\delta q}$ are of comparable
size for all quark masses, whereas $A_{\delta q\delta q}$ increases with the
mass. The behaviour of $A_{\Delta q\Delta q}$ for light and strange quarks is
qualitatively different from the one of the other functions, as is evident
from figure 18. For charm quarks, $A_{\Delta q\Delta q}$ is approximately
exponential in $y$, with a logarithmic slope similar to the one of $A_{qq}$. A
fit of the $y$ dependence of the twist-two functions for light quarks is
presented in section 5.2.

We now discuss graph $C_{2}$, for which we have data with light quarks and
with charm. For the functions $A_{\Delta q\Delta q}$, $A_{\delta q\delta q}$
and $B_{\delta q\delta q}$, the light quark data is too noisy for a meaningful
comparison with charm results, so that we focus on $A_{qq}$ and $A_{\delta
q\mskip 1.5muq}$. As is seen in figure 19, the size of both functions is
significantly smaller for charm quarks. This is plausible: as discussed in the
previous subsection, the partonic interpretation of graph $C_{2}$ always
involves a Fock state with at least two quarks and two antiquarks in the
meson, whereas for $C_{1}$ we have the regime shown in the first panel in
figure 12, which involves only the quark-antiquark Fock state. The $y$
dependence of $A_{qq}$ and $A_{\delta q\mskip 1.5muq}$ is also qualitatively
different for the two masses: for charm we observe a clear and steep
exponential falloff, whereas for light quarks, the logarithmic slope of both
functions decreases around $y\sim 0.5\operatorname{fm}$.

(a) $A_{qq}$, graph $C_{2}$ (b) $mA_{\delta q\mskip 1.5muq}$, graph $C_{2}$

Figure 19:  Mass dependence of twist-two functions for graph $C_{2}$.

###  4.3 Polarisation effects

A major aim of our study is to investigate the strength and pattern of spin
correlations between two partons in a pion. We spelled out the physical
interpretation of polarised DPDs in section 2.1. This interpretation extends
to the corresponding twist-two functions at $py=0$, provided that these are
dominated by partonic regimes associated with DPDs at $\zeta=0$. Under this
assumption, comparing $A_{\Delta q\Delta q}$ and $A_{\delta q\delta q}$ with
$A_{qq}$ indicates whether two partons prefer to have their spins aligned or
anti-aligned, with $A_{\Delta q\Delta q}$ referring to longitudinal and
$A_{\delta q\delta q}$ to transverse polarisation. We will refer to these as
“spin-spin correlations”. Note that, according to (2.2), a $q\bar{q}$ pair
with aligned spins contributes with a negative sign to $A_{qq}$ and $A_{\delta
q\delta q}$ and with a positive sign to $A_{\Delta q\Delta q}$, whereas a $qq$
pair with aligned spins contributes with a positive sign to all three
functions. The comparison of $myA_{\delta q\mskip 1.5muq}$ and
$m^{2}|y^{2}|B_{\delta q\delta q}$ with $A_{qq}$ tells us about the strength
of correlations between the transverse spin of one or both observed partons
and the distance $\boldsymbol{y}$ between these partons in the transverse
plane. We refer to this as “spin-orbit correlations” in the following. The
pre-factors $my$ and $m^{2}|y^{2}|$ in $myA_{\delta q\mskip 1.5muq}$ and
$m^{2}|y^{2}|B_{\delta q\delta q}$ follow from the decompositions (2.1) and
(2.2).

We note that the probability interpretation of polarised DPDs implies
positivity constraints Diehl:2013mla  that extend the well-known Soffer bound
for single parton distributions Soffer:1994ww . These bounds imply that
$|f_{\Delta q\Delta\bar{q}}|$, $|f_{\delta q\delta\bar{q}}|$, $|myf_{\delta
q\mskip 1.5mu\bar{q}}|$ and $|m^{2}y^{2}f^{t}_{\delta q\delta\bar{q}}|$ are
bounded by $f_{q\bar{q}}$. Corresponding bounds do not hold for the lowest
Mellin moments of DPDs because of the relative minus sign between quark and
antiquark contributions in (2.2). They hold even less for the moments of
skewed DPDs, which do not represent probabilities to start with. Nevertheless,
in a loose sense, the size of $A_{qq}$ sets a natural scale for the other
twist-two functions (multiplied with $my$ or $m^{2}|y^{2}|$ as appropriate).

(a) graph $C_{1}$, light quarks

(b) graph $C_{1}$, strange quarks

(c) graph $C_{1}$, charm quarks

Figure 20:  Effects of transverse (left) and longitudinal (right) polarisation
for graph $C_{1}$.

In the following, we consider polarisation effects separately for the
connected graphs $C_{1}$ and $C_{2}$. Their physical interpretation is rather
different, as becomes clear from figures 12 and 13 and our discussion in the
previous subsection. For graph $C_{1}$, polarisation effects reflect rather
directly correlations between the quark and antiquark in the pion valence
state, whereas for graph $C_{2}$ they are deeply connected with sea quark
degrees of freedom.

Starting our discussion with graph $C_{1}$, we see in the top panels of figure
20 that by far the strongest polarisation effect seen for light quarks is the
spin-orbit correlation for a single parton, followed by the spin-orbit
correlation involving both partons. Both the transverse and the longitudinal
spin-spin correlations are very small. This is completely different from the
simple picture of a pion as a $q\bar{q}$ pair in an $S$-wave, for which one
would obtain 100% anti-alignment of both transverse and longitudinal spins.

All spin correlations increase considerably with the quark mass. For charm
quarks, $myA_{\delta q\mskip 1.5muq}$ is almost as large as $A_{qq}$. Spin-
spin correlations are also important for charm: the spins of the quark and
antiquark are anti-aligned by about 75% for transverse and by about 50% for
longitudinal polarisation. We note that this is still quite far away from the
non-relativistic limit, in which transverse and longitudinal spin correlations
become equal.

We note that the pion mass for our simulations with light quarks,
$m_{\pi}\approx 295\operatorname{MeV}$, is quite a bit larger than the
physical value. A naive extrapolation of the polarisation patterns just
described suggests that at the physical point the spin-orbit correlation for
one polarised parton may be substantial, whilst correlations involving two
quark spins might be even smaller than the ones we see for light quarks in the
present study.

(a) graph $C_{2}$, light quarks

(b) graph $C_{2}$, charm quarks

Figure 21:  As figure 20, but for graph $C_{2}$. We have no strange quark
results for this case. The light quark data for $A_{\delta q\delta q}$ and
$m^{2}|y^{2}|B_{\delta q\delta q}$ is very noisy and not shown for the sake of
clarity.

We now turn to our results for graph $C_{2}$, which are shown in figure 21.
For light quarks, we see a substantial spin-orbit correlation of order 50% for
a single parton. The spin-spin correlation for longitudinal polarisation is
also of order 50% for $y\sim 0.35\operatorname{fm}$, but quickly decreases and
is negligible already around $y\sim 0.5\operatorname{fm}$. For all other spin
dependent correlations, the data for light quarks are too noisy to extract any
physics.

With charm quarks, we have an excellent statistical signal for all twist-two
functions. We find that all spin correlations for graph $C_{2}$ are
appreciable, apart from the one described by $m^{2}|y^{2}|B_{\delta q\delta
q}$. Notice that $A_{\Delta q\Delta q}$ has the same sign for $C_{1}$ and
$C_{2}$, unlike all other twist-two functions. If (as suggested by the sign of
$A_{qq}$) the dominant parton configuration probed by the twist-two operators
is a $c\bar{c}$ pair for graph $C_{1}$ and a $cc$ pair for graph $C_{2}$, then
the longitudinal parton spins tend to be anti-aligned in both cases.

###  4.4 Test of the factorisation hypothesis

We now test the factorisation hypothesis for $A_{ud}(y^{2},py=0)$ that we
derived in section 2.4. We restrict ourselves to the contribution from the
connected graph $C_{1}$. Taking the full combination of graphs in the first
line of (3.1) is not an option because of the huge errors in our results for
the doubly disconnected graph $D$. By contrast, we see in figure 14(a) that
$S_{1}$ is consistent with zero for $A_{qq}$ (although within errors much
larger than those on $C_{1}$). We find it plausible to expect that the
contribution from $D$ is even smaller than the one of $S_{1}$, since $D$ has
two disconnected fermion loops with one operator insertion.

The factorisation hypothesis (41) involves the vector form factor of the pion.
We have extracted this form factor from our lattice simulations, using the
full number of 2025 gauge configurations available for our lattice with
$L=40$. As we consider only the connected contribution to the two-current
correlation function, we restrict ourselves to the connected graph for the
form factor as well. We fit the form factor data to a power law

$\displaystyle F_{u,V}(t)$
$\displaystyle=-F_{d,V}(t)=\bigl{(}1-t/M^{2}\bigr{)}^{-p}\,.$ (48)

We use two fit variants, which gives us a handle on the bias of the
extrapolation to $-t>1.15\operatorname{GeV}^{2}$, where we have no data. Such
an extrapolation bias is inevitable when we Fourier transform from momentum to
position space, as is required in (41). In a monopole fit, we fix $p=1$ and
obtain $M=777(12)\operatorname{MeV}$. Leaving the power free, we obtain
$p=1.173(69)$ and $M=872(16)\operatorname{MeV}$. Both fits give a very good
description of our lattice data, as shown in figure 14a of Bali:2018nde .

With the ansatz (48), the two-dimensional Fourier transform on the r.h.s. of
(41) can be carried out analytically. We compute the remaining integral over
$\zeta$ numerically. The results obtained with the two form factor fits agree
very well for $y>0.2\operatorname{fm}$. In panel (a) of figure 22 we compare
the two sides of the factorisation hypothesis (41), and in panel (b) we show
the ratio of the r.h.s. to the l.h.s. of the equation. We see a clear
deviation from the factorised ansatz, which does however not exceed 30% in the
considered $y$ range. One may thus say that the factorised ansatz provides a
rough approximation of the two-current correlator.

(a) (b)

Figure 22:  Test of the factorisation hypothesis (41) for the invariant
function $A_{ud}$. (a): data for $A_{ud}$ (restricted to the contribution from
graph $C_{1}$) compared with the integral over form factors on the r.h.s. of
(41). The form factors are determined by a monopole or a $p$-pole fit. (b):
ratio of the form factor integral on the r.h.s. of (41) to the data for
$A_{ud}$.

###  4.5 Physical matrix elements

(a) $A_{ud}\mskip 1.5mu|_{\pi^{+}}$ (b) $myA_{\delta u\mskip 1.5mud}\mskip
1.5mu|_{\pi^{+}}$ (c) $A_{uu}\mskip 1.5mu|_{\pi^{+}}$ (d) $myA_{\delta u\mskip
1.5muu}\mskip 1.5mu|_{\pi^{+}}$

(e) $A_{uu}\mskip 1.5mu|_{\pi^{0}}$ (f) $myA_{\delta u\mskip 1.5muu}\mskip
1.5mu|_{\pi^{0}}$

Figure 23:  Twist-two functions at $py=0$ for the flavour combinations $ud$ or
$uu$ in a $\pi^{+}$ or a $\pi^{0}$. Lattice graphs are combined according to
(3.1), except for of graph $D$, which is affected by huge errors and hence
omitted. All results are for light quarks.

(a) $A_{\Delta u\Delta u}\mskip 1.5mu|_{\pi^{+}}$ (b) $A_{\Delta u\Delta
u}\mskip 1.5mu|_{\pi^{0}}$

Figure 24:  Continuation of figure 23.

We now investigate the combinations (3.1) of lattice graphs that appear in the
matrix elements of currents between charged or neutral pions. We omit the
doubly disconnected graph $D$ throughout, because its statistical errors are
much larger than the signal for any other graph. Since data for the full set
of remaining graphs is only available for light quarks, we restrict our
attention to this case. The results are shown in figures 23 and 24 for the
flavour combinations $ud$ and $uu$. The combinations $dd$ and $du$ can be
obtained from the symmetry relations (16).

As can be expected from figures 14 and 15, the statistical errors of the
physical combinations are significantly larger than those for the connected
graphs alone. Nevertheless, we see a clear negative signal for $A_{ud}$ in a
$\pi^{+}$. As discussed in section 4.2, this can be understood as a dominance
of the valence Fock state $u\bar{d}$ over Fock states that contain $ud$,
$\bar{u}\bar{d}$ or $\bar{u}d$. The function $A_{uu}$ in a $\pi^{+}$ has a
clear positive signal at small distances $y$. This reflects the behaviour of
graph $C_{2}$ and corresponds to a larger probability for finding two $u$
quarks rather than a $u\bar{u}$ pair at small separation $y$. Remarkably, the
signal at small $y$ is of comparable size for $A_{uu}$ and $A_{ud}$, which
implies that Fock states containing sea quarks do play an important role in
this region. As for polarisation effects, a clear signal for $ud$ or $uu$ in a
$\pi^{+}$ is only seen for $myA_{\delta q\mskip 1.5muq}$ and shown in the
right panels of figure 23. Comparing this with $A_{qq}$, we see that spin-
orbit correlations are appreciable for both flavour combinations.

The flavour combination $uu$ in a $\pi^{0}$ involves the sum $C_{1}+2C_{2}$.
We observe a very strong compensation between the two connected graphs, which
results in a marginal signal for $A_{uu}$ and $myA_{\delta u\mskip 1.5muu}$.
The twist-two functions for $ud$ in a $\pi^{0}$ receive no contribution from
connected graphs at all. Within errors, the corresponding results are zero for
all combinations of currents, and we do not show them here.

Among all polarised twist-two functions other than $myA_{\delta q\mskip
1.5muq}$, a marginally nonzero signal is only seen for the longitudinal spin
correlation $A_{\Delta u\Delta u}$ in a $\pi^{+}$ or a $\pi^{0}$. This is
dominated by the contribution from $C_{2}$ in both cases and shown in figure
24.

We recall that the assumption (7) going into the “pocket formula” for double
parton scattering implies that the unpolarised DPDs in a given hadron have the
same $y$ dependence for all parton combinations. If this were to hold also for
skewed DPDs, the functions $A_{ud}$ and $A_{uu}$ in a $\pi^{+}$ should have
the same $y$ dependence as well. A comparison between these two functions is
shown in figure 25. Although there is a hint for a different behaviour,
especially at low $y$, the large errors in $A_{uu}$ prevent us from making a
definitive statement.

(a) $|A_{qq}|_{\pi^{+}}$ (linear scale) (b) $|A_{qq}|_{\pi^{+}}$ (logarithmic
scale)

Figure 25:  Comparison of the twist-two functions $-A_{ud}$ and $A_{uu}$ in a
$\pi^{+}$. The data are the same as shown in figures 23(a) and 23(c). Points
with very large errors have been omitted for the sake of clarity.

##  5 Results for nonzero pion momentum

In this section, we use our data for nonzero pion momentum to study the $py$
dependence of the twist-two functions. We restrict our study to graph $C_{1}$
for light quarks on the $L=40$ lattice: only in this case do we have
simulations for a sufficient number of pion momenta. Since graph $C_{1}$
dominates the twist-two matrix elements for $ud$ in a $\pi^{+}$, we will write
$A_{ud}$, $A_{\delta u\mskip 1.5mud}$, …for twist-two functions and $I_{ud}$,
$I_{\delta u\mskip 1.5mud}$, …for Mellin moments in what follows.

###  5.1 Fit ansatz for the $py$ dependence

We start by proposing a functional ansatz for the twist-two functions, which
is based on their relation (33) with the Mellin moments of skewed DPDs. We use
this ansatz to fit the $py$ dependence of our lattice data. This will allow
for a model dependent extension of the twist-two functions to all values of
$py$, beyond the region (26) available on a Euclidean lattice. This will in
turn allow for a model dependent extraction of the Mellin moments of DPDs at
zero skewness. For ease of notation, we write $A(y^{2},py)$ to denote any of
the twist-two functions $A_{ud}$, …, $A_{\delta u\delta d}$, $B_{\delta
u\delta d}$. Likewise, we write $I(y^{2},py)$ for the Mellin moments $I_{ud}$,
…, $I_{\delta u\delta d}$, $I^{t}{}_{\\!\delta u\delta d}$.

The basis of our ansatz is the assumption that, in its support region
$-1\leq\zeta\leq 1$, the skewed moment $I(y^{2},\zeta)$ can be approximated by
a polynomial in $\zeta$,

$\displaystyle I(y^{2},\zeta)$
$\displaystyle=\pi\sum_{n=0}^{N}a_{n}(y^{2})\,\zeta^{2n}$ (49)

with some integer $N$, where we used the symmetry relation (30) to restrict
terms to even powers of $\zeta$. We write $=$ instead of $\approx$ in the
spirit of a fit ansatz, i.e. we do not claim that (49) is exact. We currently
have no guidance from theory regarding the $\zeta$ dependence of
$I(y^{2},\zeta)$, and (49) should be taken as a simple ansatz that is to be
validated by data. As there are no known constraints on the behaviour of
skewed DPDs at the edge $\zeta=1$ of their support region, we allow
$I(y^{2},\zeta)$ to be finite at that point.

Inverting the Fourier transform in (32), we obtain

$\displaystyle A(y^{2},py)$
$\displaystyle=\sum_{n=0}^{N}a_{n}(y^{2})\,h_{n}(py)\,,$ (50)

where we introduced the functions

$\displaystyle h_{n}(x)$
$\displaystyle=\frac{1}{2}\int_{-1}^{1}d\zeta\,e^{i\zeta x}\,\zeta^{2n}\,.$
(51)

A crucial property of the ansatz (50) is that its Fourier transformation (32)
has the correct support in $\zeta$.

Let us collect a few properties of the functions $h_{n}(x)$. From their
definition, one easily derives

$\displaystyle h_{n}(0)$ $\displaystyle=\frac{1}{1+2n}\,,$
$\displaystyle\frac{d^{2}h_{n}(x)}{dx^{2}}$ $\displaystyle=-h_{n+1}(x)$ (52)

and thus obtains the Taylor series

$\displaystyle h_{n}(x)$
$\displaystyle=\sum_{m=0}^{\infty}\frac{(-1)^{m}}{1+2n+2m}\;\frac{x^{2m}}{(2m)!}\;.$
(53)

An explicit representation is given by

$\displaystyle h_{n}(x)=s_{n}(x)\sin x+c_{n}(x)\cos x$ (54)

with rational functions

$\displaystyle s_{n}(x)$
$\displaystyle=\sum_{m=0}^{n}\frac{(2n)!}{(2n-2m)!}\,\frac{(-1)^{m}}{x^{1+2m}}\,,$
$\displaystyle c_{n}(x)$
$\displaystyle=\sum_{m=0}^{n-1}\frac{(2n)!}{(2n-2m-1)!}\,\frac{(-1)^{m}}{x^{2+2m}}\,.$
(55)

For $n=0$ and $n=1$, these functions read

$\displaystyle s_{0}(x)$ $\displaystyle=1/x\,,$ $\displaystyle c_{0}(x)$
$\displaystyle=0\,,$ $\displaystyle s_{1}(x)$
$\displaystyle=(x^{2}-2)/x^{3}\,,$ $\displaystyle c_{1}(x)$
$\displaystyle=2/x^{2}\,.$ (56)

In terms of the normalised quantities

$\displaystyle\widehat{A}(y^{2},py)$
$\displaystyle=\frac{A(y^{2},py)}{A(y^{2},py=0)}\,,$
$\displaystyle\hat{a}_{n}(y^{2})$
$\displaystyle=\frac{a_{n}(y^{2})}{A(y^{2},py=0)}$ (57)

our ansatz (50) reads

$\displaystyle\widehat{A}(y^{2},py)$
$\displaystyle=\sum_{n=0}^{N}\hat{a}_{n}(y^{2})\,h_{n}(py)\,.$ (58)

Using (34) and (53), we then obtain

$\displaystyle\langle\zeta^{2m}\rangle(y^{2})$ $\displaystyle=\Biggl{[}\mskip
1.5mu(-1)^{m}\,\frac{\partial^{2m}\widehat{A}(y^{2},py)}{(\partial\mskip
1.5mupy)^{2m}}\mskip
1.5mu\Biggr{]}_{py=0}=\,\sum_{n=0}^{N}\frac{1}{1+2n+2m}\;\hat{a}_{n}(y^{2})\,.$
(59)

Let us now describe our general fitting procedure. In order to achieve stable
fits, we first determine the $y^{2}$ dependence of $A(y^{2},py=0)$. This
includes the information from data with zero pion momentum and has typically
much smaller errors than the data for nonzero $py$.

In a second step, we fit the $y$ dependent coefficients $\hat{a}_{n}(y^{2})$
in the ansatz (58). To make the degrees of freedom of this fit explicit, we
consider the moments $\langle\zeta^{2m}\rangle(y^{2})$ for $m=0,\ldots,N$.
Inverting the relation (59), we obtain

$\displaystyle\hat{a}_{n}(y^{2})$
$\displaystyle=\sum_{m=0}^{N}(T^{-1})_{nm}\;\langle\zeta^{2m}\rangle(y^{2})\,,$
(60)

where $T$ is the $(N+1)\times(N+1)$ matrix with elements

$\displaystyle T_{mn}=(1+2n+2m)^{-1}\,.$ (61)

Since by definition $\langle\zeta^{0}\rangle(y^{2})=1$, we can thus fit the
$py$ dependence of the twist-two functions to (58) and (60) with $N$ fit
parameters $\langle\zeta^{2}\rangle$, …, $\langle\zeta^{2N}\rangle$ at each
value of $y^{2}$. We call this “local fits” in the following, where “local”
means “local in $y^{2}$”.

To obtain a parametrisation of both the $py$ and the $y^{2}$ dependence, we
assume an expansion

$\displaystyle\langle\zeta^{2m}\rangle(y^{2})$
$\displaystyle=\sum_{k=0}^{K}c_{mk}\,\sqrt{-y^{2}}^{\,k}\,.$ (62)

This is referred to as our “global fit”. By virtue of (58) and (60), this
corresponds to an expansion of $\widehat{A}(y^{2},py)$ in powers of
$\sqrt{-y^{2}}$. The condition $\langle\zeta^{0}\rangle(y^{2})=1$ implies
$c_{0\mskip 1.5muk}=\delta_{0\mskip 1.5muk}$.

###  5.2 Fitting the data

We recall that we have data for $p=0,1,\sqrt{2},\sqrt{3}$ and $2$ in units of
$2\pi/(La)\approx 437\operatorname{MeV}$. For a given value of $y$, this
allows for a maximum value $4\pi y/(La)\approx 6.28\,y/(20a)$ for $|py|$. We
apply the cut (47) on the angle $\theta$ to the $p=0$ data, but not to the
data with $p>0$. We then average all data points with the same values of $py$
and $y^{2}$.

We find that the twist-two functions at $py=0$ can be well described by a
superposition of two exponentials,

$\displaystyle A(y^{2},py=0)$ $\displaystyle=A_{1}\,e^{-a_{1}\mskip
1.5mu(y-y_{\text{min}})}+A_{2}\,e^{-a_{2}\mskip 1.5mu(y-y_{\text{min}})}$
$\displaystyle\text{for }y_{\text{min}}\leq y\leq y_{\text{max}}\,,$ (63)

with $y_{\text{min}}=5a=0.355\operatorname{fm}$ and
$y_{\text{max}}=20a=1.42\operatorname{fm}$. We do not include data with
$y>y_{\text{max}}$, because they have large errors and are increasingly
affected by finite size effects. The resulting fit parameters are given in
table 2. Let us emphasise that these fits are _not_ suitable for extrapolating
the twist-two functions to values significantly below $y=y_{\text{min}}$.

function | $A_{1}\,[\operatorname{fm}^{-2}]$ | $a_{1}\,[\operatorname{fm}^{-1}]$ | $A_{2}\,[\operatorname{fm}^{-2}]$ | $a_{2}\,[\operatorname{fm}^{-1}]$ | $\chi^{2}/\mathrm{dof}$
---|---|---|---|---|---
$A_{ud}$ | $-0.1163(39)$ | $2.150(68)$ | $0.0141(34)$ | $11.5\pm 2.2$ | $0.95$
$A_{\Delta u\Delta d}$ | $-0.0414(77)$ | $6.71(36)$ | $0.0326(74)$ | $4.21(52)$ | $0.94$
$A_{\delta u\mskip 1.5mud}$ | $-0.1157(62)$ | $3.786(89)$ | $0.0222(68)$ | $6.38(30)$ | $1.76$
$A_{\delta u\delta d}$ | $0.0133(24)$ | $2.11(28)$ | $-0.0018(23)$ | $7.8\pm 7.4$ | $0.34$
$B_{\delta u\delta d}$ | $-0.0491(59)$ | $4.50(20)$ | $-0.0084(64)$ | $9.4\pm 2.0$ | $0.95$
Table 2:  Parameters for the fit (63) of twist-two functions at $py=0$ in the
region $5a\leq y\leq 20a$. Throughout this section, we consider the data for
graph $C_{1}$ and light quarks on our lattice with $L=40$.

We notice a relatively high value of $\chi^{2}/\mathrm{dof}$ in the fit for
$A_{\delta u\mskip 1.5mud}$. This is due to some scatter in the data at high
$y$, which comes from points with large $p$. Repeating the fit with an upper
limit $y\leq 15a$, we find that $\chi^{2}/\mathrm{dof}$ decreases from $1.76$
to $0.9$ for $A_{\delta u\mskip 1.5mud}$. By comparison, the value of
$\chi^{2}/\mathrm{dof}$ in the fit for $A_{ud}$ decreases from $0.95$ to $0.6$
with the same reduction of the fitting range.

We then proceed and fit the $py$ dependence to (58) and (60) locally in
$y^{2}$. To have enough data in these fits, we introduce bins in $y$ and
combine all points with $(n-1/2)\mskip 1.5mua<y<(n+1/2)\mskip 1.5mua$ for
integer $n$ between $5$ and $20$. In addition, we fit the combined $y^{2}$ and
$py$ dependence of $\widehat{A}$ to (58), (60) and (62). We explored fits with
different maximum values $N$ and $K$ in the sums and find that, given the fit
range and the statistical quality of our data, an adequate choice is $N=1$ for
local fits and $N=1$, $K=1$ for the global fit. The parameters of the global
fit are given in table 3. If we take $N=2$ instead, the error bands of the fit
results for $\widehat{A}$ increase significantly, whilst the decrease of
$\chi^{2}/\text{dof}$ is minor. We hence conclude that we would over-fit the
data by choosing $N=2$ or even higher values.

function | $c_{10}$ | $c_{11}\,[\operatorname{fm}^{-1}]$ | $\chi^{2}/\mathrm{dof}$
---|---|---|---
$\widehat{A}_{ud}$ | $0.096(40)$ | $0.247(39)$ | $1.19$
$\widehat{A}_{\Delta u\Delta d}$ | $-0.43(73)$ | $0.17(75)$ | $0.68$
$\widehat{A}_{\delta u\mskip 1.5mud}$ | $0.102(50)$ | $0.111(47)$ | $1.37$
$\widehat{A}_{\delta u\delta d}$ | $-0.05(14)$ | $0.31(12)$ | $0.80$
$\widehat{B}_{\delta u\delta d}$ | $-0.023(90)$ | $0.242(92)$ | $0.99$
Table 3:  Parameters of the fit of the combined $y^{2}$ and $py$ dependence of
the normalised twist-two functions $\widehat{A}(y^{2},py)$ to (58), (60) and
(62) with $N=K=1$.

We compare our data and fits in figure 26 for different functions at $y=15a$
and in figure 27 for $\widehat{A}_{ud}$ at $y=5a$ and $10a$. We find good
agreement between the local and global fits. Note that the twist-two functions
are symmetric in $py$ due to $PT$ invariance, which is realised on the
lattice. A departure from this symmetry in the data must therefore be due to
statistical fluctuations. Many data points have admittedly large errors, which
is a consequence of at least one of $y$ or $p$ being large. Nevertheless, the
fitted parameters for all functions except $\widehat{A}_{\Delta u\Delta d}$
are in general well determined, and the corresponding error bands of the fit
results are reasonably small. As is seen in figure 26(b), the data for
$\widehat{A}_{\Delta u\Delta d}$ are much too noisy for fitting the $py$
dependence, and we exclude this function from our further discussion.

(a) $\widehat{A}_{ud}$ at $y=15a$ (b) $\widehat{A}_{\Delta u\Delta d}$ at
$y=15a$ (c) $\widehat{A}_{\delta u\mskip 1.5mud}$ at $y=15a$ (d)
$\widehat{B}_{\delta u\delta d}$ at $y=15a$

Figure 26:  Data and fits of the $py$ dependence of normalised invariant
functions. Dark points show data at $y=15a$. Light points show data in a $y$
range of $a/2$ around $15a$, which are included in the local fits. The plot
for $\widehat{A}_{\delta u\delta d}$ (not shown) is qualitatively similar to
the one for $\widehat{B}_{\delta u\delta d}$.

(a) $\widehat{A}_{ud}$ at $y=5a$ (b) $\widehat{A}_{ud}$ at $y=10a$

Figure 27:  Data and fits of the $py$ dependence of $\widehat{A}_{ud}$ for
different $y$. The meaning of dark and light points is as in figure 26.

In the data for $y=15a$, we see an indication for zero crossings around
$|py|=4$ in several twist-two functions. That this can be reproduced with a
superposition of the two functions $h_{0}(py)$ and $h_{1}(py)$ gives us some
confidence in our fit ansatz.

Using our fits, we can compute the moment $\langle\zeta^{2}\rangle(y^{2})$
associated with $I(y^{2},\zeta)$, which according to (59) follows from the
curvature of $\widehat{A}(y^{2},py)$ at $py=0$. The results are shown in
figure 28. We find again good agreement between the local and global fits. A
clear $y$ dependence of $\langle\zeta^{2}\rangle$ is observed, except for
$I_{\delta u\mskip 1.5mud}$. The values of $\langle\zeta^{2}\rangle$ are not
too large, especially for small $y$. Their size does, however, imply that
nonzero values of the skewness $\zeta$ must play some role in the integral
representation $\pi
A(y^{2},py=0)=\int_{0}^{1}\mskip-1.5mud\zeta\,I(y^{2},\zeta)$.

(a) $I_{ud}$ (b) $I_{\delta u\mskip 1.5mud}$ (c) $I_{\delta u\delta d}$ (d)
$I^{t}_{\delta u\delta d}$

Figure 28:  Values of the moment $\langle\zeta^{2}\rangle(y^{2})$ associated
with $I(y^{2},\zeta)$, extracted by local fits (data points) and the global
fit (bands).

###  5.3 Mellin moments of DPDs

We now use the global fit described in the last section to reconstruct the
lowest Mellin moments of skewed DPDs. Let us re-emphasise that such a
reconstruction is necessarily dependent on the functional ansatz we have made,
given the impossibility to constrain the full $py$ dependence of twist-two
functions with lattice simulations. We recall that the results for the spin
correlation $\Delta u\Delta d$ are too noisy and hence omitted in the
following.

We can easily derive the analytic form of the Mellin moments for our fits by
inverting the $2\times 2$ matrix $T_{mn}$ in (61). This gives

$\displaystyle I(y^{2},\zeta)$ $\displaystyle=\frac{3\pi}{4}\mskip
1.5mu\biggl{\\{}\mskip 1.5mu3-5\mskip
1.5mu\langle\zeta^{2}\rangle(y^{2})-5\mskip 1.5mu\zeta^{2}\mskip
1.5mu\Bigl{[}1-3\mskip 1.5mu\langle\zeta^{2}\rangle(y^{2})\mskip
1.5mu\Bigr{]}\biggr{\\}}\,A(y^{2},py=0)\,.$ (64)

The values of $\langle\zeta^{2}\rangle(y^{2})$ for $y\leq 20a$ are in the
range between $0$ and $0.5$ for all twist-two functions. The combination
$3-5\langle\zeta^{2}\rangle$ in (64) is therefore always positive and varies
between $3$ and $0.5$. We can hence anticipate that the dependence of the
Mellin moments $I(y^{2},\zeta=0)$ on $y$ and on the polarisation indices
should roughly follow the corresponding dependence of $A(y^{2},py=0)$. By
contrast, the coefficient of $\zeta^{2}$ in (64) has a larger variation and
can change sign as a function of $y$. Our results for the $y$ and $\zeta$
dependence of the Mellin moments are visualised in figures 29 and 30. Compared
with the data entering our fit, we have slightly extended the $y$ range from
$5a$ down to $4a$.

(a) $I_{ud}(y^{2},\zeta)$ (b) $myI_{\delta u\mskip 1.5mud}(y^{2},\zeta)$ (c)
$I_{\delta u\delta d}(y^{2},\zeta)$ (d) $m^{2}|y^{2}|I^{t}_{\delta u\delta
d}(y^{2},\zeta)$

Figure 29:  Mellin moments of skewed DPDs as a function of $y$, reconstructed
from our global fit.

(a) $I_{ud}(y^{2},\zeta)/I_{ud}(y^{2},0)$ (b) $I_{\delta u\mskip
1.5mud}(y^{2},\zeta)/I_{\delta u\mskip 1.5mud}(y^{2},0)$

Figure 30:  Mellin moments of skewed DPDs as a function of $\zeta$,
reconstructed from our global fit and normalised to their value at $\zeta=0$.

In the left panel of figure 31, we show the Mellin moments at $\zeta=0$ for
the different polarisation combinations. Comparison with the data of the
corresponding twist-two functions at $py=0$ shows the close similarity between
the two quantities. This corroborates the basic assumption of our discussion
in section 4, namely that the qualitative features of twist-two functions at
$py=0$ are representative of the Mellin moments of ordinary DPDs.

With the caveats of choosing a functional ansatz and restricting ourselves to
the connected graph $C_{1}$, we can in particular extend our discussion for
light quarks in section 4.3 to the Mellin moments of DPDs for the flavour
combination $ud$ in a $\pi^{+}$: there is a substantial spin-orbit correlation
for one transversely polarised quark or antiquark, whereas correlations
involving transverse polarisation of both partons are rather small. This is
one of the main results of our work.

(a) Mellin moments $I(y^{2},\zeta=0)$ (b) twist-two functions $A(y^{2},py=0)$

Figure 31:  (a): Mellin moments of DPDs for the flavour combination $ud$ in a
$\pi^{+}$, reconstructed from our global fit. (b): Lattice data for the
corresponding twist-two functions at $py=0$. This shows the same data as
figure 20(a), but is limited to $y\leq 17a$ for ease of comparison. Notice
that the minimum $y$ in panels (a) and (b) is slightly different.

DPDs at $\zeta=0$ satisfy sum rules, which have been proposed in Gaunt:2009re
and can be proven rigorously in QCD Gaunt:2012ths ; Diehl:2018kgr . These sum
rules express momentum and quark number conservation. The quark number sum
rule for the flavour combination $ud$ in a $\pi^{+}$ implies that

$\displaystyle 2\pi\int\limits_{y_{\text{cut}}}^{\infty}\\!dy\,y\mskip
1.5muI_{ud}(y^{2};\mu)$
$\displaystyle=-1+\mathcal{O}(\Lambda^{2}y_{\text{cut}}^{2})+\mathcal{O}\bigl{(}\alpha_{s}^{2}(\mu)\bigr{)}\,,$
(65)

where $\Lambda$ denotes a hadronic scale. The necessity of a lower cutoff on
the $y$ integral and the presence of an $\mathcal{O}(\alpha_{s}^{2})$ term on
the r.h.s. result from the singular behaviour of DPDs at perturbatively small
distances $y$, as explained in Diehl:2020xyg . To avoid large logarithms in
the $\mathcal{O}(\alpha_{s}^{2})$ term, one should take $y_{\text{cut}}\sim
1/\mu$, and a standard choice is $y_{\text{cut}}=b_{0}/\mu$, where
$b_{0}=2e^{-\gamma}\approx 1.12$ and $\gamma$ is the Euler-Mascheroni
constant. With the renormalisation scale $\mu=2\operatorname{GeV}$ of our
analysis, this gives $y_{\text{cut}}\approx 0.11\operatorname{fm}\approx
1.56\mskip 1.5mua$. Extrapolating our global fit down to this value and
evaluating the integral over $y$, we obtain

$\displaystyle 2\pi\int\limits_{b_{0}/\mu}^{\infty}\\!dy\,y\mskip
1.5muI_{ud}(y^{2};\mu)=-0.915(78)\,.$ (66)

This result is not too sensitive to the extrapolation in $y$: taking an upper
integration boundary of $20a$, we obtain $-0.908(63)$, whilst raising the
lower integration boundary by a factor $2$, we obtain $-0.885(72)$. Note that
with a larger $y_{\text{cut}}$, one expects a larger
$\mathcal{O}(\Lambda^{2}y_{\text{cut}}^{2})$ term on the r.h.s. of (65). Given
the presence of this power correction in the theory prediction, we find its
agreement with our result (66) quite satisfactory. We regard this as a strong
cross check of our analysis, and in particular of the fit ansatz we have made
in (49), (62) and (63).

###  5.4 Factorisation hypothesis for Mellin moments

With the Mellin moments reconstructed from our global fit, we can also test
the factorisation hypothesis (37), which directly follows from the
corresponding hypothesis (6) for DPDs. To evaluate the r.h.s. of (37), we use
the same two fits for the vector form factor of the pion as we did in section
4.4. The comparison of the left and right-hand sides of (37), as well as their
ratio is shown in figure 32. We see the same trend as we did in figure 22 for
$A_{ud}$ at $py=0$. At small $y$, the result of the factorised ansatz is too
large in absolute size, and at large $y$ it is too small. The discrepancy at
large $y$ is even somewhat stronger for the Mellin moment $I_{ud}$ than it is
for $A_{ud}$. We draw the same conclusion as we did in section 4.4: the
factorised ansatz for the unpolarised $ud$ flavor combination in a $\pi^{+}$
can provide a rough approximation at the level of several $10\%$. In the sense
that the factorised ansatz represents the assumption that the $u$ and the
$\bar{d}$ in a $\pi^{+}$ have independent spatial distributions, our result
for $I_{ud}$ indicates that the two partons prefer to be farther apart than if
they were uncorrelated.

(a) (b)

Figure 32:  Test of the factorisation hypothesis (37) for the lowest Mellin
moment $I_{ud}$ of the unpolarised DPD $F_{ud}$ at $\zeta=0$ in a $\pi^{+}$.
(a): Comparison of $I_{ud}$ determined by the global fit of section 5.2 with
the integral over form factors on the r.h.s. of (37). The form factors are
determined by a monopole or a $p$-pole fit as described in section 4.4. (b):
Ratio of the integral over form factors and the Mellin moment. The
factorisation hypothesis predicts this ratio to be $1$.

###  5.5 Comparison with quark model results

As we mentioned in the introduction, there are a few calculations of pion DPDs
within quark models Rinaldi:2018zng ; Courtoy:2019cxq ; Broniowski:2019rmu ;
Broniowski:2020jwk . Most of the results presented in these papers are for
distributions differential in $x_{1}$ and $x_{2}$, which are not accessible in
our lattice calculation. However, reference Courtoy:2019cxq  also gives the
lowest Mellin moment of polarised and unpolarised DPDs as a function of $y$,
as we did in the present work. Predictions for the lowest Mellin moment of the
unpolarised DPD are also shown in Broniowski:2019rmu ; Broniowski:2020jwk .
They are quite similar to the ones in Courtoy:2019cxq .

The results presented in Courtoy:2019cxq  are for the physical value of the
pion mass. The authors of that work have provided us with the numbers they
obtain when taking $m_{\pi}=300\operatorname{MeV}$ instead Courtoy:2020prc .
This parameter setting was also used in Courtoy:2020tkd , where two-current
matrix elements in the pion computed in the same model were compared with the
results of our lattice study Bali:2018nde .

In figure 33, we show the quark model results along with the Mellin moments
reconstructed from our global fit.222The notation for DPDs involving
transverse polarisation in Courtoy:2019cxq  is related to our notation in
(2.1) as $F^{v}_{u\mskip 1.5mu\delta d}=myf_{u\mskip 1.5mu\delta d}$,
$F^{s}_{\delta u\delta d}=f_{\delta u\delta d}$, and $F^{t}_{\delta u\delta
d}=m^{2}|y^{2}|f^{t}_{\delta u\delta d}$, where all functions depend on
$x_{1},x_{2}$, and $y$. Compared with our figure 31(a), the moments are
multiplied with an additional factor $-y$ so as to correspond to the curves in
figure 3 of Courtoy:2020tkd . Notice that Courtoy:2020tkd  gives results for
two different regulators of ultraviolet divergences. We only show the ones
obtained with Pauli-Villars regularisation here and note that the difference
between the two regulators in figure 3 of Courtoy:2020tkd  is quite noticeable
for $I_{u\mskip 1.5mu\delta d}$ at $y<0.4\operatorname{fm}$ and for
$I^{t}_{\delta u\delta d}$ at $y>0.2\operatorname{fm}$.

(a) $-yI_{ud}$ (b) $-m|y^{2}|I_{u\mskip 1.5mu\delta}$ (c) $-yI_{\delta u\delta
d}$ (d) $-m^{2}|y|^{3}I^{t}_{\delta u\delta d}$

Figure 33:  Comparison of our results for Mellin moments $I(y,\zeta=0)$ of
DPDs in a $\pi^{+}$ with those obtained with the Nambu Jona-Lasinio model in
Courtoy:2019cxq . The lattice results correspond to the ones in our figure
31(a) but are multiplied by an additional factor $-y$. The model curves for
$m_{\pi}^{\text{phys}}$ are the same as the ones labelled “PV” in figure 3 of
Courtoy:2019cxq . The lattice and model results refer to different
factorisation scales, as specified in the text.

The difference between the model results for $m_{\pi}=300\operatorname{MeV}$
and $m_{\pi}=140\operatorname{MeV}$ is quite small for the moments shown in
figure 33. A larger mass dependence is found for $I_{\Delta u\Delta d}$, which
we do not show because our lattice data in this channel is too noisy for a
stable reconstruction of the Mellin moment.

The quark model curves in figure 33 refer to the typical renormalisation scale
of the model, which in Courtoy:2008nf  was estimated to be
$\mu=290\operatorname{MeV}$ for evolution at LO and
$\mu=430\operatorname{MeV}$ for evolution at NLO in the strong coupling. The
Mellin moment $I_{ud}$ for unpolarised quarks involves only the vector current
and is therefore scale independent, so that the model curves can be directly
compared with our lattice values at $\mu=2\operatorname{GeV}$. We find that
the results of the two approaches are remarkably close to each other, despite
a visible difference around $y\sim 0.3\operatorname{fm}$. Let us note that the
agreement between our result and the one shown for the Spectral Quark Model in
figure 6b of Broniowski:2019rmu  is even better.

Evolution from the quark model scale to $\mu=2\operatorname{GeV}$ will reduce
the moment $I_{u\mskip 1.5mu\delta d}$ by a factor $r$ and the moments
$I_{\delta u\delta d}$ and $I^{t}_{\delta u\delta d}$ by a factor $r^{2}$. We
refrain from estimating this factor here, because it involves evolution in a
region where perturbation theory becomes quite unstable. Despite this
uncertainty, we can state that for all Mellin moments involving transverse
quark polarisation, the two approaches agree in the sign and qualitative shape
for $y>0.3\operatorname{fm}$. However, it is also clear that no value of the
evolution factor $r$ can bring the lattice and model results for all three
moments into quantitative agreement.

##  6 Summary

This paper presents the first lattice calculation that provides information
about double parton distributions in a pion. Our simulations are for a pion
mass of $m_{\pi}\approx 300\operatorname{MeV}$, a lattice spacing of $a\approx
0.07\operatorname{fm}$, and two lattice volumes with $L=32$ and $L=40$ points
in the spatial lattice directions, respectively. We also have results for the
pseudoscalar ground state made of strange or of charm quarks at their physical
masses, in a partially quenched setup.

We compute the pion matrix elements of the product of two local currents that
are separated by a space-like distance. From these tensor-valued matrix
elements, we extract Lorentz invariant functions associated with the twist-two
operators in the definition of DPDs. In the continuum and infinite volume
limits, these functions depend on the pion momentum $p^{\mu}$ and the distance
$y^{\mu}$ between the currents only via the invariant products $py$ and
$y^{2}$. This allows us to detect discretisation and finite size effects in
our data, and to devise cuts that minimise these artefacts. In particular,
most results reported here are limited to distances $y$ above $5a\approx
0.35\operatorname{fm}$. Comparing the data from our two lattice volumes, we
find only mild differences in channels that have a good statistical signal.
Comparing results obtained with different source-sink separation, we find
little evidence for contributions from excited states in our analysis. The
invariant twist-two function in the axial vector channel appears to be most
strongly affected by several of the lattice artefacts.

Comparing the importance of different Wick contractions in the twist-two
functions, we find that the connected graphs $C_{1}$ and $C_{2}$ in figure 4
are the most important ones in almost all cases. For light quarks, graph
$C_{2}$ is as important as $C_{1}$ at small distances $y$ between the two
partons, which indicates that Fock states containing sea quarks are important
in that region. As one would expect, this importance is strongly reduced for
charm quarks, but it is still visible at a level below $10\%$.

We compute matrix elements for different combinations of the vector, axial
vector and tensor currents, which respectively correspond to unpolarised
partons and partons with longitudinal or transverse polarisation. For light
quarks, we find surprisingly small correlations between the longitudinal or
transverse spins of the two partons. By contrast, a large spin-orbit
correlation is seen between the transverse component of $y$ and the transverse
polarisation of one of the partons. All spin correlations increase
considerably with the quark mass, and for charm quarks we observe large spin-
spin correlations for both longitudinal and transverse polarisation.

The invariant twist-two functions that we can determine on the lattice are not
directly related to the Mellin moments of DPDs, but rather to the moments of
what can be called “skewed” DPDs. To compute the Mellin moments of ordinary
DPDs from two-current matrix elements, one needs the dependence of the
invariant functions on the variable $py$ on the full real axis. This is
inaccessible on a Euclidean lattice. Fitting an ansatz for the $py$ dependence
to our lattice data, we can however reconstruct the lowest Mellin moments by
extrapolating this ansatz to the full $py$ range. We find that the moments
obtained in this way have a behaviour very similar to the one of the twist-two
functions at $py=0$. A valuable cross check of our procedure is the fact that
the result for the unpolarised Mellin moment is in good agreement with the
number sum rule that must be obeyed by the DPD for the flavour combination
$ud$ in a $\pi^{+}$. Comparing our results for the Mellin moments with those
obtained in quark models, we find rather close agreement for unpolarised
quarks. For moments involving transverse quark polarisation, we observe
qualitative agreement but quantitative differences. We have not reconstructed
the lowest Mellin moment for longitudinal quark polarisation, because we
consider our lattice data in this channel to be too noisy for this purpose.

A starting point of many phenomenological studies is the assumption that
unpolarised DPDs can be “factorised” into the single-particle distributions of
each parton, which would mean that the two partons are independent of each
other. We have formalised this assumption and tested it, both for the twist-
two functions directly extracted from the lattice data and for the Mellin
moments reconstructed by extrapolating a fit to these data. In both cases, we
find that the two-parton correlator deviates from the factorisation ansatz by
a few $10\%$, and that the sign of the deviation depends on the transverse
distance $y$. More specifically, the two partons tend to be farther apart from
each other than if they were independent of each other.

We see several directions into which the studies reported here should be
extended. First and foremost comes the extension from a pion to a nucleon,
which is of direct relevance for double parton scattering in proton-proton
collisions. Work in this direction is underway. On a longer time scale, one
will want to have simulations with finer lattice spacing and smaller quark
masses. Data of sufficient quality for higher hadron momenta will extend the
range in $py$ that can be probed and thus allow for a better controlled
extrapolation in this variable. Given the results obtained in the present
work, we think that the efforts required for such studies will be rewarded
with valuable physics insights.

## Acknowledgements

We gratefully acknowledge input from Sara Collins and André Sternbeck, and
discussions with Michael Engelhardt. We thank A. Courtoy, S. Noguera, and S.
Scopetta for kindly providing us with the numbers of their results in
Courtoy:2019cxq  and with the corresponding values for a pion mass close to
the one in our simulations.

We used a modified version of the Chroma Edwards:2004sx  software package,
along with the locally deflated domain decomposition solver implementation of
openQCD Luscher:2012av . The gauge ensembles have been generated by the QCDSF
and RQCD collaborations on the QPACE computer using BQCD Nakamura:2010qh ;
Nakamura:2011cd . The graphs were produced with JaxoDraw Binosi:2003yf ;
Binosi:2008ig . The simulations used for this work were performed with
resources provided by the North-German Supercomputing Alliance (HLRN). This
work was supported by the Deutsche Forschungsgemeinschaft SFB/TRR 55, by BMBF
Verbundprojekt O5P2018(ErUM-FSP T01), and by the European Union’s Horizon 2020
research and innovation programme under grant agreement No. 824093
(STRONG-2020) and Marie Skłodowska-Curie grant agreement No. 813942
(EuroPLEx).

## References

  * (1) K. Barad, M. Ogilvie and C. Rebbi, _Quark-Antiquark Charge Distributions and Confinement_ , _Phys. Lett._ 143B (1984) 222.
  * (2) K. Barad, M. Ogilvie and C. Rebbi, _Quark-antiquark charge distributions_ , _Annals Phys._ 168 (1986) 284.
  * (3) W. Wilcox and K.-F. Liu, _Charge Radii From Lattice Relative Charge Distributions_ , _Phys. Lett._ B172 (1986) 62.
  * (4) W. Wilcox, K.-F. Liu, B.-A. Li and Y.-l. Zhu, _Relative Charge Distributions for Quarks in Lattice Mesons_ , _Phys. Rev._ D34 (1986) 3882.
  * (5) W. Wilcox, _Current overlap methods in lattice QCD_ , _Phys. Rev._ D43 (1991) 2443.
  * (6) M. C. Chu, M. Lissia and J. W. Negele, _Hadron structure in lattice QCD. 1\. Correlation functions and wave functions_ , _Nucl. Phys._ B360 (1991) 31.
  * (7) M. Lissia, M. C. Chu, J. W. Negele and J. M. Grandy, _Comparison of hadron quark distributions from lattice QCD and the MIT bag model_ , _Nucl. Phys._ A555 (1993) 272.