Nucleon structure functions from lattice operator product expansion

TL;DR

This paper introduces a novel lattice QCD method using the operator product expansion to directly compute nucleon structure functions from the Compton amplitude, addressing previous renormalization and operator mixing challenges.

Contribution

It proposes a new approach to calculate nucleon structure functions directly from lattice QCD, overcoming prior renormalization and operator mixing issues.

Findings

Method successfully computes structure functions from the Compton amplitude.

Addresses renormalization and operator mixing problems in lattice calculations.

Lays groundwork for first-principles nucleon structure studies.

Abstract

Deep-inelastic scattering, in the laboratory and on the lattice, is most instructive for understanding how the nucleon is built from quarks and gluons. The long-term goal is to compute the associated structure functions from first principles. So far this has been limited to model calculations. In this Letter we propose a new method to compute the structure functions directly from the virtual, all-encompassing Compton amplitude, utilizing the operator product expansion. This overcomes issues of renormalization and operator mixing, which so far have hindered lattice calculations of power corrections and higher moments.

Equations

$$\mu (q^2)=c_2(q^2{a}^2)v_2(a)+\frac{c_4(q^2{a}^2)}{q^2}{v}_4(a)+\cdots ,$$

$$T_{\mu \nu }(p,q)={\rho }_{\lambda {\lambda }^{\prime }}\int {\mathrm{d}}^{4}x{\mathrm{e}}^{iq\cdot x}⟨p,{\lambda }^{\prime }|T{J}_{\mu }(x)J_{\nu }(0)|p,\lambda ⟩,$$

$$T_{\mu \nu }(p,q)=\left({\delta }_{\mu \nu }-\frac{q_{\mu }{q}_{\nu }}{q^2}\right){\mathcal{F}}_1(\omega ,q^2)+\left(p_{\mu }-\frac{p\cdot q}{q^2}{q}_{\mu }\right)\left(p_{\nu }-\frac{p\cdot q}{q^2}{q}_{\nu }\right)\frac{1}{p\cdot q}{\mathcal{F}}_2(\omega ,q^2),$$

$$\mathrm{Im}{\mathcal{F}}_1(\omega ,q^2)=2\pi F_1(\omega ,q^2),\mathrm{Im}{\mathcal{F}}_2(\omega ,q^2)=2\pi F_2(\omega ,q^2),$$

$$\begin{array}{rl}{T}_{\mu \nu }(p,q)=\sum _{n=2,4,\cdots }^{\infty }& \{\left({\delta }_{\mu \nu }-\frac{q_{\mu }{q}_{\nu }}{q^2}\right)4{\omega }^n{\int }_0^{1}dx{x}^{n-1}{F}_1(x,q^2)\\ & +\left(p_{\mu }-\frac{p\cdot q}{q^2}{q}_{\mu }\right)\left(p_{\nu }-\frac{p\cdot q}{q^2}{q}_{\nu }\right)\frac{8}{2p\cdot q}{\omega }^{n-1}{\int }_0^{1}dx{x}^{n-2}{F}_2(x,q^2)\}.\end{array}$$

$$\begin{array}{rl}{T}_{\mu \nu }(p,q)& =\left({\delta }_{\mu \nu }-\frac{q_{\mu }{q}_{\nu }}{q^2}\right)4\omega {\int }_0^{1}dx\frac{\omega x}{1-(\omega x{)}^2}{F}_1(x,q^2)\\ & +\left(p_{\mu }-\frac{p\cdot q}{q^2}{q}_{\mu }\right)\left(p_{\nu }-\frac{p\cdot q}{q^2}{q}_{\nu }\right)\frac{8\omega }{2p\cdot q}{\int }_0^{1}dx\frac{1}{1-(\omega x{)}^2}{F}_2(x,q^2).\end{array}$$

$$\begin{array}{rl}{\mathcal{F}}_1(\omega ,q^2)& =2\omega {\int }_1^{\infty }d\bar{\omega }\left[\frac{F_1(\bar{\omega },q^2)}{\bar{\omega }(\bar{\omega }-\omega )}-\frac{F_1(\bar{\omega },q^2)}{\bar{\omega }(\bar{\omega }+\omega )}\right]+{\mathcal{F}}_1(0,q^2),\\ {\mathcal{F}}_2(\omega ,q^2)& =2\omega {\int }_1^{\infty }d\bar{\omega }\left[\frac{F_2(\bar{\omega },q^2)}{\bar{\omega }(\bar{\omega }-\omega )}+\frac{F_2(\bar{\omega },q^2)}{\bar{\omega }(\bar{\omega }+\omega )}\right].\end{array}$$

$${\mathcal{F}}_1(\omega ,q^2)=4{\omega }^2{\int }_0^{1}dxx\frac{F_1(x,q^2)}{1-(\omega x{)}^2},{\mathcal{F}}_2(\omega ,q^2)=4\omega {\int }_0^{1}dx\frac{F_2(x,q^2)}{1-(\omega x{)}^2},$$

$$T_{33}(p,q)=\sum _{n=2,4,\cdots }^{\infty }4{\omega }^n{\int }_0^{1}dx{x}^{n-1}{F}_1(x,q^2)$$

$$T_{33}(p,q)=4\omega {\int }_0^{1}dx\frac{\omega x}{1-(\omega x{)}^2}{F}_1(x,q^2).$$

$$\mathcal{L}(x)\to \mathcal{L}(x)+\lambda {\mathcal{J}}_3(x),{\mathcal{J}}_3(x)=Z_V\cos (\vec{q}\cdot \vec{x})e_f{\bar{\psi }}_f(x){\gamma }_3{\psi }_f(x),$$

$$-2E_{\lambda}(p,q)\,\frac{\partial^{2}}{\partial\lambda^{2}}E_{\lambda}(p,q)\,\big{|}_{\lambda=0}=T_{33}(p,q)\,.$$

$$t_i=T_{33}({\omega }_i),i=1,\cdots ,N$$

$$\left(\begin{tabular}[]{c}$t_{1}$\\ $t_{2}$\\ $\vdots$\\ $t_{N}$\end{tabular}\right)\,=\,\left(\begin{tabular}[]{cccc}$4\omega_{1}^{2}$&$4\omega_{1}^{4}$&$\cdots$&$4\omega_{1}^{2M}$\\ $4\omega_{2}^{2}$&$4\omega_{2}^{4}$&$\cdots$&$4\omega_{2}^{2M}$\\ $\vdots$&$\vdots$&$\vdots$&$\vdots$\\ \;$4\omega_{N}^{2}$&\;$4\omega_{N}^{4}$&$\cdots$&$4\omega_{N}^{2M}$\end{tabular}\right)\;\left(\begin{tabular}[]{c}$\mu_{1}$\\ $\mu_{3}$\\ $\vdots$\\ \hskip 14.22636pt$\mu_{2M-1}$\end{tabular}\right)\,.$$

$$T_{33}(\omega)=4\,\big{(}\omega^{2}\mu_{1}+\omega^{4}\mu_{3}+\cdots+\omega^{2M}\mu_{2M-1}\big{)}\,.$$

$$f_i=F_1(x_i),K_{ij}=\frac{4{\omega }_i^2{x}_j}{1-({\omega }_i{x}_j{)}^2}.$$

$$t_i=ϵ\sum _{j=1}^M{K}_{ij}{f}_j,i=1,\cdots ,N.$$

$$K=U\left[\mathrm{diag}(w_1,\cdots ,w_N)\right]V^T.$$

$$f_j=\sum _{i=1}^N{K}_{ji}^{-1}{ϵ}^{-1}{t}_i,$$

$$K^{-1}=V\left[\mathrm{diag}(1/w_1,\cdots ,1/w_L,0,\cdots ,0)\right]U^T$$

$$6x{F}_1^{u-d}(x)=x\left[u(x)-d(x)\right]$$

$$\mu (s)=A(s+\alpha {)}^{-\beta },$$

Full text

– QCDSF Collaboration –

The connection between the deep-inelastic structure functions and the quark structure of the nucleon is commonly rendered by the parton model. Although providing an intuitive language, in which to interpret the deep-inelastic scattering data, the parton model is incomplete. The theoretical basis is the operator product expansion (OPE). The operators in the expansion are classified according to twist. The parton model accounts for twist-two contributions only, and cannot accommodate power corrections arising from operators of higher twist. Power corrections are inseparably connected with the leading-twist contributions, as a result of operator mixing [1]. Consider, for example, a generic moment of any deep-inelastic structure function of the nucleon regularized on the lattice by a hard cut-off $1/a$,

[TABLE]

where $v_{2}=\langle N|\mathcal{O}_{2}(a)|N\rangle$ and $v_{4}=\langle N|\mathcal{O}_{4}(a)|N\rangle$ are reduced nucleon matrix elements of local operators of twist two and four, respectively, and $c_{2}$ and $c_{4}$ are the corresponding reduced Wilson coefficients. The operator $\mathcal{O}_{4}$ mixes with $\mathcal{O}_{2}$ with mixing coefficients which diverge as $1/a^{2}$. The power divergences of $v_{4}(a)$ must be cancelled with those of $c_{2}(q^{2}a^{2})$, which demands a nonperturbative calculation of the Wilson coefficient as well. This can be accomplished by an entirely nonperturbative calculation of the structure functions only.

So far lattice calculations of nucleon structure functions have been limited to calculations on the parton level [2]. In this Letter we propose a method that goes beyond that and computes the deep-inelastic structure functions, including power corrections, directly from the product of electromagnetic currents. This approach has been called ‘OPE without OPE’ elsewhere [3]. For previous and related work on the subject see [4, 5, 6].

The starting point is the forward Compton amplitude of the nucleon [7],

[TABLE]

the time ordered product of electromagnetic currents sandwiched between nucleon states of momentum $p$ and polarization $\lambda$, where $q$ is the momentum of the virtual photon and $\rho$ is the polarization density matrix. For simplicity, we will restrict ourselves to unpolarized structure functions only with $2\rho=\mathbf{1}$. We are then left with

[TABLE]

where $\omega=2p\cdot q/q^{2}$. Euclidean metric is understood. Crossing symmetry, $T_{\mu\nu}(p,q)=T_{\nu\mu}(p,-q)$, implies that $\mathcal{F}_{1}$ is an even function of $\omega$ and $\mathcal{F}_{2}$ an odd function, $\displaystyle\mathcal{F}_{1}(-\omega,q^{2})=\mathcal{F}_{1}(\omega,q^{2})$, $\displaystyle\mathcal{F}_{2}(-\omega,q^{2})=-\mathcal{F}_{2}(\omega,q^{2})$. In the physical region $1\leq|\omega|\leq\infty$

[TABLE]

where $F_{1}$ and $F_{2}$ are the deep-inelastic structure functions of the nucleon. Using the OPE, one can express $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ in terms of moments of $F_{1}$ and $F_{2}$, which are amenable to calculation on the Euclidean lattice. Alternatively, $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ can be written as dispersion integrals over $\omega$, which leads to the same expressions.

Let us first consider the OPE of $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$. After some simple algebra we obtain [7]

[TABLE]

The series $\sum_{k\in\mathbf{N}}\,(\omega x)^{2k}$ in (5) is geometric and sums up to $\displaystyle[1-(\omega x)^{2}]^{-1}$, which leads to the alternate expression

[TABLE]

In the limit where $F_{1}(x,q^{2})$ and $F_{2}(x,q^{2})$ become independent of $q^{2}$ we have the Callan-Gross relation $\displaystyle F_{2}(x)=2xF_{1}(x)$.

Alternatively, we can express $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ directly in terms of the structure functions $F_{1}$ and $F_{2}$, circumventing the OPE. The amplitudes $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ have cuts at $-\infty\leq\omega\leq-1$ and $1\leq\omega\leq\infty$ with discontinuities (4). This leads to once subtracted dispersion relations

[TABLE]

While $\mathcal{F}_{2}(0,q^{2})=0$, the subtraction constant $\mathcal{F}_{1}(0,q^{2})$ contains information on the magnetic polarizability of the nucleon and the proton–neutron electromagnetic mass shift [8]. In the following equations we shall discard it, as it has no counterpart in $F_{1}$, nor is it accounted for by the OPE. It can be computed like any other value of $\mathcal{F}_{1}$ though and, if necessary, has to be subtracted from $\mathcal{F}_{1}(\omega,q^{2})$. (So, for example, from the data underlying Fig. 6.) Substituting $\bar{\omega}$ by $1/x$, we finally obtain

[TABLE]

where we have identified $F_{1}(\bar{\omega},q^{2})$ and $F_{2}(\bar{\omega},q^{2})$ with $F_{1}(x,q^{2})$ and $F_{2}(x,q^{2})$, respectively. If we insert (8) into (3), we obtain (6), in agreement with the OPE resummed. It should be noted that the structure functions $F_{1}(x,q^{2})$ and $F_{2}(x,q^{2})$ include higher twist contributions, as we have not made any assumptions on $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ other than on the analytic structure.

To simplify the numerical calculation, we may choose $\mu=\nu=3$ and $p_{3}=q_{3}=q_{4}=0$. We then have

[TABLE]

and, alternatively,

[TABLE]

For $|\omega|>1$ the principal value has to be taken. The matrix element $T_{33}(p,q)$ can be computed most efficiently, including singlet matrix elements, by a simple extension of existing implementations of the Feynman-Hellmann technique to lattice QCD [9]. For simplicity, we consider the local vector current only. The appropriate renormalization factor $Z_{V}$ can be computed unambiguously [10]. No further renormalization is needed. To compute the Compton amplitude from the Feynman-Hellmann relation, we introduce the perturbation to the Lagrangian

[TABLE]

where $\psi_{f}$ is the quark field of flavor $f=u,d,s,\cdots$ to which the photon is attached, and $e_{f}$ is its electric charge. Note that $\lambda$ has dimension mass. Taking the second derivative of the nucleon two-point function $\displaystyle\langle N(\vec{p},t)\bar{N}(\vec{p},0)\rangle_{\lambda}\simeq C_{\lambda}\,{\rm e}^{-E_{\lambda}(p,q)\,t}$ with respect to $\lambda$ on both sides, we obtain

[TABLE]

The derivation of (12) would go beyond the scope of this Letter and will be presented in a separate publication. Provided we compute at sufficiently large $q^{2}$, standard factorization theorems state that the Compton amplitude will be dominated by the ‘handbag’ diagram shown in the left panel of Fig. 1. Nevertheless, the amplitude does encompass all contributions, including the power-suppressed ‘cats-ears’ diagram shown in the right panel. Thereby, varying $q^{2}$ will allow us to test the twist expansion and, in particular, isolate twist-four contributions. A conventional calculation of the four-point function $\langle p,\lambda|J_{\mu}(x)J_{\nu}(0)|p,\lambda\rangle$, in contrast, would involve all-to-all quark propagators twice consecutively, which is not practicable.

Knowing $T_{33}(p,q)$ for sufficiently many values of $\omega$, we can solve (9) for the moments $\mu_{n}=\int_{0}^{1}dxx^{n}F_{1}(x,q^{2})$ and (10) directly for $F_{1}(x,q^{2})$. We want to keep the photon momentum $\vec{q}$ fixed, but vary the nucleon momentum $\vec{p}$. This amounts to a relatively cheap calculation, as all nucleon energies $E_{\lambda}(p,q)$ can be computed from a single set of background field configurations (per value of $\lambda$). For nonsinglet quantities, in which the currents couple to valence quarks only, no additional gauge field configurations will have to be generated at all.

Regarding (9), the task is to compute (say) the lowest $M$ moments $\mu$ from a finite number of sampled points

[TABLE]

with $M\leq N$. For ease of writing we have dropped the dependence on $q^{2}$. Under these conditions the OPE of $T_{33}$ can be written as a set of equations

[TABLE]

The matrix is known as Vandermonde matrix. Efficient algorithms for solving (14) can be found in the literature [11]. Alternatively, we can cast the $T_{33}(\omega_{i})$’s into a simple functional form by the interpolating polynomial

[TABLE]

The moments $\mu$ can then be determined from a fit of (15) to the sampled points (13).

To solve (10), we approximate the integral by a sum over a discrete set of $M$ points, $0<x_{1}<x_{2}<\cdots<x_{M}<1$, and write

[TABLE]

We assume the points to be equidistant with step size $\epsilon$. A generalization to adaptive step sizes is straightforward. The integral equation (10) then reduces to the set of equations

[TABLE]

In general, $N<M$. The $N\times M$ matrix $K$ can be written as the product of a $N\times N$ orthogonal matrix $U$, a $N\times N$ diagonal matrix $W$ with positive or zero eigenvalues $w_{1}<w_{2}<\cdots<w_{N}$, and the transpose of a row-orthogonal $M\times N$ matrix $V$,

[TABLE]

The matrix $W$ is singular. Singular value decomposition (SVD) is the method of choice for solving (10) under such conditions. The solution is

[TABLE]

where $K^{-1}$ is the pseudoinverse

[TABLE]

with $1/w_{l}$ being replaced by zero if $w_{l}=0$, which is assumed for $L<l\leq N$. One has to excercise some discretion at deciding at what threshold to set $1/w_{l}$ to zero. Several routines, such as PseudoInverse of Mathematica [12], solve this problem automatically.

The mathematical methods described above enable us to obtain the lower moments of the structure function, as well as $F_{1}(x)$ as a whole, rather accurately from a relatively small set of values of $T_{33}(\omega)$. To demonstrate that, we start from the experimental nonsinglet structure function $F_{1}^{u-d}(x)$, which we parameterize by

[TABLE]

with $u(x)$ and $d(x)$ taken to be the LO parton distributions at the scale $q^{2}=1\,\mbox{GeV}^{2}$ [13]. From (21) we compute $T_{33}(\omega)$ for $\omega=0.1$, $0.2$, $\cdots$, $0.9$ using (10). The result is shown in Fig. 2. This we consider our ‘data’, from which we want to retrieve the moments $\mu$ and structure function $F_{1}^{u-d}(x)$, and to which the lattice results (see Fig. 6) will aspire to.

To compute the moments, we fit the interpolating polynomial (15) to the lowest five ‘data’ points $t_{1},\cdots,t_{5}$ with $M=5$, which is equivalent to solving the set of linear equations (14) for $\mu_{1},\mu_{3},\cdots,\mu_{9}$. The result is shown in Fig. 3 for the lowest three moments, together with the target numbers. The first three moments are accurately reproduced.

To retrieve the structure function $F_{1}^{u-d}(x)$, we apply the SVD (19) to our ‘data’ [12]. Three out of nine eigenvalues of $W$ turned out to be zero. In Fig. 4 we show the result for $M=19$. It shows that the structure function $F_{1}^{u-d}$ can be well reproduced from a relatively small set of data, except perhaps for $x\lesssim 0.05$. Similar results are obtained for the singlet structure function $F_{1}^{u+d+\bar{u}+\bar{d}+\bar{s}}$. We have not made any attempts to optimize the SVD. It can be improved in several respects. A Bayesian approach [14] to alleviate overfitting, for example, might lead to particularly robust results.

There are other possibilities as well to compute the structure function from the Compton amplitude. A particularly promising approach is to fit the moments, for example in the interpolating polynomial (15), by an appropriate function $\mu(s)$ with $\mu(n)=\mu_{n}$ and employ an inverse Mellin transform on $\mu(s)$ to obtain $F_{1}(x)$. It turns out that the moments can be fitted surprisingly well by the simple expression

[TABLE]

for which the inverse Mellin transform is known analytically [15]. Starting from the moments $6\,\mu_{n}=\int_{0}^{1}dxx^{n}\left[u(x)-d(x)\right]$, the result of the Mellin transform is shown in Fig. 5.

The analysis so far has been limited to $\omega\in[0,1]$. The SVD method can be extended to larger values $\omega>1$ without problem. This will allow us to probe the small-$x$ region of $F_{1}(x)$, which is not accessible through moments of the structure function. Indeed, by extending the calculation to $\omega=2$, we were able to retrieve the singlet structure function $F_{1}^{u+d+\bar{u}+\bar{d}+\bar{s}}(x)$ [13] down to fractional momenta $x\lesssim 0.005$, which was not possible before. Odd moments of the structure functions can be obtained by also including the local axial vector current $\bar{\psi}_{f}(x)\gamma_{3}\gamma_{5}\psi_{f}(x)$ to (11) and studying the interference with the vector current. This is achievable through a simple extension of the procedure described above. The method can be generalized to nonforward Compton scattering as well. That will allow us to derive generalized parton distribution functions (GPDs).

There is the question what accuracy can be achieved with real data. It turns out that the second derivative of the nucleon energy can be computed rather accurately. In a proof-of-principle study we have computed (12) from $O(900)$ configurations generated at the SU(3) symmetric point [16] on a $32^{3}\times 64$ lattice with lattice spacing $a\approx 0.074\,\mbox{fm}$. First results are presented in Fig. 6, where the contribution from $\displaystyle\mathcal{F}_{1}(0,q^{2})=\left.-2\,E_{\lambda}(0,q)\,\partial^{2}E_{\lambda}(0,q)/\partial\lambda^{2}\,\right|_{\lambda=0}$ has been subtracted. The precision for lattice momenta $\vec{p}^{2}=1$ and $2$ is already quite impressive. We should be able to improve on the precision of the data at higher momenta by employing ‘momentum smearing’ techniques [18], which has not been attempted here. Based on this result, we consider an overall projected error of $10\%$ on the Compton amplitude $T_{33}$, marked by the shaded area in Fig. 2, a conservative estimate. The resulting errors on $\mu_{n}$ and $F_{1}^{u-d}(x)$, shown in Figs. 3, 4 and 5, have been obtained from replacing the initial values of $T_{33}(\omega_{i})$ by the corresponding numbers on the confidence envelopes. This leads us to conclude that the entire structure function, including its moments, can be reconstructed from a lattice calculation of the Compton amplitude with unprecedented accuracy, devoid of any renormalization and mixing issues.

Acknowledgement

GS thanks Akaki Rusetsky for useful discussions. The numerical calculations were carried out on the IBM BlueGene/Qs at Edinburgh and Jülich using DIRAC 2 and NIC resources, on the Cray XC30 at HLRN, Berlin and Hannover, and on the NCI National Facility at Canberra. HP and GS are supported by DFG Grant Nos. SCHI 422/10-1 and SCHI 179/8-1. PELR is supported by the STFC under contract ST/G00062X/1. RDY and JMZ are supported by the Australian Research Council Grant Nos. FT120100821, FT100100005 and DP140103067.