QCD evolution of the orbital angular momentum of quarks and gluons: Genuine twist-three part

TL;DR

This paper numerically solves the one-loop QCD evolution equations for the genuine twist-three orbital angular momentum distributions of quarks and gluons, enhancing understanding of their scale dependence in polarized nucleons.

Contribution

It provides the first numerical solution for the evolution of genuine twist-three OAM distributions and explores the scale dependence of potential angular momentum.

Findings

Evolution of twist-three OAM distributions is identical to Efremov-Teryaev-Qiu-Sterman function.

The evolution of OAM distributions is now well-understood when combined with Wandzura-Wilczek parts.

First analysis of the scale dependence of potential angular momentum.

Abstract

We present the numerical solution of the one-loop QCD evolution equation for the genuine twist-three part of the orbital angular momentum (OAM) distributions of quarks and gluons inside a longitudinally polarized nucleon. This is based on the observation that the evolution is identical to that of the Efremov-Teryaev-Qiu-Sterman function for transverse single spin asymmetry. Together with the known evolution of the Wandzura-Wilczek part, the one-loop evolution of OAM distributions is now practically under control. We also study, for the first time, the scale dependence of the potential angular momentum defined as the difference between the Ji and Jaffe-Manohar definitions of OAM.

Tables

Table 1. Table 1 : Numerical results of the sum rules.

	Initial condition 1	Initial condition 2	Initial condition 3
${\int }_0^1𝑑x{L}_{3g,F}\left(x\right)$ at ${\mu }_0^2=1$ GeV2	$0.0006$	$0.0006$	$0.0002$
${\int }_0^1𝑑x{L}_{3g,F}\left(x\right)$ at ${\mu }^2=10$ GeV2	-0.00007	0.0004	0.00016
${\int }_0^1𝑑x\left(L_{3q}\left(x\right)+L_{3g,S}\left(x\right)\right)$ at ${\mu }_0^2=1$ GeV2	$4\times {10}^{-7}$	-0.0006	-0.0002
${\int }_0^1𝑑x\left(L_{3q}\left(x\right)+L_{3g,S}\left(x\right)\right)$ at ${\mu }^2=10$ GeV2	0.00014	-0.0007	-0.0002
$L_{\mathrm{pot}}$ at ${\mu }_0^2=1$ GeV2	0.234	0.185	0.069
$L_{\mathrm{pot}}$ at ${\mu }^2=10$ GeV2	0.162	0.077	0.032

Equations

$$\frac{1}{2}=\frac{1}{2}\Delta \Sigma +\Delta G+L_q+L_g,$$

$$L_{q,g}(x)=L_{q,g}^{WW}(x)+L_{3q,3g}(x).$$

$$L_q={\int }_0^{1}dx{L}_q(x),L_g={\int }_0^{1}dx{L}_g(x).$$

$$L_q(x)$$

$$L_g(x)$$

$$\int \frac{d\lambda d\zeta }{(2\pi {)}^2}{e}^{i\frac{\lambda }{2}(x_1+x_2)+i\zeta (x_2-x_1)}⟨P^{\prime }{S}^{\prime }|\bar{\psi }(-\lambda n/2){\gamma }^{+}g{F}^{+i}(\zeta n)\psi (\lambda n/2)|PS⟩=P^{+}{ϵ}^{+i\rho \sigma }{S}_{\rho }{\Delta }_{\sigma }{\Phi }_F(x_1,x_2)+\cdots ,$$

$$\int \frac{d\lambda d\zeta }{(2\pi {)}^2}{e}^{i\frac{\lambda }{2}(x_1+x_2)+i\zeta (x_2-x_1)}⟨P^{\prime }{S}^{\prime }|\bar{\psi }(-\lambda n/2){\gamma }^{+}(-i{\gamma }_5)g{F}^{+i}(\zeta n)\psi (\lambda n/2)|PS⟩=P^{+}{S}^{+}{\Delta }^i{\tilde{\Phi }}_F(x_1,x_2)+\cdots ,$$

$$\int \frac{d\lambda d\zeta }{(2\pi {)}^2}{e}^{i\frac{\lambda }{2}(x_1+x_2)+i\zeta (x_2-x_1)}⟨P^{\prime }{S}^{\prime }|F^{+\alpha }(-\lambda n/2)g{F}^{+i}(\zeta n)F_{\alpha }^{+}(\lambda n/2)|PS⟩=(P^{+}{)}^2{ϵ}^{+i\rho \sigma }{S}_{\rho }{\Delta }_{\sigma }{M}_F(x_1,x_2)+\cdots ,$$

$$\int \frac{d\lambda d\zeta }{(2\pi {)}^2}{e}^{i\frac{\lambda }{2}(x_1+x_2)+i\zeta (x_2-x_1)}⟨P^{\prime }{S}^{\prime }|{ϵ}_{ik}{F}^{+i}(-\lambda n/2)g{F}^{+j}(\zeta n)F^{+k}(\lambda n/2)|PS⟩=-(P^{+}{)}^2{S}^{+}{\Delta }^j{\tilde{M}}_F(x_1,x_2)+\cdots ,$$

$${\Phi }_F(x_1,x_2)=-{\Phi }_F(x_2,x_1),{\tilde{\Phi }}_F(x_1,x_2)={\tilde{\Phi }}_F(x_2,x_1),$$

$$M_F(x_1,x_2)=-M_F(x_2,x_1),{\tilde{M}}_F(x_1,x_2)={\tilde{M}}_F(x_2,x_1).$$

$${\int }_0^{1}dx(L_{3q}(x)+L_{3g}(x))=0.$$

$${\int }_0^{1}dx{L}_{3g}(x)\equiv L_{\mathrm{pot}}=-{\int }_0^{1}dx{L}_{3q}(x).$$

$$L_{\mathrm{Ji}}^q-L_{\mathrm{JM}}^q=L_{\mathrm{pot}}.$$

$$L_q^n$$

$$L_g^n$$

$$\frac{\partial }{\partial t}{\left(\begin{array}{c}{L}_q^n\\ L_g^n\end{array}\right)}^{WW}$$

$$t=-\frac{2}{b_0}\ln \frac{{\alpha }_s({\mu }^2)}{{\alpha }_s({\mu }_0^2)},$$

$$\frac{\partial }{\partial t}\left(\begin{array}{c}{L}_{3q}^n\\ L_{3g}^n\end{array}\right)\stackrel{?}{=}\left(\begin{array}{cc}{A}_n& B_n\\ C_n& D_n\end{array}\right)\left(\begin{array}{c}{L}_{3q}^n\\ L_{3g}^n\end{array}\right),$$

$$\frac{\partial }{\partial t}\left(\begin{array}{c}{L}_q^n\\ L_g^n\end{array}\right)$$

$$T_{\bar{q}Fq}({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3),\Delta T_{\bar{q}Fq}({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3),T_{3F}^{+}({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3),\Delta T_{3F}^{+}({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3).$$

$$x_1={\tilde{x}}_3,x_2-x_1={\tilde{x}}_2,x_2=-{\tilde{x}}_1,{\tilde{x}}_1+{\tilde{x}}_2+{\tilde{x}}_3=0.$$

$${\Phi }_F(x_1,x_2)\leftrightarrow 2T_{\bar{q}Fq}(-x_2,x_2-x_1,x_1),$$

$${\tilde{\Phi }}_F(x_1,x_2)\leftrightarrow -2\Delta T_{\bar{q}Fq}(-x_2,x_2-x_1,x_1),$$

$$M_F(x_1,x_2)\leftrightarrow T_{3F}^{+}(-x_2,x_2-x_1,x_1),$$

$${\tilde{M}}_F(x_1,x_2)\leftrightarrow -\Delta T_{3F}^{+}({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3).$$

$${\tilde{M}}_F(x_1,x_2)=-M_F(x_2-x_1,x_2)+M_F(x_1,x_1-x_2).$$

$${\mathfrak{S}}^{\pm }({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3),{\mathfrak{F}}^{\pm }({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3)$$

$$2{\mathfrak{S}}^{+}({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3)$$

$${\mathfrak{F}}^{+}({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3)$$

$${\mathfrak{S}}^{+}(-\tilde{x})={\mathfrak{S}}^{+}(\tilde{x}),{\mathfrak{F}}^{+}(-\tilde{x})=-{\mathfrak{F}}^{+}(\tilde{x}).$$

Full text

I Introduction

After more than a decade of experimental programs at RHIC, HERMES, COMPASS and JLab Adare:2008aa ; Airapetian:2010ac ; Adamczyk:2012qj ; Alekseev:2010ub ; Adamczyk:2014ozi ; Prok:2014ltt , and concurrent theoretical efforts based on global QCD analysis deFlorian:2014yva ; Nocera:2014gqa ; Sato:2016tuz , we now know that the gluon helicity contribution $\Delta G$ to the Jaffe-Manohar sum rule of the nucleon spin Jaffe:1989jz

[TABLE]

is nonzero and can be significant. Together with the well-known quark helicity contribution $\Delta\Sigma$, the first two terms of (1) account for $\sim 70\%$ of the total spin, leaving some room for the orbital angular momentum (OAM) $L_{q,g}$ of quarks and gluons. However, this is at best a tentative conclusion. The present estimate of $\Delta G$ is based on the RHIC 200 GeV $pp$ data which can only access a limited range of the Bjorken variable $x$. Uncertainties coming from the small $x$-region $x\lesssim 0.05$ are estimated to be huge, and might completely alter the current estimates. In such circumstances, an obvious direction to proceed for the QCD spin community is to better constrain $\Delta G(x)$ in the small-$x$ region. Indeed, new data from RHIC 510 GeV $pp$ runs have just begun to appear Adam:2019aml . Moreover, the precise determination of $\Delta G$ has been designated as one of the primary goals of the future Electron-Ion Collider (EIC) Accardi:2012qut .

An equally obvious direction is to directly constrain the OAM contributions $L_{q,g}$. Unfortunately, however, this is currently not pursued experimentally due to the lack of clean observables from which one can systematically extract $L_{q,g}$ as the moment of the corresponding $x$-distributions $L_{q,g}=\int_{0}^{1}dxL_{q,g}(x)$. In fact, it is only relatively recently that the rigorous definitions of $L_{q,g}(x)$ have been derived Hatta:2012cs (following the earlier developments Lorce:2011kd ; Hatta:2011ku ; Lorce:2011ni ), although their gauge non-invariant versions have been known for quite some time Harindranath:1998ve ; Hagler:1998kg . Such a progress has finally led theorists to identify a few experimental processes that are sensitive to $L_{q,g}(x)$ Courtoy:2013oaa ; Ji:2016jgn ; Hatta:2016aoc ; Bhattacharya:2017bvs ; Bhattacharya:2018lgm . In these works, observables have been computed to leading order, whereas a realistic global analysis requires the knowledge of the next-to-leading order (NLO) corrections together with the ‘DGLAP’ equation for $L_{q,g}(x)$ to NLO, or at least to leading order.

However, so far the evolution of $L_{q,g}(x)$ has not been fully understood even to leading order. This is because $L_{q}(x)$ and $L_{g}(x)$ are not the usual twist-two parton distribution functions (PDFs). They consist of the Wandzura-Wilczek (WW) part and the genuine twist-three part

[TABLE]

The WW part is related to the twist-two unpolarized and polarized PDFs. As such, its evolution is known to the same order as that of the corresponding twist-two distributions (namely, to three-loops). The LO evolution equation derived in the early literature Harindranath:1998ve ; Hagler:1998kg can be fully understood in this way Hoodbhoy:1998yb ; Boussarie:2019icw . More recently, the small-$x$ behavior of the WW part has received a lot of attention, mainly motivated by the aforementioned uncertainties of $\Delta G$ at small-$x$ Hatta:2018itc ; More:2017zqp ; Kovchegov:2019rrz ; Boussarie:2019icw . An analysis based on the double logarithmic resummation Boussarie:2019icw has shown that there is a significant cancellation between $\Delta G(x)$ and $L_{g}(x)$ in the small-$x$ region (see, however, Kovchegov:2019rrz ).

On the other hand, the genuine twist-three part $L_{3q,3g}(x)$ is given by the matrix elements of three-parton ($q\bar{q}g,ggg$), twist-three operators. The QCD evolution of such twist-three correlators is complicated even to one-loop, and to our knowledge, it has never been discussed in the context of OAM. Fortunately, however, it has been discussed in a different context—QCD evolution for the Efremov-Teryaev-Qiu-Sterman (ETQS) function Efremov:1984ip ; Qiu:1991wg for transverse single spin asymmetry Kang:2008ey ; Vogelsang:2009pj ; Braun:2009mi . Moreover, a C++ program to numerically solve this evolution has been developed by Pirnay Pirnay:2013fra . These results can be adapted for the present purpose with appropriate modifications. Thus, the goal of this paper is to articulate the evolution equation for $L_{3q,3g}(x)$ and demonstrate the feasibility of numerically solving this equation.

The paper is organized as follows: In Sec. II we will explain the quark and gluon OAM distributions and their evolutions. Then in Sec. III the evolution of the genuine twist-three part will be discussed. We will show numerical results on the evolution of the genuine twist-three part in Sec. IV. We will discuss the cusp anomalous dimensions for the large-$n$ moments and the scale dependence of the potential angular momentum. Finally, conclusions will be drawn in Sec. V.

II OAM distributions and their evolution equation

In this section, we first review the formula (2) derived in Hatta:2012cs and then discuss its evolution equation. We employ the following normalization for the singlet quark (and antiquark) OAM distribution $L_{q}(x)=\sum_{f}(L_{f}(x)+L_{\bar{f}}(x))$ and the gluon OAM distribution $L_{g}(x)$

[TABLE]

Our normalization here is different from the one in Hatta:2012cs , but it is the same as in Hatta:2018itc ; Boussarie:2019icw . In the present normalization and notation, the result of Hatta:2012cs reads

[TABLE]

[TABLE]

where $X=\frac{x_{1}+x_{2}}{2}$ and $x^{\prime}$ in the second line of (5) means $x_{1}-x_{2}$. The first line on the right hand side is the WW part $L_{q,g}^{WW}(x)$ and the rest is the genuine twist-three part $L_{3q,3g}(x)$. $\Sigma(x)=\sum_{f}(f(x)+\bar{f}(x))$ and $G(x)$ are the usual, unpolarized quark and gluon distributions. (Ref. Hatta:2012cs used the notation $H_{g}=xG$.) $E_{q}(x)=\sum_{f}(E_{f}(x)+E_{\bar{f}}(x))$ and $E_{g}(x)$ are the helicity-flip quark and gluon generalized parton distributions (GPDs) commonly called ‘GPD $E$’.

The genuine twist-three distributions $\Phi_{F}$ and $\widetilde{\Phi}_{F}$ are defined through the matrix elements of quark-gluon correlation functions along the light-cone Hatta:2012cs ; Ji:2012ba

[TABLE]

[TABLE]

where $S^{\mu}=\delta^{\mu}_{+}S^{+}$ is the longitudinally polarized nucleon spin vector and $n^{\mu}=\delta^{\mu}_{-}$ is a fixed light-like vector. $i,j=1,2$ denote transverse indices. The Wilson lines are omitted for simplicity. The summation over quark flavors is implied. The momentum transfer $\Delta^{\mu}=P^{\prime\mu}-P^{\mu}=\delta^{\mu}_{i}\Delta^{i}$ is assumed to be small, has only transverse components, and we have kept only the linear terms in $\Delta$. The variables $1>x_{1,2}>-1$ can be interpreted as follows: The outgoing quark ($\psi$) and gluon ($F^{+i}$) carry momentum fractions $x_{1}$ and $x_{2}-x_{1}$ of the parent nucleon, respectively, and the returning quark ($\bar{\psi}$) has momentum fraction $x_{2}$, see Fig. 1. Similarly, for the three-gluon correlators $M_{F}$ and $\widetilde{M}_{F}$, one has

[TABLE]

[TABLE]

where $FFF$ means $F_{a}(T^{b})_{ac}F_{b}F_{c}=if_{abc}F_{a}F_{b}F_{c}$ in color space. These distributions have the following symmetry properties

[TABLE]

Note that the total (quark plus gluon) genuine twist-three OAM distribution integrates to zero

[TABLE]

We can write this as

[TABLE]

The quantity $L_{\mathrm{pot}}$, sometimes called the potential angular momentum Wakamatsu:2010qj , represents the difference between the kinetic (Ji) OAM Ji:1996ek and the canonical (Jaffe-Manohar) OAM Jaffe:1989jz

[TABLE]

Roughly speaking, it arises from the difference between the covariant derivative (kinetic momentum) $\vec{x}\times\vec{D}$ and the partial derivative (canonical momentum) $\vec{x}\times\vec{\partial}$ in the definition of OAM. The value of $L_{\mathrm{pot}}$ for an electron in QED has been the subject of debate in the literature Burkardt:2008ua ; Liu:2014fxa ; Ji:2015sio . More recently, $L_{\mathrm{pot}}$ has been calculated in lattice QCD and found to be nonzero Engelhardt:2017miy ; Engelhardt:2019lyy .

Let us now consider the evolution of $L_{q,g}(x)$. Taking the $n$-th moments $L_{q,g}^{n}=\int_{0}^{1}dxx^{n-1}L_{q,g}(x)$, etc., we find

[TABLE]

where $L_{3q}$ and $L_{3g}$ are the contributions from the genuine twist three part. The evolution of $\Sigma^{n},G^{n},\Delta\Sigma^{n}$ and $\Delta G^{n}$ in the renormalization scale $\mu^{2}$ is governed by the standard DGLAP anomalous dimensions $\gamma^{n}$ and $\Delta\gamma^{n}$. Moreover, the anomalous dimensions of $E^{n}_{q},E^{n}_{g}$ are the same as those for $\Sigma^{n},G^{n}$. Therefore, if one restricts oneself to the WW part, one can immediately write down the evolution equation

[TABLE]

where

[TABLE]

with $\alpha_{s}(\mu^{2})=\frac{4\pi}{b_{0}\ln(\mu^{2}/\Lambda_{QCD}^{2})}$ and $b_{0}=\frac{11N_{c}}{3}-\frac{2n_{f}}{3}$. In the WW approximation, (15) holds to all orders in perturbation theory Boussarie:2019icw . To one-loop order, it agrees with the result obtained in Harindranath:1998ve ; Hagler:1998kg , see also, Hoodbhoy:1998yb .

Inclusion of the genuine twist-three parts makes things considerably more complicated. One might naively expect that the genuine twist-three part would evolve with its own anomalous dimensions

[TABLE]

and therefore the evolution of the full OAMs would be given by

[TABLE]

However, we shall demonstrate later that this is not the case. In general, different moments (different $n$’s) mix under evolution so that the right hand side of (17) should involve a summation over all moments. This suggests that it is more convenient to study the evolution directly in the $x$-space.

III Evolution of the genuine-twist three part

At first sight, the $\mu^{2}$-evolution of $\Phi_{F}(x_{1},x_{2})$ etc., hence that of $L_{3q,3g}(x)$ seems a challenging open question. However, it is actually known in the context of transverse single spin asymmetry (SSA). There, exactly the same set of operators as in (6)-(9) appear. The difference is that in the case of SSA, one takes the forward matrix element (i.e., $\Delta=0$) in the transversely polarized nucleon state $S^{\mu}=\delta^{\mu}_{i}S_{i}$. The resulting distribution, the Efremov-Teryaev-Qiu-Sterman (ETQS) function Efremov:1984ip ; Qiu:1991wg has been extensively discussed in the literature. Because the evolution is intrinsic to the operators involved, not to external states, and here we only consider the $\Delta\to 0$ limit of the nonforward matrix element,111 In principle, $\Phi_{F}$ etc. depend on $\Delta^{2}$ but this dependence has been neglected in (6) since it is of higher order. the same evolution equation should apply. The complete derivation of the evolution equation has been given in Braun:2009mi following earlier attempts Kang:2008ey ; Vogelsang:2009pj (see also Schafer:2012ra ; Ma:2012xn ; Yoshida:2016tfh ). The result is too lengthy to be reproduced here, but fortunately, a C++ code is publicly available Pirnay:2013fra . We shall heavily rely on this code in the following.

For this purpose, first we need to clarify the difference in notations between ours and in Braun:2009mi ; Pirnay:2013fra . There the authors introduced the following distributions with three arguments

[TABLE]

To avoid confusion, we have added a tilde on momentum fractions. Their meaning is self-explanatory from Fig. 1 (beware of the direction of arrows) with the correspondence

[TABLE]

The four distributions in (19) are direct analogs of (6)-(9). The plus sign on $T_{3F}$ and $\Delta T_{3F}$ means the contraction of color indices with the $f$-symbol, cf., the comment below (9). We have carefully checked the relative normalization and found the correspondence222 This comparison is complicated by the fact that Refs. Hatta:2012cs ; Braun:2009mi use different conventions for $\gamma_{5}$ and the epsilon tensor. Ref. Hatta:2012cs used $\epsilon^{0123}=+1$ whereas Ref. Braun:2009mi used $\epsilon_{0123}=-\epsilon^{0123}=1$. The sign of $\gamma_{5}$ is also opposite in the two references.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

As observed in Braun:2009mi , $\Delta T^{+}_{3F}$ is not an independent function (see Eq. (32) there). Accordingly, $\widetilde{M}_{F}$ can be eliminated via the corresponding formula

[TABLE]

This was not noticed in Hatta:2012cs , but it can be indeed derived from the formulas in this reference.

Another complication is that in Ref. Braun:2009mi , the evolution equation has been presented in terms of

[TABLE]

which are particular linear combinations of the distributions in (19). The signs $\pm$ refer to $C$-parity, and the point is that the evolution equation does not mix the $C$-even and $C$-odd functions. From (21)-(24), we find the correspondence

[TABLE]

[TABLE]

The distributions $\mathfrak{S}^{+}$ and $\mathfrak{F}^{+}$ have the following symmetry Braun:2009mi

[TABLE]

In our case, we find from (10)333 In the three-gluon sector, there is another relation

$\displaystyle F^{+}(x_{1},x_{2})=-F^{+}(x_{2}-x_{1},x_{2})\,.$

(30) This has the same sign as the corresponding relation $\mathfrak{F}^{+}(\tilde{x}_{1},\tilde{x}_{2},\tilde{x}_{3})=-\mathfrak{F}^{+}(\tilde{x}_{1},\tilde{x}_{3},\tilde{x}_{2})$ in Braun:2009mi (see the unnumbered equation below (30)). There is no contradiction here.

[TABLE]

Notice the differences in sign, which is simply due to the fact that we are dealing with different matrix elements (longitudinally polarized nucleon, finite momentum transfer). This however causes no problem in practice because the evolution equation automatically preserves the symmetry property of the initial conditions.

Note also that in (27) and (28), we only showed the $C$-even distributions. This is because the OAM distributions are $C$-even ($PT$-even) Hatta:2011ku . Indeed, the twist-three part of (4) and (5) can be solely written in terms of $C$-even functions $S^{+}$ and $F^{+}$. It is not difficult to show that

[TABLE]

In the gluon case, let us write444Note that we can replace $\Phi_{F}(x_{1},x_{2})-\Phi_{F}(-x_{1},-x_{2})$ with $2S^{+}(x_{1},x_{2})$ because the $\widetilde{\Phi}_{F}$ part of $S^{+}$ drops out $\int dX(\widetilde{\Phi}_{F}(X,x)-\widetilde{\Phi}_{F}(-X,-x))=0$.

[TABLE]

where

[TABLE]

In terms of moments, we have

[TABLE]

[TABLE]

Eqs. (32)-(37) are the starting point of our numerical analysis which we now turn to.

IV Numerical results

In this section, we present the one-loop scale evolution of the genuine twist-three part of the OAM distributions $L_{3q,3g}(x,\mu)$. Our numerical results are based on a C++ code developed by Pirnay Pirnay:2013fra .555Here we clarify a few things and list a few typos in the code:

The code in Ref. Pirnay:2013fra is based on the evolution kernels derived in Ref. Braun:2009mi . Ref. Braun:2009mi calculates the kernels with the evolution equation defined as $\frac{\partial}{\partial\ln\mu}\cdots$. The code in Ref. Pirnay:2013fra also implements the evolution equation numerically as $\frac{\partial}{\partial\ln\mu}\cdots$. But the manuscript of Ref. Pirnay:2013fra writes the evolution equation as $\frac{\partial}{\partial\ln\mu^{2}}\cdots$. The expression (52) in Ref. Braun:2009mi has a typo: the factor $\frac{\alpha_{s}}{4\pi}$ in the third term should be $\frac{\alpha_{s}}{2\pi}$.

In the line 347 of the code file ‘t3evol.cpp’, the ‘F+_initial.txt’ should be ‘F-_initial.txt’.

In the lines 65, 90, 117 and 142 of the code file ‘fffkernels.h’, there should be no ‘Nc’ multiplying ‘beta0(nf)’.

We would like to thank Vladimir Braun and Alexander Manashov for confirming these. Although the code is intended for the ETQS function, it can be straightforwardly adapted to the present problem with almost no change. The only thing one should keep in mind is that the symmetry properties of the various distributions are different, see (29) and (30). However, this does not cause extra complications since the evolution equation preserves this symmetry, and one just needs to properly implement the relevant symmetry in the initial conditions.

We set the number of lattice sites in the interval $-1\leq x\leq 1$ to be $N_{x}=2001$ (this has to be an odd number in the code), so that the lattice size is $\Delta x=0.001$. As for the initial conditions, currently nothing is known about the functional forms of the genuine twist-three distributions $\Phi(x_{1},x_{2})$, etc., except that they must vanish at kinematic boundaries (e.g., $x_{1}=\pm 1$) and obey the symmetry properties (10). We thus introduce three different models at the initial scale $\mu_{0}^{2}$

[TABLE]

with

$f(x_{1},x_{2})=\frac{1}{10}(1-x_{1}^{2})(1-x_{2}^{2})(1-(x_{1}-x_{2})^{2})$ ; 2. 2.

$f(x_{1},x_{2})=x_{1}x_{2}(1-x_{1}^{2})(1-x_{2}^{2})(1-(x_{1}-x_{2})^{2})$ ; 3. 3.

$f(x_{1},x_{2})=x^{2}_{1}x^{2}_{2}(1-x_{1}^{2})(1-x_{2}^{2})(1-(x_{1}-x_{2})^{2})$ .

The overall sign has been fixed such that the potential angular momentum (12) calculated from these initial conditions is positive $L_{\mathrm{pot}}>0$, as suggested by an intuitive argument Burkardt:2012sd as well as a recent lattice calculation Engelhardt:2017miy . We then use (27) and (28) to convert the above initial conditions into those of $S^{+}$ and $F^{+}$. They are further converted into $\mathfrak{S}^{+}(\tilde{x}_{1},\tilde{x}_{2},\tilde{x}_{3})=S^{+}(x_{1},x_{2})$ and $\mathfrak{F}^{+}(\tilde{x}_{1},\tilde{x}_{2},\tilde{x}_{3})=F^{+}(x_{1},x_{2})$ since the code requires $\mathfrak{S}^{+}$ and $\mathfrak{F}^{+}$ as the initial conditions. In the end of the numerical calculation, we convert the final results of $\mathfrak{S}^{+}$ and $\mathfrak{F}^{+}$ back into $S^{+}$ and $F^{+}$. Finally, we compute the twist-three OAM distributions $L_{3q}(x)$ and $L_{3g}(x)$, Eqs. (32)-(35), numerically by using the Simpson’s rule. The integrands of $L_{3q}(x)$ and $L_{3g}(x)$ may have singularities when $x_{1}\to 0$ or $x^{\prime}\to 0$. Such singularities may be present already in the initial conditions, or can be dynamically generated by the evolution (see below). To improve the numerical accuracy when this happens, we first fit the integrand $f(x_{1})$ (or $f(x^{\prime})$) in the small-$x$ region (in practice, $x<0.005$) by the power function $ax_{1}^{b}$ (or $ax^{\prime b}$) with $b<0$. Then we split as

[TABLE]

The integrand of the first term has no singularities so the integral can be evaluated accurately. The second integral can be done analytically.

As a warm-up, let us first check the sum rules Hatta:2012cs :

[TABLE]

This is a more detailed version of (11) (see (33)). We numerically evolve each of the initial conditions from $\mu_{0}^{2}=1$ GeV2 to $\mu^{2}=10$ GeV2. The numerical results of the left-handed-sides of (42) and (43) for the three different initial conditions are shown in Table 1. As a reference, we also show the value of the potential angular momentum $L_{\mathrm{pot}}$ (12) in these models. We see that the sum rules are satisfied up to numerical errors666The default time step in $t$ (see (16)) in the evolution code is $\Delta t=0.01$. However, we find it necessary to make this 0.001 for the initial condition 1, because a singularity is developed at small-$x$ for this initial condition and we need a higher accuracy to stabilize the evolution and satisfy the sum rules. We also use $0.001$ as the time step for initial conditions 2 and 3..

Next, we study whether the moments of the twist-three OAM distributions $L^{n}_{3q}$ (flavor singlet) and $L^{n}_{3g}$ satisfy the homogeneous equation (17). This can be tested by assuming (17) and considering an infinitesimal evolution

[TABLE]

In practice, we use $\mu_{0}^{2}=1$ GeV2 and $\mu^{2}=1.01$ GeV2. This means that $n_{f}=3$, and the code uses the value $\Lambda_{QCD}=0.204$ GeV in this case. In order to extract $A_{n},B_{n},...$, it is not enough to consider one initial condition since (44) contains two equations but there are four unknowns. But we can combine any two sets of initial conditions and form a solvable system of linear equations. We will label the results of $A_{n},B_{n},...$ obtained from the $i$-th and $j$-th initial conditions as $A_{n}(ij),B_{n}(ij)$,… and compare these with the known anomalous dimensions of the twist-two operators

[TABLE]

We make this comparison in the (naive) hope that there may be a significant cancelation in the second line of (18) so that the evolution equation for the total OAM distributions is formally the same as that in the WW approximation. Note that we set $n^{\prime}_{f}=2$ in (45) because this is the value used in the one-loop evolution kernel in the code. It is different from $n_{f}$ used in the running coupling constant. The latter depends on the scale $\mu^{2}$.

The comparison results are shown in Fig. 2. The numerically extracted $A_{n}$, $B_{n}$, $C_{n}$ and $D_{n}$ of the twist-three operators are not universal because they depend on the choice of the initial condition. This implies that (17) is not valid and the moments of different $n$’s will mix with each other during the evolution. Interestingly, however, we find that the large-$n$ behaviors of the $A_{n}$, $B_{n}$, $C_{n}$ and $D_{n}$ are given by

[TABLE]

which agree with the $n$-dependence of $\gamma^{n+1}$’s when $n$ is large. The constant terms in the $A_{n}$, $B_{n}$, $C_{n}$ and $D_{n}$ are not universal and depend on the choice of the initial conditions. The leading, logarithmic terms are called the cusp anomalous dimensions. We have thus found that, when $n$ is large, (17) is approximately valid and the second line of (18) is small. This means that when $x\sim 1$, the total OAM distribution approximately satisfies the same evolution equation as in the WW approximation. We note that the emergence of the cusp anomalous dimension in the twist-three distributions in the large-$x$ limit was pointed out in Braun:2009mi in a different context. Incidentally, we found the correct asymptotic behavior (46) only after fixing a typo in the code mentioned in Footnote 5 (the factor of $N_{c}$).

We then study the evolution from $\mu_{0}^{2}=1$ GeV2 to $\mu^{2}=10$ GeV2. The initial and final twist-three OAM distributions are plotted in Fig. 3. We note that both $L_{3q}(x)$ (flavor singlet) and $L_{3g}(x)$ develop singular behaviors at small-$x$ during the evolution, even though their initial conditions are regular at $x=0$. We fit the singular behavior by the power function $ax^{b}$ with $b<0$ and find that $b$ depends on the choice of initial conditions. But it always satisfies $b>-1$ in our three sets of initial conditions, which guarantees that both $L_{3q}(x)$ and $L_{3g}(x)$ are integrable in the range $0\leq x\leq 1$. The situation is somewhat similar to the DGLAP evolution of the WW part Hatta:2018itc , where a singularity is developed from nonsingular initial conditions. In that case it was possible to gain some analytical insights, but in the present case a similar analytical study is difficult due to the complexity of the twist-three evolution. We also note that, curiously, the three-gluon part $L_{3g,F}(x)$ of $L_{3g}(x)$, see (33), evolves very weakly with the renormalization scale. Since this part does not contribute to the integrated OAM (42), we suspect that in general it plays a minor role in the nucleon spin decomposition. We however note that only for the initial condition 1, $L_{3g,F}(x)$ shows a sharp drop in the smallest $x$ bins that can be fitted by a power law as shown in Fig. 4. In fact, such a rapid behavior is necessary to satisfy the sum rule in this case.

Finally we study the $\mu^{2}$-dependence of the potential angular momentum (12). Here we are also interested in the flavor nonsinglet part $L_{\mathrm{pot}}^{u-d}$ because the recent lattice calculations have shown that this quantity is positive and large Engelhardt:2017miy ; Engelhardt:2019lyy . In fact, the singlet part $L^{u+d}_{\mathrm{pot}}$ is numerically much smaller due to a large cancelation between the $u$- and $d$-quark contributions. We have not mentioned the evolution in the flavor nonsinglet sector so far, but of course it is simpler than the singlet case because there is no mixing with the three-gluon part. The code can handle this as well.

We evolve the initial conditions 2 and 3 from $\mu_{0}^{2}=2$ GeV2 to $\mu^{2}=12$ GeV2 and compute how $L_{\mathrm{pot}}$ changes with the scale for both the flavor singlet and nonsinglet parts. The initial and final scales are chosen such that the number of dynamical quarks $n_{f}$ in the running coupling does not change in the course of evolution. The results are shown in the left panel of Fig. 5. In both cases, the evolution tends to suppress the magnitude of $L_{\mathrm{pot}}$, with a stronger suppression in the singlet case. This can be understood as arising from the mixing with the three-gluon correlator. We have fitted these results in the form

[TABLE]

which features an ‘anomalous dimension’ $\gamma$. Here $\Lambda_{QCD}=0.175$ GeV is used in the code for $n_{f}=4$. We find that $\gamma$ is $\mu$-dependent as shown in the right panel of Fig. 5. It also depends on the initial condition. This means that the evolution of $L_{\mathrm{pot}}$ cannot be characterized by a single anomalous dimension even in the nonsinglet sector due to the mixing between different moments. From the obtained behavior of $\gamma$, we deduce that $L_{\mathrm{pot}}\to 0$ as $\mu\to\infty$ (though this has to be checked more carefully). Therefore, the asymptotic result for the integrated OAM

[TABLE]

obtained in the WW approximation Ji:1995cu will not be affected. The subleading corrections to (48) have the $\mu$-dependence of the form (47) with $\gamma=\frac{2(16+3n_{f})}{9b_{0}}\approx 0.75$ (for $n_{f}=4$) Ji:1995cu . Since this is smaller than the value shown in Fig. 5, it seems that $L_{\mathrm{pot}}$ can be neglected in the large-$\mu$ region.

In Ji:2015sio , it has been shown that $L_{\mathrm{pot}}$ vanishes for a single electron state to one-loop order in QED. This is not necessarily in contradiction to the present result. The operator $\vec{x}\times\bar{\psi}\gamma^{+}\vec{A}\psi$ (in the light-cone gauge) relevant to $L_{\mathrm{pot}}$ mixes with a tower of operators $\vec{x}\times\bar{\psi}\gamma^{+}D^{\mu}D^{\nu}\cdots\vec{A}\psi$ and if these higher moments have nonvanishing initial values (as in our initial conditions), they will affect the renormalization of $L_{\mathrm{pot}}$.

V Conclusions

In this paper, we have studied, for the first time, the one-loop QCD evolution of the genuine twist-three part of the OAM distributions. In particular, the scale variation of the potential angular momentum has been demonstrated. As anticipated by the complicated relations between $L_{3q,3g}(x)$ and the underlying distributions $S^{+}$ and $F^{+}$ (see, e.g., (32)), different moments of $L_{3q,3g}(x)$ mix under evolution except in the large-$n$ limit. This suggests that it is more convenient to look at the evolution directly in the $x$-space. Together with the known evolution of the WW part, the one-loop evolution of the total OAM distributions $L_{q,g}(x)$ is now fully under control and ready for phenomenological applications.

In the present C++ code, the grid in $x$ is uniform and the smallest value of $x$ that we achieved is 0.001. It is not realistic to go to much lower values of $x$ because it is computationally too expensive. In view of the recent controversy regarding the small-$x$ asymptotic behavior of the OAM distributions Hatta:2018itc ; Kovchegov:2019rrz ; Boussarie:2019icw , it would be interesting to modify the code (e.g., set up a grid in $\ln 1/x$ instead of $x$) to zoom in on to the small-$x$ region as was done for the WW part Hatta:2018itc . We leave this to future work.

Acknowledgments

We thank Renaud Boussarie, Vladimir Braun, Alexander Manashov and Werner Vogelsang for discussions. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract number DE-SC0012704. It is also supported by the LDRD program of Brookhaven National Laboratory. X.Y. is supported by U.S. Department of Energy research grant DE-FG02-05ER41367 and Brookhaven National Laboratory.