Double parton distributions in the pion from lattice QCD
Gunnar S. Bali, Luca Castagnini, Markus Diehl, Jonathan R. Gaunt, Benjamin Gläßle, Andreas Schäfer, Christian Zimmermann
PPrepared for submission to JHEP
DESY 20-098, CERN-TH-2020-086
Double parton distributions in the pion from latticeQCD
Gunnar S. Bali a,b
Luca Castagnini a Markus Diehl c,d
Jonathan R. Gaunt e BenjaminGl¨aßle f Andreas Sch¨afer a and Christian Zimmermann a a Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany b Department of Theoretical Physics, Tata Institute of Fundamental Research, Homi Bhabha Road,Mumbai 400005, India c Fachbereich Physik, University of Hamburg, 22761 Hamburg, Germany d Deutsches Elektronen-Synchroton DESY, 22607 Hamburg, Germany e CERN Theory Division, 1211 Geneva 23, Switzerland f Zentrum f¨ur Datenverarbeitung, Universit¨at T¨ubingen, W¨achterstr. 76, 72074 T¨ubingen, Germany
Abstract:
We perform a lattice study of double parton distributions in the pion, using therelationship between their Mellin moments and pion matrix elements of two local currents. Agood statistical signal is obtained for almost all relevant Wick contractions. We investigatecorrelations in the spatial distribution of two partons in the pion, as well as correlationsinvolving the parton polarisation. The patterns we observe depend significantly on the quarkmass. We investigate the assumption that double parton distributions approximately factoriseinto a convolution of single parton distributions. a r X i v : . [ h e p - l a t ] J un ontents py dependence 405.2 Fitting the data 425.3 Mellin moments of DPDs 445.4 Factorisation hypothesis for Mellin moments 49 – 1 – Introduction
Matrix elements of currents in a hadron offer a variety of ways to quantify and study hadronstructure. In particular, information about correlations inside the hadron can be obtainedfrom the matrix elements of two currents that are separated by a space-like distance. Suchmatrix elements can be calculated in lattice QCD, and there has been considerable activityin this area over the years [1–11]. These studies address a broad range of physics questions,such as confinement [1, 2], the size of hadrons [3–5, 8], density correlations [6], comparisonwith quark models [7], or the non-spherical shape of hadrons with spin 1 or larger [9–11].We continued this line of investigation in a recent paper [12]. We performed a latticecomputation of the matrix elements of two scalar, pseudoscalar, vector, or axial vector cur-rents in the pion and compared our results with predictions of chiral perturbation theory. Forthe first time, we computed all Wick contractions that contribute to these matrix elements,whilst earlier work had focused on the case in which the two currents are inserted on differentquark lines between the hadron source and sink operators (see graph C in figure 4). Weobtained signals with a good statistical accuracy for almost all contractions and were thusable to study their relative importance. Our results were compared with different modelsin [13, 14].Extending our work in [12], we will in the present paper use two-current matrix el-ements from the lattice to obtain information about double parton distributions (DPDs).DPDs describe the correlated distribution of two partons inside a hadron and appear in thecross sections for double parton scattering, which occurs when there are two separate hard-scattering processes in a single hadron-hadron collision. The study of this mechanism has along history in collider physics, from early theoretical papers such as [15–21] to the detailedinvestigation of QCD dynamics and factorisation that started about ten years ago [22–33].After early experimental studies [34, 35], a multitude of double parton scattering processeshas been measured at the Tevatron and the LHC, see [36–40] and references therein. Somefinal states produced by double parton scattering are of particular interest because they area background to search channels for new physics. A prominent example are like-sign gaugeboson pairs W + W + and W − W − [40–45], the decay of which can yield like-sign lepton pairs.A wealth of further information about double parton scattering can be found in the mono-graph [46].Double parton distributions remain poorly known, and their extraction from experimentaldata is considerably more difficult than the extraction of single parton distributions (PDFs).It is therefore important to have as much theoretical guidance as possible about the propertiesand behaviour of DPDs. Apart from approaches that focus on fulfilling theoretical constraints[47–50], there exists a large number of model calculations for the DPDs of the nucleon [51–59]and a smaller number for those of the pion [60–62].A relation between the Mellin moments of DPDs and two-current matrix elements thatcan be computed on the lattice was written down in [23, 27]. This generalises the relationbetween matrix elements of one current and the Mellin moments of PDFs, which has been– 2 –xtensively exploited in lattice studies, as reviewed for instance in [63–65]. Whilst knowledgeof a few Mellin moments is insufficient for reconstructing the full DPDs, it allows one to in-vestigate crucial features of these functions, such as their dependence on the distance betweenthe two partons and on the parton polarisation. In the present paper, we pursue this ideafor the DPDs of the pion, focusing on their lowest Mellin moments. We use the same latticedata as in our study [12]. Corresponding work on the DPDs of the nucleon is in progress,and preliminary results have been presented in [66].This paper is organised as follows. In section 2, we recapitulate some basics about DPDsand then elaborate on the relation between their Mellin moments and the two-current matrixelements we compute on the lattice. This will in particular lead us to introduce the conceptof skewed DPDs. In section 3, we describe the main elements of our lattice simulations (a fullaccount is given in [12]) and investigate several lattice artefacts that are present in our data.Our results for zero pion momentum are presented and discussed in section 4. In section 5,we develop a parametrisation of the data for both zero and nonzero pion momenta, whichwill allow us to reconstruct the Mellin moments of pion DPDs, albeit in a model-dependentfashion. Our main findings are summarised in section 6. To begin with, we recall some basics about double parton distributions. An extended intro-duction to the subject can be found in [67].Factorisation for a double parton scattering process means that its cross section is givenin terms of hard-scattering cross sections at parton level and double parton distributions foreach of the colliding hadrons. For pair production of colourless particles, such as Z , W orHiggs bosons, this factorisation can be proven rigorously. A DPD gives the joint probabilityfor finding in a hadron two partons with longitudinal momentum fractions x and x at atransverse distance y from each other. The distributions for quarks and antiquarks are definedby operator matrix elements as F a a ( x , x , y ) = 2 p + (cid:90) dy − (cid:90) dz − π dz − π e i ( x z − + x z − ) p + × (cid:104) h ( p ) | O a ( y, z ) O a (0 , z ) | h ( p ) (cid:105) . (2.1)We use light-cone coordinates v ± = ( v ± v ) / √ v = ( v , v ) for any four-vector v µ . The definition (2.1) refers to a reference frame in whichthe transverse hadron momentum is zero, p = . In a frame where the hadron moves fast inthe positive z direction, x and x can be interpreted as longitudinal momentum fractions.The hadron state is denoted by h ( p ), and it is understood that an average over its polarisationis taken on the r.h.s. of (2.1) if the hadron has nonzero spin. Unless specified otherwise, theexpressions of the present section hold both for a pion and for the nucleon (and in fact forany unpolarised hadron or nucleus). – 3 –he matrix element in (2.1) involves the same twist-two operators that appear in thedefinition of ordinary PDFs. For quarks, one has O a ( y, z ) = ¯ q (cid:0) y − z (cid:1) Γ a q (cid:0) y + z (cid:1)(cid:12)(cid:12)(cid:12) z + = y + =0 , z = (2.2)with spin projectionsΓ q = γ + , Γ ∆ q = γ + γ , Γ jδq = iσ j + γ ( j = 1 , . (2.3)The analogous expressions for antiquarks can e.g. be found in [27, section 2.2]. The form (2.2)holds in light-cone gauge A + = 0, whereas in other gauges a Wilson line is to be insertedbetween the fields. Since the two fields in (2.2) have light-like separation from each other,their product requires renormalisation. This results in a scale dependence of the two operatorsand of the DPD in (2.1), which we do not indicate for the sake of brevity.Lorentz invariance implies that one can write F q q ( x , x , y ) = f q q ( x , x , y ) ,F ∆ q ∆ q ( x , x , y ) = f ∆ q ∆ q ( x , x , y ) ,F j δq q ( x , x , y ) = (cid:15) j k y k mf δq q ( x , x , y ) ,F j q δq ( x , x , y ) = (cid:15) j k y k mf q δq ( x , x , y ) ,F j j δq δq ( x , x , y ) = δ j j f δq δq ( x , x , y ) + (cid:0) y j y j − δ j j y (cid:1) m f tδq δq ( x , x , y ) (2.4)with y = y µ y µ = − y . Due to parity invariance, one has F q ∆ q = F ∆ q q = 0, and timereversal invariance implies F δq ∆ q = F ∆ q δq = 0. The hadron mass m has been introducedon the r.h.s. of (2.4) so that all scalar functions f have the dimension of an inverse area. Theoperator O jδq is a vector whose direction gives the transverse quark spin direction, and (cid:15) jk isthe two-dimensional antisymmetric tensor with (cid:15) = +1. The density interpretation of thedifferent distributions in (2.4) is then as follows: • The unpolarised distribution f q q gives the probability density to find two quarks withmomentum fractions x and x at a transverse distance y , regardless of their polarisa-tion. • f ∆ q ∆ q is the density for finding two quarks with their longitudinal polarisations alignedminus the density for finding them with their longitudinal polarisations anti-aligned. • f δq δq is the analogue of f ∆ q ∆ q for transverse quark polarisations. • f δq q describes a correlation between the transverse polarisation of the quark q andthe distance y of that quark from the unpolarised quark q . In f q δq , the first quark isunpolarised and the second quark has transverse polarisation. • f tδq δq describes a correlation between the transverse polarisations of the two quarksand their transverse distance y . – 4 –ecompositions of the same form as (2.4) can be given for the cases where one replaces oneor both of the quarks by an antiquark, with the same physical interpretation as given abovefor two quarks.Note that the polarisation dependence of DPDs is not only interesting from the point ofview of hadron structure, but can have measurable implications on double parton scattering,as was for instance shown in [27, 44, 45, 68]. Lattice calculations can give information aboutthe strength of the different spin correlations we just discussed.We note that cross sections for double parton scattering involve the product of two DPDsintegrated over the interparton distance, (cid:90) d y F a a ( x , x , y ) F b b ( x (cid:48) , x (cid:48) , y ) . (2.5)The dependence of DPDs on y can hence not be directly inferred from experimental observ-ables. If y is small, one can use perturbation theory to compute F a a in terms of PDFs andsplitting functions [27, 69]. By contrast, for large distances the y dependence is fully non-perturbative. Lattice studies can give information about this dependence, whose knowledgeis crucial for computing double parton scattering cross sections.Both unpolarised and polarised DPDs can exhibit correlations in their dependence on x , x and y . We cannot address this aspect in our present study, because the matrix elementswe compute are related to the lowest Mellin moments of DPDs, i.e. their integrals over both x and x . In principle, one could investigate higher Mellin moments, i.e. integrals weightedwith powers of x and x . This would require extending the set of currents in (2.9) to currentsthat involve covariant derivatives and is beyond the scope of the present work.Phenomenological analyses often make the assumption that in unpolarised DPDs the twopartons are independent of each other. This gives the relation F a a ( x , x , y ) ? = (cid:90) d b f a ( x , b + y ) f a ( x , b ) , (2.6)where f a ( x, b ) is an unpolarised impact parameter dependent single parton distribution. Thequestion mark above the equal sign in (2.6) indicates that this is a hypothesis. Our latticestudy allows us to test this indirectly in two different ways, as discussed in sections 2.4, 4.4and 5.4.A related but different simplifying assumption is that unpolarised DPDs can be writtenas F a a ( x , x , y ) ? = f a ( x ) f a ( x ) G ( y ) , (2.7)where f a ( x ) denotes a standard PDF and G ( y ) is a factor describing the dependence on thetransverse parton distance. This assumption leads to the so-called “pocket formula”, whichexpresses double parton scattering cross sections in terms of the cross sections for each singlescattering and a universal factor σ − = (cid:82) d y [ G ( y )] . Whilst our study cannot address thefactorisation between the x , x and y dependence assumed in (2.7), we can investigate theassumption that the y dependence is the same for all parton combinations ( a , a ) in a givenhadron. We will do this in section 4.5. – 5 – .2 Matrix elements of local currents The matrix element (2.1) involves fields at light-like distances and is hence not suitable fordirect evaluation on a Euclidean lattice. What we can study in Euclidean space-time are thematrix elements M µ ··· µ ··· q q ,i i ( p, y ) = (cid:104) h ( p ) | J µ ··· q ,i ( y ) J µ ··· q ,i (0) | h ( p ) (cid:105) , (2.8)where as in (2.1) a polarisation average is understood if the hadron h carries spin. The localcurrents J µ ··· q,i we will consider here are J µq,V ( y ) = ¯ q ( y ) γ µ q ( y ) , J µq,A ( y ) = ¯ q ( y ) γ µ γ q ( y ) , J µνq,T ( y ) = ¯ q ( y ) σ µν q ( y ) . (2.9)For spacelike distances y , which we assume throughout this work, the two currents in (2.8)commute, so that one has M µ ··· µ ··· q q ,i i ( p, y ) = M µ ··· µ ··· q q ,i i ( p, − y ) . (2.10)Together with the fact that the currents in (2.9) are Hermitian, it follows that the matrixelements (2.8) are real valued.The currents transform under charge conjugation ( C ) and under the combination of parityand time reversal ( P T ) as J µ ··· q,i ( y ) → C η i C J µ ··· q,i ( y ) , J µ ··· q,i ( y ) → P T η i PT J µ ··· q,i ( − y ) (2.11)with sign factors η i C = +1 for i = A , η i C = − i = V, T (2.12)and η i PT = +1 for i = V , η i PT = − i = A, T . (2.13)The combination of a parity and time reversal transformation gives M µ ··· µ ··· q q ,i i ( p, y ) = η i PT η i PT M µ ··· µ ··· q q ,i i ( p, − y ) (2.14)and thus relates the matrix elements for y and − y . Symmetry relations for pion matrix elements.
For pion matrix elements, one hasadditional relations due to charge conjugation and isospin invariance. For η i C η i C = 1, whichis the case for all current combinations considered in our work, one has M q q ( y, p ) (cid:12)(cid:12) π + = M q q ( y, p ) (cid:12)(cid:12) π − , (2.15)– 6 –here we indicated for which hadron the matrix element is taken but for brevity omitted theLorentz indices and the labels i , i specifying the currents. Still for η i C η i C = 1, one finds M uu ( y, p ) (cid:12)(cid:12) π c = M dd ( y, p ) (cid:12)(cid:12) π c , M ud ( y, p ) (cid:12)(cid:12) π c = M du ( y, p ) (cid:12)(cid:12) π c (2.16)for c = + , − ,
0, as well as M ud ( y, p ) (cid:12)(cid:12) π + + M uu ( y, p ) (cid:12)(cid:12) π + = M ud ( y, p ) (cid:12)(cid:12) π + M uu ( y, p ) (cid:12)(cid:12) π . (2.17)A derivation of these relations can be found in [12, section 2.1]. Tensor decomposition and extraction of twist-two functions.
The matrix elementsin (2.8) are related to the lowest Mellin moments of DPDs as (cid:90) ∞−∞ dy − M ++ q q ,V V ( p, y ) (cid:12)(cid:12)(cid:12) y + =0 , p = = 2 p + I q q ( y ) , (cid:90) ∞−∞ dy − M ++ q q ,AA ( p, y ) (cid:12)(cid:12)(cid:12) y + =0 , p = = 2 p + I ∆ q ∆ q ( y ) , (cid:90) ∞−∞ dy − M k ++ q q ,T V ( p, y ) (cid:12)(cid:12)(cid:12) y + =0 , p = = 2 p + y k mI δq q ( y ) , (cid:90) ∞−∞ dy − M + k + q q ,V T ( p, y ) (cid:12)(cid:12)(cid:12) y + =0 , p = = 2 p + y k mI q δq ( y ) , (cid:90) ∞−∞ dy − M k + k + q q ,T T ( p, y ) (cid:12)(cid:12)(cid:12) y + =0 , p = = 2 p + (cid:2) δ k k I δq δq ( y ) − (cid:0) y k y k − δ k k y (cid:1) m I tδq δq ( y ) (cid:3) (2.18)with the Mellin moments given by I a a ( y ) = (cid:90) − dx (cid:90) − dx f a a ( x , x , y )= (cid:90) dx (cid:90) dx (cid:104) f a a ( x , x , y ) + η i C f ¯ a a ( x , x , y )+ η i C f a ¯ a ( x , x , y ) + η i C η i C f ¯ a ¯ a ( x , x , y ) (cid:105) . (2.19)Here i and i refer to the currents in the matrix elements on the l.h.s. of (2.18). An analogousrelation holds between I t and the lowest moment of f t . The relations (2.18) extend the well-known connection between the Mellin moments of PDFs and the matrix elements of a singlelocal current to the case of two partons.In analogy to the case of PDFs, the matrix element (2.1) defining a DPD has supportfor both positive and negative x and x , with positive x i corresponding to a parton a i andnegative x i to its antiparton ¯ a i . On the r.h.s. of (2.19), we have limited the integration regionto positive momentum fractions. Note that if a and a are quarks and if i = V or T (butnot A ), then the quark-antiquark distributions on the r.h.s. enter with a minus sign. This is– 7 –f special importance for distributions in a pion, whose valence Fock state consists of a quarkand an antiquark. Relations analogous to (2.18) exist for higher Mellin moments in x and x and involve local currents with covariant derivatives [27], as is the case for PDFs.Contrary to Γ jδq in (2.3), the tensor current J µνq,T in (2.9) is defined without γ . As aconsequence, the vector indices k and k in (2.18) do not give the transverse quark spindirection but the transverse quark spin direction rotated by +90 ◦ in the x − y plane. Thisfollows from the relation iσ j + γ = (cid:15) jk σ k + .The relations (2.18) still refer to Minkowski space, because they involve plus-components.To make contact with matrix elements evaluated in Euclidean space, we decompose the matrixelements (2.8) in terms of basis tensors and of Lorentz invariant functions A , B , C , D , E that depend on y = y µ y µ and py = p µ y µ . We write (cid:2) M µνq q ,V V ( p, y ) + M νµq q ,V V ( p, y ) (cid:3) = t µνV V,A A q q + t µνV V,B m B q q + t µνV V,C m C q q + t µνV V,D m D q q , T M µνρq q ,T V ( p, y ) = u µνρT V,A mA δq q + u µνρT V,B m B δq q , (cid:2) M µνρσq q ,T T ( p, y ) + M ρσµνq q ,T T ( p, y ) (cid:3) = u µνρσT T,A A δq δq + u µνρσT T,B m B δq δq + u µνρσT T,C m C δq δq + u µνρσT T,D m D δq δq + u µνρσT T,E m E δq δq . (2.20)For the operator combination T V , we subtract trace terms according toT u µνρ = u µνρ + (cid:0) g µρ u ναα − g νρ u µαα (cid:1) , (2.21)where it is understood that u µνρ is antisymmetric in µ and ν . The decomposition for M q q ,AA has the same form as the one for M q q ,V V , involving the same basis tensors but differentinvariant functions A ∆ q ∆ q , . . . , D ∆ q ∆ q . The decomposition for M q q ,V T is like the one for M q q ,T V with an appropriate change in the role of the Lorentz indices. In the following, wewill not discuss the combination V T any further, because it can be traded for
T V using therelation (2.10). The basis tensors are chosen as t µνV V,A = 2 p µ p ν − g µν p ,t µνV V,B = p µ y ν + p ν y µ − g µν py ,t µνV V,C = 2 y µ y ν − g µν y ,t µνV V,D = g µν ,u µνρT V,A = 2( y µ p ν − p µ y ν ) p ρ + ( g µρ y ν − g νρ y µ ) p − ( g µρ p ν − g νρ p µ ) py ,u µνρT V,B = 2( y µ p ν − p µ y ν ) y ρ + ( g µρ y ν − g νρ y µ ) py − ( g µρ p ν − g νρ p µ ) y ,u µνρσT T,A = − (cid:0) g µρ p ν p σ − g µσ p ν p ρ ) + ( g µρ g νσ − g µσ g νρ ) p − { µ ↔ ν } ,u µνρσT T,B = − y u µνρσT T,A − y µ p ν − p µ y ν )( y ρ p σ − p ρ y σ ) + ( g µρ g νσ − g µσ g νρ ) (cid:2) p y − ( py ) (cid:3) ,u µνρσT T,C = − ( g µρ p ν y σ − g µσ p ν y ρ + g µρ y ν p σ − g µσ y ν p ρ ) + ( g µρ g νσ − g µσ g νρ ) py − { µ ↔ ν } , – 8 – µνρσT T,D = − (cid:0) g µρ y ν y σ − g µσ y ν y ρ ) + ( g µρ g νσ − g µσ g νρ ) y − { µ ↔ ν } ,u µνρσT T,E = g µρ g νσ − g µσ g νρ . (2.22)The tensor components related to twist-two matrix elements can be identified from the l.h.s.of (2.18), taking into account that y + = 0 and p = in that equation. For the basis tensors,a nonzero plus-component requires the vector p on the r.h.s. of (2.22), whilst a nonzerotransverse component requires the vector y or the metric tensor. One thus finds that theinvariant functions corresponding to operators of twist two are A q q , A ∆ q ∆ q , A δq q , A δq δq and B δq δq . We will call them “twist-two functions” in the remainder of this work. All ofthem are even functions of py due to the symmetry relation (2.14).One can project out the invariant functions by multiplying the matrix elements with suit-able linear combinations of basis tensors. For the twist-two functions, the relevant projectionsread A q q = 18 N (cid:110) y ) t µνV V,A − y py t µνV V,B + (cid:2) p y + 2( py ) (cid:3) t µνV V,C (cid:111)(cid:2) M q q ,V V (cid:3) µν ,mA δq q = 316 N (cid:110) y u µνρT V,A − py u µνρT V,B (cid:111) T (cid:2) M q q ,T V (cid:3) µνρ ,A δq δq = 164 N (cid:110) y ) u µνρσT T,A − y py u µνρσT T,C + (cid:2) p y + 2( py ) (cid:3) u µνρσT T,D (cid:111) (cid:2) M q q ,T T (cid:3) µνρσ ,m B δq δq = 164 N (cid:110) u µνρσT T,B + 6 py u µνρσT T,C − p u µνρσT T,D (cid:111) (cid:2) M q q ,T T (cid:3) µνρσ (2.23)with a normalisation factor N = p y − ( py ) . (2.24)For spacelike y µ , which we are interested in, one has N <
0, so that the projections are alwayswell defined.Using (2.18) and (2.20), one can derive the relation between Mellin moments of DPDsand integrals of twist-two functions over py : I a a ( y ) = (cid:90) ∞−∞ d ( py ) A a a ( py, y ) ,I tδq δq ( y ) = (cid:90) ∞−∞ d ( py ) B δq δq ( py, y ) , (2.25)where in the first line we have all combinations of ( a , a ) that appear on the r.h.s. of (2.18).The matrix elements (2.8) can be evaluated in Euclidean space-time at y , i.e. with thetwo current operators taken at equal Euclidean time. This entails the important restriction( py ) = ( (cid:126)p(cid:126)y ) ≤ (cid:126)p (cid:126)y , (2.26)where (cid:126)v = ( v , v , v ) denotes the spatial components of a four-vector v µ . Since the range ofaccessible hadron momenta (cid:126)p in a lattice calculation is finite, the range of the variable py is– 9 –imited, and one cannot directly evaluate the integrals in (2.25). In addition, one needs datafor nonzero hadron momentum (cid:126)p to access even a finite range in py .We note that the restriction (2.26) also applies if one computes the Mellin moments oftransverse-momentum dependent single parton distributions (TMDs) on the lattice [70–72].In that case, y µ is the distance between the quark and the antiquark field in the matrixelements that define the distributions. The same holds for lattice studies of single partondistributions in x space. There has been an enormous amount of activity in this area inrecent years; we can only cite a few papers here [73–81] and refer to the recent reviews[82, 83] for an extended bibliography. Together with the restriction (2.26), the necessity to perform an integral over all py in (2.25)presents a significant complication for relating matrix elements calculated on a Euclidean lat-tice with the Mellin moments of DPDs. This prompts us to extend the theoretical frameworkin such a way that we can discuss the physical meaning of the twist-two functions A a a and B δq δq at a given value of py .To this end, we introduce skewed double parton distributions F a a ( x , x , ζ, y ) = 2 p + (cid:90) dy − e − iζy − p + (cid:90) dz − π dz − π e i ( x z − + x z − ) p + × (cid:104) h ( p ) | O a ( y, z ) O a (0 , z ) | h ( p ) (cid:105) . (2.27)Compared with the definition (2.1) of ordinary DPDs, we have an additional exponential e − iζp + y − here. As a consequence, the partons created or annihilated by the fields ¯ q and q in O a and O a have different longitudinal momentum fractions. A sketch is given in figure 1for ( a , a ) = ( u, d ) and the case where x − ζ , x + ζ , x − ζ and x + ζ are all positive.If x − ζ becomes negative, the u quark in the wave function of | h (cid:105) becomes an antiquark ¯ u with momentum fraction − x + ζ in the wave function of (cid:104) h | . Corresponding statementshold for x + ζ , x − ζ and x + ζ .For nonzero ζ , the distributions (2.27) do not appear in cross sections for double partonscattering, but they may be regarded as a rather straightforward extension of the DPDconcept. Let us take a closer look at some of their properties. The support region of thematrix element (2.27) in the momentum fraction arguments is the same as if all four partonfields were at the same transverse position. In that case, we would have a collinear twist-four distribution. The support properties of these distributions were derived in [84], and theargument given there does not depend on the transverse position arguments of the partonfields. The result given in [84] is equivalent to the interpretation of x − ζ , x + ζ , x − ζ and x + ζ as positive or negative momentum fractions, as described in the previous paragraph. The term “skewed” refers to the parton momenta here, whilst the hadron momentum is the same in thebra and ket vector of (2.27). This is different from “skewed parton distributions”, now commonly called“generalised parton distributions”, which involve two instead of four parton fields, such that there is anasymmetry both in the parton and in the hadron momenta. – 10 – − ζ x + ζ x − ζ x + ζF ud ( x , x , ζ, y ) u d ud | h i h h | Figure 1 . Graphical representation of a skewed DPD for quark flavours u and d in the hadron h .The configuration shown is for the case where all momentum fractions given at the top of the graphare positive. For nonzero ζ there are hence different regions, in which one has either 1, 2 or 3 partons in thewave function of | h (cid:105) . With the constraints that the partons in the wave function of | h (cid:105) mustcarry the same total longitudinal momentum as those in the wave function of (cid:104) h | , and thatthis cannot be larger than the longitudinal hadron momentum, one obtains the constraint − ≤ ζ ≤ x , x ) shown in figure 2. For ζ = 0 this region becomes asquare with corners (0 , ±
1) and ( ± , ζ = ± ± , ± ).Using P T symmetry, one finds that F a a ( x , x , ζ, y ) = η i PT η i PT F a a ( x , x , − ζ, − y ) , (2.29)where η i PT = +1 for an unpolarised parton and η i PT = − ζ . The symmetry property (2.29) then implies f a a ( x , x , ζ, y ) = f a a ( x , x , − ζ, y ) (2.30)and an analogous relation for f t . Mellin moments.
We define the lowest Mellin moments of skewed DPDs as I a a ( y , ζ ) = (cid:90) − dx (cid:90) − dx f a a ( x , x , ζ, y ) (2.31)and likewise for f t , where the integration region in x , x follows from figure 2. The momentsare nonzero for ζ in the interval [ − , ζ reads I a a ( y , ζ ) = (cid:90) ∞−∞ d ( py ) e − iζpy A a a ( y , py ) , (2.32)– 11 – ζ, − ζ ) (1 − ζ, ζ ) d ¯ d | u ¯ ud | du ¯ u ¯ d | ¯ du ¯ u ud ¯ d | u ¯ ud ¯ d | ¯ u ud | duu ¯ d | ¯ du ¯ u ¯ d | ¯ d ¯ u ¯ ud | d ¯ u x + x = ≤ ζ ≤ − ζ, ζ ) (1 + ζ, − ζ ) ¯ uu | ¯ dd ¯ uud | d ¯ uu ¯ d | ¯ d u | ¯ ddu ¯ u | ¯ dd ¯ u ud | duu ¯ d | ¯ du ¯ u ¯ d | ¯ d ¯ u ¯ ud | d ¯ u x + x = − ≤ ζ ≤ x x x x Figure 2 . Support region of the distribution F ud ( x , x , ζ, y ) in the momentum fraction arguments.The notation d | du ¯ u means that one has one d quark in the wave function of | h (cid:105) and du ¯ u in thewave function of (cid:104) h | . In both panels, the triangle for the region ud | du has the corners (cid:0) | ζ | , | ζ | (cid:1) , (cid:0) | ζ | , − | ζ | (cid:1) and (cid:0) − | ζ | , | ζ | (cid:1) . Notice that the parton configuration in each of the four trianglesis the same for positive and negative ζ , whereas the configuration in each of the squares is different. which can readily be inverted for the function A a a ( y , py ). In particular, one finds A a a ( y , py = 0) = 1 π (cid:90) dζ I a a ( y , ζ ) , (2.33)where we have used the symmetry relation (2.30) to reduce the integration region to positive ζ . Rather than the Mellin moment of a DPD, a twist-two function at py = 0 is thus theaverage of the Mellin moment of a skewed DPD over the skewness parameter ζ .Quantities that characterise the ζ dependence of I a a ( y , ζ ) are the even moments in ζ , (cid:104) ζ m (cid:105) a a ( y ) = (cid:82) − dζ ζ m I a a ( y , ζ ) (cid:82) − dζ I a a ( y , ζ ) = (cid:34) ( − m A a a ( y , py ) ∂ m A a a ( y , py )( ∂ py ) m (cid:35) py =0 . (2.34)Odd moments (cid:104) ζ m +1 (cid:105) are zero because of the symmetry (2.30). To compute the moments (cid:104) ζ m (cid:105) , one needs A a a ( y , py ) in the vicinity of py = 0. According to (2.26), this can beevaluated from Euclidean data with nonzero hadron momentum (cid:126)p .Relations analogous to (2.32) to (2.34) can be written down for I tδq δq and B δq δq in theplace of I a a and A a a . We now discuss how the factorisation hypothesis (2.6) for DPDs can be formulated at thelevel of Mellin moments and twist-two functions. At this point, we specialise to the case– 12 –here the hadron h is a π + . This avoids complications due to the proton spin, which arediscussed in [27, section 4.3.1].Let us take the lowest Mellin moment in x and x of (2.6). The Mellin moment of anunpolarised impact parameter dependent parton distribution is (cid:90) − dx f q ( x, b ) = (cid:90) d ∆ (2 π ) e − i b ∆ F q,V ( − ∆ ) , (2.35)where F q,V ( t ) is the form factor of the vector current (cid:104) π + ( p (cid:48) ) | J µq,V (0) | π + ( p ) (cid:105) = ( p + p (cid:48) ) µ F q,V ( t ) with t = ( p − p (cid:48) ) . (2.36)We then obtain from (2.6) I ud ( − y ) ? = (cid:90) d ∆ (2 π ) e − i y ∆ F u,V ( − ∆ ) F d,V ( − ∆ ) . (2.37)We note that thanks to isospin invariance, one has F u,V = − F d,V . As this is not essential inthe present context, we will not use it here.Since one cannot directly determine I ud ( − y ) from Euclidean correlation functions, onecannot directly test (2.37) with lattice data. We therefore derive an analogous relation forthe twist-two function A ud ( y , py ) at py = 0.We recall from [27] that (2.6) can be obtained by inserting a complete set of intermediatestates between the operators O a ( y, z ) and O a (0 , z ) in the DPD definition (2.1) and then assuming that the dominant term in this sum is the ground state. Following exactly the samesteps for the skewed DPD (2.27), one obtains F ud ( x , x , ζ, y ) ? = (cid:90) d ∆ (2 π ) e − i y ∆ − ζ H u (cid:20) x − ζ , ζ − ζ , t ( ∆ , ζ ) (cid:21) × H d (cid:20) x − ζ , ζ − ζ , t ( ∆ , ζ ) (cid:21) (2.38)with t ( ∆ , ζ ) = − ζ m + ∆ − ζ . (2.39)Here H q ( x, ξ, t ) is the generalised parton distribution (GPD) for unpolarised quarks in a pion;its definition can be found e.g. in [85, section 3.2]. The momentum fraction arguments x and ξ of H q are defined in a symmetric way between the incoming and outgoing hadron and partonmomenta, with x referring to the sum of parton momenta and ξ to their difference, and withmomentum fractions normalised to the sum of hadron momenta in the bra and the ket state.Both x and ξ are limited to the interval [ − , x and x . The support of the l.h.s. is shown in figure 2, whereas the one of the r.h.s. is– 13 – d u u − ζx + ζ x − ζ x − ζ x + ζ (a) ζx − ζ x + ζ x + ζ x − ζd du u (b) Figure 3 . (a): Pictorial representation of the r.h.s. of the factorisation hypothesis (2.38). This isobtained by inserting a full set of intermediate states between the operators in the matrix element (cid:104) h | O u O d | h (cid:105) and then retaining only the ground state. (b): The representation obtained when insert-ing the full set of states after reordering the operators to (cid:104) h | O d O u | h (cid:105) . All momentum fractions referto the hadron h in the matrix element. We use form (a) for ζ ≥ ζ < the square delineated by − ζ ≤ x , ≤ − ζ in the ( x , x ) plane. For ζ ≥
0, this missesthe kinematic constraint | x | + | x | ≤ F ud , whereas for ζ < (cid:104) h | O u O d | h (cid:105) = (cid:104) h | O d O u | h (cid:105) . If we insert a set of intermediate states in the lattermatrix element, we obtain (2.38) with ζ replaced by − ζ on the r.h.s. This is representedin figure 3(b). In that case, the mismatch between the support regions of the two sides isless bad for ζ ≤ ζ >
0. We therefore retain (2.38) for ζ ≥ ζ → − ζ on the r.h.s. for ζ <
0. This also satisfies the symmetry in ζ required by PT invariance and stated in (2.29), which is violated if one uses (2.38) for positive and negative ζ .The mismatch of support properties just discussed also affects the case ζ = 0 and is hencenot special to the skewed kinematics we are considering here. In fact, it is well known thatthe factorisation hypothesis (2.6) for DPDs violates the momentum conservation constraint x + x ≤
1. From a theoretical point of view, inserting a full set of intermediate states betweenthe operators in the DPD definition (2.1) or its skewed analogue (2.27) is of course a legitimatemanipulation, but we see that the restriction of this set to the ground state leads to theoreticalinconsistencies such as an incorrect support region or the loss of a symmetry required by PT invariance. How the sum over all states manages to restore the correct properties is difficultto understand in an intuitive manner. We note that a similar observation was made in [84]when discussing the support properties of PDFs and of higher-twist distributions.Integrating both sides of (2.38) over their respective support regions in x and x andusing the sum rule (cid:82) dx H q ( x, ξ, t ) = F q,V ( t ), one obtains I ud ( − y , ζ ) ? = (1 − ζ ) − ζ (cid:90) d ∆ (2 π ) e − i y ∆ F u,V (cid:0) t ( ∆ , ζ ) (cid:1) F d,V (cid:0) t ( ∆ , ζ ) (cid:1) . (2.40)– 14 –sing this for ζ ≥ A ud ( − y , py = 0) ? = 1 π (cid:90) dζ (1 − ζ ) − ζ (cid:90) d ∆ (2 π ) e − i y ∆ F u,V (cid:0) t ( ∆ , ζ ) (cid:1) F d,V (cid:0) t ( ∆ , ζ ) (cid:1) . (2.41)We note that (2.41) is expressed in terms of a two-dimensional vector y . This is differentfrom the factorisation hypothesis we derived in [12, section 5.3], which involved the zero-components of currents and a three-dimensional vector (cid:126)y .Note that the two hypotheses (2.41) and (2.37) are based on the same assumption butare not equivalent to each other. Both are special cases of (2.40), obtained by either setting ζ = 0 or by integrating over ζ from 0 to 1. The assumption that the ground state dominatesthe sum over intermediate states could be a better approximation in one or the other case.Using our lattice results, we will investigate (2.41) in section 4.4 and (2.37) in section 5.4. We performed lattice simulations for the matrix elements (2.8) in a pion with the currentsgiven in (2.9). We set y = 0, so that on the lattice the two currents are inserted at thesame Euclidean time, but with a spatial separation (cid:126)y . We generated data both for zero andnonzero pion momentum (cid:126)p . The lattice techniques we employed are explained in detail in[12, section 3]. In the following, we recall only the basic steps described in that work andthen proceed to the specifics of our present analysis. To evaluate the two-current matrix elements (2.8), we compute the four-point correlationfunction of a pion source operator at Euclidean time 0, a pion sink operator at Euclideantime t , and the two currents J i and J j at Euclidean time τ . The correlation function receivescontributions from a large number of Wick contractions, which are shown in figure 4. We willalso refer to these contractions as “lattice graphs” or simply as “graphs”.The relation between pion matrix elements and lattice graphs depends on the productof C parities of the currents. Omitting Lorentz indices and the dependence on the pionmomentum p , and using the shorthand notation C = C ij ( y ) , C = (cid:2) C ij ( y ) + C ji ( − y ) (cid:3) , A = (cid:2) A ij ( y ) + A ji ( − y ) (cid:3) ,S = (cid:2) S ij ( y ) + S ji ( − y ) (cid:3) , S = S ij ( y ) , D = D ij ( y ) (3.1)for the graphs or their symmetrised combinations, we have M ud,ij ( y ) (cid:12)(cid:12) π + = C + (cid:2) S + D (cid:3) ,M uu,ij ( y ) (cid:12)(cid:12) π + = (cid:2) C + S (cid:3) + (cid:2) S + D (cid:3) ,M ud,ij ( y ) (cid:12)(cid:12) π = (cid:2) S + D (cid:3) − A ,M uu,ij ( y ) (cid:12)(cid:12) π = C + (cid:2) S + D (cid:3) + (cid:2) C + S (cid:3) + A (3.2)– 15 – i ( y ) J i ( y ) J j (0) C ij ( y ) = = η iC η jC × C ij ( y ) = A ij ( y ) = J i ( y ) J j (0) J j (0) J j (0) J i ( y ) S ij ( y ) = = η iC η jC × S ij ( y ) = J i ( y ) D ij ( y ) = J i ( y ) J j (0) J i ( y ) J j (0) J j (0) J i ( y ) J j (0) Figure 4 . Lattice graphs for the correlation functions used to extract the matrix elements (2.8) ina pion. The dependence on the pion momentum (cid:126)p is not indicated for brevity. η i C denotes the chargeconjugation parity of the current J i and is defined in (2.12). – 16 –or η i C η j C = +1, which is satisfied for all combinations of currents considered in the presentstudy. We note that this is no longer the case if one includes operators with covariantderivatives (corresponding to higher Mellin moments). One readily checks that (3.2) satisfiesthe general symmetry relation (2.17), as it must.To compute the different graphs on the lattice, we use a variety of techniques as detailedin [12, section 3.3]. We make extensive use of stochastic sources, and for graph C we use ahopping parameter expansion to reduce statistical noise for the propagation between the twocurrents.For the disconnected graphs S and D we need to subtract vacuum contributions, namelythe product of a two-point correlation function of the pion source and sink with a two-pointcorrelation function of the two currents. The latter corresponds to the vacuum expecta-tion value (cid:104) | J i ( y ) J j (0) | (cid:105) . The vacuum subtraction for the disconnected graph S involves (cid:104) | J i ( y ) | (cid:105) or (cid:104) | J j (0) | (cid:105) , which is zero because our currents carry Lorentz indices.We anticipate that the doubly disconnected graph D in general gives a good signal forthe four-point correlation function, but that there is a near-perfect cancellation between thiscorrelator and its vacuum subtraction term. The result after subtraction is consistent withzero and has huge statistical uncertainties compared with those of any other graph. We willhence not be able to report useful results for graph D . Fortunately, we encounter no suchproblem for graph S . We perform our simulations using the Wilson gauge action and n F = 2 mass degenerateflavours of non-perturbatively improved Sheikholeslami-Wohlert (NPI Wilson-clover) fermions.The gauge configurations were generated by the RQCD and QCDSF collaborations. We usetwo gauge ensembles, whose parameters are given in table 1. They have different spatial sizes, L = 32 and L = 40, which allows us to study finite volume effects in section 3.5. Despitehaving data for only a single lattice spacing, a = 0 .
071 fm, we are also able to investigatediscretisation effects, as discussed in section 3.3.For the ensemble with L = 40, we performed simulations with different κ values in theensemble β a [fm] κ L × T m π [MeV] Lm π N full N used IV 5.29 0.071 0.13632 32 ×
64 294 . ± . ×
64 288 . ± . Table 1 . Details of the gauge ensembles used in this analysis. N full is the total number of availablegauge configurations, and N used is the number of configurations used in our simulations. More detailcan be found in [86, 87]. – 17 –uark propagator:light quarks: κ = 0 . , m π = 293 MeV , strange: κ = 0 . , m π = 691 MeV , charm: κ = 0 . , m π = 3018 MeV . (3.3)Here “light quarks” refers to the κ value used for simulating the sea quarks, whereas theother two values correspond to the physical strange and charm quark masses, as determinedin [88] and [89] by tuning the pseudoscalar ground state mass to 685 . S -wave charmonium mass to 3068 . n F = 2 fermion action, the strange and charm quarks arepartially quenched.The values of m π in (3.3) are obtained from exponential fits of the pion two-point function.We quote them only for orientation and do not attempt to quantify their errors. These massesare in reasonable agreement with the value in table 1 for light quarks, and with the mass ofthe pseudoscalar ground state quoted below (3.3) for strange quarks. Pion matrix elements.
For all lattice graphs, we compute the correlation functions withzero three-momentum (cid:126)p of the pion. For the connected graphs C and C , we additionallyhave data with finite pion momenta. These data are restricted to the L = 40 lattice and tolight quarks. The pion momenta that can be realised on the lattice are given by (cid:126)p = 2 πLa (cid:126)P , (3.4)where the components of (cid:126)P are integers and 2 π/ ( La ) ≈
437 MeV in our case. For simplicitywe write P = | (cid:126)P | . Graph C is computed for all 18 nonzero momenta with P ≤
2, for 6momenta with P = 3, and for one momentum with P = 4. For C , we have results for all6 momenta with P = 1.The distance between the pion source and sink in the correlation functions is fixed to t = 15 a ≈ .
07 fm as a default. To investigate the influence of excited states, we also calculategraphs C , C and S with t = 32 a . The matrix element (2.8) is extracted from the ratiobetween the four-point correlation function around τ = t/ C and A , we measure the τ dependence of the four-point function and fit to aplateau in the τ ranges specified in [12, equation (4.1)]. The quality of the correspondingplateaus is good for matrix elements that have a nonzero value within statistical uncertainties.For the remaining graphs, we extract the matrix element from data at τ /a = 7 and 8 if t/a = 15. For the C and S data with t/a = 32, we use τ /a = 16. A comparison of datawith t = 15 a and t = 32 a is shown in section 3.4.All lattice currents are converted to the MS scheme at the renormalisation scale µ = 2 GeV . (3.5)– 18 –s described in [12, section 3.4], this is done using a combination of non-perturbative andperturbative renormalisation and includes an estimate of the quark mass dependent order a improvement term. Invariant functions.
From the matrix elements (2.8), we determine the invariant functionsfor each individual value of (cid:126)y and (cid:126)p . This is done using a minimum χ fit of the data for alltensor components to the decomposition (2.20). For invariant functions of twist two, we alsouse the projector method (2.23). In both cases, the statistical error of an invariant functionat given (cid:126)y and (cid:126)p is computed using the jackknife method. To eliminate autocorrelations, wetake the number of jackknife samples as 1 / N used of gauge configurationsgiven in table 1.For P = 0, the twist-two functions extracted with one or the other method show excellentagreement with each other and have statistical uncertainties of almost the same size. For P >
0, the values obtained with the projection method have much larger statistical errorsthan those obtained with a fit and provide only a very weak cross check. All data shown inthe following are obtained by the fit method, both for P = 0 and P > py = 0 here, because they have much smaller statistical errors than thedata for py (cid:54) = 0. We will return to the case of nonzero py in section 5.When discussing twist-two functions extracted from the correlation functions for par-ticular lattice graphs, we will generically write A qq , A ∆ q ∆ q , . . . , B δqδq , without reference tospecific quark flavors q and q . This is because the distinction between u and d quarks in apion only appears when lattice graphs are combined as specified in (3.2). The decomposition (2.20) of matrix elements in terms of basis tensors and functions of y and py assumes Lorentz invariance and thus requires both the continuum and the infinite-volume limit. If our lattice simulations are sufficiently close to these limits, then the valuesof twist-two functions extracted for individual points (cid:126)y and (cid:126)p with (cid:126)p(cid:126)y = 0 must not dependon the directions of (cid:126)y or (cid:126)p or on the size of (cid:126)p .Let us test whether this is the case in our simulations for light quarks on the lattice with L = 40. We restrict our attention to graphs C and C , for which statistical errors are smallenough to reveal the effects of interest. For the sake of legibility, we henceforth write y = | (cid:126)y | for the length of the spatial distance (cid:126)y between the two currents. We continue to use y and py to denote the products y µ y µ and p µ y µ of four-vectors in Minkowski space. Since y µ isalways spacelike in our context, this implies that y < y of order La/
2, we see a clear anisotropy with a saw-tooth pattern in all twist-two functions that have sufficiently small errors. Examples are shown in figure 5. This patternis expected on a lattice with periodic boundary conditions and can be understood in termsof “mirror images”. The same effect has been seen and discussed in previous lattice studies– 19 – 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] cos 𝜃 < 0.9 cos 𝜃 ≥ 0.9 (a) A qq , graph C −0.01−0.008−0.006−0.004−0.00200.002 1.1 1.2 1.3 1.4 1.5 1.615 16 17 18 19 20 21 22 23 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 𝛿𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] cos 𝜃 < 0.9 cos 𝜃 ≥ 0.9 (b) A δqq , graph C Figure 5 . Examples for the anisotropy in graph C at large y = | (cid:126)y | . The data shown are for L = 40,zero pion momentum, and light quarks. The results in this and all the following figures are given inthe MS scheme at the scale µ = 2 GeV. The error bars shown in the plots are statistical and obtainedwith the jackknife method. of two-current correlators [8, 11], including our study in [12] that employed the same latticedata as the present work. As shown in [8], the effect of mirror charges at a given distance y issmallest for points (cid:126)y close to one of the space diagonals, i.e. the lines given by (cid:126)z = ( z , z , z )with | z | = | z | = | z | . To quantify this, we define θ ( (cid:126)y ) as the angle between (cid:126)y and the spacediagonal in the same octant as (cid:126)y . In [12, section 4.2], we found that a cutcos θ ( (cid:126)y ) ≥ . y , whilst keeping sufficientstatistics.A different type of anisotropy in the C data is observed at small y and shown in figure 6.For A ∆ q ∆ q , A δqq and B δqδq , the data with zero pion momentum exhibit a clear discrepancybetween points (cid:126)y on a coordinate axis (i.e. with two components being zero) and all otherpoints. This discrepancy is very strong for y below 5 a and ceases to be visible above 7 a . Thedata for A δqδq (not shown in the figure) have larger errors and show only a weak anisotropyfor y < a . Only the function A qq is not affected by this phenomenon, for which we have noexplanation.By contrast, we find that the C data with nonzero pion momentum and py = 0 areisotropic in (cid:126)y . For nonzero momenta, we can hence average all data with the same valuesof y and P , which greatly decreases statistical errors. We find good agreement between the P > P = 0 data with (cid:126)y on a coordinate axis for all twist-two functions except A ∆ q ∆ q , where the agreement is only approximate.– 20 – 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 𝑞𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 ) 𝑦[𝑎]𝑃 = 0 , cos 𝜃 > √1/3𝑃 = 1𝑃 = 0 , cos 𝜃 = √1/3 (a) A qq , graph C −0.04−0.03−0.02−0.0100.010.020.030.040.050.06 0.2 0.3 0.4 0.5 0.6 0.72 3 4 5 6 7 8 9 10 𝐴 Δ 𝑞 Δ 𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 Δ𝑞Δ𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 ) 𝑦[𝑎]𝑃 = 0 , cos 𝜃 > √1/3𝑃 = 1𝑃 = 0 , cos 𝜃 = √1/3 (b) A ∆ q ∆ q , graph C −0.25−0.2−0.15−0.1−0.050 0.2 0.3 0.4 0.5 0.6 0.72 3 4 5 6 7 8 9 10 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 𝛿𝑞𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 ) 𝑦[𝑎]𝑃 = 0 , cos 𝜃 > √1/3𝑃 = 1𝑃 = 0 , cos 𝜃 = √1/3 (c) A δqq , graph C −0.25−0.2−0.15−0.1−0.0500.050.10.150.2 0.2 0.3 0.4 0.5 0.6 0.72 3 4 5 6 7 8 9 10 𝐵 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐵 𝛿𝑞𝛿𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 ) 𝑦[𝑎]𝑃 = 0 , cos 𝜃 > √1/3𝑃 = 1𝑃 = 0 , cos 𝜃 = √1/3 (d) B δqδq , graph C Figure 6 . Twist-two functions at small y for graph C , with scaled pion momenta P = 0 and P = 1as defined below (3.4). All points have py = 0 and are for L = 40 and light quarks. The data for P = √ P = √ P = 2 agree with those for P = 1 within errors but are not shown for the sakeof clarity. Data with cos θ = (cid:112) / (cid:126)y on a coordinate axis. We now turn our attention to graph C at small y . Here, we find a very strong anisotropyin the P = 0 data. This is shown in figure 7, where we distinguish points (cid:126)y on the coordinateaxes, which have cos θ = (cid:112) /
3, points with (cid:112) / < cos θ ≤ (cid:112) /
3, and points with (cid:112) / < cos θ . We note that points in a coordinate plane, i.e. with at least one component of (cid:126)y equalto zero, have (cid:112) / ≤ cos θ ≤ (cid:112) /
3. In all channels, we see a clear discrepancy between thepoints (cid:126)y on a coordinate axis and all other points. In addition, there is a significant mismatchbetween points with cos θ above or below (cid:112) / A δqq .We recall that a strong anisotropy for C at small y was also seen for the correlation functionsin our study [12]. In section 4.2 of that work, we argued that this reflects an anisotropy in the– 21 – −2 −1 | 𝐴 𝑞𝑞 | [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 𝑞𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎]√2/3 ≥ cos 𝜃 > √1/3 cos 𝜃 > √2/3 cos 𝜃 = √1/3 (a) | A qq | , graph C −4 −3 −2 −1 | 𝐴 Δ 𝑞 Δ 𝑞 | [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 Δ𝑞Δ𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎]√2/3 ≥ cos 𝜃 > √1/3 cos 𝜃 > √2/3 cos 𝜃 = √1/3 (b) | A ∆ q ∆ q | , graph C −3 −2 −1 | 𝐴 𝛿 𝑞𝑞 | [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 𝛿𝑞𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎]√2/3 ≥ cos 𝜃 > √1/3 cos 𝜃 > √2/3 cos 𝜃 = √1/3 (c) | A δqq | , graph C −4 −3 −2 −1 | 𝐵 𝛿 𝑞 𝛿 𝑞 | [ f m − ] 𝑦[ fm ]𝐶 , 𝐵 𝛿𝑞𝛿𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎]√2/3 ≥ cos 𝜃 > √1/3 cos 𝜃 > √2/3 cos 𝜃 = √1/3 (d) | B δqδq | , graph C Figure 7 . Twist-two functions at small y for graph C . All points are for L = 40, zero pionmomentum, and light quarks. For A δqδq (not shown), one finds a clear anisotropy at y < a , whilstat larger y the statistical errors are too large for drawing firm conclusions. lattice propagator between the two currents, and that points selected by the cut (3.6) shouldgive the most reliable results according to the analysis in [90].We also have P = 1 data for C , which we can compare with those for P = 0. As seen infigure 8, for A qq and A δqq the data at P = 1 are inconsistent with those at P = 0, regardlessof the value of cos θ in the latter. Since for P = 1 the condition py = 0 requires (cid:126)y to liein a coordinate plane, we can in fact not select points satisfying the cut (3.6) in this case.We therefore discard our data with nonzero P for C . Testing boost invariance of twist-twofunctions at py = 0 in the presence of the cut (3.6) would require data with at least P = √ C .At P = 0 and small y , we are now in a difficult situation. Points with large cos θ are– 22 – 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 𝑞𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 ) 𝑦[𝑎]𝑃 = 0 , cos 𝜃 ≤ √2/3𝑃 = 0 , cos 𝜃 > √2/3𝑃 = 1 (a) A qq , graph C −0.06−0.04−0.0200.020.040.060.080.10.12 0.4 0.5 0.6 0.7 0.8 0.9 16 8 10 12 14 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 𝛿𝑞𝑞 , 𝑝𝑦 = 0 ( 𝐿 = 40 ) 𝑦[𝑎]𝑃 = 0 , cos 𝜃 ≤ √2/3𝑃 = 0 , cos 𝜃 > √2/3𝑃 = 1 (b) A δqq , graph C Figure 8 . Twist-two functions at intermediate y for graph C . All points are for py = 0, L = 40and light quarks. preferred for C , while for C the points with the smallest possible value cos θ = (cid:112) / P > θ to the data for C and C would prevent us from taking linear combinations of those graphsat a given (cid:126)y . However, we regard combining data point by point in (cid:126)y as highly desirable fora transparent and consistent treatment of statistical correlations in the jackknife analysis.To avoid this problem, we choose to discard points with y < a in our further analysis,and to apply the cut in (3.6) to the P = 0 data for all lattice graphs. After this cut, datapoints with equal values of y are averaged also for P = 0. We thus avoid the regions where theanisotropy for C and C seen at P = 0 is most severe. For C , a small discrepancy betweenthe data with P = 0 and P > y ∼ a , but we consider this tobe at an acceptable level. The result of this procedure is shown for graph C in figure 9.The agreement between the data for different pion momenta is quite good, except for thefunction A ∆ q ∆ q .As an exception to the selection just described, we will in section 4.4 use the C datafor A qq down to y = 3 a , given that in this particular channel there are no indications ofanisotropy or a pion momentum dependence, as can be seen in figure 6(a). As specified in section 3.2, we have a limited set of data with a separation of t = 32 a betweenpion source and sink. Comparing this with our results for t = 15 a allows us to assess therelevance of excited state contributions in our extraction of the pion matrix elements (2.8).On our lattice with size L = 32, we have t = 32 a data for graphs S and C . Unfortu-nately, these graphs give a statistical signal consistent with zero for all twist-two functions– 23 – 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 , momentum comparison ( 𝐿 = 40 ) 𝑦[𝑎] 𝑃 = √3𝑃 = √2𝑃 = 1𝑃 = 0 , cos 𝜃 ≥ 0.9 (a) A qq , graph C −0.015−0.01−0.00500.0050.010.015 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 Δ 𝑞 Δ 𝑞 [ f m − ] 𝑦[ fm ]𝐴 Δ𝑞Δ𝑞 , momentum comparison (
𝐿 = 40 ) 𝑦[𝑎] 𝑃 = √3𝑃 = √2𝑃 = 1𝑃 = 0 , cos 𝜃 ≥ 0.9 (b) A ∆ q ∆ q , graph C −0.12−0.1−0.08−0.06−0.04−0.0200.02 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝛿𝑞𝑞 , momentum comparison ( 𝐿 = 40 ) 𝑦[𝑎] 𝑃 = √3𝑃 = √2𝑃 = 1𝑃 = 0 , cos 𝜃 ≥ 0.9 (c) A δqq , graph C −0.00200.0020.0040.0060.0080.010.0120.0140.016 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝛿𝑞𝛿𝑞 , momentum comparison ( 𝐿 = 40 ) 𝑦[𝑎] 𝑃 = √3𝑃 = √2𝑃 = 1𝑃 = 0 , cos 𝜃 ≥ 0.9 (d) A δqδq , graph C Figure 9 . Comparison of twist-two functions for graph C at different values of the scaled pionmomentum P . All points are for py = 0, L = 40 and light quarks. For B δqδq (not shown), one findsgood agreement between all points, with a similar quality as for A qq . and for both source-sink separations. We hence limit the following discussion to graph C onour L = 40 lattice.In general, we find that the results for the two source-sink separation agree reasonablywell for light quarks, as illustrated in the upper panels of figure 10. For strange quarks, thedata have smaller statistical errors and we can clearly see discrepancies between t = 15 a and 32 a , as shown in the lower panels of the figure. Except for the case of A ∆ q ∆ q , thesediscrepancies are, however, small when compared with the size of the twist-two functions.In our data for charm quarks, the statistical signal and the agreement between the twosource-sink separations is excellent for all twist-two functions, and even better than the onein figure 10(a). With the exception mentioned above, we thus find no indication for a sizeable– 24 –ontamination from excited states in our results. −0.12−0.1−0.08−0.06−0.04−0.020 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 𝑞𝑞 , 𝑝 = 0 , 𝐿 = 40𝑦[𝑎] 𝑡/𝑎 = 32𝑡/𝑎 = 15 (a) A qq , graph C , light quarks −0.05−0.04−0.03−0.02−0.0100.01 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐵 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐵 𝛿𝑞𝛿𝑞 , 𝑝 = 0 , 𝐿 = 40𝑦[𝑎] 𝑡/𝑎 = 32𝑡/𝑎 = 15 (b) B δqδq , graph C , light quarks −0.02−0.015−0.01−0.00500.005 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 Δ 𝑞 Δ 𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐴 Δ𝑞Δ𝑞 , 𝑝 = 0 , 𝐿 = 40 , strange 𝑦[𝑎] 𝑡/𝑎 = 32𝑡/𝑎 = 15 (c) A ∆ q ∆ q , graph C , strange quarks −0.02−0.010 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐵 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝐶 , 𝐵 𝛿𝑞𝛿𝑞 , 𝑝 = 0 , 𝐿 = 40 , strange 𝑦[𝑎] 𝑡/𝑎 = 32𝑡/𝑎 = 15 (d) B δqδq , graph C , strange quarks Figure 10 . Comparison of twist-two functions extracted for graph C with source-sink separations t/a = 15 or 32. All data are for L = 40, zero pion momentum, and subject to the cut (3.6). Let us finally compare our simulations for light quarks on the lattices with L = 40 and L = 32.In general, the data for the smaller lattice have larger jackknife errors. This is to be expectedfrom the parameters that determine the statistical averaging in our simulations. Details forthese are given in table 2 of [12].For twist-two functions with a small relative error, we typically find a weak volumedependence compared with the size of the functions themselves, as shown in panels (a) to (c)of figure 11. In the case of panel (b), this dependence is, however, statistically significant.For functions that have large relative errors, the volume dependence appears to be more– 25 – 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ] Volume comparison for 𝐶 , 𝐴 𝑞𝑞 , 𝑝 = 0𝑦[𝑎] 𝐿 = 32𝐿 = 40 (a) A qq , graph C −0.1−0.08−0.06−0.04−0.0200.02 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ] Volume comparison for 𝐶 , 𝐴 𝛿𝑞𝑞 , 𝑝 = 0𝑦[𝑎] 𝐿 = 32𝐿 = 40 (b) A δqq , graph C −0.05−0.04−0.03−0.02−0.0100.01 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐵 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ] Volume comparison for 𝐶 , 𝐵 𝛿𝑞𝛿𝑞 , 𝑝 = 0𝑦[𝑎] 𝐿 = 32𝐿 = 40 (c) B δqδq , graph C −0.0500.050.10.150.2 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ] Volume comparison for 𝐶 , 𝐴 𝛿𝑞𝑞 , 𝑝 = 0𝑦[𝑎] 𝐿 = 32𝐿 = 40 (d) A δqq , graph C −0.01−0.00500.0050.010.0150.02 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐵 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ] Volume comparison for A, 𝐵 𝛿𝑞𝛿𝑞 , 𝑝 = 0𝑦[𝑎] 𝐿 = 32𝐿 = 40 (e) B δqδq , graph A −0.4−0.3−0.2−0.100.10.2 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ] Volume comparison for 𝑆 , 𝐴 𝛿𝑞𝑞 , 𝑝 = 0𝑦[𝑎] 𝐿 = 32𝐿 = 40 (f) A δqq , graph S Figure 11 . Comparison of data for the two different lattice sizes in our study. All points are forzero pion momentum, light quarks, and subject to the cut (3.6). – 26 –ronounced in some cases, especially at low y . An example is figure 11(e). One may take thisas a general warning against over-interpreting statistically weak signals in our simulations. In this section, we present our results for the twist-two functions at py = 0. All data shownin the following are for zero pion momentum and have been extracted from the lattice with L = 40 with our standard source-sink separation t = 15 a . The data selection described atthe end of section 3.3 removes regions in which we see strong lattice artefacts in the form ofbroken rotational or boost symmetry.As we explained in section 2.3, twist-two functions at py = 0 are not directly related tothe Mellin moments of DPDs. Instead, they are Mellin moments of skewed DPDs, integratedover the skewness parameter ζ . As seen in figure 2, these moments receive contributions fromparton configurations that are different from those in a DPD at ζ = 0. When interpretingthe results of the present section, we will assume that these configurations are not dominant,and that the qualitative features of invariant functions at py = 0 are the same as for Mellinmoments of DPDs at ζ = 0. The results presented in section 5.3 will lend support to thisassumption.Notice that each of the lattice graphs in figure 4 can contribute to each of the partonicregimes shown in figure 2. Examples for different regimes of the connected graphs C and C are shown in figures 12 and 13. In figures 14 and 15, we compare the contributions from different lattice graphs to the twist-two functions for light quarks. The contributions from graphs S and C are multiplied witha factor 2 in the figures, since they always appear with this weight in physical matrix elementsaccording to (3.2).For all twist-two functions except A ∆ q ∆ q , graph C gives a very clear signal, which ispositive for A δqδq and negative for the other functions. By comparison, the signal for theannihilation graph A is smaller than the one for C by an order of magnitude or more, exceptfor y > a , where the statistical errors prevent us from making a clear statement. The u ¯ d u ¯ d | ¯ du u ¯ d u ¯ u ¯ d | ¯ du ¯ d ¯ u ¯ d | ¯ d ¯ u Figure 12 . Examples for the partonic regimes of graph C in a π + . The notation u ¯ u ¯ d | ¯ d is the sameas in figure 2, i.e. the partons to the left of the vertical bar belong to the pion on the left of the graph. – 27 – ¯ d ¯ d ¯ d | ¯ d ¯ du ¯ d u ¯ u | ¯ uu u ¯ d uu | uu u ¯ d u ¯ uu | uu ¯ d ¯ u ¯ u | ¯ u ¯ u Figure 13 . As figure 12, but for graph C . function A ∆ q ∆ q shows a different behaviour, with C and A being of similar size and muchsmaller than C for all other twist-two functions. We recall from section 3.3 that A ∆ q ∆ q ismore strongly affected by lattice artefacts than the other channels, see figure 9(b).A clear signal for the connected graph C is only seen for A qq and A δqq , with a signopposite to the one for graph C . This signal is most important at small y . For the graphs S and S with one disconnected fermion loop, the signal we obtain is rather noisy in all channels.For graph D with two disconnected fermion loops, the signal after vacuum subtraction is evenmore noisy and not shown.From our simulations with the strange quark mass, we have only data for graphs C and A .In all channels, we obtain an excellent signal for C , whereas for A the statistical significanceis typically not much larger than one standard deviation. In the region 5 a ≤ y ≤ a , we findthat A is smaller than C by one to two orders of magnitude, except for A ∆ q ∆ q . For A qq and A δqδq , we see in figure 16 that the behaviour of A is quite flat, unlike the one of C , so thatat large y the two graphs become more comparable in size. As in the case of light quarks, thefunction A ∆ q ∆ q behaves differently, with graph A being smaller than C at y ∼ a and thedata for both graphs having a zero crossing a bit below y = 9 a . Recall, however, that alsofor strange quarks we see stronger lattice artefacts in A ∆ q ∆ q than in other channels, as seenin figure 10(c).From our simulations with the charm quark mass, we have data for all graphs except S .A clear nonzero signal is seen for C and C up to y ∼ a to 15 a , with 2 C being smallerthan C by at least one order of magnitude. The signal for A and S is in general consistentwith zero. The only exception to this is A δqq . For this function, we see a clear signal for 2 S at y around 5 a , which is about 50 times smaller than the one for C . We also see a weak 1 σ signal for A , which we do not wish to over-interpret.By and large, we find that for all quark masses the only graphs that give signals of– 28 – 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝑆 𝐴𝐶 −0.1−0.0500.050.10.150.2 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐶 𝑆 (a) A qq , light quarks −0.06−0.04−0.0200.020.040.060.08 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 Δ 𝑞 Δ 𝑞 [ f m − ] 𝑦[ fm ]𝐴 Δ𝑞Δ𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝑆 𝐴𝐶 −0.14−0.12−0.1−0.08−0.06−0.04−0.0200.020.040.06 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 Δ 𝑞 Δ 𝑞 [ f m − ] 𝑦[ fm ]𝐴 Δ𝑞Δ𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐶 𝑆 (b) A ∆ q ∆ q , light quarks −0.1−0.08−0.06−0.04−0.0200.020.04 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝛿𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝑆 𝐴𝐶 −0.15−0.1−0.0500.050.10.150.20.25 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝛿𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐶 𝑆 (c) A δqq , light quarks Figure 14 . Contributions of the different lattice graphs to twist-two functions for light quarks. – 29 – 𝐴 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝛿𝑞𝛿𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝑆 𝐴𝐶 −0.06−0.05−0.04−0.03−0.02−0.0100.010.020.03 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝛿𝑞𝛿𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐶 𝑆 (a) A δqδq , light quarks −0.08−0.06−0.04−0.0200.020.040.060.080.1 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐵 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝐵 𝛿𝑞𝛿𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝑆 𝐴𝐶 −0.2−0.15−0.1−0.0500.050.1 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐵 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝐵 𝛿𝑞𝛿𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐶 𝑆 (b) B δqδq , light quarks Figure 15 . Continuation of figure 14. appreciable size are C and, in several cases, C . We therefore take a closer look at thesegraphs in the next subsection. The annihilation graph is negligible, except in the case of A ∆ q ∆ q for light or strange quarks, where the signal from graph C is small by itself. Disconnectedgraphs either have a negligibly small signal or large statistical errors. The contribution of graph C to the twist-two function A qq for unpolarised partons is negativefor all three quark masses in our study. We recall from (2.19) that the regime with a quarkand an antiquark in the pion contributes with a negative sign to the lowest Mellin momentof a DPD. The same holds for the Mellin moment of a skewed DPD, and hence for A qq at py = 0. A negative sign of A qq is easily understood by the dominance of the valence q ¯ q Fockstate, which is probed by graph C as shown in the first panel of figure 12.– 30 – −6 −5 −4 −3 −2 −1 | 𝐴 𝑞𝑞 | [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎] 𝐴𝐶 (a) | A qq | , graphs C and A −5 −4 −3 −2 −1 | 𝐴 𝛿 𝑞 𝛿 𝑞 | [ f m − ] 𝑦[ fm ]𝐴 𝛿𝑞𝛿𝑞 , 𝑝 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎] 𝐴𝐶 (b) | A δqδq | , graphs C and A Figure 16 . Comparison of graphs A and C for strange quarks. Here and in subsequent figures witha logarithmic scale, we stop showing data for individual quark masses at a value of y beyond whichthe error bars become so large that they would obscure the plot. The situation is different for graph C , whose partonic representation always involves ahigher Fock state of the pion. The Z -graphs in figure 13 probe the qq , ¯ q ¯ q and q ¯ q regimes in asimilar manner. We find that for all quark masses, the contribution of C to A qq is positive,which means that for a given distance y this graph gives a larger probability for finding a qq or ¯ q ¯ q rather than a q ¯ q pair.Let us now take a closer look at the mass dependence of our results for graph C . Wemultiply A δqq and B δqδq with the power of the meson mass m with which they appear in thedecomposition (2.20) of two-current matrix elements. We see in figure 17 that for all twist-twofunctions except A ∆ q ∆ q , the decrease with y becomes stronger with increasing quark mass,which simply reflects the decreasing size of the meson. At y ∼ a , the functions A qq , mA δqq and m B δqδq are of comparable size for all quark masses, whereas A δqδq increases with themass. The behaviour of A ∆ q ∆ q for light and strange quarks is qualitatively different fromthe one of the other functions, as is evident from figure 18. For charm quarks, A ∆ q ∆ q isapproximately exponential in y , with a logarithmic slope similar to the one of A qq . A fit ofthe y dependence of the twist-two functions for light quarks is presented in section 5.2.We now discuss graph C , for which we have data with light quarks and with charm.For the functions A ∆ q ∆ q , A δqδq and B δqδq , the light quark data is too noisy for a meaningfulcomparison with charm results, so that we focus on A qq and A δqq . As is seen in figure 19,the size of both functions is significantly smaller for charm quarks. This is plausible: asdiscussed in the previous subsection, the partonic interpretation of graph C always involvesa Fock state with at least two quarks and two antiquarks in the meson, whereas for C wehave the regime shown in the first panel in figure 12, which involves only the quark-antiquark– 31 – −5 −4 −3 −2 −1 − 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎] charmstrangelight (a) − A qq , graph C −5 −4 −3 −2 −1 − 𝑚 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝑚𝐴 𝛿𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎] charmstrangelight (b) − mA δqq , graph C −5 −4 −3 −2 −1 𝐴 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝛿𝑞𝛿𝑞 , 𝑝 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎] charmstrangelight (c) A δqδq , graph C −5 −4 −3 −2 −1 − 𝑚 𝐵 𝛿 𝑞 𝛿 𝑞 [ f m − ] 𝑦[ fm ]𝑚 𝐵 𝛿𝑞𝛿𝑞 , 𝑝 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎] charmstrangelight (d) − m B δqδq , graph C Figure 17 . Mass dependence of twist-two functions for graph C . −0.05−0.04−0.03−0.02−0.0100.01 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 Δ 𝑞 Δ 𝑞 [ f m − ] 𝑦[ fm ]𝐴 Δ𝑞Δ𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] charmstrangelight (a) A ∆ q ∆ q , graph C Figure 18 . As figure 17, for A ∆ q ∆ q and with a linear instead of a logarithmic scale. – 32 – −5 −4 −3 −2 −1 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎] charmlight (a) A qq , graph C −5 −4 −3 −2 −1 𝑚 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝑚𝐴 𝛿𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 , log) 𝑦[𝑎] charmlight (b) mA δqq , graph C Figure 19 . Mass dependence of twist-two functions for graph C . Fock state. The y dependence of A qq and A δqq is also qualitatively different for the twomasses: for charm we observe a clear and steep exponential falloff, whereas for light quarks,the logarithmic slope of both functions decreases around y ∼ . A major aim of our study is to investigate the strength and pattern of spin correlationsbetween two partons in a pion. We spelled out the physical interpretation of polarised DPDsin section 2.1. This interpretation extends to the corresponding twist-two functions at py = 0,provided that these are dominated by partonic regimes associated with DPDs at ζ = 0. Underthis assumption, comparing A ∆ q ∆ q and A δqδq with A qq indicates whether two partons preferto have their spins aligned or anti-aligned, with A ∆ q ∆ q referring to longitudinal and A δqδq to transverse polarisation. We will refer to these as “spin-spin correlations”. Note that,according to (2.19), a q ¯ q pair with aligned spins contributes with a negative sign to A qq and A δqδq and with a positive sign to A ∆ q ∆ q , whereas a qq pair with aligned spins contributes witha positive sign to all three functions. The comparison of myA δqq and m | y | B δqδq with A qq tells us about the strength of correlations between the transverse spin of one or both observedpartons and the distance y between these partons in the transverse plane. We refer to thisas “spin-orbit correlations” in the following. The pre-factors my and m | y | in myA δqq and m | y | B δqδq follow from the decompositions (2.4) and (2.18).We note that the probability interpretation of polarised DPDs implies positivity con-straints [91] that extend the well-known Soffer bound for single parton distributions [92].These bounds imply that | f ∆ q ∆¯ q | , | f δqδ ¯ q | , | myf δq ¯ q | and | m y f tδqδ ¯ q | are bounded by f q ¯ q . Cor-responding bounds do not hold for the lowest Mellin moments of DPDs because of the relativeminus sign between quark and antiquark contributions in (2.19). They hold even less for the– 33 –oments of skewed DPDs, which do not represent probabilities to start with. Nevertheless, ina loose sense, the size of A qq sets a natural scale for the other twist-two functions (multipliedwith my or m | y | as appropriate).Starting our discussion with graph C , we see in the top panels of figure 20 that by farthe strongest polarisation effect seen for light quarks is the spin-orbit correlation for a singleparton, followed by the spin-orbit correlation involving both partons. Both the transverseand the longitudinal spin-spin correlations are very small. This is completely different fromthe simple picture of a pion as a q ¯ q pair in an S -wave, for which one would obtain 100%anti-alignment of both transverse and longitudinal spins.All spin correlations increase considerably with the quark mass. For charm quarks, myA δqq is almost as large as A qq . Spin-spin correlations are also important for charm: thespins of the quark and antiquark are anti-aligned by about 75% for transverse and by about50% for longitudinal polarisation. We note that this is still quite far away from the non-relativistic limit, in which transverse and longitudinal spin correlations become equal.We note that the pion mass for our simulations with light quarks, m π ≈
295 MeV, isquite a bit larger than the physical value. A naive extrapolation of the polarisation patternsjust described suggests that at the physical point the spin-orbit correlation for one polarisedparton may be substantial, whilst correlations involving two quark spins might be even smallerthan the ones we see for light quarks in the present study.We now turn to our results for graph C , which are shown in figure 21. For light quarks,we see a substantial spin-orbit correlation of order 50% for a single parton. The spin-spincorrelation for longitudinal polarisation is also of order 50% for y ∼ .
35 fm, but quickly de-creases and is negligible already around y ∼ . C are appreciable, apart from the one described by m | y | B δqδq . Notice that A ∆ q ∆ q has the same sign for C and C , unlike all other twist-twofunctions. If (as suggested by the sign of A qq ) the dominant parton configuration probedby the twist-two operators is a c ¯ c pair for graph C and a cc pair for graph C , then thelongitudinal parton spins tend to be anti-aligned in both cases. We now test the factorisation hypothesis for A ud ( y , py = 0) that we derived in section 2.4.We restrict ourselves to the contribution from the connected graph C . Taking the fullcombination of graphs in the first line of (3.2) is not an option because of the huge errors inour results for the doubly disconnected graph D . By contrast, we see in figure 14(a) that S is consistent with zero for A qq (although within errors much larger than those on C ). Wefind it plausible to expect that the contribution from D is even smaller than the one of S ,since D has two disconnected fermion loops with one operator insertion.The factorisation hypothesis (2.41) involves the vector form factor of the pion. We haveextracted this form factor from our lattice simulations, using the full number of 2025 gauge– 34 – 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 𝛿𝑞𝛿𝑞 vs 𝐴 𝛿𝑞𝑞 vs 𝐵 𝛿𝑞𝛿𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝐴 𝛿𝑞𝛿𝑞 𝑚𝑦𝐴 𝛿𝑞𝑞 𝑚 |𝑦 |𝐵 𝛿𝑞𝛿𝑞 −0.12−0.1−0.08−0.06−0.04−0.0200.020.04 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 Δ𝑞Δ𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝐴 Δ𝑞Δ𝑞 (a) graph C , light quarks −0.12−0.1−0.08−0.06−0.04−0.0200.020.04 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 𝛿𝑞𝛿𝑞 vs 𝐴 𝛿𝑞𝑞 vs 𝐵 𝛿𝑞𝛿𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝐴 𝛿𝑞𝛿𝑞 𝑚𝑦𝐴 𝛿𝑞𝑞 𝑚 |𝑦 |𝐵 𝛿𝑞𝛿𝑞 −0.12−0.1−0.08−0.06−0.04−0.0200.020.04 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 Δ𝑞Δ𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝐴 Δ𝑞Δ𝑞 (b) graph C , strange quarks −0.1−0.08−0.06−0.04−0.0200.020.040.060.08 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 𝛿𝑞𝛿𝑞 vs 𝐴 𝛿𝑞𝑞 vs 𝐵 𝛿𝑞𝛿𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝐴 𝛿𝑞𝛿𝑞 𝑚𝑦𝐴 𝛿𝑞𝑞 𝑚 |𝑦 |𝐵 𝛿𝑞𝛿𝑞 −0.09−0.08−0.07−0.06−0.05−0.04−0.03−0.02−0.0100.01 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 Δ𝑞Δ𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝐴 Δ𝑞Δ𝑞 (c) graph C , charm quarks Figure 20 . Effects of transverse (left) and longitudinal (right) polarisation for graph C . – 35 – 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 𝛿𝑞𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝑚𝑦𝐴 𝛿𝑞𝑞 −0.1−0.0500.050.1 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 Δ𝑞Δ𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝐴 Δ𝑞Δ𝑞 (a) graph C , light quarks −0.00100.0010.0020.0030.0040.0050.006 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 𝛿𝑞𝛿𝑞 vs 𝐴 𝛿𝑞𝑞 vs 𝐵 𝛿𝑞𝛿𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝐴 𝛿𝑞𝛿𝑞 𝑚𝑦𝐴 𝛿𝑞𝑞 𝑚 |𝑦 |𝐵 𝛿𝑞𝛿𝑞 −0.003−0.002−0.00100.0010.0020.0030.0040.0050.006 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑞𝑞 vs 𝐴 Δ𝑞Δ𝑞 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑞𝑞 𝐴 Δ𝑞Δ𝑞 (b) graph C , charm quarks Figure 21 . As figure 20, but for graph C . We have no strange quark results for this case. The lightquark data for A δqδq and m | y | B δqδq is very noisy and not shown for the sake of clarity. configurations available for our lattice with L = 40. As we consider only the connectedcontribution to the two-current correlation function, we restrict ourselves to the connectedgraph for the form factor as well. We fit the form factor data to a power law F u,V ( t ) = − F d,V ( t ) = (cid:0) − t/M (cid:1) − p . (4.1)We use two fit variants, which gives us a handle on the bias of the extrapolation to − t > .
15 GeV , where we have no data. Such an extrapolation bias is inevitable when we Fouriertransform from momentum to position space, as is required in (2.41). In a monopole fit, wefix p = 1 and obtain M = 777(12) MeV. Leaving the power free, we obtain p = 1 . M = 872(16) MeV. Both fits give a very good description of our lattice data, as shownin figure 14a of [12]. – 36 – .4 0.6 0.8 1.0 1.2 1.4 y [fm]0.1500.1250.1000.0750.0500.0250.000 A u d [ f m ] factorization A ud vs F u F d F u F d monopole fit F u F d p -pole fit A ud y [a] (a) y [fm]0.60.70.80.91.01.11.21.31.4 F u F d / A u d ratio F u F d / A ud monopole fit p -pole fit4 6 8 10 12 14 16 18 20 y [a] (b) Figure 22 . Test of the factorisation hypothesis (2.41) for the invariant function A ud . (a): data for A ud (restricted to the contribution from graph C ) compared with the integral over form factors onthe r.h.s. of (2.41). The form factors are determined by a monopole or a p -pole fit. (b): ratio of theform factor integral on the r.h.s. of (2.41) to the data for A ud . With the ansatz (4.1), the two-dimensional Fourier transform on the r.h.s. of (2.41) canbe carried out analytically. We compute the remaining integral over ζ numerically. Theresults obtained with the two form factor fits agree very well for y > . y range. Onemay thus say that the factorised ansatz provides a rough approximation of the two-currentcorrelator. We now investigate the combinations (3.2) of lattice graphs that appear in the matrix elementsof currents between charged or neutral pions. We omit the doubly disconnected graph D throughout, because its statistical errors are much larger than the signal for any other graph.Since data for the full set of remaining graphs is only available for light quarks, we restrict ourattention to this case. The results are shown in figures 23 and 24 for the flavour combinations ud and uu . The combinations dd and du can be obtained from the symmetry relations (2.16).As can be expected from figures 14 and 15, the statistical errors of the physical combina-tions are significantly larger than those for the connected graphs alone. Nevertheless, we seea clear negative signal for A ud in a π + . As discussed in section 4.2, this can be understoodas a dominance of the valence Fock state u ¯ d over Fock states that contain ud , ¯ u ¯ d or ¯ ud .The function A uu in a π + has a clear positive signal at small distances y . This reflects thebehaviour of graph C and corresponds to a larger probability for finding two u quarks rather– 37 – 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝑢𝑑(𝜋 + ) , 𝐴 𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] (a) A ud | π + −0.07−0.06−0.05−0.04−0.03−0.02−0.0100.010.02 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝑚 𝑦 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝑢𝑑(𝜋 + ) , 𝐴 𝛿𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] (b) myA δud | π + −0.0500.050.10.150.20.25 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝑢𝑢(𝜋 + ) , 𝐴 𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] (c) A uu | π + −0.04−0.0200.020.040.060.080.10.120.14 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝑚 𝑦 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝑢𝑢(𝜋 + ) , 𝐴 𝛿𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] (d) myA δuu | π + −0.08−0.06−0.04−0.0200.020.040.060.080.1 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝑢𝑢(𝜋 ) , 𝐴 𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] (e) A uu | π −0.08−0.06−0.04−0.0200.020.040.060.080.1 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝑚 𝑦 𝐴 𝛿 𝑞𝑞 [ f m − ] 𝑦[ fm ]𝑢𝑢(𝜋 ) , 𝐴 𝛿𝑞𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] (f) myA δuu | π Figure 23 . Twist-two functions at py = 0 for the flavour combinations ud or uu in a π + or a π .Lattice graphs are combined according to (3.2), except for of graph D , which is affected by huge errorsand hence omitted. All results are for light quarks. – 38 – 𝐴 Δ 𝑞 Δ 𝑞 [ f m − ] 𝑦[ fm ]𝑢𝑢(𝜋 + ) , 𝐴 Δ𝑞Δ𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] (a) A ∆ u ∆ u | π + −0.2−0.15−0.1−0.0500.050.1 0.4 0.6 0.8 1 1.2 1.4 1.65 10 15 20 25 𝐴 Δ 𝑞 Δ 𝑞 [ f m − ] 𝑦[ fm ]𝑢𝑢(𝜋 ) , 𝐴 Δ𝑞Δ𝑞 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] (b) A ∆ u ∆ u | π Figure 24 . Continuation of figure 23. than a u ¯ u pair at small separation y . Remarkably, the signal at small y is of comparable sizefor A uu and A ud , which implies that Fock states containing sea quarks do play an importantrole in this region. As for polarisation effects, a clear signal for ud or uu in a π + is only seenfor myA δqq and shown in the right panels of figure 23. Comparing this with A qq , we see thatspin-orbit correlations are appreciable for both flavour combinations.The flavour combination uu in a π involves the sum C + 2 C . We observe a verystrong compensation between the two connected graphs, which results in a marginal signalfor A uu and myA δuu . The twist-two functions for ud in a π receive no contribution fromconnected graphs at all. Within errors, the corresponding results are zero for all combinationsof currents, and we do not show them here.Among all polarised twist-two functions other than myA δqq , a marginally nonzero signalis only seen for the longitudinal spin correlation A ∆ u ∆ u in a π + or a π . This is dominatedby the contribution from C in both cases and shown in figure 24.Comparing the functions A ud and A uu in a π + , we see a clear difference in their y dependence. This is at variance with the assumption going into (2.7) and thus into the “pocketformula” for double parton scattering, which is that the DPDs for all parton combinationshave the same y dependence in a given hadron. Of course, the twist-two functions at py = 0are not directly related to DPDs at zero skewness ζ . However, it would be remarkable ifthe strong flavour dependence we see in πA qq ( y , py = 0) = (cid:82) dζ I qq ( y , ζ ) were absent in I qq ( y , ζ = 0). In this section, we use our data for nonzero pion momentum to study the py dependence ofthe twist-two functions. We restrict our study to graph C for light quarks on the L = 40– 39 –attice: only in this case do we have simulations for a sufficient number of pion momenta.Since graph C dominates the twist-two matrix elements for ud in a π + , we will write A ud , A δud , . . . for twist-two functions and I ud , I δud , . . . for Mellin moments in what follows. py dependence We start by proposing a functional ansatz for the twist-two functions, which is based on theirrelation (2.33) with the Mellin moments of skewed DPDs. We use this ansatz to fit the py dependence of our lattice data. This will allow for a model dependent extension of the twist-two functions to all values of py , beyond the region (2.26) available on a Euclidean lattice.This will in turn allow for a model dependent extraction of the Mellin moments of DPDsat zero skewness. For ease of notation, we write A ( y , py ) to denote any of the twist-twofunctions A ud , . . . , A δuδd , B δuδd . Likewise, we write I ( y , py ) for the Mellin moments I ud ,. . . , I δuδd , I tδuδd .The basis of our ansatz is the assumption that, in its support region − ≤ ζ ≤
1, theskewed moment I ( y , ζ ) can be approximated by a polynomial in ζ , I ( y , ζ ) = π N (cid:88) n =0 a n ( y ) ζ n (5.1)with some integer N , where we used the symmetry relation (2.30) to restrict terms to evenpowers of ζ . We write = instead of ≈ in the spirit of a fit ansatz, i.e. we do not claim that(5.1) is exact. Inverting the Fourier transform in (2.32), we obtain A ( y , py ) = N (cid:88) n =0 a n ( y ) h n ( py ) , (5.2)where we introduced the functions h n ( x ) = 12 (cid:90) − dζ e iζx ζ n . (5.3)A crucial property of the ansatz (5.2) is that its Fourier transformation (2.32) has the correctsupport in ζ .Let us collect a few properties of the functions h n ( x ). From their definition, one easilyderives h n (0) = 11 + 2 n , d h n ( x ) dx = − h n +1 ( x ) (5.4)and thus obtains the Taylor series h n ( x ) = ∞ (cid:88) m =0 ( − m n + 2 m x m (2 m )! . (5.5)– 40 –n explicit representation is given by h n ( x ) = s n ( x ) sin x + c n ( x ) cos x (5.6)with rational functions s n ( x ) = n (cid:88) m =0 (2 n )!(2 n − m )! ( − m x m , c n ( x ) = n − (cid:88) m =0 (2 n )!(2 n − m − − m x m . (5.7)For n = 0 and n = 1, these functions read s ( x ) = 1 /x , c ( x ) = 0 , s ( x ) = ( x − /x , c ( x ) = 2 /x . (5.8)In terms of the normalised quantities (cid:98) A ( y , py ) = A ( y , py ) A ( y , py = 0) , ˆ a n ( y ) = a n ( y ) A ( y , py = 0) (5.9)our ansatz (5.2) reads (cid:98) A ( y , py ) = N (cid:88) n =0 ˆ a n ( y ) h n ( py ) . (5.10)Using (2.34) and (5.5), we then obtain (cid:104) ζ m (cid:105) ( y ) = (cid:34) ( − m ∂ m (cid:98) A ( y , py )( ∂ py ) m (cid:35) py =0 = N (cid:88) n =0
11 + 2 n + 2 m ˆ a n ( y ) . (5.11)Let us now describe our general fitting procedure. In order to achieve stable fits, we firstdetermine the y dependence of A ( y , py = 0). This includes the information from data withzero pion momentum and has typically much smaller errors than the data for nonzero py .In a second step, we fit the y dependent coefficients ˆ a n ( y ) in the ansatz (5.10). To makethe degrees of freedom of this fit explicit, we consider the moments (cid:104) ζ m (cid:105) ( y ) for m = 0 , . . . , N .Inverting the relation (5.11), we obtainˆ a n ( y ) = N (cid:88) m =0 ( T − ) nm (cid:104) ζ m (cid:105) ( y ) , (5.12)where T is the ( N + 1) × ( N + 1) matrix with elements T mn = (1 + 2 n + 2 m ) − . (5.13)Since by definition (cid:104) ζ (cid:105) ( y ) = 1, we can thus fit the py dependence of the twist-two functionsto (5.10) and (5.12) with N fit parameters (cid:104) ζ (cid:105) , . . . , (cid:104) ζ N (cid:105) at each value of y . We call this“local fits” in the following, where “local” means “local in y ”.– 41 –o obtain a parametrisation of both the py and the y dependence, we assume an expan-sion (cid:104) ζ m (cid:105) ( y ) = K (cid:88) k =0 c mk (cid:112) − y k . (5.14)This is referred to as our “global fit”. By virtue of (5.10) and (5.12), this corresponds to anexpansion of (cid:98) A ( y , py ) in powers of (cid:112) − y . The condition (cid:104) ζ (cid:105) ( y ) = 1 implies c k = δ k . We recall that we have data for p = 0 , , √ , √ π/ ( La ) ≈
437 MeV. Fora given value of y , this allows for a maximum value 4 πy/ ( La ) ≈ . y/ (20 a ) for | py | . Weapply the cut (3.6) on the angle θ to the p = 0 data, but not to the data with p >
0. Wethen average all data points with the same values of py and y .We find that the twist-two functions at py = 0 can be well described by a superpositionof two exponentials, A ( y , py = 0) = A e − a ( y − y min ) + A e − a ( y − y min ) for y min ≤ y ≤ y max , (5.15)with y min = 5 a = 0 .
355 fm and y max = 20 a = 1 .
42 fm. We do not include data with y > y max ,because they have large errors and are increasingly affected by finite size effects. The resultingfit parameters are given in table 2. Let us emphasise that these fits are not suitable forextrapolating the twist-two functions to values significantly below y = y min .We notice a relatively high value of χ / dof in the fit for A δud . This is due to some scatterin the data at high y , which comes from points with large p . Repeating the fit with an upperlimit y ≤ a , we find that χ / dof decreases from 1 .
76 to 0 . A δud . By comparison, thevalue of χ / dof in the fit for A ud decreases from 0 .
95 to 0 . py dependence to (5.10) and (5.12) locally in y . To haveenough data in these fits, we introduce bins in y and combine all points with ( n − / a < function A [fm − ] a [fm − ] A [fm − ] a [fm − ] χ / dof A ud − . . . . ± . . A ∆ u ∆ d − . . . . . A δud − . . . . . A δuδd . . − . . ± . . B δuδd − . . − . . ± . . Table 2 . Parameters for the fit (5.15) of twist-two functions at py = 0 in the region 5 a ≤ y ≤ a .Throughout this section, we consider the data for graph C and light quarks on our lattice with L = 40. – 42 –unction c c [fm − ] χ / dof (cid:98) A ud . . . (cid:98) A ∆ u ∆ d − . . . (cid:98) A δud . . . (cid:98) A δuδd − . . . (cid:98) B δuδd − . . . Table 3 . Parameters of the fit of the combined y and py dependence of the normalised twist-twofunctions (cid:98) A ( y , py ) to (5.10), (5.12) and (5.14) with N = K = 1. y < ( n + 1 / a for integer n between 5 and 20. In addition, we fit the combined y and py dependence of (cid:98) A to (5.10), (5.12) and (5.14). We explored fits with different maximumvalues N and K in the sums and find that, given the fit range and the statistical quality ofour data, an adequate choice is N = 1 for local fits and N = 1, K = 1 for the global fit. Theparameters of the global fit are given in table 3. If we take N = 2 instead, the error bands ofthe fit results for (cid:98) A increase significantly, whilst the decrease of χ / dof is minor. We henceconclude that we would over-fit the data by choosing N = 2 or even higher values.We compare our data and fits in figure 25 for different functions at y = 15 a and infigure 26 for (cid:98) A ud at y = 5 a and 10 a . We find good agreement between the local and globalfits. Note that the twist-two functions are symmetric in py due to P T invariance, whichis realised on the lattice. A departure from this symmetry in the data must therefore bedue to statistical fluctuations. Many data points have admittedly large errors, which is aconsequence of at least one of y or p being large. Nevertheless, the fitted parameters for allfunctions except (cid:98) A ∆ u ∆ d are in general well determined, and the corresponding error bandsof the fit results are reasonably small. As is seen in figure 25(b), the data for (cid:98) A ∆ u ∆ d aremuch too noisy for fitting the py dependence, and we exclude this function from our furtherdiscussion.In the data for y = 15 a , we see an indication for zero crossings around | py | = 4 in severaltwist-two functions. That this can be reproduced with a superposition of the two functions h ( py ) and h ( py ) gives us some confidence in our fit ansatz.Using our fits, we can compute the moment (cid:104) ζ (cid:105) ( y ) associated with I ( y , ζ ), whichaccording to (5.11) follows from the curvature of (cid:98) A ( y , py ) at py = 0. The results are shownin figure 27. We find again good agreement between the local and global fits. A clear y dependence of (cid:104) ζ (cid:105) is observed, except for I δud . The values of (cid:104) ζ (cid:105) are not too large, especiallyfor small y . Their size does, however, imply that nonzero values of the skewness ζ must playsome role in the integral representation πA ( y , py = 0) = (cid:82) dζ I ( y , ζ ).– 43 – py A u d fit on A ud at y = 15 a , M = 1 local fitglobal fitdata 14.5 a y a data y = 15 a (a) (cid:98) A ud at y = 15 a py A u d fit on A u d at y = 15 a , M = 1 local fitglobal fitdata 14.5 a y a data y = 15 a (b) (cid:98) A ∆ u ∆ d at y = 15 a py A u d fit on A ud at y = 15 a , M = 1 local fitglobal fitdata 14.5 a y a data y = 15 a (c) (cid:98) A δud at y = 15 a py B u d fit on B u d at y = 15 a , M = 1 local fitglobal fitdata 14.5 a y a data y = 15 a (d) (cid:98) B δuδd at y = 15 a Figure 25 . Data and fits of the py dependence of normalised invariant functions. Dark points showdata at y = 15 a . Light points show data in a y range of a/ a , which are included in thelocal fits. The plot for (cid:98) A δuδd (not shown) is qualitatively similar to the one for (cid:98) B δuδd . We now use the global fit described in the last section to reconstruct the lowest Mellinmoments of skewed DPDs. Let us re-emphasise that such a reconstruction is necessarilydependent on the functional ansatz we have made, given the impossibility to constrain thefull py dependence of twist-two functions with lattice simulations. We recall that the resultsfor the spin correlation ∆ u ∆ d are too noisy and hence omitted in the following.We can easily derive the analytic form of the Mellin moments for our fits by invertingthe 2 × T mn in (5.13). This gives I ( y , ζ ) = 3 π (cid:26) − (cid:104) ζ (cid:105) ( y ) − ζ (cid:104) − (cid:104) ζ (cid:105) ( y ) (cid:105)(cid:27) A ( y , py = 0) . (5.16)– 44 – .5 1.0 0.5 0.0 0.5 1.0 1.5 py A u d fit on A ud at y = 5 a , M = 1 local fitglobal fitdata 4.5 a y a data y = 5 a (a) (cid:98) A ud at y = 5 a py A u d fit on A ud at y = 10 a , M = 1 local fitglobal fitdata 9.5 a y a data y = 10 a (b) (cid:98) A ud at y = 10 a Figure 26 . Data and fits of the py dependence of (cid:98) A ud for different y . The meaning of dark and lightpoints is as in figure 25. The values of (cid:104) ζ (cid:105) ( y ) for y ≤ a are in the range between 0 and 0 . − (cid:104) ζ (cid:105) in (5.16) is therefore always positive and varies between3 and 0 .
5. We can hence anticipate that the dependence of the Mellin moments I ( y , ζ = 0)on y and on the polarisation indices should roughly follow the corresponding dependence of A ( y , py = 0). By contrast, the coefficient of ζ in (5.16) has a larger variation and canchange sign as a function of y . Our results for the y and ζ dependence of the Mellin momentsare visualised in figures 28 and 29. Compared with the data entering our fit, we have slightlyextended the y range from 5 a down to 4 a .In the left panel of figure 30, we show the Mellin moments at ζ = 0 for the different po-larisation combinations. Comparison with the data of the corresponding twist-two functionsat py = 0 shows the close similarity between the two quantities. This corroborates the basicassumption of our discussion in section 4, namely that the qualitative features of twist-twofunctions at py = 0 are representative of the Mellin moments of ordinary DPDs.With the caveats of choosing a functional ansatz and restricting ourselves to the connectedgraph C , we can in particular extend our discussion for light quarks in section 4.3 to theMellin moments of DPDs for the flavour combination ud in a π + : there is a substantialspin-orbit correlation for one transversely polarised quark or antiquark, whereas correlationsinvolving transverse polarisation of both partons are rather small. This is one of the mainresults of our work.DPDs at ζ = 0 satisfy sum rules, which have been proposed in [47] and can be proven rig-orously in QCD [93, 94]. These sum rules express momentum and quark number conservation.– 45 – .4 0.6 0.8 1.0 1.2 1.4 y [fm]0.20.00.20.40.60.81.0 for I qq , N = 1 global fitlocal fit4 6 8 10 12 14 16 18 20 y [a] (a) I ud y [fm]0.20.00.20.40.60.81.0 for I q q , N = 1 global fitlocal fit4 6 8 10 12 14 16 18 20 y [a] (b) I δud y [fm]0.20.00.20.40.60.81.0 for I q q , N = 1 global fitlocal fit4 6 8 10 12 14 16 18 20 y [a] (c) I δuδd y [fm]0.20.00.20.40.60.81.0 for I tq q , N = 1 global fitlocal fit4 6 8 10 12 14 16 18 20 y [a] (d) I tδuδd Figure 27 . Values of the moment (cid:104) ζ (cid:105) ( y ) associated with I ( y , ζ ), extracted by local fits (datapoints) and the global fit (bands). The quark number sum rule for the flavour combination ud in a π + implies that2 π ∞ (cid:90) y cut dy y I ud ( y ; µ ) = − O (Λ y ) + O (cid:0) α s ( µ ) (cid:1) , (5.17)where Λ denotes a hadronic scale. The necessity of a lower cutoff on the y integral andthe presence of an O ( α s ) term on the r.h.s. result from the singular behaviour of DPDs atperturbatively small distances y , as explained in [50]. To avoid large logarithms in the O ( α s )term, one should take y cut ∼ /µ , and a standard choice is y cut = b /µ , where b = 2 e − γ ≈ .
12 and γ is the Euler-Mascheroni constant. With the renormalisation scale µ = 2 GeV ofour analysis, this gives y cut ≈ .
11 fm ≈ . a . Extrapolating our global fit down to this– 46 – .4 0.6 0.8 1.0 1.2 y [fm]1.00.80.60.40.20.00.20.4 I u d [ f m ] Mellin Moment I ud ( y , ), N = 1, M = 1 = 0.0= 0.25= 0.5= 0.75= 1.04 6 8 10 12 14 16 y [a] (a) I ud ( y , ζ ) y [fm]0.30.20.10.00.1 m y I u d [ f m ] Mellin Moment myI ud ( y , ), N = 1, M = 1 = 0.0= 0.25= 0.5= 0.75= 1.04 6 8 10 12 14 16 y [a] (b) myI δud ( y , ζ ) y [fm]0.10.00.10.2 I u d [ f m ] Mellin Moment I u d ( y , ), N = 1, M = 1 = 0.0= 0.25= 0.5= 0.75= 1.04 6 8 10 12 14 16 y [a] (c) I δuδd ( y , ζ ) y [fm]0.100.050.000.05 m | y | I t u d [ f m ] Mellin Moment m | y | I tu d ( y , ), N = 1, M = 1 = 0.0= 0.25= 0.5= 0.75= 1.04 6 8 10 12 14 16 y [a] (d) m | y | I tδuδd ( y , ζ ) Figure 28 . Mellin moments of skewed DPDs as a function of y , reconstructed from our global fit. value and evaluating the integral over y , we obtain2 π ∞ (cid:90) b /µ dy y I ud ( y ; µ ) = − . . (5.18)This result is not too sensitive to the extrapolation in y : taking an upper integration boundaryof 20 a , we obtain − . − . y cut , one expects a larger O (Λ y ) term on ther.h.s. of (5.17). Given the presence of this power correction in the theory prediction, we findits agreement with our result (5.18) quite satisfactory. We regard this as a strong cross checkof our analysis, and in particular of the fit ansatz we have made.– 47 – .0 0.2 0.4 0.6 0.8 1.00.250.000.250.500.751.001.251.50 r a t i o Mellin Moment I ud ( y , )/ I ud ( y , 0), N = 1, M = 1 y = 5 a = 0.35fm y = 10 a = 0.71fm y = 15 a = 1.1fm (a) I ud ( y , ζ ) /I ud ( y , r a t i o Mellin Moment I ud ( y , )/ I ud ( y , 0), N = 1, M = 1 y = 5 a = 0.35fm y = 10 a = 0.71fm y = 15 a = 1.1fm (b) I δud ( y , ζ ) /I δud ( y , Figure 29 . Mellin moments of skewed DPDs as a function of ζ , reconstructed from our global fitand normalised to their value at ζ = 0. y [fm]0.60.50.40.30.20.10.00.10.2 M e lli n m o m e n t [ f m ] Mellin Moment I aa ( y , = 0) for ud , N = 1, M = 1 I ud myI ud I u d m | y | I tu d y [a] (a) Mellin moments I ( y , ζ = 0) −0.12−0.1−0.08−0.06−0.04−0.0200.020.04 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.26 8 10 12 14 16 𝐴 [ f m − ] 𝑦[ fm ]𝐴 𝑢𝑑 vs 𝐴 𝛿𝑢𝛿𝑑 vs 𝐴 𝛿𝑢𝑑 vs 𝐵 𝛿𝑢𝛿𝑑 , 𝐶 , 𝑝 = 0 ( 𝐿 = 40 ) 𝑦[𝑎] 𝐴 𝑢𝑑 𝐴 𝛿𝑢𝛿𝑑 𝑚𝑦𝐴 𝛿𝑢𝑑 𝑚 |𝑦 |𝐵 𝛿𝑢𝛿𝑑 (b) twist-two functions A ( y , py = 0) Figure 30 . (a): Mellin moments of DPDs for the flavour combination ud in a π + , reconstructed fromour global fit. (b): Lattice data for the corresponding twist-two functions at py = 0. This shows thesame data as figure 20(a), but is limited to y ≤ a for ease of comparison. Notice that the minimum y in panels (a) and (b) is slightly different. – 48 – .4 Factorisation hypothesis for Mellin moments With the Mellin moments reconstructed from our global fit, we can also test the factorisationhypothesis (2.37), which directly follows from the corresponding hypothesis (2.6) for DPDs.To evaluate the r.h.s. of (2.37), we use the same two fits for the vector form factor of thepion as we did in section 4.4. The comparison of the left and right-hand sides of (2.37), aswell as their ratio is shown in figure 31. We see the same trend as we did in figure 22 for A ud at py = 0. At small y , the result of the factorised ansatz is too large in absolute size,and at large y it is too small. The discrepancy at large y is even somewhat stronger for theMellin moment I ud than it is for A ud . We draw the same conclusion as we did in section 4.4:the factorised ansatz for the unpolarised ud flavor combination in a π + can provide a roughapproximation at the level of several 10%. In the sense that the factorised ansatz representsthe assumption that the u and the ¯ d in a π + have independent spatial distributions, our resultfor I ud indicates that the two partons prefer to be farther apart than if they were uncorrelated. y [fm]0.70.60.50.40.30.20.10.0 M e lli n m o m e n t [ f m ] factorization I ud vs F u F d F u F d monopole fit F u F d p -pole fit I ud y [a] (a) y [fm]0.40.60.81.01.21.4 F u F d / I u d ratio F u F d / I ud monopole fit p -pole fit4 6 8 10 12 14 16 18 20 y [a] (b) Figure 31 . Test of the factorisation hypothesis (2.37) for the lowest Mellin moment I ud of theunpolarised DPD F ud at ζ = 0 in a π + . (a): Comparison of I ud determined by the global fit ofsection 5.2 with the integral over form factors on the r.h.s. of (2.37). The form factors are determinedby a monopole or a p -pole fit as described in section 4.4. (b): Ratio of the integral over form factorsand the Mellin moment. The factorisation hypothesis predicts this ratio to be 1. This paper presents the first lattice calculation that provides information about double partondistributions in a pion. Our simulations are for a pion mass of m π ≈
300 MeV, a latticespacing of a ≈ .
07 fm, and two lattice volumes with L = 32 and L = 40 points in the spatiallattice directions, respectively. We also have results for the pseudoscalar ground state madeof strange or of charm quarks at their physical masses, in a partially quenched setup.– 49 –e compute the pion matrix elements of the product of two local currents that areseparated by a space-like distance. From these tensor-valued matrix elements, we extractLorentz invariant functions associated with the twist-two operators in the definition of DPDs.In the continuum and infinite volume limits, these functions depend on the pion momentum p µ and the distance y µ between the currents only via the invariant products py and y . Thisallows us to detect discretisation and finite size effects in our data, and to devise cuts thatminimise these artefacts. In particular, most results reported here are limited to distances y above 5 a ≈ .
35 fm. Comparing the data from our two lattice volumes, we find only milddifferences in channels that have a good statistical signal. Comparing results obtained withdifferent source-sink separation, we find little evidence for contributions from excited statesin our analysis. The invariant twist-two function in the axial vector channel appears to bemost strongly affected by several of the lattice artefacts.Comparing the importance of different Wick contractions in the twist-two functions, wefind that the connected graphs C and C in figure 4 are the most important ones in almostall cases. For light quarks, graph C is as important as C at small distances y between thetwo partons, which indicates that Fock states containing sea quarks are important in thatregion. As one would expect, this importance is strongly reduced for charm quarks, but itis still visible at a level below 10%. For light quarks, the combination of graphs C and C leads to a significant difference in the y dependence of the twist-two functions for the flavourcombinations ud and uu in a π + .We compute matrix elements for different combinations of the vector, axial vector andtensor currents, which respectively correspond to unpolarised partons and partons with lon-gitudinal or transverse polarisation. For light quarks, we find surprisingly small correlationsbetween the longitudinal or transverse spins of the two partons. By contrast, a large spin-orbitcorrelation is seen between the transverse component of y and the transverse polarisation ofone of the partons. All spin correlations increase considerably with the quark mass, and forcharm quarks we observe large spin-spin correlations for both longitudinal and transversepolarisation.The invariant twist-two functions that we can determine on the lattice are not directlyrelated to the Mellin moments of DPDs, but rather to the moments of what can be called“skewed” DPDs. To compute the Mellin moments of ordinary DPDs from two-current matrixelements, one needs the dependence of the invariant functions on the variable py on the fullreal axis. This is inaccessible on a Euclidean lattice. Fitting an ansatz for the py dependenceto our lattice data, we can however reconstruct the Mellin moments by extrapolating thisansatz to the full py range. We find that the moments obtained in this way have a behaviourvery similar to the one of the twist-two functions at py = 0. A valuable cross check of ourprocedure is the fact that the result for the unpolarised Mellin moment is in good agreementwith the number sum rule that must be obeyed by the DPD for the flavour combination ud in a π + .A starting point of many phenomenological studies is the assumption that unpolarisedDPDs can be “factorised” into the single-particle distributions of each parton, which would– 50 –ean that the two partons are independent of each other. We have formalised this assumptionand tested it, both for the twist-two functions directly extracted from the lattice data andfor the Mellin moments reconstructed by extrapolating a fit to these data. In both cases, wefind that the two-parton correlator deviates from the factorisation ansatz by a few 10%, andthat the sign of the deviation depends on the transverse distance y . More specifically, thetwo partons tend to be farther apart from each other than if they were independent of eachother.We see several directions into which the studies reported here should be extended. Firstand foremost comes the extension from a pion to a nucleon, which is of direct relevance fordouble parton scattering in proton-proton collisions. Work in this direction is underway. Ona longer time scale, one will want to have simulations with finer lattice spacing and smallerquark masses. Data of sufficient quality for higher hadron momenta will extend the range in py that can be probed and thus allow for a better controlled extrapolation in this variable.Given the results obtained in the present work, we think that the efforts required for suchstudies will be rewarded with valuable physics insights. Acknowledgements
We gratefully acknowledge input from Sara Collins and Andr´e Sternbeck, and discussionswith Michael Engelhardt. We used a modified version of the Chroma [95] software pack-age, along with the locally deflated domain decomposition solver implementation of open-QCD [96]. The gauge ensembles have been generated by the QCDSF and RQCD collab-orations on the QPACE computer using BQCD [97, 98]. The graphs were produced withJaxoDraw [99, 100]. The simulations used for this work were performed with resources pro-vided by the North-German Supercomputing Alliance (HLRN). This work was supportedby the Deutsche Forschungsgemeinschaft SFB/TRR 55, the European Union’s Horizon 2020research and innovation programme under grant agreement No. 824093 (STRONG-2020) andMarie Sk(cid:32)lodowska-Curie grant agreement No. 813942 (EuroPLEx).
References [1] K. Barad, M. Ogilvie and C. Rebbi,
Quark-Antiquark Charge Distributions and Confinement , Phys. Lett. (1984) 222.[2] K. Barad, M. Ogilvie and C. Rebbi,
Quark-antiquark charge distributions , Annals Phys. (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 inLattice Mesons , Phys. Rev.
D34 (1986) 3882.[5] W. Wilcox,
Current overlap methods in lattice QCD , Phys. Rev.
D43 (1991) 2443. – 51 –
6] M. C. Chu, M. Lissia and J. W. Negele,
Hadron structure in lattice QCD. 1. Correlationfunctions and wave functions , Nucl. Phys.
B360 (1991) 31.[7] M. Lissia, M. C. Chu, J. W. Negele and J. M. Grandy,
Comparison of hadron quarkdistributions 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 hadroncorrelation functions in lattice QCD , Annals Phys. (1995) 441 [ hep-lat/9406009 ].[9] C. Alexandrou, P. de Forcrand and A. Tsapalis,
Probing hadron wave functions in latticeQCD , Phys. Rev.
D66 (2002) 094503 [ hep-lat/0206026 ].[10] C. Alexandrou, P. de Forcrand and A. Tsapalis,
The Matter and the pseudoscalar densities inlattice 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 [ ].[12] G. S. Bali, P. C. Bruns, L. Castagnini, M. Diehl, J. R. Gaunt, B. Gl¨aßle et al.,
Two-currentcorrelations in the pion on the lattice , JHEP (2018) 061 [ ].[13] M. Rinaldi, Double parton correlations in mesons within AdS/QCD soft-wall models: a firstcomparison with lattice data , .[14] A. Courtoy, S. Noguera and S. Scopetta, Two-current correlations in the pion in the Nambuand Jona-Lasinio model , .[15] P. Landshoff and J. Polkinghorne, Calorimeter Triggers for Hard Collisions , Phys. Rev. D (1978) 3344.[16] R. Kirschner, Generalized Lipatov-Altarelli-Parisi Equations and Jet Calculus Rules , Phys.Lett.
B84 (1979) 266.[17] H. Politzer,
Power Corrections at Short Distances , Nucl. Phys. B (1980) 349.[18] N. Paver and D. Treleani,
Multi-Quark Scattering and Large p T Jet Production in HadronicCollisions , Nuovo Cim.
A70 (1982) 215.[19] V. P. Shelest, A. M. Snigirev and G. M. Zinovev,
The Multiparton Distribution Equations inQCD , Phys. Lett.
B113 (1982) 325.[20] M. Mekhfi,
Multiparton processes: an application to double Drell-Yan , Phys. Rev. D (1985)2371.[21] T. Sjostrand and M. van Zijl, Multiple Parton-parton Interactions in an Impact ParameterPicture , Phys. Lett. B (1987) 149.[22] B. Blok, Y. Dokshitzer, L. Frankfurt and M. Strikman,
The Four jet production at LHC andTevatron in QCD , Phys. Rev. D (2011) 071501 [ ].[23] M. Diehl and A. Sch¨afer, Theoretical considerations on multiparton interactions in QCD , Phys. Lett.
B698 (2011) 389 [ ].[24] J. R. Gaunt and W. Stirling,
Double Parton Scattering Singularity in One-Loop Integrals , JHEP (2011) 048 [ ].[25] M. G. Ryskin and A. M. Snigirev, A Fresh look at double parton scattering , Phys. Rev.
D83 (2011) 114047 [ ]. – 52 –
26] B. Blok, Yu. Dokshitser, L. Frankfurt and M. Strikman, pQCD physics of multipartoninteractions , Eur. Phys. J.
C72 (2012) 1963 [ ].[27] M. Diehl, D. Ostermeier and A. Sch¨afer,
Elements of a theory for multiparton interactions inQCD , JHEP (2012) 089 [ ].[28] A. V. Manohar and W. J. Waalewijn, A QCD Analysis of Double Parton Scattering: ColorCorrelations, Interference Effects and Evolution , Phys. Rev.
D85 (2012) 114009 [ ].[29] A. V. Manohar and W. J. Waalewijn,
What is Double Parton Scattering? , Phys. Lett. B (2012) 196 [ ].[30] M. G. Ryskin and A. M. Snigirev,
Double parton scattering in double logarithm approximationof perturbative QCD , Phys. Rev.
D86 (2012) 014018 [ ].[31] J. R. Gaunt,
Single Perturbative Splitting Diagrams in Double Parton Scattering , JHEP (2013) 042 [ ].[32] B. Blok, Y. Dokshitzer, L. Frankfurt and M. Strikman, Perturbative QCD correlations inmulti-parton collisions , Eur. Phys. J. C (2014) 2926 [ ].[33] M. Diehl, J. R. Gaunt and K. Sch¨onwald, Double hard scattering without double counting , JHEP (2017) 083 [ ].[34] Axial Field Spectrometer collaboration, T. ˚Akesson et al.,
Double Parton Scattering in pp Collisions at √ s = 63 GeV , Z. Phys. C (1987) 163.[35] UA2 collaboration, J. Alitti et al.,
A Study of multi-jet events at the CERN anti-p p colliderand a search for double parton scattering , Phys. Lett. B (1991) 145.[36]
CDF collaboration, F. Abe et al.,
Double parton scattering in ¯ pp collisions at √ s = 1 . TeV , Phys. Rev. D (1997) 3811.[37] D0 collaboration, V. M. Abazov et al., Study of double parton interactions in diphoton + dijetevents in p ¯ p collisions at √ s = 1 . TeV , Phys. Rev. D (2016) 052008 [ ].[38] LHCb collaboration, R. Aaij et al.,
Measurement of the J/ ψ pair production cross-section inpp collisions at √ s = 13 TeV , JHEP (2017) 047 [ ].[39] ATLAS collaboration, M. Aaboud et al.,
Study of the hard double-parton scatteringcontribution to inclusive four-lepton production in pp collisions at √ s = , Phys. Lett.
B790 (2019) 595 [ ].[40]
CMS collaboration, A. M. Sirunyan et al.,
Evidence for WW production from double-partoninteractions in proton–proton collisions at √ s = 13 TeV , Eur. Phys. J. C (2020) 41[ ].[41] A. Kulesza and W. Stirling, Like sign W boson production at the LHC as a probe of doubleparton scattering , Phys. Lett. B (2000) 168 [ hep-ph/9912232 ].[42] J. R. Gaunt, C.-H. Kom, A. Kulesza and W. Stirling,
Same-sign W pair production as a probeof double parton scattering at the LHC , Eur. Phys. J. C (2010) 53 [ ].[43] F. A. Ceccopieri, M. Rinaldi and S. Scopetta, Parton correlations in same-sign W pairproduction via double parton scattering at the LHC , Phys. Rev. D (2017) 114030[ ]. – 53 –
44] S. Cotogno, T. Kasemets and M. Myska,
Spin on same-sign W -boson pair production , Phys.Rev. D (2019) 011503 [ ].[45] S. Cotogno, T. Kasemets and M. Myska,
Confronting same-sign W-boson production withparton correlations , .[46] P. Bartalini and J. R. Gaunt, eds., Multiple Parton Interactions at the LHC , vol. 29. WSP,2019, 10.1142/10646.[47] J. R. Gaunt and W. Stirling,
Double Parton Distributions Incorporating Perturbative QCDEvolution and Momentum and Quark Number Sum Rules , JHEP (2010) 005 [ ].[48] K. Golec-Biernat and E. Lewandowska, How to impose initial conditions for QCD evolution ofdouble parton distributions? , Phys. Rev. D (2014) 014032 [ ].[49] K. Golec-Biernat, E. Lewandowska, M. Serino, Z. Snyder and A. M. Stasto, Constraining thedouble gluon distribution by the single gluon distribution , Phys. Lett. B (2015) 559[ ].[50] M. Diehl, J. Gaunt, D. Lang, P. Pl¨oßl and A. Sch¨afer,
Sum rule improved double partondistributions in position space , .[51] H.-M. Chang, A. V. Manohar and W. J. Waalewijn, Double Parton Correlations in the BagModel , Phys. Rev. D (2013) 034009 [ ].[52] M. Rinaldi, S. Scopetta and V. Vento, Double parton correlations in constituent quark models , Phys. Rev. D (2013) 114021 [ ].[53] W. Broniowski and E. Ruiz Arriola, Valence double parton distributions of the nucleon in asimple model , Few Body Syst. (2014) 381 [ ].[54] M. Rinaldi, S. Scopetta, M. Traini and V. Vento, Double parton correlations and constituentquark models: a Light Front approach to the valence sector , JHEP (2014) 028 [ ].[55] W. Broniowski, E. Ruiz Arriola and K. Golec-Biernat, Generalized Valon Model for DoubleParton Distributions , Few Body Syst. (2016) 405 [ ].[56] T. Kasemets and A. Mukherjee, Quark-gluon double parton distributions in the light-frontdressed quark model , Phys. Rev. D (2016) 074029 [ ].[57] M. Rinaldi, S. Scopetta, M. C. Traini and V. Vento, Correlations in Double PartonDistributions: Perturbative and Non-Perturbative effects , JHEP (2016) 063 [ ].[58] M. Rinaldi and F. A. Ceccopieri, Relativistic effects in model calculations of double partondistribution function , Phys. Rev. D (2017) 034040 [ ].[59] W. Broniowski and E. Ruiz Arriola, Double parton distributions of the pion in the NJL model , PoS
LC2019 (2020) 031 [ ].[60] M. Rinaldi, S. Scopetta, M. Traini and V. Vento,
A model calculation of double partondistribution functions of the pion , Eur. Phys. J. C (2018) 781 [ ].[61] A. Courtoy, S. Noguera and S. Scopetta, Double parton distributions in the pion in theNambu–Jona-Lasinio model , JHEP (2019) 045 [ ].[62] W. Broniowski and E. Ruiz Arriola, Double parton distribution of valence quarks in the pionin chiral quark models , Phys. Rev. D (2020) 014019 [ ]. – 54 –
63] P. H¨agler,
Hadron structure from lattice quantum chromodynamics , Phys. Rept. (2010) 49[ ].[64] H.-W. Lin et al.,
Parton distributions and lattice QCD calculations: a community white paper , Prog. Part. Nucl. Phys. (2018) 107 [ ].[65] H.-W. Lin et al.,
Parton distributions and lattice QCD calculations: toward 3D structure , .[66] C. Zimmermann, Two-current correlations and DPDs for the nucleon on the lattice , PoS
LATTICE2019 (2020) 040 [ ].[67] M. Diehl and J. R. Gaunt,
Double parton scattering theory overview , .[68] T. Kasemets and M. Diehl, Angular correlations in the double Drell-Yan process , JHEP (2013) 121 [ ].[69] M. Diehl, J. R. Gaunt, P. Pl¨oßl and A. Sch¨afer, Two-loop splitting in double partondistributions , SciPost Phys. (2019) 017 [ ].[70] P. H¨agler, B. Musch, J. Negele and A. Sch¨afer, Intrinsic quark transverse momentum in thenucleon from lattice QCD , EPL (2009) 61001 [ ].[71] B. U. Musch, P. H¨agler, J. W. Negele and A. Sch¨afer, Exploring quark transverse momentumdistributions with lattice QCD , Phys. Rev. D (2011) 094507 [ ].[72] B. Yoon, M. Engelhardt, R. Gupta, T. Bhattacharya, J. R. Green, B. U. Musch et al., NucleonTransverse Momentum-dependent Parton Distributions in Lattice QCD: RenormalizationPatterns and Discretization Effects , Phys. Rev. D (2017) 094508 [ ].[73] X. Ji, Parton Physics on a Euclidean Lattice , Phys. Rev. Lett. (2013) 262002 [ ].[74] X. Ji,
Parton Physics from Large-Momentum Effective Field Theory , Sci. China Phys. Mech.Astron. (2014) 1407 [ ].[75] A. Radyushkin, Quasi-parton distribution functions, momentum distributions, andpseudo-parton distribution functions , Phys. Rev. D (2017) 034025 [ ].[76] K. Orginos, A. Radyushkin, J. Karpie and S. Zafeiropoulos, Lattice QCD exploration of partonpseudo-distribution functions , Phys. Rev. D (2017) 094503 [ ].[77] X. Ji, L.-C. Jin, F. Yuan, J.-H. Zhang and Y. Zhao, Transverse momentum dependent partonquasidistributions , Phys. Rev. D (2019) 114006 [ ].[78] V. M. Braun, A. Vladimirov and J.-H. Zhang, Power corrections and renormalons in partonquasidistributions , Phys. Rev. D (2019) 014013 [ ].[79] C. Alexandrou, K. Cichy, M. Constantinou, K. Hadjiyiannakou, K. Jansen, A. Scapellatoet al., Systematic uncertainties in parton distribution functions from lattice QCD simulationsat the physical point , Phys. Rev. D (2019) 114504 [ ].[80] M. A. Ebert, I. W. Stewart and Y. Zhao, Towards Quasi-Transverse Momentum DependentPDFs Computable on the Lattice , JHEP (2019) 037 [ ].[81] X. Ji, Y. Liu and Y.-S. Liu, Transverse-Momentum-Dependent PDFs from Large-MomentumEffective Theory , . – 55 –
82] K. Cichy and M. Constantinou,
A guide to light-cone PDFs from Lattice QCD: an overview ofapproaches, techniques and results , Adv. High Energy Phys. (2019) 3036904[ ].[83] X. Ji, Y.-S. Liu, Y. Liu, J.-H. Zhang and Y. Zhao,
Large-Momentum Effective Theory , .[84] R. Jaffe, Parton Distribution Functions for Twist Four , Nucl. Phys. B (1983) 205.[85] M. Diehl,
Generalized parton distributions , Phys. Rept. (2003) 41 [ hep-ph/0307382 ].[86] G. S. Bali, S. Collins, B. Gl¨aßle, M. G¨ockeler, J. Najjar, R. H. R¨odl et al.,
The moment (cid:104) x (cid:105) u − d of the nucleon from N f = 2 lattice QCD down to nearly physical quark masses , Phys. Rev.
D90 (2014) 074510 [ ].[87] G. S. Bali, S. Collins, B. Gl¨aßle, M. G¨ockeler, J. Najjar, R. H. R¨odl et al.,
Nucleon isovectorcouplings from N f = 2 lattice QCD , Phys. Rev.
D91 (2015) 054501 [ ].[88]
RQCD collaboration, G. S. Bali, S. Collins, D. Richtmann, A. Sch¨afer, W. S¨oldner andA. Sternbeck,
Direct determinations of the nucleon and pion σ terms at nearly physical quarkmasses , Phys. Rev.
D93 (2016) 094504 [ ].[89] G. S. Bali, S. Collins, A. Cox and A. Sch¨afer,
Masses and decay constants of the D ∗ s (2317) and D s (2460) from N f = 2 lattice QCD close to the physical point , Phys. Rev.
D96 (2017)074501 [ ].[90] K. Cichy, K. Jansen and P. Korcyl,
Non-perturbative renormalization in coordinate space for N f = 2 maximally twisted mass fermions with tree-level Symanzik improved gauge action , Nucl. Phys.
B865 (2012) 268 [ ].[91] M. Diehl and T. Kasemets,
Positivity bounds on double parton distributions , JHEP (2013)150 [ ].[92] J. Soffer, Positivity constraints for spin dependent parton distributions , Phys. Rev. Lett. (1995) 1292 [ hep-ph/9409254 ].[93] Gaunt, Jonathan R., Double parton scattering in proton-proton collisions , Ph.D. thesis,University of Cambridge, 2012. doi:10.17863/CAM.16589.[94] M. Diehl, P. Pl¨oßl and A. Sch¨afer,
Proof of sum rules for double parton distributions in QCD , Eur. Phys. J. C (2019) 253 [ ].[95] SciDAC, LHPC, UKQCD collaboration, R. G. Edwards and B. Joo,
The Chroma softwaresystem for lattice QCD , Nucl. Phys. Proc. Suppl. (2005) 832 [ hep-lat/0409003 ].[96] M. L¨uscher and S. Schaefer,
Lattice QCD with open boundary conditions and twisted-massreweighting , Comput. Phys. Commun. (2013) 519 [ ].[97] Y. Nakamura and H. St¨uben,
BQCD - Berlin quantum chromodynamics program , PoS
LATTICE2010 (2010) 040 [ ].[98] Y. Nakamura, A. Nobile, D. Pleiter, H. Simma, T. Streuer, T. Wettig et al.,
Lattice QCDApplications on QPACE , .[99] D. Binosi and L. Theussl, JaxoDraw: A Graphical user interface for drawing Feynmandiagrams , Comput. Phys. Commun. (2004) 76 [ hep-ph/0309015 ]. – 56 – JaxoDraw: A Graphical user interface fordrawing Feynman diagrams. Version 2.0 release notes , Comput. Phys. Commun. (2009)1709 [ ].].