Systematic uncertainties in parton distribution functions from lattice QCD simulations at the physical point

TL;DR

This paper investigates the systematic uncertainties in lattice QCD calculations of nucleon parton distribution functions using the quasi-PDF approach at the physical point, providing detailed analysis and methodology.

Contribution

It offers a comprehensive analysis of systematic effects in lattice QCD quasi-PDF computations, including renormalization, excited states, and matching procedures.

Findings

Quantified systematic uncertainties in quasi-PDF calculations.

Analyzed effects of excited states and renormalization schemes.

Provided guidelines for improving lattice QCD PDF extractions.

Abstract

We present a detailed study of the helicity-dependent and helicity-independent collinear parton distribution functions (PDFs) of the nucleon, using the quasi-PDF approach. The lattice QCD computation is performed employing twisted mass fermions with a physical value of the light quark mass. We give a systematic and in-depth account of the salient features entering in the evaluation of quasi-PDFs and their relation to the light-cone PDFs. In particular, we give details for the computation of the matrix elements, including the study of the various sources of systematic uncertainties, such as excited states contamination. In addition, we discuss the non-perturbative renormalization scheme used here and its systematics, effects of truncating the Fourier transform and different matching prescriptions.

Tables

Table 1. Table 1: Parameters of the ensemble used in this work. The nucleon mass ( m N ) subscript 𝑚 𝑁 (m_{N}) , the pion mass ( m π ) subscript 𝑚 𝜋 (m_{\pi}) and the lattice spacing ( a ) 𝑎 (a) have been determined in Ref.  Alexandrou:2017xwd .

$\beta =2.10$, $c_{\mathrm{SW}}=1.57751$, $a=0.0938\left(3\right)\left(2\right)$ fm, $r_0/a=5.32\left(5\right)$
${48}^3\times 96$, $L\approx 4.5$ fm	$a\mu =0.0009$
	$m_{\pi }=0.1304\left(4\right)$ GeV
	$m_{\pi }L=2.98$
	$m_N=0.932\left(4\right)$ GeV

Table 2. Table 2: The statistics of our calculation at source-sink separation t s = 12 ​ a subscript 𝑡 𝑠 12 𝑎 t_{s}{=}12a , for each current insertion and each momentum. N conf subscript 𝑁 conf N_{\rm conf} is the number of gauge configurations, N HP subscript 𝑁 HP N_{\rm HP} ( N LP subscript 𝑁 LP N_{\rm LP} ) the number of high (low) precision measurements, and N meas subscript 𝑁 meas N_{\rm meas} the total number of measurements. A factor of 6 is included in the total measurements due to the averaging of data with the Wilson line and momentum boost in the ± x plus-or-minus 𝑥 \pm x , ± y plus-or-minus 𝑦 \pm y , ± z plus-or-minus 𝑧 \pm z -directions. An additional factor of 2 compared to all other cases is included for the transversity at P 3 = 6 ​ π L subscript 𝑃 3 6 𝜋 𝐿 P_{3}=\frac{6\pi}{L} to take into account the averages over the two possible Dirac structures, as explained in text.

	$P_3=\frac{6\pi }{L}$					$P_3=\frac{8\pi }{L}$					$P_3=\frac{10\pi }{L}$
Ins.	$N_{\mathrm{conf}}$	$N_{\mathrm{HP}}$	$N_{\mathrm{meas}}$	Ins.	$N_{\mathrm{conf}}$	$N_{\mathrm{HP}}$	$N_{\mathrm{LP}}$	$N_{\mathrm{meas}}$	Ins.	$N_{\mathrm{conf}}$	$N_{\mathrm{HP}}$	$N_{\mathrm{LP}}$	$N_{\mathrm{meas}}$
${\gamma }^3$	100	16	9600	${\gamma }^3$	425	1	16	38250	${\gamma }^3$	811	1	16	72990
${\gamma }^0$	50	16	4800	${\gamma }^0$	425	1	16	38250	${\gamma }^0$	811	1	16	72990
${\gamma }^5{\gamma }^3$	65	16	6240	${\gamma }^5{\gamma }^3$	425	1	16	38250	${\gamma }^5{\gamma }^3$	811	1	16	72990
${\sigma }^{3j}$	50	16	9600	${\sigma }^{3j}$	425	1	16	38250	${\sigma }^{3j}$	811	1	16	72990

Table 3. Table 3: Numbers of measurements of the matrix element for the unpolarized operator ( γ 0 superscript 𝛾 0 \gamma^{0} ) at each source-sink separation and P 3 = 10 ​ π / L subscript 𝑃 3 10 𝜋 𝐿 P_{3}{=}10\pi/L .

$P_3=10\pi /L$	$t_s=8a$	$t_s=9a$	$t_s=10a$	$t_s=12a$
$N_{\mathrm{conf}}$	48	98	100	811
$N_{\mathrm{meas}}$	4320	8820	9000	72990

Table 4. Table 4: Statistics used in the study of excited states for the polarized cases ( γ 5 ​ γ 3 superscript 𝛾 5 superscript 𝛾 3 \gamma^{5}\gamma^{3} : helicity, σ 3 ​ j superscript 𝜎 3 𝑗 \sigma^{3j} : transversity) at P 3 = 10 ​ π / L subscript 𝑃 3 10 𝜋 𝐿 P_{3}{=}10\pi/L .

$P_3=10\pi /L$	$t_s=8a$	$t_s=9a$	$t_s=10a$	$t_s=12a$
$N_{\mathrm{conf}}$	36	50	88	811
$N_{\mathrm{meas}}$	3240	4500	7920	72990

Table 5. Table 5: Parameters of the ensembles used to compute the Z-factors extracted in this work. For details, see Ref.  Alexandrou:2017xwd .

$\beta =2.10$,	$c_{\mathrm{SW}}=1.57751$,	$a=0.0938\left(3\right)\left(2\right)$ fm
${48}^3\times 96$	$a\mu =0.0009$	$m_{\pi }=130$ MeV
${24}^3\times 48$	$a\mu =0.003$	$m_{\pi }=235$ MeV
${24}^3\times 48$	$a\mu =0.006$	$m_{\pi }=340$ MeV

Table 6. Table 6: Numerical values for the coefficients e 𝒪 ( i ) subscript superscript 𝑒 𝑖 𝒪 e^{(i)}_{\cal O} appearing in the conversion function from MS ¯ ¯ MS \overline{\mbox{\rm MS}} to M ​ MS ¯ M ¯ MS {\rm M}\overline{\mbox{\rm MS}} in Eq. ( 26 ).

$𝒪$	$e_{𝒪}^{\left(1\right)}$	$e_{𝒪}^{\left(2\right)}$	$e_{𝒪}^{\left(3\right)}$	$e_{𝒪}^{\left(4\right)}$
$V$	$-$5	$-$3	+3	$-$3/2
$A$	$-$7	$-$3	+3	$-$3/2
$T$	$-$4	$-$4	+4	$-$4/2

Equations

$$q(x)={\int }_{-\infty }^{+\infty }d{\xi }^{-}{e}^{-ix{P}^{+}{\xi }^{-}}⟨N|\overline{\psi }({\xi }^{-})\Gamma W({\xi }^{-},0)\psi (0)|N⟩,$$

$$\tilde{q}(x,P_3)={\int }_{-z_{\max }}^{+z_{\max }}\frac{dz}{4\pi }{e}^{-ix{P}_{3}z}⟨N|\overline{\psi }(0,z)\Gamma W(z,0)\psi (0,0)|N⟩,$$

$$q(x,\mu )={\int }_{-\infty }^{\infty }\frac{d\xi }{|\xi |}C\left(\xi ,\frac{\mu }{P_3}\right)\tilde{q}\left(\frac{x}{\xi },\mu ,P_3\right)+O\left(\frac{m_N^2}{P_3^2},\frac{{\Lambda }_{QCD}^2}{P_3^2}\right),$$

$$S_F\left[\chi ,\overline{\chi },U\right]=a^4\sum _x\overline{\chi }(x)\left(D_W[U]+i{\mu }_l{\gamma }^5{\tau }^3-\frac{1}{4}{c}_{\mathrm{SW}}{\sigma }^{\mu \nu }{\mathcal{F}}^{\mu \nu }[U]\right)\chi (x),$$

$$\psi (x)\equiv e^{i\frac{\alpha }{2}{\gamma }^5{\tau }^3}\chi (x),\overline{\psi }(x)\equiv \overline{\chi }(x)e^{i\frac{\alpha }{2}{\gamma }^5{\tau }^3},$$

$$h_{\Gamma }(P_3,z)=⟨N|\overline{\psi }(0,z)\Gamma W(z,0)\psi (0,0)|N⟩,$$

$$C^{2pt}(\mathbf{P},t,0)$$

$$C^{3pt}(\mathbf{P};t_s,\tau ,0)$$

$$\mathcal{O}(\mathbf{y},\tau ;\mathbf{z})=\overline{\psi }(\mathbf{y}+\mathbf{z},\tau )\Gamma W(\mathbf{y}+z,\mathbf{y})\psi (\mathbf{y},\tau ),$$

$$h_{\Gamma }(p,z)\stackrel{0\ll \tau \ll t_s}{=}\mathcal{K}\frac{⟨C^{3pt}(p;t_s,\tau ,0){⟩}_G}{⟨C^{2pt}(p,t_s,0){⟩}_G},$$

$${\mathcal{S}}_{mom}=\frac{1}{1+6{\alpha }_G}\left(\psi (x)+{\alpha }_G\sum _j{U}_j(x)e^{i\xi \mathbf{P}\cdot \mathbf{j}}\psi (x+\hat{j})\right),$$

$$C^{imp}=\frac{1}{N_{LP}}\sum _{i=1}^{N_{LP}}{C}_{i,LP}+\frac{1}{N_{HP}}\sum _{i=1}^{N_{HP}}(C_{i,HP}-C_{i,LP}),$$

$$\mathrm{error}\propto \sqrt{2(1-r_c)+\frac{N_{HP}}{N_{LP}}},$$

$$r_{LP}\equiv \frac{|\mathrm{residue}{|}_{LP}}{|\mathrm{source}|}=2\cdot 10^{-3}$$

$$\frac{C^{3pt}(P_3;t_s,\tau )}{C^{2pt}(P_3;t_s)}=\frac{{\sum }_{n,n^{\prime }}⟨N|n^{\prime }⟩⟨n|N⟩⟨n^{\prime }|\mathcal{O}|n⟩e^{-E_{n^{\prime }}(t_s-\tau )}{e}^{-E_n\tau }}{{\sum }_n|⟨N|n⟩{|}^2{e}^{-E_n{t}_s}},$$

$$\frac{C^{3pt}(P_3;t_s,\tau )}{C^{2pt}(P_3;t_s)}=\frac{⟨0|\mathcal{O}|0⟩+f_{10}⟨1|\mathcal{O}|0⟩e^{-\Delta E(t_s-\tau )}+f_{10}^{†}⟨0|\mathcal{O}|1⟩e^{-\Delta E\tau }+|f_{10}{|}^2⟨1|\mathcal{O}|1⟩e^{-\Delta E{t}_s}+\dots }{1+|f_{10}{|}^2{e}^{-\Delta E{t}_s}+\dots },$$

$$\mathcal{R}(P_3;t_s)\equiv \sum _{\tau =a}^{t_s-a}\frac{C^{3pt}(P_3;t_s,\tau )}{C^{2pt}(P_3;t_s)}=C+\mathcal{M}{t}_s+\mathcal{O}\left(e^{-\Delta E{t}_s}\right),$$

$$C^{2pt}(P_3;t)=|A_0{|}^2{e}^{-E_{0}t}\left(1+|f_{10}{|}^2{e}^{-\Delta Et}\right),$$

$$C^{\rm 3pt}(P_{3};t_{s},\tau)=|A_{0}|^{2}e^{-E_{0}t_{s}}\Big{(}\langle 0|\mathcal{O}|0\rangle+f_{10}\langle 1|\mathcal{O}|0\rangle e^{-\Delta E\,(t_{s}-\tau)}+f_{10}^{\dagger}\langle 0|\mathcal{O}|1\rangle e^{-\Delta E\,\tau}+|f_{10}|^{2}\langle 1|\mathcal{O}|1\rangle e^{-\Delta E\,t_{s}}\Big{)}.$$

$$\displaystyle\frac{Z^{{\rm RI}^{\prime}}_{\cal O}(z,\mu_{0},m_{\pi})}{Z^{{\rm RI}^{\prime}}_{q}(\mu_{0},m_{\pi})}\frac{1}{12}{\rm Tr}\left[{\cal V}(z,p,m_{\pi})\left({\cal V}^{\rm Born}(z,p)\right)^{-1}\right]\Bigr{|}_{p^{2}{=}\mu_{0}^{2}}{=}1\,,$$

$$\displaystyle Z^{{\rm RI}^{\prime}}_{q}(\mu_{0},m_{\pi})\frac{1}{12}{\rm Tr}\left[(S(p,m_{\pi}))^{-1}\,S^{\rm Born}(p)\right]\Bigr{|}_{p^{2}=\mu_{0}^{2}}=1\,.\hskip 39.83368pt$$

$$Z_{\mathcal{O}}^{{\mathrm{RI}}^{\prime }}(z,{\mu }_0,m_{\pi })=Z_{\mathcal{O},0}^{{\mathrm{RI}}^{\prime }}(z,{\mu }_0)+m_{\pi }^2{Z}_{\mathcal{O},1}^{{\mathrm{RI}}^{\prime }}(z,{\mu }_0).$$

$${\mathcal{R}}_{\mathcal{O}}(z)\equiv \frac{\mathrm{Re}[Z_{\mathcal{O},64}^{{\mathrm{RI}}^{\prime }}(z,{\mu }_0,m_{\pi })]}{\mathrm{Re}[Z_{\mathcal{O},48}^{{\mathrm{RI}}^{\prime }}(z,{\mu }_0,m_{\pi })]},{\mathcal{I}}_{\mathcal{O}}(z)\equiv \frac{\mathrm{Im}[Z_{\mathcal{O},64}^{{\mathrm{RI}}^{\prime }}(z,{\mu }_0,m_{\pi })]}{\mathrm{Im}[Z_{\mathcal{O},48}^{{\mathrm{RI}}^{\prime }}\mathcal{O}(z,{\mu }_0,m_{\pi })]},$$

$$Z_{\mathcal{O}}^{\overline{\mathrm{MS}}}(z,\bar{\mu },{\mu }_0)={\mathcal{Z}}_{\mathcal{O},0}^{\overline{\mathrm{MS}}}(z,\bar{\mu })+(a{\mu }_0{)}^2{\mathcal{Z}}_{\mathcal{O},1}^{\overline{\mathrm{MS}}}(z,\bar{\mu }).$$

$${\mathcal{Z}}_{\mathcal{O},0}^{\mathrm{M}\overline{\mathrm{MS}}}(z,\bar{\mu })={\mathcal{Z}}_{\mathcal{O},0}^{\overline{\mathrm{MS}}}(z,\bar{\mu }){\mathcal{C}}^{\overline{\mathrm{MS}},\mathrm{M}\overline{\mathrm{MS}}},$$

$$\displaystyle{\cal C}_{\cal O}^{\overline{\rm MS},{\rm M\overline{MS}}}{=}1+\frac{C_{F}g^{2}}{16\pi^{2}}\Bigg{[}$$

$$\frac{Z^{{\rm RI}^{\prime}}_{\cal O}(z,\mu_{0},m_{\pi})}{Z^{{\rm RI}^{\prime}}_{q}(\mu_{0},m_{\pi})-A^{\infty}_{q}(\mu_{0})}\frac{1}{12}{\rm Tr}\left[{\cal V}(z,p,m_{\pi})\left({\cal V}^{\rm Born}(z,p)\right)^{-1}\right]\Bigr{|}_{p^{2}{=}\mu_{0}^{2}}{=}1\,,$$

$$h_{\Gamma }^{\mathrm{M}\overline{\mathrm{MS}}}(P_3,z,\bar{\mu })=h_{\Gamma }^{bare}(P_3,z)\cdot {\mathcal{Z}}_{\Gamma ,0}^{\mathrm{M}\overline{\mathrm{MS}}}(z,\bar{\mu }),$$

$$\tilde{q}(x,P_3)=\frac{2P_3}{4\pi }\sum _{z=-z_{\max }}^{z_{\max }}{e}^{-ix{P}_{3}z}{h}_{\Gamma }(P_3,z).$$

$$q_{{\gamma }^0}(x,\Lambda )=q_{bare}(x)+\frac{{\alpha }_s}{2\pi }{\int }_x^1\frac{dy}{y}\left({\Pi }_{{\gamma }^0}\left(\Lambda \right)\delta \left(1-\frac{x}{y}\right)+{\Gamma }_{{\gamma }^0}\left(\frac{x}{y},\Lambda \right)\right)q_{bare}(y),$$

$${\Gamma }_{{\gamma }^0}^{\overline{\mathrm{MS}}}\left(\xi ,\frac{\bar{\mu }}{{\mu }_F}\right)=\frac{1+{\xi }^2}{1-\xi }\left(-\frac{1}{{ϵ}_{IR}}+\ln \left(\frac{{\bar{\mu }}^2}{{\mu }_F^2}\right)\right),$$

Full text

DESY 19-001

**Systematic uncertainties in parton distribution functions from lattice QCD simulations at the physical point **

Constantia Alexandrou(a,b), Krzysztof Cichy*(c), Martha Constantinou(d)*,

Kyriakos Hadjiyiannakou*(b), Karl Jansen(e), Aurora Scapellato(a,f), Fernanda Steffens(g)* **

(a) Department of Physics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus

(b) Computation-based Science and Technology Research Center, The Cyprus Institute, 20 Kavafi Str., Nicosia 2121, Cyprus

(c) Faculty of Physics, Adam Mickiewicz University, Umultowska 85, 61-614 Poznań, Poland

(d) Temple University,1925 N. 12th Street, Philadelphia, PA 19122-1801, USA

(e) NIC, DESY, Platanenallee 6, D-15738 Zeuthen, Germany

(f) University of Wuppertal, Gaußstr. 20, 42119 Wuppertal, Germany

(g) Institut für Strahlen- und Kernphysik, Rheinische Friedrich-Wilhelms-Universität Bonn, Nussallee 14-16, 53115 Bonn

We present a detailed study of the helicity-dependent and helicity-independent collinear parton distribution functions (PDFs) of the nucleon, using the quasi-PDF approach. The lattice QCD computation is performed employing twisted mass fermions with a physical value of the light quark mass. We give a systematic and in-depth account of the salient features entering in the evaluation of quasi-PDFs and their relation to the light-cone PDFs. In particular, we give details for the computation of the matrix elements, including the study of the various sources of systematic uncertainties, such as excited states contamination. In addition, we discuss the non-perturbative renormalization scheme used here and its systematics, effects of truncating the Fourier transform and different matching prescriptions. Finally, we show improved results for the PDFs and discuss future directions, challenges and prospects for evaluating precisely PDFs from lattice QCD with fully quantified uncertainties.

I Introduction

Quantum Chromodynamics (QCD), being the fundamental theory of the strong interactions, describes the interactions among quarks and gluons both in the perturbative regime, as well as in the non-perturbative regime with the emergence of complex structures like hadrons. As the latter requires a non-perturbative approach, the use of a discretized Euclidean spacetime to numerically solve QCD has been proven to be extremely successful in the computation of physical properties of hadrons, such as their masses and decay constants. The rich and diverse structure of hadrons is described in terms of observables like various kinds of form factors, revealing different structural aspects, such as the average distribution of electric charge or parton (quark and gluon) momentum and spin inside the hadron. Even more information is contained in distribution functions, such as parton distribution functions (PDFs), generalized parton distributions (GPDs) and transverse momentum dependents PDFs (TMDs), which decompose the aggregate observables into probability densities for partons with specific longitudinal momentum, transverse momentum and momentum transfer within a scattering event. The simplest of these are the PDFs, which have been studied intensively and continuously in experimental facilities over the last few decades, most notably for the proton, and also provide input in collider experiments. It is well-established that the main source of information on PDFs are global QCD analyses, providing accurate results due to theoretical advances and new data emerging from accelerators, covering different kinematical regions (see e.g. Refs. Perez:2012um ; DeRoeck:2011na ; Alekhin:2011sk ; Ball:2012wy ; Aidala:2012mv ; Forte:2013wc ; Jimenez-Delgado:2013sma ; Rojo:2015acz ; Butterworth:2015oua ; Accardi:2016ndt ; Gao:2017yyd ; Lin:2017snn ). Such analyses rely on the factorization framework, in which scattering cross sections, generally obtaining contributions from all energy scales, are written as convolutions of perturbatively computable coefficient functions describing the high-energy scales, and the low-energy PDFs. Despite the tremendous progress in phenomenological parametrizations of PDFs, the procedures are not without ambiguities Jimenez-Delgado:2013sma , as there are still kinematical regions that are not easily accessible experimentally, e.g. the large Bjorken-$x$ region. Uncertainties on the PDFs in the large-$x$ range propagate to smaller-$x$ regions when QCD evolution is applied.

A calculation of PDFs from first principles is a valuable addition to the global fitting analyses and of crucial importance for the deeper understanding of the inner structure of hadrons. Apart from providing a fundamental framework for the study of these quantities, it can also serve as input for experimental analyses in collision experiments. The non-perturbative nature of PDFs makes lattice QCD an ideal ab initio formulation to determine them, utilizing large-scale simulations. A novel method to extract parton distribution functions from lattice QCD was proposed five years ago by Ji Ji:2013dva . It is based on considering matrix elements probing purely spatial correlations, making them accessible in Euclidean lattice QCD.

The first studies on quasi-PDFs within lattice QCD Lin:2014zya ; Alexandrou:2014pna ; Alexandrou:2015rja have appeared soon after their proposal by Ji. These papers and the follow-up ones Chen:2016utp ; Alexandrou:2016jqi focused mainly on the feasibility of the approach only calculating bare nucleon matrix elements of boosted nucleons for nonsinglet operators. The first computations applying a proper procedure for renormalization and matching to $\overline{\mathrm{MS}}$ light-cone PDFs and using simulations with physical pion mass were published in Ref. Alexandrou:2018pbm for the unpolarized isovector PDFs, $u(x)-d(x)$, and helicity, $\Delta u(x)-\Delta d(x)$ and in Ref. Alexandrou:2018eet for the isovector transversity PDF, $h_{1}^{u}(x)-h_{1}^{d}(x)$. These works works were published in shorter letters, where many of the technical and methodological details of the calculations needed to be left out.

It is the purpose of this paper to fill this gap and provide the before mentioned details in a comprehensive way, give the present understanding and control of systematics effects appearing in the computations and to present the status of the calculations of PDFs for our twisted mass setup. In particular, we extend the work of Refs. Alexandrou:2018pbm ; Alexandrou:2018eet and give a systematic and in-depth account of the salient features entering in the evaluation of quasi-PDFs and their relation to the light-cone PDFs. We explain in detail the computation of the matrix elements including the study of the various sources of systematic uncertainties, for example excited states contamination. An improved renormalization procedure is presented using three ensembles for the non-perturbative extraction of the renormalization functions, allowing to take the chiral limit, as well as an investigation of other sources of systematic uncertainties related to the renormalization, such as the finite volume and scale dependence. In addition, different prescriptions for applying the Fourier transform and the matching prescription are employed and critically compared.

On the physics side, in order to compare with existing phenomenological estimates, we convert our results in the RI*′* scheme into the $\overline{\mathrm{MS}}$ scheme following the procedure given in Ref. Constantinou:2017sej . We then apply appropriate matching to relate quasi-PDFs in the $\overline{\mathrm{MS}}$ scheme to light-cone PDFs extracted from global fits. With the detailed analysis of systematic effects, as decribed in this work, having a physical value of the pion mass and the improved renormalization and matching prescriptions, we are finally in position to make indeed contact to the phenomenological analyses of PDFs and we present final results of our lattice PDF calculations in Sec. VI below. In addition to our work, results on the quasi-distribution approach were reported in Refs. Zhang:2017bzy ; Chen:2017mzz ; Lin:2017ani ; Chen:2017gck ; Chen:2018xof ; Chen:2018fwa ; Lin:2018qky ; Liu:2018uuj ; Fan:2018dxu ; Liu:2018hxv ; Petreczky:2018jqc that include, besides the nucleon PDFs, exploratory studies of pion distribution amplitudes and PDFs as well as of gluon PDFs. For complete references and an overview of the theoretical and numerical developments, see Ref. Cichy:2018mum .

On more general grounds, it is very important to realize that quasi-PDFs and light-cone PDFs have been shown to share the same infrared physics Xiong:2013bka ; Briceno:2017cpo , which is the fundamental observation that allows one to relate both quantities using perturbation theory, provided that the hadron is moving with a large, although necessarily finite, momentum in a chosen spatial direction. It has also been proven that quasi-PDFs can be extracted from lattice QCD in Euclidean spacetime Briceno:2017cpo and that they do not suffer from power-divergent mixings with lower-dimensional operators Ji:2017rah ; Radyushkin:2018nbf ; Karpie:2018zaz . A factorization formula makes it possible to extract the PDFs from the quasi-PDFs, an operation called matching Xiong:2013bka ; Ji:2015qla ; Xiong:2015nua ; Wang:2017qyg ; Stewart:2017tvs ; Izubuchi:2018srq ; Alexandrou:2018pbm ; Alexandrou:2018eet ; Liu:2018uuj ; Liu:2018hxv . In general, this procedure is based on a newly developed large-momentum effective theory (LaMET) Ji:2014gla , and it is renormalizable to all orders in perturbation theory Ishikawa:2017faj ; Ji:2017oey ; Zhang:2018diq ; Li:2018tpe . Other approaches for a direct computation of the $x$-dependence of PDFs include the hadronic tensor Liu:1993cv ; Liu:1998um ; Liu:1999ak , fictitious heavy quark Detmold:2005gg , auxiliary light quark Braun:2007wv , good lattice cross sections Ma:2014jla ; Ma:2017pxb (closely related to the auxiliary light quark method), “OPE without OPE” Chambers:2017dov and pseudo-PDFs Radyushkin:2016hsy ; Radyushkin:2017ffo ; Radyushkin:2017cyf ; Radyushkin:2017lvu , where the latter can be seen as a generalization of PDFs off the light-cone. These alternative approaches have been explored in lattice QCD, and recent results can be found in Refs. Liang:2017mye ; Detmold:2018kwu ; Bali:2017gfr ; Bali:2018spj ; Chambers:2017dov ; Orginos:2017kos ; Karpie:2017bzm ; Radyushkin:2018cvn ; Karpie:2018zaz . The new formulation and its successful implementation within lattice QCD has also led to a wider interest on phenomenological studies using models and toy theories of QCD Gamberg:2014zwa ; Gamberg:2015opc ; Jia:2015pxx ; Jia:2018qee ; Bacchetta:2016zjm ; Radyushkin:2016hsy ; Radyushkin:2017ffo ; Broniowski:2017wbr ; Broniowski:2017gfp ; Bhattacharya:2018zxi . A detailed overview of the current status of lattice QCD calculation of PDFs and other partonic distributions can be found in the recent reviews of Refs. Monahan:2018euv ; Cichy:2018mum .

The remainder of the paper is organized as follows: In Sec. II, we provide the general theoretical aspects, lattice QCD action and parameters. We discuss the computation of the bare matrix elements, along with the lattice techniques necessary for attaining high statistical accuracy. In Sec. III, we present our results on bare matrix elements, together with various analyses methods that are required to control excited states contamination. Sec. IV describes in detail the non-perturbative renormalization program, which investigates pion mass dependence, volume effects and renormalization scale dependence. The matching procedure is presented in Sec. V and we demonstrate that a modified $\overline{\mathrm{MS}}$ scheme is preferable. Different prescriptions are presented and we compare the effect on the final PDFs. Sec. V also includes an investigation of systematic uncertainties resulting from the Fourier transform. In Sec. VI, we present our final results with discussion on the systematic uncertainties, while in Sec. VII we conclude.

II Methodology

II.1 From PDF to quasi-PDF

The original definition of quark momentum distribution within a hadron can be derived from the operator product expansion (OPE) of hadronic deep-inelastic scattering and is given by

[TABLE]

where $x$ is the momentum fraction of the quark, $|N\rangle$ the hadron state at rest assumed to be the nucleon, the momentum $P^{+}{=}(P^{0}+P^{3})/\sqrt{2}{=}M/\sqrt{2}$, $\xi^{-}{=}(\xi^{0}-\xi^{3})/\sqrt{2}$, $W(\xi^{-},0)$ is the Wilson line connecting $\xi^{-}$ and [math], and $\Gamma$ is a Dirac matrix whose structure defines the momentum distributions of quarks having spin parallel or perpendicular to the nucleon momentum. The factorization scale is implicit. Eq. (1) is light-cone dominated, as it receives contributions only in the region where $\xi^{2}{=}t^{2}-\vec{r}^{\,2}{=}0$. Consequently, a direct implementation of this equation is impossible on the lattice, where only a single point can satisfy the Euclidean light-cone condition $\xi^{2}{=}t^{2}+\vec{r}^{\,2}{=}0$. The alternative approach proposed by Ji overcomes this difficulty by reconstructing light-cone distributions from correlation functions where the bilinear fermionic fields, $\psi$ and $\bar{\psi}$, are separated by a purely spatial distance $z$ and the nucleon is boosted with a finite momentum. Quasi-PDFs are defined as

[TABLE]

where $|N\rangle$ represents a nucleon, which is boosted in the $z$-direction with momentum $P{=}(P_{0},0,0,P_{3})$. On the lattice, quasi-PDFs are thus computed using matrix elements of non-local operators containing a straight Wilson line of finite length $z$, that varies from 0 to some maximum value, $z_{\rm max}$. This aspect will be further discussed in the following sections. The standard light-cone PDFs can be recovered from quasi-PDFs using the Large Momentum Effective Theory Ji:2014gla , developed specifically for these operators. Eventually, quasi-PDFs are matched to the physical PDFs through the factorization

[TABLE]

where $q(x,\mu)$ is the light-cone PDF at the scale $\mu$ and we follow the usual choice of common factorization and renormalization scales, $\mu_{F}=\mu_{R}\equiv\mu$. $C$ is the matching kernel, which can be computed perturbatively and has so far been evaluated to one-loop level. For sufficiently large momenta, ideally much larger than the nucleon mass, $m_{N}$, and $\Lambda_{QCD}$, higher order corrections in $\alpha_{s}$ are very small and the extraction of quark distributions is fully robust. Presently, such a situation cannot be accomplished and it is part of this paper to investigate whether practically reachable momenta are sufficient.

II.2 Lattice gauge ensemble

The analysis is performed using a gauge ensemble of two dynamical degenerate light quarks ($N_{f}{=}2$), generated by the ETM 111Effective from this year, the European Twisted Mass Collaboration has officially changed its name to Extended Twisted Mass Collaboration, as it comprises now members also from non-European institutions. Along with the name change, there is also a new logo. Collaboration Abdel-Rehim:2015pwa . The light quark mass has been tuned to reproduce the physical value of the pion mass, henceforth referred to as the physical point. The lattice volume is $48^{3}\times 96$, the lattice spacing $a=0.0938(2)(3)$ fm Alexandrou:2017xwd , the spatial lattice extent $L\approx 4.5$ fm and $\,\,m_{\pi}L=2.98$. Based on a model prediction of excited states effects in certain hadron structure quantities Hansen:2016qoz , the dependence on $m_{\pi}L$ is mild (see, e.g., Fig. 3 of Ref. Green:2018vxw ). The complete list of parameters of this ensemble is reported in Table 1.

The gauge configurations were generated with the Iwasaki improved gauge action Iwasaki:2011np . In the fermionic sector, the twisted mass fermion action Frezzotti:2000nk ; Frezzotti:2003ni with a clover term Sheikholeslami:1985ij was employed. The fermion action is given by:

[TABLE]

where $D_{W}$ is the Wilson-Dirac operator, $\mu_{l}$ is the bare twisted mass for the light quarks, $\tau^{3}={\rm diag}(1,-1)$ is the third Pauli matrix in flavor space and the last term includes the field strength tensor $\mathcal{F}^{\mu\nu}[U]$ weighted by $c_{\rm SW}$, known as the Sheikoleslami-Wohlert (clover) coefficient. In the gauge ensemble used in this work, the clover parameter is taken to be $c_{\rm SW}{=}1.57551$ Aoki:2005et . In Eq. (4) $\chi(x){=}(u,d)^{T}$ denotes the light quark doublet in the “twisted basis” at maximal twist. The fermion fields in the “physical basis”, denoted by $\psi(x)$, are recovered by the following chiral rotation

[TABLE]

where $\alpha{=}\pi/2$ at maximal twist. In the next sections of this paper, the interpolating fields and the nucleon matrix elements have to be understood with quark fields in the physical basis, unless otherwise specified.

The introduction of the twisted mass term in the lattice action has a series of advantages in hadron structure calculations, such as excluding zero eigenvalues from the spectrum of the Wilson Dirac operator allowing a speed-up in the numerical simulations. Moreover, it also simplifies renormalization properties of operators Frezzotti:2003ni ; Jansen:2005cg ; Farchioni:2004us ; Farchioni:2005bh and at maximal twist it provides an automatic $O(a)$ improvement. However, the isospin symmetry breaking of the twisted mass formulation could lead to instabilities when simulations are carried out at light quark masses close to the physical value. The isospin breaking, which manifests itself in the neutral pion being lighter than the charged one, is reduced with the addition of the clover term in the fermion action Abdel-Rehim:2015pwa .

II.3 Nucleon bare matrix elements

In this work, we evaluate the unpolarized, helicity and transversity quasi-PDFs using three values for the nucleon momentum, that is, $P_{3}{=}6\pi/L$, $P_{3}{=}8\pi/L$ and $P_{3}{=}10\pi/L$, corresponding in physical units to 0.83 GeV, 1.11 GeV and 1.38 GeV, respectively. The matrix element of interest is given by

[TABLE]

for a straight Wilson line, $W$, with varying length from $z{=}0$, corresponding to the standard ultra-local operators, up to half of the spatial extension, $L/2$. $\Gamma$ is the Dirac structure leading to different types of PDFs, as discussed below. The matrix elements of Eq. (6) are extracted from a ratio of two- and three-point functions, averaged over the gauge field ensemble. The two-point and three-point functions are given by

[TABLE]

where $N_{\alpha}(x)$ is the proton interpolating field, $t$ ($t_{s}$) is the time separation of the sink relative to the source in the two-point (three-point) function222We use different symbols to emphasize that the three-point function is computed only for selected values of the source-sink separation, $t_{s}/a{=}8,9,10,12$ in this work, while the two-point function is evaluated for the full range of time separations to constrain more precisely the gap between the ground state and the first excited state in two-state fits, see Sec. III.3., $\tau$ is the insertion time of the $\mathcal{O}$ operator and ${\cal P}_{\alpha\beta}$, $\tilde{\cal P}_{\alpha\beta}$ are the parity projectors for the two- and three-point functions, respectively. For the two-point functions, we use the plus and minus parity projectors ${\cal P}^{\pm}=\frac{1\pm\gamma^{0}}{2}$ and average over the forward and backward correlators. The projector $\tilde{\cal P}_{\alpha\beta}$ depends on the operator under study. The operator $\mathcal{O}(\mathbf{y},\tau;z)$ has the following form

[TABLE]

where the Wilson line is aligned to the direction of the nucleon momentum and $\psi$ is the fermion doublet of flavors $(u,d)$. We compute the isovector combination $u{-}d$, inserting a $\tau_{3}$ in flavor space. With this choice, the disconnected diagrams cancel and the total contribution is obtained from the connected diagram shown in Fig. 1 .

As mentioned in Section II.1, the choice for the Dirac structure in the operator of Eq. (1) gives access to different quark distribution functions. In fact, to extract a given PDF, different $\Gamma$ matrices may be used, for example $\gamma^{3}$ (parallel to the Wilson line direction) or $\gamma^{0}$ (temporal direction) for the unpolarized PDF. However, some choices are preferable (e.g., $\gamma^{0}$) as they avoid finite mixing under renormalization, that was found to be allowed in lattice regularization Constantinou:2017sej . In the absence of mixing, the operators renormalize multiplicatively, and one avoids systematic uncertainties related to the elimination of mixing. In this work, we compute matrix elements of the operators with

• •

${\Gamma}=\gamma^{3}\,,\gamma^{0}$ for the unpolarized distribution, $\tilde{q}(x)=\tilde{q}(x)_{\uparrow}+\tilde{q}_{\downarrow}(x)$ ,

• •

${\Gamma}=\gamma^{5}\gamma^{3}$ for the helicity distribution, $\Delta\tilde{q}(x)=\tilde{q}_{\uparrow}(x)-\tilde{q}_{\downarrow}(x)$ ,

• •

${\Gamma}=\sigma^{3j}$ for the transversity distribution, \delta\tilde{q}(x)=\tilde{q}_{\perp}(x)+\tilde{q}_{\scriptsize{\rotatebox[origin={c}]{-180.0}{\perp}}}(x) .

$q_{\uparrow}$ ($q_{\downarrow}$) and ${q}_{\perp}$ ({q}_{\scriptsize{\rotatebox[origin={c}]{-180.0}{\perp}}}) indicate quarks with helicity aligned (anti-aligned) with that of a longitudinally and transversely polarized proton, respectively. The index $j$ entering the matrix element of the transversity distribution is purely spatial and denotes the direction of the quark spin, orthogonal to the proton momentum. Each choice of an operator requires an appropriate parity projector $\tilde{\cal P}_{\alpha\beta}$, that is, $\frac{1+\gamma^{0}}{2}$ for the unpolarized, $i\gamma^{3}\gamma^{5}\frac{1+\gamma^{0}}{2}$ for the helicity, and $i\gamma^{k}\frac{1+\gamma^{0}}{2}$ for the transversity PDF. The desired matrix elements of Eq. (6) are finally extracted from fitting the ratio of three- over two-point functions to isolate the ground state. One choice is a constant fit as a function of the time insertion of the operator over which the ratio is time independent (plateau region). Other choices are two-state fits that take into account the first excited state contributions, as discussed in Sec. III with a comparison among them. Thus, the desired matrix element in the plateau fit is given by

[TABLE]

where $\mathcal{K}$ is a kinematic factor equal to $\mathcal{K}\,{=}\,iE/p$ (with $E{=}\sqrt{m_{N}^{2}+p^{2}}$) for the vector current $\gamma^{3}$, whereas $\mathcal{K}{=}1$ for all other choices of the Dirac structure used. A gauge average ($\langle\cdots\rangle_{G}$) is performed on the two- and three-point functions prior to taking the ratio.

II.4 Lattice techniques

To evaluate the correlation functions of Eqs. (7),(8), a state with the same quantum numbers as the nucleon is created at the source and annihilated at the sink. In our calculation, we employ point sources generated with the proton interpolating field $N_{\alpha}(x){=}\epsilon^{abc}u^{a}_{\alpha}(x)\left(d^{b^{T}}(x)\mathcal{C}\gamma^{5}u^{c}(x)\right)$, with $\mathcal{C}{=}\gamma^{0}\gamma^{2}$ the charge conjugation matrix. In general, the overlap of the states generated with such interpolating fields with the desired ground state is improved by employing smearing. Gaussian smearing Gusken:1989qx ; Alexandrou:1992ti on quark fields, used in combination with APE smearing on the gauge field, is a very effective approach to improve ground state dominance. For nucleon states at rest, previous studies Abdel-Rehim:2015owa performed on the gauge ensemble employed in this work, have extracted the optimal parameters for both Gaussian and APE smearing, tuned to approximately reproduce a nucleon state with root mean square (rms) radius of about 0.5 fm and maximize the overlap with the ground state. These parameters are $(N_{G},\alpha_{G})=(50,4)$ for Gaussian smearing and $(N_{APE},\alpha_{APE}=0.5)$ for APE smearing.

However, in the computation of PDFs, we need to optimize the overlap to a boosted nucleon state keeping the statistical noise minimal. We, thus, modify the standard Gaussian smearing by using momentum-smeared interpolating fields Bali:2016lva . This has proven extremely advantageous when the nucleon is boosted with a large momentum, as is the case for quasi-PDFs. The momentum smearing modifies the standard Gaussian smearing function by including a complex phase factor that affects only the gauge links along the direction of the boost, namely

[TABLE]

where $U_{j}$ denotes a gauge link in the $j$-direction. The value of the free momentum smearing parameter $\xi$ depends, in general, on the parameters of the gauge ensemble and on the nucleon momentum. Tuning of $\xi$ is necessary in order to improve the overlap of our interpolator with the boosted proton. We, thus, have optimized the parameter $\xi$ for each value of the momentum employed by minimizing the statistical errors of the nucleon two-point functions. In Fig. 2, we demonstrate the effect of the momentum smearing, by plotting the scaling of the error of the two-point correlator and the effective energy. The results shown have been extracted using $100$ measurements for the nucleon boost $P_{3}{=}8\pi/L$. For comparison, we include also the results using the standard Gausian smearing ($\xi{=}0$). As can be seen, the errors in the correlation functions reduce dramatically as the value of $\xi$ increases and convergence is observed in the range $\xi\in[0.6{-}0.75]$. Thus, any value of $\xi$ in this window of values leads to a similar signal-to-noise ratio. The tuning procedure for the other two momenta, $P_{3}{=}6\pi/L$ and $P_{3}{=}10\pi/L$, leads to the same conclusions and we fix $\xi{=}0.6$ through out this work.

The three-point functions are computed using the sequential method Martinelli:1988rr , which has the advantage of summing automatically the sink spatial volume to produce the sequential propagator for all insertion points. Compared to the stochastic method Alexandrou:2013xon , the sequential method does not introduce stochastic noise, but has the drawback that new inversions of the Dirac matrix have to be carried out for different source-sink separations and nucleon momenta. For details on the comparison between the stochastic method (with pure Gaussian smearing) and the sequential method (with and without momentum smearing), we refer to Refs. Alexandrou:2014pna ; Alexandrou:2016jqi . In addition, the sequential method is preferable when the momentum smearing technique on the quark fields is employed.

To increase the number of measurements of the three-point functions, we use multiple source positions ($N_{src}$) for each configuration. To confirm that the data extracted from different source positions on the same configuration are statistically independent, we study the scaling of the statistical error with the increase of $N_{src}$. In Fig. 3 (left), we show the error on the matrix elements for $z{=}0$, obtained from the analysis of a sample of gauge configurations for nucleon boost $6\pi/L$. For comparison purposes, the absolute error is normalized to the one obtained using only one source position per configuration. As can be seen, the errors for the three choices of $\Gamma$ structure decrease approximately as $1/\sqrt{N_{src}}$, which is the scaling expected if the measurements are not correlated.

To further decrease the statistical uncertainties, we boost the nucleon along different spatial directions and orientations, that is, $\pm x$, $\pm y$ and $\pm z$, with the Wilson line taken always along the axis in which the spatial component of the momentum is nonzero. The correlation functions obtained from these $N_{dir}{=}6$ possible directions lead to the same physical results due to the spatial rotational symmetry on the lattice, and therefore can be averaged within each configuration. The statistical error reduction fluctuates is close to the ideal behavior $1/\sqrt{N_{dir}}$, as shown in the right panel of Fig. 3.

Despite the use of the momentum and APE smearing, the exponential decrease of the signal-to-noise ratio persists as the nucleon momentum and source-sink separation increase. To increase statistics at a reduced cost, we employ, together with the momentum and APE smearing, the Covariant Approximation Averaging (CAA) Blum:2012uh that belongs to the class of truncated solver methods with a bias correction. For each configuration, $N_{LP}$ low-precision (LP) inversions of the Dirac matrix are carried out from a set of random source positions $\left\{X_{src}\right\}$ and the bias from the measurements is removed using a small number $N_{HP}$ of high-precision (HP) inversions. Denoting with $C_{LP}$ and $C_{HP}$ the correlation functions produced with low and high precision of the solver, respectively, the improved correlation functions $C^{imp}$ for each configuration are defined by

[TABLE]

where, in the second sum, $C_{LP}$ and $C_{HP}$ are computed on the same source position, otherwise the bias cannot be corrected. The error for a given observable scales with the ratio $N_{HP}/N_{LP}$ as

[TABLE]

where $r_{c}$ is the correlation coefficient among nucleon correlators computed at high and low precision. A compromise is needed to have $r_{c}\simeq 1$, while keeping the inversions as fast as possible and $N_{LP}\gg N_{HP}$. Thus, a tuning of the precision of the solver has to be carried out. To invert the Dirac operator, we use the adaptive multigrid solver with twisted mass fermion support Alexandrou:2016izb and require the residual to be $r_{HP}{=}10^{-10}$ for HP inversions. After testing different values of the residual for LP inversions, we find that the stopping criterion

[TABLE]

guarantees a correlation coefficient $r_{c}\geq 0.999$ with a considerable speed-up in the inversion time. Moreover, taking 15 HP inversions as the reference setup, a comparison of the HP and CAA estimates for the two- and three-point correlators verified that the bias introduced from LP inversions is negligible compared to the gauge noise of our measurements when one HP inversion for each configuration is performed. Thus, to extract the nucleon matrix elements for quasi-PDFs at momenta $8\pi/L$ and $10\pi/L$, we use the CAA setup with $(N_{HP},N_{LP}){=}(1,16)$.

The quasi-PDFs are computed at different source-sink time separations $t_{s}$ in order to investigate excited states effects. In particular, we use $t_{s}/a{=}8,9,10,12$ for the unpolarized and $t_{s}/a{=}8,10,12$ for the helicity and transversity cases. For a detailed discussion on the excited states effects, we refer to Sec. III.3. In the remaining part of this section and unless otherwise stated, we focus on the results extracted from $t_{s}{=}12a$, which is the one where excited states are found to be suppressed for all three PDFs when the statistical precision is around $10\%$. The number of configurations and the total statistics for each operator and momentum is reported in Table 2. In the number of total measurements, we include a factor of $6$, coming from the average over correlators computed from the boost aligned along $\pm x$, $\pm y$, $\pm z$-directions, and a factor 16 (15) from the source positions for $P_{3}=6\pi/L$ ($P_{3}=[8,10]\pi/L$). This translates into an additional factor of $96$ for $P_{3}{=}6\pi/L$ and $90$ for $P_{3}{=}8,10\pi/L$, where the CAA is employed. Moreover, we note that only for the transversity PDF at $P_{3}{=}6\pi/L$, we have averaged over the matrix elements computed for the two possible choices of Dirac matrices (for example, $\sigma^{31}$ and $\sigma^{32}$ for a boost in the $z$-direction). This contributes with an additional factor of 2 in the total number of measurements at this momentum, as well as in the computational cost (due to different parity projectors needed for $\sigma^{31}$ and $\sigma^{32}$).

Although large nucleon momenta are needed to approach light-cone PDFs, high values of the momentum on a Euclidean lattice may lead to substantial cut-off effects if the condition $P_{3}\ll 1/a$ is not satisfied. One possible check of cutoff effects is via the computation of the dispersion relation. To this end, we compute the nucleon energies for momenta {0,2,4,6,8,10}$\pi/L$ using the momentum smearing method, and check whether the relativistic dispersion relation $E^{2}{=}m_{N}^{2}c^{4}+P_{3}^{2}c^{2}$ is satisfied for our results. As can be seen in Fig. 4, no deviations from the continuum energy-momentum relation are observed for the values employed in this work (up to 1.38 GeV), which is well below the inverse lattice spacing ($10\pi/L\approx 0.65/a$). Moreover, by performing a two-parameter fit to the lattice data, we obtain $c^{2}{=}1.00(3)$ for the squared speed of light and $a^{2}m_{N}^{2}c^{4}{=}0.207(4)$ for the nucleon mass in lattice units, which is compatible with the value $am_{N}{=}0.4436(11)$ extracted at zero-momentum Alexandrou:2017xwd . Whether this finding of only small cutoff effects in the dispersion relation also holds in case of the PDFs we are interested in, can eventually only be answered when results at several lattice spacing are available.

III Techniques for the evaluation of bare quasi-PDF nucleon matrix elements in lattice QCD

In this section, we discuss crucial aspects related to the lattice QCD computation of bare quasi-PDF nucleon matrix elements, such as dependence on the number of stout smearing iterations used in the operator, the choice of the Dirac structure and identification of excited states contamination.

III.1 Stout smearing

We apply three-dimensional stout smearing to the gauge links of the Wilson line of the operator, following the prescription of Ref. Morningstar:2003gk . This reduces the power divergence that is present in the matrix elements of fermion operators with Wilson lines (see, e.g. Refs. Dotsenko:1979wR ; Brandt:1981kf , as well as Ref. Cichy:2018mum for a review of recent investigations of the power divergence). The application of smearing was especially crucial in the first calculations of quasi-PDFs Lin:2014zya ; Alexandrou:2015rja ; Chen:2016utp ; Alexandrou:2016jqi when the renormalization procedure was not yet developed. Even though the complete renormalization procedure was developed recently Alexandrou:2017huk , a few iterations of stout smearing is still useful for noise reduction in the renormalized matrix elements. In Fig. 5, we show examples of the effect of stout smearing for the case of the bare matrix elements of the unpolarized (insertion $\gamma^{0}$), helicity and transversity quasi-PDFs for nucleon momentum $P_{3}{=}8\pi/L$ without stout smearing, using 0, 5, 10 and 15 stout steps. As expected, the stout smearing modifies the matrix elements, increasing the values of the real and imaginary parts at each $z$. We find convergence of the matrix elements after a few iterations of stout smearing and, thus, in what follows we discuss in detail results obtained with $5$ steps of stout. The renormalized matrix elements are expected to no longer show any dependence on the smearing. We discuss the details of the renormalization in Sec. IV.

III.2 Choice of the Dirac structure

Although the natural choice to extract the unpolarized PDF would be $\gamma^{+}{=}(\gamma^{0}{+}\gamma^{3})/\sqrt{2}$, a recent study Constantinou:2017sej has revealed that the matrix element with $\gamma^{3}$ exhibits mixing with the twist-3 scalar operator (later confirmed by symmetry properties Chen:2017mzz ) and the computation of a mixing renormalization matrix is required Alexandrou:2017huk . However, the twisted mass formulation has the advantage that the mixing is between the vector and pseudoscalar operator. The latter vanishes in the continuum limit and only contributes as a discretization effect increasing the gauge noise. Consequently, it may be neglected in the first approximation, since eventually one is interested in the results for $a{\rightarrow}0$. In this work, we compute both matrix elements, but use the data for $\gamma^{0}$ to extract our final results for the unpolarized quasi-PDF. Even though the mixing is a purely lattice artifact for twisted mass fermions, we expect increased noise contamination for the case of $\gamma^{3}$, which is clearly visible in the lattice data, as shown in Fig. 6. This behavior is momentum-independent and we compare the matrix elements for the two $\gamma$-structure at momentum $8\pi/L$. The statistical errors for $\gamma^{3}$ are twice larger as compared to the ones for $\gamma^{0}$ with the same number of measurements. We, thus, conclude that the insertions $\gamma^{3}$ or $\gamma^{+}$ are not optimal for extracting the unpolarized PDFs within a lattice QCD calculation.

We compute the matrix elements for the three values of the nucleon momentum given in Table 2 with the associated number of measurements. The momentum dependence of the resulting matrix elements is shown in Fig. 7 for the three operators considered in this work. We find that with increasing momentum, the matrix elements decay faster to zero and, for the highest momentum employed, both real and imaginary parts are compatible with zero for $z\geq 11a\simeq 1.03$ fm.

III.3 Excited States Contamination

Identification of the nucleon state is crucial in order to extract the correct nucleon matrix elements from lattice QCD measurements. This requires a careful analysis of excited states. An additional challenge is the need to boost the nucleon with a relatively large momentum, something that it is not needed in e.g. studies of nucleon charges and form factors. Requiring large nucleon momentum in combination with using simulations at the physical point leads to a more severe contamination due to the excited states, since the spectrum is denser. Therefore, a thorough assessment of the excited states is even more essential for the reliable extraction of PDFs. We extend the study of excited states first presented in Ref. Alexandrou:2018yuy in order to eliminate, as much as possible, one of the major systematic uncertainties in the evaluation of nucleon matrix elements from which the PDFs are extracted within lattice QCD. In what follows, we review the analysis methods that we employ to isolate the ground state. As will be discussed in this subsection, the excited states contributions are milder for the matrix elements of the unpolarized operator as compared to the matrix elements of the helicity and transversity operators.

To extract the ground state contribution in the correlation functions, we employ three analysis methods, briefly described below.

1. Plateau method (single-state fit). In this method, one seeks for a region where the ratio of Eq. (10) becomes independent of the insertion time, $\tau$, which indicates suppression of excited states (possibly partly). The ratio in this region is fitted to a constant value, yielding the matrix element of the ground state. Indeed, inserting in the two- and three-point functions two sets of complete eigenstates of the QCD Hamiltonian ($|n(P_{3})\rangle$ and $|n^{\prime}(P_{3})\rangle$, where we indicated that the states are momentum-dependent, but we keep this dependence implicit below to simplify the notation) with quantum numbers of the nucleon, the ratio of Eq. (10) can be written as333All correlation functions in this section should be understood as averaged over the gauge field configurations ensemble.

[TABLE]

where $|N\rangle$ is the nucleon state, $E_{n,n^{\prime}}$ are energies of the eigenstates $|n\rangle$, $|n^{\prime}\rangle$ (with $|0\rangle$ – ground state, $|1\rangle$ – first excited state, etc.) and the timeslice of the source is set to zero. Isolating in the ratio the contribution of the ground state and expanding the sum up to the first excited state, Eq.(15) reads

[TABLE]

where we introduced the constant $f_{10}{=}{\langle N|1\rangle}/{\langle N|0\rangle}$ and $\Delta E=E_{1}-E_{0}$. As can be seen, the excited states contributions fall off exponentially and for $\Delta E\,t_{s}\gg 0$, $\Delta E\,\tau\gg 0$ and $\Delta E\,(t_{s}-\tau)\gg 0$, the first time-independent term dominates, yielding the matrix element of the nucleon state, $\langle 0|\mathcal{O}|0\rangle$. Thus, fitting the ratio within the plateau region to a constant value yields the desired nucleon matrix element. Since statistical errors grow exponentially with $t_{s}$, the challenge is to identify the smallest value of $t_{s}$ that suppresses the contributions of excited states in the ratio to a level negligible as compared to the statistical accuracy.

2. Summation method. This method was introduced in Ref. Maiani:1987by and entails the sum of the ratio of Eq. (16) over the insertion time $\tau$, excluding the timeslices of the source and the sink. The sum is a geometric series in the terms involving the excited states contributions and yields the expression

[TABLE]

where $C$ is a constant and the desired matrix element $\mathcal{M}{=}\langle 0|\mathcal{O}|0\rangle$ can be extracted from a linear two-parameter fit. As a consequence, this procedure yields, in general, results with larger uncertainties as compared to the plateau method, but has the advantage of suppressing the excited states contamination by a faster decaying factor $e^{-\Delta E\,t_{s}}$, as compared to the leading exponential factors in Eq. (16).

3. Two-state fits. In this method, one retains the first excited state contributions in the two-point and three-point correlation functions. We use two types of fits: (a) either a fit to the two-point correlator followed by a fit to the three-point function, or (b) a simultaneous fit to both the two- and three-point correlators.

We refer to the (a) type as a sequential fit and parameterize the nucleon two-point function with momentum $P_{3}$ as

[TABLE]

with fit parameters being the ground state amplitude, $|A_{0}|^{2}=|\langle N|0\rangle|^{2}$, the ground state energy $E_{0}$, $|f_{10}|^{2}$ and $\Delta E$. The three-point correlator can be written as

[TABLE]

The amplitudes $|A_{0}|^{2}$ and $|f_{10}|^{2}$ and energies $E_{0}$ and $\Delta E$ are determined from first fitting the two-point function in Eq. (18) and used as inputs into the three-point function, which is fitted separately for the real and imaginary parts. The fit parameters are $\textrm{Re}\langle 0|\mathcal{O}|0\rangle$, $\textrm{Re}\widetilde{f}_{10}{\equiv}\textrm{Re}f_{10}\langle 1|\mathcal{O}|0\rangle$, $\textrm{Re}\widetilde{f}_{01}{\equiv}\textrm{Re}f_{10}^{\dagger}\langle 0|\mathcal{O}|1\rangle$, $\textrm{Re}\langle 1|\mathcal{O}|1\rangle$ and $\textrm{Im}\langle 0|\mathcal{O}|0\rangle$, $\textrm{Im}\widetilde{f}_{10}$, $\textrm{Im}\widetilde{f}_{01}$, $\textrm{Im}\langle 1|\mathcal{O}|1\rangle$, respectively. Through this fitting procedure that entails 4 fit parameters for the two-point function fit and additional 4 parameters for the three-point function, we extract the desired matrix element, $\langle 0|\mathcal{O}|0\rangle$.

The type (b) fit we call simultaneous and it is a combined fit to Eqs. (18) and (19), using all 8 parameters. We expect that both fit types lead to very similar results. Proper analysis requires the evaluation of two- and three-point correlation functions on exactly the same gauge field configurations and for the same set of source positions. Only in such a case, the correlations among the fit parameters can be probed in a fully consistent manner. In practice, this is a severe problem for the observables we are interested in, which is due to the exponential increase of the noise with increasing source-sink time separation. When the time separation grows by one lattice spacing, the necessary statistics to suppress noise to the same level as before increases by a factor 2-3 in our data. Thus, small source-sink time separations yield a precise signal already for $\mathcal{O}(10^{3})$ measurements, while correlators at $t_{s}{=}12a$ require $\mathcal{O}(10^{5})$ measurements to reach a similar precision. For a given amount of computational resources, one, therefore, has to make a choice between using the same statistics for all $t_{s}$ values and thus obtaining very precise data at small values of $t_{s}$ and considerably less precise ones for larger values of $t_{s}$, or using different statistics achieving similar precision for all values of $t_{s}$. The severe drawback of the former choice is that the two-state fits are dominated by the precise results at small $t_{s}$ that can lead to biased results for the ground state matrix element. This bias becomes more severe as the boost increases due to the denser spectrum prohibiting the robust identification of the effects of excited states.

III.3.1 Two-point correlator

Before we compare different analysis methods for the extraction of quasi-PDF matrix elements, we discuss the two-state fit of the two-point correlator at our largest momentum, using Eq. (18). We obtain the boosted nucleon ground state energy of $aE_{0}^{\textrm{2-state}}{=}0.773(18)$, which is in line with a single-state fit, leading to a value $aE_{0}^{\textrm{1-state}}{=}0.788(9)$. These values are illustrated in the left panel of Fig. 8. As we showed in Fig. 4, the obtained ground state energy satisfies the continuum dispersion relation. The energy difference between the ground state and the first excited state is $a\Delta E{=}0.36(7)$, which in physical units is 0.76(15) GeV. This can be compared with the expectation from the approximation of noninteracting stable hadrons in a box, which leads, at our volume and the physical pion mass ($m_{\pi}L\approx 3$), to the following splittings (for $N\pi\pi$, $N\pi$, $N\pi\pi\pi\pi$, $N\pi\pi\pi$) with respect to the ground state ($N$): $\approx$ 0.27 GeV, 0.36 GeV, 0.54 GeV, 0.63 GeV, respectively. The above values hold for a nucleon at rest, whereas the spectrum becomes denser with increasing nucleon momentum. In the interacting case, the spectrum obviously changes, however its general features are unchanged; for a comprehensive discussion on the spectrum of excitations, we refer to the recent review of Ref. Green:2018vxw .

The energy extracted for the first excited state does not match the energy expected for the two- and three-particle states lying below. Furthermore, it is known that a two-state fit will model, in one exponential, all the states above the ground state that have an overlap with the interpolating field. Thus, one expects a bias in the determination of the first excited state. Consequently, the conclusions from two-state fits need to be checked against the plateau values or three-state fits, if the accuracy of the data allows for this. If consistency is found, this indicates ground state dominance up to the reached statistical precision.

III.3.2 Matrix elements of the unpolarized PDF

We first present the analysis of excited states in the nucleon matrix elements for the unpolarized operator using the $\gamma^{0}$ structure and the largest momentum, $P_{3}{=}10\pi/L$, where excited states effects are expected to be most severe. We use four values of the source-sink time separation, namely $t_{s}/a{=}8,9,10,12$ or in physical units $t_{s}{\simeq}0.75,0.84,0.94,1.13$ fm. We opt to increase the number of measurements as we increase $t_{s}$, to have statistical errors that are approximately the same for consequent values of $t_{s}$. This enables us to perform a reliable analysis at each value of $t_{s}$. We use the same CAA setup for all separations, as explained above. In Table 3, we collect the statistics for each value of $t_{s}$.

The plateau method is employed to analyze the data from each $t_{s}$ value and the results are shown in Fig. 9. For $t_{s}{=}8a$, the real part of the matrix element shifts towards smaller values at each $z/a$ compared to larger $t_{s}$, while the effect in the imaginary part is less prominent. Within statistical uncertainties, a convergence between the results at $t_{s}{=}10a$ and $t_{s}{=}12a$ is observed in both the real and imaginary parts. Comparing the data for the four $t_{s}$ values, we observe signs of a nonmonotonic behavior that affects the real and imaginary parts differently, depending on the value of $z$. This can introduce a complicated effect in the determination of the PDF. Ideally, one must compute the matrix elements for several large enough values of $t_{s}$ with equally small statistical errors and demonstrate convergence to one value. However, the exponential increase of the noise-to-signal ratio seen in Fig. 2 and the need for new sets of sequential inversions for each $t_{s}$, require computational resources beyond what is currently available, placing limitations on the maximum value of $t_{s}$. Nevertheless, if consistency can be found between the converged plateau values of the single-state fits and other methods, ground state dominance can be reliably established.

We apply the summation method for two sets of $t_{s}$ values, namely: (a) $t_{s}{=}8a,\,9a,\,10a,\,12a$ and (b) $t_{s}{=}9a,\,10a,\,12a$, with the results shown in Fig. 10. The first set yields results that, although consistent with those from the second fit, are systematically higher. Furthermore, for large values of $z/a$, the slope obtained using the (b) set is significantly larger than from the (a) set and $\chi^{2}/{\rm d.o.f.}>2$ for $z/a{\approx}10{-}15$, indicating that the (a) fits do not provide good description of the data. This is a consequence of the ratio at $t_{s}{=}8a$ being considerably below the ones at other source-sink separations. Fig. 11 shows examples of summation fits for our data. In the right panel, only the (b) fit yields an acceptable $\chi^{2}/{\rm d.o.f.}$ of around 0.16, while the full fit to all $t_{s}$ values has $\chi^{2}/{\rm d.o.f.}\approx 2.4$. Thus, the fitting ansatz of Eq. (17) does not provide a correct description of the data and the large-$z$ values in Fig. 10 for $t_{s}{=}[8a{-}12a]$ are not reliable. The effect is visible also for small $z/a$ (see the left panel of Fig. 11 for $z/a{=}0$ fits) and may result in overestimating the real part of the matrix elements. At $z/a{=}0$, the extracted value of the matrix element is 1.21(16) for the (a) set and 1.02(24) for the (b), after applying in both cases the renormalization factor $Z_{V}{=}0.7565$ Alexandrou:2015sea . We note that the matrix element should be equal to 1 upon renormalization, which is satisfied only by the set (b), while the value from the summation method that includes the smallest $t_{s}$ in the fit is considerably larger. All the above lead to the conclusion that using too small source-sink separations gives incorrect results. We will thus omit $t_{s}{=}8a$ in our summation method analysis, and take fits (b) as our final estimates from the summation method. All of these fits, both in the real and imaginary part, have $\chi^{2}/{\rm d.o.f.}{\lesssim}1$, indicating that the exponential contributions from higher excited states in Eq. (17) are sufficiently suppressed, yielding a good estimate of the ground state matrix element ${\cal M}$. However, the precision of the summation method estimates is much worse than the one from the plateau or two-state fits (see below).

We present now the two-state fits and compare results from the sequential and simultaneous two-state fits. In addition, we check the robustness of these fits by using three or four values of $t_{s}$. The fits using only two $t_{s}$ lead to considerably larger errors (with compatible central values) and hence, we do not show them. The different choices of fits with three $t_{s}$ values lead to very similar results in terms of the central values and their errors. Thus, in Fig. 12, we present results from one choice of three source-sink separations, namely $t_{s}{=}8a,\,9a,\,10a$, and from the fit to all $t_{s}{=}8a,\,9a,\,10a,\,12a$. We conclude that the simultaneous and sequential fits lead to statistically identical results for both sets of $t_{s}$ values. Including the results for the largest $t_{s}$ yields consistent fits, but with errors that are slightly smaller (up to 10-15%) for most $z$ values. This is visible particularly for the simultaneous fits at all $z$ and for sequential fits at small $z$ values in the real part. Note that the central values are essentially the same, suggesting that excited states are sufficiently suppressed.

In Fig. 13, we collect the results of the analyses using the three aforementioned procedures, for all $z/a$ values, and in Fig. 14 we present part of the same data for two selected values of the Wilson line length, $z/a{=}5$ and $z/a{=}10$, for better visibility. Excited states effects are clearly visible for $t_{s}{=}8a$, especially in the real part when using the plateau method. For $t_{s}{=}9a$ and $10a$, we observe small tension in the imaginary part – the plateau values are consistently lower for $t_{s}{=}9a$ and consistently higher for $t_{s}{=}10a$, with respect to the two-state fits. We note that, in particular, there is tension between the two-state fit using $t_{s}{=}8a{-}10a$ and the plateau values for $t_{s}{=}10a$, at $z/a{\approx}10$, which fails to satisfy our criterion of consistency between two-state fits and plateau fits. Hence, data at source-sink separations below $12a$ are likely still contaminated by excited states, although statistical fluctuation as the source of the trend cannot be excluded. We find that the results extracted from the plateau method using $t_{s}{=}10a$ and $t_{s}{=}12a$ are in agreement with each other and the ones from $t_{s}{=}12a$ are compatible with those extracted using the all-$t_{s}$ two-state fits for all values of $z/a$ and for both the real and the imaginary parts. The errors in the summation method are too large to draw any meaningful conclusions. Given that this investigation is carried out for the largest value of the momentum, we can take as our final value for the nucleon matrix elements of the unpolarized operator the data at $t_{s}{=}12a$ for all the momentum values. We note that increasing the nucleon boost would need a similarly thorough re-analysis on the effects of the excited states. Likewise, increased statistical precision could reveal excited states contamination even at this nucleon momentum. Thus, we emphasize that the attained conclusion about ground state dominance is valid only within the present statistical uncertainties of $\mathcal{O}(10)\%$.

III.3.3 Matrix elements of the helicity and transversity PDFs

In this subsection, we discuss the excited states effects for the two other Dirac structures used in our work – axial and tensor, associated with the helicity and transversity PDFs, respectively. We perform a similar analysis as for the unpolarized case at the largest momentum, $P_{3}{=}10\pi/L$, and using again four values of the source-sink separation, namely $t_{s}/a{=}8,9,10,12$, or in physical units $t_{s}\simeq\{0.75,0.84,0.94,1.13\}$ fm. The numbers of measurements are listed in Table 4. The methodology of our investigation is the same as for the unpolarized case, i.e. for each $t_{s}$ value, we perform single-state fits within a plateau region, and we combine data at different $t_{s}$ using the summation method approach, as well as two-state sequential and simultaneous fits.

In Figs. 15 and 16, we plot the full $z$-dependence of bare matrix elements for helicity and transversity PDFs, respectively. Fig. 17 displays a zoom for two selected values of $z/a$, i.e. $z/a{=}5$ and $z/a{=}10$. In this case, the fits using the summation method have good quality ($\chi^{2}/{\rm d.o.f.}{\lesssim}1$) even when including $t_{s}{=}8a$ and hence, we use all source-sink separations for this approach. For the real part, the two-state fits and the summation method yield results compatible with plateau fits for $t_{s}{=}12a$. The only observed tensions in the real part are rather small, possibly statistical fluctuations, and are visible for the plateau fits using $t_{s}/a{=}8,9$ versus the summation method, for intermediate and large $z/a$ in the helicity case. However, for the imaginary part, we consistently observe a rather strong dependence of the plateau values extracted when using $t_{s}{=}12a$ as compared to those from $t_{s}/a{=}8,9,10$, with 2-$\sigma$ to 3-$\sigma$ tensions. The latter are also incompatible with the results extracted using two-state fits to all $t_{s}$ and the summation method, both of which are consistent with the plateau results for $t_{s}{=}12a$. Moreover, the two-state fits using $t_{s}/a{=}8,9,10$ are incompatible with plateau fits for $t_{s}{=}10a$ (see right panel of Fig. 15), which again violates our criterion for consistency, necessitating the use of $t_{s}{=}12a$. This behavior reinforces the conclusion reached in the case of the bare matrix elements for the unpolarized quasi-PDF, namely that the tensions observed between results from source-sink separations below 1 fm are indeed manifestations of excited states contamination and ground state dominance is achieved, to around 10% statistical accuracy, only at $t_{s}{=}12a$. Hence, to extract PDFs, we take the plateau results for $t_{s}{=}12a$ as our preferred ones, since they are more precisely determined as compared to both the results extracted using the two-state and the summation approaches, but show full consistency with them. The more severe excited states effects observed in the cases of helicity and transversity are in accordance with observations connected to the extraction of the nucleon axial and tensor charges, where excited states contamination is more significant as compared to the case of the vector operator. Thus, studies that aim at a higher precision in the determination of quasi-PDF at increasing values of the nucleon boost would need large statistics to account for the increased errors resulting from having a nucleon state with large boost, but also for investigating excited states.

III.3.4 Remarks and computational costs

This concludes our investigation of excited states effects. We emphasize that the spectrum of nucleon excitations is rich, particularly for a boosted nucleon with quarks of physical masses. Thus, one method of extracting bare matrix elements can be misleading, as the fitted energy gap between the ground state and the explicitly modeled first excited state suggests there are tens of excited states. In such a situation, excited states need to be suppressed by going to large enough source-sink separations and robust statements can be made only when they are based on compatible results between different methods – in our case the plateau fits at the largest source-sink separation of around 1.1 fm, the two-state fits and the summation method. We note that the largest $t_{s}{=}12a$ is crucially needed to establish this compatibility. Having only $t_{s}{=}8a,9a,10a$ in the two-state fits and comparing to plateau fits at $t_{s}{=}10a$, one would conclude significant tensions between the former and the latter. This is best illustrated in the imaginary part of bare matrix elements for the helicity quasi-PDF, see the right panels of Fig. 15.

Finally, we would like to make some remarks on the computational resources that go well beyond what one needs for typical hadron structure calculations. The quasi-PDFs will reproduce the light-cone PDFs in the limit of large boosts. How large the boost should be, needs further investigation. Within the lattice QCD formulation, as already explained, one cannot increase the momentum of the nucleon to arbitrarily large values. The reasons are:

As the momentum increases, the signal-to-noise ratio rapidly deteriorates, despite the utilization of special methods, such as smearing techniques that reduce the noise. We find that to increase the momentum from $6\pi/L\approx 0.83$ GeV to $8\pi/L\approx 1.11$ GeV, we need to increase the statistics by a factor of around 3.7 (for unpolarized PDFs) and around 6 (helicity and transversity PDFs) and to increase to $10\pi/L\approx 1.38$ GeV by an additional factor of 3.7 or 6, respectively, in order to keep the statistical error approximately the same at $t_{s}{=}12a$. 2. 2.

A careful study of excited states must be carried out and becomes increasingly more difficult as the momentum increases, because of the denser spectrum. As the source-sink separation increases, the statistical errors grow exponentially, e.g. for $t_{s}{=}12a$, we needed a factor of about 10-12 (similar for all Dirac structures) more statistics for the same error as for $t_{s}{=}10a$ (at the largest boost). It is imperative to have large enough source-sink separations for at least three values of $t_{s}$ with comparable errors to perform a reliable analysis of excited states effects and extract the ground state matrix element.

Therefore, in order to reach a nucleon momentum of e.g. 2 GeV at $t_{s}{=}12a$, for unpolarized PDFs, we estimate that one would need $\mathcal{O}(100)$ million corehours (Mch) on a typical supercomputer, as compared to $\mathcal{O}(5)$ Mch at $P_{3}\approx 1.38$ GeV studied in this work and for momentum 3 GeV, we would need $\mathcal{O}(10000)$ MCh. For the polarized PDFs, the projected estimate for $P_{3}\approx 3$ GeV reads $\mathcal{O}(10^{6})$ Mch. In addition, one may have to increase the source-sink separation to account for the increased excited states contamination. If $t_{s}{=}14a$ is used, then the computational resources for a boost of 3 GeV would be $\mathcal{O}(10^{5})$/$\mathcal{O}(10^{7})$ million corehours for the unpolarized/polarized case, which is prohibitively expensive, given the computers available presently. These requirements may be alleviated with possible development of better algorithms, enhancing the signal-to-noise ratio at large nucleon boosts and large source-sink separations. Another way to milden the need for huge computational resources is the derivation of two-loop matching and conversion factors, foreseen in the near future. In this way, quasi-PDFs may be robustly connected to light-cone PDFs already at lower momenta.

We emphasize that a way to go is not to increase the nucleon boost uncontrollably, relying on precise data only at low source-sink separations. Failing to keep the statistical errors approximately constant as one increases $t_{s}$ may introduce uncontrolled systematic errors, since the fits will be determined mostly by the more accurate data at small $t_{s}$. As we have shown, it is essential that all source-sink separations have approximately equal errors for a reliable extraction of the matrix elements and this is the criterion we adopt. Unlike other studies Chen:2018xof ; Lin:2018qky ; Liu:2018hxv , we do not rely solely on two-state fits of data from source-sink separations with widely varied errors, which are thus dominated by the precise data at the small values of $t_{s}$. Such an approach can be uncontrolled and lead to a systematic bias in the final results.

IV Renormalization

Renormalization is needed in order to relate the bare lattice QCD matrix elements to physical results removing ultra-violet divergences, as well as finite dependence on the lattice action 444Elimination of residual dependence on the lattice formulation requires continuum extrapolation.. In the case of quasi-PDFs, one needs to also eliminate additional divergences arising due to the utilization of operators with a finite Wilson line. Renormalization of Wilson loops has been addressed long time ago using dimensional regularization (DR) for smooth contours Dotsenko:1979wR , as well as for contours containing singular points Brandt:1981kf . Based on arguments valid to all orders in perturbation theory, it was demonstrated that smooth Wilson loops in DR are finite functions of the renormalized coupling, while the presence of cusps and self-intersections introduces logarithmically divergent multiplicative renormalization factors, referred to as $Z$-factors. More importantly, it was shown that other regularization schemes are expected to lead to further $Z$-factors, which are power-law divergent with respect to the dimensionful ultraviolet cutoff. This also appears in the lattice formulation, where a divergence arises as a function of the lattice spacing, that increases exponentially with the length of the Wilson line as ${\sim}\,e^{{z}/{a}}$. Such a divergence must be removed prior to the extrapolation to the continuum limit. We describe in this section our renormalization program that includes the removal of the exponential divergence. A recent work on the perturbative renormalization of Wilson-line fermion operators of the type given in Eq. (9) has identified the mixing pattern among non-local straight Wilson-line operators, and led to the development of an appropriate renormalization prescription for both the multiplicative renormalization and the mixing coefficient Constantinou:2017sej . The non-perturbative renormalization program that we developed is a generalization of the RI′-scheme Martinelli:1994ty , appropriate for operators incluing a Wilson line Alexandrou:2017huk  555For alternative approaches, using, e.g. an auxiliary field method Green:2017xeu , see the recent review of Ref. Cichy:2018mum .. The $Z$-factors are extracted by imposing the following conditions

[TABLE]

We use the general notation $Z_{\cal O}$ and consider ${\cal O}{=}V_{0},\,A,\,T$ corresponding to the unpolarized, helicity and transversity operators, respectively. $Z_{\cal O}$ and $Z_{q}$ are the renormalization functions of the operator and the quark field, respectively. Both $Z_{\cal O}$ and $Z_{q}$ are scheme and scale dependent, and are expected to have some dependence on the pion mass. Also, $Z_{\cal O}$ is a function of the length of the Wilson line, $z$. ${\cal V}$ is the amputated vertex function of the operator and $S$ the fermion propagator, while ${\cal V}^{{\rm Born}}$ and $S^{{\rm Born}}$ are the corresponding tree-level values. Note that this condition is applied independently for each value of $z$. The RI′ renormalization scale, $\mu_{0}$, is chosen to be democratic in the spatial directions, that is $a\mu_{0}{=}\frac{2\pi}{L_{s}}(n_{t}+\frac{\pi}{2},n,n,n)$, which minimizes the ratio $P4{\equiv}{\sum_{i}\mu_{i}^{4}}/{(\sum_{i}\mu_{i}^{2})^{2}}$ and suppresses discretization effects Constantinou:2010gr ; Alexandrou:2015sea .

For the calculation of renormalization functions, we employ the momentum source method Gockeler:1998ye ; Alexandrou:2015sea that has the advantage of yielding results of high statistical accuracy. This method requires new inversions for each momentum used, but significant reduction in the gauge noise is observed, which by far outweighs the additional computational cost. Data are produced using the three $N_{f}{=}2$ ensembles given in Table 5 that have a different value of the pion mass. The twisted mass parameters of the light quarks in the sea and valence sectors have been set equal (isospin limit and unitary setup). Even though the ensemble simulated with the smallest value of the pion mass has a larger volume, the $Z$-factors are obtained at the same values of $(a\mu_{0})^{2}$ and $P4$. In the current work, we improve our previous analysis on the $Z$-factors of one-derivative operators Alexandrou:2017huk by performing:

• $\bullet$

a chiral extrapolation using the three ensembles at different value of the pion mass;

• $\bullet$

a fit on the chiral data to eliminate residual dependence on the RI′ scale. We use several values of the RI′ scale that cover the range $(a\,\mu_{0})^{2}\in[1-4]$. An extensive study on the choice of the renormalization scale and the corresponding systematic uncertainties can be found in Ref. Alexandrou:2017huk .

IV.1 Pion mass dependence

The pion mass dependence for the $Z$-factors for operators with a finite Wilson line has never been explored, and is expected that small values of $z$ will have weak dependence on $m_{\pi}$, as observed in the case of local operators computed within the same setup. However, it is not known how the $Z$-factors will behave when the length of the Wilson line is large. To obtain the $Z$-factors in the chiral limit, we fit the data from the three ensembles at each value of $z$. This fit is applied to the real and imaginary parts independently. The data for $Z^{{\rm RI}^{\prime}}_{\cal O}$ are expected to have a quadratic dependence on the pion mass (equivalently linear with respect to the twisted mass parameter) and are fitted using

[TABLE]

The chirally extrapolated value is given by the fit parameter ${Z}^{{\rm RI}^{\prime}}_{{\cal O},0}(z,\mu_{0})$. Since the fitted data are obtained on different ensembles, we use the super-jackknife method (see, e.g., Ref. AliKhan:2001xoi ) to correctly calculate the statistical error. This method is applicable to both correlated and uncorrelated data.

In Fig. 18, we show the pion mass dependence of the real and imaginary parts of the $Z$-factors for all three operators. For a clear presentation, we focus on $z/a{=}1,\,5$ and we plot against $m_{\pi}^{2}$ for the scale $a\,\mu_{0}{=}2.5$. We find that the dependence is almost constant in $m_{\pi}^{2}$ and the chiral fit yields a slope consistent with zero for both the real and imaginary parts for small values of $z$ (see, e.g., $z/a{=}1$). As $z$ increases, a non-zero slope is observed, with the dependence being linear, as expected. The dashed line corresponds to the chiral fit of Eq. (22), and the open symbols are the extrapolated values ${Z}^{{\rm RI}^{\prime}}_{{\cal O},0}(z,\mu_{0})$.

IV.2 Volume effects

Another source of systematic uncertainty entering the determination of matrix elements is due to the finite lattice extent. Finite volume effects are expected to be suppressed as $\exp(-m_{\pi}L)$ and based on empirical studies, $m_{\pi}L\geq 4$ is considered sufficient in most applications. However, most lattice QCD studies deal with matrix elements of local operators, and as discussed in Ref. Briceno:2018lfj , the length of the Wilson line entering the operator may enhance finite volume effects. To date, this systematic uncertainty has not been investigated due to the absence of lattice QCD computations using ensembles with different volumes, keeping the rest of the parameters constant. Such a study requires significant computational resources and is one of our future goals. Volume effects in the renormalization functions can be easily examined due to the reduced computational resources compared to the matrix elements. Therefore, we compute the $Z$-factors using the $48^{3}{\times}96$ ensemble of Table 5 and a $64^{3}{\times}128$ ensemble with the same pion mass. In Fig. 19, we plot the ratios

[TABLE]

as a function of the length of the Wilson line, for all three types of operators. The ratios are defined at the pion mass $m_{\pi}{=}$130 MeV and a renormalization scale $(a\,\mu_{0})^{2}{=}2.5$. The additional subscript 64 and 48 indicates the lattice volume. We find that both the real and imaginary parts do not show a statistically significant dependence on the volume, as the ratios take a maximum value of 1.02 and 1.03, respectively, for $z/a$ up to 15, which is well within the range of interest. In addition, we find an almost linear increase of ${\cal R}_{\cal O}(z)$, while ${\cal I}_{\cal O}(z)$ has an oscillatory dependence on the volume.

Hence, we conclude that volume effects in the $Z$-factors are small and have little effect in the renormalized matrix elements. This partly results from the fact that the bare matrix elements go to zero in the large-$z$ region. In determining the final values for the $Z$-factors, we do not include the $64^{3}{\times}128$ ensemble, because the ratio $P4$ is large for most of the scales $(a\mu_{0})^{2}\in[1{-}4]$ compared to the smaller volumes, leading to contamination from finite-$a$ effects, as seen in Ref. Alexandrou:2015sea . However, for $(a\,\mu_{0})^{2}{=}2.5$ which was used in the comparison, $P4$ is the same as the one obtained from the smaller volumes, which allows one to isolate the volume effects. Finite-$a$ effects are in fact under investigation for this class of non-local operators in lattice perturbation theory MC_HP_artifacts .

IV.3 Conversion to the standard and the modified $\overline{\mbox{\rm MS}}$ scheme

In order to compare renormalized lattice QCD matrix elements with phenomenological results extracted from global analyses, the $Z$-factors must be converted to the same scheme and evolved to the same scale as those used in the phenomenological analyses. Traditionally, the chosen scheme is the $\overline{\rm MS}$ and the scale $\bar{\mu}$ is typically set to 2 GeV. The appropriate conversion factors for the non-local operators with a straight Wilson line are taken from Ref. Constantinou:2017sej , where a calculation was carried out to one-loop level in perturbation theory, using dimensional regularization. Technical complications related to the non-locality of the operators under study make it very hard to extend such a calculation to higher loops, as done for local operators, usually known to three and four loops. As a consequence, it is expected that the $Z$-factors will have a residual dependence on the initial RI′ scale $\mu_{0}$. Thus, one typically computes the $Z$-factors at several values of the RI′ scale, as set by Eq. (20), and then uses an appropriate conversion to bring each ${Z}^{{\rm RI}^{\prime}}_{{\cal O},0}(z,\mu_{0})$ to ${Z}^{\overline{\mathrm{MS}}}_{{\cal O},0}(z,2{\rm GeV})$. The remaining dependence on $\mu_{0}$ can be studied with a linear fit in $(a\,\mu_{0})^{2}$, which is the leading order of a more complicated scale dependence. The final values are obtained by taking $(a\,\mu_{0})^{2}\to 0$ using

[TABLE]

For ${Z}^{\overline{\mathrm{MS}}}_{\cal O}(z,\bar{\mu},\mu_{0})$, we use the chirally extrapolated $Z$-factors converted to $\overline{\mbox{\rm MS}}$ at the scale $\bar{\mu}{=}$2 GeV. For simplicity, in the notation we dropped the subscript “0” appearing in the fit of Eq. (22). The desired quantity is the fit parameter ${\cal Z}^{\overline{\mathrm{MS}}}_{{\cal O},0}(z,\bar{\mu})$.

Our renormalization program entails a new element, namely the use of a modified $\overline{\mbox{\rm MS}}$ scheme (${\rm M}\overline{\mbox{\rm MS}}$). The development of such a scheme was motivated by the fact that the existing matching formulae do not satisfy particle number conservation (e.g. Ref. Izubuchi:2018srq ). The matching using the ${\rm M}\overline{\mbox{\rm MS}}$ scheme was already presented in our recent work Alexandrou:2018pbm ; Alexandrou:2018eet . Here we complement the previous analyses by giving the appropriate conversion of the $Z$-factors to the ${\rm M}\overline{\mbox{\rm MS}}$ scheme, instead of $\overline{\mbox{\rm MS}}$. We find that the resulting modification is numerically very small, but moves the final values of the PDFs towards the phenomenological ones. Details on the extraction of the ${\rm M}\overline{\mbox{\rm MS}}$ are given in Sec. V. In a nutshell, an additional conversion factor is needed to bring ${\cal Z}^{\overline{\mathrm{MS}}}_{{\cal O},0}(z,\bar{\mu})$ to ${\cal Z}^{{\rm M}\overline{\mathrm{MS}}}_{{\cal O},0}(z,\bar{\mu})$ via

[TABLE]

which has been computed perturbatively in dimensional regularization to one-loop level and is presented in the next section. The expression for the conversion ${\cal C}_{\cal O}^{\overline{\mathrm{MS}},{\rm M\overline{MS}}}$ is different for each operator under study and its general form is given by:

[TABLE]

where $\mu_{F}$ is the factorization scale that is taken equal to the $\overline{\mbox{\rm MS}}$ scale, that is, $\mu_{F}{=}\bar{\mu}{=}2$ GeV. The expression also contains the special functions ${\rm Ci}$ (cosine integral), ${\rm Si}$ (sine integral) and ${\rm Ei}$ (exponential integral), as well as the sign function (${\rm sgn}$). The coefficients $e^{(i)}_{\cal O}$ depend on the operator and their numerical values are given in Table 6.

In Fig. 20, we show ${\cal Z}^{{\rm M}\overline{\mathrm{MS}}}_{{\cal O},0}(z,\bar{\mu})$ after multiplying ${\cal Z}^{\overline{\mathrm{MS}}}_{{\cal O},0}(z,\bar{\mu})$ with the conversion factor given in Eq. (25). The data are shown against $(a\,\mu_{0})^{2}$ to demonstrate the dependence of the $Z$-factors on the initial RI′ scale. As done for the previous figures, we choose representative values of the length of the Wilson line, namely $z/a{=}1$ and $z/a{=}5$. We find that the imaginary part has a stronger dependence on the $\mu_{0}$ value. Note, however, that the imaginary part is an order of magnitude smaller than the real part of the matrix element and thus this dependence is still suppressed in the total matrix element, especially in the region $(a\,\mu_{0})^{2}\in[1-2]$, which is an indication of non-perturbative effects. Such behavior is also observed in local operators for $(a\,\mu_{0})^{2}{<}1$. This tendency seems to affect the region $(a\,\mu_{0})^{2}{<}2$ for non-local operators with $z/a{\geq}5$. Since we are using perturbative expressions for the conversion to the ${\rm M}\overline{\mbox{\rm MS}}$ scheme, it is important to choose a region where perturbation theory is valid, and thus, we choose the values obtained from $(a\,\mu_{0})^{2}\in[2-4]$ for both the real and imaginary parts. There is a systematic uncertainty attached to the $Z$-factor due to the choice of the fit range, and here we used the ranges $(a\,\mu_{0})^{2}\in[1-3],[1-4],[2-3],[2-4]$ to estimate the systematic uncertainty. Even though the real and imaginary parts of the $Z$-factors do have mild dependence on the range, the final values are consistent. Therefore, we do not give any systematic uncertainty from the $(a\,\mu_{0})^{2}{\to}0$ extrapolation. In the same figure, we also show the improvement of the non-perturbative results when subtracting ${\cal O}(g^{2}\,a^{\infty})$ effects computed perturbatively in Ref. Alexandrou:2015sea on the same ensembles. They are obtained by replacing Eq. (20) with

[TABLE]

where the denominator has been modified by subtracting the ${\cal O}(g^{2}\,a^{\infty})$ artifacts in $Z_{q}$, denoted by $A^{\infty}_{q}(\mu_{0})$. As expected, the differences between the subtracted results of Eq. (27) and the unsubtracted results of Eq. (20) are small, and mostly affect the real part of the $Z$-factors. In addition, subtraction of the artifacts in $Z_{q}$ is not sufficient to eliminate the dependence on $a\,\mu_{0}$. For the later to be achieved, subtraction of the artifacts in ${\cal V}(z,p,m_{\pi})$ is crucial, as has been discussed in Ref. Alexandrou:2015sea for the local operators. Such a computation for the non-local operators is not only more complicated technically, but requires the addition of stout smearing in the links of the operator resulting in lengthy expressions. This calculation is under way and will be presented in a separate publication MC_HP_artifacts .

The final values for the $Z$-factors are extracted by extrapolating $(a\,\mu_{0})^{2}{\to}0$ linearly the ${\rm M}\overline{\mbox{\rm MS}}$ data in the region $(a\,\mu_{0})^{2}\in[2-4]$ (Eq. (24)). They are shown in Fig. 21 for the vector operator. As can be seen, the statistical errors of the $Z$-factors are much smaller as compared to those of the matrix elements, because of the momentum source method used in the calculation of the vertex functions. The real part increases significantly with $z$, which is due to the power-law divergence of the Wilson line. Note that for $z/a{>}10$, the bare matrix elements decay to zero, and thus, the renormalized matrix elements, given by the complex multiplication

[TABLE]

are also zero but with large statistical uncertainty originating from the large values of the real part of the $Z$-factors.

Both the matrix elements and the renormalization functions of a given operator share similar properties with respect to $z$ (symmetric real part and antisymmetric imaginary part), with the difference that the $Z$-factor decreases when a smearing technique is applied to the Wilson line, while the matrix element increases. Of course, the dependence in the smearing is non-monotonic due to the complex nature of the matrix elements and $Z$-factors, but the stout smearing dependence is expected to cancel out in the renormalized matrix elements. This is demonstrated in Fig. 22, where we compare the renormalized matrix elements for the helicity operator extracted using 0, 5, 10 and 15 stout smearing steps, in the $\overline{\mbox{\rm MS}}$ scheme at 2 GeV and for momentum $\frac{6\pi}{L}$. This confirms that the elimination of the power divergences is correctly realized via the renormalization program, yielding compatible results for the renormalized matrix elements for different stout iterations. We find that this hold also for the renormalized matrix elements of the unpolarized and transversity operators and therefore any stout step may be used without changing the final physical result. It is worth mentioning that the agreement is more striking upon the $(a\,\mu_{0})^{2}{\to}0$ extrapolation of the $Z$-factors.

V Matching to light-cone PDFs

V.1 Derivation of the matching formulae

In this section, we discuss, in detail, the matching procedure that relates the quasi-PDFs, renormalized in some scheme, to light-cone PDFs in the same or other renormalization scheme, typically chosen to be the $\overline{\mathrm{MS}}$ scheme. We derive new matching formulae that relate $\overline{\mathrm{MS}}$-renormalized quasi-PDFs to $\overline{\mathrm{MS}}$-renormalized light-cone PDFs and conserve the particle number. To satisfy this condition, we introduce a modification of the $\overline{\mathrm{MS}}$ scheme, i.e. the ${\rm M}\overline{\mathrm{MS}}$ scheme, which was already partially discussed in the previous section, since it requires also the modification of conversion of renormalization functions.

Quasi-PDFs can be obtained as a Fourier transform (FT) of renormalized matrix elements, $h_{\Gamma}(P_{3},z)$,

[TABLE]

To relate the quasi-PDFs $\tilde{q}(x,P_{3})$ to light-cone PDFs, one relies on a perturbative matching procedure Xiong:2013bka ; Ji:2014gla ; Wang:2017qyg ; Stewart:2017tvs ; Izubuchi:2018srq . To one-loop order, and in the Feynman gauge, one needs to compute self-energy corrections, which include the usual self-energy plus the virtual ”sail” and ”tadpole” diagrams, and the vertex corrections, with the corresponding the real ”sail” and ”tadpole” diagrams. We use operators with four Dirac structures, namely $\gamma^{0}$ and $\gamma^{3}$ for the unpolarized distribution, $\gamma^{3}\gamma^{5}$ for the helicity, and $\gamma^{3}\gamma^{j}$, $j=1,2$ for the transversity case.

Since the extraction of the matching formula follows a similar process for all four operators, we use the $\gamma^{0}$ structure as an example and we only present the final results for the other cases. As already mentioned, we take the nucleon momentum in the third direction, $P{=}(P^{0},0,0,P^{3})$. It is assumed that, before gluon emission, the quark momentum $p$ is collinear to the nucleon momentum, i.e. $p{=}(p^{0},0,0,p^{3})$. It also obeys the Dirac equation $\not{p}u(p){=}0$, and carries a fraction of momentum $y$ of the parent hadron. After gluon emission, the quark has momentum $\xi p_{3}{=}\xi yP_{3}$. When the bare quark distributions are dressed taking the $p_{3}{\rightarrow}\infty$ limit, we have the usual quark distributions in the infinite momentum frame or, equivalently, on the light cone. Defining $x{\equiv}\xi y$, the nonsinglet quark distribution at one-loop is given by

[TABLE]

where $\Pi_{\gamma^{0}}$ denotes the self-energy corrections, $\Gamma_{\gamma^{0}}$ is the vertex corrections, and $\Lambda$ stands for the infrared (IR) and UV regulators in any given scheme. Using DR to regulate both the IR and UV divergences, the resulting one-loop correction for $0<\xi<1$ in the $\overline{\mathrm{MS}}$ scheme is

[TABLE]

where the pole from the UV divergence has already been subtracted, and $\bar{\mu}$ is the corresponding renormalization scale in $\overline{\mathrm{MS}}$. The pole from the soft IR, ${1}/{\epsilon_{IR}}$, can be absorbed in the bare distribution, at the factorization scale $\mu_{F}$, but it is explicitly written in Eq. (31), as it must cancel a similar pole arising in the one-loop correction to the quasi-PDF. The fact that the poles in $1/\epsilon_{IR}$ are the same for the light-cone PDF and the quasi-PDF, is a crucial observation that allows one to match them using perturbation theory. There are remaining IR divergences, which are located at $\xi{=}1$ and have their origin in the emission of soft gluons. However, they cancel between the vertex and self-energy one-loop corrections, which are related by $\Pi_{\gamma^{0}}(\Lambda){=}-\int_{0}^{1}d\xi\,\Gamma_{\gamma^{0}}(\xi,\Lambda)$.

To obtain the quasi-PDF, the quark is dressed with finite $p_{3}$. To simplify notation, we write the one-loop correction for positive $x$ only, that is

[TABLE]

Because $\tilde{\Gamma}_{\gamma^{0}}\neq 0$ for $\xi>1$, the lower limit of integration in Eq. (32) goes to zero, and the quasi-PDF has support for $x>1$. Using DR to perform the integrals, the one-loop correction to the vertex is given by:

[TABLE]

Note that the vertex corrections have no poles when DR is used, but they are non-zero outside the physical region $0{<}\xi{<}1$, a behavior that is different from the one-loop correction to the light-cone PDF. Thus, in addition to the IR divergences at $\xi{=}1$, there are UV divergences associated with the $\xi\rightarrow\pm\infty$ limits when Eq. (33) is integrated. So far, these UV divergences have not been subtracted, and this is reflected by the $\Lambda$-dependence in the functional form of $\tilde{\Pi}_{\gamma^{0}}(\Lambda,p_{3}/\mu_{F})$. DR is again used to isolate the poles, in which case $d\xi\rightarrow d\xi(\xi p_{3})^{d-1}(\bar{\mu}^{2}e^{\gamma_{E}}/4\pi)^{1/2-d/2}$. For the region $\xi{>}1$, the self-energy is:

[TABLE]

where $d{=}1{-}2\epsilon$, with $\epsilon{>}0$. Using the plus prescription for the collinear divergence at $\xi{=}1$, Eq. (34) can be written as

[TABLE]

Similar computations can be done for the other regions, and the renormalized one-loop self-energy in the $\overline{\mathrm{MS}}$ scheme for the quasi-PDF is given by

[TABLE]

with the corresponding renormalization function, $\tilde{Z}_{\Pi}^{\overline{\mathrm{MS}}}{=}1+\frac{\alpha_{s}}{2\pi}C_{F}\frac{3}{2}\left(\frac{1}{\epsilon}+\gamma_{E}-\ln\left({4\pi}\right)\right)$. Upon integration of Eq. (33) and employing the $\overline{\mathrm{MS}}$ scheme, the Ward identity of the resulting local current is naturally respected.

For the computation of the $x$-dependence of the distributions, however, the convolution involving the vertex correction prevents the aforementioned cancellation to occur, and one needs to impose a prescription to ensure conservation of the quark number for the nonsinglet distributions. From the equations for the quark distributions, Eq. (30), and quasi-distributions, Eq. (32), and including the negative $x$ regions we have to one-loop

[TABLE]

where $C_{\gamma^{0}}^{\overline{\mathrm{MS}}}(\xi)=\delta(1-\xi)-\alpha_{s}C_{F}[(\tilde{\Pi}_{\gamma^{0}}^{\overline{\mathrm{MS}}}-\Pi_{\gamma^{0}}^{\overline{\mathrm{MS}}})\delta(1-\xi)+\tilde{\Gamma}_{\gamma^{0}}\left(\xi\right)-\Gamma_{\gamma^{0}}^{\overline{\mathrm{MS}}}\left(\xi\right)]/2\pi$. The dependence on $p_{3}$, $\bar{\mu}$, and $\mu_{F}$ is implicit in $C^{\overline{\mathrm{MS}}}(\xi)$. $\tilde{\Gamma}(\xi)$ is the bare function, hence the absence of an $\overline{\mathrm{MS}}$ superscript. The poles of $1/\epsilon_{IR}$, on the other hand, automatically cancel in $C_{\gamma^{0}}^{\overline{\mathrm{MS}}}$, which is explicitly written as

[TABLE]

For the rest of this paper, we will employ the usual choice of renormalization and factorization scales, $\mu_{F}{=}\bar{\mu}$. Eq. (V.1) has an imbalance when integrated, because $\tilde{\Gamma}_{\gamma^{0}}$ outside the physical region picks up a logarithmic divergence in the momentum fraction as $\xi\rightarrow\pm\infty$, while this divergence has already been removed from the self-energy part in Eq. (V.1). This implies an anomalous non-conservation of the quark number, even if the classical current is conserved. In Ref. Izubuchi:2018srq , $C_{\gamma^{0}}^{\overline{\mathrm{MS}}}$ has also been computed and a plus prescription at infinity was used to cancel the UV divergences between the vertex and self-energy corrections. This explains the source of the difference between our result in Eq. (V.1) and that in Eq. (68) of Ref. Izubuchi:2018srq . Namely, in our case, the term proportional to the Dirac $\delta$-function depends on the ratio of the factorization and renormalization scales, while in the case of Ref. Izubuchi:2018srq , it depends on the quark momentum $p_{3}$. Both prescriptions, however, do not conserve quark number.

To treat the unbalanced divergence, we introduce a modified $\overline{\mathrm{MS}}$ scheme (${\rm M}\overline{\mathrm{MS}}$), which has already been discussed in the previous section and used in our previous work Alexandrou:2018pbm ; Alexandrou:2018eet . In this scheme, an extra subtraction is made outside the physical region of the unintegrated vertex corrections, which, in practice, renormalizes the $\xi$-dependence for $\xi>1$ and $\xi<0$ and removes the potential divergences,

[TABLE]

One can write Eq. (41) in $z$-space, noticing that $\xi$ outside the physical region is now renormalized, at the scale $\bar{\mu}$. The inverse Fourier transform is then from $\xi$ to $z\bar{\mu}$, and the extra subtraction in $z$-space is written as

[TABLE]

For consistency, Eq. (42) has been also applied to the renormalization functions to bring them to the ${\rm M}\overline{\mbox{\rm MS}}$ scheme, as described in Sec. IV.3. Thus, these renormalization functions are obtained as follows: the $Z$-factors calculated in the RI*′* scheme are converted to the $\overline{\mathrm{MS}}$ scheme according to the perturbative formulae of Ref. Constantinou:2017sej and then multiplied by the factor given in Eq. (26) to convert them to the ${\rm M}\overline{\mathrm{MS}}$ scheme. An important self-consistency check is that the expression of Eq. (42) must cancel the $z\rightarrow 0$ divergence in $\ln\left(z^{2}\right)$ present in the $\overline{\mathrm{MS}}$ scheme Constantinou:2017sej ; Ishikawa:2017faj . Indeed, in the limit $z\rightarrow 0$, one has:

[TABLE]

i.e. the $Z$-factor of the vertex corrections at $z{=}0$ is the same in our scheme and in the “ratio” scheme introduced in Ref. Izubuchi:2018srq and both cancel the divergence. The latter scheme was proposed as an alternative to our solution of the current conservation problem when using the pure $\overline{\mathrm{MS}}$ expression of Eq. (68) of Ref. Izubuchi:2018srq . The “ratio” scheme is another modification of the $\overline{\mathrm{MS}}$ scheme that also needs an additional conversion factor in order to bring the renormalization functions in the $\overline{\mathrm{MS}}$ scheme into the “ratio” scheme. The difference with respect to our solution is that the form of the matching kernel (or of $Z_{\Gamma_{\gamma^{0}}}^{\rm ratio}(z\bar{\mu})$ after doing the Fourier transform to $z$-space) implies that the $\xi$-dependence of the matching equation in the physical region is also renormalized, while our approach does not modify this region. Thus, the “ratio” scheme can induce potentially large modifications to the matched PDF, as we show below numerically.

We apply $\tilde{Z}_{\Gamma_{\gamma^{0}}}^{{\rm M}\overline{\mathrm{MS}}}(\xi)$ to Eq. (V.1) to obtain a matching kernel that keeps the norm of the nonsinglet distributions unchanged. We also include the results for the $\gamma^{3}$ and $\gamma^{3}\gamma^{5}$ Dirac structures. Because quarks are taken to be massless, the one-loop corrections are the same for the $\gamma^{3}$ and $\gamma^{3}\gamma^{5}$ cases. Compared to the $\gamma^{0}$ case, a shift of $+2(1-\xi)$ inside the plus prescription in the physical region of Eq. (V.1) is needed, and also, the factor of 5/2 in the last line of Eq. (V.1) becomes 7/2. The matching equation is then written as:

[TABLE]

with the matching kernels, for the different Dirac structures, given by

[TABLE]

where $\iota{=}0$ for $\gamma^{0}$, and $\iota{=}1$ for $\gamma^{3}$ and $\gamma^{3}\gamma^{5}$. We note that an alternative procedure to get $\overline{\mathrm{MS}}$-renormalized light-cone PDFs is to directly match to them from RI-renormalized quasi-PDFs. The relevant formulae for this were derived in Ref. Stewart:2017tvs for the operators with Dirac structures $\gamma^{3}$ and $\gamma^{3}\gamma^{5}$ and in Ref. Liu:2018uuj for $\gamma^{0}$. Such a procedure is equivalent, to one-loop order, to the one described above, but with different higher-order effects. This choice is investigated numerically below, by comparing the direct matching from a variant of the RI scheme to $\overline{\mathrm{MS}}$ with a two-step procedure that first brings the renormalization functions to the ${\rm M}\overline{\mathrm{MS}}$ scheme and then performs the ${\rm M}\overline{\mathrm{MS}}$ to $\overline{\mathrm{MS}}$ matching.

For the transversity case, the computation is similar, albeit simpler, because the usual one-loop vertex correction is zero in the Feynman gauge. The quark self-energy is not zero, but the unbalanced terms cancel each other in the matching equation, since they are the same for $\Pi_{\gamma^{3}\gamma^{j}}$ and $\tilde{\Pi}_{\gamma^{3}\gamma^{j}}$. Thus, $Z_{{\Gamma}_{\gamma^{3}\gamma^{j}}}^{{\rm M}\overline{\mathrm{MS}}}(z\bar{\mu})$ is given by the same functional form as Eq. (42), with the replacement of the numerical factors 3/2 and 5/2 by 2. $C_{\gamma^{3}\gamma^{j}}^{{\rm M}\overline{\mathrm{MS}}}$ is given by:

[TABLE]

The matching kernel for the transversity has been calculated previously, with a hard cut-off and a quark mass to regularize the UV and the IR divergences, respectively, but without renormalization Xiong:2013bka . More recently, following the idea of Ref. Stewart:2017tvs to match directly from RI-renormalized quasi-PDFs to light-cone PDFs in the $\overline{\mathrm{MS}}$ scheme, the matching for the transversity was also derived in Ref. Liu:2018hxv .

V.2 Comparison of matching in the $\overline{\mathrm{MS}}$, ${\rm M}\overline{\mathrm{MS}}$ and ratio schemes

Our choice for the matching to extract the light-cone PDFs from quasi-PDFs is to use a minimal modification of the $\overline{\mathrm{MS}}$ scheme that ensures current conservation. Thus, the matching formula is applied on a quasi-PDF renormalized in a modified $\overline{\mathrm{MS}}$ scheme, the ${\rm M}\overline{\mathrm{MS}}$ scheme, to yield the corresponding light-cone PDF in the pure $\overline{\mathrm{MS}}$ scheme. The difference between $\overline{\mathrm{MS}}$ and ${\rm M}\overline{\mathrm{MS}}$ schemes was also taken into account at the level of the renormalization functions in Sec. IV. Note that in our previous work Alexandrou:2018pbm ; Alexandrou:2018eet , we did not modify the conversion and, thus, our quasi-PDFs were renormalized in the standard $\overline{\mathrm{MS}}$ scheme. The formulae for the conversion modification, derived in the current paper, were not available at the time of our previous work and we argued that the numerical effect of the modification is subleading with respect to other uncertainties (as the modification in the $\overline{\mathrm{MS}}$ is minimal and only in the unphysical regions). Moreover, it disappears in the infinite momentum limit. Having now derived the relevant formulae to convert renormalization functions from the $\overline{\mathrm{MS}}$ scheme to ${\rm M}\overline{\mathrm{MS}}$, we test this and validate that it is indeed a small effect.

In Fig. 24, we compare matched PDFs obtained from quasi-PDFs renormalized in the $\overline{\mathrm{MS}}$ scheme, where the additional conversion of Eq. (26) is not included, and in the modified ${\rm M}\overline{\mathrm{MS}}$ scheme. As anticipated, the numerical effect is very small, particularly in the unpolarized case. Nevertheless, it is important to take the conversion modification into account to have a self-consistent procedure. We note that the modification brings the matched PDFs slightly towards phenomenological extractions, which is reassuring.

Another possibility to obtain a matching formulae with current conservation is to use the procedure proposed in Ref. Izubuchi:2018srq , the so-called “ratio” scheme, already mentioned in Sec. V.1. It consists in a different way of changing the $\overline{\mathrm{MS}}$ scheme to achieve current conservation, including a modification of the physical region. Thus, the effect on the matched PDFs is expected to be larger numerically than when using the ${\rm M}\overline{\mathrm{MS}}$ scheme. We show the comparison of the matching from our previous procedure and from the “ratio” scheme in Fig. 25. Indeed, the “ratio” scheme has a noticeable effect on the matched PDF, which is particularly large in the small-$x$ region. Theoretically, both the “ratio” and ${\rm M}\overline{\mathrm{MS}}$ schemes are valid modifications of the $\overline{\mathrm{MS}}$ scheme. Obviously, they are equivalent at the one-loop level, but they have significantly different higher-order effects. Hence, they can be both used to match quasi-PDFs to light-cone PDFs and the difference between them can serve as an estimate of truncation effects at one loop. For the remainder of the paper, we use the ${\rm M}\overline{\mathrm{MS}}$ scheme, as the one for which the higher-order effects are likely to be smaller, since this scheme is a milder modification of the $\overline{\mathrm{MS}}$ scheme, i.e. one that does not affect the physical region in the matching.

V.3 Comparison of $\overline{\mathrm{MS}}\rightarrow\overline{\mathrm{MS}}$ and RI$\rightarrow\overline{\mathrm{MS}}$ matching

As we discussed above, our renormalization and matching procedure proceeds in two steps. First, the renormalization functions computed in the RI*′* scheme are converted to the ${\rm M}\overline{\mathrm{MS}}$ scheme and evolved to a chosen scale, $\bar{\mu}{=}2$ GeV, using one-loop perturbative formulae derived in Ref. Constantinou:2017sej , and including the additional conversion from the $\overline{\mathrm{MS}}$ scheme to the ${\rm M}\overline{\mathrm{MS}}$ scheme, derived in this work. The resulting ${\rm M}\overline{\mathrm{MS}}$ $Z$-factors are used to renormalize the bare matrix elements, which are then Fourier-transformed from $z$-space to $x$-space, yielding the renormalized quasi-PDF in the ${\rm M}\overline{\mathrm{MS}}$ scheme. In the second step, a matching procedure is applied that brings the ${\rm M}\overline{\mathrm{MS}}$-renormalized quasi-PDF to the light-cone PDF in the $\overline{\mathrm{MS}}$ scheme, at the same renormalization scale $\bar{\mu}$. We will refer to this procedure as a two-step procedure.

In Ref. Stewart:2017tvs , a one-step procedure was proposed. It consists in applying the renormalization functions computed in the RI scheme to the bare matrix elements and taking the Fourier transform of RI-renormalized matrix elements to obtain a quasi-PDF renormalized in the RI scheme. The matching procedure then serves to simultaneously bring the quasi-PDFs to the light-cone PDFs and to convert from the RI scheme to the $\overline{\mathrm{MS}}$ scheme, evolving them from the given RI scale to the reference scale of 2 GeV. For this procedure, we use the formulae derived in Ref. Liu:2018uuj and consider the operator with the $\gamma_{0}$ Dirac structure. We choose the variant with the $\not{p}$ projection that differs from the projection that we used in our variant of the RI′ scheme. We refer to this choice of the projection as RI*／p*.

The comparison of the one-step and two-step procedures is presented in Fig. 26. In the left panel, we show our quasi-PDFs, renormalized in the RI*／p* scheme and in the ${\rm M}\overline{\mathrm{MS}}$ scheme. The corresponding renormalization functions are applied to the same lattice data for bare matrix elements. The results of application of the one-step and the two-step procedures are given in the right panel. We observe that the results of both procedures are fully compatible for small positive and negative $x$. For large positive $x$, the two-step procedure gives a PDF that goes to zero more slowly, while the PDF from the one-step procedure crosses zero at $x\approx 0.6$ and remains negative until $x\approx 1$ (reaching a minimum of around -0.16 at $x\approx 0.8$). In the large negative $x$ region, the situation is analogous – the one-step procedure leads to the PDF approaching zero from below and the two-step from above.

Both the one-step and the two-step procedures use a one-loop computation in continuum perturbation theory. Obviously, both contain higher-order contributions that cannot be quantified unless the two-loop formulae are available. Hence, neither of them can be considered to be the preferred one. It might happen that one of the procedures evinces smaller higher-loop effects without any theoretical arguments to support this, but it is not possible to say which one. The comparison between them is, thus, useful, as it reveals systematic effects due to higher-order terms, which may differ for different $x$ regions. The present study suggests that the large-$|x|$ regions suffer more from such higher-order effects. As Fig. 26 indicates, there may be significant two-loop effects in the conversion/evolution/matching procedures. Thus, a two-loop computation is indeed necessary. The two-step procedure allows for separating the conversion, the evolution and the matching and hence, a two-loop computation even for one of these, which is simpler than the full computation of the one-step procedure, can provide some insight.

V.4 Truncation of the Fourier transform

We turn now to the investigation of systematic effects related to the Fourier transform (FT) that is applied on the renormalized matrix elements in $z$-space, giving the quasi-PDFs in $x$-space. In principle, the FT integrates over all Wilson line lengths from zero to infinity, whereas lattice provides data for a finite number of discrete lengths $z/a$. Hence, one needs to decide about the maximum value of $z/a$ taken in the discretized FT integral. Ideally, it should be a value for which both the real part and the imaginary part of the matrix elements have decayed to zero. For unrenormalized matrix elements, the choice of such a value poses no problem. As can be seen in Fig. 27, the bare matrix elements are zero for any $|z/a|\gtrsim 15$. However, renormalized matrix elements involve multiplication of the bare ones by complex $Z$-factors that mix bare real and imaginary parts, according to

[TABLE]

[TABLE]

As a result, the real part of renormalized matrix elements can be non-zero even if the real part of bare matrix elements has already decayed to zero. This stems from an unphysical truncation effect in the perturbative conversion between the intermediate RI′ scheme $Z$-factors and the $\overline{\mathrm{MS}}$ ones. Analogously, the imaginary part of renormalized matrix elements gets contaminated by ${\rm Im}[Z]$ that multiplies the real part of bare matrix elements and leads to non-physical contributions. Moreover, the large values of $Z$-factors at large Wilson line lengths amplify the bare matrix elements and even if the latter are compatible with zero, this amplification introduces very large noise to the data, which propagates through the FT and matching to the final PDFs. These effects call for a truncation of renormalized matrix elements at some justified value of $z/a$, that we call $z_{\rm max}/a$. In Fig. 28, we illustrate the $z_{\rm max}/a$-dependence of the resulting final PDFs for all operators and for our largest momentum. The smallest (largest) value of $z_{\rm max}/a$ is chosen according to where the real (imaginary) part of the renormalized matrix elements is compatible with zero.

For the transversity case (lower plot of Fig. 28), we observe very little dependence on the choice of $z_{\rm max}/a$ (apart from increased statistical noise for larger values of $z_{\rm max}/a$) and we choose $z_{\rm max}/a{=}12$ as the value for our final plots. In this case, both the real and the imaginary parts of renormalized matrix elements are compatible with zero. Hence, the unphysical effects of truncating the Fourier transform at finite $z_{\rm max}/a$ are minimal. For the unpolarized case, a value of $z_{\rm max}/a$ where ${\rm Re}[h^{\rm ren}]$ and ${\rm Im}[h^{\rm ren}]$ are both zero does not exist. Hence, we observe that the variation of $z_{\rm max}/a$ leads to a visible effect in the upper left plot of Fig. 28. However, matched PDFs from all considered values of the cutoff are compatible with one another. In the end, we choose the middle value, $z_{\rm max}/a{=}10$, as the final one. In this way, the unphysical effect of truncation is split between the real and the imaginary parts. For the former, the matched PDF gets a contribution from unphysical negative values of ${\rm Re}[h^{\rm ren}]$ (resulting from the imaginary part of $Z$-factors being non-zero due to the truncation of the perturbative conversion between renormalization schemes). In turn, for the latter, part of the contribution from the imaginary part of renormalized matrix elements is missing. In the helicity case, ${\rm Re}[h^{\rm ren}]$ decays to zero at a similar value of $z/a$ as for the other cases and, contrary to the unpolarized matrix elements, does not go below zero. However, $|{\rm Im}[h^{\rm ren}]|$, reaches a minimum around $z/a\approx 12-13$ and then increases again, to become compatible with zero at $z/a\approx 17-18$, within huge errors. As our preferred value of the cutoff, we choose $z_{\rm max}/a{=}12$, where the real part is compatible with zero and the absolute value of the imaginary part is locally minimal. Varying $z_{\rm max}/a$ (upper right plot of Fig. 28) leads to changes in the matched helicity PDF that are statistically insignificant.

Apart from the non-physical effects appearing when the renormalized matrix elements are not compatible with zero at $z_{\rm max}/a$, further non-physical effects are introduced by the FT if the decay of $h^{\rm ren}$ is not fast enough. Namely, the periodicity of the Fourier transform leads to unphysical oscillations in the quasi-PDFs and finally in the matched PDFs. This effect is visible in all of our final distributions as a distorted approach of PDFs to zero for $x$ between 0.5 and 1, as a mild oscillatory behavior for large negative $x$, and as an unphysical minimum in the antiquark part at $-x\approx 0.1-0.2$. Decomposing the FT into cosine and sine transforms, one can observe that the quasi-PDF is bound to be negative if matrix elements entering the transforms do not decay to zero fast enough. The rate of decay of these matrix elements depends on the nucleon boost, as was demonstrated in Fig. 7. For larger momenta, the decay becomes faster, i.e. zero is reached at a smaller value of $z/a$, both for the real and the imaginary parts. Thus, the unphysical oscillations that we observe in the data may be attributed to insufficiently high momentum and therefore be considered higher-twist contaminations. In Fig. 29, we display the results for the unpolarized quasi-PDF for $z_{\rm max}{=}10a$, comparing bare, RI′ and ${\rm M}\overline{\mathrm{MS}}$ renormalized cases. As can be seen, the oscillations do not emerge when Fourier-transforming bare matrix elements that decay to zero faster. The amplification of the matrix elements by the $Z$-factors is responsible for their relatively slow decay. At this stage, one needs to remember that the $Z$-factors, as calculated now, are subject to two kinds of non-physical effects: lattice artifacts in their non-perturbative evaluation on a hypercubic lattice in the RI′ scheme, and truncation effects in the perturbative conversion to the ${\rm M}\overline{\mathrm{MS}}$ scheme and evolution to the reference scale of 2 GeV. After curing these effects, by computing higher-order conversion formula and by subtracting lattice artifacts computed in lattice perturbation theory MC_HP_artifacts , we expect that the oscillatory behavior will be considerably less prominent. This is supported also by the behavior observed for the quasi-PDF renormalized using the RI′ scheme, which has no conversion truncation effects, where we find that the oscillation is reduced with respect to the ${\rm M}\overline{\mathrm{MS}}$-renormalized quasi-PDF, as seen in Fig. 29. However, even in the RI′ scheme, the bare matrix elements are still amplified by the $Z$-factors, whose real part is exponentially increasing due to the presence of the power divergence related to the Wilson line. Thus, the oscillatory behavior in renormalized PDFs can only be avoided if the bare matrix elements decay faster, which can only happen for larger nucleon boosts.

An approach that may address non-physical oscillations was proposed in Ref. Lin:2017ani . The idea is to rewrite the FT that yields the quasi-PDF, using integration by parts as follows:

[TABLE]

where $h^{\prime}(z)$ is the derivative of the matrix elements with respect to the Wilson line length $z$. This step is exact, however, the proposed approximation is to neglect the surface term, even if it is non-vanishing. The latter happens when the matrix elements have not decayed to zero at $z{=}z_{\rm max}$. The resulting FT, given by the second term of Eq. (54), thus, involves the derivatives of renormalized matrix elements, which decay to zero faster than the matrix elements themselves. As argued above, a fast enough decay of the transformed matrix elements is sufficient to avoid oscillations. In this so-called “derivative” method, the oscillations are effectively transferred to the neglected surface term. We note, however, that this procedure is dangerous, due to the $1/x$ factor of the surface term, which leads to an uncontrolled effect for small values of $x$, where also the exponential of the surface term has a large real part. This aspect is illustrated in Fig. 30 for all three types of PDFs. The explicit $1/x$ factor amplifies the final PDF for small $|x|$, leading to its huge values, particularly in the small positive-$x$ region. This is very visible in the case of the unpolarized PDF, where the results are statistically the most precise. In the polarized cases, in the small-$|x|$ regions, one observes that the statistical errors are significantly enhanced by the “derivative” method. All three final PDFs evince unphysical divergences to $-\infty$ as $x\rightarrow 0$ in the antiquark part. Furthermore, the oscillatory behavior is very mildly affected, undermining the very aim of using the “derivative” procedure.

The difference between our results obtained using the standard FT and the “derivative” method, vanishes if $h(|z|\geq z_{\rm max}){=}0$. Only in this case, one can claim control over different regions of $x$. The problem associated with renormalized matrix elements not decaying fast enough to zero or, equivalently, having a higher-twist contaminations in the PDF, cannot be artificially hidden by using the “derivative” method. This is because the problem does not disappear, since the Fourier transform that defines the quasi-PDF is still truncated and part of the truncation encoded in the surface term may be hiding physical information relevant for reconstructing the full PDF. In other words, for a fully reliable extraction of PDFs, the surface term has to go to zero when applying actual lattice QCD data, and not by setting it to zero by hand. Another disadvantage of using the “derivative” method is that it introduces additional discretization effects, which need to be controlled. Moreover, as argued very recently in Ref. Karpie:2019eiq , the “derivative” method does not solve the problem with the Fourier transform, shown to be ill-conditioned if truncation is necessary and the data are limited to $\mathcal{O}(10)$ values of $z/a$. For all the above reasons, we opt not to use the “derivative” method and treat the residual oscillations as an indication that the nucleon boost needs to be increased in a controlled manner to extract reliable results on the PDFs. A possible way to improve this aspect for the currently attained nucleon momenta is advocated in Ref. Karpie:2019eiq , using advanced reconstruction methods, such as Backus-Gilbert or neural networks. It is an interesting direction for further study and we plan to explore these ideas in the future.

VI Final results - Discussion

In this section, we present our final results for the collinear non-singlet quark PDFs, whose determination was guided by the extensive investigation of several systematic effects presented in the previous sections. Although it is not possible to quantify all of the systematics, this investigation clearly points to the steps that one must follow in the future in order to arrive at precise extractions of PDFs from lattice QCD, with systematic uncertainties under control. We will be discussing these steps together with the presentation of the final results.

In Fig. 31, we show the quasi-PDFs, the matched PDFs and the final PDFs that take into account the nucleon mass corrections (NMCs) for the largest value of the momentum. The nucleon mass corrections are performed according to the formulae of Ref. Chen:2016utp . They are $\mathcal{O}(m_{N}^{2}/\Lambda_{\rm QCD}^{2})$ corrections to the PDF resulting from the finite mass of the nucleon and they can be derived in a closed form to all orders. We note that the one-loop matching that connects the quasi-PDFs to the light-cone PDFs is still relatively large for our largest nucleon boost. Increasing the nucleon boost while excited states contamination is fully under control is one aspect that has to be addressed in future computations. Another is the two-loop perturbative conversion between renormalization schemes and evolution to the chosen reference scale. This will further reduce the dependence on the initial RI*′* scale, and may be used to quantify truncation effects in the conversion factor. NMCs, on the other hand, are very small and negligible in most regions of $x$ and thus not expected to contribute to large systematic effects.

In Fig. 32 we show the dependence of the final PDFs on the nucleon boost. For the unpolarized PDF, we observe a strong effect when increasing the momentum from $6\pi/L$ to $8\pi/L$, which seems to converge after that yielding values that are compatible with those obtained for nucleon boost of $10\pi/L$ over almost the whole range of $x$. In addition, the oscillatory behavior becomes milder as the momentum increases. As we discussed above, this behavior results from the fact that for larger boosts, renormalized matrix elements become consistent with zero at smaller Wilson line lengths $z$, and affects the region $x{\approx}1$, as can be seen from the plots. For the helicity PDF, on the other hand, increasing the momentum from $6\pi/L$ to $8\pi/L$ is a minor effect, i.e. the final PDFs are compatible with each other for both of these momenta, over almost the whole range of $x$. However, increasing the momentum from $8\pi/L$ to $10\pi/L$, a significant shift is observed, especially for small values of $x$. Thus, it is likely that higher-twist effects are more pronounced for the helicity PDF. The oscillatory behavior is also milder in this case as the momentum increases, and the final PDF is compatible with zero at $x{=}1$. For the transversity PDF, the momentum dependence is the mildest and results using momentum $6\pi/L$ and $8\pi/L$ are compatible over a wide range of $x$-values. Similarly to the unpolarized and helicity PDFs, the unphysical oscillations are significantly reduced and the transversity PDF vanishes at $x{=}1$ for momentum $10\pi/L$.

We compare our final PDFs at the largest momentum with the phenomenologically extracted ones in Fig. 33. For the phenomenological results, we use the ones by CJ15 Accardi:2016qay , ABMP16 Alekhin:2017kpj and NNPDF3.1 Ball:2017nwa for the unpolarized PDF, DSSV08 deFlorian:2009vb , NNPDF1.1pol Nocera:2014gqa and JAM17 Ethier:2017zbq for the helicity PDF and two parameterizations for the PDFs transversity, one extracted from experimental SIDIS data and one where the tensor charge computed on lattice QCD was used as input in the phenomenological analysis Lin:2017stx . We stress that the comparison with phenomenological PDFs is intended to be qualitative, since we only include statistical uncertainties. The unpolarized PDF has a similar slope to the phenomenological PDFs, but lies above them in the positive $x$-region. Increasing the momentum should bring the PDF closer to the phenomenological values, as demonstrated in Fig. 32. Our results for the helicity PDF are compatible with phenomenological ones for $x\lesssim 0.4-0.5$. While this may indicate faster convergence of the quasi-PDF, it can also be the result of different systematic effects, with some of them possibly canceling out. Likewise, our result for the transversity PDF is in agreement with both the phenomenologically extracted one, as well as the one using the lattice determination of the tensor charge as input, for $x\lesssim 0.4-0.5$. An interesting feature is that the precision of our results is better than in both the phenomenological determinations. Thus, lattice QCD holds the prospect of impacting significantly our knowledge of the transversity PDF, in particular as we expect that in the next generation of lattice QCD computations of PDFs a better control of the systematic effects will be achieved.

The errors shown in the above figures are statistical, and significant effort is needed to properly quantify systematic uncertainties, present in the various steps of the analysis. Currently, these consist the dominating source of uncertainty, and addressing them will allow one to draw final conclusions from the lattice results. Based on the current studies, it is not possible to quantify and hierarchize them. Typical systematics related to the lattice calculation include discretization effects, finite volume effects and the role of the pion mass value. The pion mass of the ensemble used here is already at its physical value and hence, there is no need for a chiral extrapolation. The latter would be an important source of uncertainty, as the fit function is not known. We remark that, in practice, the effect of a non-physical pion mass can be large, as we demonstrated in Ref. Alexandrou:2018pbm , comparing the physical pion mass PDF with one at $m_{\pi}\approx 375$ MeV.

The parameters of the ensembles are expected to satisfy certain criteria for the range of values of the pion mass, the volume and the lattice spacing, to study uncertainties such as:

$\bullet$ Cutoff effects: A reliable control of cutoff effects requires at least three values of the lattice spacing smaller than 0.1 fm. Normally, cutoff effects are found to be relatively small in lattice hadron structure calculations. In the quasi-PDF computation, the nucleon is boosted to momenta for which $P_{3}$ becomes significant in comparison to the inverse lattice spacing and this may lead to increased cutoff effects. We note that for our largest momentum, we have $aP_{3}{=}0.65$ which is below the lattice cutoff (unlike Refs. Chen:2018xof ; Lin:2018qky ; Liu:2018hxv where the employed nucleon momenta are significantly above the lattice cutoff), and the continuum dispersion relation is still satisfied, as shown in Fig. 4. Still, it is unclear to what extent the good quality of the dispersion relation translates into discretization effects of the matrix elements considered here.

$\bullet$ Finite volume effects: Similarly to discretization effects, finite volume effects are also usually found to be rather small in hadron structure observables. The situation with quasi-PDFs is likely to be somewhat more complicated, since we use operators with Wilson lines of significant length. The volume behavior of such extended operators was considered by Briceño et al. Briceno:2018lfj within a model using current-current correlators in a scalar theory. Despite the fact that the model is not directly applicable to our investigation, it does provide a warning that the suppression of finite volume effects for matrix elements of spatially extended operators may change from the standard $\exp(-m_{\pi}L)$ to $(L^{m}/|L-z|^{n})\exp(-m_{N}(L-z))$, where $m$ and $n$ are undetermined exponents, which may become dominating for large $z$. Thus, finite volume effects may turn out to be a significant source of systematics and their investigation is crucial in the future.

$\bullet$ Systematic uncertainties in the determination of the renormalization functions: Uncertainties also arise in the computation of renormalization functions due to the breaking of rotational invariance. We have partly improved our work by subtracting lattice artifacts in the renormalization factor of the quark field, computed in lattice perturbation theory (see Sec. IV). A similar subtraction for the $Z$-factors of local operators was very successful Constantinou:2010gr ; Constantinou:2014fka ; Alexandrou:2015sea , and we intend to perform such a subtraction in the future MC_HP_artifacts . We note that pion mass dependence and finite volume effects in the renormalization functions studied here are found to be insignificant.

$\bullet$ Uncertainties specific to the quasi-PDF approach: The most significant uncertainty of this type is the requirement of having boosted nucleon states with values much larger than the nucleon mass. Although we observe convergence when increasing the boost, there are still undesirable effects from the finite value of the momentum that cause contamination from higher-twist effects (HTE), currently not known. HTE can be suppressed by going to larger boosts or alternatively explicitly computed and subtracted. Moreover, relatively small values of the momentum are likely to be responsible for unphysical oscillations, by making the matrix elements decay too slowly to zero, i.e. at rather large values of $z$. Faster decay of these matrix elements will naturally suppress the oscillations and the only way to achieve this reliably is to increase the momentum. Other ad hoc prescriptions, such as the “derivative” method, do not actually remove the source of the oscillations. In this technique, it is unclear how to handle the small-$|x|$ region, since the ignored surface term contains an explicit $1/x$ factor. In addition, the need for discretizing the derivative of the matrix elements introduces additional cutoff effects. Increasing the nucleon boost is difficult, as the signal is exponentially decaying at larger momenta and the latter also significantly increase excited states contributions, particularly when using simulations at the physical pion mass. As we argued, good control over excited states is necessary for extracting reliable physical results. Thus, a detailed study is required to maintain ground state dominance, following the approach explained in detail in this work. In this paper, we showed that excited states are suppressed below ${\sim}$10% if the source-sink separation is 1.1 fm at $P_{3}{\approx}1.4$ GeV, as established with a detailed analysis based on three methods. We note that at the physical pion mass and at small source-sink separations, the contamination comes from tens of excited states, and thus, a one-state or a two-state fit alone are not reliable. Compatibility between one-state and two-state fits provides a strong argument that ground state dominance has been achieved. Another possibly significant source of systematic errors are the perturbative truncation effects in the conversion to the $\overline{\mathrm{MS}}$ scheme, evolution to the reference scale and matching, all currently known to one-loop level. We observe that the effect of the latter is significant at our largest momentum and hence two-loop effects are likely to be substantial. Therefore, a two-loop computation would lead to better connection to light-cone PDFs with smaller values of the nucleon momentum.

VII Conclusions

This works presents a detailed investigation of the formalism employed to extract $x$-dependent PDFs within lattice QCD. The analysis was carried out using one ensemble of $N_{f}{=}2$ twisted mass fermions with quark mass fixed to its physical value. Results on the collinear unpolarized, helicity and transversity PDFs were obtained for the isovector flavor combination. The quasi-PDF method used here relies on a computation of equal-time correlators involving boosted nucleons, with the momentum increased to large enough values, so that the Large Momentum Effective Theory is applicable. In this work, we investigated three values of the momentum, reaching a maximum of 1.38 GeV. All necessary components to obtain light-cone PDFs have been considered, namely, calculation of bare nucleon matrix elements, renormalization, matching to light-cone PDFs and subtraction of finite nucleon mass corrections. The work presented here is an extension of our earlier work Alexandrou:2018pbm ; Alexandrou:2018eet , and includes details on technical and theoretical aspects, as well as improvements in various aspects of the calculation and examination of systematic effects. In particular, we provide the lattice techniques used in the evaluation of the matrix elements, such as momentum smearing and methods to eliminate (within statistical accuracy) excited states contributions. Theoretical developments associated with the non-perturbative renormalization are presented, where a chiral extrapolation and a modification of the conversion to the ${\rm M}\overline{\mathrm{MS}}$ scheme have been employed. The latter is a variant of the $\overline{\mathrm{MS}}$ scheme that ensures particle number conservation in the matching procedure. Issues related to the Fourier transform as well as different matching prescriptions were also explored.

Particular attention was given to discussing the role of systematic effects in the various steps of the analysis. Typical systematics related to the lattice calculation that are also common in other hadron structure calculations include discretization effects and finite volume effects. These can be addressed and eliminated by simulations using additional gauge field configuration ensembles. At present, these are part of the unquantified uncertainties, that need to be addressed in future studies. Another typical systematic enters lattice data if the pion mass of the ensemble is at an unphysical value. In this work, we use simulations at the physical point and chiral extrapolation is not needed. Analyses using non-physical values of the quark mass would introduce an important source of uncertainty, as the fit function to extrapolate the results to the physical point is not known.

Apart from the aforementioned systematics that are common in other hadron structure investigations, there are additional ones that are specific to the quasi-PDF approach. As discussed in the previous sections, increasing the nucleon boost can lead to further systematic uncertainties. High values of the momenta reduce the appearance of oscillations, but at the same time, increase the number of contributing excited states. Systematics related to the truncation of the conversion factor and the matching formula to one-loop level are non-negligible. In particular, having a matching formula to two-loops may lead to better convergence to light-cone PDFs at smaller nucleon momenta. Hence, a two-loop computation is strongly desired and would help to establish better connection to light-cone PDFs with smaller values of the nucleon momentum.

Despite these uncertainties, this work demonstrates a tremendous progress in the determination of PDFs from the quasi-distribution approach Cichy:2018mum . Lattice QCD results confirm the feasibility of extracting PDFs from first-principle calculations. The success in the quasi-PDF approach for the nucleon has also resulted in studies of other hadrons and alternative approaches, as summarized in Ref. Cichy:2018mum . The theoretical and technical aspects are now well understood for the nucleon studies and addressing the systematic uncertainties is the next step that will be enabled by using the methodology developed, in combination with the availability of increased computational resources. In fact, our preliminary results for an $N_{f}{=}2{+}1{+}1$ twisted mass ensemble with physical values of the quark masses, lattice spacing 0.082 fm and a larger volume of $64^{3}{\times}128$ have been shown in Ref. Alexandrou:2018own and demonstrate the future direction and progress of such computations. Furthermore, the production of another ensemble at a finer lattice spacing is already on its way.

Acknowledgements.

We would like to thank all members of ETMC for their constant and pleasant collaboration. We also thank Y. Zhao for useful discussions. This work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 642069 (HPC-LEAP). K.C. was supported by National Science Centre (Poland) grant SONATA BIS no. 2016/22/E/ST2/00013. F.S. was funded by DFG project number 392578569. M.C. 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. This research used computational resources of the Oak Ridge Leadership Computing Facility (OLCF), which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. Additional computational resources are provided by John-von-Neumann-Institut für Computing (NIC) for the JURECA supercomputer (allocation ECY00), Prometheus supercomputer at the Academic Computing Centre Cyfronet AGH in Cracow (grant ID quasipdfs), Eagle supercomputer at the Poznan Supercomputing and Networking Center (grant no. 346), Okeanos supercomputer at the Interdisciplinary Centre for Mathematical and Computational Modelling in Warsaw (grant IDs gb70-17, ga71-22).