# Two-current correlations in the pion on the lattice

Gunnar S. Bali c Peter C. Bruns a Luca Castagnini d,e Markus Diehl f Jonathan
R. Gaunt g Benjamin Gläßle a Andreas Schäfer h André Sternbeck 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, IndiaNuclear Physics
Institute, 25069 Řež, Czech RepublicFachbereich 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,
GermanyTheoretisch-Physikalisches Institut, Friedrich-Schiller-Universität
Jena, 07743 Jena, Germany

###### Abstract

We perform a systematic study of the correlation functions of two quark
currents in a pion using lattice QCD. We obtain good signals for all but one
of the relevant Wick contractions of quark fields. We investigate the quark
mass dependence of our results and test the importance of correlations between
the quark and the antiquark in the pion. Our lattice data are compared with
predictions from chiral perturbation theory.

††preprint:  DESY 18-105, CERN-TH-2018-147

##  1 Introduction

The study of hadronic form factors has a long history in QCD and remains an
active field of study. Defined from the matrix elements of single currents,
form factors contain information about the spatial distribution of charge
inside a hadron. Going beyond this, the matrix elements of _two_ currents at
different points in space are sensitive to charge _correlations_ and thus
yield qualitatively new information about QCD bound states and their many-body
structure.

Hadronic matrix elements of one or two currents can be calculated in lattice
QCD. In fact, it was realised long ago that two-current correlators on the
lattice can be regarded as gauge invariant probes of hadron wave functions.
Following earlier work Barad:1984px ; Barad:1985qd ; Wilcox:1986ge ;
Wilcox:1986dk ; Wilcox:1990zc ; Chu:1990ps ; Lissia:1991gv ; Burkardt:1994pw ,
detailed computations for the vector, scalar and pseudoscalar currents were
presented in Alexandrou:2002nn ; Alexandrou:2003qt  and extended studies for
the vector current in Alexandrou:2008ru . These studies focused on a broad
range of physics and observables, such as confinement Barad:1984px ;
Barad:1985qd , the size of hadrons Wilcox:1986ge ; Wilcox:1986dk ;
Wilcox:1990zc ; Burkardt:1994pw , comparison with quark models Lissia:1991gv ,
or the non-spherical shape of hadrons with spin $1$ or larger
Alexandrou:2002nn ; Alexandrou:2003qt ; Alexandrou:2008ru .

The computation of two-current correlators on the lattice involves many
different Wick contractions of the quark fields, including disconnected graphs
and graphs with all-to-all propagators. This presents major challenges, for
the sheer amount of calculations and for obtaining statistically significant
results. Whilst the work in Barad:1984px ; Barad:1985qd ; Wilcox:1986ge ;
Wilcox:1986dk ; Wilcox:1990zc ; Chu:1990ps ; Lissia:1991gv ; Burkardt:1994pw ;
Alexandrou:2002nn ; Alexandrou:2003qt ; Alexandrou:2008ru  focused on the
graph denoted by $C_{1}$ in figure 2, we study all relevant contractions for a
meson here (leaving the case of baryons for future work). We also extend the
set of currents mentioned above by the axial current. Last but not least,
increased computing power has allowed us to use larger lattices and a smaller
pion mass than previous studies. Indeed, we will see that finite volume
effects can be appreciable for the quantities we are interested in.

At large distances between the two currents, the correlation functions we
consider can be evaluated in chiral perturbation theory, which describes the
low-energy limit of QCD in terms of mesons and their interactions. A leading-
order calculation, supplemented by an estimate of higher orders using
resonance exchange graphs, has been performed in Bruns:2015yto , and we will
compare our lattice results with these predictions.

An extension of the work presented here allows us to compute matrix elements
that, as explained in Diehl:2011yj , can be connected with the Mellin moments
of double parton distributions. These distributions are necessary to compute
coherent hard scattering on two spatially separated partons inside a hadron.
An additional challenge in this case is that one must compute correlation
functions for all vector or tensor components of the inserted currents and
then extract their twist-two part. Corresponding results for the pion will be
presented in a forthcoming paper. The ultimate goal is to extend such studies
to double parton distributions in the nucleon, which is of acute interest for
the precise description of high-multiplicity final states at the LHC and at
possible future hadron colliders.

This paper is organised as follows. In section 2, we derive a number of
general properties of two-current matrix elements of the pion and recall the
predictions of chiral perturbation theory relevant to our study. Properties of
the different Wick contractions and their computation on the lattice are
described in section 3. The general quality of our data and the presence of
lattice artefacts is investigated in section 4. Results for individual lattice
contractions are shown in section 5, including the quark mass dependence and
the relevance of correlations between the two currents. Physical matrix
elements are presented in section 6 and compared with chiral perturbation
theory. A summary of our work is given in section 7.

##  2 Correlation functions of two currents

The object of our study are correlation functions of two currents in a pion,

$\displaystyle\bra{\pi^{k}(p)}\mathcal{O}^{q_{1}q_{2}}_{i}(y)\mskip
1.5mu\mathcal{O}^{q_{3}q_{4}}_{j}(0)\ket{\pi^{k^{\prime}}(p)}\,,$ (1)

where the pion charges $k,k^{\prime}=+,-,0$ may differ between the bra and ket
state, whereas the four-momentum $p$ is the same for both. The matrix elements
(1) are understood to be fully connected, with disconnected contributions like
$\braket{\pi}{\pi}\cdot\bra{0}\mathcal{O}_{i}(y)\mskip
1.5mu\mathcal{O}_{j}(0)\ket{0}$ or
$\bra{\pi}\mathcal{O}_{i}(y)\ket{0}\cdot\bra{0}\mathcal{O}_{j}(0)\ket{\pi}$
removed. The currents we consider are

$\displaystyle\mathcal{O}^{qq^{\prime}}_{i}(y)$
$\displaystyle=\bar{q}(y)\mskip 1.5mu\Gamma_{i}\mskip 1.5muq^{\prime}(y)\,,$
(2)

where $q$ and $q^{\prime}$ are $u$ or $d$ quark fields. The full set of Dirac
matrices $\Gamma_{i}$ corresponds to the scalar, pseudoscalar, vector, axial
and tensor currents:

$\displaystyle S_{qq^{\prime}}$ $\displaystyle=\bar{q}\mskip
1.5muq^{\prime}\,,$ $\displaystyle P_{qq^{\prime}}$ $\displaystyle=i\mskip
1.5mu\bar{q}\mskip 1.5mu\gamma_{5}\mskip 1.5muq^{\prime}\,,$ $\displaystyle
V_{qq^{\prime}}^{\mu}$ $\displaystyle=\bar{q}\gamma^{\mu}\mskip
1.5muq^{\prime}\,,$ $\displaystyle A_{qq^{\prime}}^{\mu}$
$\displaystyle=\bar{q}\gamma^{\mu}\gamma_{5}\mskip 1.5muq^{\prime}\,,$
$\displaystyle T_{qq^{\prime}}^{\mu\nu}$
$\displaystyle=\bar{q}\sigma^{\mu\nu}\mskip 1.5muq^{\prime}\,.$ (3)

Note that for equal quark flavours all currents are hermitian.

In the present work, we investigate lattice data for correlation functions of
two equal currents $V^{0}$, $A^{0}$, $S$ or $P$. In a pion at rest, these
currents are associated with the vector, axial, scalar or pseudoscalar
charges. We abbreviate the corresponding matrix elements (1) as

$\displaystyle\langle V^{0}V^{0}\rangle,\langle A^{0}A^{0}\rangle,\langle
SS\rangle,\langle PP\rangle,$ (4)

respectively. For later use, we shall however consider the full set of
currents (2) in the theoretical exposition below.

The matrix elements $\langle SS\rangle$ and $\langle PP\rangle$ are boost
invariant and thus depend only on the products $y^{2}$ and $yp$ of four-
vectors, given that $p^{2}=m_{\pi}^{2}$ is fixed. On the lattice, we can
compute the matrix elements for $y^{0}=0$ and for different pion momenta. The
rotation to Euclidean time and the method to extract hadronic matrix elements
from a lattice calculation single out a particular reference frame. It is
hence a valuable check for lattice artefacts to verify that for a given
$\vec{\mskip 0.5muy}\mskip 1.5mu^{2}$, the matrix elements with scalar and
pseudoscalar currents are independent of the pion momentum when $\vec{\mskip
0.5muy}\mskip 1.5mu\cdot\vec{\mskip 0.5mup}\mskip 1.5mu=0$, a condition that
can be achieved for both zero and nonzero $\vec{\mskip 0.5mup}\mskip 1.5mu$.

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

$\displaystyle\mathcal{O}_{i}^{qq^{\prime}}(y)$
$\displaystyle\,\underset{C}{\to}\,\eta_{C}^{i}\,\mathcal{O}_{i}^{q^{\prime}\mskip-1.5muq}(y)\,,$
$\displaystyle\mathcal{O}_{i}^{qq^{\prime}}(y)$
$\displaystyle\,\underset{PT}{\to}\,\eta_{PT}^{i}\,\mathcal{O}_{i}^{q^{\prime}\mskip-1.5muq}(-y)\,,$
(5)

with sign factors

$\displaystyle\eta_{C}^{i}$ $\displaystyle=+1$
$\displaystyle\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ i$
$\displaystyle=S,P,A,$ $\displaystyle\eta_{C}^{i}$ $\displaystyle=-1$
$\displaystyle\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ i$
$\displaystyle=V,T,$ (6) and $\displaystyle\eta_{PT}^{i}$ $\displaystyle=+1$
$\displaystyle\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ i$
$\displaystyle=S,P,V,$ $\displaystyle\eta_{PT}^{i}$ $\displaystyle=-1$
$\displaystyle\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ i$
$\displaystyle=A,T.$ (7)

We also define the products $\eta_{C}^{ij}=\eta_{C}^{i\phantom{j}\\!\\!}\mskip
1.5mu\eta_{C}^{j}$ and $\eta_{PT}^{ij}=\eta_{PT}^{i\phantom{j}\\!\\!}\mskip
1.5mu\eta_{PT}^{j}$.

###  2.1 Isospin decomposition and constraints

The matrix elements (1) are not all independent due to constraints from
isospin symmetry and from the discrete symmetries just discussed. To exploit
isospin symmetry, it is convenient to use linear combinations with isospin $0$
and $1$,

$\displaystyle\mathcal{O}^{\text{s}}_{i}$
$\displaystyle=\mathcal{O}^{uu}_{i}+\mathcal{O}^{dd}_{i}\,,$
$\displaystyle\mathcal{O}^{\text{ns}}_{i}$
$\displaystyle=\mathcal{O}^{uu}_{i}-\mathcal{O}^{dd}_{i}\,,$ (8)

as well as the isotriplet current

$\displaystyle\mathcal{O}^{a}_{i}$ $\displaystyle=\bar{Q}\mskip
1.5mu\tau^{a}\Gamma_{i}\mskip 1.5muQ\,,$ $\displaystyle a$
$\displaystyle=1,2,3$ (9)

with the Pauli matrices $\tau^{a}$ and the isodoublet $Q=(u,d)$ of quark
fields. One then has

$\displaystyle\mathcal{O}^{\text{ns}}_{i}$
$\displaystyle=\mathcal{O}^{3}_{i}\,,$ $\displaystyle\mathcal{O}^{ud}_{i}$
$\displaystyle=\bigl{(}\mathcal{O}^{1}_{i}+i\mathcal{O}^{2}_{i}\bigr{)}/{2}\,,$
$\displaystyle\mathcal{O}^{du}_{i}$
$\displaystyle=\bigl{(}\mathcal{O}^{1}_{i}-i\mathcal{O}^{2}_{i}\bigr{)}/{2}\,.$
(10)

Expressing the pion states in the isospin basis, we have

$\displaystyle|\pi^{+}\rangle$
$\displaystyle=\bigl{(}|\pi^{1}\rangle+i|\pi^{2}\rangle\bigr{)}/\sqrt{2}\,,$
$\displaystyle|\pi^{-}\rangle$
$\displaystyle=\bigl{(}|\pi^{1}\rangle-i|\pi^{2}\rangle\bigr{)}/\sqrt{2}\,,$
$\displaystyle|\pi^{0}\rangle$ $\displaystyle=|\pi^{3}\rangle\,.$ (11)

We can now decompose

$\displaystyle\bra{\pi^{d}}\mathcal{O}_{i}^{\text{s}}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{s}}(0)\ket{\pi^{c}}$
$\displaystyle=\delta^{cd}\mskip 1.5muF_{0}(y)\,,$
$\displaystyle\bra{\pi^{d}}\mathcal{O}_{i}^{a}(y)\mskip
1.5mu\mathcal{O}_{j}^{b}(0)\ket{\pi^{c}}$
$\displaystyle=\delta^{ab}\delta^{cd}\mskip
1.5muF_{1}(y)+\bigl{(}\delta^{ac}\delta^{bd}+\delta^{ad}\delta^{bc}\bigr{)}\mskip
1.5muF_{2}(y)+\bigl{(}\delta^{ac}\delta^{bd}-\delta^{ad}\delta^{bc}\bigr{)}\mskip
1.5muiF_{3}(y)\,,$ $\displaystyle\bra{\pi^{d}}\mathcal{O}_{i}^{b}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{s}}(0)\ket{\pi^{c}}$
$\displaystyle=i\epsilon^{bcd}\mskip 1.5muG_{1}(y)\,,$
$\displaystyle\bra{\pi^{d}}\mathcal{O}_{i}^{\text{s}}(y)\mskip
1.5mu\mathcal{O}_{j}^{b}(0)\ket{\pi^{c}}$ $\displaystyle=i\epsilon^{bcd}\mskip
1.5muG_{2}(y)\,,$ (12)

where for brevity we have suppressed the dependence of $F_{n}$, $G_{n}$ on the
pion momentum $p$ and on the indices $i,j$ that specify the operators. Taking
the complex conjugate of (2.1), one readily finds that the isospin amplitudes
$F_{n}$ and $G_{n}$ are real valued.

For the matrix elements in the quark flavor basis, we obtain

$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{\text{s}}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{s}}(0)\ket{\pi^{+}}=\bra{\pi^{0}}\mathcal{O}_{i}^{\text{s}}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{s}}(0)\ket{\pi^{0}}$ $\displaystyle=F_{0}\,,$
$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{\text{ns}}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{s}}(0)\ket{\pi^{+}}=-\sqrt{2}\mskip
1.5mu\bra{\pi^{0}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{s}}(0)\ket{\pi^{+}}$ $\displaystyle=G_{1}\,,$
$\displaystyle\bra{\pi^{0}}\mathcal{O}_{i}^{\text{ns}}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{s}}(0)\ket{\pi^{0}}$ $\displaystyle=0$ (13)

and analogous relations for the matrix elements parameterised by $G_{2}$, now
also suppressing the $y$-dependence of $F_{n}$ and $G_{n}$. For matrix
elements with two isotriplet operators, we have

$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{\text{ns}}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{ns}}(0)\ket{\pi^{+}}$ $\displaystyle=2\mskip
1.5mu\bra{\pi^{0}}\mathcal{O}_{i}^{ud}(y)\mskip
1.5mu\mathcal{O}_{j}^{du}(0)\ket{\pi^{0}}=2\mskip
1.5mu\bra{\pi^{0}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{ud}(0)\ket{\pi^{0}}=F_{1}\,,$
$\displaystyle\bra{\pi^{0}}\mathcal{O}_{i}^{\text{ns}}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{ns}}(0)\ket{\pi^{0}}$
$\displaystyle=F_{1}+2F_{2}\,,$
$\displaystyle\bra{\pi^{-}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{du}(0)\ket{\pi^{+}}$ $\displaystyle=F_{2}\,,$
$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{ud}(y)\mskip
1.5mu\mathcal{O}_{j}^{du}(0)\ket{\pi^{+}}$
$\displaystyle=(F_{1}+F_{2}-iF_{3})/2\,,$
$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{ud}(0)\ket{\pi^{+}}$
$\displaystyle=(F_{1}+F_{2}+iF_{3})/2\,,$
$\displaystyle\bra{\pi^{0}}\mathcal{O}_{i}^{\text{ns}}(y)\mskip
1.5mu\mathcal{O}_{j}^{du}(0)\ket{\pi^{+}}$
$\displaystyle=(F_{2}-iF_{3})/\sqrt{2}\,,$
$\displaystyle\bra{\pi^{0}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{ns}}(0)\ket{\pi^{+}}$
$\displaystyle=(F_{2}+iF_{3})/\sqrt{2}\,.$ (14)

From (2.1) one also finds that the matrix elements with two isosinglet or with
two isotriplet operators are even under the simultaneous exchange111According
to (10) and (11) this corresponds to changing the sign of $|\pi^{a}\rangle$
and $\mathcal{O}^{a}_{i}$ for $a=2$ in the isospin basis.

$\displaystyle|\pi^{+}\rangle$
$\displaystyle\leftrightarrow|\pi^{-}\rangle\,,$
$\displaystyle\mathcal{O}^{du}_{i}$
$\displaystyle\leftrightarrow\mathcal{O}^{ud}_{i}\,,$ (15)

whereas the matrix elements with one isosinglet and one isotriplet operator
are odd under that exchange. The charge conjugate of a matrix element is
obtained by applying (15) and multiplying with the product $\eta_{C}^{ij}$ of
intrinsic $C$ parities of the two operators. One readily sees that $G_{n}=0$
if $\eta_{C}^{ij}=+1$ and that $F_{n}=0$ if $\eta_{C}^{ij}=-1$. The charge
conjugation constraints give relations between matrix elements in the quark
flavor basis. For $\eta_{C}^{ij}=+1$ the vanishing of $G_{1}$ and $G_{2}$
implies in particular

$\displaystyle\bra{\pi^{k}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{k}}$
$\displaystyle=\bra{\pi^{k}}\mathcal{O}_{i}^{dd}(y)\mskip
1.5mu\mathcal{O}_{j}^{dd}(0)\ket{\pi^{k}}\,,$
$\displaystyle\bra{\pi^{k}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{dd}(0)\ket{\pi^{k}}$
$\displaystyle=\bra{\pi^{k}}\mathcal{O}_{i}^{dd}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{k}}$ (16) and thus
$\displaystyle{\textstyle\frac{1}{2}}\mskip
1.5mu\bra{\pi^{k}}\mathcal{O}_{i}^{\text{s}}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{s}}(0)\ket{\pi^{k}}$
$\displaystyle=\bra{\pi^{k}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{k}}+\bra{\pi^{k}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{dd}(0)\ket{\pi^{k}}\,,$
$\displaystyle{\textstyle\frac{1}{2}}\mskip
1.5mu\bra{\pi^{k}}\mathcal{O}_{i}^{\text{ns}}(y)\mskip
1.5mu\mathcal{O}_{j}^{\text{ns}}(0)\ket{\pi^{k}}$
$\displaystyle=\bra{\pi^{k}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{k}}-\bra{\pi^{k}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{dd}(0)\ket{\pi^{k}}$ (17)

for $k=+,-,0$.

###  2.2 Predictions from chiral perturbation theory

Consider pions at rest and take the matrix elements $\langle SS\rangle,\langle
PP\rangle,\langle V^{0}V^{0}\rangle$ and $\langle A^{0}A^{0}\rangle$ at
$y^{0}=0$ and large $|\vec{\mskip 0.5muy}\mskip 1.5mu|$. These can then be
computed in chiral perturbation theory. The chiral expansion, on which this
theory is based, requires the pion mass and momenta $p$ to be much smaller
than $4\pi F\sim 1\operatorname{GeV}$, where $F$ is the pion decay constant.
In position space, one should hence require $|\vec{\mskip 0.5muy}\mskip
1.5mu|\gg 0.2\operatorname{fm}$. At leading order in the chiral expansion, the
matrix elements can be computed from the tree-level graphs in figure 1(a), (b)
and (c). The corresponding calculation is detailed in Bruns:2015yto . As an
estimate for higher-order contributions, the same work evaluated the resonance
exchange graphs (d) and (e) from the appropriate leading-order Lagrangian in
the approximation of vanishing pion four-momentum (thus setting $m_{\pi}$ to
zero). Resonance exchange graphs with the topology of figure 1(c) would
involve a vertex between pions and resonances, which vanishes in this
approximation. Considered were the lowest-lying resonances $\rho$, $a_{1}$,
$a_{0}$, $\eta$ and $\sigma$.

Figure 1:  Graphs for the matrix elements $\langle SS\rangle,\langle
PP\rangle,\langle V^{0}V^{0}\rangle,\langle A^{0}A^{0}\rangle$ in chiral
perturbation theory. Single lines denote pions, double lines resonances, and
crossed circles indicate current insertions. The leading-order chiral
Lagrangian gives rise to graph (a) for $\langle SS\rangle,\langle
V^{0}V^{0}\rangle$ and to graphs (b) and (c) for $\langle PP\rangle,\langle
A^{0}A^{0}\rangle$. The resonance exchange graphs (d) and (e) can contribute
to any of the four matrix elements, depending on the quantum numbers of the
resonance.

The corresponding results are given in sections II.A and III.B of
Bruns:2015yto . For convenience, we recast them into the notation used here,
making use of the simplifications that arise when setting $y^{0}=0$. The
isospin amplitudes $F_{i}$ defined here are related to those in Bruns:2015yto
as $F_{0}=C^{00}$, $F_{1}=C_{1}$, $F_{2}=C_{1}+C_{2}$ and $iF_{3}=C_{3}$. We
also note that the isospin currents $V^{a}_{\mu}$ and $A^{a}_{\mu}$ in
Bruns:2015yto  have an extra factor of $1/2$ compared with our currents (9).

Vector and axial currents were only considered for the isovector case in
Bruns:2015yto , so that no predictions are available for $F_{0}$ in the
$\langle V^{0}V^{0}\rangle$ and $\langle A^{0}A^{0}\rangle$ channels. Under
the conditions stated above, the isotriplet amplitude $F_{3}$ is found to
vanish for all channels $\langle SS\rangle$, $\langle PP\rangle$, $\langle
V^{0}V^{0}\rangle$, $\langle A^{0}A^{0}\rangle$.

For the remaining matrix elements, we obtain the following results from
Bruns:2015yto , writing $y=|\vec{\mskip 0.5muy}\mskip 1.5mu|$ for simplicity.
The leading order chiral Lagrangian gives rise to the amplitudes

$\displaystyle F_{0\mskip 1.5mu(\mathrm{LO})}^{{SS}}$
$\displaystyle=\frac{2B^{2}m_{\pi}}{\pi^{2}y}\,K_{1}(m_{\pi}\mskip 1.5muy)\,,$
$\displaystyle F_{0\mskip 1.5mu(\mathrm{LO})}^{{PP}}$ $\displaystyle=0\,,$
$\displaystyle F_{1\mskip 1.5mu(\mathrm{LO})}^{{SS}}$ $\displaystyle=0\,,$
$\displaystyle F_{1\mskip 1.5mu(\mathrm{LO})}^{{PP}}$
$\displaystyle=-\frac{B^{2}m_{\pi}^{2}}{2\pi^{2}}\,K_{0}(m_{\pi}\mskip
1.5muy)\,,$ $\displaystyle F_{2\mskip 1.5mu(\mathrm{LO})}^{{SS}}$
$\displaystyle=0\,,\phantom{\frac{1}{1}}$ $\displaystyle F_{2\mskip
1.5mu(\mathrm{LO})}^{{PP}}$
$\displaystyle=\frac{B^{2}m_{\pi}^{2}}{2\pi^{2}}\left[K_{0}(m_{\pi}y)-\frac{2K_{1}(m_{\pi}\mskip
1.5muy)}{m_{\pi}\mskip 1.5muy}\right]$ (18)

and

$\displaystyle F_{1\mskip 1.5mu(\mathrm{LO})}^{{VV}}$
$\displaystyle=\frac{2m_{\pi}^{3}}{\pi^{2}y}\left[K_{1}(m_{\pi}\mskip
1.5muy)-\frac{K_{2}(m_{\pi}\mskip 1.5muy)}{m_{\pi}\mskip 1.5muy}\right]\,,$
$\displaystyle F_{1\mskip 1.5mu(\mathrm{LO})}^{{AA}}$
$\displaystyle=\frac{m_{\pi}^{3}}{2\pi^{2}y}\left[K_{1}(m_{\pi}\mskip
1.5muy)+\frac{4K_{2}(m_{\pi}\mskip 1.5muy)}{m_{\pi}\mskip 1.5muy}\right]\,,$
$\displaystyle F_{2\mskip 1.5mu(\mathrm{LO})}^{{VV}}$
$\displaystyle=-\frac{1}{2}F_{1\mskip
1.5mu(\mathrm{LO})}^{\mathrm{VV}}\,,\phantom{\left[\frac{0}{0}\right]}$
$\displaystyle F_{2\mskip 1.5mu(\mathrm{LO})}^{{AA}}$
$\displaystyle=-\frac{m_{\pi}^{3}}{2\pi^{2}y}\left[K_{1}(m_{\pi}\mskip
1.5muy)+\frac{2K_{2}(m_{\pi}\mskip 1.5muy)}{m_{\pi}\mskip 1.5muy}\right]\,,$
(19)

where $K_{n}$ denotes the Macdonald functions (modified Bessel functions of
the second kind). The pion decay constant $F$ and the chiral symmetry breaking
parameter $B$ are defined as usual in chiral perturbation theory, see section
I.A in Bruns:2015yto . The resonance exchange contributions read

$\displaystyle F_{0\mskip 1.5mu(\mathrm{R})}^{{SS}}$
$\displaystyle=-\frac{32B^{2}(c_{m}^{\sigma})^{2}\,m_{\sigma}}{\pi^{2}F^{2}\mskip
1.5muy}\,K_{1}(m_{\sigma}\mskip 1.5muy)\,,$ $\displaystyle F_{0\mskip
1.5mu(\mathrm{R})}^{{PP}}$ $\displaystyle=\frac{16B^{2}}{\pi^{2}F^{2}\mskip
1.5muy}\left[c_{m}^{2}\mskip 1.5mum_{a_{0}}\,K_{1}(m_{a_{0}}\mskip
1.5muy)-2\mskip 1.5mud_{\eta}^{2}\mskip 1.5mum_{\eta}\,K_{1}(m_{\eta}\mskip
1.5muy)\right]\,,$ $\displaystyle F_{1\mskip 1.5mu(\mathrm{R})}^{{SS}}$
$\displaystyle=0\,,\qquad\quad F_{1\mskip
1.5mu(\mathrm{R})}^{{PP}}=0\,,\qquad\quad F_{2\mskip
1.5mu(\mathrm{R})}^{{SS}}=-\frac{1}{2}\mskip 1.5muF_{0\mskip
1.5mu(\mathrm{R})}^{{PP}}\,,\qquad\quad F_{2\mskip
1.5mu(\mathrm{R})}^{{PP}}=-\frac{1}{2}\mskip 1.5muF_{0\mskip
1.5mu(\mathrm{R})}^{{SS}}\,,$ $\displaystyle F_{1\mskip
1.5mu(\mathrm{R})}^{{VV}}$ $\displaystyle=-F_{1\mskip
1.5mu(\mathrm{R})}^{{AA}}=\sum_{X=\rho,\mskip
1.5mua_{1}}\\!s_{X}\,\frac{2f_{X}^{2}\mskip 1.5mum_{X}^{5}}{\pi^{2}F^{2}\mskip
1.5muy}\left[K_{1}(m_{X}\mskip 1.5muy)+\frac{K_{2}(m_{X}\mskip
1.5muy)}{m_{X}\mskip 1.5muy}\right]\,,$ $\displaystyle F_{2\mskip
1.5mu(\mathrm{R})}^{{VV}}$ $\displaystyle=-F_{2\mskip
1.5mu(\mathrm{R})}^{{AA}}=-\frac{1}{2}\mskip 1.5muF_{1\mskip
1.5mu(\mathrm{R})}^{{VV}}$ (20)

with signs $s_{\rho}=1,s_{a_{1}}=-1$. In the numerical evaluation of section
6.2, we will take the resonance parameters estimated in Bruns:2015yto :

$\displaystyle m_{\rho}$ $\displaystyle=0.8\operatorname{GeV}\,,$
$\displaystyle f_{\rho}$ $\displaystyle=0.2\,,$ $\displaystyle m_{a_{1}}$
$\displaystyle=1.25\operatorname{GeV}\,,$ $\displaystyle f_{a_{1}}$
$\displaystyle=0.1\,,$ $\displaystyle m_{a_{0}}$
$\displaystyle=1\operatorname{GeV}\,,$ $\displaystyle c_{m}$
$\displaystyle=50\operatorname{MeV}\,,$ $\displaystyle m_{\eta}$
$\displaystyle=600\operatorname{MeV}\,,$ $\displaystyle d_{\eta}$
$\displaystyle=15\operatorname{MeV}\,,$ $\displaystyle m_{\sigma}$
$\displaystyle=0.5\operatorname{GeV}\,,$ $\displaystyle c_{m}^{\sigma}$
$\displaystyle=35\operatorname{MeV}\,.$ (21)

For the parameters in the chiral sector we will take
$m_{\pi}=300\operatorname{MeV}$, $F=100\operatorname{MeV}$ and
$B=2.4\operatorname{GeV}$, which are rounded values of what we extract from
our lattice simulations with $L=40$, see (5.5).

##  3 Lattice techniques

The lattice computation of the matrix elements (1) involves a considerable
number of different Wick contractions between the quark fields in the two
currents and in the pion source and sink. This is to be contrasted with, e.g.,
the computation of single-current matrix elements, where there is only one
connected and one disconnected graph. In the present section, we give details
about the lattice contractions for two-current correlators, their relation
with the physical pion matrix elements, and their implementation in the
lattice simulation.

###  3.1 Lattice contractions and their symmetry properties

The different lattice contractions are pictorially represented in figure 2. We
have two fully connected graphs $C_{1}$ and $C_{2}$, a graph $A$ in which the
pion in annihilated by one current and created again by the second current,
two graphs $S_{1}$ and $S_{2}$ with one disconnected quark loop, and a doubly
disconnected graph $D$ with two quark loops.

The symbols $C_{1}^{ij}(y)$, $C_{2}^{ij}(y)$, etc. denote contributions to the
physical matrix elements (1), i.e. it is understood that

  * •

the contractions shown are averaged over the gauge ensemble and divided by the
ensemble average of the pion two-point function,

  * •

the two currents are inserted at the same time slice $\tau$ (i.e. $y^{0}=0$),
and the ratio of four-point and two-point functions is evaluated at a plateau
in $\tau$. To use time reversal invariance, the plateau must be symmetric
between source and sink times,

  * •

both source and sink are projected on definite momentum $p$. By translation
invariance one can thus shift their spatial positions by a common amount.

Details are given in section 3.3 below. For brevity, we do not indicate the
dependence of $C_{1}^{ij}(y)$ etc. on the pion momentum $p$.

Figure 2:  Lattice contractions for the two-current correlation function (1)
in a pion. The dependence on the pion momentum $p$ is not indicated for
brevity. The product $\smash{\eta_{C}^{ij}}$ of charge conjugation parities of
the two operators is defined below (7).

From translation invariance one readily finds

$\displaystyle M^{ij}(y)$ $\displaystyle=M^{ji}(-y)$ $\displaystyle\text{for
$M=S_{2},D$}.$ (22) Charge conjugation invariance gives $\displaystyle
M^{ij}(y)$ $\displaystyle=\eta_{C}^{ij}\,M^{ji}(-y)$ $\displaystyle\text{for
$M=C_{1}$}.$ (23)

If $\eta_{C}^{ij}=-1$ one furthermore obtains $M^{ij}(y)=0$ for $M=A,S_{2},D$.
From $PT$ invariance one gets

$\displaystyle M^{ij}(y)$ $\displaystyle=\eta_{PT}^{ij}\mskip 1.5muM^{ij}(-y)$
$\displaystyle\text{for $M=C_{1},S_{1},S_{2},D$},$ (24) $\displaystyle
M^{ij}(y)$ $\displaystyle=\eta_{PT}^{ij}\mskip 1.5muM^{ji}(y)$
$\displaystyle\text{for $M=A,C_{2}$},$ (25)

where for the formulation of time reversal in the Euclidean path integral we
refer to Hasenfratz:2005ch .

Finally, a useful relation for the contractions is obtained by taking the
complex conjugate of the contraction, followed by a parity transformation and
charge conjugation. This gives

$\displaystyle\bigl{[}M^{ij}(y)\bigr{]}^{*}$ $\displaystyle=\eta_{PT}\mskip
1.5muM^{ij}(-y)$ for all contractions. (26)

Combining (24) with (26) we find that $C_{1},S_{1},S_{2}$ and $D$ are real
valued, whereas combining (25) with (26) gives

$\displaystyle\operatorname{Re}M^{ij}(y)$
$\displaystyle={\textstyle\frac{1}{2}}\bigl{[}M^{ij}(y)+M^{ji}(-y)\bigr{]}\,,$
$\displaystyle i\operatorname{Im}M^{ij}(y)$
$\displaystyle={\textstyle\frac{1}{2}}\bigl{[}M^{ij}(y)-M^{ji}(-y)\bigr{]}$
(27)

for $M=A,C_{2}$.

###  3.2 Physical matrix elements

The form of the relation between pion matrix elements and lattice contractions
depends on the product of $C$ parities. 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)\,,$ (28)

we have for $\eta_{C}^{ij}=+1$

$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{dd}(0)\ket{\pi^{+}}$
$\displaystyle=C_{1}+\bigl{[}2S_{1}+D\bigr{]}\,,$
$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{+}}$
$\displaystyle=\bigl{[}2C_{2}+S_{2}\bigr{]}+\bigl{[}2S_{1}+D\bigr{]}\,,$
$\displaystyle\bra{\pi^{0}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{dd}(0)\ket{\pi^{0}}$
$\displaystyle=\bigl{[}2S_{1}+D\bigr{]}-A\,,$
$\displaystyle\bra{\pi^{0}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{0}}$
$\displaystyle=C_{1}+\bigl{[}2S_{1}+D\bigr{]}+\bigl{[}2C_{2}+S_{2}\bigr{]}+A\,,$
$\displaystyle\bra{\pi^{0}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{ud}(0)\ket{\pi^{0}}$
$\displaystyle={}-C_{1}+\bigl{[}2C_{2}+S_{2}\bigr{]}\,,$
$\displaystyle\bra{\pi^{-}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{du}(0)\ket{\pi^{+}}$ $\displaystyle=2C_{1}+2A\,,$
$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{ud}(0)\ket{\pi^{+}}$
$\displaystyle=2C_{2}^{ij}(y)+S_{2}+A^{ij}(y)\,,$ $\displaystyle\sqrt{2}\mskip
1.5mu\bra{\pi^{0}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{+}}$
$\displaystyle=C_{1}+\bigl{[}C_{2}^{ij}(y)-C_{2}^{ji}(-y)\bigr{]}+A^{ij}(y)\,.$
(29)

One can verify that this satisfies the isospin relations in section 2.1. We
have three independent real valued combinations, which can for instance be
chosen as

$\displaystyle{\textstyle\frac{1}{2}}F_{0}(y)$
$\displaystyle=\bra{\pi^{+}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{+}}+\bra{\pi^{+}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{dd}(0)\ket{\pi^{+}}$ $\displaystyle=C_{1}+2\mskip
1.5mu\bigl{[}2S_{1}+D\bigr{]}+[2C_{2}+S_{2}]\,,$
$\displaystyle{\textstyle\frac{1}{2}}F_{1}(y)$
$\displaystyle=\bra{\pi^{+}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{+}}-\bra{\pi^{+}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{dd}(0)\ket{\pi^{+}}$
$\displaystyle={}-C_{1}+\bigl{[}2C_{2}+S_{2}\bigr{]}\,,$
$\displaystyle{\textstyle\frac{1}{2}}F_{2}(y)$
$\displaystyle={\textstyle\frac{1}{2}}\bra{\pi^{-}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{du}(0)\ket{\pi^{+}}$ $\displaystyle=C_{1}+A\,,$ (30)

where $F_{0}$, $F_{1}$ and $F_{2}$ are the isospin amplitudes defined in
(2.1). These combinations are real valued owing to (27), as the isospin
amplitudes must be. The imaginary combination $iF_{3}$ of matrix elements can
e.g. be isolated by

$\displaystyle iF_{3}(y)$
$\displaystyle=\bra{\pi^{+}}\mathcal{O}_{i}^{du}(y)\mskip
1.5mu\mathcal{O}_{j}^{ud}(0)\ket{\pi^{+}}-\bra{\pi^{+}}\mathcal{O}_{i}^{ud}(y)\mskip
1.5mu\mathcal{O}_{j}^{du}(0)\ket{\pi^{+}}$
$\displaystyle=2\bigl{[}C_{2}^{ij}(y)-C_{2}^{ji}(-y)\bigr{]}+\bigl{[}A^{ij}(y)-A^{ji}(-y)\bigr{]}\,.$
(31)

In the present work, we only study the real valued combinations in (3.2).

For $\eta_{C}^{ij}=-1$ we have three nonzero lattice contractions, $C_{1}$,
$C_{2}$ and $S_{1}$, and two independent matrix elements, which may be taken
as

$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{dd}(0)\ket{\pi^{+}}$
$\displaystyle=C_{1}+S_{1}^{ij}(y)-S_{1}^{ji}(y)\,,$
$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{j}^{uu}(0)\ket{\pi^{+}}$
$\displaystyle=2C_{2}+S_{1}^{ij}(y)+S_{1}^{ji}(y)\,.$ (32)

Among all contractions, $C_{1}$ has typically the smallest statistical
uncertainties in the lattice simulation and thus plays a special role. The
above equations show that it does not appear in isolation in any physical pion
matrix element. However, $C_{1}$ can be regarded as “approximately physical”
in the following sense. Consider QCD with $n_{F}=4$ (instead of $n_{F}=2$)
mass degenerate quarks, where SU(4) flavour symmetry is exact. It is easy to
see that in this theory

$\displaystyle\bra{\pi^{+}}\mathcal{O}_{i}^{uc}(y)\mskip
1.5mu\mathcal{O}_{j}^{sd}(0)\ket{D_{s}^{+}}$ $\displaystyle=C_{1}^{ij}(y)\,,$
(33)

where the quark content of the $D_{s}^{+}$ meson is $c\bar{s}$. The
correlation function $C_{1}$ computed in our study does not exactly correspond
to this matrix element, because our lattice action has $n_{F}=2$ and not
$n_{F}=4$ sea quark flavours. We therefore must interpret (33) within the
partially quenched approximation.

###  3.3 Simulation Details

#### Lattice action and quark mass values.

We have performed $N_{f}=2$ lattice simulations with the Wilson gauge action
and non-perturbatively improved Sheikholeslami-Wohlert (NPI Wilson-clover)
fermions. The gauge configurations were generated by the RQCD and QCDSF
collaborations. As is discussed for instance in Bali:2014nma ; Bali:2014gha ,
there are eleven standard ensembles with pion masses down to
$150\operatorname{MeV}$. Two of these are used in the present study, namely
ensembles IV and V, with a reduced number of gauge configurations as indicated
in table 1.

ensemble | $\beta$ | $a\,[\operatorname{fm}]$ | $\kappa$ | $L^{3}\times T$ | $m_{\pi}\,[\operatorname{MeV}]$ | $Lm_{\pi}$ | $N_{\text{full}}$ | $N_{\text{used}}$ | $N_{\text{sm}}$
---|---|---|---|---|---|---|---|---|---
IV | 5.29 | 0.071 | 0.13632 | $32^{3}\times 64$ | 294.6(14) | 3.42 | 2023 | 960 | 400
V | 5.29 | 0.071 | 0.13632 | $40^{3}\times 64$ | 288.8(11) | 4.19 | 2025 | 984 | 400
Table 1:  Details of the ensembles used in this analysis. $N_{\text{full}}$ is
the total number of available gauge configurations and $N_{\text{used}}$ the
number of configurations used in the present study. $N_{\text{sm}}$ indicates
the number of Wuppertal smearing iterations in the pion source. The error on
the pion mass combines statistical and systematic errors, see Bali:2014nma .
The time difference between pion source and sink is $t=15a$ or $t=32a$ for
both ensembles.

We will also study the dependence of the pion matrix elements (1) on the quark
mass. To this end, we have performed simulations with ensemble V and different
values of $\kappa$ in the quark propagator, namely

light quarks: $\displaystyle\kappa$ $\displaystyle=0.13632$ $\displaystyle
am_{q}=0.00291(3)\,,$ strange: $\displaystyle\kappa$ $\displaystyle=0.135616$
$\displaystyle am_{q}=0.02195(3)\,,$ charm: $\displaystyle\kappa$
$\displaystyle=0.125638$ $\displaystyle am_{q}=0.31475(3)\,,$ (34)

where we also give the values of the bare quark masses $m_{q}$ in lattice
units. Here “light quarks” refers to the $\kappa$ value of ensemble V, 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.

#### Correlation functions.

To compute the pion matrix element in (1) on a lattice in Euclidean space-
time, we need the four-point correlators

$\displaystyle C^{ij,\vec{\mskip 0.5mup}\mskip
1.5mu}_{4\mathrm{pt}}(\vec{\mskip 0.5muy}\mskip 1.5mu,\tau,t)$
$\displaystyle=a^{6}\mskip 1.5mu\sum_{\vec{\mskip 0.5mux}\mskip
1.5mu,\vec{\mskip 0.5muz}\mskip 1.5mu}e^{-i\vec{\mskip 0.5mup}\mskip
1.5mu\cdot(\vec{\mskip 0.5mux}\mskip 1.5mu-\vec{\mskip 0.5muz}\mskip
1.5mu)}\,\bra{0}\Pi(\vec{\mskip 0.5mux}\mskip
1.5mu,t)\,\mathcal{O}_{i}(\vec{\mskip 0.5muy}\mskip
1.5mu,\tau)\,\mathcal{O}_{j}(\vec{\mskip 0.5mu0}\mskip
1.5mu,\tau)\,\Pi^{\dagger}(\vec{\mskip 0.5muz}\mskip 1.5mu,0)\ket{0}\,,$ (35)

where $\Pi(x)$ is a pion interpolator. We use interpolators with Dirac
structure $\gamma_{5}$. Pion matrix elements are extracted by calculating the
following ratio at time slices where excited states should be suppressed:

$\bra{\pi(p)}\mathcal{O}_{i}(y)\mskip
1.5mu\mathcal{O}_{j}(0)\ket{\pi(p)}=R_{ij}^{\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5muy}\mskip 1.5mu)=2E_{\vec{\mskip 0.5mup}\mskip
1.5mu}\mskip 1.5muV\,\frac{C^{ij,\vec{\mskip 0.5mup}\mskip
1.5mu}_{4\mathrm{pt}}(\vec{\mskip 0.5muy}\mskip 1.5mu,\tau,t)}{C^{\vec{\mskip
0.5mup}\mskip 1.5mu}_{2\mathrm{pt}}(t)}\;\biggl{|}_{0\ll\tau\ll t}\,.$ (36)

Here $V=L^{3}\mskip 1.5mua^{3}$ is the spatial volume and

$\displaystyle C^{\vec{\mskip 0.5mup}\mskip 1.5mu}_{2\mathrm{pt}}(t)$
$\displaystyle=a^{6}\mskip 1.5mu\sum_{\vec{\mskip 0.5mux}\mskip
1.5mu,\vec{\mskip 0.5muz}\mskip 1.5mu}e^{-i\vec{\mskip 0.5mup}\mskip
1.5mu\cdot(\vec{\mskip 0.5mux}\mskip 1.5mu-\vec{\mskip 0.5muz}\mskip
1.5mu)}\bra{0}\Pi(\vec{\mskip 0.5mux}\mskip
1.5mu,t)\,\Pi^{\dagger}(\vec{\mskip 0.5muz}\mskip 1.5mu,0)\ket{0}$ (37)

the usual pion two-point function. In analogy to (35) and (36), we use an
appropriate ratio of the three-point function $\smash{C^{i,\vec{\mskip
0.5mup}\mskip 1.5mu}_{3\mathrm{pt}}(\tau,t)}$ and $C^{\vec{\mskip
0.5mup}\mskip 1.5mu}_{2\mathrm{pt}}(t)$ to extract the matrix elements
$\bra{\pi(p^{\prime})}\mathcal{O}_{i}(0)\ket{\pi(p)}$ of the vector and scalar
currents. We thus obtain the vector and scalar form factors for the pion,
which are used in sections 3.4 and 5.3.

The pion energy in (36) is computed using the continuum dispersion relation
$E_{\vec{\mskip 0.5mup}\mskip 1.5mu}=\sqrt{\vec{\mskip 0.5mup}\mskip
1.5mu^{2}+m_{\pi}^{2}}$, which is well satisfied for the momenta investigated
here. The value of $m_{\pi}$ is obtained from an exponential fit of the two-
point function, which gives

$\displaystyle 293\operatorname{MeV}$
$\displaystyle\quad(\text{light},L=40)\,,$ $\displaystyle
299\operatorname{MeV}$ $\displaystyle\quad(\text{light},L=32)\,,$
$\displaystyle 691\operatorname{MeV}$
$\displaystyle\quad(\text{strange},L=40)\,,$ $\displaystyle
3018\operatorname{MeV}$ $\displaystyle\quad(\text{charm},L=40)\,.$ (38)

The statistical errors of the fits are less than 1%, and since we do not aim
at a high-precision analysis, we do not attempt to quantify the uncertainty
due to excited states. We observe that the values in (3.3) agree reasonably
well with the pion masses given in table 1 for light quarks and quoted after
(3.3) for strange quarks.

The source-sink distance is fixed to $t=15a\approx 1.07\operatorname{fm}$ as a
default. To investigate the influence of excited states, we have also
calculated graphs $C_{1}$, $C_{2}$ and $S_{1}$ with $t/a=32$. To extract the
desired matrix elements, we calculate the ratio (36) by fitting or averaging
over the plateau around $\tau=t/2$. For the contractions $C_{1}$ and $A$, we
measure the $\tau$ dependence of $C_{4\mathrm{pt}}$ and fit to a plateau in
the $\tau$ ranges specified in (4.1). For $t/a=15$ we have $S_{2}$ data for
$\tau/a=7$ and $8$, which we average, whereas for the remaining contractions
$C_{2},S_{2}$ and $D$ we take $\tau/a=7$ or $8$ for each gauge configuration
on a statistical basis. We recall that the plateau extraction must be
symmetric around $\tau=t/2$ in order to respect the time reversal invariance
relations of section 3.1. The data for $C_{2}$ and $S_{1}$ with $t/a=32$ is
restricted to $\tau/a=16$.

#### Details on the contractions.

We now give details on the calculation of the different lattice contractions
for $C_{4\mathrm{pt}}$. Let us start with a broad overview, which is
pictorially represented in figure 3. Most graphs are calculated using the one-
end trick Foster:1998vw  at the source, which for $A$ and $C_{1}$ we are also
able to use at the sink. $C_{2}$ and $S_{1}$ have a similar structure, for
which the sequential source technique Martinelli:1988rr  is suitable. In
general, loops are obtained using stochastic insertions, except for the
double-insertion loop of $S_{2}$, where usual point-to-all propagators are
used. The $D$ graph is calculated using point sources only. For the evaluation
of the two-point graph $G_{2}$ and the connected part $G_{3}$ of the three-
point graph (see figure 3) we alternatively use point sources or the one-end
trick.

Figure 3:  Lattice methods employed for computing the different graphs of the
four-point function. Colours are used to distinguish the different parts of a
graph. Two stochastic sources within a pion source or sink lead always to the
application of the one-end trick or two-hand-trick, respectively.

We use stochastic $\mathbb{Z}_{2}\otimes i\mathbb{Z}_{2}$ wall sources, i.e. a
set of complex random vectors $\eta_{t}^{(\ell)}$ carrying space-time, spin
and colour indices (not explicitly written here). $\eta_{t}^{(\ell)}$ is
nonzero only on the time slice $t$, where it has components $(\pm 1\pm
i)/\sqrt{2}$ that are random for each gauge configuration. Averaged over all
realisations, one then has

$\displaystyle\frac{1}{N^{\mathrm{st}}_{\mathrm{src}}}\sum_{\ell}^{N^{\mathrm{st}}_{\mathrm{src}}}\eta_{t}^{(\ell)}\otimes\eta_{t}^{\dagger(\ell)}$
$\displaystyle\to\;\mathbb{1}_{t}\,,$ (39)

where the matrix $\mathbb{1}_{t}$ is unity if both time indices are equal to
$t$ and zero otherwise. Using this identity as an expression for the pion
source represents the one-end trick.

The structure of graphs $C_{1}$ and $A$ allows us to perform a second one-end
trick at the sink, which we refer to as the “two-hand trick”. For $C_{1}$ this
trick was applied already in Alexandrou:2008ru . In a different context it was
used in Yang:2015zja , where it was dubbed the “stochastic sandwich method”.
Including the sign of the Wick contraction (easily obtained from the number of
fermion loops) we explicitly have222Note that $C_{1}^{ij,\vec{\mskip
0.5mup}\mskip 1.5mu}(\vec{\mskip 0.5muy}\mskip 1.5mu,\tau,t)$ is a
contribution to the four-point correlator $C_{4\mathrm{pt}}$, whereas
$C_{1}^{ij}(y)$ introduced in section 3.1 is a contribution to the pion matrix
element (36). The same holds for the other contractions.

$\displaystyle C_{1}^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5muy}\mskip 1.5mu,\tau,t)$
$\displaystyle=\frac{a^{3}}{VN^{\mathrm{st}}_{\mathrm{src}}}\sum_{\ell}^{N^{\mathrm{st}}_{\mathrm{src}}}\sum_{\vec{\mskip
0.5mux}\mskip 1.5mu}\biggl{\langle}\left[\Psi_{t}^{\dagger(\ell),\vec{\mskip
0.5mu0}\mskip 1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu+\vec{\mskip 0.5muy}\mskip
1.5mu,\tau)\mskip 1.5mu\gamma_{5}\mskip
1.5mu\Gamma_{i}\Psi_{0}^{(\ell),-\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5mux}\mskip 1.5mu+\vec{\mskip 0.5muy}\mskip 1.5mu,\tau)\right]$
$\displaystyle\quad\times\left[\Psi_{0}^{\dagger(\ell),\vec{\mskip
0.5mu0}\mskip 1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\mskip
1.5mu\gamma_{5}\mskip 1.5mu\Gamma_{j}\Psi_{t}^{(\ell),\vec{\mskip
0.5mup}\mskip 1.5mu}(\vec{\mskip 0.5mux}\mskip
1.5mu,\tau)\right]\biggr{\rangle}\,,$ $\displaystyle A^{ij,\vec{\mskip
0.5mup}\mskip 1.5mu}(\vec{\mskip 0.5muy}\mskip 1.5mu,\tau,t)$
$\displaystyle=-\frac{a^{3}}{V}\sum_{\vec{\mskip 0.5mux}\mskip
1.5mu}\left\langle B_{0}^{i,-\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5mux}\mskip 1.5mu+\vec{\mskip 0.5muy}\mskip 1.5mu,\tau)B_{t}^{j,\vec{\mskip
0.5mup}\mskip 1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\right\rangle$
$\displaystyle\quad+\frac{a^{3}}{V}\sum_{\vec{\mskip 0.5mux}\mskip
1.5mu}\left\langle B_{0}^{i,-\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5mux}\mskip 1.5mu+\vec{\mskip 0.5muy}\mskip
1.5mu,\tau)\right\rangle\left\langle B_{t}^{j,\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\right\rangle\,,$ $\displaystyle
B_{t}^{j,\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip 0.5mux}\mskip
1.5mu,\tau)$
$\displaystyle=\frac{1}{N^{\mathrm{st}}_{\mathrm{src}}}\sum_{\ell}^{N^{\mathrm{st}}_{\mathrm{src}}}\left[\Psi_{t}^{\dagger(\ell),\vec{\mskip
0.5mu0}\mskip 1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\mskip
1.5mu\gamma_{5}\mskip 1.5mu\Gamma_{j}\mskip 1.5mu\Psi_{t}^{(\ell),\vec{\mskip
0.5mup}\mskip 1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\right]\,,$ (40)

where again $V=L^{3}\mskip 1.5mua^{3}$ is the spatial volume. Here and in the
following, the notation $\langle\ldots\rangle$ denotes the average over the
gauge ensemble and $[\ldots]$ indicates a closed spin-colour structure. The
vector $\Psi^{(\ell),\vec{\mskip 0.5mup}\mskip 1.5mu}_{t}$ is obtained from an
inversion of the Wilson-clover Dirac operator $\mathcal{D}$ on the random
source,

$\displaystyle\mathcal{D}\mskip 1.5mu\Psi^{(\ell),\vec{\mskip 0.5mup}\mskip
1.5mu}_{t}$ $\displaystyle=\Phi\mskip 1.5mu\mathcal{E}^{\vec{\mskip
0.5mup}\mskip 1.5mu}\mskip 1.5mu\eta^{(\ell)}_{t}\,,$ (41)

where

$\displaystyle\bigl{(}\,\mathcal{E}^{\vec{\mskip 0.5mup}\mskip
1.5mu}\,\bigr{)}_{\vec{\mskip 0.5mux}\mskip 1.5mu\mskip-1.5mu,\mskip
1.5mu\vec{\mskip 0.5muy}\mskip 1.5mu}$ $\displaystyle=e^{-i\vec{\mskip
0.5mup}\mskip 1.5mu\cdot\vec{\mskip 0.5mux}\mskip 1.5mu}\,\delta_{\vec{\mskip
0.5mux}\mskip 1.5mu\vec{\mskip 0.5muy}\mskip 1.5mu}$ (42)

is a diagonal matrix in the spatial coordinates. For better legibility, the
$(\vec{\mskip 0.5mux}\mskip 1.5mu,t)$-component of $\Psi$ is written as
$\Psi(\vec{\mskip 0.5mux}\mskip 1.5mu,t)$ rather than $\Psi_{\vec{\mskip
0.5mux}\mskip 1.5mut}$ in (3.3). We will use the same notation for other
quantities that are vectors in space-time.

Source smearing is implemented by $\Phi$, which is a hermitian matrix acting
on spatial coordinates and colour indices and consists of 400 Wuppertal
smearing iterations Gusken:1989qx . It turned out that taking $\Phi\mskip
1.5mu\mathcal{E}^{\vec{\mskip 0.5mup}\mskip 1.5mu}$ instead of
$\mathcal{E}^{\vec{\mskip 0.5mup}\mskip 1.5mu}\mskip 1.5mu\Phi$ in (41)
greatly improves the signal for nonzero momenta. This observation contributed
to introducing momentum smearing in Bali:2016lva . The contractions for which
we have data with nonzero $\vec{\mskip 0.5mup}\mskip 1.5mu$ are $C_{1}$ and
$A$ on the lattice with $L=40$. Specifically, $C_{1}$ was computed for all 24
nonzero momenta satisfying $p^{2}\leq 3$ and $A$ for all 6 nonzero momenta
with $p^{2}=1$. Here $p$ is given as a multiple of the smallest non-trivial
lattice momentum $2\pi/(La)$, which is equal to $437\operatorname{MeV}$ for
$L=40$.

As we see in (3.3), the calculation of $A$ involves the subtraction of vacuum
contributions in order to give the fully connected pion matrix element. For
symmetry reasons, these subtractions are only nonzero for the currents $A$ and
$P$.

The connected part of graph $S_{1}$ is obtained by applying the sequential
source method either to stochastic sources with the one-end trick (denoted by
“oet”) or to point sources (denoted by “pt”). The loop appearing in $S_{1}$ is
calculated by using a stochastic source $\eta_{\tau}^{(\ell)}$ at time slice
$\tau$ with the corresponding solution $\chi_{\tau}^{(\ell)}$ of the Dirac
equation,

$\mathcal{D}\chi_{\tau}^{(\ell)}=\eta_{\tau}^{(\ell)}\,.$ (43)

Specifically, we have

$\displaystyle S_{1}^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5muy}\mskip 1.5mu,\tau,t)$ $\displaystyle=-\frac{a^{3}}{V}\sum_{\vec{\mskip
0.5mux}\mskip 1.5mu}\left\langle G_{3}^{i,\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu+\vec{\mskip 0.5muy}\mskip
1.5mu,\tau,t)\,L_{1}^{j}(\vec{\mskip 0.5mux}\mskip
1.5mu,\tau)\right\rangle+\frac{a^{3}}{V}\sum_{\vec{\mskip 0.5mux}\mskip
1.5mu}\left\langle G_{3}^{i,\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5mux}\mskip 1.5mu,\tau,t)\right\rangle\left\langle\\!\left\langle
L_{1}^{j}(\tau)\right\rangle\\!\right\rangle\,,$ $\displaystyle
L_{1}^{j}(\vec{\mskip 0.5muy}\mskip 1.5mu,\tau)$
$\displaystyle=\frac{1}{N^{\mathrm{st}}_{\mathrm{ins}}}\sum_{\ell}^{N^{\mathrm{st}}_{\mathrm{ins}}}\left[\eta_{\tau}^{\dagger(\ell)}(\vec{\mskip
0.5muy}\mskip 1.5mu,\tau)\mskip 1.5mu\Gamma_{j}\mskip
1.5mu\chi_{\tau}^{(\ell)}(\vec{\mskip 0.5muy}\mskip 1.5mu,\tau)\right]$ (44)

and

$\displaystyle\left\langle\\!\left\langle
L_{1}^{j}(\tau)\right\rangle\\!\right\rangle$
$\displaystyle=\frac{a^{3}}{V}\sum_{\vec{\mskip 0.5mux}\mskip
1.5mu}\left\langle L_{1}^{j}(\vec{\mskip 0.5mux}\mskip
1.5mu,\tau)\right\rangle\,.$ (45)

Note that for symmetry reasons, the vacuum subtraction for $S_{1}$ in (3.3) is
only nonzero for the scalar current $S$. If the one-end trick is used, then
$G_{3}^{i,\vec{\mskip 0.5mup}\mskip 1.5mu}$ is given by

$\displaystyle G_{3,\mathrm{oet}}^{i,\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau,t)$
$\displaystyle=\frac{1}{N^{\mathrm{st}}_{\mathrm{src}}}\sum_{\ell}^{N^{\mathrm{st}}_{\mathrm{src}}}\left[{X}_{0\mskip
1.5mut,\mathrm{oet}}^{\dagger(\ell),-\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\mskip 1.5mu\gamma_{5}\mskip
1.5mu\Gamma_{i}\mskip 1.5mu\Psi_{0}^{(\ell),\vec{\mskip 0.5mu0}\mskip
1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\right]$ (46)

with the sequential propagator ${X}_{0\mskip 1.5mut}^{(\ell),\vec{\mskip
0.5mup}\mskip 1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)$ obtained by
inversion of

$\displaystyle\Bigl{(}\mathcal{D}{X}_{0\mskip
1.5mut,\mathrm{oet}}^{(\ell),\vec{\mskip 0.5mup}\mskip
1.5mu}\Bigr{)}(\vec{\mskip 0.5mux}\mskip 1.5mu^{\prime},t^{\prime})$
$\displaystyle=\delta_{tt^{\prime}}\Bigl{(}\Phi\,\mathcal{E}^{\vec{\mskip
0.5mup}\mskip 1.5mu}\,\Phi\mskip 1.5mu\gamma_{5}\mskip
1.5mu\Psi_{0}^{(\ell),-\vec{\mskip 0.5mup}\mskip 1.5mu}\Bigr{)}(\vec{\mskip
0.5mux}\mskip 1.5mu^{\prime},t^{\prime})\,.$ (47)

We refrain from writing out the corresponding expressions for
$G_{3,\mathrm{pt}}$ and $\mathrm{X}_{0t,\mathrm{pt}}$ with point sources,
given that they are identical to those used in standard computations of the
pion form factor, see e.g. Brommel:2006ww . We find very good agreement
between $S_{1}$ computed with the one-end trick and with point sources, with
slightly larger statistical errors for the latter. In later sections, only the
one-end trick results are used. By contrast, we take the point-source version
of $G_{3}$ to evaluate the connected contributions to the vector and scalar
form factors.

The contraction $C_{2}$ has a similar structure as the connected part of
$S_{1}$, but it requires the calculation of an additional propagator between
the two current insertions. For its evaluation we again use a stochastic
source at the insertion time slice $\tau$. To reduce statistical noise, we
make the following improvements. We consider the hopping parameter expansion
of the propagator Thron:1997iy ; Bali:2005fu ; Bali:2009hu , writing
$\mathcal{D}=(\mathbb{1}-H)/(2\kappa)$ with the hopping term $H$, and make use
of the geometric series:

$\mathcal{M}=\mathcal{D}^{-1}=2\kappa\left(\mathbb{1}-H\right)^{-1}=2\kappa\sum_{n=0}^{\infty}H^{n}=2\kappa\\!\\!\sum_{n=0}^{n(\vec{\mskip
0.5muy}\mskip 1.5mu)-1}H^{n}+2\kappa\\!\\!\sum_{n=n(\vec{\mskip 0.5muy}\mskip
1.5mu)}^{\infty}H^{n}\,.$ (48)

Since for the Wilson-clover action, $H$ involves at most nearest neighbours on
the lattice, one needs at least

$n(\vec{\mskip 0.5muy}\mskip
1.5mu)=\sum_{i=1}^{3}\min\biggl{(}\frac{|y_{i}|}{a},L-\frac{|y_{i}|}{a}\biggr{)}$
(49)

hopping terms to obtain a non-vanishing contribution to the propagator from a
point $\vec{\mskip 0.5mux}\mskip 1.5mu$ to $\vec{\mskip 0.5mux}\mskip
1.5mu+\vec{\mskip 0.5muy}\mskip 1.5mu$ on a periodic lattice of spatial size
$La$. Hence, the first sum on the r.h.s. of (48) can be omitted, and we get

$\displaystyle\mathcal{M}$ $\displaystyle=2\kappa\\!\\!\sum_{n=n(\vec{\mskip
0.5muy}\mskip 1.5mu)}^{\infty}H^{n}=H^{n(\vec{\mskip 0.5muy}\mskip
1.5mu)}\;2\kappa\sum_{n=0}^{\infty}H^{n}=H^{n(\vec{\mskip 0.5muy}\mskip
1.5mu)}\mathcal{M}$ for propagation from $\vec{\mskip 0.5mux}\mskip 1.5mu$ to
$\vec{\mskip 0.5mux}\mskip 1.5mu+\vec{\mskip 0.5muy}\mskip 1.5mu$, (50)

where in the last step we used (48) again. Taking $H^{n(\vec{\mskip
0.5muy}\mskip 1.5mu)}\mathcal{M}$ instead of $\mathcal{M}$ itself, we
implicitly omit those terms that do not contribute to the propagation but may
add to the stochastic noise. The expression to be evaluated for the $C_{2}$
graph is then given by

$\displaystyle C_{2}^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5muy}\mskip 1.5mu,\tau,t)$
$\displaystyle=\frac{a^{3}}{VN^{\mathrm{st}}_{\mathrm{src}}}\sum_{\ell}^{N^{\mathrm{st}}_{\mathrm{src}}}\sum_{\vec{\mskip
0.5mux}\mskip 1.5mu}\biggl{\langle}\left[{X}_{0\mskip
1.5mut,\mathrm{oet}}^{\dagger(\ell),-\vec{\mskip 0.5mup}\mskip
1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\mskip 1.5mu\gamma_{5}\mskip
1.5mu\Gamma_{j}\,\xi_{\tau}^{(\ell),n(\vec{\mskip 0.5muy}\mskip
1.5mu)}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\right]$
$\displaystyle\quad\times\left[\eta_{\tau}^{\dagger(\ell)}(\vec{\mskip
0.5mux}\mskip 1.5mu+\vec{\mskip 0.5muy}\mskip 1.5mu,\tau)\mskip
1.5mu\Gamma_{i}\mskip 1.5mu\Psi_{0}^{(\ell),\vec{\mskip 0.5mu0}\mskip
1.5mu}(\vec{\mskip 0.5mux}\mskip 1.5mu+\vec{\mskip 0.5muy}\mskip
1.5mu,\tau)\right]\biggr{\rangle}$ (51)

with

$\xi_{\tau}^{(\ell),n}=H^{n}\chi_{\tau}^{(\ell)}=H^{n}\mskip
1.5mu\mathcal{D}^{-1}\eta^{(\ell)}_{\tau}\,.$ (52)

For $D$ and $S_{2}$, one needs the pion two-point graph $G_{2}$ and two single
insertion loops $L_{1}$ or one double insertion loop $L_{2}$, respectively.
For $G_{2}$ we again use either one-end trick or point sources. The double
insertion loop is evaluated using a point source $S_{x}$ at $x=(\vec{\mskip
0.5mu0}\mskip 1.5mu,\tau)$, from which the point-to-all propagator
$\mathcal{M}_{x}$ is obtained by inverting

$\displaystyle\mathcal{D}\mathcal{M}_{x}$ $\displaystyle=S_{x}\,.$ (53)

Note that $S_{x}$ and hence also $\mathcal{M}_{x}$ is a vector in space-time
but a matrix in spin and colour. The space-time components of the point source
are $S_{x}(x^{\prime})=\mathbb{1}\mskip 1.5mu\delta_{xx^{\prime}}$. We then
have

$\displaystyle L_{2}^{ij}(\tau,\vec{\mskip 0.5muy}\mskip 1.5mu)$
$\displaystyle=\operatorname{tr}\left\\{\gamma_{5}\mskip
1.5mu\mathcal{M}^{\dagger}_{\vec{\mskip 0.5mu0}\mskip 1.5mu\tau}(\vec{\mskip
0.5muy}\mskip 1.5mu,\tau)\mskip 1.5mu\gamma_{5}\mskip
1.5mu\Gamma_{i}\,\mathcal{M}_{\vec{\mskip 0.5mu0}\mskip 1.5mu\tau}(\vec{\mskip
0.5muy}\mskip 1.5mu,\tau)\,\Gamma_{j}\right\\}$ (54)

with the trace referring to spin and colour indices, and

$\displaystyle S_{2}^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5muy}\mskip 1.5mu,\tau,t)$ $\displaystyle=-\left\langle
G_{2,\mathrm{oet}}^{\vec{\mskip 0.5mup}\mskip
1.5mu}(t)\,L_{2}^{ij}(\vec{\mskip 0.5muy}\mskip
1.5mu,\tau)\right\rangle+\left\langle G_{2,\mathrm{oet}}^{\vec{\mskip
0.5mup}\mskip 1.5mu}(t)\vphantom{L_{2}^{ij}}\right\rangle\left\langle
L_{2}^{ij}(\vec{\mskip 0.5muy}\mskip 1.5mu,\tau)\right\rangle\,,$
$\displaystyle D^{ij,\vec{\mskip 0.5mup}\mskip 1.5mu}(\vec{\mskip
0.5muy}\mskip 1.5mu,\tau,t)$ $\displaystyle=\frac{a^{3}}{V}\sum_{\vec{\mskip
0.5mux}\mskip 1.5mu}\,\biggl{\\{}\left\langle G_{2,\mathrm{pt}}^{\vec{\mskip
0.5mup}\mskip 1.5mu}(t)\,L_{1}^{i}(\vec{\mskip 0.5mux}\mskip 1.5mu+\vec{\mskip
0.5muy}\mskip 1.5mu,\tau)\,L_{1}^{j}(\vec{\mskip 0.5mux}\mskip
1.5mu,\tau)\right\rangle$ $\displaystyle\qquad-\left\langle
G_{2,\mathrm{pt}}^{\vec{\mskip 0.5mup}\mskip
1.5mu}(t)\vphantom{L_{2}^{ij}}\right\rangle\left\langle L_{1}^{i}(\vec{\mskip
0.5mux}\mskip 1.5mu+\vec{\mskip 0.5muy}\mskip
1.5mu,\tau)\,L_{1}^{j}(\vec{\mskip 0.5mux}\mskip 1.5mu,\tau)\right\rangle$
$\displaystyle\qquad-\left\langle G_{2,\mathrm{pt}}^{\vec{\mskip 0.5mup}\mskip
1.5mu}(t)\,L_{1}^{j}(\vec{\mskip 0.5mux}\mskip
1.5mu,\tau)\right\rangle\left\langle\\!\left\langle
L_{1}^{i\phantom{j\\!\\!}}(\tau)\right\rangle\\!\right\rangle-\left\langle
G_{2,\mathrm{pt}}^{\vec{\mskip 0.5mup}\mskip
1.5mu}(t)\,L_{1}^{i\phantom{j\\!\\!}}(\vec{\mskip 0.5mux}\mskip
1.5mu,\tau)\right\rangle\left\langle\\!\left\langle
L_{1}^{j}(\tau)\right\rangle\\!\right\rangle\biggr{\\}}$
$\displaystyle\quad+2\left\langle G_{2,\mathrm{pt}}^{\vec{\mskip 0.5mup}\mskip
1.5mu}(t)\vphantom{L_{2}^{ij}}\right\rangle\left\langle\\!\left\langle
L_{1}^{i\phantom{j\\!\\!}}(\tau)\right\rangle\\!\right\rangle\left\langle\\!\left\langle
L_{1}^{j}(\tau)\right\rangle\\!\right\rangle$ (55)

with $\langle\langle L_{1}\rangle\rangle$ defined in (45). When evaluated with
the one-end trick, the two point graph reads

$\displaystyle G_{2,\mathrm{oet}}^{\vec{\mskip 0.5mup}\mskip 1.5mu}(t)$
$\displaystyle=\frac{1}{N^{\mathrm{st}}_{\mathrm{src}}}\sum_{\ell}^{N^{\mathrm{st}}_{\mathrm{src}}}\biggl{[}\Psi_{0}^{\dagger(\ell),\vec{\mskip
0.5mu0}\mskip 1.5mu}\mskip 1.5mu\Phi\mskip 1.5mu\mathcal{E}^{\vec{\mskip
0.5mup}\mskip 1.5mu}\mskip 1.5mu\Phi\mskip 1.5mu\Psi_{0}^{(\ell),-\vec{\mskip
0.5mup}\mskip 1.5mu}\biggr{]}\,.$ (56)

Note that $\langle L_{2}^{ij}\rangle$ and $\langle
L_{1}^{i\phantom{j\\!\\!}}\,L_{1}^{j}\mskip 1.5mu\rangle$ depend on the
distance $\vec{\mskip 0.5muy}\mskip 1.5mu$ between the currents and can be
nonzero for most combinations of Dirac matrices $\Gamma_{i}$ and $\Gamma_{j}$.
This makes vacuum subtractions necessary for $S_{2}$ and $D$ in most channels.

As indicated in (3.3), we use only stochastic sources for $S_{2}$ and only
point sources for $D$. For the pion two-point function
$C_{2\mathrm{pt}}(t)=\langle G_{2}(t)\rangle$, we compare both methods and
find excellent agreement.

| $L=40$ | $L=32$
---|---|---
graph | $N_{\mathrm{src}}^{\mathrm{st}}$ | $N_{\mathrm{ins},i}^{\mathrm{st}}$ | $N_{\mathrm{ins},j}^{\mathrm{st}}$ | $N_{\mathrm{src}}$ | $N_{\mathrm{ins}}$ | $N_{\mathrm{src}}^{\mathrm{st}}$ | $N_{\mathrm{ins},i}^{\mathrm{st}}$ | $N_{\mathrm{ins},j}^{\mathrm{st}}$ | $N_{\mathrm{src}}$ | $N_{\mathrm{ins}}$
$C_{1}$ | 1 | - | - | 1 | $L^{3}$ | 1 | - | - | 1 | $L^{3}$
$C_{2}$ | 1 | 10 | - | 5 | $L^{3}$ | 1 | 20 | - | 1 | $L^{3}$
$A$ | 1 | - | - | 4 | $L^{3}$ | 1 | - | - | 1 | $L^{3}$
$S_{1}$ | - | - | 120 | 4 | $L^{3}$ | - | - | 120 | 3 | $L^{3}$
$S_{2}$ | 1 | - | - | 16 | 16 | 1 | - | - | 1 | 2
$D$ | - | 60 | 60 | 1 | $L^{3}$ | - | 60 | 60 | 1 | $L^{3}$
Table 2:  Numbers $N^{\text{st}}$ of stochastic noise vectors and numbers of
sources and current insertions used for each graph in our simulations for the
two lattice volumes. For most graphs, stochastic propagators are connected to
the insertion, which implies an average over the entire spatial lattice
volume. This is indicated by $N_{\mathrm{ins}}=L^{3}$.

To increase statistics, we repeat the calculations at $N_{\mathrm{src}}$ time
positions $t_{s}$ of the source (keeping $\tau$ and $t$ fixed relative to
$t_{s}$) and for several spatial insertion positions $N_{\mathrm{ins}}$ of the
current. For most graphs, the latter is implicitly realised by using
stochastic wall sources. The corresponding numbers, as well as the size of the
stochastic noise vector sets, are given in table 2. Note that in some cases
these numbers differ between the lattices with $L=40$ and $L=32$ (ensembles V
and IV of table 1).

###  3.4 Renormalisation

We convert all lattice currents, which are defined in the lattice scheme at a
lattice spacing $a$, into the $\overline{\mathrm{MS}}$ scheme at the
renormalisation scale $\mu=2\operatorname{GeV}$. The corresponding
renormalisation constants depend on the gauge coupling $g^{2}=6/\beta$, where
in our case $\beta=5.29$. For currents with nonzero anomalous dimension, such
as in the scalar and pseudoscalar cases, there is also a dependence on $a\mu$.
The conversion factors used for the correlator (36) of currents $i$ and $j$
read

$\displaystyle R_{ij}^{\scalebox{0.6}{$\overline{\mathrm{MS}}$}}$
$\displaystyle=\tilde{Z}_{i}\mskip 1.5mu\tilde{Z}_{j}\mskip
1.5muR_{ij}^{\mathrm{lat}}$ $\displaystyle\text{ \leavevmode\nobreak\
with\leavevmode\nobreak\ }\tilde{Z}_{i}$
$\displaystyle=Z^{\scalebox{0.6}{$\overline{\mathrm{MS}}$}}_{i}\left(1+am_{q}\mskip
1.5mub_{i}\right)\,,$ (57)

see e.g. Bali:2014nma . Here we have included the coefficients $b_{i}$ for the
mass-dependent order $a$ improvement. They become particularly relevant at the
charm quark mass. In the vector and axial vector cases, there are additional
mass-independent order $a$ improvement terms (accompanied by $c_{V}$ and
$c_{A}$, respectively), which we ignore in the present work. Note that some of
the matrix elements listed in section 3.2 receive contributions from flavour
singlet current combinations. In these cases our renormalisation procedure
should be regarded as approximate, because we do not take into account
operator mixing.

The renormalisation factors $Z_{i}^{\scalebox{0.6}{$\overline{\mathrm{MS}}$}}$
have been obtained in Gockeler:2010yr  (updated in Bali:2014nma ) in a two-
step procedure: first the currents were matched non-perturbatively from the
lattice scheme to the RI’/MOM scheme and then from there to the
$\overline{\mathrm{MS}}$ scheme in continuum perturbation theory. The
coefficients $b_{i}=1+\mathcal{O}(g^{2})$ have been computed in lattice
perturbation theory Capitani:2000xi ; Sint:1997jx ; Taniguchi:1998pf .
Furthermore, in Fritzsch:2010aw  the coefficient $b_{S}$ was determined non-
perturbatively as

$b_{S}^{\mathrm{np}}=\bigl{(}1+0.19246\mskip
1.5mug^{2}\bigr{)}\,\frac{1-0.3737\mskip 1.5mug^{10}}{1-0.5181\mskip
1.5mug^{4}}$ (58)

with an uncertainty of about 5%. We can determine $b_{V}$ non-perturbatively
by evaluating the vector form factor at zero momentum transfer with charm
quarks ($am_{q}=0.3147$). Multiplied by
$Z_{V}^{\scalebox{0.6}{$\overline{\mathrm{MS}}$}}(1+am_{q}\mskip 1.5mub_{V})$
this must give $1$, owing to charge conservation. Using the result
$F_{V}(0)=0.9068(2)$ from our $L=40$ lattice and the value of
$Z_{V}^{\scalebox{0.6}{$\overline{\mathrm{MS}}$}}$ from Bali:2014nma , we
extract

$\displaystyle b_{V}^{\mathrm{np}}$ $\displaystyle=1.586(32)\,,$ (59)

where the uncertainty is dominated by the uncertainty on
$Z_{V}^{\scalebox{0.6}{$\overline{\mathrm{MS}}$}}$.

| $S$ | $P$ | $V$ | $A$ | $T$
---|---|---|---|---|---
$Z^{\scalebox{0.6}{$\overline{\mathrm{MS}}$}}$ | 0.6153(25) | 0.476(13) | 0.7356(48) | 0.76487(64) | 0.8530(25)
$b^{\text{pert}}$ | 1.3453 | 1.2747 | 1.2750 | 1.2731 | 1.2497
$b^{\text{np}}$ | 1.091(55) |  | 1.586(32) |  |
$b^{\text{resc}}$ | 1.673 | 1.586 | 1.586 | 1.584 | 1.555
Table 3:  Renormalisation factors
$Z_{i}^{\scalebox{0.6}{$\overline{\mathrm{MS}}$}}$ from Bali:2014nma  for the
different currents (3) in the $\overline{\mathrm{MS}}$ scheme at scale
$\mu=2\operatorname{GeV}$. We also list the perturbative estimates of the
coefficients $b_{i}$ from Bali:2014nma , the non-perturbative determinations
(58) and (59), and the rescaled values (60). For the latter, we estimate an
error of 34% as explained in the text.

Table 3 gives the values for
$Z_{i}^{\scalebox{0.6}{$\overline{\mathrm{MS}}$}}(2\operatorname{GeV})$, as
well as estimates $b_{i}^{\text{pert}}$ using one-loop perturbation theory
with an “improved” coupling constant (for details see (26) and (27) in
Bali:2014nma ). In the next row, we give the non-perturbative estimates
$b_{i}^{\text{np}}$ from (58) and (59). We see that the non-perturbative value
for $b_{V}$ is 24% larger than the perturbative estimate, while $b_{S}$ from
Fritzsch:2010aw  is about 19% smaller than the perturbative result. We remark
that different non-perturbative determinations of $b_{i}$ only need to agree
up to order $a$ effects.

Based on our result for $b_{V}$, we make the naive assumption that all non-
perturbative coefficients are larger than the perturbative estimates and
define rescaled coefficients

$b^{\text{resc}}_{i}=b^{\text{pert}}_{i}\,\frac{b^{\text{np}}_{V}}{b^{\text{pert}}_{V}}\,,$
(60)

which are also listed in the table. As is clear from the scalar channel, there
is a huge uncertainty in this procedure. We find
$|1-b_{S}^{\text{np}}/b_{S}^{\text{resc}}|=34\%$ and take this as our
uncertainty of the coefficients $b_{i}^{\text{resc}}$ with $i\neq V$. The
results in the following chapters are obtained with the coefficients
$b_{i}^{\mathrm{np}}$ for $i=S,V$ and with $b_{i}^{\mathrm{resc}}$ for all
other currents.

##  4 Data quality and lattice artefacts

In this section we investigate the quality of our data and discuss a number of
lattice artefacts. We concentrate on the data for light quarks here, for which
we have the largest set of simulations. We verified that for heavy quarks, the
lattice artefacts are not worse than for light ones.

###  4.1 Plateaux and excited state contributions

For the contractions $C_{1}$ and $A$, we have data for the ratio (36) of four-
point and two-point functions in a wide range of the current insertion time
$\tau$ and can hence verify the existence of a plateau. Figure 4 shows this
ratio for two cases that are representative of channels with a good signal. In
the case of figure 4(a) we have data for sink-source separations of $t/a=15$
and $32$, which we plot in such a way that the midpoints $\tau=t/2$ coincide
in both cases. Throughout this paper, error bars are purely statistical and
determined by 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 (see table 1).

The time separation $t/a=32$ is half the temporal extent of our lattices, so
that the pion two-point function $C_{2\mathrm{pt}}(t)$ at that point has equal
contributions from pions that propagate forward or backward from the source at
time zero. To extract the pion matrix element, we must hence multiply
$C_{4\mathrm{pt}}(\vec{\mskip 0.5muy}\mskip 1.5mu,\tau,t)/C_{2\mathrm{pt}}(t)$
in (36) with a normalisation factor $N(32\mskip 1.5mua)=2$, which is also done
in the figure. For $t/a=15$, the contribution from backward propagating pions
in the two-point function is below $3\%$. We do not correct for this
contamination and take $N(15\mskip 1.5mua)=1$ in this case.

(a) $C_{1},\langle A^{0}A^{0}\rangle,\vec{\mskip 0.5muy}\mskip
1.5mu/a=(0,3,1),p=0,L=40$ (b) $A,\langle SS\rangle,\vec{\mskip 0.5muy}\mskip
1.5mu/a=(0,3,1),p=0,L=40,t=15\mskip 1.5mua$

Figure 4:  Dependence of the four-point function on the current insertion time
$\tau$ at a fixed distance $\vec{\mskip 0.5muy}\mskip 1.5mu$ between the two
currents. Shown are the contributions of the contraction $C_{1}$ to the
correlator $\langle A^{0}A^{0}\rangle$ and of the contraction $A$ to $\langle
SS\rangle$. Here and in all following plots, errors are purely statistical.
The normalisation factor $N(t)$ is specified in the text. $p$ denotes the pion
momentum and $L$ the spatial lattice extent. Unless mentioned otherwise, all
data shown in this paper are for light quarks.

The data in figure 4 show a clear plateau around $\tau=t/2$. In figure 4(a)
the plateaux for the two source-sink separations are in good agreement with
each other, which indicates that contributions from excited states are not
very large. Based on the plots and their analogues for other channels, we
extract plateau values by fitting to a constant in the ranges

$\displaystyle 6\leq\tau/a\leq 9$ $\displaystyle\text{ \leavevmode\nobreak\
for\leavevmode\nobreak\ }t/a=15\,,$ $\displaystyle 14\leq\tau/a\leq 18\text{
\leavevmode\nobreak\ and\leavevmode\nobreak\ }46\leq\tau/a\leq 50$
$\displaystyle\text{ \leavevmode\nobreak\ for\leavevmode\nobreak\ }t/a=32\,.$
(61)

For $t/a=32$ we thus combine the two plateaux around $\tau/a=16$ and $48$.

Examples for the extracted two-current matrix elements as a function of the
distance between the currents are shown in figure 5. Here and in the
following, we write $y=|\vec{\mskip 0.5muy}\mskip 1.5mu|$. We see that the
data for $t/a=15$ and $32$ are consistent within their statistical
uncertainties. In these plots and in their analogues for the other channels,
we do not find any indication for large excited state contamination.

(a) $C_{1},\langle V^{0}V^{0}\rangle,p=0,L=40$ (b) $C_{1},\langle
PP\rangle,p=0,L=40$ (c) $C_{2},\langle SS\rangle,p=0,L=32$ (d) $S_{1},\langle
PP\rangle,p=0,L=32$

Figure 5:  Comparison between data with $t/a=15$ and $32$ for different
contractions and different lattice sizes. Shown are only data points
satisfying the cut (62) discussed in the next section. The correlation
functions in this and all following figures have been converted to the
$\overline{\mathrm{MS}}$ scheme at the scale $\mu=2\operatorname{GeV}$.

###  4.2 Anisotropy effects

In the continuum, the correlation functions $\langle V^{0}V^{0}\rangle,\langle
A^{0}A^{0}\rangle,\langle SS\rangle$ and $\langle PP\rangle$ in a pion at rest
can only depend on the modulus $y=|\vec{\mskip 0.5muy}\mskip 1.5mu|$. Lattice
regularisation breaks rotational invariance, so that an anisotropy in the
$\vec{\mskip 0.5muy}\mskip 1.5mu$ dependence of our correlators is a clear
indication of lattice artefacts.

At large distances $y/a\sim L/2$ we see clear signs of anisotropy, as shown in
figure 6(a) to (c). As explained in detail in Burkardt:1994pw , this can be
traced back to the use of a finite lattice with periodic boundary conditions:
the two currents are sensitive to the “images” of the pion in lattice cells
adjacent to the elementary cell of size $L^{3}$. At a given distance $y$, this
effect is smallest for vectors $\vec{\mskip 0.5muy}\mskip 1.5mu$ with
directions close to space diagonals $(\pm 1,\pm 1,\pm 1)$ and largest if
$\vec{\mskip 0.5muy}\mskip 1.5mu$ points along one of the coordinate axes.
This gives rise to the distinct “sawtooth” pattern seen in our plots, which
was also observed in earlier studies Burkardt:1994pw ; Alexandrou:2008ru . An
analytic investigation of the same phenomenon using a low-energy effective
field theory can be found in Briceno:2018lfj . A method to correct for these
periodic images was proposed in Burkardt:1994pw  and employed in
Alexandrou:2008ru , where simulations were performed on $L=24$ lattices with
$a\approx 0.077\operatorname{fm}$. The resulting physical length $La$ was thus
0.81 times the length of the smallest lattice used in our study.

(a) $\langle V^{0}V^{0}\rangle,p=0,L=40$ (b) $\langle SS\rangle,p=0,L=40$ (c)
$\langle PP\rangle,p=0,L=40$ (d) $\langle PP\rangle,p=0,L=40$

Figure 6:  Anisotropy effects in our data at large and intermediate $y$ (light
points) and their removal by the cut (62) (dark points).

Benefiting from the larger size of our lattices, we take a less sophisticated
approach here and consider only data for distances $\vec{\mskip 0.5muy}\mskip
1.5mu$ close to the space diagonals. To this end, we use exact lattice
symmetries to map the vector $\vec{\mskip 0.5muy}\mskip 1.5mu$ into the octant
where all its components are non-negative and define $\vartheta(\vec{\mskip
0.5muy}\mskip 1.5mu)$ as the angle between that vector and the diagonal
$(1,1,1)$. We then retain only data points satisfying

$\displaystyle\cos\vartheta(\vec{\mskip 0.5muy}\mskip 1.5mu)$
$\displaystyle\geq 0.9\,.$ (62)

After this cut, we statistically average the correlators over points
$\vec{\mskip 0.5muy}\mskip 1.5mu$ of the same length $y$, which greatly
decreases statistical errors. The points thus combined are not necessarily
related by an exact lattice symmetry operation, and we checked that
restricting the data combination to points that are equivalent on the lattice
gives results consistent with those of the more inclusive combination.

As seen in figure 6(a) to (c), the above procedure efficiently removes
anisotropy effects while retaining enough data points. It also removes a broad
anisotropic structure that we observe for $\langle PP\rangle$ in the $S_{2}$
graph, figure 6(d), for which we have no interpretation.

(a) $\langle SS\rangle,p=0,L=40$ (b) $\langle PP\rangle,p=0,L=40$

Figure 7:  Anisotropy effects for the contraction $C_{2}$ at small $y$ (light
points) and their removal by the cut (62) (dark points).

A different type of anisotropy effect appears at small $y$ in our data for the
contraction $C_{2}$: in channels with small statistical errors, we observe a
huge dependence of the correlators on the direction of $\vec{\mskip
0.5muy}\mskip 1.5mu$, as shown in figure 7. To understand this, we note that
in $C_{2}$ the two current insertions are directly connected by a fermion
propagator, which is hence evaluated for a fixed space-time direction. The
effect we see may thus directly reflect a discretisation effect, namely the
anisotropy of the fermion propagator on the lattice. This anisotropy was
studied in Cichy:2012is  (see figure 1 in that work), where it was shown that
lattice artefacts are smallest for distances $y$ close to the diagonal
$n=(1,1,1,1)$ in Euclidean space-time. A cut on $yn$ advocated in Cichy:2012is
is equivalent to our cut on $\cos\vartheta(\vec{\mskip 0.5muy}\mskip 1.5mu)$,
given that we always have $y^{0}=0$. As seen in our figure 7, the cut (62)
indeed removes the anisotropy effects in $C_{2}$ and yields points that can be
connected by a reasonably smooth curve. One may expect a similar situation for
the contraction $S_{2}$, where the two currents are connected by two lattice
propagators, but our statistical errors for $S_{2}$ are too large to
unambiguously identify this effect.

To summarise, we apply the cut (62) and the averaging procedure described
above to all lattice data shown in the rest of this work. One must of course
bear in mind that the filtered data may still be affected by discretisation
effects at small $y$ and by finite size effects at large $y$. Regarding the
former, one should be most careful in the interpretation of data with $y/a$
below, say, $2$ or $3$. For the latter, we have an additional handle from the
comparison of our two lattice volumes, which we discuss next.

###  4.3 Finite volume effects

With lattices of two spatial volumes and otherwise identical parameters, we
can investigate finite volume effects on our observables. Figure 8 compares
the results of the two lattices with $L=32$ and $L=40$ for a representative
selection of contractions and currents. Here and in the rest of this paper, we
show correlators on a linear scale, except in cases where their values cover
such a wide range that a logarithmic scale is more advantageous.

As is to be expected, finite size effects are more pronounced at larger
distances $y$ between the two currents. However, we see that, depending on the
contraction and on the operator, volume effects may persist down to small $y$.
Clearly, one needs to keep these in mind when comparing our lattice data with
predictions e.g. from chiral perturbation theory in an infinite volume. We
will come back to this aspect in section 6.1.

(a) $|\mskip 1.5muC_{1}|,\langle V^{0}V^{0}\rangle,p=0$ (b) $A,\langle
A^{0}A^{0}\rangle,p=0$ (c) $C_{2},\langle SS\rangle,p=0$ (d) $A,\langle
SS\rangle,p=0$ (e) $C_{1},\langle PP\rangle,p=0$ (f) $S_{1},\langle
PP\rangle,p=0$

Figure 8:  Comparison of correlators for the two spatial volumes used in our
simulations.

We notice in figure 8 that for the contractions $A$ and $S_{1}$, the
statistical errors are significantly larger for the $L=32$ lattice. For $A$
this is plausible from the number of source time insertions, which is
$N_{\text{src}}=4$ for $L=40$ but only $N_{\text{src}}=1$ for $L=32$ (see
table 2). For $S_{1}$ the difference in $N_{\text{src}}$ between the two
lattices is less pronounced, and we conclude that statistical fluctuations for
this disconnected graph are particularly sensitive to the simulation volume.

###  4.4 Finite momentum and boost invariance

The correlation functions we are studying in this work are matrix elements of
two currents at equal time in a pion at rest. As discussed after equation (4),
we can test that our lattice computation yields boost invariant results using
the correlation functions $\langle SS\rangle$ and $\langle PP\rangle$. These
are Lorentz invariant and hence depend only on the invariant scalar products
$y^{2}$ and $py$. Taking always $y^{0}=0$, this means that the correlation
functions in these channels must coincide for equal $\vec{\mskip 0.5muy}\mskip
1.5mu^{2}$ if we select distance vectors with $\vec{\mskip 0.5mup}\mskip
1.5mu\cdot\vec{\mskip 0.5muy}\mskip 1.5mu=0$. This comparison is shown for the
contraction $C_{1}$ in figure 10 with momenta satisfying $\vec{\mskip
0.5mup}\mskip 1.5mu^{2}=0,1,2$, where $\vec{\mskip 0.5mup}\mskip 1.5mu$ is
given in multiples of $2\pi/(La)\approx 437\operatorname{MeV}$. We see a good
signal for the higher momenta and full agreement between data for different
momenta within statistical uncertainties. This agreement persists for
$\vec{\mskip 0.5mup}\mskip 1.5mu^{2}=3$, which is not shown in the figure
because the larger errors would obscure the clarity of the plot.

(a) $C_{1},\langle SS\rangle,L=40$ (b) $C_{1},\langle PP\rangle,L=40$

Figure 9:  Comparison of the correlator $C_{1}$ for $\langle SS\rangle$ and
$\langle PP\rangle$ with different pion momenta. We select distances
satisfying $\vec{\mskip 0.5mup}\mskip 1.5mu\cdot\vec{\mskip 0.5muy}\mskip
1.5mu=0$ and combine lattice data for all $\vec{\mskip 0.5mup}\mskip 1.5mu$
with the same modulus $p=|\vec{\mskip 0.5mup}\mskip 1.5mu|$, which is given in
multiples of $2\pi/(La)\approx 437\operatorname{MeV}$.

(a) $F_{2},\langle SS\rangle,L=40$ (b) $F_{2},\langle PP\rangle,L=40$

Figure 10:  As figure 10, but for the isospin amplitude $F_{2}=2(C_{1}+A)$.

We also have finite momentum data for the contraction $A$, restricted in this
case to $\vec{\mskip 0.5mup}\mskip 1.5mu^{2}=1$. We compare this with zero
momentum data in figure 10 for the combination $F_{2}=2(C_{1}+A)$ of
contractions, which according to (3.2) represents a physical matrix element in
a definite isospin channel. Again, the data for different momenta are in good
agreement with each other.

##  5 Results for individual lattice contractions

In this section we present results for the individual lattice contractions
shown in figure 2. Having compared zero and finite pion momenta in section
4.4, we restrict our attention to zero momentum data from now on.

###  5.1 Light quarks

(a) $\langle V^{0}V^{0}\rangle,p=0,L=40$ (b) $\langle
A^{0}A^{0}\rangle,p=0,L=40$ (c) $\langle SS\rangle,p=0,L=40$ (d) $\langle
PP\rangle,p=0,L=40$

Figure 11:  Correlators of the four current combinations considered in this
work. The contributions of all graphs are shown separately, except for $D$,
which has huge statistical errors.

An overview of the different contractions for each current combination is
given in figure 11. We show all contractions except for the doubly
disconnected graph $D$, where we find huge errors in all channels. Not only do
these errors prevent us from seeing a signal for $D$, but in addition they are
large compared with the signal of all other contractions, so that even the
upper bounds on $D$ that we could extract from our data would not be very
useful. The large statistical fluctuations of $D$ can be traced back to the
vacuum subtraction terms in (3.3). For all currents, graph $D$ has a fully
connected contribution $\langle G_{2}\,L_{1}^{i\phantom{j\\!\\!}}\mskip
1.5muL_{1}^{j}\mskip 1.5mu\rangle$ and a vacuum subtraction term $\langle
G_{2}\rangle\,\langle L_{1}^{i\phantom{j\\!\\!}}\mskip 1.5muL_{1}^{j}\mskip
1.5mu\rangle$ (plus additional subtractions involving $\langle\langle
L_{1}\rangle\rangle$ in the case of $\langle SS\rangle$). Whilst we see
relatively good signals for the separate terms in some cases, at least for
smaller distances $y$, the connected and subtraction terms turn out to be
almost equal in size and opposite in sign. Adding them, we lose any signal in
the statistical noise. We will discuss in section 6.1 how to deal with this
situation for physical matrix elements that receive contributions from $D$.

Let us return to figure 11. Given the strong decrease of most correlators with
$y$, we start the plots at $y=6a\approx 0.43\operatorname{fm}$ to give a
clearer picture of the situation at larger distances. Compared with what is
seen in the figure, no qualitative changes occur at smaller $y$.

We observe a marked difference between the different currents as regards the
relative importance of different contractions. For the vector current, the
connected graphs $C_{1}$ and $C_{2}$ are clearly dominant, whilst the
disconnected contributions $S_{1}$, $S_{2}$ and the annihilation graph $A$ are
much smaller and in fact consistent with zero. In the axial current case,
$C_{1}$ is clearly dominating over all other contractions at smaller
distances; for $y$ above $1\operatorname{fm}$ the errors on $C_{2}$ and
$S_{2}$ become too large to draw strong conclusions. For both $\langle
SS\rangle$ and $\langle PP\rangle$ we obtain a very clear signal for all
contractions over a wide $y$ range and see that apart from $C_{1}$ and $C_{2}$
the annihilation graph $A$ is of appreciable size, along with $S_{1}$ and
$S_{2}$ in the case of $\langle PP\rangle$.

We note that the graph $C_{1}$ is only dominant in the case of $\langle
A^{0}A^{0}\rangle$ and at relatively small distances. This is among the most
striking results of our study, given that one readily associates this graph
with the “valence content” of a pion and might have expected it to be most
prominent in general.

###  5.2 Quark mass dependence

Let us see how the situation described in the previous subsection changes with
increasing quark mass. The quark and pion mass values of our study are given
in (3.3) and (3.3). With the strange quark mass, we have data only for the
contractions $C_{1}$ and $A$, whilst for charm we have results for all
contractions except $S_{2}$. The errors on $D$ are again huge in this case,
which we will discuss no further here.

(a) $\langle V^{0}V^{0}\rangle,|\mskip 1.5muC_{1}|,p=0,L=40$ (b) $\langle
A^{0}A^{0}\rangle,|\mskip 1.5muC_{1}|,p=0,L=40$ (c) $\langle SS\rangle,|\mskip
1.5muC_{1}|,p=0,L=40$ (d) $\langle PP\rangle,|\mskip 1.5muC_{1}|,p=0,L=40$

Figure 12:  Comparison of $|\mskip 1.5muC_{1}|$ for the three different quark
masses of our study. For clarity, we have omitted data points for which the
statistical error is too large compared with the central value. The $\langle
PP\rangle$ correlator undergoes a sign change for strange and charm quarks.

We compare the correlator $C_{1}$ for our three quark masses in figure 12. As
is to be expected, the behaviour for charm quarks differs more strongly from
the two other cases than light and strange quarks do from each other. The
steeper decrease with $y$ for charm quarks is also expected and indicates that
a “pion” made from heavy quarks is more compact. We will return to this in
section 5.4. In the case of $\langle PP\rangle$, we observe that $C_{1}$
crosses zero (seen as a sharp dip in the logarithmic plot). This happens
around $0.95\operatorname{fm}$ for strange quarks and around
$0.35\operatorname{fm}$ for charm.

(a) $\langle V^{0}V^{0}\rangle,p=0,L=40$, charm mass (b) $\langle
A^{0}A^{0}\rangle,p=0,L=40$, charm mass (c) $\langle SS\rangle,p=0,L=40$,
charm mass (d) $\langle PP\rangle,p=0,L=40$, charm mass

Figure 13:  As figure 11 but for charm quarks. The contribution from the
annihilation graph $A$ is not shown for $\langle A^{0}A^{0}\rangle$ and
$\langle PP\rangle$, where its relatively large errors would obscure the
plots.

Not shown in our plots is the annihilation contribution for strange quarks. It
is very small compared with $C_{1}$, except in the $\langle SS\rangle$
channel, where the ratio $-A/C_{1}$ is about $40\%$ at $y\sim
0.5\operatorname{fm}$ and decreases with $y$. The annihilation mechanism thus
significantly loses importance when going from light to strange quark masses.
With charm quarks, we find $A$ to be negligible in all channels.

In figure 13 we show different contractions for charm quarks. Compared with
the case of light quarks in figure 11, we see that not only the importance of
$A$ but also the one of $S_{1}$ is considerably reduced for charm. In turn,
$C_{2}$ remains important compared with $C_{1}$ for both $\langle
A^{0}A^{0}\rangle$ and $\langle PP\rangle$. The expectation that $C_{1}$ is
the dominant contraction in each channel is thus not even realised for charm
quarks.

Looking at the largest contraction in each current combination, we observe
that for charm quarks $\langle A^{0}A^{0}\rangle$ is small compared to
$\langle V^{0}V^{0}\rangle$ and $\langle PP\rangle$ relatively small compared
to $\langle SS\rangle$. By contrast, for light quarks $\langle
A^{0}A^{0}\rangle$ is comparable in size to $\langle V^{0}V^{0}\rangle$, and
$\langle PP\rangle$ is comparable in size to $\langle SS\rangle$. Our results
for charm quarks confirm that in the heavy-quark limit, matrix elements of the
spin dependent currents $A^{\mu}$ and $P$ are suppressed compared with their
spin independent counterparts $V^{\mu}$ and $S$. Seen from this angle, the
contraction $C_{1}$ is indeed dominant for charm because it is the most
prominent contribution to $\langle V^{0}V^{0}\rangle$ and $\langle SS\rangle$.

###  5.3 Test of factorisation

The two-current matrix elements $\langle V^{0}V^{0}\rangle$ and $\langle
SS\rangle$ quantify the distribution of vector or scalar charges at two
separate points in the pion. By contrast, the corresponding single-current
matrix elements $\bra{\pi(p^{\prime})}V^{0}\ket{\pi(p)}$ and
$\bra{\pi(p^{\prime})}S\ket{\pi(p)}$ give the vector or scalar charge
distribution at one point in the pion if one performs a Fourier transform from
the momentum transfer $\vec{\mskip 0.5mup}\mskip 1.5mu-\vec{\mskip
0.5mup}\mskip 1.5mu^{\prime}$ to three-dimensional space.333As is well known
and discussed e.g. in Burkardt:2000za , this simple interpretation receives
relativistic corrections, which become relevant for distances of order
$1/m_{\pi}$ or below. Our treatment here is fully relativistic and not
affected by this limitation. In the absence of correlations, the two-current
distribution can be computed from the single-current one. A detailed analysis
of this has already been given in Burkardt:1994pw , with some focus on the
non-relativistic limit. A way to express the expectation for the two-current
distribution in the absence of correlations is to insert a complete set of
intermediate states between the two currents and to keep only the pion ground
state Burkardt:1994pw ; Diehl:2011yj . In momentum space, we have

$\displaystyle\mathcal{M}_{ii}(\vec{\mskip 0.5muq}\mskip 1.5mu^{2})$
$\displaystyle=\int d^{3}{y}\;e^{i\vec{\mskip 0.5muq}\mskip 1.5mu\vec{\mskip
0.5muy}\mskip 1.5mu}\,\bra{\pi^{+}(p)}\mathcal{O}_{i}^{uu}(y)\mskip
1.5mu\mathcal{O}_{i}^{dd}(0)\ket{\pi^{+}(p)}$ $\displaystyle\overset{?}{=}\int
d^{3}{y}\;e^{i\vec{\mskip 0.5muq}\mskip 1.5mu\vec{\mskip 0.5muy}\mskip
1.5mu}\int\frac{d^{3}{p}^{\prime}}{(2\pi)^{3}\mskip 1.5mu2p^{\prime\mskip
1.5mu0}}\bra{\pi^{+}(p)}\mathcal{O}_{i}^{uu}(y)\ket{\pi^{+}(p^{\prime})}\,\bra{\pi^{+}(p^{\prime})}\mathcal{O}_{i}^{dd}(0)\ket{\pi^{+}(p)}$
$\displaystyle=\frac{\eta^{i}_{C}}{2E_{q}}\;\bigl{|}\bra{\pi^{+}(E_{q},-\vec{\mskip
0.5muq}\mskip 1.5mu)}\mathcal{O}_{i}^{uu}(0)\ket{\pi^{+}(p)}\bigr{|}^{2}\,,$
(63)

where $E^{2}_{q}=m_{\pi}^{2}+\vec{\mskip 0.5muq}\mskip 1.5mu^{2}$ and it is
understood that $\vec{\mskip 0.5mup}\mskip 1.5mu=\vec{\mskip 0.5mu0}\mskip
1.5mu$. The question mark above the equality sign indicates an assumption,
which will be tested in the following. In the last step of (5.3) we used the
relation
$\bra{\pi^{+}}\mathcal{O}_{i}^{dd}\ket{\pi^{+}}=\eta^{i}_{C}\,\bra{\pi^{+}}\mathcal{O}_{i}^{uu}\ket{\pi^{+}}$
with the $C$ parity $\eta^{i}_{C}$ defined in (6). Defining vector and scalar
form factors as

$\displaystyle\bra{\pi^{+}(p^{\prime})}V_{uu}^{\mu}(0)\ket{\pi^{+}(p)}$
$\displaystyle=(p^{\prime}+p)^{\mu}\mskip 1.5muF_{V}(Q^{2})\,,$
$\displaystyle\bra{\pi^{+}(p^{\prime})}S_{uu}(0)\ket{\pi^{+}(p)}$
$\displaystyle=F_{S}(Q^{2})$ (64)

with $Q^{2}=-(p-p^{\prime})^{2}$, we get

$\displaystyle{}-\mathcal{M}_{V^{0}V^{0}}(\vec{\mskip 0.5muq}\mskip
1.5mu^{2})$
$\displaystyle\overset{?}{=}{}\frac{(m_{\pi}+E_{q})^{2}}{2E_{q}}\,\Bigl{[}\mskip
1.5muF_{V}\bigl{(}2m_{\pi}E_{q}-2m_{\pi}^{2}\bigr{)}\mskip
1.5mu\Bigr{]}^{2}\,,$ $\displaystyle\mathcal{M}_{SS}(\vec{\mskip 0.5muq}\mskip
1.5mu^{2})$ $\displaystyle\overset{?}{=}{}\frac{1}{2E_{q}}\,\Bigl{[}\mskip
1.5muF_{S}\bigl{(}2m_{\pi}E_{q}-2m_{\pi}^{2}\bigr{)}\mskip
1.5mu\Bigr{]}^{2}\,.$ (65)

In momentum space, the absence of correlations thus implies that the two-
current matrix element factorises into the product of two single-current form
factors (and a kinematic prefactor). To compare the matrix elements in
position space, we analytically Fourier transform the r.h.s. of (5.3) using a
parameterisation of the form factors. We note in passing that in the limit
$\vec{\mskip 0.5muq}\mskip 1.5mu^{2}\ll m_{\pi}^{2}$ the argument of the form
factors in (5.3) becomes $\vec{\mskip 0.5muq}\mskip 1.5mu^{2}$, as appropriate
for a non-relativistic treatment. For our study with $m_{\pi}\approx
300\operatorname{MeV}$, this limit is however not relevant.

It follows from (3.2) that the two-current correlators in (5.3) receive
contributions from different contractions in the combination $C_{1}+2S_{1}+D$.
Given the poor quality of our data for $D$, we formulate a version of the
factorisation hypothesis that requires only $C_{1}$ for the two-current
correlator and the connected three-point graph $G_{3}$ for the elastic form
factors. To this end, we use the partially quenched scenario of $n_{F}=4$
mass-degenerate quarks $u$, $d$, $s$ and $c$, described at the end of section
3.2. In this case, one starts from the matrix element
$\bra{\pi^{+}}\mathcal{O}_{i}^{uc}(y)\mskip
1.5mu\mathcal{O}_{i}^{sd}(0)\ket{D_{s}^{+}}$ and inserts a full set of
intermediate states with quark content $c\bar{d}$. Retaining only the ground
state term $\ket{D^{+}}\bra{D^{+}}$ and using SU(4) flavour symmetry to relate
the two single-current matrix elements to each other, one obtains the analogue
of (5.3). Disconnected contributions are in this case excluded by the quantum
numbers of the $uc$ and $sd$ currents. Our lattice computation can be
understood as giving the corresponding matrix elements, up to the effect of
partial quenching due to the lattice action with $n_{F}=2$.

#### Factorisation at $\boldsymbol{q=0}$.

Let us first test the factorisation hypothesis (5.3) at $\vec{\mskip
0.5muq}\mskip 1.5mu=\vec{\mskip 0.5mu0}\mskip 1.5mu$, where it reads

$\displaystyle{}-(2m_{\pi})^{-1}\int d^{3}{y}\;\langle
V^{0}V^{0}\rangle(\vec{\mskip 0.5muy}\mskip 1.5mu)$
$\displaystyle\overset{?}{=}1\,,$ $\displaystyle 2m_{\pi}\int
d^{3}{y}\;\langle SS\rangle(\vec{\mskip 0.5muy}\mskip 1.5mu)$
$\displaystyle\overset{?}{=}\bigl{[}F_{S}(0)\mskip 1.5mu\bigr{]}^{2}\,.$ (66)

For the vector current we used the normalisation condition $F_{V}(0)=1$. We
evaluate the integrals on the l.h.s. of (66) as discrete sums over all lattice
sites.444In this case we do not use the cut in (62). As discussed in section
4.2, this leaves us with lattice artefacts from images at large $y$, but this
region is not important in the sum over $\vec{\mskip 0.5muy}\mskip 1.5mu$. The
result for the vector correlator is

$\displaystyle{}-(2m_{\pi})^{-1}\,a^{3}\,\smash{\sum_{\vec{\mskip
0.5muy}\mskip 1.5mu}}\langle V^{0}V^{0}\rangle(\vec{\mskip 0.5muy}\mskip
1.5mu)$ $\displaystyle=0.975(33)\text{ \leavevmode\nobreak\
\leavevmode\nobreak\ (light quarks), }$ $\displaystyle=0.985(31)\text{
\leavevmode\nobreak\ \leavevmode\nobreak\ (strange), }$
$\displaystyle=1.053(36)\text{ \leavevmode\nobreak\ \leavevmode\nobreak\
(charm)}$ (67)

for the three quark masses used in our study. In all cases, the agreement with
factorisation is excellent, despite the possibility of substantial lattice
spacing effects for the charm quark.

In the scalar channel we obtain

$\displaystyle\biggl{[}2m_{\pi}\mskip 1.5mua^{3}\sum_{\vec{\mskip
0.5muy}\mskip 1.5mu}\;\langle SS\rangle(\vec{\mskip 0.5muy}\mskip
1.5mu)\biggr{]}^{1/2}$ $\displaystyle=3.14(4)\operatorname{GeV}$ (68)

for light quarks, which is to be compared with the value
$F_{S}(0)=2.2\operatorname{GeV}$ we extract from our scalar form factor data
(see table 4 below). While this is of the same order of magnitude,
factorisation clearly fails in this case.

Let us give a heuristic argument why the factorisation hypothesis works so
well in the vector channel. $V_{qq}^{\mu}(y)$ is a conserved current in QCD,
so that the associated Noether charge $\rho_{q}=\int
d^{3}{y}\;V_{qq}^{\mu}(y)$ corresponds to a good quantum number. For a
positive pion one has $\bra{\pi^{+}}\rho_{u}=\bra{\pi^{+}}$ and thus obtains

$\displaystyle\int d^{3}{y}\,\bra{\pi^{+}}V_{uu}^{0}(y)\mskip
1.5muV_{dd}^{0}(0)\ket{\pi^{+}}$ $\displaystyle=\bra{\pi^{+}}\rho_{u}\mskip
1.5muV_{dd}^{0}(0)\ket{\pi^{+}}=\bra{\pi^{+}}V_{dd}^{0}(0)\ket{\pi^{+}}={}-2m_{\pi}$
(69)

for zero pion momentum. To turn this argument into a theorem, one would need
to discuss possible short-distance singularities of the integrated two-current
matrix element at $\vec{\mskip 0.5muy}\mskip 1.5mu=\vec{\mskip 0.5mu0}\mskip
1.5mu$. We shall not attempt this here. It is clear that the argument just
given cannot be extended to the scalar charge, which is not the Noether charge
of a conserved current.

#### Factorisation as a function of ${\boldsymbol{y}}$.

We have extracted the vector form factor and the connected part of the scalar
form factor from our lattice simulations, using in this case the full number
of 2025 gauge configurations available for our lattice with $L=40$. To Fourier
transform the r.h.s. of (5.3) to position space, we fit the form factor data
to a power law

$\displaystyle F(Q^{2})$
$\displaystyle=\frac{F(0)}{\bigl{(}1+Q^{2}/M^{2}\bigr{)}^{n}\rule{0.0pt}{10.50005pt}}\,.$
(70)

For each form factor, we take two fit variants with different powers $n$, so
as to have a handle on the bias of the extrapolation to large $Q^{2}$, where
we have no data. Such an extrapolation bias is inevitable when we Fourier
transform to position space.

fit | form factor | $F(0)$ | $M[\operatorname{GeV}]$ | $n$
---|---|---|---|---
$1$ | $F_{V}$ | $1\text{\leavevmode\nobreak\ (fixed)}$ | $0.777(12)$ | $1\text{\leavevmode\nobreak\ (fixed)}$
$2$ |  | $1\text{\leavevmode\nobreak\ (fixed)}$ | $0.872(16)$ | $1.173(69)$
$1$ | $F_{S}$ | $2.222(19)\operatorname{GeV}$ | $1.314(39)$ | $1\text{\leavevmode\nobreak\ (fixed)}$
$2$ |  | $2.212(19)\operatorname{GeV}$ | $2.023(50)$ | $2\text{\leavevmode\nobreak\ (fixed)}$
Table 4:  Parameters obtained in fits of our form factor data to (70).

The extracted form factor data and the results of the fits are shown in figure
15, and the fitted parameters are given in table 4. We see that the fits
describe the $F_{V}$ data very well. For $F_{S}$, the quality of the data is
less good, and so is the agreement between our fits and the data points, but
we consider this sufficient for our present study.

(a) fits of $F_{V}(Q^{2})$, $L=40$ (b) fits of $F_{S}(Q^{2})$, $L=40$

Figure 14:  Data for the vector and scalar pion form factors extracted from
our simulations, together with the fits specified in table 4. $F_{S}$ is given
in the $\overline{\mathrm{MS}}$ scheme at scale $\mu=2\operatorname{GeV}$.

(a) factorisation test for ${}-\langle V^{0}V^{0}\rangle,L=40$ (b)
factorisation test for $\langle SS\rangle,L=40$

Figure 15:  Test of the factorisation hypothesis for the two-current
correlators $\langle V^{0}V^{0}\rangle$ and $\langle SS\rangle$. The lattice
data corresponds to l.h.s. of (5.3) in position space, and the curves with
error bands are obtained by Fourier transforming the r.h.s. of (5.3)
accordingly, using the form factor parameterisations specified in (70) and
table 4.

Having determined the form factors, we can evaluate the Fourier transforms of
(5.3) and compare the result with our lattice data for $C_{1}$ as a function
of $y$. We see in figure 15 that the Fourier transforms of the two form factor
parameterisations agree very well in the $y$ range shown. For both currents,
the factorisation hypothesis fails very clearly over a wide range of $y$,
indicating significant correlation effects for the distribution of both vector
and scalar charge in a pion of mass $300\operatorname{MeV}$.

###  5.4 Root mean square radii

In figure 12 we noted a clear difference between the decrease of the two-
current correlators for light and heavy quarks, reflecting the different size
of the pseudoscalar bound state in the two cases. In the present section, we
make this observation more quantitative. We focus on the current combinations
$\langle V^{0}V^{0}\rangle$ and $\langle A^{0}A^{0}\rangle$ and on graph
$C_{1}$, which can be interpreted as a physical matrix element in the sense
discussed at the end of section 3.2. Following Wilcox:1986ge  and later work,
we quantify the characteristic length scale of the two-current correlators via
root mean square (rms) radii $r_{VV}$ and $r_{AA}$, defined by

$\displaystyle r^{2}_{VV}$ $\displaystyle=\frac{\int
d^{3}y\;{y}^{2}\,C_{1}^{\mskip 1.5muVV}(\vec{\mskip 0.5muy}\mskip 1.5mu)}{\int
d^{3}y\;C_{1}^{\mskip 1.5muVV}(\vec{\mskip 0.5muy}\mskip 1.5mu)}$ (71)

for $\langle V^{0}V^{0}\rangle$ and in analogy for $\langle
A^{0}A^{0}\rangle$. For the sake of legibility, we omit superscripts $0$ in
the labels $VV$ and $AA$ in the remainder of this section.

One might think of evaluating (71) from the ratio

$\displaystyle\sum_{\vec{\mskip 0.5muy}\mskip 1.5mu\in E(\vec{\mskip
0.5mu0}\mskip 1.5mu)}y^{2}\,C_{1}^{\mskip 1.5muVV}(\vec{\mskip 0.5muy}\mskip
1.5mu)\Bigg{/}\sum_{\vec{\mskip 0.5muy}\mskip 1.5mu\in E(\vec{\mskip
0.5mu0}\mskip 1.5mu)}C_{1}^{\mskip 1.5muVV}(\vec{\mskip 0.5muy}\mskip 1.5mu)$
(72)

of lattice data, with the sum running over all $\vec{\mskip 0.5muy}\mskip
1.5mu$ in the elementary lattice cell $E(\vec{\mskip 0.5mu0}\mskip 1.5mu)$
defined by $-L/2<y^{i}/a\leq L/2$ for $i=1,2,3$. However, the weight $y^{2}$
in the numerator of (72) enhances the importance of large distances to the
extent that for our $L=40$ lattice, the contribution of points with $y/a>L/2$
to the sum over all $\vec{\mskip 0.5muy}\mskip 1.5mu$ amounts to 39% for light
and to 16% for strange quarks. At such distances, the effect of periodic
images discussed in section 4.2 is considerable. Moreover, the slow decrease
of $y^{2}\,C_{1}^{\mskip 1.5muVV}(\vec{\mskip 0.5muy}\mskip 1.5mu)$ with $y$
implies that the numerator of (72) misses important contributions from
distances $y/a>\sqrt{3}\mskip 1.5muL/2$ that are not included in the
elementary cell.

In the following we present a method to deal with this situation. Let us
assume that the correlation function computed on the lattice has the form

$\displaystyle C_{1}^{\text{lat}}(\vec{\mskip 0.5muy}\mskip 1.5mu)$
$\displaystyle=\sum_{\vec{\mskip 0.5mun}\mskip
1.5mu\in\mathbb{Z}_{3}}C_{1}(\vec{\mskip 0.5muy}\mskip 1.5mu+\vec{\mskip
0.5mun}\mskip 1.5muLa)\,,$ (73)

where $C_{1}(\vec{\mskip 0.5muy}\mskip 1.5mu)$ is the correlation function in
the physical limit and the terms with $\vec{\mskip 0.5mun}\mskip
1.5mu\neq\vec{\mskip 0.5mu0}\mskip 1.5mu$ are due to the periodic boundary
conditions on the lattice Burkardt:1994pw . The rms radii can be computed from
the Fourier transform

$\displaystyle\widetilde{C}_{1}(\vec{\mskip 0.5muq}\mskip 1.5mu)$
$\displaystyle=\int d^{3}y\;e^{i\vec{\mskip 0.5muq}\mskip 1.5mu\vec{\mskip
0.5muy}\mskip 1.5mu}\,C_{1}(\vec{\mskip 0.5muy}\mskip 1.5mu)$ (74)

as

$\displaystyle r^{2}$ $\displaystyle={}-6\mskip
1.5mu\bigl{[}\widetilde{C}_{1}(\vec{\mskip 0.5mu0}\mskip
1.5mu)\bigr{]}^{-1}\,\frac{d\mskip 1.5mu\widetilde{C}_{1}(\vec{\mskip
0.5muq}\mskip 1.5mu)}{d\mskip 1.5muq^{2}}\biggl{|}_{q^{2}=0}\,,$ (75)

where we used that $\widetilde{C}_{1}$ depends on $\vec{\mskip 0.5muq}\mskip
1.5mu$ only via $q=|\vec{\mskip 0.5muq}\mskip 1.5mu|$. On the lattice we can
evaluate the discrete Fourier transform

$\displaystyle\widetilde{C}_{1}^{\text{lat}}(\vec{\mskip 0.5muq}\mskip 1.5mu)$
$\displaystyle=a^{3}\\!\\!\sum_{\vec{\mskip 0.5muy}\mskip 1.5mu\in
E(\vec{\mskip 0.5mu0}\mskip 1.5mu)}e^{i\vec{\mskip 0.5muq}\mskip
1.5mu\vec{\mskip 0.5muy}\mskip 1.5mu}\,C_{1}^{\text{lat}}(\vec{\mskip
0.5muy}\mskip 1.5mu)$ $\displaystyle\text{ for }\vec{\mskip 0.5muq}\mskip
1.5mu=\vec{\mskip 0.5muk}\mskip 1.5mup_{0}$ (76)

with $\vec{\mskip 0.5muk}\mskip 1.5mu\in\mathbb{Z}_{3}$ and $p_{0}=2\pi/(La)$.
We have

$\displaystyle\widetilde{C}_{1}^{\text{lat}}(\vec{\mskip 0.5muq}\mskip 1.5mu)$
$\displaystyle=a^{3}\sum_{\vec{\mskip 0.5mun}\mskip
1.5mu\in\mathbb{Z}_{3}}\;\sum_{\vec{\mskip 0.5muy}\mskip 1.5mu\in
E(\smash{\vec{\mskip 0.5mu0}\mskip 1.5mu})}e^{i\vec{\mskip 0.5muq}\mskip
1.5mu\vec{\mskip 0.5muy}\mskip 1.5mu}\,C_{1}(\vec{\mskip 0.5muy}\mskip
1.5mu+\vec{\mskip 0.5mun}\mskip 1.5muLa)=a^{3}\sum_{\vec{\mskip 0.5mun}\mskip
1.5mu\in\mathbb{Z}_{3}}\;\sum_{\vec{\mskip 0.5muy}\mskip 1.5mu\in
E(\vec{\mskip 0.5mun}\mskip 1.5mu)}e^{i\vec{\mskip 0.5muq}\mskip
1.5mu\vec{\mskip 0.5muy}\mskip 1.5mu-i\vec{\mskip 0.5muq}\mskip
1.5mu\vec{\mskip 0.5mun}\mskip 1.5muLa}\,C_{1}(\vec{\mskip 0.5muy}\mskip
1.5mu)$ $\displaystyle=a^{3}\\!\\!\sum_{\vec{\mskip 0.5muy}\mskip
1.5mu/a\in\mathbb{Z}_{3}}e^{i\vec{\mskip 0.5muq}\mskip 1.5mu\vec{\mskip
0.5muy}\mskip 1.5mu}\,C_{1}(\vec{\mskip 0.5muy}\mskip 1.5mu)\,,$ (77)

where $E(\vec{\mskip 0.5mun}\mskip 1.5mu)$ denotes the shifted lattice cell
defined by $-L/2<y^{i}/a+n^{i}\leq L/2$ for $i=1,2,3$, and in the last step we
have used the condition $\vec{\mskip 0.5muq}\mskip 1.5mu=\vec{\mskip
0.5muk}\mskip 1.5mup_{0}$. We thus find that for $a\to 0$ the discrete Fourier
transform (76) becomes equal to the infinite-volume expression (74). The
periodic images included in $\widetilde{C}_{1}^{\text{lat}}(\vec{\mskip
0.5muq}\mskip 1.5mu)$ provide the contribution of the infinite-volume
correlator $C_{1}(\vec{\mskip 0.5muy}\mskip 1.5mu)$ at distances outside the
elementary cell $E(\vec{\mskip 0.5mu0}\mskip 1.5mu)$.

To evaluate the rms radii from (75), we need to construct a smooth function of
$\vec{\mskip 0.5muq}\mskip 1.5mu$ out of $\widetilde{C}_{1}(\vec{\mskip
0.5muq}\mskip 1.5mu)$ at the points $\vec{\mskip 0.5muq}\mskip
1.5mu=\vec{\mskip 0.5muk}\mskip 1.5mup_{0}$. Whilst direct computation of (71)
would require us to _extrapolate_ $C_{1}(\vec{\mskip 0.5muy}\mskip 1.5mu)$ to
large values of $y$ and to remove the contributions from periodic images, we
now need to _interpolate_ between discrete values of
$\widetilde{C}_{1}(\vec{\mskip 0.5muq}\mskip 1.5mu)$ in the vicinity of
$\vec{\mskip 0.5muq}\mskip 1.5mu=\vec{\mskip 0.5mu0}\mskip 1.5mu$. The
reliability of this interpolation is of course higher for a higher density of
points $\vec{\mskip 0.5muk}\mskip 1.5mup_{0}=2\pi\vec{\mskip 0.5muk}\mskip
1.5mu/(La)$ and thus for larger physical lattice size $La$.

In figure 16 we show $\widetilde{C}_{1}^{\mskip 1.5muVV}(\vec{\mskip
0.5muq}\mskip 1.5mu)$ obtained from our lattice data for light and charm
quarks. The results for strange quarks and those for
$\widetilde{C}_{1}^{\mskip 1.5muAA}(\vec{\mskip 0.5muq}\mskip 1.5mu)$ look
qualitatively similar. At values $q\sim\pi/a$ we see a clear anisotropy in the
Fourier transform (which in the continuum limit depends on the length but not
on the direction of $\vec{\mskip 0.5muq}\mskip 1.5mu$). This is to be expected
and reflects discretisation effects in $C_{1}(\vec{\mskip 0.5muy}\mskip
1.5mu)$ at distances $y$ of a few lattice units $a$. With decreasing $q$, the
anisotropy gradually disappears and we have a very clear and smooth signal.

Figure 16:  The Fourier transform of the correlator $C_{1}(\vec{\mskip
0.5muy}\mskip 1.5mu)$, evaluated using (76). The values for charm quarks are
scaled by $0.05$ for the sake of clarity. Statistical uncertainties on the
data points are not visible at this scale. We divide the Fourier transform by
$-2m_{\pi}$, so that according to our findings in section 5.3 we have a simple
normalisation $-\widetilde{C}_{1}^{\mskip 1.5muVV}(\vec{\mskip 0.5mu0}\mskip
1.5mu)/(2m_{\pi})=1$.

(a) $\langle V^{0}V^{0}\rangle,p=0$ (b) $\langle A^{0}A^{0}\rangle,p=0$

Figure 17:  Volume comparison of the Fourier transform $\widetilde{C}_{1}$ for
light quarks.

In figure 17 we focus on the small $q$ region and compare our results for
light quarks on the lattices with $L=32$ and $L=40$. Although there is some
indication for a systematic shift of $\widetilde{C}_{1}$ between the two
lattices, we find the agreement very satisfactory. This indicates that the
main finite size effect in the two-current correlators is indeed due to
periodic images, which was the hypothesis underlying our arguments following
(73).

(a) $\langle V^{0}V^{0}\rangle,p=0,L=40$, light quarks (b) $\langle
V^{0}V^{0}\rangle,p=0,L=40$, strange and charm (c) $\langle
A^{0}A^{0}\rangle,p=0,L=40$

Figure 18:  The Fourier transform $\widetilde{C}_{1}$ at low values of
$q^{2}$. The curves (without error bands) show the series of fits discussed in
the text. In panel (a) the results for $-\widetilde{C}_{1}/(2m_{\pi})$ are
labelled by “4pt FT light”, whereas the points and curve labelled by “3pt
light” show the r.h.s. of (79) evaluated with our data for the vector form
factor $F_{V}$ and fit 1 of table 4.

(a) $\langle V^{0}V^{0}\rangle,p=0,L=40$ (b) $\langle
A^{0}A^{0}\rangle,p=0,L=40$

Figure 19:  Squared rms radii $r^{2}$ extracted from fitting
$\widetilde{C}_{1}$ to the form (78) in the interval $0\leq q^{2}\leq
q^{2}_{\text{max}}$, together with their jackknife errors. The entries
labelled “avg” give the averaged values corresponding to table 5, with
statistical and systematic errors added in quadrature.

In order to compute the derivative in $q^{2}$ around $q=0$, we fit our results
to the form

$\displaystyle\widetilde{C}_{1}(\vec{\mskip 0.5muq}\mskip 1.5mu)$
$\displaystyle=\frac{\widetilde{C}_{1}(\vec{\mskip 0.5mu0}\mskip
1.5mu)}{\bigl{[}1+r^{2}q^{2}/(6n)\bigr{]}^{n}}$ (78)

in the range $0\leq q^{2}\leq q^{2}_{\text{max}}$, varying the value of
$q^{2}_{\text{max}}$ between $2\mskip 1.5mup_{0}^{2}$ and $6\mskip
1.5mup_{0}^{2}$. For the lowest value of $q^{2}_{\text{max}}$ we have as many
data points as parameters in (78), whereas the highest value is motivated by
the fact that there is no point with $q^{2}=7p_{0}^{2}$ because $7$ cannot be
written as the sum of squares of three integers. Fits extending to even higher
$q^{2}$ put less and less emphasis on the fit quality around $q^{2}=0$ and are
therefore less well suited for determining the derivative at that point.
Figure 18 shows the low $q^{2}$ data together with the different fit curves.
The latter are in general very close to each other and can barely be
distinguished individually. In figure 19 we show the fitted values of
$r^{2}_{VV}$ and $r^{2}_{AA}$. The dependence of $r^{2}$ on
$q^{2}_{\text{max}}$ is in general mild, although for light or strange quarks
there seems to be a systematic increase of $r^{2}_{AA}$ with
$q^{2}_{\text{max}}$. We do not regard the fits with lower or higher
$q^{2}_{\text{max}}$ as intrinsically more reliable around $q^{2}=0$ and
therefore determine an overall value and error for each radius in the
following way. The averaged value of $r^{2}$ is obtained as the arithmetic
mean of the fitted values $r^{2}(q^{2}_{\text{max}})$, and the root of the
variance of these values is used as systematic error on $r^{2}$ due to our
fitting procedure. The statistical error on $r^{2}$ is taken as the arithmetic
mean of the jackknife errors on the fit results $r^{2}(q^{2}_{\text{max}})$.
The outcome of this procedure is given in table 5, converted to values and
errors for $r$ instead of $r^{2}$.

radius [fm] | light quarks | strange | charm
---|---|---|---
$r_{VV}$ | 1.046(40)(09) | 0.826(32)(09) | 0.460(60)(25)
$r_{AA}$ | 0.580(15)(09) | 0.576(21)(08) | 0.380(38)(12)
Table 5:  Rms radii $r$ for the correlator $C_{1}(y)$. The first error is
statistical and the second one due to the fitting procedure, as specified in
the text.

We see that $r_{VV}$ decreases with increasing quark mass, as expected for a
quantity that characterises the size of the pion. The ratio $r_{AA}/r_{VV}$
increases with the quark mass, being clearly below $1$ for light and strange
quarks and consistent with $1$ for charm. A dynamical interpretation of this
interesting finding is beyond the scope of the present work.

The factorisation hypothesis for $C_{1}^{\mskip 1.5muVV}$ discussed in section
5.3 gives the relation

$\displaystyle{}-\frac{\widetilde{C}_{1}(\vec{\mskip 0.5muq}\mskip
1.5mu)}{2m_{\pi}}$
$\displaystyle\overset{?}{=}\frac{(m_{\pi}+E_{q})^{2}}{4m_{\pi}E_{q}}\,\Bigl{[}\mskip
1.5muF_{V}\bigl{(}2m_{\pi}E_{q}-2m_{\pi}^{2}\bigr{)}\mskip 1.5mu\Bigr{]}^{2}$
(79)

according to (5.3). We see in figure 18a that this is clearly ruled out, which
confirms our finding in position space in figure 15a.

For rms radii, the hypothesis (79) implies $r_{VV}^{2}\overset{?}{=}2\mskip
1.5mur_{V}^{2}$, where $r_{V}$ is the radius associated with the vector form
factor $F_{V}$, i.e. the conventional charge radius. Interestingly, the
kinematic prefactor on the r.h.s. of (79) does not contribute to $d/dq^{2}$ at
$q^{2}=0$, so that for the radii one obtains the same result as in the non-
relativistic approximation, where the prefactor is $1$. Quantitatively, we
obtain $\sqrt{2}\,r_{V}=0.878(14)\operatorname{fm}$ with fit 1 and
$\sqrt{2}\,r_{V}=0.849(12)\operatorname{fm}$ with fit 2 of table 4, which is
clearly below the value $r_{VV}=1.046(40)(09)\operatorname{fm}$ for light
quarks in table 5. A detailed discussion of the difference between $r_{VV}$
and $\sqrt{2}\,r_{V}$ can be found in Burkardt:1994pw .

###  5.5 Subtraction term for the annihilation graph

The vacuum subtraction term for the annihilation graph $A$ in (3.3) involves
the same three-point functions $\langle B_{t}\rangle$ that are necessary to
compute the matrix elements of the axial or pseudoscalar current between a
pion state and the vacuum. As a by-product of our simulations, we can thus
extract these matrix elements. We can then compute the pion decay constant
from the relation

$\displaystyle\bigl{|}\mskip 1.5mu\bra{\pi^{+}}A^{0}_{ud}(0)\ket{0}\bigr{|}$
$\displaystyle=\sqrt{2}\mskip 1.5mum_{\pi}\mskip 1.5muF_{\pi}$ (80)

for a pion at rest. Using $\partial_{\mu}A^{\smash{\mu}}_{ud}(x)=2m_{q}\mskip
1.5muP_{ud}(x)$, where $m_{q}$ is the average of $u$ and $d$ quark masses, we
have for a pion at rest

$\displaystyle\biggl{|}\frac{\bra{\pi^{+}}P_{ud}(0)\ket{0}}{\bra{\pi^{+}}A^{0}_{ud}(0)\ket{0}}\biggr{|}$
$\displaystyle=\frac{m_{\pi}}{2m_{q}}=\frac{B}{m_{\pi}}\,,$ (81)

where the last relation holds at leading order in chiral perturbation theory.
The quark mass $m_{q}$ and the chiral symmetry breaking parameter $B$ are
understood to be renormalised in the $\overline{\mathrm{MS}}$ scheme at the
scale $\mu=2\operatorname{GeV}$, corresponding to the renormalisation
constants given in table 3. From the data of our two lattices, we obtain

$\displaystyle F_{\pi}$ $\displaystyle=(100.2\pm 0.4)\operatorname{MeV}\,,$
$\displaystyle B$ $\displaystyle=(2.41\pm 0.02)\operatorname{GeV}$
$\displaystyle\text{ \leavevmode\nobreak\ for\leavevmode\nobreak\ }L=40$
$\displaystyle(m_{\pi}=293\operatorname{MeV})\,,$ $\displaystyle F_{\pi}$
$\displaystyle=(98.2\pm 1.0)\operatorname{MeV}\,,$ $\displaystyle B$
$\displaystyle=(2.50\pm 0.05)\operatorname{GeV}$ $\displaystyle\text{
\leavevmode\nobreak\ for\leavevmode\nobreak\ }L=32$
$\displaystyle(m_{\pi}=299\operatorname{MeV})\,,$

where we also recall the pion masses fitted from our two-point correlation
functions. The quoted errors are purely statistical.

The above values of $F_{\pi}$ for $m_{\pi}\approx 300\operatorname{MeV}$ agree
within a few percent with chiral perturbation theory at NLO, see figure 10 in
Durr:2014oba . From equations (46), (65) and table 14 in the FLAG review
Aoki:2013ldr , we obtain $B=\Sigma/F^{2}\approx 2.6\text{\leavevmode\nobreak\
to\leavevmode\nobreak\ }2.7\operatorname{GeV}$ for $n_{F}=2$ quark flavours in
the chiral limit. This agrees reasonably well with our values in (5.5), given
that we use lowest-order chiral perturbation theory for their extraction and
have not attempted to correct our data for lattice artefacts.

##  6 Results for isospin amplitudes

In this section we present our results for the isospin amplitudes $F_{0}$,
$F_{1}$ and $F_{2}$ defined in (2.1). We compare them with the predictions of
chiral perturbation theory that were computed in Bruns:2015yto  and recalled
in section 2.2 of the present paper. We consider the case of light quarks
throughout.

###  6.1 Isospin amplitudes

In figures 20 and 21, we show our results for the isospin amplitudes in all
channels except $\langle PP\rangle$ for $F_{0}$ and $\langle SS\rangle$ for
$F_{1}$, where the statistical errors are too large to obtain a nonzero
signal. The $y$ range shown is selected such that regions where we see no
signal have been omitted.

(a) $F_{1},\langle V^{0}V^{0}\rangle,p=0$ (b) $F_{1},\langle
A^{0}A^{0}\rangle,p=0$ (c) $F_{0},\langle V^{0}V^{0}\rangle,p=0$ (d)
$F_{2},\langle V^{0}V^{0}\rangle,p=0$ (e) $F_{2},\langle
A^{0}A^{0}\rangle,p=0$ (f) $F_{2},\langle A^{0}A^{0}\rangle,p=0$

Figure 20:  Lattice results for isospin amplitudes of $\langle
V^{0}V^{0}\rangle$ and $\langle A^{0}A^{0}\rangle$. The channel $\langle
A^{0}A^{0}\rangle$ for $F_{0}$ has very large statistical errors and is not
shown.

(a) $F_{0},\langle SS\rangle,p=0$ (b) $F_{0},\langle PP\rangle,p=0$ (c)
$F_{1},\langle PP\rangle,p=0$ (d) $F_{2},\langle SS\rangle,p=0$ (e)
$F_{2},\langle PP\rangle,p=0$

Figure 21:  Lattice results for isospin amplitudes of $\langle SS\rangle$ and
$\langle PP\rangle$. The channel $\langle SS\rangle$ for $F_{1}$ has very
large statistical errors and is not shown.

Comparison of the data for the $L=40$ and $L=32$ lattices shows a mixed
situation regarding finite-volume effects. These are significant for $\langle
PP\rangle$ in all three isospin combinations, for $\langle A^{0}A^{0}\rangle$
in $F_{1}$ and for $\langle SS\rangle$ is $F_{2}$. In the latter case, even
the sign of the signal changes when going from $L=32$ to $L=40$. For all other
cases, volume effects are moderate, except for $\langle SS\rangle$ in $F_{0}$
at small $y$ and for $\langle V^{0}V^{0}\rangle$ in $F_{2}$ at large $y$.

We must now return to the doubly disconnected graph $D$, which contributes to
$F_{0}$ (but not to $F_{1}$ or $F_{2}$). In all plots shown, the contribution
from $D$ is omitted since the result would have too large errors to detect a
signal. We must therefore ask whether there is any reason to omit $D$ in
$F_{0}=C_{1}+2S_{1}+D$. Given the structure of the corresponding graphs in
figure 2, it seems plausible to assume that the ratio $D:S_{1}$ is similar in
size to the ratio $S_{1}:C_{1}$. If one is willing to make this hypothesis,
then one may justify the omission of $D$ in $F_{0}$ for cases in which the
data shows that $|S_{1}|\ll|C_{1}|$. In figure 22 we see that this condition
is satisfied for $\langle SS\rangle$ up to about $y\sim 1\operatorname{fm}$,
but not at all for $\langle PP\rangle$. For $\langle V^{0}V^{0}\rangle$ in the
$y$ region shown in figure 20(c), we find a very small ratio
$|S_{1}/C_{1}|<0.035$. The correlator for $\langle A^{0}A^{0}\rangle$ is too
noisy for $F_{0}$ even if we only consider the sum $C_{1}+2S_{1}$.

(a) $\langle SS\rangle,p=0,L=40$ (b) $\langle PP\rangle,p=0,L=40$

Figure 22:  Comparison of $F_{0}$ with $C_{1}$ and $S_{1}$ for $\langle
SS\rangle$ and $\langle PP\rangle$. The contribution from $D$ to $F_{0}$ is
omitted, as explained in the text.

###  6.2 Comparison with chiral perturbation theory

In figures 23 and 24, we compare our lattice results with the predictions of
chiral perturbation theory discussed in section 2.2. We remind the reader that
in leading-order chiral perturbation theory, it is admissible to replace the
pion decay constant $F$ in the chiral limit with its value $F_{\pi}$ at the
considered pion mass, which is the value we have extracted in (5.5). Let us
recall that predictions for the isospin amplitude $F_{0}$ are not available
for the channels $\langle V^{0}V^{0}\rangle$ and $\langle A^{0}A^{0}\rangle$
in Bruns:2015yto . Our lattice data in both channels is compatible with zero
within errors for $y>0.65\operatorname{fm}$ and not shown here.

(a) $F_{1},\langle V^{0}V^{0}\rangle,p=0,L=40$ (b) $F_{1},\langle
A^{0}A^{0}\rangle,p=0,L=40$ (c) $F_{2},\langle V^{0}V^{0}\rangle,p=0,L=40$ (d)
$F_{2},\langle A^{0}A^{0}\rangle,p=0,L=40$

Figure 23:  Comparison of our lattice results for $\langle V^{0}V^{0}\rangle$
and $\langle A^{0}A^{0}\rangle$ with chiral perturbation theory. The results
obtained from the leading-order chiral Lagrangian are denoted by “LO” and
those including resonance exchange by “LO $+$ Res”.

(a) $F_{0},\langle SS\rangle,p=0,L=40$ (b) $F_{0},\langle PP\rangle,p=0,L=40$
(c) $F_{1},\langle SS\rangle,p=0,L=40$ (d) $F_{1},\langle PP\rangle,p=0,L=40$
(e) $F_{2},\langle SS\rangle,p=0,L=40$ (f) $F_{2},\langle PP\rangle,p=0,L=40$

Figure 24:  As figure 23 but for the correlators $\langle SS\rangle$ and
$\langle PP\rangle$. The resonance exchange contribution is zero for $\langle
SS\rangle$ and $\langle PP\rangle$ in $F_{1}$.

Since chiral perturbation theory requires large distances to be valid, we
start our plots at $y\approx 0.65\operatorname{fm}\approx 0.95\mskip
1.5mum_{\pi}^{-1}\approx 4.0/(4\pi F_{\pi})$. Comparing the leading-order
chiral result with the resonance exchange contribution, we find three
qualitatively different cases: for $\langle SS\rangle$ and $\langle PP\rangle$
in $F_{1}$ the resonance terms are zero, for $\langle PP\rangle$ in $F_{0}$
and $\langle SS\rangle$ in $F_{2}$ the leading-order terms vanish, whereas in
all other channels both contributions are nonzero and there is a large
difference between them in the $y$ range under study. In the latter case, one
may doubt whether the chiral expansion is stable, given that the resonance
terms are formally subleading in the low-energy expansion. We also recall from
the discussion in Bruns:2015yto  that the resonance terms are only meant to be
an estimate of higher-order contributions, and that the resonance parameters
in the scalar and pseudoscalar channels are not well known.

With these caveats in mind, we now compare the chiral predictions with our
lattice data. As we see in figure 23, the agreement between lattice data and
theoretical results is quite good for $\langle V^{0}V^{0}\rangle$ in $F_{1}$
and $F_{2}$ and for $\langle A^{0}A^{0}\rangle$ in $F_{2}$. These are the
cases for which the volume effects seen in figure 20 are moderate. The
agreement for $\langle A^{0}A^{0}\rangle$ in $F_{1}$ is still fair, and we
recall that in this channel the volume dependence is somewhat larger. For the
vector and axial currents, we thus obtain a rather satisfactory picture.

The situation is quite different for the scalar and pseudoscalar currents. For
$\langle SS\rangle$ in $F_{0}$, we find no large volume effects on the lattice
and have argued that the neglect of graph $D$ should not have dramatic
effects, at least in the lower $y$ range. The comparison with the chiral
prediction including resonance exchange is very bad in this case. For the four
channels shown in figure 24(b), (d), (e) and (f), there is also a strong
disagreement between lattice data and theory. We recall from figure 21 that in
these cases the lattice results for $L=40$ and $L=32$ point to significant
finite volume effects. For $\langle SS\rangle$ in $F_{1}$, our lattice signal
is consistent with zero within errors, and the chiral prediction is zero as
well at the order considered in Bruns:2015yto . In the sense that $\langle
SS\rangle$ in $F_{1}$ is small (compared for instance with $\langle PP\rangle$
in the same isospin amplitude) we can state agreement between our lattice
computation and chiral perturbation theory in this one case.

##  7 Summary

This paper contains a detailed study of two-current correlation functions in
the pion in lattice QCD. We use two gauge ensembles with a pion mass of
$m_{\pi}\approx 300\operatorname{MeV}$, a lattice spacing of $a\approx
0.07\operatorname{fm}$ and spatial extensions of $L=32$ and $L=40$,
respectively. Additional results for heavier quarks (strange or charm) are
obtained in a partially quenched setup.

We derive several general theoretical results concerning symmetry properties
and the connection between lattice graphs and physical amplitudes. In the
remainder of our study, we focus on correlation functions of two currents
$V^{0}$, $A^{0}$, $S$ or $P$ in a pion at rest. In our lattice simulations, we
make extensive use of stochastic sources, which allows us to obtain
satisfactory (and in many cases excellent) signals for all lattice
contractions except for the double disconnected graph $D$, where vacuum
subtractions entail large cancellations and a loss of the signal in
statistical noise.

We study a number of lattice artefacts in our data. The comparison of results
for two different time differences between the pion source and sink gives no
indication for large contributions from excited states. Comparing the Lorentz
invariant correlation functions $\langle SS\rangle$ and $\langle PP\rangle$
for different pion momenta, we find that our lattice results are fully
consistent with boost invariance. We observe several types of anisotropic
behaviour, indicating a loss of rotational invariance due to either the
periodic boundary conditions or to discretisation effects. We find that these
effects are significantly reduced by selecting distances $\vec{\mskip
0.5muy}\mskip 1.5mu$ between the two currents that are close to the spatial
diagonals of the lattice. Finally, the comparison of data for two lattices
shows moderate volume effects in several channels, but large ones in others.

In the connected graph $C_{1}$, the “valence quark” and “valence antiquark” in
the pion are probed individually by the two currents. One might naively expect
this contraction to be dominant. An important result of our study is that this
is not true: we find that several lattice contractions are important for
$m_{\pi}\approx 300\operatorname{MeV}$. This concerns especially the connected
graph $C_{2}$ with two current insertions on the same quark line, but in the
case of $\langle SS\rangle$ and $\langle PP\rangle$ also the disconnected
graph $S_{1}$ and the annihilation contribution $A$. For heavier quark masses,
the importance of these graphs is reduced, although $C_{2}$ remains prominent
in $\langle A^{0}A^{0}\rangle$ and $\langle PP\rangle$ even for charm quarks.
Of course, not every contraction contributes to every physical matrix element,
as shown in (3.2).

For the connected graph $C_{1}$, we test the hypothesis that the two-current
density can be computed in terms of the single-current density, assuming the
absence of correlation effects between the quark and antiquark in the pion. We
find that this hypothesis clearly fails in a large range of current
separations $\vec{\mskip 0.5muy}\mskip 1.5mu$, both for the vector and for the
scalar current.

We also extract rms radii $r_{VV}$ and $r_{AA}$ of the correlator
$C_{1}(\vec{\mskip 0.5muy}\mskip 1.5mu)$ for the vector and for the axial
current. We find that a Fourier transform w.r.t. $\vec{\mskip 0.5muy}\mskip
1.5mu$ mitigates finite-size effects. The radius $r_{VV}$ shows a clear
decrease with the quark mass. Interestingly, it turns out that $r_{AA}<r_{VV}$
for light and strange quarks.

Combining the different lattice graphs to physical amplitudes in an isospin
basis, we can compare our lattice results with the computation in chiral
perturbation theory performed in Bruns:2015yto , which includes the
contributions from the leading-order chiral Lagrangian and an estimate of
higher-order contributions using resonance exchange graphs. We find rather
good agreement between the lattice data and the chiral calculation for vector
and axial currents, whereas the comparison for scalar and pseudoscalar
currents is poor in five out of six channels. Since the exchanged resonances
are $\rho$ and $a_{1}$ in the former case and $\sigma$, $a_{0}$, $\eta$ in the
latter, our findings are consistent with the hypothesis that the exchange of
the lowest-mass vector and axial vector mesons often gives a good estimate of
higher orders in the chiral expansion, whereas the same does not necessarily
hold for the exchange of spin-zero mesons. This is in line with the
conclusions drawn in Donoghue:1988ed ; Ecker:1988te  and the general success
of models based on vector meson dominance.

The lattice methods presented in this work are suitable for the study of two-
current correlators that can be related to double parton distributions in the
pion, thus providing a connection with an active field of research in collider
physics. This will be presented in a forthcoming paper.

## Acknowledgements

We gratefully acknowledge input from Sara Collins and Philipp Wein and
discussions with Michael Engelhardt. 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 discrete
Fourier transforms in section 5.4 were computed using the package FFTW3
Frigo:2005 , and the graphs in figures 1 and 2 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.

## References

  * (1) K. Barad, M. Ogilvie and C. Rebbi, _Quark - Anti-quark Charge Distributions and Confinement_ , _Phys. Lett._ 143B (1984) 222.
  * (2) K. Barad, M. Ogilvie and C. Rebbi, _Quark-anti-quark 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.
  * (8) M. Burkardt, J. M. Grandy and J. W. Negele, _Calculation and interpretation of hadron correlation functions in lattice QCD_ , _Annals Phys._ 238 (1995) 441 [hep-lat/9406009].
  * (9) C. Alexandrou, P. de Forcrand and A. Tsapalis, _Probing hadron wave functions in lattice QCD_ , _Phys. Rev._ D66 (2002) 094503 [hep-lat/0206026].
  * (10) C. Alexandrou, P. de Forcrand and A. Tsapalis, _The Matter and the pseudoscalar densities in lattice QCD_ , _Phys. Rev._ D68 (2003) 074504 [hep-lat/0307009].
  * (11) C. Alexandrou and G. Koutsou, _A Study of Hadron Deformation in Lattice QCD_ , _Phys. Rev._ D78 (2008) 094506 [0809.2056].
  * (12) P. C. Bruns, _Pion correlation functions in position space from chiral perturbation theory with resonance exchange_ , 1512.06650.
  * (13) M. Diehl, D. Ostermeier and A. Schäfer, _Elements of a theory for multiparton interactions in QCD_ , _JHEP_ 03 (2012) 089 [1111.0910].
  * (14) P. Hasenfratz and M. Bissegger, _CP, T and CPT in the non-perturbative formulation of chiral gauge theories_ , _Phys. Lett._ B613 (2005) 57 [hep-lat/0501010].
  * (15) G. S. Bali, S. Collins, B. Gläßle, M. Göckeler, J. Najjar, R. H. Rödl et al., _Nucleon isovector couplings from $N_{f}=2$ lattice QCD_, _Phys. Rev._ D91 (2015) 054501 [1412.7336].
  * (16) G. S. Bali, S. Collins, B. Gläßle, M. Göckeler, J. Najjar, R. H. Rödl et al., _The moment $\langle x\rangle_{u-d}$ of the nucleon from $N_{f}=2$ lattice QCD down to nearly physical quark masses_, _Phys. Rev._ D90 (2014) 074510 [1408.6850].
  * (17) RQCD collaboration, G. S. Bali, S. Collins, D. Richtmann, A. Schäfer, W. Söldner and A. Sternbeck, _Direct determinations of the nucleon and pion $\sigma$ terms at nearly physical quark masses_, _Phys. Rev._ D93 (2016) 094504 [1603.00827].