Two-current correlations in the pion on the lattice

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.


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.
Hadronic matrix elements of one or two currents can be calculated in lattice QCD. After the pioneering work [1] on two-current correlators, detailed computations for the vector, scalar and pseudoscalar currents were presented in [2] and extended studies for the vector current in [3]. 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 [1][2][3] 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. A leading-order calculation, supplemented by an estimate of higher orders using resonance exchange graphs, has been performed in [4], 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 [5], can be connected with double parton distributions, which are necessary to compute coherent hard scattering on two spatially separated partons inside a hadron. 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.

Correlation functions of two currents
The object of our study are correlation functions of two currents in a pion, where the pion charges k, k = +, −, 0 may differ between the bra and ket state, whereas the four-momentum p is the same for both. The matrix elements (2.1) are understood to be fully connected, with disconnected contributions like π|π · 0| O i (y)O j (0) |0 or π| O i (y) |0 · 0| O j (0) |π removed. The currents we consider are O qq i (y) =q(y)Γ i q (y) , (2.2) where q and q are u or d quark fields. The full set of Dirac matrices Γ i corresponds to the scalar, pseudoscalar, vector, axial and tensor currents: S qq =q q , P qq = iq γ 5 q , V µ qq =qγ µ q , A µ qq =qγ µ γ 5 q , T µν qq =qσ µν q . (2.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 (2.1) as V 0 V 0 , A 0 A 0 , SS , P P , (2.4) respectively. For later use, we shall however consider the full set of currents (2.2) in the theoretical exposition below. The matrix elements SS and P P are boost invariant and thus depend only on the products y 2 and yp of four-vectors, given that p 2 = m 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 y 2 , the matrix elements with scalar and pseudoscalar currents are independent of the pion momentum when y · p = 0, a condition that can be achieved for both zero and nonzero p.
The operators (2.2) transform under charge conjugation (C) and under the combination of parity and time reversal (P T ) as We also define the products η ij C = η i C η j C and η ij P T = η i P T η j P T .

Isospin decomposition and constraints
The matrix elements (2.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, with the Pauli matrices τ a and the isodoublet Q = (u, d) of quark fields. One then has Expressing the pion states in the isospin basis, we have We can now decompose 0) |π c = δ ab δ cd F 1 (y) + δ ac δ bd + δ ad δ bc F 2 (y) + δ ac δ bd − δ ad δ bc iF 3 (y) , 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.12), 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 (2.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 (2.14) From (2.12) one also finds that the matrix elements with two isosinglet or with two isotriplet operators are even under the simultaneous exchange 1 Figure 1. Graphs for the matrix elements SS , P P , V 0 V 0 , A 0 A 0 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 SS , V 0 V 0 and to graphs (b) and (c) for P P , A 0 A 0 . The resonance exchange graphs (d) and (e) can contribute to any of the four matrix elements, depending on the quantum numbers of the resonance.
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 (2.15) and multiplying with the product η ij C of intrinsic C parities of the two operators. One readily sees that G n = 0 if η ij C = +1 and that F n = 0 if η ij C = −1. The charge conjugation constraints give relations between matrix elements in the quark flavor basis. For η ij C = +1 the vanishing of G 1 and G 2 implies in particular and thus for k = +, −, 0.

Predictions from chiral perturbation theory
Consider pions at rest and take the matrix elements SS , P P , V 0 V 0 and A 0 A 0 at y 0 = 0 and large | y |. 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πF ∼ 1 GeV, where F is the pion decay constant. In position space, one should hence require | y | 0.2 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 [4]. 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 fourmomentum (thus setting m π 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 ρ, a 1 , a 0 , η and σ.
The corresponding results are given in sections II.A and III.B of [4]. 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 [4] 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 µ and A a µ in [4] have an extra factor of 1/2 compared with our currents (2.9). Vector and axial currents were only considered for the isovector case in [4], so that no predictions are available for F 0 in the V 0 V 0 and A 0 A 0 channels. Under the conditions stated above, the isotriplet amplitude F 3 is found to vanish for all channels SS , P P , For the remaining matrix elements, we obtain the following results from [4], writing y = | y | for simplicity. The leading order chiral Lagrangian gives rise to the amplitudes 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 [4]. The resonance exchange contributions read with signs s ρ = 1, s a 1 = −1. In the numerical evaluation of section 6.2, we will take the resonance parameters estimated in [4]: For the parameters in the chiral sector we will take m π = 300 MeV, F = 100 MeV and B = 2.4 GeV, which are rounded values of what we extract from our lattice simulations with L = 40, see (5.11).

Lattice techniques
The lattice computation of the matrix elements (2.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.

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 ij 1 (y), C ij 2 (y), etc. denote contributions to the physical matrix elements (2.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 τ (i.e. y 0 = 0), and the ratio of four-point and two-point functions is evaluated at a plateau in τ . 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 ij 1 (y) etc. on the pion momentum p. From translation invariance one readily finds Lattice contractions for the two-current correlation function (2.1) in a pion. The dependence on the pion momentum p is not indicated for brevity. The product η ij C of charge conjugation parities of the two operators is defined below (2.7).

Charge conjugation invariance gives
If η ij C = −1 one furthermore obtains M ij (y) = 0 for M = A, S 2 , D. From P T invariance one gets where for the formulation of time reversal in the Euclidean path integral we refer to [6]. 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 Combining (3.3) with (3.5) we find that C 1 , S 1 , S 2 and D are real valued, whereas combining (3.4) with (3.5) gives

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 we have for η ij 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 where F 0 , F 1 and F 2 are the isospin amplitudes defined in (2.12). These combinations are real valued owing to (3.6), as the isospin amplitudes must be. The imaginary combination iF 3 of matrix elements can e.g. be isolated by In the present work, we only study the real valued combinations in (3.9).
For η ij C = −1 we have three nonzero lattice contractions, C 1 , C 2 and S 1 , and two independent matrix elements, which may be taken as 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 where the quark content of the D + s meson is cs. 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 (3.12) within the partially quenched approximation.

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 [7,8], there are eleven standard ensembles with pion masses down to 150 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.
We will also study the dependence of the pion matrix elements ( Table 1. Details of the ensembles used in this analysis. N full is the total number of available gauge configurations and N used the number of configurations used in the present study. N sm indicates the number of Wuppertal smearing iterations in the pion source. The error on the pion mass combines statistical and systematic errors, see [7]. The time difference between pion source and sink is t = 15a or t = 32a for both ensembles.
where we also give the values of the bare quark masses m q in lattice units. Here "light quarks" refers to the κ value of ensemble V, whereas the other two values correspond to the physical strange and charm quark masses, as determined in [9] and [10] by tuning the pseudoscalar ground state mass to 685.8 MeV in the first case and the spin-averaged S-wave charmonium mass to 3068.5 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 (2.1) on a lattice in Euclidean space-time, we need the four-point correlators where Π(x) is a pion interpolator. We use interpolators with Dirac structure γ 5 . Pion matrix elements are extracted by calculating the following ratio at time slices where excited states should be suppressed: Here V = L 3 a 3 is the spatial volume and the usual pion two-point function. In analogy to (3.14) and (3.15), we use an appropriate ratio of the three-point function C i, p 3pt (τ, t) and C p 2pt (t) to extract the matrix elements π(p )| O i (0) |π(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 (3.15) is computed using the continuum dispersion relation E p = p 2 + m 2 π , which is well satisfied for the momenta investigated here. The value of m π is obtained from an exponential fit of the two-point function, which gives 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.17) agree reasonably well with the pion masses given in table 1 for light quarks and quoted after (3.13) for strange quarks. The source-sink distance is fixed to t = 15a ≈ 1.07 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 (3.15) by fitting or averaging over the plateau around τ = t/2. For the contractions C 1 and A, we measure the τ dependence of C 4pt and fit to a plateau in the τ ranges specified in (4.1). For t/a = 15 we have S 2 data for τ /a = 7 and 8, which we average, whereas for the remaining contractions C 2 , S 2 and D we take τ /a = 7 or 8 for each gauge configuration on a statistical basis. We recall that the plateau extraction must be symmetric around τ = 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 τ /a = 16.
Details on the contractions. We now give details on the calculation of the different lattice contractions for C 4pt . Let us start with a broad overview, which is pictorially represented in figure 3. Most graphs are calculated using the one-end trick [11] 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 [12] 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.
We use stochastic Z 2 ⊗ iZ 2 wall sources, i.e. a set of complex random vectors η ( ) t carrying space-time, spin and colour indices (not explicitly written here). η ( ) t is nonzero only on the time slice t, where it has components (±1 ± i)/ √ 2 that are random for each gauge configuration. Averaged over all realisations, one then has where the matrix 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 [3]. In a different context it was used in [13], where it was dubbed the "stochastic sandwich method". Including the sign of the Wick contraction (easily obtained from the number of  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. fermion loops) we explicitly have 2 is a contribution to the four-point correlator C4pt, whereas C ij 1 (y) introduced in section 3.1 is a contribution to the pion matrix element (3.15). The same holds for the other contractions.
where again V = L 3 a 3 is the spatial volume. Here and in the following, the notation . . . denotes the average over the gauge ensemble and [. . .] indicates a closed spin-colour structure. The vector Ψ ( ), p t is obtained from an inversion of the Wilson-clover Dirac operator D on the random source, is a diagonal matrix in the spatial coordinates. For better legibility, the ( x, t)-component of Ψ is written as Ψ( x, t) rather than Ψ xt in (3.19). We will use the same notation for other quantities that are vectors in space-time. Source smearing is implemented by Φ, which is a hermitian matrix acting on spatial coordinates and colour indices and consists of 400 Wuppertal smearing iterations [14]. It turned out that taking ΦE p instead of E p Φ in (3.20) greatly improves the signal for nonzero momenta. This observation contributed to introducing momentum smearing in [15]. The contractions for which we have data with nonzero p are C 1 and A on the lattice with L = 40. Specifically, C 1 was computed for all 24 nonzero momenta satisfying p 2 ≤ 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π/(La), which is equal to 437 MeV for L = 40.
As we see in (3.19), 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 η Specifically, we have Note that for symmetry reasons, the vacuum subtraction for S 1 in (3.23) is only nonzero for the scalar current S. If the one-end trick is used, then G i, p 3 is given by with the sequential propagator X We refrain from writing out the corresponding expressions for G 3,pt and X 0t,pt with point sources, given that they are identical to those used in standard computations of the pion form factor, see e.g. [16]. 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 τ . To reduce statistical noise, we make the following improvements. We consider the hopping parameter expansion of the propagator [17][18][19], writing D = (1 − H)/(2κ) with the hopping term H, and make use of the geometric series: (3.27) One needs at least hopping terms to obtain a non-vanishing contribution to the propagator from a point x to x + y on a periodic lattice of spatial size La. Hence, the first sum on the r.h.s. of (3.27) can be omitted, and we get Taking H n( y ) M instead of 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 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 = ( 0, τ ), from which the point-to-all propagator M x is obtained by inverting (3.32) Note that S x and hence also 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 ) = 1δ xx . We then have with the trace referring to spin and colour indices, and with L 1 defined in (3.24). When evaluated with the one-end trick, the two point graph reads Note that L ij 2 and L i 1 L j 1 depend on the distance y between the currents and can be nonzero for most combinations of Dirac matrices Γ i and Γ j . This makes vacuum subtractions necessary for S 2 and D in most channels.
As indicated in (3.34), we use only stochastic sources for S 2 and only point sources for D. For the pion two-point function C 2pt (t) = G 2 (t) , we compare both methods and find excellent agreement. Table 2. Numbers N st of stochastic noise vectors and numbers of sources and current insertions used for each graph in our simulations for the L = 40 lattice. 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 ins = L 3 .
To increase statistics, we repeat the calculations at N src time positions t s of the source (keeping τ and t fixed relative to t s ) and for several spatial insertion positions N 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 for the L = 40 lattice (ensemble V of table 1). These numbers are lower for the L = 32 lattice (ensemble IV) in some cases.

Renormalisation
We convert all lattice currents, which are defined in the lattice scheme at a lattice spacing a, into the MS scheme at the renormalisation scale µ = 2 GeV. The corresponding renormalisation constants depend on the gauge coupling g 2 = 6/β, where in our case β = 5.29. For currents with nonzero anomalous dimension, such as in the scalar and pseudoscalar cases, there is also a dependence on aµ. The conversion factors used for the correlator (3.15) of currents i and j read see e.g. [7]. 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 MS i have been obtained in [20] (updated in [7]) in a twostep procedure: first the currents were matched non-perturbatively from the lattice scheme to the RI'/MOM scheme and then from there to the MS scheme in continuum perturbation theory. The coefficients b i = 1 + O(g 2 ) have been computed in lattice perturbation theory [21][22][23]. Furthermore, in [24] the coefficient b S was determined non-perturbatively  .3) in the MS scheme at scale µ = 2 GeV. We also list the perturbative estimates of the coefficients b i from [7], the non-perturbative determinations (3.37) and (3.38), and the rescaled values (3.39). For the latter, we estimate an error of 34% as explained in the text.
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 MS V (1 + am q b 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 where the uncertainty is dominated by the uncertainty on Z MS V . Table 3 gives the values for Z MS i (2 GeV), as well as estimates b pert i using one-loop perturbation theory with an "improved" coupling constant (for details see (26) and (27) in [7]). In the next row, we give the non-perturbative estimates b np i from (3.37) and (3.38). We see that the non-perturbative value for b V is 24% larger than the perturbative estimate, while b S from [24] 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

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.

Plateaux and excited state contributions
For the contractions C 1 and A, we have data for the ratio (3.15) of four-point and two-point functions in a wide range of the current insertion time τ 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 τ = 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 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 2pt (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 4pt ( y , τ, t)/C 2pt (t) in (3.15) with a normalisation factor N (32a) = 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 (15a) = 1 in this case.
The data in figure 4 show a clear plateau around τ = 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 6 ≤ τ /a ≤ 9 for t/a = 15 , 14 ≤ τ /a ≤ 18 and 46 ≤ τ /a ≤ 50 for t/a = 32 . For t/a = 32 we thus combine the two plateaux around τ /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 = | y |. 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.

Anisotropy effects
In the continuum, the correlation functions V 0 V 0 , A 0 A 0 , SS and P P in a pion at rest can only depend on the modulus y = | y |. Lattice regularisation breaks rotational      invariance, so that an anisotropy in the y dependence of our correlators is a clear indication of lattice artefacts.
At large distances y/a ∼ L/2 we see clear signs of anisotropy, as shown in figure 6(a) to (c). As explained in detail in [1], 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 y with directions close to space diagonals (±1, ±1, ±1) and largest if y 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 [1,3]. An analytic investigation of the same phenomenon using a low-energy effective field theory can be found in [25]. A method to correct for these periodic images was proposed in [1] and employed in [3], where simulations were performed on L = 24 lattices with a ≈ 0.077 fm. The resulting physical length La was thus 0.81 times the length of the smallest lattice used in our study.  Benefiting from the larger size of our lattices, we take a less sophisticated approach here and consider only data for distances y close to the space diagonals. To this end, we use exact lattice symmetries to map the vector y into the octant where all its components are non-negative and define ϑ( y ) as the angle between that vector and the diagonal (1, 1, 1). We then retain only data points satisfying cos ϑ( y ) ≥ 0.9 . (4.2) After this cut, we statistically average the correlators over points y 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 P P in the S 2 graph, figure 6(d), for which we have no interpretation.
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 y , 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 [26] (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 [26] is equivalent to our cut on cos ϑ( y ), given that we always have y 0 = 0. As seen in our figure 7, the cut (4.2) 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 (4.2) 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.

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.
We notice in figure 8 that for the contractions A and S 1 , the statistical errors are significantly larger for the L = 32 lattice. This is in part due to a smaller number of source time insertions N src used for these channels in our simulations.

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 (2.4), we can test that our lattice computation yields boost invariant results using the correlation functions SS and P P . 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 y 2 if we select distance vectors with p · y = 0. This comparison is shown for the contraction C 1 in figure 9 with momenta satisfying p 2 = 0, 1, 2, where p is given in multiples of 2π/(La) ≈ 437 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 p 2 = 3, which is not shown in the figure because the larger errors would obscure the clarity of the plot.
We also have finite momentum data for the contraction A, restricted in this case to p 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.9) represents a physical matrix element in a definite isospin channel. Again, the data for different momenta are in good agreement with each other.

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.

Light quarks
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.34). For all currents, graph D has a fully connected contribution G 2 L i 1 L j 1 and a vacuum subtraction term G 2 L i 1 L j 1 (plus additional subtractions involving L 1 in the case of SS ). 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 ≈ 0.43 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 fm the errors on C 2 and S 2 become too large to draw strong conclusions. For both SS and P P 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 P P .
We note that the graph C 1 is only dominant in the case of A 0 A 0 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.

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.13) and (3.17). 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.
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. In the case of P P , we observe that C 1 crosses zero (seen as a sharp dip in the logarithmic plot). This happens around 0.95 fm for strange quarks and around 0.35 fm for charm.
Not shown in our plots is the annihilation contribution for strange quarks. It is very small compared with C 1 , except in the SS channel, where the ratio −A/C 1 is about 40% at y ∼ 0.5 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 A 0 A 0 and P P . 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 A 0 A 0 is small compared to V 0 V 0 and P P relatively small compared to SS . By contrast, for light quarks A 0 A 0 is comparable in size to V 0 V 0 , and P P is comparable in size to SS . Our results for charm quarks confirm that in the heavy-quark limit, matrix elements of the spin dependent currents A µ and P are suppressed compared with their spin independent counterparts V µ and S. Seen from this angle, the contraction C 1 is indeed dominant for charm because it is the most prominent contribution to V 0 V 0 and SS .

Test of factorisation
The two-current matrix elements V 0 V 0 and SS quantify the distribution of vector or scalar charges at two separate points in the pion. By contrast, the corresponding singlecurrent matrix elements π(p )| V 0 |π(p) and π(p )| S |π(p) give the vector or scalar charge distribution at one point in the pion if one performs a Fourier transform from the momentum transfer p − p to three-dimensional space. 3 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 [1], 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 [1,5]. In momentum space, we have where E 2 q = m 2 π + q 2 and it is understood that p = 0. The question mark above the equality sign indicates an assumption, which will be tested in the following. In the last step of (5.1) we used the relation π + | O dd i |π + = η i C π + | O uu i |π + with the C parity η i C defined in (2.6). Defining vector and scalar form factors as 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 q 2 m 2 π the argument of the form factors in (5.3) becomes q 2 , as appropriate for a non-relativistic treatment. For our study with m π ≈ 300 MeV, this limit is however not relevant.
It follows from (3.8) that the two-current correlators in (5.1) 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 π + | O uc i (y)O sd i (0) |D + s and inserts a full set of intermediate states with quark content cd. Retaining only the ground state term |D + 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 q = 0. Let us first test the factorisation hypothesis (5.3) at q = 0, where it reads For the vector current we used the normalisation condition F V (0) = 1. We evaluate the integrals on the l.h.s. of (5.4) as discrete sums over all lattice sites. 4 The result for the vector correlator is  Table 4. Parameters obtained in fits of our form factor data to (5.8).
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 2m π a 3 y SS ( y ) 1/2 = 3.14(4) GeV (5.6) for light quarks, which is to be compared with the value F S (0) = 2.2 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 (y) is a conserved current in QCD, so that the associated Noether charge ρ q = d 3 y V µ qq (y) corresponds to a good quantum number. For a positive pion one has π + | ρ u = π + | and thus obtains 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 y = 0. 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 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 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. The extracted form factor data and the results of the fits are shown in figure 14, and the fitted parameters are given in table 4. We see that the fits describe the F V data very  Test of the factorisation hypothesis for the two-current correlators V 0 V 0 and SS . 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 (5.8) and table 4. 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.
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 MeV.

Subtraction term for the annihilation graph
The vacuum subtraction term for the annihilation graph A in (3.19) involves the same three-point functions B t 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 for a pion at rest. Using ∂ µ A µ ud (x) = 2m q P ud (x), where m q is the average of u and d quark masses, we have for a pion at rest 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 π for m π ≈ 300 MeV agree within a few percent with chiral perturbation theory at NLO, see figure 10 in [28]. From equations (46), (65) and table 14 in the FLAG review [29], we obtain B = Σ/F 2 ≈ 2.6 to 2.7 GeV for n F = 2 quark flavours in the chiral limit. This agrees reasonably well with our values in (5.11), given that we use lowest-order chiral perturbation theory for their extraction and have not attempted to correct our data for lattice artefacts.

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.12). We compare them with the predictions of chiral perturbation theory that were computed in [4] and recalled in section 2.2 of the present paper. We consider the case of light quarks throughout.

Isospin amplitudes
In figures 16 and 17, we show our results for the isospin amplitudes in all channels except P P for F 0 and SS 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. Comparison of the data for the L = 40 and L = 32 lattices shows a mixed situation regarding finite-volume effects. These are significant for P P in all three isospin combinations, for A 0 A 0 in F 1 and for SS 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 SS in F 0 at small y and for V 0 V 0 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 | |C 1 |. In figure 18 we see that this condition is satisfied for SS up to about y ∼ 1 fm, but not at all for P P . For V 0 V 0 in the y region shown in figure 16(c), we find a very small ratio |S 1 /C 1 | < 0.035. The correlator for A 0 A 0 is too noisy for F 0 even if we only consider the sum C 1 + 2S 1 .

Comparison with chiral perturbation theory
In figures 19 and 20, 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 π at the considered pion mass, which is the value we have extracted in (5.11). Let us recall that predictions for the isospin amplitude F 0 are not available for the channels V 0 V 0 and A 0 A 0 in [4]. Our lattice data in both channels is compatible with zero within errors for y > 0.65 fm and not shown here.
Since chiral perturbation theory requires large distances to be valid, we start our plots at y ≈ 0.65 fm ≈ 0.95m −1 π ≈ 4.0/(4πF π ). Comparing the leading-order chiral result with the resonance exchange contribution, we find three qualitatively different cases: for SS and P P in F 1 the resonance terms are zero, for P P in F 0 and SS in F 2 the leadingorder 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. We also recall from the discussion in [4] that the resonance contributions are meant as an estimate for the full next-to-leading order chiral results, 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 19, the agreement between lattice data and theoretical results is quite good for V 0 V 0 in F 1 and F 2 and for A 0 A 0 in F 2 . These are the cases for which the volume effects seen in figure 16 are moderate. The agreement for A 0 A 0 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 SS 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 As figure 19 but for the correlators SS and P P . The resonance exchange contribution is zero for SS and P P in F 1 .
four channels shown in figure 20(b), (d), (e) and (f), there is also a strong disagreement between lattice data and theory. We recall from figure 17 that in these cases the lattice results for L = 40 and L = 32 point to significant finite volume effects. For SS 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 [4]. In the sense that SS in F 1 is small (compared for instance with P P in the same isospin amplitude) we can state agreement between our lattice computation and chiral perturbation theory in this one case.

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 π ≈ 300 MeV, a lattice spacing of a ≈ 0.07 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 SS and P P 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 y 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 π ≈ 300 MeV. This concerns especially the connected graph C 2 with two current insertions on the same quark line, but in the case of SS and P P 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 A 0 A 0 and P P even for charm quarks. Of course, not every contraction contributes to every physical matrix element, as shown in (3.8).
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 y , both for the vector and for the scalar current.
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 [4], 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 ρ and a 1 in the former case and σ, a 0 , η 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 [30,31] 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.