Abstract
Within the theory of Quantum Chromodynamics (QCD), the rich structure of hadrons can be quantitatively characterized, among others, using a basis of universal nonperturbative functions: parton distribution functions (PDFs), generalized parton distributions (GPDs), transverse momentum dependent parton distributions (TMDs), and distribution amplitudes (DAs). For more than half a century, there has been a joint experimental and theoretical effort to obtain these partonic functions. However, the complexity of the strong interactions has placed severe limitations, and first-principle information on these distributions was extracted mostly from their moments computed in Lattice QCD. Recently, breakthrough ideas changed the landscape and several approaches were proposed to access the distributions themselves on the lattice. In this paper, we review in considerable detail approaches directly related to partonic distributions. We highlight a recent idea proposed by X. Ji on extracting quasidistributions that spawned renewed interest in the whole field and sparked the largest amount of numerical studies within Lattice QCD. We discuss theoretical and practical developments, including challenges that had to be overcome, with some yet to be handled. We also review numerical results, including a discussion based on evolving understanding of the underlying concepts and the theoretical and practical progress. Particular attention is given to important aspects that validated the quasidistribution approach, such as renormalization, matching to light-cone distributions, and lattice techniques. In addition to a thorough discussion of quasidistributions, we consider other approaches: hadronic tensor, auxiliary quark methods, pseudodistributions, OPE without OPE, and good lattice cross-sections. In the last part of the paper, we provide a summary and prospects of the field, with emphasis on the necessary conditions to obtain results with controlled uncertainties.
1. Introduction
Among the frontiers of nuclear and particle physics is the investigation of the structure of hadrons, the architecture elements of the visible matter. Hadrons consist of quarks and gluons (together called partons), which are governed by one of the four fundamental forces of nature, the strong force. The latter is described by the theory of Quantum Chromodynamics (QCD). Understanding QCD can have great impact on many aspects of science, from the subnuclear interactions to astrophysics, and, thus, a quantitative description is imperative. However, this is a very challenging task, as QCD is a highly nonlinear theory. This led to the development of phenomenological tools such as models, which have provided important input on the hadron structure. However, studies from first principles are desirable. An ideal ab initio formulation is Lattice QCD, a space-time discretization of the theory that allows the study of the properties of fundamental particles numerically, starting from the original QCD Lagrangian.
Despite the extensive experimental program that was developed and evolved since the first exploration of the structure of the proton [1, 2], a deep understanding of the hadrons’ internal dynamics is yet to be achieved. Hadrons have immensely rich composition due to the complexity of the strong interactions that, for example, forces the partons to exist only inside the hadrons (color confinement), making the extraction of information from experiments very difficult.
Understanding internal properties of the hadrons requires the development of a set of appropriate quantities that can be accessed both experimentally and theoretically. The QCD factorization provides such formalism and can relate measurements from different processes to parton distributions. These are nonperturbative quantities describing the parton dynamics within a hadron and have the advantage of being universal, that is, do not depend on the process used for their extraction. The comprehensive study of parton distributions can provide a wealth of information on the hadrons, in terms of variables defined in the longitudinal direction (with respect to the hadron momentum) in momentum space, and two transverse directions. The latter can be defined either in position or in momentum space. These variables are as follows: the longitudinal momentum fraction carried by the parton, the longitudinal momentum fraction obtained via the longitudinal momentum transferred to the hadron, and the momentum transverse to the hadron direction of movement. Parton distributions can be classified into three categories based on their dependence on , , and the momentum transferred to the hadron, , as described below.
Parton distribution functions (PDFs) are one-dimensional objects and represent the number density of partons with longitudinal momentum fraction while the hadron is moving with a large momentum.
Generalized parton distributions (GPDs) [3–7] depend on the longitudinal momentum fractions and and, in addition, on the momentum transferred to the parent hadron, . They provide a partial description of the three-dimensional structure.
Transverse momentum dependent parton distribution functions (TMDs) [8–12] describe the parton distribution in terms of the longitudinal momentum fraction and the transverse momentum . They complement the three-dimensional picture of a hadron from GPDs.
As is clear from the above classification, PDFs, GPDs, and TMDs provide complementary information on parton distributions, and all of them are necessary to map out the three-dimensional structure of hadrons in spatial and momentum coordinates. Experimentally, these are accessed from different processes, with PDFs being measured in inclusive or semi-inclusive processes such as deep-inelastic scattering (DIS) and semi-inclusive DIS (SIDIS); see e.g., [13], for a review of DIS. GPDs are accessed in exclusive scattering processes such as Deeply Virtual Compton Scattering (DVCS) [14], and TMDs in hard processes in SIDIS [10, 11]. Most of the knowledge on the hadron structure is obtained from DIS and SIDIS data on PDFs, while the GPDs and TMDs are less known. More recently, data emerge from DVCS and Deeply Virtual Meson Production (DVMP) [15]. This includes measurements from HERMES, COMPASS, RHIC, Belle and Babar, E906/SeaQuest, and the 12 GeV upgrade at JLab. A future Electron-Ion Collider (EIC), that was strongly endorsed by the National Academy of Science, Engineering and Medicine [16], will be able to provide accurate data related to parton distributions and will advance dramatically our understanding on the hadron tomography. Together with the experimental efforts, theoretical advances are imperative in order to obtain a complete picture of hadrons. First, to interpret experimental data, global QCD analyses [17–26] are necessary that utilize the QCD factorization formalism and combine experimental data and theoretical calculations in perturbative QCD. Note that these are beyond the scope of this review and we refer the interested Reader to the above references and a recent community white paper [27]. Second, theoretical studies are needed to complement the experimental program and, in certain cases, provide valuable input. This is achieved using models of QCD and more importantly calculations from first principles. Model calculations have evolved and constitute an important aspect of our understanding of parton structure. An example of such a model is the diquark spectator model [28] that has been used for studies of parton distributions (for more details, see Section 4). The main focus of the models discussed in Section 4 is the one-dimensional hadron structure (-dependence of PDFs), but more recently the interest has been extended to the development of techniques that are also applicable to GPDs and TMDs (some aspects are discussed in this review). Let us note that there have been studies related to TMDs from the lattice, and there is intense interest towards that direction (see, e.g., [29–31], and references therein).
Despite the tremendous progress in both the global analyses and the models of QCD, parton distributions are not fully known, due to several limitations: global analysis techniques are not uniquely defined [22]; certain kinematic regions are difficult to access, for instance, the very small -region [32–34]; and models cannot capture the full QCD dynamics. Hence, an ab initio calculation within Lattice QCD is crucial, and synergy with global fits and model calculations can lead to progress in the extraction of distribution functions.
Lattice QCD provides an ideal formulation to study hadron structure and originates from the full QCD Lagrangian by defining the continuous equations on a discrete Euclidean four-dimensional lattice. This leads to equations with billions of degrees of freedom, and numerical simulations on supercomputers are carried out to obtain physical results. A nonperturbative tool, such as Lattice QCD, is particularly valuable at the hadronic energy scales, where perturbative methods are less reliable, or even fail altogether. Promising calculations from Lattice QCD have been reported for many years with the calculations of the low-lying hadron spectrum being such an example. More recently, Lattice QCD has provided pioneering results related to hadron structure, addressing, for instance, open questions, such as the spin decomposition [35] and the glue spin [36] of the proton. Another example of the advances of numerical simulations within Lattice QCD is the calculation of certain hadronic contributions to the muon , for example, the connected and leading disconnected hadronic light-by-light contributions (see recent reviews of [37, 38]). Direct calculations of distribution functions on a Euclidean lattice have not been feasible due to the time dependence of these quantities. A way around this limitation is the calculation on the lattice of moments of distribution functions (historically for PDFs and GPDs) and the physical PDFs can, in principle, be obtained from operator product expansion (OPE). Realistically, only the lowest moments of PDFs and GPDs can be computed (see, e.g., [39–44]) due to large gauge noise in high moments, and also unavoidable power-divergent mixing with lower-dimensional operators. Combination of the two prevents a reliable and accurate calculation of moments beyond the second or third, and the reconstruction of the PDFs becomes unrealistic.
Recent pioneering work of X. Ji [45] has changed the landscape of lattice calculations with a proposal to compute equal-time correlators of momentum boosted hadrons, the so-called quasidistributions. For large enough momenta, these can be related to the physical (light-cone) distributions via a matching procedure using Large Momentum Effective Theory (LaMET) (see Sections 3.1 and 8). This possibility has opened new avenues for direct calculation of distribution functions from Lattice QCD and first investigations have revealed promising results [46, 47] (see Section 3.2). Despite the encouraging calculations, many theoretical and technical challenges needed to be clarified. One concern was whether the Euclidean quasi-PDFs and Minkowski light-cone PDFs have the same collinear divergence, which underlies the matching programme. In addition, quasi-PDFs are computed from matrix elements of nonlocal operators that include a Wilson line. This results in a novel type of power divergences and the question whether these operators are multiplicatively renormalizable remained unanswered for some time. While the theoretical community was addressing such issues, the lattice groups had to overcome technical difficulties related to the calculation of matrix elements of nonlocal operators, including how to obtain reliable results for a fast moving nucleon, and how to develop a nonperturbative renormalization prescription (see Section 7). For theoretical and technical challenges, see Sections 5-6. Our current understanding on various aspects of quasi-PDFs has improved significantly, and lattice calculations of quasi-PDFs have extended to quantities that are not easily or reliably measured in experiments (see Sections 9-10), such as the transversity PDF [48, 49]. This new era of LQCD can provide high-precision input to experiments and test phenomenological models.
The first studies on Ji’s proposal have appeared for the quark quasi-PDFs of the proton (see Sections 3.2 and 9). Recently, the methodology has been extended to other hadrons, in particular mesonic PDFs and distribution amplitudes (DAs). Progress towards this direction is presented in Section 10. Other recent reviews on the -dependence of PDFs from Lattice QCD calculations can be found in [27, 50, 51]. The quasi-PDFs approach is certainly promising and can be generalized to study gluon quasi-PDFs, quasi-GPDs, and quasi-TMDs. In such investigations, technical difficulties of different nature arise and must be explored. First studies are presented here. Apart from the quasidistribution approach, we also review other approaches for obtaining the -dependence of partonic functions, both the theoretical ideas underlying them (see Section 2) and their numerical explorations (Section 11).
The central focus of the review is the studies of the -dependence of PDFs. We present work that appears in the literature until November 10, 2018 (published, or on the arXiv). The discussion is extended to conference proceedings for recent work that has not been published elsewhere. The presentation is based on chronological order, unless there is a need to include follow-up work by the same group on the topic under discussion. Our main priority is to report on the progress of the field, but also to comment on important aspects of the described material based on theoretical developments that appeared in later publications, or follow-up work. To keep this review at a reasonable length, we present selected aspects of each publication discussed in the main text and we encourage the interested Reader to consult the referred work. Permission for reproduction of the figures has been granted by the Authors and the scientific journals (in case of published work).
The rest of the paper is organized as follows. In Section 2, we introduce methods that have been proposed to access the -dependence of PDFs from the lattice, which include a method based on the hadronic tensor, auxiliary quark field approaches, quasi- and pseudodistributions, a method based on OPE, and the good lattice cross-sections approach. A major part of this review is dedicated to quasi-PDFs, which are presented in more detail in Section 3, together with preliminary studies within Lattice QCD. The numerical calculations of the early studies have motivated an intense theoretical activity to develop models of quasidistributions, which are presented in Section 4. In Section 5, we focus on theoretical aspects of the approach of quasi-PDFs, that is, whether a Euclidean definition can reproduce the light-cone PDFs, as well as the renormalizability of operators entering the calculations of quark and gluon quasi-PDFs. The lattice techniques for quasi-PDFs and difficulties that one must overcome are summarized in Section 6. Recent developments on the extraction of renormalization functions related to logarithmic and/or power divergences are explained in Section 7, while Section 8 is dedicated to the matching procedure within LaMET. Lattice results on the quark quasi-PDFs for the nucleon are presented in Section 9. The quasi-PDFs approach has been extended to gluon distributions, as well as studies of mesons, as demonstrated in Section 10. In Section 11, we briefly describe results from the alternative approaches presented in Section 2. We close the review with Section 12 that gives a summary and future prospects. We discuss the -dependence of PDFs and DAs, as well as possibilities to study other quantities, such as GPDs and TMDs. A glossary of abbreviations is given in the Appendix.
2. -Dependence of PDFs
In this section, we briefly outline different approaches for obtaining the -dependence of partonic distribution functions, in particular collinear PDFs. We first recall the problem by directly employing the definitions of such functions on the lattice, using the example of unpolarized PDFs. The unpolarized PDF, denoted here by , is defined on the light cone:where is the hadron state with momentum , in the standard relativistic normalization1, the light-cone vectors are taken as , and is the Wilson line connecting the light-cone points and , while the factorization scale is kept implicit. Such light-cone correlations are not accessible on a Euclidean spacetime, because the light-cone directions shrink to one point at the origin. As discussed in the Introduction, this fact prevented lattice extraction of PDFs for many years, apart from their low moments, reachable via local matrix elements and the operator product expansion (OPE). However, since the number of moments that can be reliably calculated is strongly limited, alternative approaches were sought for to yield the full Bjorken- dependence.
The common feature of all the approaches is that they rely to some extent on the factorization framework. For a lattice observable that is to be used to extract PDFs, one can generically write the following:2where is a perturbatively computable function and is the desired PDF. In the above expression we distinguish between the factorization scale, , and the renormalization scale, . These scales are usually taken to be the same and, hence, from now on we will adopt this choice and take . Lattice approaches use different observables that fall into two classes: (1)Observables that are generalizations of light-cone functions such that they can be accessed on the lattice; such generalized functions have direct -dependence, but does not have the same partonic interpretation as the Bjorken-.(2)Observables in terms of which hadronic tensor can be written; the hadronic tensor is then decomposed into structure functions like and , which are factorizable into PDFs.
Below, we provide the general idea for several proposals that were introduced in recent years.
2.1. Hadronic Tensor
All the information about a DIS cross-section is contained in the hadronic tensor [52–56], defined bywhere is the hadron state labeled by its momentum and polarization , is virtual photon momentum, and is electromagnetic current at point . The hadronic tensor can be related to DIS structure functions and hence, in principle, PDFs can be extracted from it. is the imaginary part of the forward Compton amplitude and can be written as the current-current correlation functionwhere the subscript denotes averaging over polarizations. The approach has been introduced as a possible way of investigating hadronic structure on the lattice by K.-F. Liu and S.-J. Dong already in 1993. They also proposed a decomposition of the contributions to the hadronic tensor according to different topologies of the quark paths, into valence and connected or disconnected sea ones. In this way, they addressed the origin of Gottfried sum rule violation.
A crucial aspect for the implementation in Lattice QCD is the fact that the hadronic tensor , defined in Minkowski spacetime, can be obtained from the Euclidean path-integral formalism [54–58], by considering ratios of suitable four-point and two-point functions. In the limit of the points being sufficiently away from both the source and the sink, where the hadron is created or annihilated, the matrix element receives contributions from only the ground state of the hadron. Reconstructing the Minkowski tensor from its Euclidean counterpart is formally defined by an inverse Laplace transform of the latter and can, in practice, be carried out using, e.g., the maximum entropy method or the Backus-Gilbert method. Nevertheless, this aspect is highly nontrivial and improvements thereof are looked for. As pointed out in [59], a significant role may be played by power-law finite volume effects related to the matrix elements being defined in Euclidean spacetime. A similar phenomenon was recently observed also in the context of mixing [60]. Another difficulty of the hadronic tensor approach on the lattice is the necessity to calculate four-point correlation functions, which is computationally more intensive than for three-point functions, the standard tools of hadron structure investigations on the lattice. However, the theoretical appeal of the hadronic tensor approach recently sparked renewed interest in it [61–63]. We describe some exploratory results in Section 11.1.
2.2. Auxiliary Scalar Quark
In 1998, a new method was proposed to calculate light-cone wave functions (LCWFs) on the lattice [64]. This finds its motivation from the fact that LCWFs enter the description of many processes, such as electroweak decays and meson production. LCWF is the leading term of the full hadronic wave function in the expansion, where is the hadron momentum. For concreteness, we write the defining expression for the most studied LCWF, the one of the pion, , where is the momentum fraction:with being boosted pion state, being vacuum state, and being pion decay constant. is the Wilson line that ensures gauge invariance of the matrix element.
The essence of the idea is to “observe” and study on the lattice the partonic constituents of hadrons instead of the hadrons themselves [65, 66]. As shown in [64], the pion LCWF can be extracted by considering a vacuum-to-pion expectation value of the axial vector current with quark fields separated in spacetime. Gauge invariance is ensured by a scalar quark propagator with color quantum numbers of a quark, and at a large momentum transfer. The relation between the Fourier transform of this matrix element, computed on the lattice, and the pion LCWF, , is given by the following formula:where is momentum transfer, is scalar colored propagator, and is discrete set of partonic momentum fractions (allowed by the discretized momenta in a finite volume). The spacetime points are explained in Figure 1, which shows the three-point function that needs to be computed. The interval needs to be large to have an on-shell pion. To extract the LCWF, several conditions need to be satisfied: injected pion momentum needs to be large (to have a “frozen” pion and see its partonic constituents), the scalar quark needs to carry large energy, the time (time of momentum transfer and “transformation” of a quark to a scalar quark) has to be small (to prevent quantum decoherence and hadronization), and the lattice volume has to be large enough (to minimize effects of discretizing parton momenta). We refer to the original papers for an extensive discussion of these conditions. An exploratory study of the approach was presented in [65, 66] and later in [67], both in the quenched approximation. Naturally, the conditions outlined above are very difficult to satisfy simultaneously on the lattice, due to restrictions from the finite lattice spacing and the finite volume. However, the knowledge of the full hadronic wave function from first principles would be very much desired and further exploration of this approach may be interesting. In particular, integrals of hadronic wave functions over transverse momenta yield distribution amplitudes and PDFs.
2.3. Auxiliary Heavy Quark
In 2005, another method was proposed [58] to access hadron structure on the lattice, including PDFs. The idea relies on simulating on the lattice the Compton scattering tensor, using currents that couple physical light quarks of a hadron with a purely valence fictitious heavy quark. In the continuum limit, one can extract matrix elements of local operators in the OPE in the same renormalization scheme in which the Wilson coefficients are calculated. In this way, one gets the moments of PDFs. The crucial difference with respect to standard lattice calculations of moments is that the approach removes power divergent mixings with lower-dimensional operators, unavoidable in the lattice formulation for fourth and higher moments due to the breaking of the rotational invariance into the discrete hypercubic subgroup . This is derived and discussed in detail in [58]. Thus, in principle, any PDF moment can be extracted and the whole PDF can be reconstructed. Moreover, the heavy fictitious quark suppresses long-range correlations between the currents and also removes many higher-twist contaminations (twist=dimension−spin). The results are independent of the mass of the heavy quark, , as long as it satisfies a lattice window restriction; that is, it should be much larger than , but much smaller than the ultraviolet cutoff . In practice, this means a requirement of rather small lattice spacings.
The considered heavy-light current is defined as follows:with denoting the light quark field and the fictitious heavy quark field. The Dirac structure can be general and is typically chosen according to the desired final observable. The Euclidean Compton scattering tensor is then constructed:where are hadron states with momentum and spin and the spins are summed over. The momentum transfer between the hadron and the scattered lepton in a DIS process should be smaller or at most of the order of the heavy quark mass. Furthermore, the momenta should satisfy the constraint , where the momenta with subscript are the Minkowski counterparts of the Euclidean ones. If this condition is satisfied, analytic continuation of the hadronic tensor to Euclidean spacetime is straightforward and can be achieved by relating to . Expanding this tensor, in the continuum, using OPE, one can relate it to moments of PDFs. On the lattice, one needs to compute four-point functions to access the Compton tensor for PDFs. Analogous procedure may be applied to compute distribution amplitudes, for which the hadronic interpolating operator needs to be applied only once to the vacuum state and hence a computation of only three-point functions is required.
Numerical exploration, in the quenched approximation, is in progress [68], aimed at extracting the moments of the pion DA. Since the matching to OPE has to be performed in the continuum, at least three values of the lattice spacing need to be employed for a reliable extrapolation. Preliminary results are presented in Section 11.2 demonstrating the feasibility of this method.
2.4. Auxiliary Light Quark
Another possibility for extraction of light-cone distribution functions appeared in 2007 by V. Braun and D. Müller [69] and is based on the lattice calculation of exclusive amplitudes in coordinate space. It is similar to the fictitious heavy quark approach, but the heavy quark is replaced by an auxiliary light quark. One considers a current-current correlator:with being the electromagnetic current, but other choices of currents are also possible. For a discussion on extracting partonic distributions at light-like separations from such Euclidean correlators, we refer to Section II of [70]. On the lattice, can be computed as an appropriate ratio of a three-point and a two-point function. If the separation between currents is small, the correlator can be computed perturbatively (using OPE) and in such a case Equation (9) yields:at leading order and leading twist. Equation (10) is proportional to the Fourier transform, , of the pion DA, , where is the quark momentum fraction. The perturbative expression for the correlator was also derived in [69] to NNLO, including twist-4 corrections. The LO and leading twist expression for the case of scalar-pseudoscalar densities in Equation (9) was given in [71]. It has been emphasized that the pion boost plays a different role than in some other approaches, as it does not suppress higher-twist contributions, but rather enters the Ioffe time . Thus, going to large boosts is important to have the full information on the coordinate space pion DA, , which can allow disentanglement between phenomenological models considered in the literature, that disagree in the regime of large Ioffe times. Advantages of the approach include the possibility of having arbitrary direction of with respect to the boost direction, which may make it possible to minimize discretization effects. Moreover, one avoids complications related to the renormalization in the presence of a Wilson line (see Section 7); that is, one only needs renormalization of standard local operators which is at most logarithmically divergent. Finally, different possible Dirac structures may give the possibility of better control of higher-twist contamination. Obviously, the approach can also be generalized to extract PDFs, which, however, would necessitate the computation of four-point functions (see also Section 2.8).
The first numerical investigation of this approach is under way by the Regensburg group [70, 71] and is aimed at computing the pion DA, using multiple channels. The results fully prove the feasibility of this method and establish its status as a promising way of studying hadron structure; see also Section 11.3. Nevertheless, the requirement of calculation of four-point functions for extracting PDFs may prove to be a serious restriction and an exploratory study for, for example, nucleon PDFs, is not yet available.
2.5. Quasidistributions
In 2013, X. Ji proposed a new approach to extracting the -dependence of structure functions [45]. Although historically it was not the first idea, it can be presently judged that it has been a breakthrough in the community’s thinking about -dependence from numerical simulations on a Euclidean lattice. In particular, it clearly renewed the interest also in approaches proposed earlier and described above. Ji’s approach, obviously, bears similarities with the earlier methods and is also based on the factorization framework, in which a lattice computable function is factorized into a hard coefficient and a nonperturbative object like a PDF or a DA. The main difference is another type of object that is used to connect a quark and an antiquark separated by some distance and that ensures gauge invariance. In earlier proposals, different types of auxiliary quark propagators were used for this: scalar, heavy, or light quark propagators. In Ji’s technique, this role is played by a Wilson line, that is, the same object that is used in definitions of PDFs and other distribution functions. Thus, in general, the quasidistribution approach is the closest transcription of a light-cone definition to Euclidean spacetime, effectively boiling down to replacing light-cone correlations by equal-time correlators along the direction of the Wilson line.
We illustrate the idea using the example of PDFs, while analogous formulations can be used to define DAs, GPDs, etc. It is instructive to see the direct correspondence between the light-cone definition (Equation (1)) and the definition of quasi-PDFs. As pointed out above, since light-cone correlations cannot be accessed on a Euclidean lattice, Ji proposed to evaluate on the lattice the following distribution, now termed the quasidistribution:where , is the quark momentum in the -direction, and is the Wilson line in the boost direction3. The light-cone definition corresponds to the above expression at infinite momentum boost, in line with Feynman’s original parton model [72, 73]. Since the momentum of the nucleon on the lattice is obviously finite, the partonic interpretation is formally lost and some quarks can carry more momenta than the whole nucleon () or move in the opposite direction to it ().
The quasidistribution differs from the light-cone one by higher-twist corrections suppressed with and , where is the nucleon mass; see Section 3.1 for more details. A vital observation of Ji was that the difference between the two types of distributions arises only in the UV region; that is, their structure in the IR is the same. This means that the UV difference can be computed in perturbation theory and subtracted from the result, which comes under the name of matching to a light-cone distribution or Large Momentum Effective Theory (LaMET) [74]. The possibility of correcting the higher-twist effects by LaMET is an important difference with respect to previously mentioned approaches. However, explicit computation of such effects is also possible in them, as demonstrated already in the original paper for the auxiliary light quark approach [69].
The quasidistribution approach received a lot of interest in the community and sparked most of the numerical work among all the direct -dependence methods. In further sections, we discuss in more detail its various aspects and the plethora of numerical results obtained so far.
2.6. Pseudodistributions
The approach of quasidistributions was thoroughly analyzed by A. Radyushkin [75–77] in the framework of virtuality distribution functions introduced by the same Author [78, 79] and straight-link primordial TMDs. In the process, he discovered another, but strongly related, type of distribution that is accessible on the lattice and can be related to light-cone distributions via factorization. It can be extracted from precisely the same matrix element that appears in Equation (1), , viewed as a function of two Lorentz invariants, the “Ioffe time” [52], and . Thus, has been termed the Ioffe time distribution (ITD). As in Ji’s approach the vector can be chosen to be purely spatial, on a Euclidean lattice. Then, one defines a pseudodistribution:Thus, the variation with respect to a quasi-PDF is the Fourier transform that is taken over the Ioffe time (at fixed ), as opposed to being over the Wilson line length (at fixed momentum ). A consequence of this difference is that pseudo-PDFs have considerably distinct properties from quasi-PDFs. In particular, the distribution has the canonical support, .
We briefly mention here the issue of power divergences induced by the Wilson line, to be discussed more extensively in Sections 5.2.1 and 7. In the pseudodistribution approach, a convenient way of eliminating these (multiplicative) divergences is to take the ratio [77, 80]. The reduced ITD, , can then be perturbatively matched to a light-cone Ioffe time PDF [81–86], as demonstrated in Section 8. The (inverse) length of the Wilson line plays the role of the renormalization scale and can be related to, e.g., the scale.
Numerical investigation of the pseudodistribution approach has proceeded in parallel with the theoretical developments and promising results are being reported [81, 87–90] (see also Section 11.4).
2.7. OPE without OPE
Yet another recent proposal to compute hadronic structure functions was suggested in [91]. It is closely related to known ideas introduced around 20 years ago, dubbed “OPE without OPE” by G. Martinelli [92] and applied in, e.g., flavor physics [93]. The name originates from the fact that one directly computes the chronologically ordered product of two currents rather than matrix elements of local operators. In addition, one works in the regime of small spacetime separations between currents (to use perturbation theory to determine the expected form of the OPE), but large enough to avoid large discretization effects. The idea is also an ingredient of the proposal to compute LCWFs with the aid of a fictitious scalar quark [64].
The starting point is the forward Compton amplitude of the nucleon, defined similarly as in Equation (3). It can be decomposed in terms of DIS structure functions and . With particular choice of kinematics, one can obtain the following relations between the -component of the Compton amplitude and :where . Being able to access for large enough number of values of , one can extract the moments of or even the whole function.
Another important ingredient of the method proposed in [91] is the efficient computation of , that is, one that avoids the computation of four-point functions. It relies on the Feynman-Hellmann relation [94]. One extends the QCD Lagrangian with a perturbationwhere is the electric charge of the -th flavor and is a parameter with dimension of mass. Evaluating the derivative of the nucleon energy with respect to , which requires dedicated simulations at a few values, leads to estimates of :The Authors also showed first results obtained in this framework and point to directions of possible improvements and to prospects of computing the entire structure function based on this method (see Section 11.5).
2.8. Good Lattice Cross Sections4
A novel approach to extracting PDFs or other partonic correlation functions from ab initio lattice calculations was proposed by Y.-Q. Ma and J.-W. Qiu [85, 86, 95]. They advocate for a global fit of “lattice cross-sections” (LCSs), i.e., appropriate lattice observables defined below, to which many of the ones described above belong. The logic is that standard phenomenological extractions of PDFs rely on an analogous fit to hadronic cross-sections (HCSs) obtained in experiments and a global fit approach can average out some of the systematics and yield ultimately good precision.
Good LCSs, i.e., ones that can be included in such a global fit, are the ones that have the following properties:(1)They are calculable in Euclidean Lattice QCD(2)Have a well-defined continuum limit(3)Have the same and factorizable logarithmic collinear divergences as PDFs.
All of these properties are crucial and nontrivial. The first one excludes the direct use of observables defined on the light cone. In practice, the second one requires the observables to be renormalizable. Finally, the third property implies that the analogy with global fits to HCSs is even more appropriate; both strategies need to rely on the factorization framework: LCSs and HCSs are then written as a convolution of a perturbatively computable hard coefficient with a PDF.
Ma and Qiu constructed also a class of good LCSs in coordinate space that have the potential of being used in the proposed global fits, demonstrating that the three defining properties of LCSs are satisfied [86]. The considered class is very closely related to the one proposed by Braun and Müller (see Section 2.4), but the latter Authors concentrated on the pion DA, while the analysis of Ma and Qiu deals with the case of hadronic PDFs. In general, the relevant matrix element can be written aswhere stands for different possible operators that can be shown to be factorizable into the desired PDF. is the hadron momentum, and is the largest separation of fields from which the -th operator is constructed (), . One suggested choice for are the current-current correlators:where stands for the dimension of the renormalized current , with being the renormalization function of the current . Different possible options for the currents were outlined and, then, factorization was demonstrated for this whole class of LCSs. Renormalizability of these objects is straightforward, as they are constructed from local currents. Also, the feasibility of a lattice calculation is easy to establish if has no time component. Thus, this class of matrix elements belongs to the set of good LCSs. It was also shown in [86] that three of the observables discussed above, quasi-PDFs, pseudo-PDFs, and the Compton amplitude , are also examples of good LCSs.
An explicit numerical investigation of the current-current correlators is in progress by the theory group of Jefferson National Laboratory (JLab) and first promising results for pion PDFs, using around 10 different currents, have been presented. For more details see Section 11.6.
3. Quasi-PDFs: More Details and Early Numerical Studies
We discuss now, in more detail, the quasidistribution approach which is the main topic of this review. The focus of this section is on the theoretical principles of this method and we closely follow the original discussion in Ji’s first papers. Since these were soon followed by numerical calculation within Lattice QCD exploring the feasibility of the approach, we also summarize the progress on this side. We also identify the missing ingredients in these early studies and aspects that need significant improvement.
3.1. Theoretical Principles of Quasi-PDFs
Ji’s idea of quasi-PDFs [45] relies on the intuition that if light-cone PDFs can be equivalently formulated in the infinite momentum frame (IMF), then the physics of a hadron boosted to a large but finite momentum has to have much in common with the physics of the IMF. Moreover, the difference between a large momentum frame and the IMF should vanish when the hadron momentum approaches infinity. These intuitions were formalized by Ji in his original paper and we reproduce here his arguments.
Consider a local twist-2 operatorwhere parentheses in superscript indicate symmetrization of indices and the subtracted trace terms include operators of dimension with at most Lorentz indices. The matrix element of such an operator in the nucleon state readswhere is a symmetric rank- tensor [96] and the coefficients are moments of PDFs, i.e., with even . Taking all indices , one recovers the light-cone, time-dependent correlation that defines the PDF. We now consider a different choice of indices, without any temporal component, :with the trace terms containing operators with again at most Lorentz indices. Because of Lorentz invariance, matrix elements of the trace terms in the nucleon state are at most multiplied by . On the other hand [96],where is a combinatorial coefficient and is the nucleon mass. As a consequence, we find thatThe form of this expression implies that using an operator with the Wilson line in a spatial direction, in a nucleon state with finite momentum, leads to the light-cone PDF up to power-suppressed corrections in the inverse squared momentum. The corrections are of two kinds: generic higher-twist corrections and ones resulting from the nonzero mass of the nucleon. As we will discuss below, the latter can be calculated analytically and subtracted out. However, the former can only be overcome by simulating at a large enough nucleon boost and by using a matching procedure.
In the original paper that introduced the quasidistribution approach [45], Ji pointed out an intuitive way to understand the above result: “(…) consider the Lorentz transformation of a line segment connecting with the origin of the coordinates. As the boost velocity approaches the speed of light, the space-like line segment is tilted to the light-cone direction. Of course, it cannot literally be on the light-cone because the invariant length cannot change for any amount of boost. However, this slight off-light-cone-ness only introduces power corrections which vanish asymptotically.” This intuition is schematically represented in Figure 2.
We turn now to discussing how to match results obtained on the lattice, with a hadron momentum that is finite and relatively small, to the IMF. The subtlety of this results from the fact that regularizing the UV divergences does not commute with taking the infinite momentum limit. When defining PDFs, the latter has to be taken first, i.e., before removing the UV cutoff, whereas on the lattice one is bound to take all scales, including the momentum boost of the nucleon, much smaller than the cutoff, whose role is played by the inverse lattice spacing. To overcome this difficulty, one needs to formulate an effective field theory, termed Large Momentum Effective Theory (LaMET) [74], which takes the form of matching conditions that take the quasidistribution to the IMF, or light-cone, distribution. LaMET is an effective theory of QCD in the presence of a large momentum scale , in a similar sense as Heavy Quark Effective Theory (HQET) [97] is an effective theory of QCD in the presence of a heavy quark, that can have a mass larger than the lattice UV cutoff.
The parallels of LaMET with HQET are more than superficial. We again follow Ji’s discussion [74]. In HQET, a generic observable depends on the heavy mass and a cutoff . The matching with an observable defined in the effective theory, in which the heavy quark has infinite mass, can be written in the following way, due to asymptotic freedom:where is renormalized at a scale in the effective theory. Additionally, renormalization of the full theory translates the cutoff scale to a renormalization scale . The crucial aspect is that and have the same infrared physics. Thus, the matching coefficient, , is perturbatively computable as an expansion in the strong coupling constant. Apart from the perturbative matching, there are power-suppressed corrections, which can also be calculated.
Using the same ideas, one can write the relation between an observable in the lattice theory, , dependent on the analogue of a heavy mass, i.e., a large momentum (and on the cutoff scale), and an observable in a theory in the IMF, , thus corresponding to Feynman’s parton model or to a light-cone correlation. This is again valid because of asymptotic freedom. The matching reads as follows:We have, therefore, established the close analogy between HQET and the IMF parton model and the latter plays the role of an effective theory for a nucleon moving with a large momentum, just as HQET is an effective theory for QCD with a heavy mass. The infrared properties are, again, the same in both theories and the matching coefficient, , can be computed in perturbation theory. There are power-suppressed corrections in inverse powers of , vs. inverse powers of in HQET.
To summarize, the need for LaMET when transcribing the finite boost results to light-cone parton distributions is the consequence of the importance of the order of limits. Parton physics corresponds to taking in the observable first, before renormalization. On the lattice, in turn, UV regularization is necessarily taken first, before the infinite momentum limit, since no scale in the problem can be larger than the UV cutoff. However, interchanging the order of limits does not influence infrared physics and, hence, only matching in the ultraviolet has to be carried out and can be done perturbatively. The underlying factorization can be proven order by order in perturbation theory. It is important to emphasize that any partonic observable can be accessed within this framework, with the same universal steps:(1)Construction of a Euclidean version of the light-cone definition. The Euclidean observable needs to approach its light-cone counterpart in the limit of infinite momentum(2)Computation of the appropriate matrix elements on the lattice and renormalizing them(3)Calculation of the matching coefficient in perturbation theory and use of LaMET, Equation (24), to extract the light-cone distribution.
There is complete analogy also with accessing parton physics from scattering experiments, using factorization theorems and, thus, separating the nonperturbative (low-energy) and perturbative (high-energy) scales. To have similar access to partonic observables from lattice computations, LaMET plays the role of a tool for scale separation. Moreover, just as parton distributions can be extracted from a variety of different scattering processes, they can also be approached with distinct lattice operators.
We continue the discussion of LaMET by considering now the matching process in more detail. In the first paper devoted to the matching in the framework of LaMET, the nonsinglet PDF case was discussed [98]. We remind here the definition of the quasi-PDF:taking the original choice of the Dirac structure, i.e., , for the unpolarized case (see discussion about mixing for certain Dirac structures in Section 7). The matching condition should take the formwhere the quasi-PDF, , is renormalized at a scale . The calculation of the matching is performed in a simple transverse momentum cutoff scheme, regulating the UV divergence, and, later in Section 8, we will consider further developments, including matching from different schemes to the scheme. The motivation behind using the transverse momentum cutoff scheme is to take trace of the linear divergence related to the presence of the Wilson line, which would not be possible when using dimensional regularization.
The tree level of both the quasi- and the light-cone distributions is the same, i.e., a Dirac delta . At one-loop level, two kinds of contributions appear: the self-energy diagram (left one in Figure 3) and the vertex diagram (right one in Figure 3). The quasidistribution receives, hence, the following one-loop correction: where are one-loop self-energy corrections (wave function corrections) and are the one-loop vertex corrections. Expressions for an explicit form of and are given in [98]. Their crucial aspect is that they are nonzero not only in the canonical range , but also outside of it, for any positive and negative . This corresponds to the loss of the standard partonic interpretation mentioned above. An important aspect is the particle number conservation, . Different kinds of singularities appear:(i)Linear (UV) divergences due to the Wilson line, taking in this scheme the form (ii)Collinear (IR) divergences, only in , expected to be the same as in the light-cone distribution(iii)Soft (IR) divergences (singularities at ), canceling between the vertex and the self-energy corrections (“plus prescription”)(iv)Logarithmic (UV) divergences in self-energy corrections, regulated with another cutoff5.
We turn now to the light-cone distribution. It can be calculated in the same transverse momentum cutoff scheme by taking the limit à la Weinberg [99]. We do not write here the final formulae, which can be found again in [98]. The result is the same as obtained from the light-cone definition. Crucially, the collinear divergence is the same as in the quasi-PDF, as anticipated based on physical arguments that the construction of quasidistributions should not modify the IR properties. Obviously, the diagrams in this case are nonzero only for ; i.e, has a partonic interpretation.
Having computed the one-loop diagrams, one is ready to calculate the matching coefficient in Equation (26). Its perturbative expansion can be written as with the following one-loop function: Note that the matching process effectively trades the dependence on the large momentum for renormalization scale dependence (the term with the logarithm of ), another characteristic feature of effective field theories. The antiquark distribution and the -factor satisfy ; hence including antiquarks is straightforward. Similar matching formulae were also derived for the case of helicity and transversity distributions [98].
The early papers of [45, 74, 98] provided the systematic framework for defining quasidistributions and matching them with their light-cone counterparts. Since then, there have been several improvements of many aspects of this programme, including renormalization, matching, target mass corrections, and other theoretical aspects, as well as developments for distributions other than the nonsinglet quark PDF of the nucleon discussed here. Before we turn to them, we report the early numerical efforts in Lattice QCD that illustrate the state-of-the-art calculations of that time.
3.2. Early Numerical Investigations
Ji’s proposal for a novel approach of extracting partonic quantities on the lattice, in particular PDFs, sparked an enormous wave of interest, including numerical implementation and model investigations (see Section 4).
The first lattice results were presented in 2014 in [46] by H.-W. Lin et al. and later in [47, 100] by the ETM Collaboration6. Lin et al. used a mixed action setup of clover valence quarks on a HISQ sea, lattice volume , fm, and pion mass () around 310 MeV, while ETMC used a unitary setup with maximally twisted mass quarks, lattice volume , fm, and MeV. Both papers implemented the bare matrix elements of the isovector unpolarized PDF ( flavor structure, Dirac structure ). The statistics for Lin et al. is 1383 measurements, while ETMC used a larger statistics of 5430 measurements. The employed nucleon boosts were in both cases the three lowest multiples of , i.e., 0.43, 0.86, and 1.29 GeV (Lin et al.) and 0.47, 0.94, and 1.42 GeV (ETMC), with noticeable increase of noise for the larger boosts, resulting in larger statistical errors. In view of the missing renormalization programme, both collaborations used HYP smearing [101] to bring the renormalization functions closer to their tree-level values (ETMC also applied the renormalization factor to correctly renormalize the local matrix element, i.e., one without the Wilson line). ETMC presented a study of the bare matrix elements dependence on the number of HYP smearing iterations, finding large sensitivity to this number especially for the imaginary part (the matrix elements are real only in the local case). Furthermore, ETMC tested the contamination by excited states by using two source-sink separations () of fm and fm, finding compatible results, but within large uncertainties. The source-sink separation in the study of Lin et al. was not reported. We note that separations below 1 fm are more susceptible to excited states contamination. However, the goal of these preliminary studies is to explore the approach of quasi-PDFs, postponing the investigation of excited states for later calculations. Having the bare matrix elements, the Fourier transform was taken to obtain the corresponding quasi-PDFs. The quasi-PDFs were matched to light-cone PDFs using the formulae of [98] and nucleon mass corrections were also applied. The obtained final PDFs are shown in Figure 4 for each study. One observes a similar picture from both setups and certain degree of qualitative agreement with phenomenological PDFs [102–104], shown for illustration purposes. Lin et al. also computed the helicity PDF (Dirac structure in the matrix elements) and quoted the value of the sea quark asymmetry, but without showing the quasi- or final distributions.
The two earliest numerical investigations of Ji’s approach showed the feasibility of lattice extraction of PDFs. However, they also identified the challenges and difficulties. On one side, these were theoretical, like the necessity of development of the missing renormalization programme and the matching from the adopted renormalization scheme to the desired scheme. On the other side, it became also clear that the computation is technically challenging, in particular because of the decreasing signal-to-noise ratio when increasing the nucleon boost. The computational cost also increases with the source-sink separation, for which a large value (typically above 1 fm) is needed to suppress excited states. In addition, full control over typical lattice systematics, e.g., cutoff effects, finite volume effects, or the pion mass dependence, was also missing. At this stage, some difficulties were still unidentified, for example, the mixing between certain Dirac structures due to the chiral symmetry breaking in the used lattice discretizations, first identified by Constantinou and Panagopoulos [105, 106].
Further progress was reported in the next two papers by the same groups (with new members), early in 2016 by Chen et al. [107] and later in the same year by Alexandrou et al. (ETMC) [108]. Both groups used the same setups as in [46, 100] but implemented a number of improvements and considered all three types of collinear PDFs: unpolarized, helicity, and transversity. Chen et al. [107] considered two source-sink separations, fm and fm, and performed measurements on 449 gauge field configuration ensembles with 3 source positions on each configuration, using the same set of nucleon momenta as in [46], 0.43, 0.86, and 1.29 GeV. They also derived and implemented nucleon mass corrections (NMCs, also called target mass corrections, TMCs7) for all three cases of PDFs. The NMCs will be discussed below in Section 6. In the work of ETMC [108], a large-statistics study was performed with 30000 measurements for each of the three momenta, 0.47, 0.94, and 1.42 GeV, at an increased source-sink separation of fm. In the course of this work, the method of momentum smearing was introduced [109] (see Section 6 for details) to overcome the difficulty of reaching larger nucleon boosts. The technique was implemented by ETMC and results were presented for additional momenta, 1.89 and 2.36 GeV with small statistics of 150 and 300 measurements, respectively. Moreover, a test of compatibility between standard Gaussian smearing, applied in the earlier work of both groups, and the momentum smearing was performed at GeV for the unpolarized case. This revealed a spectacular property that similar statistical error as for Gaussian smearing with 30000 measurements can be obtained with only 150 measurements employing momentum smearing.
As an illustration, we show the final helicity PDFs in Figure 5. Direct visual comparison between the two results is not possible, since the plot by Chen et al. shows the PDF multiplied by . Nevertheless, the qualitative picture is similar, revealing that no striking differences occur due to different lattice setups. The much smaller uncertainty in the plot by ETMC results predominantly from over 20 times larger statistics. Analogous plots for the unpolarized and transversity cases can be seen in [107, 108].
This concludes our discussion of the early explorations of the quasi-PDF approach. References [46, 100, 107, 108] proved its feasibility on the lattice and initiated identification of the challenges, already mentioned above. Further progress was conditioned on theoretical and practical improvements that will be described in later sections.
4. Quasidistributions: Model Investigations
Apart from theoretical analyses and numerical investigations in Lattice QCD, insights about the quasidistribution approach were obtained also from considerations in the framework of phenomenological or large- models. In this section, we take a closer look at quasi-PDFs, quasi-GPDs, and quasi-DAs in such models and review the conclusions obtained in setups where direct access to analytical forms of quasi- and light-cone distributions is possible.
4.1. Diquark Spectator Model
The aim of the early (2014) work of L. Gamberg et al. [110] and I. Vitev et al. [111] was to provide guidance about nucleon momenta needed for a reliable approach to the light-cone PDFs, for all collinear twist-2 distributions, i.e., unpolarized, helicity, and transversity PDFs. The Authors considered the diquark spectator model (DSM) [28], a phenomenological model that captures many of the essential features of the parton picture. The central idea of the DSM is to insert a completeness relation with intermediate states in the operator definition of PDFs (or some quark-quark correlation functions) and then truncate them to a single state with a definite mass. Such state is called the diquark spectator. This procedure boils down to making a particular ansatz for the spectral decomposition of the considered observable. The diquark spectator, in the simplest picture, can have spin-0 (scalar diquark) or spin-1 (axial-vector diquark). Finally, the nucleon is viewed as a system consisting of a constituent quark of some mass and a scalar or axial-vector diquark. The basic object in this approximation is the nucleon-quark-diquark interaction vertex, which contains a suitably chosen form factor, taken in the so-called dipolar form in [110].
With such setup, one can derive the model expressions for all kinds of collinear quasi-PDFs, combining the expressions for scalar and axial-vector diquarks. The obtained relations can be used to study the approach to the light-cone PDFs, also calculated in the DSM. Gamberg et al. [110] got 3 couplings, fixed by normalization and 9 parameters of the model that were fixed by fitting to experimental data, with good quality of fits. Then, they considered the quasi-PDFs for different boosts from 1 to 4 GeV. It was found that the shape of quasi-PDFs approaches the PDF for GeV. The agreement is especially good in the small to intermediate- regime, while large- needs significantly larger boost for a satisfactory agreement. The Authors also studied the Soffer inequality [112], stating that the transversity distribution should not be larger than the average of unpolarized and helicity ones. It holds for the standard PDFs. For quasi-PDFs, it was found that the inequality is always satisfied for the quark, while it is violated for the quark in the entire range of for small momenta of around 0.5 GeV.
Further model study of quasi-PDFs was presented in [113] in 2016. Bacchetta et al. confirmed the conclusions of [110] and, motivated by the conclusion that the large- region of quasi-PDFs converges much more slowly to the appropriate light-cone PDF, they proposed a procedure to aid the large- extraction of PDFs from the quasidistribution approach. This is a relevant aspect for computations of quasi-PDFs in Lattice QCD. The main idea is to combine the result of a quasi-PDF and that of the corresponding moments of PDFs. One divides the whole -region into two intervals, with a “matching” point . For , one assumes that the computed quasi-PDF is already a good approximation to the standard PDF. In turn, for , a parametrization is used, with parameters fixed by conditions of smoothness at (for the value of the quasi-PDF and for its derivative) and by the available moments of PDFs. The procedure to reconstruct unpolarized and helicity PDFs was tested numerically in the DSM for two matching points, , and two nucleon momenta, GeV. In all cases, there is significant improvement of agreement with respect to the standard PDF, especially at . Excellent agreement was observed for the quark even for the lower nucleon momentum, while the quark seems to require a larger value of , which is due to the worse agreement of the quasi-PDF and the standard PDF at the matching point. Overall, the procedure was proven to be successful in the DSM and the Authors are hopeful that it can be effectively applied also for actual Lattice QCD data.
The DSM (with scalar diquarks) was also employed as a framework for studying quasi-GPDs [114] in 2018. S. Bhattacharya, C. Cocuzza, and A. Metz calculated twist-2 unpolarized quasi-GPDs (the so-called and functions; for some definitions of GPDs variables and functions, see Section 8.2), using two Dirac structures, and , motivated by the discovery of mixing for one of them [105, 106]. They verified that, in the forward limit, their expressions reduce to the ones for the respective quasi-PDFs and, in the infinite momentum limit, their quasi-GPDs approach the appropriate GPDs. They found that all results for quasidistributions are continuous and argued that this feature should hold also at higher twist in the DSM. The Authors also checked the reliability of the cut-diagram approach, widely used in spectator models, and concluded it does not reproduce certain terms appearing from handling the calculation exactly. Thus, this approach is a simplification that should be avoided when dealing with quasidistributions. Having the final analytical expressions for quasi-GPDs and quasi-PDFs, they studied numerically the approach to infinite nucleon boost. They found that of order 2 GeV and larger yields quasifunctions within of their light-cone counterparts in a wide range of . The problematic region, as for quasi-PDFs, is the large regime, and the discrepancies increase for larger skewness . Interestingly, the derivation of matching for GPDs [115], described shortly in Section 8.2, indicates that no matching is required for the function (at leading order). However, in this model study, no significant differences in the convergence of the and functions were seen. In the ERBL region, , the agreement with standard GPDs is good, provided that is not too small. The results for both Dirac structures were found to be very similar at large enough momentum ( GeV). To verify and strenghten the conclusions, the Authors also checked the sensitivity to parameter variations (constituent mass and spectator mass) and found no significant differences. As numerical exploration of quasi-GPDs on the lattice is still missing, the DSM results can provide very useful guidance to such attempts.
4.2. Virtuality Distribution Functions
A model investigation of quasi-PDFs was performed also by A. Radyushkin in 2016-17 [75, 76]. He used his formalism of virtuality distribution functions (VDFs) [78, 79]. In the VDF formalism, a generic diagram for a parton-hadron scattering corresponds to a double Fourier transform of the VDF, , where is related to the parton virtuality (giving the name to the VDF) and is the hadron mass. The variables conjugates in the double Fourier transform are ( is hadron momentum and is separation of fields) and . The VDF representation holds for any and , but the case relevant for PDFs is the light-cone projection, . Then, one can define primordial (straight-link) TMDs and derive relations between VDFs, TMDs, and quasi-PDFs. Working in a renormalizable theory, one can represent the VDF as a sum of a soft part, i.e., generating a nonperturbative evolution of PDFs, and a hard tail, vanishing with the inverse of .
Numerical interest in these papers was in the investigation of the nonperturbative evolution generated by the soft part of the VDF or, equivalently, the soft part of the primordial TMD. Radyushkin considers two models thereof, with a Gaussian-type dependence on the transverse momentum (“Gaussian model”) and a simple non-Gaussian model (“ model”). These models are two extreme cases of a family of models, one with a too fast and one with a too slow fall-off in the impact parameter. In the numerical part of [75], the formalism was applied to a simple model PDF, . Both TMD models give similar evolution patterns, implying that one observes some universal features related to the properties of quasi-PDFs. It was also observed that the approach to the limiting PDF is not uniform for different and it can even be nonmonotonic for small nucleon momenta. These conclusions can provide very useful guidance to Lattice QCD calculations, meaning, e.g., that simple extrapolations in the inverse squared momentum might not be justifiable.
In the work of [76], Radyushkin considered target mass corrections (TMCs) in quasi-PDFs using the same framework. In both TMD models, it was found that TMCs become negligible already significantly before the quasi-PDF approaches the standard PDF. The Author suggested that, given the realistic precision of lattice simulations, TMCs can be neglected for nucleon boosts larger than around twice the nucleon mass.
4.3. Chiral Quark Models and Modeling the Relation to TMDs
Further model studies of the quasidistribution approach were performed in 2017 by W. Broniowski and E. Ruiz-Arriola [116, 117].
In the first paper [116], the pion quasi-DA and quasi-PDFs were computed in the framework of chiral quark models, namely, the Nambu-Jona-Lasinio (NJL) [118, 119] and the spectral quark model (SQM) [120–122]. The NJL model is a well-known toy model of QCD, which is a low-energy approximation to it and encompasses a mechanism of spontaneous chiral symmetry breaking from the presence of strong four-quark interactions. The SQM model, in turn, is a spectral regularization of the chiral quark model based on the introduction of the Lehmann representation of the quark propagator. The Authors derived analytical expressions for the quasi-DA and the quasi-PDF, together with their underlying unintegrated versions dependent on the transverse momentum, as well as the ITDs. They also verified the relations between different kinds of distributions found by Radyushkin [75–77]. This allowed them also to study the approach of the quasi-DA and quasi-PDF towards their light-cone counterparts and they found clear convergence for pion momenta in the range of a few GeV. Moreover, a comparison to lattice data [123] was made. For the NJL model, very good agreement was found with the lattice results at both considered pion momenta, GeV. In the case of the SQM model, similar agreement was observed at GeV and at the smaller momentum there were some differences between the model and the lattice data, but the agreement was still satisfactory. This implies that both models are able to capture the essential physical features.
In the second paper [117], Broniowski and Ruiz-Arriola explored further the relations between nucleon quasi-PDFs, PDFs and TMDs, following the work of Radyushkin [75, 76]. They derived certain sum rules, e.g., relating the moments of quasi-PDFs, PDFs, and the width of TMDs. Furthermore, Broniowski and Ruiz-Arriola modeled the factorization separating the longitudinal and transverse parton dynamics. They applied this model to study the expected change of shape of ITDs and reduced ITDs, for both quarks and gluons. They also considered the breakdown of the longitudinal-transverse factorization induced by the evolution equations, in the context of consequences for present-day lattice simulations, finding that the effects should be rather mild in quasi-PDFs, but could be visible in ITDs. Finally, they also performed comparisons to actual lattice data of the ETM Collaboration [108] for isovector unpolarized PDFs of the nucleon. The model quasi-PDFs, resulting from the assumed factorization, do not agree well with the ETMC data at 4 values of the nucleon boost between 0.94 and 2.36 GeV (see Section 3.2 for more details about these results). The discrepancy was attributed to the large pion mass, shifting the distributions to the right of the phenomenological ones, and to other lattice systematics (see also discussion about the role of the pion mass in Section 9.2.1 and [124]). More successful was the test of the aforementioned sum rule, predicting linearly increasing deviation of the second central moment of the quasi-PDF from that of the PDF with increasing , with the slope giving the TMD width. Using the ETMC data, they indeed observed the linear behavior and, moreover, extrapolating to infinite momentum, they found the second central moment to be compatible with a phenomenological analysis. In the last part of the paper, the Authors offered considerations for the pion case, presenting predictions for the valence-quark quasi-PDFs and ITDs.
4.4. Quasidistributions for Mesons in NRQCD and Two-Dimensional QCD
Meson DAs were first considered in the quasidistribution formalism in 2015 by Y. Jia and X. Xiong [125]. They calculated the one-loop corrections to quasi-DAs and light-cone DAs employing the framework of nonrelativistic QCD (NRQCD). This resulted in analytical formulae for quasi- and light-cone DAs for three -wave charmonia: the pseudoscalar and both the longitudinally and the transversely polarized vector . They checked analytically the convergence of quasi-DAs to standard DAs and performed also a numerical investigation of the rate of convergence. A function was introduced, called degree of resemblance, that quantifies the difference between the quasi- and the standard DAs. In general, momentum of around 3 times the meson mass is needed to bring the quasidistribution to within 5% of the light-cone one. The Authors also considered first inverse moments of quasi- and light-cone DAs, concluding that their rate of convergence is somewhat smaller and the difference at equal to three times the hadron mass may still be of order 20%, with 5% reached at six times larger than the meson mass.
Following the NRQCD investigation, Y. Jia and X. Xiong continued their work related to model quasidistributions of mesons. In 2018, together with S. Liang and R. Yu [126], they presented results on meson quasi-PDFs and quasi-DAs in two-dimensional QCD in its large- limit, often referred to as the ’t Hooft model [127]. The Authors used the Hamiltonian operator approach and Bars-Green equations in equal-time quantization [128], instead of the more standard diagrammatic approach in light-cone quantization. They performed a comprehensive study comparing the quasidistributions and their light-cone counterparts, studying the approach of the former to the latter at increasing meson momentum. Among the most interesting conclusions is the observation that the approach to the standard distributions is slower for lighter mesons than for heavier quarkonia of [125]. This observation was illustrated with numerical studies of the derived analytical equations for the different distributions. It was found that, for the pion, even momentum 8 times larger than the pion mass leads to significant discrepancies between the shapes of quasidistributions and light-cone ones. For () meson, in turn, momentum of five (two) times the meson mass already leads to the two types of distributions almost coinciding. An analogous phenomenon is also beginning to emerge in lattice studies and provides a warning that, e.g., pion PDFs might be more difficult to study than nucleon PDFs, i.e., require relatively larger momentum boosts.
Additionally, Jia et al. studied both types of distributions in perturbation theory, thus being able to consider the matching between quasi- and light-cone PDFs/DAs. The very important aspect of this part is that they were able to verify one of the crucial features underlying LaMET, that quasi- and light-cone distributions share the same infrared properties at leading order in . This is interesting, because the two-dimensional model has a more severe IR divergence than standard QCD.
As such, this work in two-dimensional QCD provides a benchmark for lattice studies of quasidistributions in four-dimensional QCD. It is expected that many of the obtained conclusions regarding the ’t Hooft model hold also in standard QCD. Moreover, the setup can also be used to study other proposals for obtaining the -dependence of light-cone distributions, in particular pseudo-PDFs and LCSs.
5. Theoretical Challenges of Quasi-PDFs
In this section, we summarize the main theoretical challenges related to quasi-PDFs, that have been identified early on. Addressing and understanding these challenges was very critical in order to establish sound foundations for the quasidistribution method. We concentrate on two of them, the role of the Euclidean signature (whether an equal-time correlator in Euclidean spacetime can be related to light-cone parton physics in Minkowski) and renormalizability. The latter is not trivial due to the power-law divergence inherited from the Wilson line included in the nonlocal operator. It is clear that problems related to either challenge could lead to abandoning the whole programme for quasi-PDFs. Therefore, it was absolutely crucial to prove that both of these aspects do not hide insurmountable difficulties.
5.1. Euclidean vs. Minkowski Spacetime Signature
One of the crucial assumptions of the quasidistribution approach is that these distributions computed on the lattice with Euclidean spacetime signature are the same as their Minkowski counterparts. In particular, they should share the collinear divergences, such that the UV differences can be matched using LaMET. In [129], C. Monahan and K. Orginos considered the Mellin moments of bare PDFs and bare quasi-PDFs in the context of smeared quasidistributions that differ from the standard ones only in the UV region, by construction (see Section 7 for more details about the smeared quasi-PDFs). They found that the Wick rotation from the bare Euclidean quasi-PDF to the light-cone PDF is simple.
However, in [130], a perturbative investigation was performed by C. Carlson and M. Freid, who discovered that there are qualitative differences between loop corrections in Euclidean and Minkowski spacetimes. In particular, it seemed that the IR divergence of the light-cone PDF is absent in the Euclidean quasi-PDF, which would be a problem at the basic root of LaMET. The complication emerged in certain diagrams, because the integration contours along real and imaginary axes of the complex loop temporal momentum plane could not be linked by smooth deformation, with physical observables being related to the integration along real and the lattice objects being extracted from integration along the imaginary axis. The Authors gave also a physical intuition justifying this finding. The IR divergence in Minkowski spacetime comes from collinear configurations of nearly on-shell quarks and gluons, with quark mass preventing exactly parallel configuration. Such a parallel situation is not possible in Euclidean spacetime and, hence, the Authors argued that no divergence can appear, invoking, thus, also mismatch of IR regions that could not be corrected for, perturbatively.
The serious doubts about the importance of spacetime signature were addressed in [131] by R. Briceño, M. Hansen, and C. Monahan. They formulated a general argument that for a certain class of matrix elements computed on the lattice, observables from the Euclidean-time dependence and from the LSZ reduction formula in Minkowski spacetime coincide. The class of (single-particle) matrix elements requires currents local in time, but not necessarily in space, and it includes the matrix elements needed to obtain quasi-PDFs. More precisely, the correlation functions depend on the spacetime signature, but matrix elements do not. The central issue was illustrated with a computation in a toy model, without the added, irrelevant from the point of view of the argument, complications of QCD. The focus was on a Feynman diagram directly analogous to the problematic one in [130]. The Authors, using the LSZ reduction formula, calculated its contribution to the quasi-PDF. They indeed found that the result depends on the contour of integration along the axis but pointed out that the contour along the imaginary axis does not coincide with what is done on the lattice. Instead, the connection can be made by computing the diagram contribution to a Euclidean correlator in a mixed time-momentum representation. At large Euclidean times, the result is dominated by a term which is exactly the same one as in the Minkowski calculation. After establishing the perturbative connection in the toy model for some specific kind of diagram, the proof was extended to all orders in perturbation theory. The Authors concluded their paper with a general statement about the proper prescription that yields the same result from Euclidean and Minkowski diagrams: the chosen contour must be an analytic deformation of the standard, Minkowski-signature definition of the diagram.
Thus, the apparent contradiction pointed out in [130] was fully resolved. As its Authors identified, the problem lied in the definition of the integration contour in the plane. However, the contour along the imaginary axis does not correspond to the perturbative contribution to Euclidean matrix elements, as shown in [131]. Even though the arguments of [130] turned out to be misplaced, they certainly discussed an interesting problem and they induced very valuable insights and a general proof in [131]. To our knowledge, the arguments were accepted by the Authors of [130] and no further arguments were given that would question the connection between Euclidean and Minkowski signatures in the context of quasidistributions.
5.2. Renormalizability of Quasi-PDFs
One of the indispensable components of the quasi-PDFs approach is the ability to match equal-time correlation functions (calculable on the lattice) to the light-cone PDFs using LaMET. For this approach to be successful, it is crucial that the quasi-PDFs can be factorized to normal PDFs to all orders in QCD perturbation theory, and this requires that quasi-PDFs can be multiplicatively renormalized [85]. However, the renormalization programme of quasi-PDFs is not straightforward due to the UV power divergences and, for quite some time, was not well understood (see Section 7 for recent progress).
One of the main concerns is whether the nonlocal operators are renormalizable. For example, the nonlocality of the operators does not guarantee that all divergences can be removed, due to the additional singularity structures compared to local operators and also the divergences with coefficients that are nonpolynomial. Due to the different UV behavior of quasi-PDFs and light-cone PDFs, the usual renormalization procedure is not ensured. Based on the work of [98], this originates from the different definition of the momentum fraction; that is (where () is plus momentum for the quark in the loop (initial quark)) for light-cone PDFs and (: space-like vector) in quasi-PDFs. In addition, the momentum fraction for the light-cone PDFs is restricted to , while for the quasi-PDFs can extend to . As a consequence of the above, the vertex correction (see diagrams in the second row of Figure 6) has different behavior. In fact, the UV divergences appear in the self-energy, while the vertex correction is UV finite. As pointed out in [132], the renormalizability of nonlocal operators has been proven up to two loops in perturbation theory by the analogy to the static heavy-light currents.
Thus, it is of utmost importance for the renormalizability to be confirmed to all orders in perturbation theory. This issue has been addressed independently by two groups [133–135], concluding that the Euclidean spacelike correlation functions leading to the quasi-PDFs are indeed renormalizable. These are based on two different approaches: the auxiliary heavy quark method [133, 135–137] and the diagrammatic expansion method [134, 138], employed for both quark and gluon quasi-PDFs. Below we highlight their main findings.
5.2.1. Renormalizability of Quark Quasi-PDFs
X. Ji and J.-H. Zhang in one of their early works [133] have studied renormalization of the unpolarized nonsinglet distribution. They performed an analytic calculation in the Feynman gauge to one-loop level using dimensional regularization (DR) to extract the full contributions of the diagrams entering the calculation, as shown in Figure 6. This includes the self-energy diagrams (top row) and the vertex correction diagrams (bottom row). We point out that the self-energy diagrams require integration over all components of the loop momentum, while the vertex correction diagrams have the component of the loop momentum that is parallel to the Wilson line unintegrated. This has implications on the discussion about renormalizability. This work exhibits how in this gauge all UV divergences in the vertex correction do not alter the renormalization of the quasi-PDFs, as they are removed by counterterms for subdiagrams from the interaction. Based on this, the renormalization of the quark quasidistribution reduces to the renormalization of two quark fields in the axial gauge. This study was also extended to two loop corrections, assuming one-to-one correspondence between the two-loop diagrams, as well as equivalence between the UV divergences of the two-loop self-energy in the quasiquark PDFs and of the two-loop corrections of the heavy-light quark current. Thus, multiplicative renormalizability was proven to hold up to two loops in perturbation theory. The arguments presented in this work can be generalized to include helicity and transversity PDFs.
The renormalizability of quasi-PDFs to all orders in perturbation theory has been proven for the first time by T. Ishikawa et al. in [134]. They performed the complete one-loop calculation of the quasi-PDFs in coordinate space and in the Feynman gauge, which is convenient because the renormalization of the QCD Lagrangian is known in this gauge. The one-loop calculation shows explicitly the renormalizability (to that order) of the quasi-PDFs.
More interestingly, the Authors have studied all sources of UV divergences for the nonlocal operators that enter the quasi-PDFs calculation using a primitive basis of diagrams (see Figures 3-6 in [134]). These diagrams were used to construct all possible higher-order Feynman diagrams that are presented schematically in Figure 7, and the Authors explained in great detail the proof of both power-law and logarithmic divergences being renormalized multiplicatively to all orders. This can be summarized in the calculation of the diagrams shown in Figure 7.
(a)
(b)
(c)
Diagrams of the topology shown in Figure 7(a) can be reordered in terms of one-particle-irreducible (1PI) diagrams and, therefore, one can derive all corresponding linear UV power divergences explicitly into an exponential. The latter may be removed by a mass renormalization of a test particle moving along the gauge link [139]. In addition, these diagrams have logarithmic UV divergences that can be removed by a “wave function” renormalization of the test particle [140]. The second type of diagrams (Figure 7(b)) has only logarithmic UV divergences, which can be absorbed by the coupling constant renormalization of QCD [140]. The last type of UV divergent diagrams is shown in Figure 7(c) that differs from types (a) and (b), because the loop momentum goes through an external quark, leading to divergences from higher-order loop corrections to the quark-gauge-link vertex. It was concluded that the UV divergent term of diagrams (c) is proportional to the tree level of the operator, and, therefore, a constant counterterm is sufficient to remove it. All the above constitute a concrete proof that all remaining perturbative UV divergences of the quark quasi-PDFs can be removed by introducing multiplicative renormalization factors. Exact calculations to one-loop level show that quasi-PDFs of different types do not mix under renormalization, which completes the proof of the renormalizability in coordinate space [134].
The study of the renormalizability of quark quasi-PDFs has been complemented with the work of Ji et al. in [135], in which the auxiliary heavy quark field formalism was employed. The approach shows renormalizability to all orders in perturbation theory and was confirmed in both the dimensional and the lattice regularizations, the latter for the first time. As in other studies, the focus was on the unpolarized PDFs and it was shown explicitly that the procedure mimics the renormalization of two heavy-light quark currents; the latter is valid to all orders in perturbation theory. We note that the conclusions hold for all types of PDFs and can be confirmed following the same procedure as the unpolarized one.
The introduction of a heavy quark auxiliary field, , modifies the QCD Lagrangian by including an additional term. This allows us to replace the nonlocal straight Wilson line operator by a composite operator, which is the product of two auxiliary heavy quark fieldsThus, the question of the renormalizability of the nonlocal operator can be addressed based on the renormalization of the above operator in the extended QCD theory. This has been demonstrated in DR and we urge the interested Reader to see the proof in [135]. Here we discuss the case of the lattice regulator which is particularly interesting for the numerical simulations in Lattice QCD. Unlike the case of DR, in lattice regularization (LR) the self-energy of the auxiliary quark introduces a divergence beyond leading order in perturbation theory. This may be absorbed as an effective mass counterterm, that is [141]Using the above, and for spacelike correlators, the linear divergence of Equation (31) can be factorized in the renormalized operatorwhere the remaining divergence is at most logarithmic and can be canceled to all orders in perturbation theory.
5.2.2. Renormalizability of Gluon Quasi-PDFs
For completeness, we also address the renormalizability of the gluon quasi-PDFs, which are more complicated to study compared to nonsinglet quark PDFs due to the presence of mixing. Their renormalizability was implied using arguments based on the quark quasi-PDFs [134, 135], but more recently there are direct studies for the renormalization of gluon quasi-PDFs [136–138].
The first investigation appeared in 2017 by W. Wang and S. Zhao [136], using the auxiliary field approach to study the renormalization of gluon nonlocal operators, and, in particular, the power divergences. The mixing under renormalization was also addressed. This follows their work on the matching between the quasi- and the normal gluon PDFs [142], as described in Section 8.2. The light-cone gluon PDFs are nonlocal matrix elements of the formwhere is the field strength tensor. Based on this, gluon quasidistribution can be defined by nonlocal spacelike matrix elementin which the sum over is in all directions except the direction of the Wilson line. This definition is slightly modified from the definition used in [45, 86, 142], where the sum is over the transverse directions. Despite this modification, Equation (34) is still a proper definition of a gluon quasi-PDF, as demonstrated in [136] based on the energy-momentum tensor decomposition. This is also confirmed numerically, as the one-loop matching to the light-cone PDFs coincides for the two definitions.
Reference [136] presented the complete one-loop calculation for the gluon operator of Equation (34), introducing a UV cutoff on the transverse momentum. The calculation was performed in Feynman gauge and in the adjoint representation. The relevant one-loop diagrams can be separated into two categories: (1) diagrams in which the vertex from the operator does not include gluons from the Wilson line (shown in Figure 8) and diagrams that have at least one gluon from the Wilson line in the vertex of the operator, as shown in Figure 9. This calculation identified all divergences including the linear divergence, and, unlike the case of the quark quasi-PDFs, the Wilson line self-energy (right diagram of Figure 9) is not the only source of linear divergence in the gluon distributions. As a consequence, it is not possible to absorb all linear divergences in the renormalization of the Wilson line, but a more complicated renormalization is needed. However, as argued in [137], this is due to the choice of a nongauge invariant regulator.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
One approach to study the renormalization of the quasigluon PDFs is to introduce an auxiliary heavy quark field, as adopted in the renormalization of the quark distributions. This auxiliary field is in the adjoint representation of and does not have spin degrees of freedom. Therefore, this approach allows one to study local operators instead of the nonlocal operator of Equation (34). Mixing between these new operators and the gauge invariant gluon field strength tensor is permitted. In addition, it was shown at the level of one-loop corrections that the power divergence can be absorbed in the matrix elements of the local operators, which is expected to hold to all orders in perturbation theory. This work by W. Wang and S. Zhao has contributed to understanding the renormalization of the gluon quasi-PDFs, but there were a number of issues to be addressed. This included, but was not limited to, (a) the study of gauge fields in the fundamental representation and the corresponding mixing and (b) the study of the renormalization in the lattice regularization, preferably nonperturbatively. The latter is highly nontrivial and technically more complicated than in the case of quarks.
The renormalizability of both the unpolarized and the helicity gluon PDFs has been studied by J.-H. Zhang et al. in [137], including possible mixing that is permitted by the symmetries of the theory. The auxiliary field formalism was employed in a similar fashion as the studies presented above [133, 135–137]. Explicit results were given for the unpolarized quasigluon PDFs in the dimensional and gauge-invariant cutoff regularizations.
In the auxiliary field formalism, the operator presented in Equation (34) ( summed over the transverse directions) can be replaced by a new operator; that is, where , . denotes the auxiliary adjoint “heavy quark” field. For a proof, see [137]. Based on symmetry properties, such a composite operator can mix with lower-dimensional operators that are gauge-invariant, BRST variations or vanish by the equations of motion. The identified mixing pattern helps to construct the proper operators for the gluon quasi-PDFs that are multiplicatively renormalizable. In particular, three (four) operators are identified for the unpolarized (helicity) gluon PDFs, that do not suffer from mixing. Here we provide the operators for the unpolarized case: where the index represents the temporal direction and the direction of the Wilson line. In addition, runs over all Lorentz components, while over the transverse components only (). In a similar way, it was found that three operators related to the gluon helicity distributions can be renormalized multiplicatively. For details, see Section III C of [137]. This work provides crucial guidance for numerical simulations in Lattice QCD and the development of a nonperturbative renormalization prescription. Based on the mixing pattern, the Authors provided a renormalization prescription suitable for lattice simulations, and a factorization for gluon and quark quasi-PDFs.
Z.-Y. Li, Y.-Q. Ma, and J.-W. Qiu have studied renormalizability of gluon quasi-PDFs in [138], a work that appeared simultaneously with [137]. Their work is based on diagrammatic expansion approach, as studied for the quark quasi-PDFs [134, 138]. By studying the UV divergence of gluon operators, it was demonstrated that appropriate combinations can be constructed, so that their renormalization is multiplicative to all orders in perturbation theory. Such operators are related to gluon quasi-PDFs. The demonstration is based on a quasigluon operator that has a general form where is the Wilson line with gauge links in the adjoint representation.
The procedure followed in this work is based on a one-loop calculation of the Green’s functions which is performed in DR. It was demonstrated that the linear UV divergences of the gluon-gauge-link vertex are canceled explicitly. This was extended to all loops in perturbation theory by investigating all possible UV divergent topologies of higher-order diagrams, showing that the corresponding linear UV divergences are canceled to all orders in perturbation theory. It was also discussed in detail that the UV divergences of all 36 pure quasigluon operators (including the antisymmetry of gluon field strength) can be multiplicatively renormalized. This work, thus, constitutes a powerful proof of the renormalizability of gluon quasi-PDFs.
6. Lattice Techniques and Challenges for Quasi-PDFs
Apart from theoretical challenges of the quasidistribution approach, discussed in the previous section, also the lattice implementation and efficiency of computations are a major issue for the feasibility of the whole programme. In this section, we discuss these aspects in some detail, showing that tremendous progress has been achieved also on this side. In addition, we discuss challenges for the lattice that need to be overcome for a fully reliable extraction of PDFs.
6.1. Lattice Computation of Matrix Elements
To access quasi-PDFs of the quarks in the nucleon, one needs to compute the following matrix elements:where the Dirac structure determines the type of quasi-PDF (see below), is the boosted nucleon state with momentum , and is a Wilson line of length along the spatial direction of the boost. To obtain the above matrix elements, one constructs a ratio of three-point and two-point functions:where is a kinematic factor that depends on the Dirac structure, and the correlation functions are computed according towith the proton interpolating operator, , the current insertion time, parity plus projector for the two-point function, , and parity projector for the three-point functions, , dependent on the Dirac structure of the current.
The Wick contractions for the three-point function lead, in general, to a quark-connected and a quark-disconnected diagram. Since the evaluation of the latter is far more demanding than that of the former, the numerical efforts were so far restricted to connected diagrams only. One uses the fact that disconnected diagrams cancel when considering the flavor nonsinglet combination in the formulation of Lattice QCD with degenerate light quarks. The connected diagram that contributes to the three-point function is shown in Figure 10.
Special attention has to be paid to the Dirac structure of the insertion operator, because mixing appears among certain structures, as discovered in [106]. In particular, the originally suggested for the unpolarized PDF mixes with the scalar operator; see Section 7 for details. Such mixing can be taken into account by explicitly computing a mixing matrix of renormalization functions and matrix elements for both Dirac structures. However, in practice, this leads to much worse signal and finally to a much less precise estimate of the PDFs. For this reason, the strongly preferred choice is for the unpolarized quasi-PDF. Similar mixing occurs for the polarized cases for certain Dirac structures, with the choice of and or for helicity and transversity, respectively, guaranteeing that no mixing is present [106].
We now turn to describing the lattice computation in more detail. For the two-point function, Wick contractions lead to standard point-to-all propagators that can be obtained from inversions of the Dirac operator matrix on a point source. The computation of the three-point function is more complicated. Apart from the point-to-all propagator, it requires the knowledge of the all-to-all propagator. Two main techniques exist to evaluate this object: the sequential method [143] and the stochastic method [144]. In the former, one constructs a so-called sequential source from a suitable point-to-all propagator. Inverting the Dirac matrix on this source, the sequential propagator is obtained that enters in the three-point function. The other method employs stochastic noise sources on a single time slice, leading to a stochastic estimate of the all-to-all propagator upon inversion of the Dirac matrix. In principle, the second method is more flexible, as it allows for obtaining results for all Dirac structures and all momenta with the same inversions, the most costly part of the computation. The price to pay is the introduction of stochastic noise, but the overhead introduced by the necessity to suppress this noise is still more than compensated by the gain from flexibility, in principle. Using the sequential method or, more precisely, its fixed sink variant, implies that the momentum at the sink has to be fixed and separate inversions are needed for each nucleon boost, as well as for each Dirac structure due to different projectors.
In the early studies, both approaches were tested by ETMC [47], with the conclusion that they yield compatible results and the additional noise from the stochastic method can be suppressed by using 3-5 stochastic noise vectors. Given the flexibility of the stochastic method, ETMC decided to pursue studies with this approach in [100, 108]. In [108], the technique was changed to one involving the sequential propagator for reasons explained in the next subsection. The method for computing the all-to-all propagator was not revealed by the Authors of the other exploratory numerical study of quasi-PDFs in [46, 107].
Having computed the three-point and two-point functions, the relevant matrix elements can be obtained. The crucial issue that has to be paid special attention to is the contamination of the desired ground state matrix elements by excited states. Three major techniques are available: single-state (plateau), multistate, and summation fits. We briefly describe all of them below.(i) Plateau method. The most straightforward way of obtaining the matrix element from the three-point and two-point functions is to identify a region where their ratio is independent of the insertion time and fitting to a constant, which is the matrix element of the ground state. As can be seen from the spectral decomposition of the three-point function, excited states manifest themselves as curvature in the ratio of Equation (88) and also in the shift of its central value. Under realistic statistical uncertainties, it is, therefore, not always clear whether an actual plateau has been reached and, thus, it is not advisable to use this method as the sole method of extracting the ground state properties.(ii) Summation method. This approach [145, 146] consists in summing the ratios of three-point and two-point functions over the insertion time . By decomposing the correlators into sums of exponential terms, one obtains a geometric series, leading finally to where the source and sink timeslices are excluded avoiding contact terms and is a constant. The ground state matrix element, , is then extracted from a linear two-parameter fit to data at sufficiently large source-sink separations . The method has the advantage that excited states are suppressed by a faster-decaying exponential with respect to the plateau fits, but the statistical uncertainties are, typically, much larger.(iii) Multistate fits. A natural generalization of the plateau method is to include higher-order exponential terms in the decomposition of the two-point and three-point functions, typically the first excited state (two-state fits) or the lowest two excited states (three-state fits). In general, the two-point correlator can be written as with amplitudes and energies of subsequent states . The three-point function reads with matrix elements of the suitable operator in addition to parameters in the two-point correlator. Note that, in practice, it is difficult to consistently go beyond one or two excited states, as the number of fitting parameters is increasing faster than linearly with increased number of excited states taken into account, due to the presence of the mixed matrix elements, , for a growing number of pairs .
In principle, the multistate method (realistically two-state method) allows for a better control of excited states contamination. However, in realistic lattice situations, the interpolating operators used to create the nucleon from the vacuum excite numerous states with the same quantum numbers. This contamination increases with pion mass decreasing towards the physical point; see, e.g., [39], for an illustration. Moreover, the number of excited states increases with larger nucleon boosts. All the above imply that it is unlikely to achieve a regime of source-sink separations where precisely two states play a role. Thus, also relying solely on two-state fits should not be used as the only method. Instead, ground state dominance should be established by aiming at compatibility between all three methods of extracting the ground state matrix elements. Such compatibility ensures that the probed regime of values has enough suppression of excited states and excludes that many excited states mimic a single excited state. Numerically, it is hard to disentangle several excited states and the manifestation of many of them appearing would be clear incompatibility of two-state fits with plateau fits. Note also that, with exponentially decaying signal-to-noise ratio at larger source-sink separations, the danger of the two-state approach is that the fits may easily be dominated by data at the lower values, heavily contaminated by excited states. Thus, we are led to conclude that the most reliable estimates of ground state matrix elements ensue from compatible results obtained using all of the three above methods (with the summation method being, in most cases, inconclusive due to large statistical uncertainty). It is still important to bear in mind that excited states are never fully eliminated, but only exponentially suppressed. This means that the ground state dominance is always established only to some level of precision. Aiming at increased statistical precision, the previously reliable source-sink separation(s) may prove to be insufficient. Obviously, when targeting larger momenta and/or smaller pion masses, conclusions for the role of excited states at a smaller boost or a larger pion mass do not apply; hence, a careful analysis is always needed at least in the setup most prone to excited states.
Having extracted the relevant matrix elements, one is finally ready to calculate the quasi-PDF. We rewrite here the definition of quasi-PDFs with a discretized form of the Fourier transform:where the factor of ensures correct normalization for different momenta and is to be chosen such that the matrix elements have decayed to zero both in the real and in the imaginary part (see also Section 6.3).
6.2. Optimization of the Lattice Computation
In the previous subsection, we have established the framework for the computation of quasi-PDF matrix elements on the lattice. Now, we describe some more techniques that are usually used to perform the calculation as effectively as possible.
The first technique, commonly employed in lattice hadron structure computations, serves the purpose of optimizing the overlap of the interpolating operator that creates and annihilates the nucleon with the ground state. This can be achieved by employing Gaussian smearing [147, 148] of fermionic fields, which reflects the fact that hadrons are not point-like, but are extended objects. Moreover, the smearing is further optimized by combining it with a technique for reducing short-range (UV) fluctuations of the gauge fields; gauge links used in the quark fields smearing are subjected to APE smearing [149]. The procedure involves optimizing four parameters, the parameters regulating the “strength” of the Gaussian and APE smearing, and , respectively, and the number of Gaussian and APE smearing iterations. The typical criterion of optimization is that the root mean square (rms) radius of the proton should be around 0.5 fm.
Smearing techniques are used also to decrease UV fluctuations in gauge links entering the Wilson line in the operator insertion. In principle, any kind of smearing can be used for this purpose, with practical choices employed so far of HYP smearing [101] and stout smearing [150]. The smearing of the Wilson line has the additional effect of reducing the UV power divergence related to the Wilson line, i.e., shifting the values of renormalization factors towards their tree-level values, and thus suppressing the power-like divergence. Thus, a relatively large number of smearing iterations was used in the early works, which was necessary due to the absence of the renormalization. In principle, the renormalized matrix elements should not depend on the amount of smearing applied to the operator and it is an important consistency check to confirm this. We note that, before the advent of the full nonperturbative renormalization programme for quasi-PDFs [151], the role played by the Wilson line links smearing was somewhat different. Without explicit renormalization, the results were contaminated by the power divergence and the smearing had the task of subtracting a possibly large part of this divergence in the hope of obtaining preliminary results at not too small lattice spacings. After this premise lost its significance, this kind of smearing is applied only to reduce gauge noise to a certain extent. Alternatively, smearing of gauge links can also be applied to the whole gauge field, i.e., enter both the Wilson line and also the Dirac operator matrix. Note, however, that in this way it is not possible to check explicitly that the renormalized results are independent of the smearing level, at least without costly additional Dirac matrix inversions for different numbers of smearing iterations.
All the above techniques are rather standard and have been employed in the quasi-PDFs computations already in the very first exploratory studies. However, the recent progress that we review in Section 9 would not have been possible without the technique of so-called momentum smearing [109, 152]. It is a relatively simple extension of the quark fields smearing described above. The crucial observation is that a Gaussian-smeared nucleon state has maximal overlap with a nucleon at rest; i.e., it is centered around zero momentum in momentum space. It is, hence, enough to move this center to the desired momentum to obtain an improved signal for a boosted nucleon. The modification is the addition of a phase factor in the position space definition of the smearing, where is the desired nucleon momentum and a tunable parameter8. Explicitly, the modified Gaussian momentum smearing function readswhere are gauge links in the -direction. For optimal results, the parameter should be tuned separately for every momentum and every ensemble. In the context of quasi-PDFs, momentum smearing has first been applied in the ETMC study reported in Section 3.2 [108, 153]. By now, it has become a standard technique for enhancing the signal. We note, however, that momentum smearing does not fully solve the exponentially hard problem of decaying signal at large boosts, but rather it moves it towards larger momenta. Therefore, accessing highly boosted nucleon on the lattice, necessary for reliable matching to light-cone PDFs via LaMET, remains a challenge.
To finalize this subsection, we mention one more useful technique that is applied nowadays to decrease statistical uncertainties at fixed computing time. The most expensive part of the calculation of the correlation functions is the computation of the quark propagators, i.e., the inversion of the Dirac operator matrix on specified sources. This is typically done using specialized iterative algorithms, often tailored to the used fermion discretization. The iterative algorithm is run until the residual, quantifying the distance of the current solution with respect to the true solution, falls below some tolerance level, . The standard way is to set to a very small number, of order . However, obviously that may need iterating the solver for a long time. To save some considerable fraction of computing time, truncated solver methods have been invented, where the precision is relaxed to . Naturally, relaxed precision of the solver leads, in general, to a bias introduced in the solution. Hence, the second ingredient of these methods is bias correction. Below, we shortly describe one of such methods, the Covariant Approximation Averaging (CAA) [154]. One performs a certain number of low-precision (LP) inversions, , accompanied by a smaller number of standard, high-precision (HP) inversions, . The final correlation functions are defined as follows:where and denote correlation functions obtained from LP and HP inversions, respectively. To correct the bias properly, HP and LP inversions have to be done for the same source positions. The choice of the numbers of LP and HP inversions has to be tuned in such a way to maintain a large correlation coefficient (typically 0.99-0.999) between LP and HP correlators, which guarantees that the bias is properly subtracted.
6.3. Lattice Challenges
In this section, we discuss the challenges for lattice computations of quasi-PDFs. On the one side, this includes “standard” lattice challenges, like control over different kinds of systematic effects, some of them enhanced by the specifics of the involved observables. On the other side, the calculation of quasi-PDFs offered new challenges that had to or have to be overcome for the final reliable extraction of light-cone distributions. Below, we discuss these issues in considerable detail, starting with the “standard” ones and going towards more specific ones.(1) Discretization effects. Lattice simulations are, necessarily, performed at finite lattice spacings. Nevertheless, the goal is to extract properties or observables of continuum QCD. At finite lattice spacing, these are contaminated by discretization (cutoff) effects, which need to be subtracted in a suitable continuum limit extrapolation. Obviously, prior to taking the continuum limit, the observables need to be renormalized and we discuss this issue in Section 7. Assuming divergences have been removed in a chosen renormalization scheme, the continuum limit can be taken by simulating at three or more lattice spacings and fitting the data to an appropriate ansatz, typically linear in the leading discretization effects, of order or . In most Lattice QCD applications, -improved fermionic discretizations or observables are used. In many cases this, however, requires calculation of observable-specific improvement coefficients (e.g., for Wilson-clover fermions). It remains to be shown how to obtain improvement of quasi-PDFs at least for some of the fermionic discretizations. Up to date, quasi-PDFs studies have been performed for a single lattice spacing in a given setup and, hence, discretization effects have not been reliably estimated. Going to smaller lattice spacings remains a challenge for the future. It is not a problem in principle, but obviously it requires huge computational resources, especially at the physical pion mass. However, there are indirect premises that discretization effects are not large. Firstly, they have been relatively small in general lattice hadron structure calculations. Secondly, indirect evidence for the smallness of cutoff effects is provided by checks of the dispersion relation, i.e., the relation between energy of a boosted nucleon and its momentum (see Section 9). In the absence of large discretization effects, the continuum relativistic dispersion relation holds. Note, however, that cutoff effects can be enhanced if the nucleon boost becomes larger than the lattice UV cutoff; i.e., if . In principle, no energy scale on the lattice should exceed . Precisely for this reason, lattice calculations involving the heavy quark need its special treatment; with typical lattice spacings, the bottom quark mass exceeds and a reliable computation must involve an effective theory treatment, such as the one provided by HQET or by NRQCD.(2) Finite volume effects. Apart from finite lattice spacing, also the volume of a numerical simulation is necessarily finite. Thus, another lattice systematic uncertainty may stem from finite volume effects (FVE). FVE become important if the hadron size becomes significant in comparison with the box size. The hadron size is to a large extent dictated by the inverse mass of the lightest particle in the theory. Hence, leading-order FVE are related to the pion mass of the simulation and smaller pion masses require larger lattice sizes in physical units to suppress FVE. Usually, FVE are exponentially suppressed as , where is the spatial extent of the lattice. The typical rule adopted in lattice simulations is that this suppression is enough if . At nonphysical pion masses of order 300-400 MeV, this corresponds to a box size of 2-2.5 fm, which is easy to reach with typically used lattice spacings, 0.05-0.1 fm. When simulating at the physical pion mass, the minimal box size that yields is 6 fm and, thus, finer lattice spacings require huge lattices. Nevertheless, lattice hadron structure calculations have usually evinced rather small FVE already with . Still, an explicit check of FVE is highly advisable when aiming at a fully reliable computation. Above, the main source of FVE that we considered was related to the size of hadrons. However, it was pointed out in [155] that, for quasi-PDFs, a relevant source of FVE may be the size of the Wilson line in the operator inserted in the matrix elements defining quasidistributions. The Authors studied perturbatively a toy scalar model with a light degree of freedom (mimicking the pion in QCD) and a heavy one (corresponding to the nucleon). The studied matrix element involved a product of two currents displaced by a vector of length and they found two kinds of FVE: one decaying with and the other one with , where is the mass of the heavy state. Moreover, both exponentials have prefactors scaling as (with some exponents and ), that can further enhance FVE for larger displacements . In the case of pion matrix elements, the FVE may be particularly enhanced by . Even though the studied case concerned a product of two currents, not quark fields connected by a Wilson line, some enhancement of FVE may also occur for the latter case. In view of this, investigation of FVE in matrix elements for quasi-PDFs, especially ones with larger lengths of the Wilson line, is well motivated. It is also important to mention that finite lattice extent in the direction of the boost, , imposes a limit on the minimal Bjorken- that can be reached. The parton momentum is , which determines its correlation length to be of order . This value should be smaller than the physical size of the boost direction, . At the same time, the boost should be smaller than the lattice UV cutoff; i.e., . Replacing “” symbols in the above inequalities with “” signs, one arrives at the minimal accessible on the lattice: ; i.e., . Note that it is the number of sites in the boost direction that determines , not its physical size.(3) Pion mass dependence. The computational cost of Lattice QCD calculations depends on the pion mass. Hence, exploratory studies are usually performed with heavier-than-physical pions, as was also the case for quasi-PDFs (see Section 3.2). Obviously, this introduces a systematic effect. If no physical pion mass calculations are available, one can extrapolate to the physical point, if the fitting ansatz for this extrapolation is known (e.g., from chiral perturbation theory). However, the cleanest procedure is to simulate directly with pions of physical mass. Recently, quasi-PDFs computations with physical pions have become available; see Section 9 for their review including a direct comparison between ensembles with different pion mass [124].(4) Number of flavors, isospin breaking. QCD encompasses six flavors of quarks. However, due to the orders of magnitude difference between their masses, only the lightest two, three, or four flavors are included in lattice simulations. Moreover, the up and down quarks are often taken to be degenerate; i.e, one assumes exact isospin symmetry. One then speaks of a , , or setup, respectively. Differences among these setups are observable dependent, but usually smaller than other systematic uncertainties and the statistical errors. For examples of the small dependence on the number of dynamical quarks in various observables, see, e.g., the FLAG review [156]. Hence, for most applications, all these setups can be considered to be equivalently suitable. Only when aiming at total uncertainty, well beyond the current precision of the field of lattice PDFs, it may be necessary to include dynamical strange and charm quarks. Similar or smaller effects are expected from isospin breaking by the different up and down quark masses (QCD effect) and their different electric charges (QED effect). The order of magnitude of these effects can be deduced from the difference of proton and neutron masses, less than two per mille. Note that the setup with degenerate light quarks is very useful in lattice hadron structure calculations also for practical reasons; in such a setup, the disconnected contributions cancel in the flavor combination and, moreover, isovector PDFs do not mix under matching and renormalization. Thus, it is clear that all the effects discussed in this point are currently subleading but may become important in the future, when aiming at very precise extractions of PDFs.(5) Source-sink separation and excited states contamination. As already discussed in Section 6.1, a significant systematic effect may emerge in lattice matrix elements due to excited states contamination. From the correlation functions decomposition, one can see that excited states are suppressed with the source-sink separation, . Hence, a careful analysis of a few separations is needed to establish ground state dominance; see Section 6.1 for more details. The issue of reaching large values is nontrivial from the computational point of view, as the signal-to-noise ratio decays exponentially with increasing . For this reason, a compromise is needed to keep the computational cost under control. Yet, the compromise must not affect the reliability of the results.(6) Momentum boost and higher-twist effects. Contact with the IMF via LaMET is established at large nucleon momenta. Hence, it is desirable to use large nucleon boosts on the lattice. However, this is highly nontrivial for several reasons. First, the signal-to-noise ratio decays exponentially with increasing hadron momentum, necessitating increase of statistics to keep similar statistical precision at larger boosts. Second, excited states contamination increases considerably at larger momenta, calling for an increase of the source-sink separation to maintain suppression of excited states at the same level. As argued in the previous point, the increase of further decays the signal, enlarging the required statistics. Third, large hadron momenta may induce enhanced discretization effects, in particular when the boost becomes similar to or larger than the lattice UV cutoff, i.e., the inverse lattice spacing. Thus, momenta larger than the UV cutoffs of the currently employed lattice spacings, of order 2-2.5 GeV, should only be aimed at with ensembles at finer lattice spacings. We now consider effects that may appear if the nucleon momentum is too small. Looking at the formulation of LaMET, it is clear that higher-twist effects (HTE), suppressed as , may become sizable and hinder the extraction of leading-twist PDFs. In principle, one can compute the HTE explicitly and subtract them. This would be an interesting direction of further studies, especially that HTE are of interest in their own right. Alternatively, one may compute the leading functional dependence of HTE and extrapolate them away. An example of such computation was presented in [157], based on the study of renormalons in coefficient functions within the bubble-chain approximation. The result for quasi-PDFs is an correction. Note, however, that the matrix elements underlying the quasi-PDFs in this analysis are normalized to unity at zero momentum, as done in the pseudo-PDF approach (see Section 2.6). This suppresses HTE at small- at the price of enhancement for large-. Clearly, the renormalization programme employed for quasi-PDFs, e.g., based on a variant of RI/MOM (see Section 7), can lead to different functional form of HTE. Moreover, the Authors of [157] put another warning that a perturbative analysis might not see all sources of HTE and their results should rather be considered as a minimal model that may miss nonperturbative features. Note also that knowing the functional form of leading-order HTE (with unknown prefactors) does not clarify what the range of hadron momenta is where these terms are indeed leading. At too small momenta, it may still be that higher-order HTE are sizable and even change the overall sign of the correction, rendering the extrapolation unreliable. Another type of HTE is nucleon mass corrections (NMCs). These, in turn, can be exactly corrected by using the formulae derived by Chen et al. [107]. The calculation presented in this reference allowed us to obtain closed expressions for the mass corrections relevant for all types of quasi-PDFs. An important feature of these NMCs is that the particle number is conserved. We note that NMCs are already small at momenta not much larger than the nucleon mass, as also argued by Radyushkin [76] from a model calculation. It is important to remark that the NMCs for quasi-PDFs (also commonly referred to as TMCs) are different from TMCs in phenomenological analyses for standard PDFs (see, e.g., [158] for a review). NMCs in quasi-PDFs result from the nonzero ratio of the nucleon mass to its momentum (while this ratio is zero in the IMF), whereas TMCs in phenomenological analyses refer to corrections needed because of a nonzero mass of the target in a scattering experiment. At the level of matrix elements, the momentum dependence is manifested, inter alia, by the physical distance at which they decay to zero. This distance, entering in the limits of summation for the discretized Fourier transform in Equation (49), becomes smaller for larger values of . If it is too large, periodicity of the Fourier transform will induce nonphysical oscillations in the quasi-PDFs, especially at large . We note that these oscillations do not appear because of the truncation at finite , but rather because of a too large value of at low momenta. This effect can be naturally suppressed by simulating at larger nucleon boosts and indeed, as we show in Section 9, oscillations are dampened at larger . The uncertainty induced by this behavior can also result from uncertainties related to the renormalization of bare matrix elements. The large values of -factors amplify both the real and the imaginary part and, for complex -factors, also mix them with each other. The -factors should be purely real, but this feature holds only if conversion between the intermediate lattice renormalization scheme and the scheme is done to all orders in perturbation theory. Together with lattice artifacts appearing in the estimate of the intermediate scheme renormalization functions, this effectively induces a slower decay of matrix elements with the Wilson line length and shifts the value of where matrix elements become zero to larger distances. Hence, the combination of too small boost and uncertainties in -factors manifests itself in the oscillations. Note also that the problem may be more fundamental. It is not presently clear how reliable is a distribution reconstruction procedure from a set of necessarily limited data. This issue is being investigated in the context of pseudodistributions [159]. The reconstruction techniques, mentioned in the context of the hadronic tensor in Section 2.1, may be crucial for control of this aspect. It was also speculated [51] that the Fourier transform may be a fundamental limitation of the quasi- and pseudodistribution approaches and the fundamental object may be the renormalized matrix element, or ITD. A method to remove the nonphysical oscillations was proposed in [160] and was termed the derivative method. One rewrites the Fourier transform using integration by parts: where the derivative of the matrix elements with respect to the Wilson line length gives the name to the method. The integration by parts is exact and this definition of the Fourier transform is equivalent to the standard one if the matrix elements have decayed to zero at and up to discretization effects induced by the need to lattice size the continuous derivative. Otherwise, one neglects the surface term in Equation (52), which effectively absorbs oscillations. However, it is debatable whether the procedure is safe and the neglected surface term does not hide also physical contributions at a given nucleon boost. Also, the presence of an explicit factor in the surface term leads to an uncontrolled approximation for small values of . Other proposed methods to remove the oscillations are a low-pass filter [160], including a Gaussian weight in the Fourier transform [161]. However, they have not been used with real lattice data. Ideally, the nucleon momentum needs to be large enough to remove oscillations in a natural way, instead of attempting to suppress them artificially.(7) Other effects in the PDFs extraction procedure. For the sake of completeness, we mention other effects that can undermine the precision of lattice extraction of PDFs, although they are not challenges for the lattice per se. In the previous point, we have already mentioned uncertainties related to renormalization. In RI/MOM-type schemes, they manifest themselves in the dependence of -factors on RI scales from which they were extracted, even after evolution to a common scale. This can be traced back to the breaking of continuum rotational invariance () to a hypercubic subgroup . A way to overcome this problem is to subtract lattice artifacts computed in lattice perturbation theory, which can be done to all orders in the lattice spacing at the one-loop level; see [162] for more details about this method and an application to local -factors. Another renormalization-related issue is the perturbative conversion from the intermediate lattice scheme to the scheme and evolution to a reference scale. Although not mandatory, the aim of the whole programme is to provide PDFs in the scheme of choice for phenomenological applications, i.e., the scheme. The conversion and evolution are currently performed using one-loop formulae and, hence, subject to perturbative truncation effects. A two-loop calculation of these steps will shed light on the magnitude of truncation effects. Similarly, truncation effects emerge also in the matching of quasi-PDFs to light-cone PDFs, currently done to one-loop level; see Section 8 for a more thorough discussion on matching.(8) Finite and power-divergent mixings. A general feature of quantum field theory is that operator mixing under renormalization is bound to appear among operators that share the same symmetry properties. On the lattice, some continuum symmetries, that otherwise prevent mixing, are broken. For operators of the same dimension, the mixing is finite. Important example of such mixing was mentioned above; for some fermionic discretizations, operator with the Dirac structure (for unpolarized PDF) has the same symmetries as the analogous scalar operator [106] and hence mixes with it, while the structure has different symmetry properties and avoids the mixing. This mixing is a lattice effect stemming from chiral symmetry breaking by the lattice discretization and does not appear for lattice fermion formulations that preserve this symmetry, e.g., overlap fermions. We discuss this finite mixing in more detail in the next section. If the dimension of the operator with the same symmetries is lower, then the mixing will be power divergent in the lattice spacing; i.e., it will contribute a term , where is the difference in the dimension. The possibility that such mixings occur for quasi-PDFs, as well as pseudo-PDFs and LCSs, was considered by G.C. Rossi and M. Testa in [163, 164]. They considered a toy model, devoid of QCD complications irrelevant in the context of their argument, and showed that moments of quasi-PDFs evince power-divergent mixings with lower-dimensional operators coming from the trace terms. They argued that, for a proper lattice extraction, such mixings would have to be computed and subtracted. However, it was counterargued in three papers [84, 90, 165] that the problem actually does not exist. In [165], it was pointed out that indeed all moments, except for the zeroth, do not converge. However, the light-cone PDFs are extracted from the nonlocal quasidistributions that avoid the power divergence problem; i.e, moments of quasi-PDFs are never intended to be computed. They trace it back to the much simpler ultraviolet physics in the nonlocal formulation, where, apart from the Wilson-line-induced power divergence, shown to be renormalizable (see Sections 5 and 7), there are only logarithmic divergences. All of these divergences can be properly renormalized on the lattice, e.g., in a RI/MOM-type scheme. It was also argued by Rossi and Testa that divergent moments of quasi-PDFs, , necessarily imply divergent moments of extracted light-cone PDFs, , since the latter are proportional to the former. However, this argument ignores the presence of moments of the matching function, : It is exactly the matching function that makes the moments of standard PDFs finite after the subtraction of the UV differences between the two types of distributions. In other words, the divergence in the moments is exactly canceled by the divergence of moments , yielding finite moments of light-cone PDFs. Further explanations were provided in [84]. Radyushkin pointed out that Rossi and Testa rely on a Taylor expansion in . This expansion may be justified in the very soft case when all derivatives with respect to exist at . However, in the general case, the use of the Taylor expansion for the hard logarithm “amounts to just asking for trouble.” Crucially, it is the part that contributes slowly decreasing terms into the large- part of quasi-PDFs and these terms lead to the divergence of the quasidistribution moments. These terms are not eliminated by just taking the infinite momentum limit, but they disappear upon the matching procedure. As a result, one can calculate the moments of light-cone PDFs from the quasi-PDF data. Finally, J. Karpie, K. Orginos, and S. Zafeiropoulos demonstrated [90] explicitly that the problem does not appear for pseudo-PDFs, refuting claim thereof by Rossi and Testa. The reduced ITDs were OPE expanded in lattice-regularized twist-2 operators, which indeed have power divergences on the lattice, due to the breaking of the rotational symmetry. However, the Wilson coefficients in the OPE have exactly the same power divergences and cancel the power divergences of the matrix elements, order by order in the expansion, and the final series is finite to all orders. The Authors also provided an explicit numerical demonstration for the first two moments, obtaining compatibility within errors with an earlier lattice calculation in the same quenched setup; see Section 11.4 for more details. With all these developments, it has been convincingly established that the problem advocated by Rossi and Testa does not hinder the lattice extraction of light-cone PDFs. Thus, power-divergent mixings only manifest themselves in certain quantities, like moments of quasi-PDFs, which are nonphysical. In turn, finite mixings can be avoided by a proper Dirac structure of matrix elements.
We finalize this section with a schematic flowchart (Figure 11), prepared by C. Monahan, representing the various steps on the way from bare matrix elements to final light-cone PDFs. Some of the discussed above challenges for the lattice computations are indicated. We refer also to [51] for another discussion of systematic effects.
7. Renormalization of Nonlocal Operators
The renormalization of nonlocal operators that include a Wilson line is a main component of the lattice calculation related to quasi-PDFs. Lattice results from the numerical simulations can only be related to physical quantities upon appropriate renormalization and only then comparison with experimental and phenomenological estimates becomes a real possibility. As discussed in Section 5.2, the renormalizability of the straight Wilson line bilinear operators has been investigated early on by Ji and Zhang [133] to one loop in perturbation theory, concluding that such operators are multiplicatively renormalizable. The argument was also extended to two-loop level. Ishikawa et al. showed in [132] the feasibility of the subtraction of the power divergence present in the operators under study to achieve a well-defined matching between the quasi-PDFs with the light-cone PDFs. These studies were later expanded to prove renormalizability of the operators to all orders in perturbation theory [134, 135], including the lattice regularization.
Since the proposal of Ji in 2013, several aspects of quasi-PDFs have been investigated, such as the feasibility of a calculation from Lattice QCD. This includes algorithmic developments [109, 152, 154] that lead to simulations at the physical point, and the matching between the quasi- and light-cone PDFs [48, 83, 124, 142, 166]. Thus, the lattice calculations have progressed with a missing ingredient: its renormalization. It is not until 2017 that a proper renormalization prescription has been proposed, despite the theoretical developments on the renormalizability of the nonlocal operators of interest. This has proven to be a challenging and delicate process due to the presence of the Wilson line that brings in additional power divergences, the nonlocality and the complex nature of the matrix elements. As a consequence, the first studies of quasi-PDFs on the lattice were either neglecting renormalization [46] or multiplying naively the matrix elements with the renormalization function of the corresponding local operators [100, 107, 108], a procedure understood as normalization.
7.1. Power Divergence
Among the first attempts to understand the renormalization of nonlocal operators was to address the power divergence inherited from the Wilson line in the static potential approach, as described in this subsection. Eliminating the power divergence results in a well-defined matching between the quasi-PDFs and the light-cone PDFs.
The renormalization of nonlocal operators in gauge theories has been investigated long time ago [167–170], and later in the 1980’s and 1990’s [140, 141, 171–176]. In these seminal works, it was identified that Wilson loops along a smooth contour with length , computed in dimensional regularization (DR), are finite functions of the renormalized coupling constant, while other regularization schemes may lead to additional renormalization functions, that isIn the above expression, the subscript indicates the distance between the end points of the Wilson line, whereas is mass renormalization of a test particle that moves along . Also, the logarithmic divergences can be factorized within , and the power divergence is included in . In particular, in the lattice regularization (LR), the latter divergence manifests itself in terms of a power divergence with respect to the UV cutoff, the inverse of the lattice spacing ,where is dimensionless. This is analogous to the heavy quark approach, where similar characteristics are observed. For instance, a straight Wilson line may represent a static heavy quark propagator, and a corresponding mass shift. Inspired by this analogy, Chen et al. [177] and Ishikawa et al. [132] proposed a modification such that spacelike correlators do not suffer from any power divergence. In their work, the matrix element appearing in Equation (25) can be replaced bythat has only logarithmic divergence in the lattice regulator. Note that, in the above expression, a general Dirac structure appears, as this methodology is applicable for all types of PDFs. A necessary component of this improved matrix elements is the calculation of the mass counterterm . First attempts to obtain appear in the literature [132, 177] and are based on adopting a static potential [178].
Following the notation of [132], we define the static potential for separation , via an Wilson loop: where is large. The Wilson loop is renormalized as where is due to the corners of the chosen rectangle. Using Equations (57)-(58), one can relate the desired mass counterterm to the static potential, An additional condition is necessary to determine , which can be fixed using where the choice of depends on the scheme of choice. The appearance of an arbitrary scale is not a surprise and is in accordance with the work of R. Sommer [179], suggesting a further finite dimensionful scale that appears in the exponential of Equation (55), based on arguments from Heavy Quark Effective Theory.
A proper determination of requires a nonperturbative evaluation on the same ensemble used for the calculation of the quasi-PDF. This is essential in order to eliminate a source of systematic uncertainty related to the truncation of a perturbative evaluation, which is limited to typically one to two loops. Furthermore, such a quantity can be used for a purely nonperturbative matrix element. Nevertheless, this quantity has been computed to one-loop level in perturbation theory [132, 177] in an effort to qualitatively understand the structure of the power divergence. Within these works, it was demonstrated that such a mass counterterm removes the power divergence to all orders in perturbation theory.
7.2. Lattice Perturbation Theory
The promising results from the first exploratory studies of the quasi-PDFs [46, 47] have led to an interest in developing a renormalization prescription appropriate for Lattice QCD. Several features of the quasi-PDFs have been studied in continuum perturbation theory (see, e.g., Section 5.2), but more recently there appeared also calculations in lattice perturbation theory. Such development is highly desirable, as the ultimate goal is to relate quasi-PDFs extracted from numerical simulations in Euclidean spacetime to standard PDFs in continuum Minkowski space. In this subsection, we highlight three calculations that provided important insights into quasi-PDFs.
7.2.1. IR Divergence in Lattice and Continuum
X. Xiong et al. have computed in [180] the unpolarized quasi-PDF in lattice perturbation theory using clover fermions and Wilson gluons. The calculation was performed in Feynman gauge and included a nonzero quark mass. This allowed the study of the matching between lattice and continuum, and, as discussed in that work, the massless and continuum limits do not commute, leading to different IR behaviors. The calculation contained the one-loop Feynman diagrams shown in Figure 12, where the quark () and internal gluon () momenta are shown explicitly.
(a)
(b)
(c)
(d)
Here, we do not provide any technical details and focus only on the qualitative conclusions, but we encourage the interested Reader to consult [180] for further details. The one-loop results show that a correct recovery of the IR divergence of the continuum quasi-PDFs can only be achieved for and ; the complete continuum quasi-PDF is obtained from . However, it is stressed that this is necessary only for perturbative calculations, as the nonperturbative ones do not contain collinear divergence. This is encouraging and serves as a proof of the matrix elements in Euclidean space being the same as the ones in Minkowski space. Same conclusions have been obtained from [131, 165].
7.2.2. Renormalization of Nonlocal Operators in Lattice Perturbation Theory
M. Constantinou and H. Panagopoulos have calculated in [106] the renormalization functions for the nonlocal operators in perturbation theory in lattice regularization (LR). The calculation was performed to one-loop level, in which the diagrams shown in Figure 12 were evaluated for clover fermions and a variety of Symanzik-improved gluon actions, including Iwasaki and tree-level Symanzik. Note that the schematic representation of the diagrams shown in Figure 12 appears to be the same for dimensional and lattice regularizations, but a calculation in LR is by far more complicated numerically. This is a consequence of the QCD Lagrangian discretization, coupled with the additional divergences that depend on the lattice regulator. The renormalization functions were computed for massless fermions in the scheme for general values of the action parameters and general gauge. The latter has served as a cross-check for gauge-independent quantities. In addition to the calculation in LR, Green’s function of the nonlocal operators has been obtained in dimensional regularization (DR), which, in combination with the corresponding lattice results, gives a direct definition of the renormalization functions in the scheme.
The operator under study includes a straight Wilson line in the direction and has the general form In the above operator, only is to be considered, due to the appearance of contact terms beyond tree level, making the limit nonanalytic. Green’s functions of the above operators are evaluated for all independent combinations of the Dirac matrices, ; that is, Note that the above includes twist-2 operators as well as higher twist. We will later see that it is important to distinguish between the cases in which the spin index is in the same direction as the Wilson line (), or perpendicular to it ().
One of the main findings of this work is the difference between the bare lattice Green’s functions and the -renormalized ones (from DR). This contributes to the renormalization of the operator in LR and was found to receive two contributions, one proportional to the tree level of the operator (), and one that has a different Dirac structure (); that is, where are numerical coefficients that depend on the action parameters. Note that the term proportional to is the one-loop counterpart of the power divergence discussed in the previous subsection, and its numerical coefficient has been computed in [106]. Perturbation theory is not reliable in providing numerical values for mixing and power coefficients; nevertheless it provides crucial qualitative input for the quantities under study.
The conclusion from Equation (63) is that the operator with structure will renormalize multiplicatively only if the anticommutator between and vanishes. This is true for the axial and tensor operators that are used for the helicity () and transversity () PDFs, which have one index in the direction of the Wilson line. On the contrary, the vector current turns out to mix with the scalar, a higher-twist operator. This finding impacted significantly all numerical simulations of the unpolarized PDFs, as they were using this particular operator [46, 100, 107, 108, 181], unaware of the aforementioned mixing. With this work, the Authors proposed to use a vector operator with the spin index perpendicular to the Wilson line and the ideal candidate is the temporal direction () for reasons beyond the mixing. First, the matching procedure between the quasi-PDFs and normal-PDFs also holds for , as it belongs to the same universality class of [182]. In addition, the temporal vector operator offers a faster convergence to the light-cone PDFs, as discussed in [75]. Note, however, that and do not share the same matching formula, with the latter being calculated much later. Detailed discussion on the matching to light-cone PDFs is provided in Section 8.
The work of [106] has led to a number of useful information not only on the renormalization pattern of nonlocal operators, but also on the conversion from a renormalization scheme of choice () to the scheme. This is extracted from the ratio of renormalization functions in the two schemes computed in DR, and multiplies nonperturbative estimates of in order to obtain . Due to the mixing found in the lattice calculation, a convenient scheme which is applicable nonperturbatively is an RI-type [183]. A well-defined prescription within RI-type schemes exists for both the multiplicative and the mixing coefficients, as described in Section 7.3. The is a natural choice for nonperturbative evaluations of renormalization functions, because it does not require to separate finite contributions with tensor structures which are distinct from those at tree level (typically denoted by that appears in the local vector and axial operators in the limit of zero quark mass). The conversion factor has been computed and was used in the renormalization program of the ETMC [151]. The conversion factor shares certain features with the matrix elements; that is, it is complex and symmetric/antisymmetric in the real/imaginary part. A representative example is shown in Figure 13 for the vector, axial, and tensor operators that have a Dirac index in the same direction as the Wilson line.
7.2.3. Nonlocal Operators for Massive Fermions in Dimensional Regularization
G. Spanoudes and H. Panagopoulos [184] have extended the work of [106] presented in the previous paragraph, by examining the effect of nonzero quark masses on the renormalization functions and conversion factors between the and schemes, as obtained in DR at one-loop level. This work was motivated by the fact that lattice simulations are not performed exactly at zero renormalized mass. Of course, one expected that the correction will be very small for the light quarks, but not necessarily for the heavier quarks, which are typically used in dynamical simulations ( and ). In principle, one should adopt a zeromass renormalization scheme for all quarks that requires dedicated production of ensembles with all degenerate flavors (e.g., and ), as typically done for local operators, but this entails additional complications.
Including massive quarks requires a proper modification of the RI-type renormalization conditions, as developed in [184]. Also, in addition to the fermion field renormalization function, the quark mass renormalization is required as well (see Equations (4)-(6)). More interestingly, the conditions for the nonlocal operators must be generalized to account for the more complicated structure of Green’s functions. In particular, it is found that the mixing revealed in [106] extends beyond the anticommutator (: direction of the Wilson line), which still holds for operators with the same flavor in the external quark fields. However, operators with a different flavor give rise to additional mixing, affecting, among other operators, (mixes with ), (mixes with ), and (mixes with ). Depending on the size of the mixing and the simulation setup, a nonnegligible effect may occur in numerical simulations, as all these operators are used in the quasi-PDFs calculations. This is more likely to impact the results extracted using strange and charm in the sea. This includes the first studies (e.g., [46, 47]), but also the more recent work of , in which the Authors use a single ensemble for the renormalization functions and an extrapolation to the chiral limit is not possible for the -factors. Unlike the case of the local operators, where quark mass dependence is negligible, the nonlocal operators exhibit quite visible mass dependence for Wilson lines of length larger than 0.8 fm [185]. However, the mixing is expected to be at most finite and thus not present in the scheme.
As a consequence of the additional mixing, the conversion factors are matrices usually determined in DR, as they are regularization-independent quantities. In [184], the and renormalization functions were obtained by using appropriate conditions on the bare Green’s functions. This is a complicated process and the results can be found in Section III of [184]. Here, we present in Figure 14 the conversion factor for the mixing coefficient between the pseudoscalar and the axial () operators. As can be seen, the mixing is small, but nonnegligible, especially if the flavors involved have mass difference above 100 MeV. The action parameters are given in [184] and the notation is .
7.3. Nonperturbative Renormalization
The progress in the renormalization of the nonlocal operators from lattice perturbation theory has encouraged investigations of nonperturbative calculations. This was supported by theoretical developments proving the renormalizability of the operators under study to all orders in perturbation theory (see Section 5.2). The full development of a proper nonperturbative prescription has been a natural evolution of the knowledge gained from the perturbative calculations, and in particular the pattern identified in [106]. The Authors of this work have proposed an RI-type scheme that was employed by ETMC [151] giving, for the first time, properly renormalized quasi-PDFs. The approach was also adopted by [181] with a slight variation due to a different projection entering the renormalization prescription. The latter was motivated by the fact that the matrix elements of the vector and axial operators have additional tensor structure different than the tree level. We close the discussion on the nonperturbative renormalization with a presentation of an alternative prescription based on the auxiliary field formalism [186].
7.3.1. Scheme
C. Alexandrou et al. [151] have employed a renormalization scheme that is of similar nature as the scheme [183] that is widely used for local operators. Using the renormalization pattern of [106], the Authors developed a nonperturbative method for computing the renormalization functions of nonlocal operators that include a straight Wilson line. In this scheme, one imposes the condition that Green’s functions of the operator must coincide with the corresponding tree-level values at each value of . This approach is also applicable in the presence of mixing, via matrices (: number of operators that mix with each other). The proposed program has the advantage that it eliminates both power and logarithmic divergences at once, without the need to introduce another approach to calculate the power divergence. This is due to the fact that the vertex functions of the operator that enter the RI-type prescription have the same divergences as the matrix elements. The prescription can be summarized as follows for a pair of nonlocal operators, and , assuming they mix under renormalization: According to the above mixing, the renormalized matrix element of , , is related to the bare matrix elements of the two operators via where and are computed in the scheme and then are converted to the scheme, at an energy scale GeV. The renormalization factors can be computed followingwhere the elements of the vertex function matrix are given by the trace In the above equation, is the tree-level expression of the operator . Thus, all matrix elements of can be extracted by a set of linear equations, which can be written in the following matrix form: The complication of the mixing is not relevant for recent calculations of quasi-PDFs, as the vector operator has become obsolete and was replaced by . In the absence of mixing, the above equations simplify significantly and reduce to where is related to the inverse of the vertex function of the operator. Let us repeat that the prescription is applied on each value of independently.
In Figure 15, we show a representative example of the renormalization function of the axial nonlocal operator (), using an ensemble of twisted mass fermions with a clover term () and lattice size . In the left panel of the plot, we overlay the results for the (open symbols) and the (filled symbols) schemes, for the real and imaginary part of the -factor. The momentum source technique [162, 187] was employed that offers high statistical accuracy with a small number of measurements. The scale was set to . As can be seen from the plot, the imaginary part of is smaller than . It is worth mentioning that the perturbative renormalization function in the , as extracted in DR, is a real function to all orders in perturbation theory. Therefore, it is expected that the imaginary part of the nonperturbative estimates should be highly suppressed.
In the aforementioned work, the Authors used several values of the renormalization scales, and each one was converted to the and evolved to 2GeV. Residual dependence on the initial scale was eliminated by an extrapolation to , and the results can be seen in the right panel of Figure 15. An investigation of systematic uncertainties was presented in [151] and upper bounds for uncertainties were estimated.
7.3.2. RI/MOM Scheme
A modification of the RI-type prescription that was first proposed by Constantinou and Panagopoulos [105] was presented by J.-W. Chen et al. in [181]. The main motivation for the modification was the intent to employ a matching procedure that relates the quasi-PDFs in RI/MOM scheme to the light-cone PDFs in , which was developed by Stewart and Zhao [166]. However, both and RI/MOM prescriptions can be used to obtain an appropriate RI-type formula without complication.
Based on the RI/MOM prescription, the vertex function of the operator under study was projected by (instead of the tree level) in order to account for the extra tensor structure included in the vertex function. We note in passing that the difference between and RI/MOM is finite and should be removed by appropriate modification in the conversion factor to the scheme. Besides the different choice in the projector appearing in the RI/MOM prescription, the rest of the setup is equivalent to that of [151]. A minor exception is the fact that the definition of the renormalization functions of [181] is inverse to the one used in [151], which, however, has no implications in the extracted physics. For example, the RI/MOM prescription for the operator that has mixing is given byThe renormalization matrix can be extracted via
In the calculation of the renormalization factors, the Authors used an clover on HISQ ensemble with a volume [188]. The momentum source method was used that leads to high statistical accuracy, and a single RI/MOM renormalization scale () was employed for each nucleon momentum, which corresponds to .
In Figure 16, we show the multiplicative renormalization factor of (red squares) and the mixing coefficient (blue circles). As expected, it is found that the size of the mixing coefficient is about an order of magnitude smaller than the renormalization factor in the large- region. However, the mixing coefficient should multiply the matrix element of the scalar operator that has very large numerical values, leading to a nonnegligible contribution. The mixing is ignored in the rest of the analysis of [181].
Another investigation of [181] aimed at understanding the mixing discussed in [106] using symmetry properties. This was based on the invariance under parity , time reversal , and charge conjugation , where parity and time reversal are generalized into any Euclidean direction. The operator used was which has the advantage that it is either Hermitian or anti-Hermitian. Taking as an example the vector current in the direction of the Wilson line, , one can see how the mixing arises: the “+” (“−”) combination of Equation (76) is anti-Hermitian (Hermitian). In this case, the transformation properties allow mixing with the unity (scalar) operator. Unlike the case of , the other directions of the vector operators do not suffer from mix even for formulations that break chiral symmetry. This conclusion is fully compatible with the findings of [106].
The study of the symmetries was extended to include nonlocal operators including a covariant derivative () or a power of mass () [189]. It was found that operators with and without a covariant derivative may mix even if axial or chiral symmetry is preserved in the formulation under study. In addition, operators may also mix with operators regardless of chiral symmetry breaking. Further details on this analysis can be found in Tables 3 and 4 of [189].
In closing, let us add that a proper determination of the renormalization functions computed nonperturbatively in an RI-type scheme (e.g., the works presented in Sections 7.3.1 and 7.3.2) requires a few improvements. For once, dedicated calculations are needed on ensembles with all degenerate flavor quarks. These ensembles should correspond to the same value of the coupling constants as the ensembles used for the calculation of the hadron matrix elements. For instance, matrix elements obtained on or ensembles should be renormalized using and , respectively. Renormalization functions should then be computed on multiple ensembles with different quark masses, so the chiral limit can be taken. Several values of the RI scale () should be employed for each ensemble in order to take the limit upon conversion to the and at a common scale. This will eliminate residual dependence on the initial RI scale and give more reliable estimates for the renormalization functions. We note that the extrapolation has been performed in the work of [151] (see left panel of Figure 15). Potential improvements could also be a two-loop conversion factor from an RI-type to scheme, and also a subtraction technique of finite effects using one-loop perturbation theory. This method was successfully employed for local operators of different lattice formulations [162, 190]. It is anticipated that both aforementioned improvements will be available in the near future.
7.4. Auxiliary Field Formalism
An alternative proposal for the renormalization of nonlocal operators is based on an auxiliary field method, a formulation also adopted to prove the renormalizability of the operators under study [135]. The use of this approach is not new but originates from other studies in the continuum [172, 173], adopting an auxiliary scalar field results to a pair of operators in an extended theory instead of the usual nonlocal operators. In this case, a renormalization prescription reduces to a three-parameter equation instead of a single equation for each value, which characterizes the RI-type renormalization. J. Green et al. presented in [186] this nonperturbative approach and employed the twisted mass formulation on two ensembles that have pion mass of around 370 MeV and different lattice spacings ( fm), in order to determine the three parameters of the auxiliary field renormalization scheme.
The auxiliary scalar color triplet field () is defined on the line , where is the length of the Wilson line in physical units. The main component of the approach is the replacement of correlation functions with ones from the extended theory including the field, which involve the local color singlet bilinear . The introduction of the auxiliary field requires modification of the action (for details, see [186]), which yields a bare propagator in a fixed gauge background:In the above expression is a straight Wilson line between points and , and the exponent with the mass is an counterterm. One obtains for the operator including the Wilson line, whose renormalization we are seeking,Besides the counterterm , the renormalization functions of the bilinear () and the operator () must be calculated. Due to mixing allowed by the breaking of chiral symmetry, a proper renormalization is in this caseA different basis of operators may be employed to achieve diagonal renormalization in a mixing matrix; that is,As can be seen from the equations above, the renormalization of requires knowledge of the linearly divergent , the log-divergent , and the finite . Note that is of similar nature as the mixing identified in [106]. In addition, this approach is not applicable for , in which case is a local operator and its renormalization can be extracted from standard RI-type techniques.
In the work of [186], the Authors renormalized nucleon matrix elements obtained from two ensembles of twisted mass fermions. For extracting the renormalization functions, they used ensembles of four degenerate quarks () as expected for mass independent renormalization schemes. However, the chiral limit is yet to be taken for this approach. In summary, the three parameters and the auxiliary field renormalization, , are determined by the RI-xMOM conditionswhere is the usual fermion field renormalization obtained from a standard RI-type prescription. As in the case of the nonperturbative schemes described in the previous paragraphs, a prescription is needed to bring the renormalized quasi-PDFs into the scheme and a conversion factor is necessary. This has been computed in DR to one-loop level in perturbation theory and the formula is given in [186].
Here, we present selected results from [186] and Figure 17 showing the quantity in Equation (81) (left panel) for the two ensembles discussed above with fm () and fm (). As can be seen, smearing of the gauge links reduces the statistical noise but more importantly reduces the difference between the two ensembles. This is an evidence of reduction of the linear divergence. In case of no mixing (axial operator ), is not relevant and only and need to be determined. is shown in the right panel of Figure 17 upon conversion to the scheme and evolution to the scale 2 GeV. The one-loop conversion factor removes the bulk of the dependence on the scheme parameter , and the two-loop evolution removes most of the dependence on the scale .
7.5. Other Developments
In this subsection, we review some other developments related to the renormalization of PDF-related operators, in particular the Wilson-line-induced power divergence.
In 2016, the idea of removing such divergence by smearing was proposed by C. Monahan and K. Orginos [129]. It takes advantage of the properties of the gradient flow (GF), introduced by M. Lüscher a few years ago [191, 192] and applied to many problems in Lattice QCD [193]. As shown by Lüscher and P. Weisz in [194], it defines a 4+1-dimensional field theory, wherein the extra dimension is the flow time. The crucial property of GF, proven to all orders in perturbation theory, is that the correlation functions defined at nonzero flow time are finite after usual renormalization of the 4-dimensional theory. As such, GF defines a renormalization scheme with the flow time, , being the renormalization scale. Thus, using GF, one can define smeared quasidistributions, which are finite, in particular devoid of the power divergence from the presence of the Wilson line. Note that GF only regulates the UV behavior, leaving the IR structure intact, which is a prerequisite for factorization. The quasi-PDF results in the GF scheme can be converted perturbatively to other renormalization schemes or directly matched to light-cone PDFs, e.g., in the scheme. The Authors of [129] demonstrated a simple relation between the moments of the smeared quasi-PDF and the renormalized moments of the light-cone PDF. For this relation to be valid and to allow matching to light-cone PDFs, the scales in the problem have to satisfy . Apart from usual higher-twist corrections, , there are also corrections of . The explicit one-loop perturbative analysis of smeared quasi-PDFs was performed, in 2017, by C. Monahan [195]. It was shown that indeed the IR divergences of smeared quasi- and light-cone PDFs are the same. The perturbative computation led in the end to establishing the matching equation that could be used to extract light-cone PDFs from a lattice computation. An interesting aspect also shown by Monahan is that the smeared matrix element satisfies a relation akin to a usual renormalization group equation. This could, in principle, allow a nonperturbative step-scaling procedure to be defined, which would connect lattice-extracted matrix elements to high scales at which matching could be performed with much reduced truncation effects.
Smeared operators are the fundament of another method, introduced in 2012 by Z. Davoudi and M. Savage [196]. It aims at calculations of arbitrarily many moments of PDFs or other structure functions, that could, in principle, allow us to reconstruct the full distributions. The main idea is to avoid the power-divergent mixings with lower-dimensional operators in higher moments by removing their source, the breaking of rotational invariance by the lattice. The paper considers a mechanism for the restoration of this symmetry in the continuum limit of lattice field theories, in particular the theory and QCD. In general, the interpolating operators that are used to excite a hadron do not have definite angular momentum; i.e., it is not possible to assign a well-defined angular momentum to a lattice state and the latter is a linear combination of infinitely many different angular momentum states. The essence of the approach of [196] is to construct appropriate operators on the hypercubic lattice with maximum overlap with states with definite angular momentum in the continuum. Such operators are constructed on multiple lattice sites using smearing that renders the contributions of both lower- and higher-dimensional operators subleading and totally suppressed in the continuum limit. The Authors performed detailed calculations in the theory demonstrating the mechanism. For the QCD case, things are complicated by the gauge symmetry. However, Davoudi and Savage showed that the idea can also be applied for this case, relevant for moments of partonic functions. Apart from smearing of the operators, the essential ingredient is tadpole improvement. Recently, this approach has been revisited and exploratory numerical results were recently presented [197].
We finalize by shortly discussing one more method of dealing with the power divergence related to the Wilson line in quasidistributions. In 2016, H.-n. Li proposed [198] to modify the definition of such distributions by using “nondipolar” gauge links, i.e., two pieces of links oriented in orthogonal directions. He showed with an explicit calculation of one-loop corrections that the linear divergence of the standard quasiapproach (with dipolar links) is absent in such a case and the IR region is untouched. In general, the hadron boost direction needs to differ from the direction of the Wilson line links to avoid the linear divergence. However, due to developments in the renormalization of the power divergence (and other divergences present in quasidistributions), in particular the full nonperturbative renormalization, this interesting idea of Li has not been implemented in numerical calculations. Clearly, the implementation itself is possible, but much less practical than with straight links for the standard definition. It is also worth mentioning that [198] discussed also a potential problem with two-loop factorization of standard quasi-PDFs (but absent in the nondipolar ones). The power divergence in such setup induces an additional collinear divergence at the two-loop order, rendering the matching kernel IR divergent and breaking the factorization. However, the problem does not appear if the power divergence is properly renormalized [165].
8. Matching of Quasi-PDFs to Light-Cone PDFs
In this section, we focus on the matching from quasi-PDFs to light-cone PDFs. Since the inception of LaMET, there has been a lot of effort devoted to understanding many aspects of this procedure. In particular, the first matching paper [98], discussed in Section 3, considered the nonsinglet quark quasi- and light-cone PDFs in a simple transverse momentum cutoff scheme. Later work concentrated on matching from different renormalization schemes to the scheme, on the issue of particle number conservation and on observables different than nonsinglet quark PDFs, in particular gluon PDFs, singlet quark PDFs, GPDs, TMDs, and meson DAs. We review all of the below and we also include a discussion on the matching of pseudo-PDFs/ITDs.
For convenience, we repeat here the general factorization formula for the matching:where is the common factorization and renormalization scale, and the second argument of the matching kernel emphasizes that the relevant momentum is that of a parton.
Let us first briefly revisit the early attempt to remove the Wilson-line-related power divergence, discussed in Section 7, from the point of view of the matching process. We will then move to the presentation of the matching of -renormalized and RI-renormalized quasi-PDFs. T. Ishikawa et al. discussed in [132] that the counterterm that subtracts this divergence to all orders in the coupling can be provided by an independent lattice observable that shares the same power divergence as the nonlocal operator defining the quasi-PDF. It was noted that a natural and simple choice for such an observable is the static potential. The Authors calculated the matching (to PDFs in the UV cutoff scheme) in one-loop lattice perturbation theory for the case of naive fermions. The idea was also followed in [177], where J.-W. Chen, X. Ji, and J.-H. Zhang defined improved quasi-PDFs with the power divergence, calculated from, e.g., the static potential, subtracted. The modification amounts to multiplication of bare matrix elements by an exponential factor, . The matching formulae of [98] are then modified by ignoring the terms containing the cutoff .
8.1. Matching of Nonsinglet Quark Quasi-PDFs to the Scheme PDFs
One of the possibilities of renormalizing the quasi-PDF is to obtain it in the scheme. Obviously, this scheme cannot be directly applied on the lattice and, hence, nonperturbative renormalization of lattice matrix elements proceeds via an intermediate scheme, like a variant of RI (see Section 7). Having renormalization functions in such an intermediate scheme, one then converts them perturbatively to the scheme and evolves to some reference scale, like 2 GeV. The last step is the Fourier transform that yields the quasi-PDF in the scheme.
The first paper that considered the matching from quasi-PDF was [142] by W. Wang, S. Zhao, and R. Zhu. The Authors presented complete matching for quarks and gluons, that we discuss more below in Section 8.2. For the case of nonsinglet quark PDFs (with or Dirac structure), it was found that the change with respect to [98] is simple; the terms with the transverse momentum cutoff do not appear and there is a modified polynomial dependence in the physical region of quasi-PDFs. Explicitly, the matching kernel readswhere we now use the notation with plus functions at over some domain of integration , defined as
However, one more issue remained unresolved for the to matching. Namely, the self-energy corrections have a UV divergence in the limit (cf. Equation (85)). Thus, the form of the matching kernel in [142] still needs a cutoff for the -integration. The issue was addressed by T. Izubuchi et al. in [83]. The aforementioned divergence can be canceled by adding a term (for ) or (for ) to the self-energy corrections. In the scheme, another term arises from this modification, outside of the integral sign, and finally the matching kernel readswhere for and for or . Note that the polynomial term in the physical interval agrees with the one of [142] when . Equation (87) is the pure expression for the matching kernel. However, it violates particle number conservation; i.e., after the matching process. Moreover, the violation grows for increasing . To satisfy particle number conservation, the Authors proposed a modified scheme, the so-called “ratio scheme.” It is a modification of the scheme, in which the problem is avoided by using pure plus functions:In this scheme, all regions in the -integration of the plus functions contain the same term and no additional term appears. Formally, this is a different renormalization scheme and, hence, the quasi-PDF used in the matching procedure needs to be renormalized in this scheme. This requires a relatively simple modification of the perturbative conversion from the intermediate renormalization scheme to :This factor simply multiplies the conversion factor or the -factors.
Alternative procedure was used in [124] by ETMC. Similarly to the ratio scheme, the matching kernel contains only pure plus functions:It amounts to the kernel of Equation (87) but without the terms outside the plus functions and without the additional -dependent term outside of the integral and, thus, satisfies the particle number conservation requirement by construction. Similarly to the ratio scheme, it is a modification of the scheme (that we denote by in the superscript), thus requiring modification of conversion. In the procedure used by ETMC [124], this conversion modification was not taken into account, on grounds that the modification of is done only in the unphysical region and it disappears in the infinite momentum limit. After the publication of [124], ETMC has calculated the required conversion modification that was presented in the results of [51]. More details can be found in [185]. As anticipated, the effect is numerically very small and the ensuing light-cone PDFs are compatible with the ones obtained from the simplified procedure. This is in contrast with the ratio scheme, wherein the modification of the physical region in the matching kernel, combined with the factor of Equation (89), brings about large modification of the quasi-PDF and the final light-cone PDF [185].
The matching kernel for transversity PDFs () for the matching has been calculated by ETMC in [48], following the same method to preserve particle number. Thus, it also needs the conversion modification that will be shown in [185]. Explicitly, it readsThe formula for the transversity case is not the same as the unpolarized and helicity distributions due to the different splitting function, different polynomial dependence in the physical region, and different term added in the nonphysical regions to renormalize the UV divergence in the self-energy corrections.
An alternative way of bringing the results from the intermediate RI renormalization scheme to the scheme is to match directly the RI-renormalized quasi-PDFs onto the light-cone PDFs. This way was advocated by I. Stewart and Y. Zhao [166], including the derivation of the relevant formulae. Such one-step procedure can be used as the sole means of obtaining light-cone PDFs or compared to the two-step procedure (first conversion to and evolution to a reference scale and matching as the second step), with differences in the final PDFs taken as a measure of systematic uncertainty. Both procedures have been derived to one-loop order in perturbation theory, but clearly they can differ in the magnitude of neglected higher-order contributions. The derivation of the RI matching is somewhat more complicated than that for the case. Stewart and Zhao presented it (for the or Dirac structures) in the general covariant gauge, including the practically relevant case of the Landau gauge, typically implemented on the lattice. They also showed a detailed numerical study of the dependence on the choice of the gauge and on the initial and final scales. While the matching has only one scale involved, the RI case depends on three scales: the final scale and the two scales of the RI scheme: the overall scale and the scale defined by the momentum in the -direction. Explicit checks showed that, when aiming at a result at some reference scale, the dependence on the intermediate RI scales is rather small. It is important to note that the RI conserves the particle number and also that the problem with the UV divergence in self-energy corrections does not appear, since the RI scheme introduces a counterterm to the quasi-PDF that cancels this divergence. Results for the case and for transversity PDFs matching were presented in [49, 199] by Y.-S. Liu et al. For final RI matching formulae, we refer to the original publications.
8.2. Matching of Other Quasidistributions and Pseudodistributions
In this section, we review other developments in the matching of quasidistributions to their light-cone counterparts. We also shortly discuss the matching process for the pseudo-PDFs/ITDs.
GPDs. Apart from PDFs, also other kinds of parton distributions can be accessed on the lattice via LaMET. Already in 2015, matching was worked out for GPDs, for (nonsinglet) unpolarized and helicity in [115] and for transversity in [200]. In both papers, the transverse momentum cutoff was used, as in the first paper for the matching of PDFs. The lattice matrix elements are extracted in a similar way as for quasi-PDFs, but there is momentum transfer in the boost direction between the source and the sink, . The obtained quasi-GPD can be decomposed into two functions, and , for the unpolarized case, and , , for helicity (chiral-even), and four functions , , , , for transversity (chiral-odd), where the additional variables with respect to standard PDFs are and . For , the , , and quasifunctions become the standard quasi-PDFs and all the functions and have no quasi-PDF counterparts. After matching, the -integrals of unpolarized and give the Dirac and Pauli form factors and , respectively. The -integrals of helicity and give the generalized axial and pseudoscalar form factors and . Finally, the first moments of transversity GPDs give the generalized tensor form factors , , , and , for , , , and , respectively. In the papers [115, 200], it was shown that the matching is nontrivial for the functions , , and and reduces to the matching for the corresponding quasi-PDFs in the forward limit, as expected. In turn, the matching kernel for all the functions is a trivial -function at leading order in the coupling. The fourth transversity quasi-GPD, , is power suppressed by the hadron momentum and omitted at leading power accuracy. We refer to the original publications for the final matching formulae. It is worth mentioning that, for quasi-PDFs, the formulae decompose into three intervals in , the physical one and two nonphysical ones outside of the partonic support for , whereas for quasi-GPDs there are, in general, four intervals with different matching functions for the physical ERBL () and DGLAP ( and ) regions.
Complete Matching for Quark and Gluon PDFs. In 2017, the first calculation of the matching of gluon quasi-PDFs to light-cone PDFs was done by W. Wang, S. Zhao, and R. Zhu [142]. This paper was already discussed in the previous subsection in the context of matching for nonsinglet quark PDFs. However, the aim of the paper was more broad, to consider the complete matching for quark and gluon quasi-PDFs. The Authors used two ways to regulate UV divergences: the UV cutoff scheme and DR, and also two ways for IR divergences: finite gluon mass and offshellness. The gluon quasi-PDF was defined as the Fourier transform of a boosted nucleon matrix element of two gluon field strength tensors displaced by length and connected with a Wilson line in the adjoint representation, , with a sum over transverse directions, . The Authors calculated the one-loop expressions for gluon quasi-PDFs and light-cone PDFs in the UV cutoff scheme and confirmed they have the same infrared structure, with IR divergences present only in the physical region, as expected. They also noted the presence of a linear divergence in the quasidistribution, however, related in the gluon case not only to the presence of the Wilson line. Having performed the calculation in the UV cutoff scheme, they pointed to a difficulty in the self-energy diagram, coming from the breaking of gauge invariance by this scheme. This motivated further computations in DR, performed for both quark and gluon quasi- and light-cone PDFs. They considered all possible cases of quark-in-quark, gluon-in-gluon, gluon-in-quark, and quark-in-gluon distribution functions, required for the complete matching. The quark-in-quark matching, the only one relevant for the nonsinglet quark distributions, has already been described above; see Equation (85). Together with the derived equations for the other cases, one is ready to write the final matching formula:where () denotes the light-cone (quasi-) distribution. The indices and the four cases mentioned above correspond to matching kernels , , , and , respectively. The matching equation implies mixing under matching between quark and gluon distributions, which can only be avoided in nonsinglet quark distributions. Finally, the Authors derived evolution formulae for quasidistributions that turned out to be the DGLAP evolution equations of light-cone PDFs.
In a follow-up work [136], W. Wang and S. Zhao considered in more detail the issue of the power divergence in quasigluon PDFs; see Section 5.2.2 for more details from the point of view of renormalizability. As remarked above, linear divergences exist also in one-loop diagrams without a Wilson line, which means that the divergence cannot be absorbed into the renormalization of the Wilson line. The adopted definition of the gluon quasi-PDF was slightly modified with respect to [142] by extending the sum in from the transverse directions to all directions except the direction of the boost; i.e., . The calculation of one-loop corrections to quasigluon distributions was performed in a UV cutoff scheme, with the cutoff interpreted as the lattice cutoff. The Authors included diagrams arising in lattice perturbation theory (counterterm from the measure in the path integral and quark and ghost tadpoles) that preserve the gauge invariance, broken in the naive cutoff scheme. The main result of the paper, derived in the auxiliary field formalism, is that the linear divergences can be renormalized by considering the contribution from operator mixing (only with certain gluonic operators, i.e., no mixing with quark quasi-PDFs occurs) and the mass counterterm of the Wilson line. This allowed the Authors to define an improved quasigluon PDF with matrix elements multiplied by , where the mass counterterm can be determined nonperturbatively, and with a subtraction of the mixing calculated in perturbation theory. In addition to the one-loop calculation, they discussed two-loop corrections and conjectured that they hold to all orders. Finally, they provided the formula for the one-loop matching of the improved gluon quasi-PDF, which is IR finite and free from the linear UV divergence.
The proof of renormalizability to all orders was indeed provided a few months later by the same Authors, together with J.-H. Zhang, X. Ji, and A. Schäfer [137] (see Section 5.2 for more details on this paper and another proof of renormalizability of gluon quasi-PDFs [138]). From the point of view of matching, the important contribution of this paper was to confirm that the conclusions of [136] hold when using gauge-invariant DR instead of the UV cutoff scheme. Moreover, it was pointed out that one can construct gluonic operators that are multiplicatively renormalizable; i.e., they evince no mixing under renormalization. However, mixing still occurs at the level of matching, as in Equation (92). The Authors wrote schematic matching equations for the proposed nonperturbatively renormalized gluon quasi-PDFs in the RI/MOM scheme, postponing the calculation of the matching kernels RI to a forthcoming publication. The latter computation, as well as the alternative possibility of RI conversion and matching, will open the prospect of obtaining light-cone gluon PDFs in the scheme from the lattice.
TMDs. Yet another important class of partonic distributions that can, in principle, be accessed on a Euclidean lattice is TMDs. The quasi-TMDs were considered already in 2014 by X. Ji and collaborators in [201]. They performed a one-loop perturbative calculation of quasi-TMDs in the Drell-Yan process. The crucial subtlety that makes the TMD case much more cumbersome than the PDF case is the subtraction of the soft term. It needs to be constructed in such a way to make it computable on the lattice. It is related to the presence of a light-cone singularity in TMDs. The unsubtracted matrix element for quasi-TMDs, , is defined as the correlation between a quark and an antiquark in a boosted nucleon, with the quark fields spatially separated by a distance and connected by two gauge links: one going from the quark field to infinity (for Drell-Yan) and the second one from infinity to the antiquark (in the covariant gauge; in the axial gauge an explicit link at infinity is additionally needed). The TMD depends on the longitudinal momentum fraction, , and the transverse momentum, , where the latter is often exchanged for the impact parameter, , via a two-dimensional Fourier transform. Having defined the quasi-TMD, the Authors proposed a lattice-calculable subtraction of the soft factor. The latter was conjectured to also play an important role in the two-loop matching for quasi-PDFs, where it could be handled similarly. Further, they proceeded with the derivation of the one-loop formulae and demonstrated the one-loop factorization. Finally, they also considered the TMD evolution (Collins-Soper evolution) in the scale related to the hadron momentum or the hard scale of the scattering process.
Early in 2018, X. Ji et al. reinvestigated quasi-TMDs [202]. They considered gauge links of finite lengths (staples), instead of infinite ones. Moreover, they redefined the subtraction of the soft factor, since the one defined in [201] could have practical implementation difficulties on the lattice. With this modified subtraction and finite-link TMDs, the Authors could show that the so-called pinch pole singularities are regulated. The new subtraction leads to an additional term in the one-loop computation of the quasi-TMD. Before establishing the matching formula, resummation of large logarithms needed to be performed to avoid scheme dependence in regulating light-cone singularities. This could be done using the Collins-Soper evolution derived in [201]. Finally, the matching equation was given to the TMDs in the standard TMD scheme.
Very recently, a third paper considering quasi-TMDs appeared by M. Ebert, I.W. Stewart, and Y. Zhao [203]. TMDs depend on two scales, the virtuality scale and the scale introduced above. Evolution in the former can usually be done fully perturbatively, á la DGLAP. For the latter, the (Collins-Soper) evolution involves the impact parameter dependent anomalous dimension, ( – parton flavor index), and becomes nonperturbative for transverse momenta of the order of , even for . The focus of this paper was on this aspect. The Authors proposed a method of a first-principle nonperturbative determination of , using the quasi-TMD formalism. They defined the quasibeam function (unsubtracted quasi-TMD) with finite-length () gauge links that can be related to the corresponding light-cone beam function. However, for the soft function that provides the subtraction of the soft term, they argued that a straightforward definition of a quasianalogue is not possible, since the Wilson lines of the soft function involve both light-cone directions and would require opposite boosts to be recovered from Wilson lines in the spatial directions. A detailed study of this aspect was postponed to a forthcoming publication. For this paper, the Authors introduced a function that describes the missing IR physics and -dependence, . This function removes the linear divergences in the Wilson line self-energy and an explicit form that cancels all divergences in may be used in the form proposed in [202]. The crucial aspect for the extraction of is that the factor cancels in the ratios of quasi-TMDs defined at different nucleon boosts. The matching between quasi-TMDs and light-cone TMDs can be spoiled by the issue in the soft function and the Authors introduced a function expressing the mismatch between quasi- and light-cone soft functions and allowed it to be nonperturbative. They expressed the quasi-TMD in terms of the light-cone TMD for the nonsinglet case via the perturbative kernel (matching between quasi- and light-cone beam functions), the unknown and the Collins-Soper anomalous dimension. Knowing the that matches the IR physics of the light-cone soft function, the could also be calculated perturbatively. The interpretation of the matching equation differs from the analogous one in [202], wherein the analogue of is assumed to be fully perturbative, which is claimed to be incorrect due to missing the nonperturbative physics when is of order . Taking the ratio of two matching equations, the -independent factor drops out and one can extract the anomalous dimension based on the perturbative matching relation between the quasi- and the standard beam functions. The method was illustrated by an explicit one-loop computation. It was also remarked that it is restricted to the nonsinglet quark channel, because of mixings between singlet quarks and gluons under matching (see the previous paragraph).
Meson DAs. Another type of observables that can be considered in the framework of LaMET is meson DAs. They are defined as vacuum-to-meson matrix elements of the same operator as for PDFs, quark, and antiquark connected with a Wilson line, with -structure of, e.g., , for pseudoscalar mesons. They are easier to calculate, since they require only two-point functions, as the pion is not annihilated in the matrix element. The matching can be extracted as a limit of the matching formula for GPDs by crossing the initial quark to the final state and it was extracted for the first time (for the pseudoscalar case) in the paper [115] commented on above. We refer to this paper for explicit matching formulae in the transverse momentum cutoff scheme.
Further, (pseudoscalar) meson mass corrections were calculated analytically in [123], yielding an infinite series in which the few first terms are enough to take into account for practical application.
The heavy quarkonium case was considered in [125] by Y. Jia and X. Xiong, with the one-loop corrections to both quasi- and light-cone DAs computed in the framework of NRQCD factorization. The matching for meson DAs and PDFs was also analyzed by Jia et al. within two-dimensional QCD [126]. In both papers, the UV divergences were regulated with a transverse momentum cutoff, interpreted as a renormalization scale. For more details about these two papers, see Section 4.4.
The matching for vector meson DAs was also considered [204], by J. Xu, Q.-A. Zhang, and S. Zhao. They derived the formulae in the UV cutoff scheme and in DR (with subtraction), both for longitudinally and transversely polarized mesons.
Recently, the matching for meson DAs was also obtained for the case of RI-renormalized quasi-DAs to bring them into -renormalized light-cone DAs [205], by Y.-S. Liu et al. They considered the cases of pseudoscalar, as well as longitudinally and transversely polarized vector mesons. The quasi-DA can be renormalized with the same renormalization factors as the corresponding quasi-PDF, in a variant of the RI/MOM scheme. The one-loop calculation of the matching relation proceeded along the lines of analogous computations for quasi-PDFs, first done in [166], and we refer to the original paper for the final formulae.
Pseudo-PDFs. The one-loop corrections to pseudo-PDFs were first considered in [87] by K. Orginos et al., in the leading logarithmic approximation (LLA), appropriate to study the dependence. In the LLA, pseudo-PDFs are simply related to the PDFs at the scale : , where plays the role of the renormalization scale for the pseudodistribution. The full one-loop corrections to pseudo-PDFs were calculated by X. Ji, J.-H. Zhang, and Y. Zhao [165] and also by A. Radyushkin [80]. Further insights into the structure of these corrections were given in [81] of Radyushkin, which also contains the explicit matching between reduced ITDs and standard light-cone PDFs in the scheme. The matching was also simultaneously computed by two other independent studies: J.-H. Zhang, J.-W. Chen, C. Monahan [82], and T. Izubuchi et al. [83], and the preprints were made available for all three papers almost simultaneously. After initial discrepancies due to finite terms, all three results agree with one another.
The matching of pseudo-PDFs is, to some extent, simpler than that for quasi-PDFs, since there are no complications related to the nonperturbative renormalization of the pseudo-PDF when taking the ratio of matrix elements to construct the reduced ITD. Crucially, taking the ratio does not alter the IR properties and the factorization framework can be applied, as in the case of matching quasidistributions. We write here the final matching formula in the notation of [81]:where is the light-cone ITD, at Ioffe time and renormalized at the scale in the scheme, is the pseudo-ITD at the scale , and is the Altarelli-Parisi kernel. The first term under the integral corresponds to the LLA result used in [87], i.e., the invoked above multiplicative difference between the pseudo-ITD and scales. The term containing leads to a large negative contribution and causes that the -dependence of vertex diagrams involving the gauge link is generated by an effective scale smaller than . This can be seen by rewriting the matching equation in such a way that the logarithmic term has an argument ; i.e., it involves instead of . In this way, the evolution is governed by this combined logarithm instead of simple , which leads to rescaling with a coefficient , numerically found to be relatively large, around 4 for the setup of [87] (cf. its LLA value of approximately 1.12).
In [83], the relation between quasi-PDFs, pseudo-PDFs, and ITDs was emphasized. This relation implies that their matching involves a unique factorization formula that involves small distances and large nucleon boosts. For these reasons, Izubuchi et al. claim that LaMET and pseudo-/Ioffe-time distribution approaches are, in principle, equivalent. However, it should be noted that the structure of one-loop corrections is different between them and, obviously, the lattice systematics are not equivalent. Because of this, in the absence of all-order (or nonperturbative) matching formulae and under realistic lattice situations, it seems more proper to view them as complementary approaches that aim at the same physical observables.
9. Quark Quasi-PDFs of the Nucleon
The preliminary studies presented in Section 3.2 have evolved based on the progress on various aspects of PDFs, including simulations with improved parameters, renormalization, and choice for the matching procedure. It is the goal of this section to present the advances in the numerical simulations, including a critical discussion on the systematic uncertainties outlined in Section 6. We first present results on ensembles with the quark masses tuned to produce a pion mass larger than its physical value, and we extend the discussion for the simulations with physical values for the quark masses (physical point). To avoid repetition, let us point out that all the works presented here correspond to the isovector flavor combination , which receives contributions only from the connected diagram (up to cutoff effects).
9.1. Simulations at Unphysical Quark Masses
Once the nonperturbative renormalization of the nonlocal operators with straight Wilson line has been developed and presented to the community [105] (see Section 7.3), the first implementation for the quasi-PDFs appeared in the literature in 2017, by ETMC [151] in the scheme, and a modification of the proposal in the RI/MOM scheme by the collaboration [181].
In the original proposal for the nonperturbative renormalization [151], Alexandrou et al. (ETMC) applied the renormalization prescription on their previous work of [108] for an ensemble with MeV (see Section 3.2 for a discussion on the simulations). This employed large-statistics results for nucleon momentum 1.42 GeV and source-sink separation of about 0.98 fm, to demonstrate the effect of the renormalization for the helicity PDFs, which has a multiplicative renormalization, 9. The renormalization function was extracted in the scheme at a scale of 2 GeV, and the remaining dependence on the RI scale () was reduced by an extrapolation: where is the desired quantity. In the work of [151], the fit was performed in the range . One technical consequence of the renormalization is the behavior of the renormalized matrix element in the large- region: while the real (imaginary) part of the bare matrix element decays to zero for (), the renormalization function grows exponentially due to the power divergence. This leads to the unwanted effect of enhanced values for the matrix elements that are almost compatible with zero within uncertainties. This effect is propagated to the quasi-PDF (with the truncation of the integration limits of the Fourier transform), as well as the final extraction of the PDFs. Let us also add that, in the RI-type renormalization prescription, each value of is renormalized independently. More discussion on this systematic effect can be found in Section 6.3.
A Fourier transform is applied on renormalized matrix elements leading to -dependent quasi-PDFs, followed by the matching procedure and target mass correction to finally extract the light-cone PDFs. The obtained helicity PDF from the aforementioned work is shown in Figure 18. To demonstrate the effect of a proper renormalization, we compare the PDF computed with either fully renormalized matrix elements (blue band) or matrix elements renormalized with the local axial vector current renormalization function for all values (magenta band) that was previously used in [108]. As can be seen from the figure, the blue band has a form that is closer to the phenomenological PDFs, as compared to the magenta band. There is also a visual improvement for the antiquark region (), which, however, should not be considered conclusive, as this region is, to date, unreliably extracted.
Despite the improvement from previous works on quasi-PDFs, a number of further improvements were still necessary at this point, as described in Section 6.3. Let us also add that [151] employed the matching formula of Xiong et al. [98] that was obtained in the transverse momentum cutoff scheme and was later replaced by matching formulae calculated in dimensional regularization (see Section 8).
In the work of J.-W. Chen et al. () presented in [181], a nonperturbative renormalization was also applied, using the RI/MOM scheme. They focused on results for the unpolarized PDF, which however uses the “” vector operator. The mixing present in this operator was ignored under the assumption that it is small. Indeed, the mixing coefficient is smaller than the multiplicative factor (see Figure 16), but the scalar operator (that mixes with “”) is expected to be sizable. This can also be seen from the extraction of the scalar charge using the same ensemble, which has the bare value [232]. The RI/MOM renormalization scale is fixed to the nucleon momentum, , which also appears in the matching and, thus, cancels to leading order. Even so, residual dependence on this scale can be nonnegligible (estimated to up to 10% based on studies with ultralocal operators ()) and an extrapolation would be desirable; otherwise this systematic uncertainty cannot be assessed.
The renormalization function of this work was used on the results obtained in [107] (for an ensemble with MeV) for the unpolarized PDF, together with the matching obtained by I. Stewart and Y. Zhao [166] (see also Section 8.1). The matching formula of the latter work was the first one obtained for renormalized quasi-PDFs in the RI scheme matched to the light-cone PDFs in the scheme. A different kind of comparison between lattice and phenomenological data is presented in Figure 19. The renormalized matrix elements for the unpolarized case are compared to phenomenological data [214] on which an inverse Fourier transform and matching have been applied to bring them to coordinate space. This procedure was applied on the central values and, thus, statistical and systematic uncertainties are absent. It is found that the lattice data have a narrower peak around (real part) and are not compatible with the CJ15 data for large values of the Ioffe time, . Note, however, that the lattice data carry very large uncertainties for the large- region that prevents proper comparison. In addition, there are concerns on whether such a comparison is meaningful due to higher-twist effects.
A recent effort to quantify systematic uncertainties was presented by Y.-S. Liu et al. () in [199], using an ensemble with pion mass value of about MeV [188]. Clover valence fermions were employed on an HISQ ensemble [233]. The lattice spacing is fm, and the volume has a spatial extent fm. In this work, Liu et al. computed the unpolarized PDF with nucleon momentum 1.7, 2.15, and 2.6 GeV, and source-sink separations that correspond to 0.60, 0.72, 0.84, 0.96, and 1.08 fm. The main goal of this work was to study uncertainties related to excited states, the nonperturbative renormalization, and the matching to light-cone PDFs. For the Fourier transform to momentum () space, the Authors used the derivative method [160], which is based on an integration by parts, instead of the standard Fourier transform. In this procedure, the corresponding surface term is neglected (see Section 6.3 for details), which carries systematic uncertainties; the latter is not addressed in this work.
Possibly the largest systematic effect comes from the excited states contamination, which is sensitive to the pion mass (worsens for simulations at the physical point) [39], an issue that also appears in matrix elements of local operators. In fact, the situation for the nonlocal operators entering the quasi-PDFs calculation is more severe, as the number of excited states increases with an increase of the nucleon momentum. The effect of excited states can be understood using different analysis methods, as presented in Section 6.1, with the single- and two-state fits being crucial for identifying the ground state of the nucleon. This is particularly important for nonlocal operators that are limitedly studied and are less understood than other hadron structure quantities. Ideally, one should perform a combined analysis with source-sink separations higher than 1 fm. The need for two different analysis techniques is to ensure that the dominant excited states are eliminated by achieving convergence between different techniques. In addition, a single-state fit (applied on each source-sink separation separately) gives important information on the statistical uncertainties of the lattice data. Such information is not to be underestimated, as multistate fits will be driven by the most accurate data. Since statistical noise increases exponentially with the source-sink separation, the most accurate data typically correspond to small separations, which are severely affected by excited states contamination. In the work of Liu et al., the analysis is exclusively based on two-state fits using either all five separations, or four/three largest ones. The Authors did not provide any details on the statistics used in this work, nor the statistical accuracy of the data on each separation, leading to an inadequacy in the quality of their analysis procedure.
The work of [199] addressed systematic uncertainties related to a convenient choice of an RI-type scheme, by examining two possible projectors in the renormalization prescription. This was motivated by the fact that Green’s function, , of the unpolarized operator has additional tensor structures; that is,where ’s are form factors. The minimal projection only projects out , while an alternative choice for the projection is [166], which we call the projection, leading to the conditions An appropriate matching formula to the scheme had been also derived for each RI scheme, and it was concluded that the minimal projector leads to better controlled final estimates, shown in Figure 20, compared with global analysis PDFs [210–212]. The Authors reported reasonable agreement with global analyses in small- and large- regions, while the slope of the lattice data at intermediate -values is different, possibly due to uncertainties in the derivative method for the Fourier transform. The pion mass of the ensemble used for the data is 310 MeV, making the comparison only qualitative.
9.2. Simulations at Physical Quark Masses
One of the highlights of the current year is the appearance of lattice results on quasi-PDFs using simulations at the physical point 10 by ETMC [48, 124] and [49, 218, 219]. Unlike previous studies, these results include proper nonperturbative renormalization and an appropriate matching procedure, for the unpolarized, helicity, and transversity PDFs. We note that, in these works, the use of the Dirac structure parallel to the Wilson line, , has been abandoned due to the mixing discussed in Section 7.2.2 and is replaced by the vector operator with the Dirac structure in the temporal direction, . Here we outline the most important results from each work.
9.2.1. Unpolarized and Helicity PDFs
The work by C. Alexandrou et al. (ETMC) presented in [124] is the first complete calculation of ETMC with several of the systematic uncertainties under control: simulations at the physical point, nonperturbative renormalization, matching to light-cone PDFs computed in dimensional regularization in the scheme. The ensemble corresponds to twisted mass fermions (at maximal twist) with a clover improvement [234]. The ensemble has a lattice spacing of 0.093 fm, lattice spatial extent of 4.5 fm (), and a pion mass of 130 MeV. The nucleon matrix elements of the nonlocal vector and axial operator were computed for three values of the momentum, 0.83, 1.11, and 1.38 GeV, and employ momentum smearing on the nucleon interpolating field [109], that leads to a better signal for the high momenta at a reasonable computational cost (see also Section 6.2 for more details about optimization of the lattice setup). In addition, stout smearing [150] was applied to the links of the Wilson line entering the operator, that reduces the power divergence, and it was checked that different numbers of steps for the stout smearing lead to compatible (almost equivalent) renormalized matrix elements.
A large number of configurations is necessary to keep the statistical uncertainties under control, in particular, as the nucleon momentum increases. The work of [124] analyzed 9600, 38250 and 58950 independent correlators for the momenta 0.83, 1.11, and 1.38 GeV, respectively, so that statistical uncertainties are at the same level. A first study of excited states contamination was presented using only two values of the source-sink separation, 0.93 and 1.12 fm, and demonstrating that within statistical uncertainties the matrix elements are compatible. Nevertheless, a dedicated study of excited states is missing from the presentation and was recently completed [216], concluding that the separation 1.12 fm is sufficient for a nucleon momentum of about 1.5 GeV. We will discuss this investigation below.
The renormalization was performed according to the procedure outlined in Section 7.3.1 and the quasi-PDFs were extracted by the standard Fourier transform. The matching formula used in the work of ETMC was a modified expression with respect to the one suggested in [83] (see discussion is Section 8.1), that preserves the normalization of the distribution functions. However, there is a small mismatch in the renormalization procedure and the matching process, as the conversion factor brings the quasi-PDFs to the scheme, while the matching assumes that the quasi-PDFs are given in the scheme. Preliminary investigation showed a small effect, but this mismatch adds to the overall systematic uncertainties. A follow-up work by ETMC eliminated this uncertainty by computing the quasi-PDFs in the proper scheme [51, 185]. Nucleon mass corrections were applied according to the formulae of [107].
In Figure 21, we show the final results for the unpolarized (left) and helicity (right) distributions for the three values of the nucleon boost. For qualitative comparison, we also include the phenomenological determinations: CJ15 [214], ABMP16 [213], NNPDF3.1 [211], DSSV08 [206], NNPDF1.1pol [207], and JAM17 [215]. The Authors reported that the increase of the nucleon momentum shifts the lattice data towards the phenomenological results. For the unpolarized PDF, the two largest momenta give compatible result, while it is not the case for the helicity PDF. For the latter, there is better agreement with phenomenology, compared to the unpolarized case. As seen from the plots, the large- region suffers from the so-called oscillations that are unphysical. This result from the fact that the bare matrix element does not decay to zero fast enough for large (due to finite momentum), while the renormalization grows exponentially. It is worth mentioning that the oscillations become milder as the momentum increases from 0.83 GeV to 1.38 GeV. It is clear that there are several aspects of the current studies to be improved and the removal of the oscillations is one of them. For this to be achieved while systematic uncertainties are under control, different directions must be pursued, for instance, new techniques that can contribute to a reduction of the gauge noise in the correlators.
An interesting discussion presented in [124] is the comparison between results at the physical point and results from an ensemble with pion mass of about 370 MeV [108] (labeled as B55), as shown in Figure 22. The nucleon momentum is the same for both ensembles (≈1.4 GeV) and a clear pion mass dependence is observed. This is not surprising, as similar pion mass dependence is found in the first moment, , computed with other techniques in Lattice QCD.
A follow-up study by ETMC was presented recently [216] and focused on understanding systematic uncertainties originating from excited states contamination. This study used a high-statistics analysis for the physical point ensemble used in [124]. Four (three) source-sink separations () were used for the unpolarized (helicity and transversity) case, corresponding to 0.75, 0.84, 0.93, and 1.12 fm (0.75, 0.93, and 1.12 fm) in physical units. All three analysis techniques described in Section 6.1, that is a single-state fit for each separation , a two-state fit, and the summation method, were used. For a reliable analysis, it is absolutely critical to keep the statistical uncertainties at the same level for all separations, and this is achieved with 4320, 8820, 9000, and 72990 measurements for the unpolarized PDFs at the four separations. For the helicity and transversity PDFs, the number of measurements is 3240, 7920, and 72990 for the separations , respectively.
A comparison of the three methods for the helicity is presented in Figure 23, where one clearly observes the discrepancy between separations 0.75 and 1.12 fm for both the real and the imaginary parts. In addition, the real part (left plot) obtained for 0.93 fm is compatible with both 0.75 and 1.12 fm. The most striking effect of excited states for this particular study can be seen in the imaginary part (right plot), where separations 0.75 and 0.93 fm are compatible, but in huge disagreement with fm, indicating that excited states are severe and one should focus on separations above 1 fm. The two-state fit is compatible with the results from the largest separation , but not with the two lower separations in the imaginary part. The summation method has large statistical uncertainties and is not providing any useful information. Based on these findings, the Authors concluded that a source-sink separation of 1.12 fm for nucleon momentum up to ~1.5 GeV is sufficient for isolating the ground state dominance within statistical uncertainties. We would like to stress the importance of having raw lattice data with similar statistical precision to avoid bias in the various analysis techniques.
We now continue the discussion with a presentation of the work of for the unpolarized distribution of [218]. The calculation was carried out using a mixed action setup of clover fermions in the valence sector on a HISQ ensemble that has lattice spacing fm, with spatial lattice extent fm and a pion mass ≈135 MeV [188]. A single step of hypercubic smearing (HYP) was employed to improve discretization effects, but also to possibly address a delicate issue: the mixed action setup of clover on HISQ is nonunitary and suffers from exceptional configurations as the quark masses approach their physical value for a fixed lattice spacing [235, 236]. As a consequence, the results would be biased in the presence of exceptional configurations. Based on the empirical evidence of [235, 236] for local operators, it is expected that, for physical value of the pion mass, the ensembles with lattice spacing above 0.09 fm could be vulnerable to exceptional configurations. However, this problem is not addressed in the work of for the quasi-PDFs, and a more concrete investigation is imperative to eliminate possible bias in the results.
In this work, the Gaussian momentum smearing [109] was employed, and the nucleon was boosted with momenta 2.2, 2.6, and 3 GeV. As pointed out by the Authors, one should be particularly cautious in the investigation of excited states contamination, which are expected to worsen with momentum boost, as the energy states come closer to each other. Thus, four variations of two-state fits were tested using source-sink separation of 0.72, 0.81, 0.90, and 1.08 fm giving compatible results. Despite the effort to employ different analysis techniques with the intention to eliminate excited states contamination, we believe that it is unlikely for this procedure to be conclusive, as the two-state fit alone does not guarantee reliability and the different variations used in the work of [218] are correlated. In addition, the success of the fits relies on having all correlators with similar accuracy; otherwise the fit is biased by the accurate data (typically at small values of the separation). Note that [218] does not report any measurements for the nucleon matrix elements. We stress that the statistical accuracy for the data should be verified from plots of the ratio on each separation that enters the fit.
The lattice data were properly renormalized using an RI-type scheme [181], as described in Section 7.3.2, and the quasi-PDFs were obtained using the “derivative” method. Finally, a matching appropriate for the choice of renormalization was applied [166, 181] to bring the final estimates in the scheme. This is an alternative to the procedure of ETMC in which a two-step process is used in order to bring the renormalized quasi-PDFs in the and then match using a proper matching formula. Both processes are equivalent to a one-loop correction, which is currently the level at which both the conversion and the matching formula are available. It is yet to be identified which process brings the final results closer to a two-loop correction; this will be possible once the two-loop expressions are extracted.
The final result for the unpolarized PDF is shown in the left plot of Figure 24, together with global fit data from CT14 [210], with agreement between the two within uncertainties. The same setup was applied for the helicity PDF presented in [219], where two values for the source-sink separations were added, giving fm. The number of measurements for each separation is 16000, 32000, 32000, 64000, 64000, and 128000, respectively; these data are used exclusively for a two-state fit, but it would certainly be critical to compare with plateau values and the summation method. The renormalization program includes various choices for the scales appearing in the RI and schemes, and we refer the Reader to [219] for details. The final estimates are given in the scheme at 3 GeV and are shown in the right plot of Figure 24 (red curve with grey band for reported systematics). The lattice data have similar behavior as the phenomenological estimates [207, 215, 217].
9.2.2. Transversity PDF
Extracting the transversity PDF is a powerful demonstration of the advances in the quasi-PDFs approach using Lattice QCD simulations. Preliminary studies can be found in the literature already in 2016 [107, 108]. However, these lack two major components that prevent comparison with global analysis fits: proper renormalization and matching procedures. Complete studies of the transversity quasi-PDFs appeared this year by ETMC [48] and by [49] using the same lattice setup as their work for the unpolarized and helicity PDFs described above.
The main motivation for first-principle calculations of the transversity PDF is the fact that it is less known experimentally [220, 237–241], because it is chirally odd, and totally inclusive processes cannot be used. In particular, one may extract information on the transversity PDF from annihilation into dihadrons with transverse momentum [242–244] and semi-inclusive deep-inelastic scattering (SIDIS) TMD data for single hadron production [245–247]. This method requires disentanglement of the dependence on the momentum fraction from the transverse momentum on TMD form factors and TMD PDFs. Alternatively, dihadron SIDIS cross-section data can be analyzed to obtain the transversity distribution directly from the measured asymmetry [248–250]. However, this analysis leads to large uncertainties, as the available data are less precise, and the collinear factorization at large is problematic [251].
The ETM Collaboration presented the first computation of the -dependence for the transversity PDF in [48] in Lattice QCD which includes a nonperturbative renormalization in lattice regularization (), and a matching procedure similar to the scheme of [124]. The latter was recalculated using the appropriate tensor nonlocal operator. We remind the Reader of the parameters for the ensemble at a pion mass of 130 MeV [234], which has the lattice spacing and the volume of . As in the case of the unpolarized and helicity PDFs, the nucleon was boosted with momenta 0.83, 1.11, and 1.38 GeV, while the source-sink separation was fixed to fm for the final results. This value has been chosen after a thorough investigation of excited states [216]. The statistics increases with the nucleon momentum, that is, 9600, 38250, and 72990 measurements for momenta 0.83, 1.11, and 1.38 GeV, respectively.
The final lattice data for the transversity isovector PDF, , are shown in Figure 25 in the scheme and at a scale of GeV, so that they can be compared to phenomenological fits extracted at the same scale. In the left plot, we show the dependence on the nucleon momentum, which is found to be small for most values of , with the highest momentum having milder oscillatory behavior. In the right panel, we present the lattice data for the highest momentum and compare with phenomenological fits on SIDIS data without [220] or with [220] constraints from lattice estimates of the tensor charge (“SIDIS+lattice”). The difference in the statistical accuracy between the global fit and the lattice data is impressive, with the data of [48] being more accurate than both the constrained and the unconstrained SIDIS results. One way to check for systematic uncertainties is to compare the tensor charge as extracted: (a) directly from the local tensor operator, and (b) by integrating over within the interval of PDFs. This consistency check reveals that both results are well compatible within uncertainties and both give a value of (the exact matching of the two numbers is to some degree accidental). Even though the agreement is nontrivial, as the steps leading to both values are different, it is, obviously, not sufficient for a complete quantitative understanding of systematic effects.
The latest work of on the quasi-PDFs was very recently extended to the transversity distribution [49], using the same ensemble with clover valence quarks on a HISQ sea [188], physical pion mass, the lattice spacing fm, and the volume of . For the lattice setup, we refer the Reader to Section 9.2.1 and [218, 219]. Six source-sink separations were used with the highest at 1.08 fm and the same statistics as in [219]. These were analyzed based on different variations of a two-state fit, and the extracted matrix elements are shown in Figure 26 for the three momenta employed in this work, that is 2.2, 2.6, and 3 GeV. It is observed that the dependence on the nucleon momentum is weak within the uncertainties, which also holds for the matched PDFs. This can be seen in Figure 3 of [49], with the exception of the very small- region. However, this is not conclusive, as lattice calculations have limitations on the reliability for this region. The observed convergence could be partly due to limitations in the matching formula, which is available to one-loop level only. Given the latter, a convergence can be possibly achieved at smaller nucleon momentum, which has the advantage that excited states can be better controlled. Evidence of nonnegligible excited states contamination for momenta as high as 3 GeV can be seen in Figure 26, particularly in the real part where the matrix element becomes negative for large values of . The latter is a clear evidence of excited states and it has been observed in other works that increasing source-sink separation (thus decreasing the contamination) brings the real part of large- bare matrix elements to values compatible with zero; see, e.g., the upper left plot of Figure 1 in [216] and, to a lesser extent, the left panel of Figure 23.
Final estimates for the transversity PDF are given in Figure 27, where the lattice results (blue curve) underestimate the global fits from LMPSS17 [220] for and are slightly higher in the region . Note that the results of ETMC shown in Figure 25 overlap with the fit from LMPSS17 (“SIDIS+lattice” in Figure 25) [220] for and overestimate it for , possibly due to the oscillatory behavior. We believe that the difference in the behavior of the data from ETMC and has its origin in the employment of the derivative method by instead of the standard Fourier transform, which, as argued in Section 6.3, may lead to uncontrolled systematic uncertainties.
10. Other Results from the Quasidistribution Approach
In the previous section, we have concentrated on numerical results for the isovector quark PDFs in the nucleon. Now, we review other results obtained with the quasidistribution method, for mesonic DAs and PDFs, as well as first exploratory results for gluon PDFs.
10.1. Meson DAs
Arguably the simplest partonic functions are distribution amplitudes (DAs) of mesons. The interest in them is at least for two reasons. First, being very simple, they can be used for investigating and comparing different techniques. Many exploratory studies were or are performed focusing on the pion DA. Second, mesonic DAs are of considerable physical interest as well. They represent probability amplitudes of finding a configuration in the final meson state, with the quark carrying fraction of the total momentum and the antiquark fraction . In phenomenology, they serve as nonperturbative inputs in analyses of hard exclusive processes with mesons, most notably the pion, in the final state. The shape of the pion DA is well known at large momentum transfers, where it follows an asymptotic form . However, for smaller momentum transfers, different models lead to different functional forms and, hence, a first-principle investigation on the lattice could shed light on this issue and eliminate the theoretical uncertainty in analyses requiring DA as an input.
The first lattice computation of the pion quasi-DA was presented early in 2017 by J.-H. Zhang et al. [123]. They used a setup of clover valence quarks on an HISQ sea with pion mass of 310 MeV, lattice spacing fm, and lattice volume that yields . The measurements were done on 986 gauge field configurations with 3 source positions, averaging over two directions of boost. The employed pion momenta were and , which corresponds to around 0.86 and 1.32 GeV, respectively. The matrix elements defining the quasi-DA can be accessed with two-point correlation functions and, after taking the Fourier transform, the distribution can be matched to its light-cone counterpart. At this stage, only matching formulae in the transverse momentum cutoff scheme were available from [115]. The Authors calculated the pion mass correction of along the lines of their earlier derivation of NMCs for nucleon quasi-PDFs [107]. They also parametrized the higher-twist corrections by extrapolating linearly in to zero after employing the matching and the mass correction. The results were presented, first, without any renormalization of the Wilson-line-related power divergence and, next, with the latter being subtracted by multiplication of the matrix elements by (“improved” pion DA), with extracted from the static potential. The latter computation was performed only on one lattice spacing and, hence, the obtained value, MeV, was attributed a large uncertainty of 200 MeV.
The final result for the improved DA, after matching and mass corrections, is shown in Figure 28. In the left panel, the curves correspond to GeV for the transverse momentum cutoff and the central value of without uncertainty (error bands correspond to statistical uncertainties). One can see significant dependence on the pion momentum and the nonphysical nonzero values outside of . In the right panel, the uncertainty in the determination of is included and dominates the total error. Within this large uncertainty, there is reasonable agreement with various models and parametrizations. However, the precision is clearly not enough to disentangle between different possibilities suggested from phenomenology. Naturally, that was not the aim of an exploratory study, where several systematic uncertainties are yet to be addressed (see Section 6.3 for a general discussion of such systematics). The main result of the paper is, thus, establishing the feasibility of the computation and the qualitative agreement with phenomenology can certainly be considered as reassuring.
The above study was extended by the Collaboration [224] to include also the kaon and mesons, with the view of studying the flavor symmetry breaking and testing predictions of chiral perturbation theory (PT). Further extension with respect to [123] was to include momentum smearing to improve the signal for the boosted meson and access one more unit of lattice momentum, i.e., , corresponding to around 1.74 GeV. The used gauge field configurations ensemble was the same as in [123]
Technically, the computation of the kaon DA amounts to changing the mass of one valence quark to represent the strange quark mass. For the meson, things are more subtle, because of the ensuing quark-disconnected diagrams and mixing with the singlet state. The Authors argued that the mixing is small and can be safely neglected, while the effect from using only connected diagrams (corresponding to the unphysical meson) can be taken into account and the final result for can be approximated as . They again used the “improved” pion DA definition but employed three additional ensembles, with fm, all at the physical pion mass, to determine precisely the mass counterterm , the dominating source of uncertainty in their previous work. The computation yielded the value -253 MeV. The final DAs show that the data at the two largest momenta are compatible with each other in most regions of , while there are also regions where the behavior is nonmonotonic in . Hence, the Authors did not attempt the extrapolation to . The data for are rather close to the ones for ; hence the result for is also close to the two.
A comparison of the pion and kaon DAs (at the largest meson boost) with models and parametrizations is shown in Figure 29. At the attained meson momenta, there are still sizable contributions outside of the physical region. Since the distributions are normalized to 1, the central regions of the LaMET DAs are significantly below all other results. The Authors concluded that larger momenta are needed, together with higher-order matching. Moreover, most of the standard lattice systematics is yet to be addressed (see Section 6.3). The Authors also converted their results on to the data for the pseudoscalar-scalar current correlator, to compare to the auxiliary light quark approach of [71], and found compatible behavior (see also Section 11.3). Finally, first attempt at testing the flavor symmetry breaking was made, with indications of agreement with PT. The effect manifests itself mostly as the difference between the DAs of and , predicted to be by PT. For a more complete study, simulations at additional light quark masses are needed.
10.2. Meson PDFs
Apart from DAs of mesons, the interest is, obviously, also in their PDFs, particularly for the pion. Phenomenological extraction of the pion PDF uses predominantly experimental data from the Drell-Yan process in the pion-nucleon scattering. This established that the large- behavior of the pion PDF is [228], corroborated by certain models. However, other models indicate rather a decay. A first-principle computation could solve this discrepancy.
The first lattice extraction of the pion PDF based on LaMET was shown in [226] by the Collaboration. They used again the same ensemble as for the pion DA (see previous subsection) and applied boosts of 0.86, 1.32, and 1.74 GeV to the pion. The (isovector) quasi-PDF is defined analogously to the nucleon case and the Dirac structure was chosen to be to avoid the mixing discovered in [106]. The Authors used four source-sink separations, ranging from to (0.72 to 1.08 fm), to investigate excited states contamination. They demonstrated that different two-state fits lead to consistent results in the real part of the matrix elements, at their intermediate pion momentum. The effects in the imaginary part were, unfortunately, not shown. As we argued in Section 6.1, the two-state method is, by itself, not enough to check excited states effects. Much stronger conclusions can be drawn from comparison of two-state fits with the plateau method. Else, the danger is that two-state fits are dominated by the lowest source-sink separations and/or many excited states mimic one excited state. Moreover, it is not clear what happens in this study at the largest pion boost, where the excited states contamination is bound to be larger.
For renormalization, followed two procedures. They used a variant of RI/MOM but also decided to apply the procedure of removing the power divergence by the mass counterterm determined from the static potential for comparison. The RI-renormalized quasi-PDF results were matched directly to the scheme using the kernel of [199] and mass corrections were applied [107]. To reduce the oscillations in the large- region, the Authors used the derivative method. They investigated the momentum dependence of the final results and for RI results they also varied the renormalization scale . Comparison between the RI and Wilson line renormalizations revealed large differences, attributed by the Authors to possibly large higher-order corrections in the matching.
The final results for the -renormalized pion PDF, taken from lattice quasi-PDFs renormalized in the RI scheme and matched to at GeV, are shown in Figure 30. The result is contrasted with a model calculation based on Dyson-Schwinger equations (DSE at a different scale of GeV) [227] and with the ASV fit to experimental Drell-Yan data [228]. Within the reported uncertainty, coming from statistical errors and comparing results for two values of the RI intermediate scale, the Authors observed compatibility with the ASV fit for small , where the ASV fit disagrees with the Dyson-Schwinger analysis. For large , the phenomenological fit agrees with DSE, but the extraction lies significantly above the two. The reliability of the computation (in particular the large- region) is expected to increase when using larger pion boosts and decreasing the pion mass towards its physical value, as well as when taking higher-order matching into account. Obviously, other systematics, such as cutoff effects and FVE, need to be addressed too; see Section 6.3.
10.3. Gluon PDFs
Very recently, the first investigation of quasi-gluon PDFs appeared [229], by Z.-Y. Fan et al. Needless to say, gluon PDFs are relevant for many analyses, especially in the small- region, where they become the dominating partons. Phenomenologically, they are determined from DIS and jet-production cross-sections. The employed lattice setup consisted of valence overlap quarks on an domain-wall sea with lattice spacing fm, lattice volume , and pion mass of 330 MeV. The Authors used two valence pion masses, one slightly larger than the sea quark mass (340 MeV) and one corresponding to light quarks having the strange quark mass (pion mass 678 MeV). The computations were performed on 203 gauge field configurations with many smeared point sources, yielding total measurements for the two-point functions. The bare matrix elements were extracted using the method proposed in [252], based on the derivative of the summed ratio of three-point and two-point functions, grounded on the Feynman-Hellmann theorem.
Fan et al. employed the following definition of gluon quasi-PDF:with the bare matrix element being the boosted proton state expectation value of the Euclidean operator:where and the gluon operator is subject to HYP smearing to improve the signal. This operator was shown not to be multiplicatively renormalizable by the Authors of [137] (see also discussion in Section 5.2.1 about the renormalizability of gluon quasi-PDFs). However, in this exploratory study, the Authors did not perform a rigorous renormalization procedure, but only tried to eliminate the power divergence by taking the ratiowith equal to . This was justified by an empirical observation from unpolarized quark quasi-PDFs, where an analogous ratio reproduces the RI-renormalized matrix elements with deviation.
In their numerical investigation, Fan et al. compared the -dependence of bare and ratio-renormalized matrix elements for different levels of HYP smearing, using nucleon momenta of 0, 0.46 and 0.92 GeV (without momentum smearing). At this level of precision, not much sensitivity to could be seen. The bare matrix elements are significantly enhanced by the removal of the power divergence. Since the lattice computation is very noisy in the gluon sector, the signal extends only to fm. The Authors also compared results from the operator to three other operators that can be used to define gluon quasi-PDFs, finding that the other ones either suffer from large mixing with higher-twist operators or provide a worse signal. They also plotted the results for the ratio-renormalized matrix elements together with two phenomenological gluon PDFs inverse Fourier transformed to coordinate space, observing compatibility within large uncertainties for their smaller valence pion mass; see Figure 31. Finally, matrix elements were shown also for gluon quasi-PDF in the pion.
The Authors concluded that, at the present level of precision, their study could not constrain gluon PDFs, which would require taking the Fourier transform and performing the matching to the light-cone PDF. Due to the fact that the magnitude of the gluon PDF is significant predominantly for small , the distribution in coordinate space is very broad, necessitating reaching large values of (while in the current study only could be reached). Thus, significant improvements are needed to obtain a reliable gluon PDF from the quasidistribution approach. The challenge is further extended by the mixing between the gluon quasi-PDF and the singlet quark quasi-PDFs (see Section 8.2), which have not been yet explored on the lattice and would require calculations involving quark-disconnected diagrams.
11. Results from Other Approaches
The last two sections were devoted to reviewing results obtained for the -dependence of nonsinglet quark PDFs, gluon PDFs, and meson DAs/PDFs from the quasidistribution method. In the present one, we discuss some other results obtained in the last few years from alternative approaches, shortly described in Section 2. We review them in the order of discussion in Section 2.
11.1. Hadronic Tensor
Despite being proposed in the early 1990s, the hadronic tensor approach [54–56] (see also Section 2.1) has not led to many numerical applications, because it requires the computation of difficult four-point correlators and faces the inverse Laplace transform problem. However, recently there is renewed interest in it, due to hugely increased computational powers and new reconstruction techniques to tackle the inverse problem. In [63], J. Liang, K.-F. Liu, and Y.-B. Yang presented preliminary results obtained using the classical Backus-Gilbert technique [253]. They used an ensemble of clover fermions on an anisotropic lattice with pion mass 640 MeV and lattice spacing of 0.1785 fm, performing measurements on 500 gauge field configurations.
The preliminary results are shown in Figure 32. The Euclidean hadronic tensor (left plot) vs. the current separation is shown for nucleon at rest () with momentum transfer and corresponds to connected sea antiup and antidown partons. The reconstructed Minkowski tensor , where is conjugate to in the inverse Laplace transform, is shown in the right plot. The first peaks are elastic and correspond to the energy transfer invoked by the momentum transfer. The less pronounced second peaks are quasielastic and are related to nucleon excitations. Unfortunately, with these kinematics, the DIS region is inaccessible, as it would require both ( GeV in this case) and at the same time much larger than the one corresponding to the quasielastic peaks (extending to , which yields GeV). This could be achieved on lattices with much smaller lattice spacings. The Authors, nevertheless, concluded that the observation of both elastic and quasielastic peaks is encouraging.
The investigations are continued and further results were presented in the Lattice 2018 Symposium, using other reconstruction methods and an ensemble with much finer lattice spacing, fm (in the temporal direction), lattice size , and lower pion mass of 380 MeV; see upcoming proceedings [254] for more details.
11.2. Auxiliary Heavy Quark
The approach with auxiliary heavy quark [58] (see also Section 2.3) was also recently revived by its Authors, W. Detmold and C.-J. D. Lin, in collaboration with I. Kanamori, S. Mondal, and Y. Zhao [68]. Their study is aimed at extracting the pion DA and the current investigations employed three quenched ensembles (Wilson plaquette action discretization), with lattice spacings of 0.05 fm, 0.06 fm, and 0.075 fm and fixed physical spatial extent of fm, . The valence pion mass is 450 MeV, and the auxiliary heavy quark mass 1.3 or 2 GeV.
The calculation proceeds via evaluating the vacuum-to-pion matrix elements of the product of two heavy-light currents separated in spacetime. The spatial Fourier transform of such matrix elements, for large enough temporal separation of the three points in the correlator, gives a quantity called , where is the pion momentum, the momentum transfer, and the separation of currents. is then an input to a temporal Fourier transform yielding the Euclidean hadronic tensor , which, in the continuum limit, gives access to moments of the structure function by varying . As an illustration, the integrand of this Fourier transform is shown in Figure 33 (left), for , pion at rest, and with minimal spatial momentum transfer of in the -direction. The heavy quark mass is 1.3 GeV and two lattice spacings and two values of are shown. The signal is clear, but lattice cutoff effects are not negligible, as also evidenced in the right plot of Figure 33, showing the full quantity for three lattice spacings and three choices of . Since the extraction of moments requires reliable extrapolation to the continuum limit, the Authors prefer first to analyze smaller lattice spacings. To this end, they already have quenched ensembles with lattice spacings down to 0.025 fm. Furthermore, the momentum smearing technique will be employed to enhance the signal for a moving pion. Preliminary investigation of this case was also shown in [68].
11.3. Auxiliary Light Quark
Instead of an auxiliary heavy quark, one can also use an auxiliary light quark [69] (see also Section 2.4). The wave of renewed interest in light-cone distribution functions on the lattice in recent years sparked also revival of numerical studies of this approach, by the Regensburg group [70, 71]. Their aim is to extract the pion DA. In their exploratory study, they employed one gauge field ensemble of clover fermions, with lattice spacing fm, lattice volume , and pion mass 295 MeV. The auxiliary light quark has the same mass as the physical quarks. Relatively large momenta were reached, up to around 2 GeV, thanks to the momentum smearing technique introduced by the same group. It is clear that going much beyond 2 GeV is currently impossible on the lattice, if aiming at a reliable analysis, in particular large enough temporal separations between points in the three-point correlator.
As in the auxiliary heavy quark approach, the lattice part consists in calculating the vacuum-to-pion matrix element of two currents, separated spatially by . In [71], the pion DA was extracted from the scalar-pseudoscalar channel. The Authors paid particular attention to discretization effects from the breaking of rotational invariance that leads to very different behavior of points with the same , but different choices of its components. In particular, the “democratic” points, like (1,1,1), tend to behave better than “nondemocratic” ones, e.g. (1,0,0). This is a well-known effect in coordinate space and it can be seen already in the free theory (cf., e.g., [255]). To improve the behavior, one can discard points that are too “nondemocratic” and also define a tree-level improvement coefficient. Renormalization (involving only local operators) was performed in the RI/MOM scheme, with a three-loop conversion to the scheme. The data at different renormalization scales and different Ioffe times were compared to continuum perturbation theory predictions for three different phenomenological models, at leading twist and with twist-4 corrections. The Authors concluded that there are indications of deviating from the asymptotic form of the pion DA, , in the large Ioffe time region; however, for reliable conclusions one needs to access this region at larger pion boosts, to keep in the perturbative region. Larger pion boosts should be accompanied by computations at smaller lattice spacings, to keep the momenta sufficiently away from the cutoff. At small Ioffe times, one would need significantly larger statistics to disentangle between the three models.
The follow-up work of [70], by the same group and using the same lattice ensemble, concentrated on exploring higher-twist effects (HTE) and comparing results from six channels: vector-vector (VV), axial-axial (AA), vector-axial (VA), axial-vector (AV), scalar-pseudoscalar (SP), and pseudoscalar-scalar (PS). Other channels, like scalar-vector, although possible in principle, may suffer from enhanced HTE. For the employed channels, the Authors calculated the leading HTE in the framework of three phenomenological models. Results from some channels can be combined to eliminate certain effects; e.g., imaginary parts cancel in SP+PS. In the end, three linear combinations were formed: VV+AA, VA+AV, and SP+PS. On the lattice side, the Regensburg group also tested another technique to calculate the all-to-all propagator, using stochastic estimators instead of the sequential source method. This technique allowed them to take a volume average at a smaller computational cost and hence it was considered superior to the previously employed one. They chose six momentum vector choices at five different boost magnitudes and to avoid enhanced lattice artifacts observed at very small distances.
Example results for the Ioffe-time dependence of the pion DA are shown in Figure 34 (left). They correspond to two of the linear combinations, VV+AA and SP+PS, and one spatial distance of fm. The lattice data are compared to tree-level and one-loop-corrected continuum perturbative results, with and without leading HTE. The Authors observed that the sign and magnitude of the predicted splitting are in good agreement with the data, but quantitative differences emerge. Obviously, no quantitative agreement was expected, since lattice data have their systematics and the phenomenological models may not be correct and/or are subject to unknown higher-order corrections. To investigate the final shape of DA, Bali et al. performed a global fit to all channels and all data at different separations and momenta, using three different parametrizations of the leading-twist DA and different fitting ranges. An example result (with only statistical errors), for two parametrizations and one selected fitting range, is shown in Figure 34 (right). Both DAs describe the lattice data equally well, having similar second Gegenbauer coefficients , which is the only relevant parameter for the description of available data. With data extending to larger Ioffe times, the next Gegenbauer coefficient should become accessible and allow us to disentangle between the two parametrizations. The Authors concluded that these results are very promising and the dominating uncertainty is the systematic one, which can be reliably improved, in particular, by using smaller lattice spacings, larger pion boosts, and higher-order perturbative corrections and HTE.
11.4. Pseudodistributions
The first numerical investigation of the pseudodistribution approach [75–77] (see also Section 2.6) was performed by J. Karpie, K. Orginos, A. Radyushkin, and S. Zafeiropoulos in 2017. The computation proceeded using a quenched ensemble with lattice spacing fm, lattice volume , and clover fermions in the valence sector, with pion mass around 600 MeV. The employed momenta for the nucleon boost reached up to , i.e., approximately 2.5 GeV. The matrix elements (lattice ITDs) were obtained using the methodology of [252]. From these, reduced matrix elements, , were formed and they require no further renormalization. After plotting vs. the Ioffe time, the Authors noticed a significant -dependence of the results and applied the one-loop LLA evolution for all points with , i.e., MeV. When using and evolving to , this led to all points collapsing close to a universal line, for both the real part and the imaginary part. Clearly, it is difficult to imagine one-loop perturbative formula to work rigorously at scales down to 500 MeV. Hence, the LLA evolution should rather be treated as a model of evolution. The model was further extended to check the behavior of data under LLA for even lower scales . Around , the evolution was observed to stop. Hence, points for were treated as if they corresponded to the scale . The result of this procedure is shown in the left panel of Figure 35. The evolved data were fitted to cosine Fourier transforms of -type functions ( – normalization, , ), which yielded the blue band in the plot. The corresponding PDFs at two scales are shown in Figure 35 (right) and compared to three sets of phenomenological PDFs. Obviously, no quantitative agreement was expected, but the general shape of the ensuing PDF evinces features of the experimental distributions and the evolution from the original scale of GeV to 2 GeV moves the lattice-extracted PDFs closer to phenomenology.
As argued by Radyushkin in [80, 81] (see Section 8.2), the LLA is only an approximation appropriate for studying the dependence. To obtain the full PDF, one should perform the matching procedure based on factorization [81–83, 85, 86], taking into account all one-loop corrections. The matching equation (93) has the outcome of effectively changing the relation between the lattice scale and the scale, as discussed in Section 8.2. Radyushkin [81] applied the matching to the data of [87] and found that the matched ITD, denoted by , is approximately equal to the reduced ITD and thus the rescaling factor is close to 4, as opposed to the LLA value of about 1.12. The matched ITD is shown in the left plot of Figure 36 and the resulting PDF is in its right panel. Both plots also contain, for comparison, data from phenomenological parametrizations, inverse Fourier transformed for the ITD plot. As could be expected, the matched ITDs lie close to a universal curve and the curve corresponds to a fit to the same model as in [87], with parameters and . The fitted curve lies significantly below the phenomenological ITD. Correspondingly, the final PDF deviates from phenomenology, especially for small and intermediate . The Author pointed out that alternative fitting ansatzes lead to a similar curve as in the left panel, but the final PDF may significantly differ. The reason for this is that the ITD is unknown in the whole region and, having a limited set of Ioffe times, one needs to add assumptions about the behavior of the ITD outside the region or about the functional form of the PDF. Radyushkin also compared the present result to the one from LLA in [87]. The final PDF is changed to a large extent and is further away from phenomenology. He pointed out that this is because the LLA analysis assumes that the final scale differs from by only the factor 1.12, while the full one-loop formula implies that the true scale is in fact around , i.e., about 4 GeV. Thus, the evolution to the reference scale of phenomenological PDFs, 2 GeV, should proceed downwards from 4 GeV to 2 GeV and not upwards from 1 GeV to 2 GeV.
The final result that we report from the pseudodistribution approach is the computation of the two lowest moments of the isovector unpolarized PDF, erroneously claimed to be impossible due to fatal flaws in the approach in [164]. We refer to Section 6.3 for more details about this argument and its refutation. In [90], J. Karpie, K. Orginos, and S. Zafeiropoulos used the same quenched ensemble as in [87] and demonstrated that the two lowest moments agree with an earlier explicit computations thereof by the QCDSF collaboration [230]; see Figure 37.
Further progress was reported in the Lattice 2018 Symposium, including first calculations with dynamical flavors [256] and the issue of reconstruction of distributions from a limited set of data 11.
11.5. OPE without OPE
The approach dubbed “OPE without OPE” was first investigated numerically in [91] (see also Section 2.7) by the QCDSF collaboration. The Authors took an exemplary parametrization of a nonsinglet PDF and applied the proposed method. They showed that the parametrized PDF can be reconstructed from computed moments with very promising agreement already using a very limited set of data points; see Figure 38 (left). Moreover, they performed an exploratory study with real lattice data, employing an ensemble of clover fermions, with lattice spacing fm and lattice volume . They computed the Compton amplitude for 10 spatial momenta and one momentum transfer . The result is shown in the right plot of Figure 38. For low momenta, the precision was found to be already very good and for larger ones the usage of the momentum smearing technique is planned. Further exploration is in progress, at three lattice spacings and a pion mass of 470 MeV, and results were reported in the Lattice 2018 Symposium; see upcoming proceedings [257].
11.6. Good Lattice Cross-Sections
This approach, suggested in [85, 86, 95] (see Section 2.8) and closely related to the auxiliary light quark method, is being pursued by the theory group at JLab, aiming at meson PDFs [231]. They use clover fermions with lattice spacing fm, pion mass of 430 MeV, and the largest momentum employed is about 1.5 GeV. Preliminary results are illustrated in Figure 3912. It shows the vector-vector () current-current matrix element for the pion PDF calculation vs. the Ioffe time , where is the pion boost and the separation of currents. Different colors correspond to different separations (in lattice units). The higher-twist effects are visible at large separations and the Authors are calculating the NLO perturbative kernel that will give a correction in . For more results, see [231].
12. Summary and Future Prospects
In this paper, we give an overview of several approaches to obtain the Bjorken- dependence of partonic distribution functions from ab initio calculations in Lattice QCD. A major part of this review is dedicated to a discussion on the state-of-the-art of the field, demonstrated with modern numerical simulations. We considered different theoretical ideas that were proposed over the last years to access parton distribution functions (PDFs) and parton distribution amplitudes (DAs), as well as more complex generalized parton distributions (GPDs) and transverse momentum dependent PDFs (TMDs). Even though their -dependence was believed to be practically impossible to calculate on the lattice, breakthrough ideas were conceived and sparked renewed interest in these difficult observables. Arguably, the single most seminal idea was the one of X. Ji, who developed a general framework for accessing light-cone quantities on a Euclidean lattice, the quasidistribution approach. This framework itself has been heavily studied and has led to very encouraging results, but, moreover, it has prompted also the rediscovery of previously proposed ideas, like the hadronic tensor, and approaches with auxiliary heavy/light quarks. It has spawned also new or related concepts, such as pseudodistributions, OPE without OPE, and good lattice cross-sections.
As a summary, we would like to offer the Reader a flowchart (Figure 40) with an overview of how progress of the different approaches has been evolving. For all these new methods, we distinguish four general stages in the evolution of our understanding.(i)Starting with the proposed theoretical idea (e.g., quasidistributions, good lattice cross-sections, pseudodistributions, etc.), several challenges (theoretical and technical) must be studied and overcome to achieve a successful implementation of the method. Theoretical analyses of the idea may lead to additional challenges on the lattice.(ii)The second stage is exploratory studies aiming at a demonstration of the feasibility of the method. During this stage, further technical difficulties can be revealed, as well as possible additional theoretical challenges.(iii)The next stage consists of more advanced studies focusing on a more thorough investigation of the method and first estimation of certain systematic effects. Before precision calculations can be carried out with full systematics taken into account, usually further technical difficulties must be overcome. During this evolution of knowledge, additional theoretical challenges may arise, as well as subleading systematic uncertainties.(iv)The final desired outcome is an accurate and reliable Lattice QCD estimate of the observable of interest. For this to be achieved, the various sources of uncertainties must be quantified and brought under control.
Based on Figure 40, we comment on the status of the different approaches presented in this paper. Most of the methods are still at an exploratory stage, or toward the third phase of advanced studies. Notable exceptions are, in our view, the isovector quark quasi-PDFs, as the numerical exploration began immediately after Ji’s proposal. As we argued, the exploratory studies of 2014-2016 (see Section 3.2) showed the feasibility of the method and identified theoretical and lattice challenges. Among the former, we discussed the role of the spacetime signature and renormalizability (Section 5), renormalization studies (Section 7), and matching onto light-cone PDFs (Sections 3.1 and 8). The lattice challenges were of various origins and we described them in detail (Section 6). The most recent results of 2018 are, undoubtedly, in the advanced stage, using ensembles at physical pion masses and optimized lattice techniques, as well as reliable renormalization and matching procedures (Section 9). However, reaching into the precision era is still extremely demanding and will require overcoming further challenges, most of them classifiable as lattice ones. Careful investigation of systematic uncertainties is imperative and this will necessitate additional simulations employing ensembles with finer lattice spacings, larger volumes, accessing larger nucleon boosts, etc., as thoroughly reviewed in Section 6.3. This will require tremendous amount of computing time but is, in principle, possible. The difficult part of this programme is to reliably access large nucleon momenta and the main obstacle is the exponential signal-to-noise problem when increasing the boost and, at the same time, increasing the source-sink separation to avoid excited states contamination. We have highlighted the latter, since, in our view, this is an essential feature, if quasi-PDFs are to give reliable results. The present results are highly encouraging and steady increase of convergence towards phenomenologically extracted PDFs is being observed, even with partial agreement within uncertainties in some Bjorken- regions. However, fully reliable results are still to be obtained. Nevertheless, it is highly conceivable that these lattice-extracted results may have extensive phenomenological impact, in particular the transversity PDF, which is much less constrained experimentally.
The quasidistribution approach has also been applied to other kinds of distributions (besides the isovector flavor combination) and notable progress has recently been achieved. We discussed the exploratory studies concerning quark DAs/PDFs for mesons and gluonic PDFs (Section 10). These results are promising for prospective reliable calculations that will also have an impact on phenomenological studies. However, as Figure 40 indicates, there are already challenges to go to the advanced stage, especially in the gluonic sector, which is characterized by noisy signal and mixings under matching with singlet quark PDFs, the latter requiring computation of noisy quark-disconnected diagrams. Yet other quasidistributions that are accessible, in principle, are quasi-GPDs and quasi-TMDs (Section 8.2). These are, obviously, much more difficult to compute, given the fact that they involve additional variables such as momentum transfer or transverse momentum. Both are receiving considerable theoretical attention and continuous progress, but numerical explorations are still absent and, in the case of quasi-TMDs, important theoretical challenges are yet to be overcome.
Even though quasidistributions are currently the most explored, other approaches are beginning to yield very interesting results as well. Several exploratory studies have been reported for quark PDFs and DAs of nucleons and pions (Section 11). These methods are in different phases of exploratory studies, but steadily pushing towards more advanced investigations. Theoretical and lattice challenges are beginning to be clear. We note that many of them are common to all approaches, such as cutoff effects, other typical lattice systematics, or the need for precise signal extraction for highly boosted hadrons. However, some of them are more specific to certain approaches, such as the renormalization of nonlocal operators for quasidistributions. The level of numerical difficulty may also vary. For example, some approaches require the computation of three-/two-point functions for PDFs/DAs (e.g., quasidistributions), while other ones necessitate the use of four-/three-point correlators (e.g., hadronic tensor, auxiliary quark methods). It is also clear that all these approaches, even though aiming at the same physical observables, may have very different systematics in practice. Hence, it can be expected that a global fitting strategy, combining results from various methods, can prove in the end to be the optimal one. Thus, all the efforts of the lattice community, with the aid of experts in phenomenology, can contribute to obtaining reliable first-principle determinations of partonic distributions.
Appendix
Abbreviations
| 1PI: | 1-Particle irreducible |
| APE: | Array processor experiment |
| CAA: | Covariant Approximation Averaging |
| PT: | Chiral perturbation theory |
| DA: | Distribution amplitude |
| DIS: | Deep-inelastic scattering |
| DR: | Dimensional regularization |
| DSM: | Diquark spectator model |
| DVCS: | Deeply Virtual Compton Scattering |
| DVMP: | Deeply Virtual Meson Production |
| EIC: | Electron-Ion Collider |
| ETMC: | Extended Twisted Mass Collaboration13 |
| FVE: | Finite volume effects |
| GPD: | Generalized parton distribution |
| HCS: | Hadronic cross-section |
| HISQ: | Highly improved staggered quarks |
| HP: | High precision |
| HQET: | Heavy Quark Effective Theory |
| HYP: | Hypercubic |
| HTE: | Higher-twist effects |
| IMF: | Infinite momentum frame |
| IR: | Infrared |
| ITD: | Ioffe-time distribution |
| JLab: | Jefferson Laboratory |
| LaMET: | Large Momentum Effective Theory |
| LCS: | Lattice cross-section |
| LCWF: | Light-cone wave function |
| LLA: | Leading logarithmic approximation |
| LP: | Low precision |
| : | Lattice Parton Physics Project |
| LR: | Lattice regularization |
| NJL: | Nambu-Jona-Lasinio |
| NLO: | Next-to-leading order |
| NMC: | Nucleon mass correction |
| NRQCD: | Nonrelativistic Quantum Chromodynamics |
| OPE: | Operator product expansion |
| PDF: | Parton distribution function |
| RI: | Regularization independent |
| RI/MOM: | Regularization-independent momentum subtraction |
| rms: | Root mean square |
| QCD: | Quantum Chromodynamics |
| QED: | Quantum Electrodynamics |
| SIDIS: | Semi-inclusive deep-inelastic scattering |
| SQM: | Spectral quark model |
| TMC: | Target mass correction |
| TMD: | Transverse momentum dependent parton distribution function |
| UV: | Ultraviolet |
| VDF: | Virtuality distribution function. |
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
We first want to thank the editors of the special issue “Transverse Momentum Dependent Observables from Low to High Energy: Factorization, Evolution, and Global Analyses” for the invitation to prepare this review and the guidance they provided throughout the process. We are also grateful to all Authors for giving permission for the figures we used to illustrate the progress of the field. We are indebted to several people with whom we had discussions over the years and who helped to shape our view on different aspects discussed in this review. Their names are, in alphabetical order, C. Alexandrou, G. Bali, R. Briceño, W. Broniowski, J.-W. Chen, I. C. Cloet, W. Detmold, V. Drach, M. Engelhardt, L. Gamberg, E. García-Ramos, K. Golec-Biernat, J. Green, K. Hadjiyiannakou, K. Jansen, X. Ji, P. Korcyl, G. Koutsou, P. Kotko, K. Kutak, C.-J.D. Lin, K.-F. Liu, S. Liuti, W. Melnitchouk, A. Metz, Z.-E. Meziani, C. Monahan, K. Orginos, H. Panagopoulos, A. Prokudin, J.-W. Qiu, A. Radyushkin, G. C. Rossi, N. Sato, M. Savage, A. Scapellato, R. Sommer, F. Steffens, I. W. Stewart, R. Sufian, J. Wagner, Ch. Wiese, J. Wosiek, Y.-B. Yang, F. Yuan, S. Zafeiropoulos, J.-H. Zhang, and Y. Zhao. We also thank all members of the TMD Topical Collaboration for enlightening discussions. Krzysztof Cichy is supported by the National Science Centre (Poland) Grant SONATA BIS no. 2016/22/E/ST2/00013. Martha Constantinou acknowledges financial support by the U.S. Department of Energy, Office of Nuclear Physics, within the framework of the TMD Topical Collaboration, as well as by the National Science Foundation under Grant no. PHY-1714407.