Ratio of strange to $u/d$ momentum fraction in disconnected insertions

TL;DR

This study calculates the ratio of strange to light quark momentum fractions in disconnected insertions using lattice QCD, providing a nonperturbative constraint to improve global parton distribution function fits.

Contribution

It presents the first lattice QCD determination of the strange to light quark momentum fraction ratio in disconnected insertions with comprehensive renormalization and mixing effects.

Findings

Ratio of strange to light quark momentum fractions: 0.795(79)(77) at 2 GeV

Provides a new nonperturbative constraint for global PDF analyses

Connects lattice QCD results with phenomenological global fits

Abstract

The ratio of the strange quark momentum fraction $\langle x\rangle_{s+\bar{s}}$ to that of light quark $u$ or $d$ in disconnected insertions (DI) is calculated on the lattice with overlap fermions on four domain wall fermion ensembles. These ensembles cover three lattice spacings, three volumes and several pion masses including the physical one, from which a global fitting is carried out. A complete nonperturbative renormalization and the mixing between the quark and glue operators are taken into account. We find the ratio to be $\langle x\rangle_{s+\bar{s}}/\langle x\rangle_{u+\bar{u}} ({\rm DI})=0.795(79)(77)$ at $\mu = 2$ GeV in the $\overline{{\rm MS}}$ scheme. This ratio can be used as a constraint to better determine the strange parton distribution especially in the small $x$ region in the global fittings of PDFs when the connected and disconnected sea are fitted and evolved separately, demonstrating a new way that connects lattice calculations with global analyses.

Tables

Table 1. Table 1: Naming scheme in the path-integral formulation of the Euclidean hadronic tensor.

“connected”	connected insertions of the Green’s function
“disconnected”	disconnected insertions
“parton”	forward propagating quark between currents
“antiparton”	backward propagating quark
“valence”	quarks from the interpolating field operators
“sea”	quarks and antiquarks from gluons

Table 2. Table 2: Parameters of the RBC/UKQCD ensembles: spatial/temporal size, lattice spacing, sea strange quark mass under MS ¯ ¯ MS \overline{\textrm{MS}} scheme at 2 GeV, pion mass with degenerate light sea quark, and the number of configurations.

Symbol	$L^3\times T$	$a$(fm)	$m_s^s$(MeV)	$m_{\pi }^s$(MeV)	$N_{\mathrm{cfg}}$
32ID	${32}^3\times 64$	0.1431(7)	89.4	171	200
24I	${24}^3\times 64$	0.1105(3)	120	330	203
48I	${48}^3\times 96$	0.1141(2)	94.9	139	81
32I	${32}^3\times 64$	0.0828(3)	110	300	309

Table 3. Table 3: The details of the overlap simulation in the valence sector for the CI case, including the name of the lattice, the grid type of source 𝒢 src subscript 𝒢 src \mathcal{G}_{{\rm src}} (the notations such as 12-12-12 denote the intervals of the grid in the three spatial directions; see reference (Yang  et al. , 2016a ) for more details), the number of noises for the source grids N src subscript 𝑁 src N_{{\rm src}} , the time positions of sources t src subscript 𝑡 src t_{{\rm src}} , the grid type of sink 𝒢 sink subscript 𝒢 sink \mathcal{G}_{{\rm sink}} , the number of noises for the sink grids N sink subscript 𝑁 sink N_{{\rm sink}} , the source-sink separations ( t sink − t src ) subscript 𝑡 sink subscript 𝑡 src (t_{{\rm sink}}-t_{{\rm src}}) , and the bare valence quark masses m q v ​ a superscript subscript 𝑚 𝑞 𝑣 𝑎 m_{q}^{v}a .

Lattice	${𝒢}_{\mathrm{src}}$	$N_{\mathrm{src}}$	$t_{\mathrm{src}}$	${𝒢}_{\mathrm{sink}}$	$N_{\mathrm{sink}}$	$\left(t_{\mathrm{sink}}-t_{\mathrm{src}}\right)$	$m_q^{v}a$
					5	0.88 fm
24I	12-12-12	1	(0, 32)	2-2-2	5	1.11 fm	(0.0102, 0.0135, 0.0160, 0.0203)
					5	1.33 fm
					3	0.99 fm
32I	16-16-16	1	(0, 32)	1-1-1	3	1.16 fm	(0.00765, 0.00885, 0.0112, 0.0152)
					3	1.24 fm
					4	1.29 fm
32ID	16-16-16	6	(0, 32)	1-1-1	5	1.43 fm	(0.0042, 0.0060, 0.011, 0.014, 0.017, 0.022)
					12	1.57 fm
					4	0.88 fm
48I	12-12-12	5	(0, 32, 64)	1-1-1	8	1.11 fm	(0.0024, 0.0030, 0.00809, 0.0102, 0.0135, 0.0160, 0.0203)
					12	1.33 fm

Table 4. Table 4: The bare valence strange quark mass parameters and mass renormalization constants ( MS ¯ ¯ MS \overline{\rm MS} at 2 GeV) used in DI.

	32I	24I	48I	32ID
$m_{s}a$	0.04454	0.06347	0.06548	0.08500
$Z_m$	0.9467(57)	0.8872(68)	0.8872(68)	0.8094(26)

Table 5. Table 5: The nonperturbative renormalization constants on different ensembles at MS ¯ ¯ MS \overline{\textrm{MS}} 2 GeV . The 24I and 48I ensembles share the same renormalization constants due to the same lattice spacing.

Symbol	$Z_{QQ}^{\overline{\text{MS}}}$	$\delta Z_{QQ}^{\overline{\text{MS}}}$	$Z_{QG}^{\overline{\text{MS}}}$
32ID	1.25(0)(2)	0.018(2)(2)	0.017(17)
24I/48I	1.24(0)(2)	0.012(2)(2)	0.007(14)
32I	1.25(0)(2)	0.008(2)(2)	0.000(14)

Table 6. Table 6: Setups of the two-state fits in the DI part. The source-sink separations used in the fits, the number of points dropped on the source side and the sink side, and the prior value and width of δ ​ m 𝛿 𝑚 \delta m are listed for each lattice.

lattice	separations (a)	source drop	sink drop	prior $\delta ma$
32ID	6, 7, 8	1	1	0.4($\infty$)
24I	8, 9, 10, 11	2	1	0.3($\infty$)
48I	6, 7, 8	1	1	0.3(0.6)
32I	8, 9, 10, 11	2	3	0.2($\infty$)

Table 7. Table 7: Setups of the two-state fits in the CI case. The source-sink separations used in the fits, the number of points dropped on the source side and the sink side and the prior value and width of δ ​ m 𝛿 𝑚 \delta m are listed for each lattice.

lattice	separations (a)	source drop	sink drop	prior $\delta ma$
32ID/$u$	9, 10, 11	3	2	0.4(0.2)
24I/$d$	8, 10, 12	2	1	0.3(0.3)
48I/$d$	8, 10, 12	2	1	0.3(0.3)
32I/$d$	12, 14, 15	2	3	0.2(0.1)

Table 8. Table 8: The systematic error budget of the ratio ℛ ℛ {\cal R} .

source	absolute error	relative error
two-state fit (including the use of prior on $\delta m$)	0.050	6.3%
finite lattice spacing	0.003	0.4%
finite lattice volume	0.027	3.4%
pion mass extrapolation	0.045	5.7%
mixed action	0.016	2.0%
strange quark mass	0.008	1.0%
renormalization and mixing	0.016	2.0%
lack of mixing from charm quark	0.004	0.5%
all (combined in quadrature)	0.077	9.7%

Equations

$$\mathcal{R}\equiv \frac{⟨x{⟩}_{s+\bar{s}}}{⟨x{⟩}_{u+\bar{u}}(\mathrm{DI})},$$

$$⟨x{⟩}_{q,g}\equiv -\frac{⟨N|\frac{4}{3}{\overline{T}}_{44}^{q,g}|N⟩}{M_N⟨N|N⟩},$$

$$C_3^{q,g}(t_f,\tau )=\sum _{\vec{x},\vec{y}}⟨\chi (t_f,\vec{y}){\overline{T}}_{44}^{q,g}(\tau ,\vec{x})\bar{\chi }(0,\mathcal{G})⟩,$$

$$\begin{array}{rl}& {\Pi }^{q,g}(t_f,\tau )=\frac{\mathrm{Tr}\left[{\Gamma }_e{C}_3^{q,g}(t_f,\tau )\right]}{\mathrm{Tr}\left[{\Gamma }_e{C}_2(t_f)\right]}\\ =& ⟨N|{\overline{T}}_{44}^{q,g}|N⟩+c_1^{q,g}{e}^{-\delta m\left(t_f-\tau \right)}+c_2^{q,g}{e}^{-\delta m\tau }+c_3^{q,g}{e}^{-\delta m{t}_f}.\end{array}$$

$$⟨x{⟩}_{u,d,s}^{R,\mathrm{DI}}$$

$$+$$

$$\mathcal{R}$$

$${m^{\prime }}_{\pi ,vs}^2=\frac{1}{2}\left(m_{\pi ,vv}^2+m_{\pi ,ss}^2\right)+a^2{\Delta }_{\mathrm{mix}},$$

$$s(x)+\bar{s}(x)=\frac{1}{\mathcal{R}}\left(\bar{u}(x)+\bar{d}(x)\right)-c\left(\bar{u}(x)-\bar{d}(x)\right),$$

Full text

I Introduction

Understanding the structure of the nucleon in terms of quarks and gluons from QCD is one of the most challenging aspects of modern nuclear and particle physics Aprahamian et al. (2015) and is of great importance in learning about how the visible Universe is built. Parton distribution functions (PDFs), which describe the number density of a parton with a certain longitudinal momentum fraction $x$ and at a particular energy scale $Q^{2}$ inside a nucleon, reveal a lot of pertinent and essential information about the nucleon structure. In general, PDFs are determined by global analyses of deep inelastic scattering (DIS) and Drell-Yan experiments under the framework of QCD factorization theorems.

For the extensively studied unpolarized PDFs, recent attention is focused on the less-known flavor structure, which is believed to implicate the nonperturbative nature of the parton distributions due to confinement. A typical example is the strange parton distribution which is the most uncertain among the unpolarized PDFs. Three recent global fittings Ball et al. (2017); Dulat et al. (2016); Harland-Lang et al. (2015) with NNLO analysis show that $x(s(x)+\bar{s}(x))$ has large errors, $\sim 50\%$ or more at $x=10^{-3}$, and the central values of the three fits differ by $\sim 30\%$, at $Q^{2}=4\,{\rm GeV}^{2}$. On the other hand, as a nonperturbative approach of solving QCD from first principles, lattice QCD could also play an important role in the study of nucleon structure. Although there exist several pioneering approaches aiming to directly calculate the $x$ dependent PDFs on the lattice (e.g., (Liu and Dong, 1994; Liu et al., 1999; Liu, 2000; Ji, 2013; Radyushkin, 2017; Ma and Qiu, 2018)), it is maturer and more straightforward to calculate the moments of PDFs on the lattice, which provides constraints to the PDFs.

Recent lattice calculations can already determine several quantities, e.g. the strange quark magnetic moment Sufian et al. (2017) and strange quark spin contribution Liang et al. (2018), to a higher accuracy than experiments have done to date, but difficulties still exist in constraining the unpolarized strange parton distribution. The direct difficulty is that the lattice signals of the strange quark momentum fraction $\langle x\rangle_{s+\bar{s}}$ that involves only the disconnected insertions (DIs) are not good enough (Yang et al., 2018a) to provide strong constraint to the global fittings. On the other hand, lattice ratios of correlated quantities like

[TABLE]

where $\langle x\rangle_{u+\bar{u}}({\rm DI})$ stands for the light quark momentum fraction in disconnected insertions only (using $d$ quark makes no difference based on present lattice setup) usually have much smaller statistical uncertainty compared with the momentum fractions themselves due to the cancellation of the statistical fluctuations of the numerator and the denominator. However, the ratio ${\cal R}$ cannot be directly connected to the current global fittings without further theoretical insight, as the DI component on the lattice is part of the $u/d$ sea partons in global fittings.

The parton degrees of freedom have been rigorously defined in the path-integral formulation of the hadronic tensor and classified according to the topologically distinct connected insertions (CIs) and DIs Liu and Dong (1994); Liu et al. (1999); Liu (2000). This precise definition is natural for the lattice community and it is advocated to be adopted and accommodated in global fittings Liu et al. (2012); Liu (2017). Upon this basis, the ratio $\cal R$ has clear physical meaning and can be used as a strong constraint to better determine the strange quark distributions in future global fittings employing the path-integral classification. In this manuscript, we report a complete lattice calculation of the ratio $\cal R$ at three lattice spacings and several pion masses including the physical one. Nonperturbative renormalization and the mixing from the glue momentum fraction are considered.

II Theoretical background

The parton classification is revealed in the path-integral formulation of the Euclidean hadronic tensor $\langle p|\int\frac{d^{3}\vec{x}}{2\pi}e^{-i\vec{q}\cdot\vec{x}}J_{\mu}(\vec{x},t_{2})J_{\nu}(\vec{0},t_{1})|p\rangle$ Liu and Dong (1994); Liu et al. (1999); Liu (2000), where $|p\rangle$ is a nucleon state, $J_{\mu}$ and $J_{\nu}$ are two currents inserted at $t_{2}$ and $t_{1}$ and $\vec{q}$ is the momentum transfer. Three gauge invariant and topologically distinct path-integral diagrams of the 4-point functions of the Euclidean hadronic tensor, which entail leading twist contributions, are illustrated in Fig. 1. The solid lines represent quark propagators. The Minkowski hadronic tensor is the inverse Laplace transform of its Euclidean counterpart and is used to extract PDFs experimentally, thus these diagrams of insertions (Fig. 1a and 1b are CIs and 1c is DI) classify parton degrees of freedom. A complete naming scheme is listed in Table 1. Following this scheme, we denote Fig. 1a as $q^{v+cs}$ since in addition to the obvious valence contribution, there is also the connected-sea (CS) contribution coming from the higher Fock-state components in the Z-graph of quark lines between the two currents. Similarly, we also have the CS antipartons ($\bar{q}^{cs}$) in Fig. 1b and the disconnected-sea (DS) partons and antipartons (${q}^{ds}+\bar{q}^{ds}$) in Fig. 1c.

This is a general classification of partons based on the continuum path-integral formulation of QCD, which is applicable to the lattice. The nomenclature of CS and DS follows those in time-ordered perturbation theory sterman1993introduction; schweber2011introduction. These two sources of sea quarks have interesting flavor dependence. While $u$ and $d$ have both the CS and the DS, $s$ and $c$ have only the DS. Several experimental results demonstrate the necessity of this classification, e.g. the Gottfried sum rule violation is explained by the existence of CS Liu and Dong (1994). Also, under this classification, the net valence contribution is defined as $q^{v}\equiv q^{v+cs}-\bar{q}^{cs}$ which is not the same as the usual definition of $q^{v}\equiv q-\bar{q}=q^{v+cs}+q^{ds}-\bar{q}^{cs}-\bar{q}^{ds}$ if one does not assume ${q}^{ds}=\bar{q}^{ds}$. Actually, this definition of $q^{v}$ from QCD path-integral avoids some ambiguities like having a “valence” strange quark distribution when the NNLO evolution equations are involved which makes $s(x)\neq\bar{s}(x)$ Liu (2000).

Before we can study all these parton components explicitly on the lattice by directly calculating the hadronic tensor, this classification extends the way of using common lattice calculation of 3-point functions to study the CS and DS. It is shown Liu (2000) that, upon short distance expansion, Fig. 1a together with Fig. 1b becomes the CIs of 3-point functions in Fig. 2a for a series of local operators $\sum_{n}O_{q}^{n}$, from which the CI moments of PDFs are obtained. By the same token, the disconnected 4-point functions in Fig. 1c become the DIs of 3-point functions in Fig. 2b to obtain the DI moments. The DI ratio $\cal R$ therefore represents the DS ratio of the strange quark to light quark, containing the information needed in global fittings to separate CS and DS and improve the strange parton results.

III Numerical Details

We use overlap fermions (Neuberger, 1998) as valence quarks on four $2+1$-flavor RBC/UKQCD gauge ensembles with domain wall fermions (Aoki, 2009; Blum et al., 2016). The parameters of the ensembles are listed in Table 2. We have three different lattice spacings and lattice volumes respectively, and four values of sea pion mass with one at the physical point. For the valence sector, multiple partially-quenched valence quark masses are used, owing to the multi-mass algorithm. We choose four valence quark masses ranging from $\sim 250$ to $\sim 400$ MeV on the 24I and 32I ensembles and 7/6 quark masses in the range $\left[130,400\right]$ MeV on the 48I/32ID ensemble. Combining these ensembles and valence pion masses in a global analysis helps to control the lattice systematic uncertainties and leads to our final result at the physical limit.

The quark and glue momentum fractions in the nucleon can be defined by the matrix element of the traceless diagonal part of the energy-momentum tensor (EMT) in the rest frame (Horsley et al., 2012),

[TABLE]

with $\overline{T}_{44}^{q}=\int d^{3}x\overline{\psi}(x)\frac{1}{2}\left(\gamma_{4}\overleftrightarrow{D}_{4}-\frac{1}{4}{\displaystyle\sum_{i=0,1,2,3}}\gamma_{i}\overleftrightarrow{D}_{i}\right)\hat{\psi}(x)$ and $\overline{T}_{44}^{g}=\int d^{3}x\frac{1}{2}\left[E(x)^{2}-B(x)^{2}\right].$ Here $\hat{\psi}=(1-\frac{1}{2}D_{{\rm ov}})\psi$ is for giving rise to the effective quark propagator $\left(D_{c}+m\right)^{-1}$, where $D_{c}$ satisfying $\left\{D_{c},\gamma_{5}\right\}=0$ is exactly chiral and can be defined from the original overlap operator $D_{{\rm ov}}$ as $D_{c}=\frac{\rho D_{{\rm ov}}}{1-D_{{\rm ov}}/2}$ (Liu, 2005). More details regarding the calculation of the overlap operator and eigenmodes deflation in the inversion of the fermion matrix can be found in (Li et al., 2010). To calculate the matrix elements, we need first to construct 3-point correlation functions

[TABLE]

where $\chi$ is the nucleon interpolation field and $\mathcal{G}$ denotes the source grid. Then, we make a ratio of the 3-point correlation function to the nucleon 2-point function and extract the matrix element by fitting the ratio using the so-called two-state form

[TABLE]

Here $\Gamma_{e}$ is the non-polarized projector, $C_{2}(t_{f})=\sum_{\vec{x}}\langle\chi(t_{f},\vec{x})\bar{\chi}(0,\mathcal{G})\rangle$, $c$’s are fitting coefficients, and $\delta m$ is the effective energy difference between the ground state and the excited states. To better use this formula, multiple source-sink separations $t_{f}$ ranging from $\sim 0.7$ fm to $\sim 1.5$ fm are constructed for ${\Pi}^{q,g}(t_{f},\tau)$ on each ensemble for all the current positions $\tau$ between the source and sink.

As mentioned above, the 3-point correlation functions have two kinds of current insertions, CI and DI, as illustrated in Fig. 2. Since both the CI and glue matrix elements mix to DI through the renormalization of bare quantities under lattice regularization Yang et al. (2018a), the calculation of the ratio ${\cal R}$ under $\overline{\rm MS}$ scheme involves also the CI and glue contributions.

For the CI calculations, we use the stochastic sandwich method (SSM) (Yang et al., 2016a) with low-mode substitution (LMS) Li et al. (2010) to better control the statistical uncertainty. $Z_{3}$-noise grid sources with Gaussian smearing (for the 48I, 24I and 32I lattices) or block smearing (Liang et al., 2016) (for the 32ID lattice) are placed coherently at $t_{{\rm{src}}}=0$ and $t_{{\rm{src}}}=32$ ($t_{{\rm{src}}}=64$ also for 48I) in one inversion. Nucleon sinks are located at different positions with different separations in time from the source. Technical details regarding the LMS of a random $Z_{3}$ grid source and the use of SSM with LMS for constructing 3-point functions can be found in Refs. (Gong et al., 2013; Yang et al., 2016a; Liang et al., 2016). Setups regarding the valence sector of the CI case are listed in Table 3. Due to the fact that the multi-mass inversion algorithm is applicable to the overlap fermion with eigenvector deflation, we calculate four to seven valence masses for each of the four lattices.

For the DI calculations, we use the low-mode average (LMA) technique to calculate the quark loops which improves the signal-to-noise ratios. The low-mode part of the quark loops is calculated exactly since we have solved the low-lying eigenvectors of the overlap Dirac operator on all these lattices, The high-mode part is estimated with 8 sets of $Z_{4}$-noise on a 4-4-4-2 space-time grid with even-odd dilution and additional time shift (32 inversions in total). The same smeared $Z_{3}$-noise grid sources as used in the CI case are used in the production of the nucleon propagators. We make multiple measurements by shifting the source along the time direction to improve statistics. The spatial position of the center of the grid is randomly chosen for each source at different times to reduce autocorrelation. References (Gong et al., 2013, 2017; Yang et al., 2016b) contain more details regarding the DI calculation. When constructing quark loops, we include more valence quark masses to cover the strange region. The bare valence strange quark masses are determined on each lattice by the global-fit value at 2 GeV in the $\overline{{\rm MS}}$ scheme calculated in our previous study (Yang et al., 2015) and the nonperturbative mass renormalization constants calculated in (Liu et al., 2014); corresponding numbers are collected in Table 4. For all the 24I, 48I and 32I lattices the renormalized strange quark mass is around 100.5 MeV and for the 32ID lattice the number is around 95 MeV, which are all consistent with our global-fit value 101(3)(6) MeV (Yang et al., 2015) within error. We used the clover definition of the glue operator (Yang et al., 2018b) for the DI calculation of the glue momentum fraction. The cluster-decomposition error reduction (CDER) technique is applied to improve the signal (Liu et al., 2018; Yang et al., 2018b).

IV Renormalization

As demonstrated in (Liang et al., 2018), the renormalization can be processed separately for CI and DI and we will focus on the DI part in this work. The general form of the renormalized momentum fractions in DI $\langle x\rangle^{R,{\rm DI}}$ in the $\overline{\textrm{MS}}$ scheme at scale $\mu$ reads

[TABLE]

where $\langle x\rangle_{u,d,s}^{{\rm DI/CI}}$ is the bare quark momentum fraction in the DI/CI sector under lattice regularization, $\langle x\rangle_{g}$ is the glue momentum fraction, $Z_{QQ}^{\overline{\textrm{MS}}}(\mu)$ is the renormalization constant and $\delta Z_{QQ}^{\overline{\textrm{MS}}}(\mu)$ and $Z_{QG}^{\overline{\textrm{MS}}}(\mu)$ account for the mixing. To renormalize a lattice-regularized quantity with ${\overline{\textrm{MS}}}$ scheme, we first use the RI/MOM scheme to renormalize it at a scale $\mu_{R}$ nonperturbatively. And then, we convert the RI/MOM renormalized quantity to the ${\overline{\textrm{MS}}}$ scheme using a perturbatively calculated matching coefficient and evolve it to certain scale $\mu$. The complete renormalization is a combination of these two steps and can be expressed formally as $Z^{\overline{{\rm MS}}}(\mu)=\left[\left(Z(\mu_{R})R(\mu/\mu_{R})\right)|_{a^{2}\mu_{R}^{2}\rightarrow 0}\right]^{-1}$, where $Z(\mu_{R})$ and $R(\mu/\mu_{R})$ denote the RI/MOM renormalization and matching respectively.

In the nonperturbative renormalization procedure, one needs to carry out the RI/MOM renormalization several times with several quark masses on each lattice and extrapolate the results to the massless limit before the matching from RI/MOM to $\overline{\rm MS}$ since the massless renormalization scheme is used. Example plots for the isovector RI/MOM renormalization constants of the traceless diagonal piece of the EMT $Z_{QQ}$ as a function of the bare valence quark mass are shown in Fig. 3. Different colors denote different $a^{2}p^{2}$ scales. Fig. 3a is for the 24I lattice and Fig. 3b is for the 32ID lattice. It can be observed from the figures that the quark mass dependence is mild for both lattices at all scales, and linear fits (the solid lines in the figure) can be used to extrapolate the results to the massless limit.

In principle, this kind of extrapolation needs to be done also for the mixing coefficients. However, in practice, it is not necessary since from Fig. 3a we learned that the finite quark mass effect is quite small compared with the large statistical errors of the mixing coefficients. Also, we did not extrapolate the sea quark masses to zero. But our previous work Liu et al. (2014), where we used three 24I and three 32I lattices with different sea quark masses and several valence quark masses on each lattice to take both the valence and sea quark mass to zero for the massless renormalization scheme, shows that the sea quark mass effects are usually very weak for other renormalization constants (e.g., $Z_{A}$, $Z_{P}$ and $Z_{S}$). Accordingly, we can estimate the systematic uncertainty of the renormalization constants of this work due to the non-zero sea masses. Actually, our previous work Yang et al. (2018a) has estimated all the systematic uncertainties of the renormalization constants of EMT on these lattices and a full error budget can be found in its supplementary materials. All those systematic uncertainties of the renormalization constants are included in the present calculations.

We use the 3-loop result for the isovector matching coefficient (Gracey, 2003) while only 1-loop results exist for the others (Yang et al., 2016c). More detailed discussion about the calculation of nonperturbative renormalization and mixing is beyond the scope of this paper, and can be found in our previous works (Yang et al., 2018b, a). The renormalization constants used in this work are listed in Table 5. Although the mixing coefficients are of order $0.01$ or less, the mixing effects of this study are significant ($\sim 10\%$) since the DI bare values are themselves smaller than those of CI and glue.

V Results

The two bare matrix elements of the strange and light quarks are fitted using the two-state formula (Eq. (4)) in a joint correlated fit, such that the correlation between the two matrix elements is properly kept. This ensures the cancellation of the fluctuations of the two matrix elements in the ratio and leads to statistically more stable results.

Two-state fits are employed to handle the excited-state contaminations at finite source-sink separations. Fig. 4 shows example plots of the momentum fractions in DI on the 24I lattice at its unitary point. Different colors denote different source-sink separations. The colored lines in Panels (a) and (b) show the two-state fittings for each separation and their ratios are plotted in Panel (c). The blue bands indicate the final results at infinite source-sink separations. We use a joint fit involving both the light and strange quark, such that the ratio of the strange quark momentum fraction to that of light quark has much smaller relative error than the momentum fraction themselves due to the cancellation of the statistical fluctuations. All the fittings result in good $\chi^{2}/d.o.f.$ (around or less than one) and the final errors represented by the height of the bands are similar to the errors of the data points at the largest source-sink separation, which is reasonable and reassuring. Detailed fitting setups can be found in Table 6.

Another way to look at the fitting of the ratios is plotted in Fig. 5. There, the data points are modified such that the finite source-sink separation effects are removed using the results of the two-state fits. We can see that the new data points at different source-sink separations are all consistent with each other within errors and they all coincide with the fitted results (the blue bands), which means the two-state fits can successfully track the excited-state effects.

Although the main topic of this work is to calculate the ratio of DI, the CI part contributes too through the mixing. Example plots of the momentum fractions in CI for the 48I and 24I lattice are shown in Fig. 6. As in the DI case, all the fittings result in good $\chi^{2}/d.o.f.$ (around or less than one). Detailed fitting setups can be found in Table 7.

After renormalization, the final ${\cal R}$ ratios on different ensembles with different valence pion masses are fitted by the following form to track the pion mass, lattice spacing and volume dependence

[TABLE]

where the $C$’s are free parameters, $m_{\pi}^{v}$/$m_{\pi}^{s}$ is the valence/sea pion mass, and $m_{\pi}^{0}$ is the physical pion mass. The third term is to account for the partial quenching effect. A total of 21 data points are used for this global analysis. The extrapolated result to the physical limit is ${\cal R}^{\overline{\rm MS}}(2\,{\rm GeV})=\langle x\rangle_{s+\bar{s}}^{R}/\langle x\rangle_{u+\bar{u}}^{R}({\rm DI})=0.795(79)(77)$ with $\chi^{2}/d.o.f.=0.16$, where the first error is the statistical error and the second error is the total systematic one. A complete breakdown of the systematic uncertainties can be found in detail in Table 8. Details are discussed as follows.

The systematic uncertainty coming from the two-state fits is estimated in a way shown in Fig. 5. If the two-state fits work well, after the finite source-sink separation effect is removed from the data points by using the results of the two-state fits, all the points at different source-sink separations should lie on the same horizontal line with respect to the current insertion time. Fig. 5 shows that this is true within errors. We then use the difference between the modified data points and the two-state fitting results to estimate the corresponding systematic uncertainty to be $\sim 6\%$ of the central value of the final result.

The systematic uncertainties related to the global extrapolation, including finite lattice spacing effect, finite lattice volume effect, pion mass extrapolation with mixed action and strange quark mass effect, are estimated following Ref. Yang et al. (2016b). We use the difference between the continuum-extrapolated result and the result on our finest lattice, i.e. 32I, at the physical pion mass point and at the infinite volume limit for the systematic uncertainty of the finite lattice spacing effect. We have two coefficients for the finite lattice spacing effects in Eq. (6) since we are using two kinds of gauge actions. The results, $C_{3}^{I}=-0.4(5.3)$ and $C_{3}^{ID}=0.8(3.3)$, have no statistical significance for both cases, which is consistent with the behavior our data (the upper panel of Fig. 7). Numerically, the corresponding uncertainty is $|C_{3}^{I}|\times a^{2}({\rm 32I})\sim 0.003$. It is also consistent with some other calculations with the overlap fermion where there are no discernible discretization effects within statistics Liang et al. (2018). The lower panel of Fig. 7 shows the $\cal R$ dependence on lattice volume. The band shows the fitted volume dependence at the physical $m_{\pi}$ and at the continuum limit. From the band, visible finite volume effect can be seen at $m_{\pi}L<4$. Numerically, the coefficient of the volume effect $C_{4}$ in Eq. (6) is $-1.3(1.3)$. It has only a one-sigma signal but is more statistically significant than the coefficients of the lattice spacing effect. Using the same strategy, the corresponding systematic uncertainty is estimated to be $|C_{4}|\times e^{-m_{\pi,{\rm phy}}L(\rm{48I})}\sim 0.027$.

Some details about the chiral extrapolation are discussed as follows. For a partially quenched chiral extrapolation, one needs to separate the valence pion mass $m_{\pi,vv}$ and the sea pion mass $m_{\pi,ss}$ dependence. Additionally, for a mixed action case where the valence is a chiral fermion (as in our case where we use the overlap fermions), the LO chiral Lagrangian involves only one more low-energy constant Chen et al. (2007) and the correction vanishes as ${\cal O}(a^{2})$:

[TABLE]

where the term $a^{2}\Delta_{\rm mix}$ gauges the difference between the valence and sea lattice actions at a finite lattice spacing $a$. Also, we found that, when using overlap valence on domain wall sea, $\Delta_{\rm mix}$ is much smaller than those of other mixed action combinations. For example, it is an order of magnitude smaller than that of domain wall on staggered fermion Orginos and Walker-Loud (2008); Aubin et al. (2008) and only shifts $m_{\pi,vs}$ from the would be unitary mass $m_{\pi}$ by $\sim 10$ MeV for $m_{\pi}\sim 300$ MeV at $a=0.11$ fm Lujan et al. (2012). In the chiral extrapolation of the ratio $\cal R$, the errors of the data are around 10% or even larger and one cannot fit complicated pion mass dependence except for the leading linear terms in $m_{\pi,vv}^{2}$ and $m_{\pi,ss}^{2}$, in which case the delta mix effects are already included in the $a^{2}$ extrapolation. Additional and higher order terms (including the ${m^{\prime}}_{\pi,vs}$ terms appearing in their third power) show no statistical significance which means the difference caused by ignoring higher order terms is much smaller than the current statistical uncertainty. For the strange quark mass part, as mentioned before, the bare valence strange quark masses are determined by the global-fit value at 2 GeV in the $\overline{{\rm MS}}$ scheme calculated in our previous study (Yang et al., 2015). For all the 24I, 48I and 32I lattices the renormalized strange quark mass is around 100.5 MeV and for the 32ID lattice the number is around 95 MeV, which are all consistent with our global-fit value 101(3)(6) MeV (Yang et al., 2015) within error. As discussed in Ref. Yang et al. (2016b), the valence strange quark mass effect is very weak and we assume so for the sea strange quark too. So combining both the valence and sea effects, we estimate the systematic uncertainty due to strange quark mass to be around 1%.

The systematic uncertainties corresponding to the renormalization are discussed in detail in the supplemental materials of Ref. Yang et al. (2018a). They contribute in total $\sim 2\%$ systematic uncertainty of ${\cal R}$. In addition, we did not consider the charm quark contribution by mixing. However, the charm momentum fraction itself is small. More importantly, the sea quark-to-quark mixing coefficients are of order one percent as given in our paper, so the neglect of the heavy flavor contributions is safe and we estimate that it can only leads to $\sim 0.5\%$ uncertainty.

We use the traceless diagonal part of the EMT $\bar{T}_{44}$ to carry out our calculation. Actually, in our previous papers discussing the glue momentum fraction renormalization Yang et al. (2018b) and nucleon mass decomposition Yang et al. (2018a), we found that the assumed rotational symmetry breaking effect of the renormalization constant between the off-diagonal EMT $T_{4i}$ and the traceless diagonal part $\bar{T}_{44}$ which are in different irreducible representations of the cubic point group $O_{h}$ is much smaller than the statistical uncertainties for both the glue and quark case. So we do not add one more systematic uncertainty due to the rotational symmetry breaking. In total, as shown in the table, the systematic uncertainty of ${\cal R}$ is about 9.7% which is close to the statistical uncertainty.

The bands in Fig. 8 show our final prediction on the unitary pion mass dependence of the ratio in the continuum and infinite volume limits. The width of the blue band indicates the statistical error and the width of the wider cyan band the total error. The data points from different ensembles with the partially quenching effect subtracted are also plotted in the figure. The agreement of the bands and the data points shows that the finite volume and lattice spacing effects are small.

VI Discussion and Summary

To manifest that the CS and DS partons have distinct small $x$ behaviors, the ratio of the distribution $(s(x)+\bar{s}(x))/(\bar{u}(x)+\bar{d}(x))$ from 3 global fittings at NNLO Ball et al. (2017); Dulat et al. (2016); Harland-Lang et al. (2015) at $Q^{2}=4\,{\rm GeV^{2}}$ is plotted in Fig. 9. We see that the ratios are kind of flat at small $x\lesssim 10^{-2}$ with large errors. On the other hand, it is conspicuous that they all have a characteristic shoulder with a fall off around $x\sim 10^{-2}$ toward larger $x$. Since under our classification, $s(x)+\bar{s}(x)$ contains only the DS contribution while $\bar{u}(x)+\bar{d}(x)$ includes both CS and DS, this reflects the fact that the small $x$ behavior of $\bar{q}^{ds}(x)$ is more singular than that of $q^{v}$ and $\bar{q}^{cs}$ such that at $x\lesssim 10^{-2}$ the DS dominates, so the ratio stays roughly constant Liu (2000); Liu et al. (2012). When $x$ is larger than $\sim 10^{-2}$, the CS $\bar{u}^{cs}(x)+\bar{d}^{cs}(x)$ component in $\bar{u}(x)+\bar{d}(x)$ (N.B. $\bar{u}(x)+\bar{d}(x)=\bar{u}^{cs}(x)+\bar{d}^{cs}(x)+\bar{u}^{ds}(x)+\bar{d}^{ds}(x)$) sets in to make the ratio smaller. This can also be understood in Regge theory where $q^{v},\bar{q}^{cs}_{\stackrel{{\scriptstyle\longrightarrow}}{{x\rightarrow 0}}}x^{-1/2}$ since the CS partons are in the connected insertion which is flavor non-singlet as are the valence partons and their small $x$ behaviors reflect the leading Reggeon exchanges Brodsky and Schmidt (1991). On the other hand, $q^{ds},\bar{q}^{ds}_{\stackrel{{\scriptstyle\longrightarrow}}{{x\rightarrow 0}}}x^{-1}$ since the DS is flavor singlet and can have Pomeron exchanges Kuti and Weisskopf (1971); Reya (1981). Also, we find numerically that in the small $x$ region ($10^{-4}$ to $10^{-2}$) of the global fittings for PDF Ball et al. (2017); Dulat et al. (2016); Harland-Lang et al. (2015); Buckley et al. (2015) at $Q^{2}=4$ GeV2, the power of the small $x$ behavior for $\bar{q}(x)$ is in the range $[-1.22,-1.15]$ and for $q^{v}(x)$ in the range [-0.6, -0.2], which are close to those prescribed in Regge theory and consistent with our argument that DS dominates the small $x$ behavior.

The dominance of the DS partons at small $x$ reveals the possibility that lattice QCD can help to constrain the strange quark distribution in the small $x$ region using the DI ratio $\cal R$ as defined in Eq. (1). To do so, it is advocated Liu (2017) to completely separate the CS and DS partons with their corresponding evolutions in new global fittings. No approximation such as $u^{ds}(x)=\bar{u}^{ds}(x)$ is needed in this approach. Once done, close connection can be built between lattice computations of moments and those from global fittings.

In addition, before new global fittings that separate CS and DS are carried out, an easy-to-implement first trial is to use the following approximate form for the strange PDF

[TABLE]

where $c$ is a coefficient to be fixed, to better control the statistical uncertainty of the global fittings. The form is based on the approximation that $s(x)+\bar{s}(x)$ is proportional to $\bar{u}^{ds}(x)+\bar{d}^{ds}(x)$ by the factor of $\cal R$ if $u^{ds}(x)=\bar{u}^{ds}(x)$ and isospin symmetry are assumed. Since the first term $\bar{u}(x)+\bar{d}(x)$ in Eq. (8) contains both $\bar{u}^{ds}(x)+\bar{d}^{ds}(x)$ and $\bar{u}^{cs}(x)+\bar{d}^{cs}(x)$, we use the second term to subtract the latter which is taken to be proportional to $\bar{u}(x)-\bar{d}(x)$ due to the fact that it equals $\bar{u}^{cs}(x)-\bar{d}^{cs}(x)$ in the isospin limit. Both $\bar{u}(x)+\bar{d}(x)$ and $\bar{u}(x)-\bar{d}(x)$ in Eq. (8) are obtained in global fittings, thus this form should be easy to implement. This form serves as an explicit example of how a lattice result enters directly the global fittings of PDFs. Further lattice calculation of the fourth moment $\langle x^{3}\rangle$ of the DI will serve to gauge the validity of this approach and suggest possible modification of the fitting function. We also plot our $\cal R$ in Fig. 9 up to $x=10^{-2}$ to visually show that the uncertainty of the strange PDF can be reduced significantly by using this lattice constraint as indicated by the lattice error as compared to those from the NNLO analyses.

Since lattice calculations with low quark and glue moments are getting mature and complete, the QCD path-integral classification that separates CS and DS extends the ability that lattice calculations can serve as meaningful constraints for the global analysis of PDFs. The present result of the ratio $\langle x\rangle_{s+\bar{s}}/\langle x\rangle_{u+\bar{u}}({\rm DI})=0.795(79)(77)$ at $\overline{\rm MS}$ scale $\mu=2$ GeV is the first such calculation that can constrain the global fittings in the small $x$ region. Our results will have important impact on future global fittings and experiments from EIC and LHC.

Acknowledgements.

We thank the RBC and UKQCD Collaborations for providing their DWF gauge configurations. This work is supported in part by the U.S. DOE Grant No. DE-SC0013065 and DOE Grant No. DE-AC05-06OR23177 which is within the framework of the TMD Topical Collaboration. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This work used Stampede time under the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1053575. We also thank the National Energy Research Scientific Computing Center (NERSC) for providing HPC resources that have contributed to the research results reported within this paper. We acknowledge the facilities of the USQCD Collaboration used for this research in part, which are funded by the Office of Science of the U.S. Department of Energy.