Impact of DCSB and dynamical diquark correlations on proton GPDs

TL;DR

This paper computes the proton's generalized parton distributions using a QCD-inspired NJL model that incorporates dynamical chiral symmetry breaking and diquark correlations, providing insights into proton structure and form factors.

Contribution

It introduces a consistent calculation of proton GPDs within the NJL model including DCSB and diquark correlations, ensuring Lorentz covariance and polynomiality sum rules.

Findings

D-term of -1.08 with dressed correlator

Proton stability linked to dressed constituent quarks

Vanishing anomalous gravitomagnetic moment

Abstract

We calculate the leading-twist, helicity-independent generalized parton distributions (GPDs) of the proton, at finite skewness, in the Nambu--Jona-Lasinio (NJL) model of quantum chomodynamics (QCD). The NJL model reproduces low-energy characteristics of QCD, including dynamical chiral symmetry breaking (DCSB). The proton bound-state amplitude is solved for using the Faddeev equation in a quark-diquark approximation, including both dynamical scalar and axial vector diquarks. GPDs are calculated using a dressed non-local correlator, consistent with DCSB, which is obtained by solving a Bethe-Salpeter equation. The model and approximations used observe Lorentz covariance, and as a consequence the GPDs obey polynomiality sum rules. Extractions of electromagnetic and gravitational form factors are performed. We find a D-term of $-1.08$ when the non-local correlator is properly dressed, and $0.85$ when the bare correlator is used instead, suggesting that within this framework proton stability requires the constituent quarks to be dressed consistently with DCSB. We also find that the anomalous gravitomagnetic vanishes, as required by Poincar\'{e} symmetry.

Equations

$$V_{\lambda {\lambda }^{\prime }}^q$$

$$V_{\lambda {\lambda }^{\prime }}^q=\bar{u}(p^{\prime },{\lambda }^{\prime })\left[/nH(x,\xi ,t)+\frac{i{\sigma }^{n\Delta }}{2M_N}E(x,\xi ,t)\right]u(p,\lambda ),$$

$$H_{X,i}(x,\xi ,t)=\sum _j\int \frac{\mathrm{d}y}{|y|}{h}_{Y/X,ij}(y,\xi ,t)H_{Y,j}\left(\frac{x}{y},\frac{\xi }{y},t\right).$$

$$\tau (p)\approx {\tau }_{\mathrm{pole}}(p)=\frac{1}{p^2-M_{\mathrm{dq}}^2+i0},$$

$$v=p^2-M_{\mathrm{dq}}^2,v^{\prime }=p^{\prime 2}-M_{\mathrm{dq}}^2.$$

$$\tau (p^{\prime })\Lambda (x,\xi ,t,v,v^{\prime })\tau (p).$$

$$\tau (p^{\prime })\Lambda (x,\xi ,t,v,v^{\prime })\tau (p)\approx {\tau }_{\mathrm{pole}}(p^{\prime })\Lambda (x,\xi ,t,0,0){\tau }_{\mathrm{pole}}(p),$$

$$V_s^q=H_s^q(x,\xi ,t).$$

$$\begin{gathered} V_{a,{\lambda }^{\prime }\lambda }^q=-(ϵ{ϵ}^{\prime \ast })H_{1a}^q+\frac{({ϵ}^{\prime \ast }n)(ϵ\Delta )-(ϵn)({ϵ}^{\prime \ast }\Delta )}{2(Pn)}{H}_{2a}^q \\ +\frac{(ϵ\Delta )({ϵ}^{\prime \ast }\Delta )}{2M^2}{H}_{3a}^q-\frac{({ϵ}^{\prime \ast }n)(ϵ\Delta )+(ϵn)({ϵ}^{\prime \ast }\Delta )}{2(Pn)}{H}_{4a}^q \\ +\left[\frac{M^2(ϵn)({ϵ}^{\prime \ast }n)}{(Pn{)}^2}+\frac{1}{3}(ϵ{ϵ}^{\prime \ast })\right]H_{5a}^q, \end{gathered}$$

$$\begin{gathered} V_{a,{\lambda }^{\prime }\lambda }^{q,\mu \nu }=-g^{\mu \nu }{H}_{1a}^q+\frac{n^{\nu }{\Delta }^{\mu }-n^{\mu }{\Delta }^{\nu }}{2(Pn)}{H}_{2a}^q+\frac{{\Delta }^{\mu }{\Delta }^{\nu }}{2M^2}{H}_{3a}^q \\ -\frac{n^{\nu }{\Delta }^{\mu }+n^{\mu }{\Delta }^{\nu }}{2(Pn)}{H}_{4a}^q+\left[\frac{M^2{n}^{\mu }{n}^{\nu }}{(Pn{)}^2}+\frac{1}{3}{g}^{\mu \nu }\right]H_{5a}^q. \end{gathered}$$

$$V_{sa}^q=\frac{1}{(Pn)}\frac{1}{M_s+M_a}i{ϵ}^{P\Delta nϵ}{H}_{sa}^q(x,\xi ,t).$$

$$H_{I=0}(x,\xi ,t)$$

$$H_{I=1}(x,\xi ,t)$$

$$H_I(x,\xi ,t)$$

$$E_I(x,\xi ,t)$$

$$H_{I=0,1}^{\prime }(x,\xi ,t)$$

$$D_Q(x,\xi ,t)$$

$${\mathcal{E}}_1\left(z,t\right)=E_1\left(\frac{4M^2-t(1-z^2)}{4{\Lambda }_{\mathrm{UV}}^2}\right)-E_1\left(\frac{4M^2-t(1-z^2)}{4{\Lambda }_{\mathrm{IR}}^2}\right),$$

$$-|\xi |\le x\le \max (1,|\xi |),$$

$$\xi =\frac{(np)-(n{p}^{\prime })}{(np)+(n{p}^{\prime })}.$$

$$\int \mathrm{d}x{H}^q(x,\xi ,t)=F_1^q(t),\int \mathrm{d}x{E}^q(x,\xi ,t)=F_2^q(t),$$

$$\int \mathrm{d}xx{H}^q(x,\xi ,t)$$

$$\int \mathrm{d}xx{E}^q(x,\xi ,t)$$

Full text

I Introduction

Generalized parton distributions (GPDs) appear in the calculation of hard exclusive reactions such as deeply virtual Compton scattering (DVCS) and deeply virtual meson production (DVMP). Factorization Collins et al. (1997); Collins and Freund (1999); Ji and Osborne (1998) allows the amplitudes of these processes to be broken down (up to power-suppressed corrections) into the convolution of a hard scattering kernel and a soft matrix element of quark and/or gluon fields which contains the GPDs.

GPDs are Lorentz-invariant functions of four variables, including the renormalization scale, that describe many interesting properties of hadrons. They encode spatial light cone distributions of partons through two-dimensional Fourier transforms Burkardt (2003). Additionally, their Mellin moments encode the electromagnetic and gravitational properties of hadrons—allowing, in the latter case, for such properties to be studied through hard exclusive reactions, in lieu of the impossibility of graviton-exchange experiments. The gravitational properties shed light on the way that mass and angular momentum are distributed among the quarks and gluons within the hadron, thus directly addressing deep questions about the mass Lorcé et al. (2019); Hatta et al. (2018) and spin Leader and Lorcé (2014) decompositions of the proton.

In light of these considerations, calculations of proton GPDs are to be highly desired. It is vital that any model calculation observe the symmetries and low-energy properties of quantum chromodynamics (QCD), so that the qualitative and quantitative effects of these phenomena manifest in the GPDs themselves. For instance, the relationship between Mellin moments and electromagnetic and gravitational form factors is a consequence of Lorentz covariance Ji (1997a), and the magnitude of the gravitational form factors has an intimate relationship with dynamical chiral symmetry breaking (DCSB) Polyakov and Schweitzer (2018); Freese and Cloët (2019). For this reason, we use the Nambu–Jona-Lasinio (NJL) model of QCD to perform calculations of proton GPDs.

The NJL model is an effective model of quark interactions based on a four-fermi contact interaction Vogl and Weise (1991); Klevansky (1992); Hatsuda and Kunihiro (1994). It successfully incorporates several low-energy aspects of QCD, most notably (approximate) chiral symmetry and its dynamical breaking. The breaking of chiral symmetry dresses the quarks, causing them to propagate with a large effective mass $M\sim 400$ MeV, as is described by the gap equation. Moreover, confinement can be simulated in the NJL model through use of proper time regularization Ebert et al. (1996); Hellstern et al. (1997); Cloët et al. (2014). The NJL model has been used to successfully describe many properties of both mesons Vogl and Weise (1991); Klevansky (1992); Cloët et al. (2014); Ninomiya et al. (2017); Freese and Cloët (2019) and baryons Ishii et al. (1993a, b, 1995); Cloët et al. (2014).

A particular aspect of proton structure we will emphasize is the presence of diquark correlations. Quark-diquark correlations have had considerable success in modeling the properties of baryons Cahill et al. (1989); Anselmino et al. (1993); Roberts et al. (2011); Cloët et al. (2014), and the presence of diquark correlations is also borne out by such evidence as the $Q^{2}$ dependence of flavor-separated form factors Cates et al. (2011) and an approximate meson-baryon supersymmetry Brodsky et al. (2016). We will thus calculate the proton’s GPDs in a dynamical quark-diquark model.

II Proton GPDs in a dynamical quark-diquark model

Generalized parton distributions (GPDs) are defined through the matrix elements of bilocal light cone correlators. The leading-twist, helicity-independent quark GPDs of a hadron are defined through:

[TABLE]

where $[x,y]$ is a Wilson line from $y$ to $x$, $p$ and $p^{\prime}$ are the initial and final momenta, and $\lambda$ and $\lambda^{\prime}$ are the initial and final helicities (if applicable). The GPDs themselves are obtained by decomposing $V_{\lambda\lambda^{\prime}}^{q}$ in terms of linearly independent Lorentz structures. For a spin-half hadron such as the proton, we have:

[TABLE]

where $P=\frac{p+p^{\prime}}{2}$, $\Delta=p^{\prime}-p$, $\xi=-2(\Delta n)/(Pn)$, $t=\Delta^{2}$, and $n$ is a lightlike vector defining the light front. The GPDs $H(x,\xi,t)$ and $E(x,\xi,t)$ are Lorentz-invariant functions of the Lorentz-invariant arguments $x$, $\xi$, and $t$. Similar correlators are defined for helicity-dependent and helicity-flip GPDs, as well as for gluon GPDs.

The proton can be considered as a bound state of three constituent quarks, which are dressed amalgamations of more elementary current quarks (and, in QCD, gluons). The bound state amplitude can be found by solving the Faddeev equation. It has been found in many model calculations Cahill et al. (1989) that two of the three quarks are often bound in a diquark correlation that is either isoscalar and Lorentz scalar, or isovector and Lorentz axial vector. We will consider these configurations specifically in the calculations to follow.

A set of “body GPDs” (named in analogy to the “body form factors” in Ref. Cloët et al. (2014)) can be found by calculating the distribution of these dressed quarks within the proton. In momentum space, an operator of the form $\Gamma\delta(n[xP-k])$ is placed on each of the dressed quark lines, with $\Gamma=\not{n}$ (times an isospin structure) for the helicity-independent body GPDs. Within the quark-diquark approximation, the quark can be within or outside of the diquark correlation, or possibly within the interaction kernel that binds the proton. Within the approximations considered in this work, the latter does not contribute to the GPDs. The diagrams that do contribute are depicted in Fig. 1.

Since GPDs describe the structure of hadrons in terms of current quarks, the body GPDs are not by themselves sufficient. Dressed quarks are not current quarks, but contain current quarks as a more elementary substructure. One can address this by dressing the bilocal operator defining the light cone correlator, or, equivalently, one can combine the body GPDs of the proton with the GPDs of the dressed constituent quarks using a convolution formula. Such a convolution equation would also have applicability to describing the non-elementary vertex in the diquark diagram of Fig. 1.

Since it is of central importance to this work, we will first consider how GPD convolution is to be done. We will then explore the (body) GPDs and transition GPDs of the diquarks, and subsequently the GPDs of the dressed constituent quarks.

II.1 The convolution equation

Let us consider a hadron $X$ to contain a composite constituent $Y$. The insertion of a bilocal operator onto $Y$ can be expanded in terms of the substructure of $Y$, as depicted in Fig. 2. Since $Y$ is in general off its mass shell, the composite operator (rightmost diagram in Fig. 2) is a function of the initial and final virtuality of $Y$. In specific cases where there is no functional dependence on virtuality (or when the dependence on virtuality can be safely neglected), a convolution formula holds for the GPDs $H_{X,i}(x,\xi,t)$ of $X$:

[TABLE]

The indices $i$ and $j$ label the multiplicity of GPDs appearing in front of the available Lorentz structures for $X$ and $Y$, respectively. $H_{Y,j}(x,\xi,t)$ are the GPDs of an on-shell $Y$. The functions $h_{Y/X,ij}(y,\xi,t)$ are the body GPDs, which are defined by using a Lorentz structure $\Gamma_{j}$ associated with the GPD $H_{Y,j}(x,\xi,t)$ in place of $\not{n}$ in the elementary bilocal operator.

When the constituents of $Y$ (which we call $Z$ for concreteness) also have non-elementary substructure, one can simply apply Eq. (3) consecutively. It is possible to show that GPD convolution is associative, that is, ${\left(h_{Y/X}\otimes h_{Z/Y}\right)\otimes H_{Z}=h_{Y/X}\otimes\left(h_{Z/Y}\otimes H_{Z}\right)}$, so the consecutive applications of Eq. (3) can be done in either order.

In Sec. II.3, we will observe that the dressed quark GPD has no functional dependence on virtuality, allowing use of Eq. (3) without caveats. The GPD of an off-shell diquark does depend on virtuality in general, but the use of a pole approximation for the diquark propagator requires replacing vertices sandwiched between the propagators by their on-shell form for consistency.

II.2 Diquark body GPDs

In order to determine the diquark diagram contribution to the proton body GPDs, we must determine the body GPDs of the diquarks themselves. We have both scalar and axial vector diquarks to consider, in addition to the transition GPD between the two diquark species.

In this work, we use a pole approximation for the diquark propagators, i.e., we take

[TABLE]

which is the dominant contribution to the propagator. Such approximations are ubiquitous in the baryon modeling literature Mineo et al. (1999); Cloet et al. (2005a, b); Mineo et al. (2005); Cloet et al. (2006, 2008a); Eichmann et al. (2008); Cloet et al. (2008b); Eichmann et al. (2009); Nicmorus et al. (2009, 2010); Matevosyan et al. (2012); Roberts et al. (2011); Wilson et al. (2012); Cloet and Miller (2012); Cloet et al. (2013); Segovia et al. (2014a, b, 2015); Carrillo-Serrano et al. (2016). The operators $\Lambda(x,\xi,t,v,v^{\prime})$ defining the diquark GPDs in general depend on the initial and final virtualities:

[TABLE]

However, this operator is always found sandwiched between propagators in the form:

[TABLE]

The pole approximation is in effect a truncated Laurent series expansion, in $v$ (or $v^{\prime}$) which is truncated above $v^{-1}$ ($v^{\prime-1}$). When multiplying truncated series expansions, the product must be truncated at the order of the approximation, meaning that any order-$v$ (order-$v^{\prime}$) terms in $\Lambda(x,\xi,t,v,v^{\prime})$ must be discarded. In other words, consistency with the pole approximation requires neglecting the functional dependence of the GPD vertex on virtuality. Thus, we replace the expression in Eq. (6) with:

[TABLE]

and proceed to consider on-shell diquark GPDs.

II.2.1 Scalar diquarks

An on-shell scalar diquark has only a single GPD111 In principle, an off-shell scalar diquark has an additional $\sf{T}$-odd GPD, but it is suppressed by the virtuality and, within the pole approximation, should be neglected.

which is identical for up and down quarks since the scalar diquark is isoscalar. If one evaluates the far-right triangle diagram in Fig. 2 with a scalar Bethe-Salpeter vertex, the matrix element decomposition is simply:

[TABLE]

The corresponding body GPD for the distribution of scalar diquarks within the proton is then found by evaluating the right diagram in Fig. 1 with $(kn)\delta(n[xP-k])$ as the vertex. This body GPD can be combined with $H_{s}^{q}(x,\xi,t)$ via Eq. (3) to obtain the scalar diquark diagram contribution to the proton body GPDs.

The scalar diquark is isoscalar, so makes equal contributions to the up and down body GPDs. Additionally, since it contains the proton’s valence down quark, diagrams with a spectator scalar diquark only contribute to the up quark body GPD (which can still contain a down current quark).

II.2.2 Axial vector diquarks

The axial vector diquark has five on-shell GPDs. The relevant Lorentz decomposition is Berger et al. (2001):

[TABLE]

where the functional dependence on $x$, $\xi$, and $t$ has been suppressed for compactness.

One can obtain an off-shell version of this by not including the polarization vectors in the calculation of $V_{a,\lambda^{\prime}\lambda}^{q}$. There appears to be an ambiguity in this, since prior to stripping the polarization vectors, the identities $(\epsilon p)=(\epsilon^{\prime*}p^{\prime})=0$ can be used to rewrite (for instance) $(\epsilon\Delta)$ in terms of $(\epsilon P)$. However, within the pole approximation, the axial vector diquark propagators are transverse, thus enforcing the similarity relations $p^{\mu}\sim 0$ and $p^{\prime\nu}\sim 0$. We may use these similarity relations to rewrite the uncontracted correlator $V_{a,\lambda^{\prime}\lambda}^{q,\mu\nu}$ with $\Delta$ as the only uncontracted momentum, giving us:

[TABLE]

The Lorentz structures above can be used to calculate the body GPDs for the distribution of axial vector diquarks within the proton, provided the substitution $P\mapsto k$ is made, and the structures are then multiplied by $(kn)\delta(n[xP-k])$.

The axial vector diquark is isovector, meaning it comes in $uu$, $ud$, and $dd$ varieties. The isospin algebra necessary to determine the weights with which diagrams involving axial vector diquarks enter into up and down quark body GPDs has previously been done, with the results in Eqs. (102,103) of Ref. Cloët et al. (2014).222 Although Ref. Cloët et al. (2014) is about form factors, the isospin algebra is the same, and the relevant isospin factors are also the same.

II.2.3 Diquark transition GPDs

Lastly, there are scalar-to-axial and an axial-to-scalar transition GPDs. For scalar $\rightarrow$ axial vector transitions, the correlator takes the form:

[TABLE]

Hermiticity and time reversal properties can be used to show that $H_{sa}^{q}(x,\xi,t)=H_{sa}^{q}(x,-\xi,t)$ and that the axial vector $\rightarrow$ scalar transition GPD satisfies $H_{as}^{q}(x,\xi,t)=-H_{sa}^{q}(x,\xi,t)$. Neglecting virtuality dependence (as required by the pole approximation), there remains one GPD in the off-shell case, since stripping the polarization vector $\epsilon$ from $V_{sa}^{q}$ can still only produce a single unique Lorentz structure.

As for the axial vector diquark, the isospin weights for the transition diagrams can be found in Eqs. (102,103) of Ref. Cloët et al. (2014).

II.3 Dressed quark GPDs

In order to calculate any hadronic matrix element within the NJL model, one must dress the operator in question. This dressing is a result of DCSB, and is just as necessary for bilocal light cone correlators and GPDs as it is for the electromagnetic current and form factors. Failing to dress the operator defining GPDs will result in its Mellin moments reproducing the incorrect electromagnetic and gravitational form factors Freese and Cloët (2019); Shi et al. . The need for dressing arises because the GPDs appear in a bilocal correlator of current quark fields, while a hadron in the NJL model is constructed from composite dressed quarks. Dressing of the operator essentially amounts to describing the dressed quark in terms of an elementary current quark substructure.

The bilocal operator defining the leading-twist GPDs is dressed according to a Bethe-Salpeter equation (BSE), which is depicted in the Hartree-Fock approximation333 This approximation is standard for NJL model calculations Vogl and Weise (1991); Klevansky (1992); Hatsuda and Kunihiro (1994), and excludes diagrams with more than one loop. by Fig. 3. One can see this by considering Mellin moments of the bilocal operator, which define local operators that each obey a BSE that is also depicted by Fig. 3. At leading twist, these local operators are traceless and accordingly do not receive contributions from operators inserted on the four-point vertex. Due to the uniqueness of the inverse Mellin transform, the bilocal operator is likewise dressed according to Fig. 3.

This BSE holds for both up and down quarks, and the NJL interaction kernel mixes these equations. It is most straightforward to solve the decoupled BSEs for the isoscalar and isovector GPDs $H_{I}(x,\xi,t)$, defined as:

[TABLE]

and analogously for $E_{I}(x,\xi,t)$, which appear in Lorentz decompositions of the operators $\not{n}\delta(n[xP-k])$ and $\not{n}\tau_{3}\delta(n[xP-k])$, respectively.

We find the following solutions for the dressed quark GPDs:

[TABLE]

and $E_{1}(z)=\int_{1}^{\infty}\mathrm{d}t\,t^{-1}e^{-zt}$ is the exponential integral function. These solutions are exact, and do not contain any explicit functional dependence on the quark virtuality, despite the quark being off-shell in general. This is a consequence of the interaction that produces the dressing being a contact interaction.

It is worth remarking on the limit $\xi\rightarrow 0$ in Eq. (13). It can be shown that the integral of $H^{\prime}_{I}(x,\xi,t)$ over $x$ is independent of $\xi$ when $\xi>0$, and that all higher Mellin moments contain an overall factor $\xi$, and thus vanish when $\xi\rightarrow 0$. From the uniqueness of the Mellin transform, we infer that when $\xi\rightarrow 0$, $H_{I}^{\prime}(x,\xi,t)$ is proportional to a Dirac delta distribution. This was also found in Ref. Shi et al. , where it was remarked that even in the zero-skewness limit, the GPD contains a “hidden ERBL region” at $x=0$. This hidden ERBL region persists through GPD convolution, meaning that numerical results at $\xi=0$ for hadron GPDs in the NJL model will necessarily be missing the delta distribution in the hidden ERBL region.

II.4 Support region

The support region in $x$ of the GPDs and body GPDs is determined during the course of evaluating the relevant diagrams. Each diagram’s contribution is non-zero only when the delta function $\delta(n[xP-k])$ is picked up by the integration over $k$. We find in particular that:

[TABLE]

except in the case of the GPDs of dressed quarks, for which $|x|<|\xi|$, as is explicitly noted by the presence of step functions in Eqs. (13). The condition $|\xi|\leq 1$ holds for any on-shell particle by virtue of kinematics, which would strengthen Eq. (14) to $x\in[-|\xi|,1]$. However, it is possible to have $|\xi|>1$ for off-shell particles. To see this, consider that:

[TABLE]

For an on-shell particle, $(np)\propto E_{p}+p_{z}$ is strictly non-negative, so the constraint $|\xi|\leq 1$ follows from the triangle inequality. On the other hand, off-shell particles are not required to satisfy any such constraint, and in fact $(np)$ can be negative. Thus $\xi$ is not constrained in general.

Since we consider off-shell particles in using the convolution formula Eq. (3), we leave the support region in Eq. (14) as general as possible. In particular, since $y\sim 0$ is present in the integral, the off-shell constituent $Y$ can have arbitrarily large skewness. We also find in numerical calculations using Eq. (3) that having support at $x>1$ for $H_{Y}(x,\xi,t)$ is necessary for $H_{X}(x,\xi,t)$ to satisfy polynomiality (see Sec. III.1).

III Proton GPD results

Several variations of the NJL model exist. In Ref. Cloët et al. (2014), the electromagnetic properties of the proton were found to be well-described within a two-flavor variant of the model. We thus use the model variant described in Ref. Cloët et al. (2014), including the numerical values for the model parameters and the approximations described therein, to calculate the helicity-independent, leading-twist proton GPDs.

We begin by presenting results for the GPDs $H^{q}(x,\xi,t)$ and $E^{q}(x,\xi,t)$ at two skewness values, and at a model renormalization scale of $\mu^{2}=M^{2}=0.16$ GeV2. In Fig. 4, we have $\xi=0$, while in Fig. 5, we consider a moderate $\xi=0.5$.

To help understand the results, we first give the relationships that the GPDs have to more familiar observables. For instance, we have the forward limit relation $H_{q}(x,0,0)=q(x)$. Additionally, Mellin moments of the GPDs are related to the electromagnetic form factors (EMFFs):

[TABLE]

where $F_{i}(t)=\sum_{q}e_{q}F_{i}^{q}(t)$, and to the gravitational form factors (GFFs):

[TABLE]

An angular momentum form factor $J^{q}(t)=\frac{1}{2}\left(A^{q}(t)+B^{q}(t)\right)$ can also be defined, and is related to Ji’s sum rule Ji (1997a).

The $x$ and $t$ dependence of the GPDs in both Figs. 4,5 can be seen to differ between the up and down quarks. This is due to the presence of multiple isospin-dependent effects. Among these is the presence of diquark correlations, whose effects can be most easily seen in the non-skewed GPDs.

At zero skewness ($\xi=0$), there is a peak in $H^{q}(x,0,t)$ for each fixed-$t$ slice. (See top panel of Fig. 4.) The location for this peak can be quantified by the average momentum fraction $\langle\langle x_{q}(t)\rangle\rangle=A^{q}(t)/F_{1}^{q}(t)$. In the forward limit, $\langle\langle x_{u}(0)\rangle\rangle=0.34\approx\langle\langle x_{d}(0)\rangle\rangle=0.32$. Each down quark thus carries about the same momentum on average than each up quark. Were only scalar diquarks present, we would expect $\langle\langle x_{u}(0)\rangle\rangle\gg\langle\langle x_{d}(0)\rangle\rangle$, since a (dressed) down quark would only be found within the diquark, thus giving the down quark a lower effective mass. However, axial vector diquark configurations are also present, and (due to how the recoupling coefficients work out) the down quark is more often alone than in the axial vector diquark. Thus, the difference between $\langle\langle x_{u}(0)\rangle\rangle$ and $\langle\langle x_{d}(0)\rangle\rangle$ is softened by the presence of axial vector diquarks.

Finite $t$ is a novel aspect of GPDs that elaborates the roles of different diquark species further. High $-t$ acts as a filter that selects for configurations where the probed parton was already carrying most of the hadron’s momentum. One can accordingly observe in Fig. 4 that increasing $-t$ moves the peaks of both GPDs to higher $x$. The down quark peak migrates further than the up quark, with $\langle\langle x_{u}(-2~{}\mathrm{GeV}^{2})\rangle\rangle=0.45$ and $\langle\langle x_{d}(-2~{}\mathrm{GeV}^{2})\rangle\rangle=0.50$. This occurs because axial vector diquark configurations fall more slowly with $-t$, so at large $-t$ the down quark becomes sampled more often outside a diquark.

Despite $\langle\langle x_{d}(t)\rangle\rangle$ exceeding $\langle\langle x_{u}(t)\rangle\rangle$ at large $-t$, the down quark GPD still falls faster than the up quark GPD with increasing $-t$, as has previously been seen in measurements of the flavor-separated electromagnetic form factors Cates et al. (2011). We define the ratio $S_{q}(x,t)=H_{q}(x,0,t)/H_{q}(x,0,0)$, which characterizes how quickly an $x$ slice of a quark GPD falls with $-t$. The super-ratio $S_{d}(x,t)/S_{u}(x,t)$ then characterizes how much faster the down quark GPD falls with $-t$ than the up quark GPD.

Such a super-ratio is plotted in Fig. 6 for several values of $x$. We find at all $x$ values that the down quark GPD falls more steeply than the up quark GPD, but that the steepness itself varies as a function of $x$. In particular, the relative steepness of the down quark falloff is less extreme at large $x$ than at small $x$. Since the steeper down quark falloff is a consequence of diquark correlations, this suggests that diquark correlations dominate at small $x$, and one is more likely to observe a direct quark outside the diquark correlation when $x$ is close to 1.

The other GPD, $E^{q}(x,\xi,t)$, does not correspond to a familiar observable in the forward limit, but can be understood by its relationships to the form factors $F_{2}^{q}(t)$, $B^{q}(t)$, and $J^{q}(t)$. Unlike with $H^{q}(x,\xi,t)$, the peaks for $E^{q}(x,\xi,t)$ are at lower $x$ for the up quarks than down quarks for all values of $t$. (This can be clearly seen in the lower panel of Fig. 4.) This occurs due to axial vector diquarks having a larger impact on $E^{q}(x,\xi,t)$ than scalar diquarks.

Another novel aspect of GPDs is finite skewness. We emphasize that a fully covariant calculation is necessary for a correct description at finite skewness, particularly in the ERBL region, where $|x|<|\xi|$. A popular non-covariant method to calculate GPDs is the overlap representation using a truncated light front basis expansion Diehl et al. (2001). The truncation in particular is not invariant under the full Lorentz group de Melo et al. (1998); Li et al. (2018), and moreover leads to missing terms in the overlap representation of the ERBL region, since the GPD in the ERBL region is generated by the overlap of Fock states with different numbers of partons.

Since the calculations done in this work are fully covariant, we are able to obtain self-consistent results at finite skewness, including in the ERBL region. In Fig. 5, GPD results are shown for $\xi=0.5$. The most immediately striking feature are the jump discontinuities at $x=\pm\xi$. These discontinuities are inherited from the dressed quark GPDs given in Eq. (13), and are not present if the quark GPDs are not dressed. These jump discontinuities are characteristic of effective model calculations with constant dressed quark mass, and have previously been observed in other model GPD calculations Petrov et al. (1998); Polyakov and Weiss (1999); Theussl et al. (2004). These discontinuities appear to present a problem for factorization of the DVCS amplitude, but factorization is not expected to apply at the model scale of $\mu^{2}=0.16$ GeV2 and evolution of the model GPDs to a scale where QCD factorization is expected to apply removes the discontinuities (as we will see in Fig. 7).

Another striking feature of Fig. 5 is the drastic difference in the $x$ dependence of the up and down GPDs in the ERBL region, where $|x|<|\xi|$. In this region, the dressed quark GPD is no longer trivial, and accordingly it is possible to find (for instance) a down current quark inside a dressed up quark.

At large $-t$, the down quark GPD $H_{d}(x,\xi,t)$ is dominated by $D_{Q}(x,\xi,t)$ in the ERBL region. This happens because the down quark body GPD falls faster, due to the dominance of scalar diquark configurations, causing $H_{d}(x,\xi,t)$ to be dominated by the current down quarks found within dressed up quarks. Moreover, the GPD for a current quark to appear within a dressed quark of the opposite flavor goes as $\frac{1}{2}(H_{I=0}(x,\xi,t)-H_{I=1}(x,\xi,t)+D_{Q}(x,\xi,t))$, where $H_{I=0}(x,\xi,t)\approx H_{I=1}(x,\xi,t)$ because of the nearly degenerate masses of the $\rho$ and $\omega$ mesons. Thus, $D_{Q}(x,\xi,t)$, which is an odd function of $x$ and is negative for $x>0$, dominates the ERBL region of $H_{d}(x,\xi,t)$ at large $-t$.

Many of the peculiar model features are washed out by GPD evolution. In Fig. 7, we present results of evolving the model scale GPDs to $\mu^{2}=4$ GeV2. Leading-order evolution kernels Ji (1997b) were used along with a zero-mass variable flavor number scheme. The jump discontinuities at $x=\pm\xi$ are removed, and the shape of $D_{Q}(x,\xi,t)$ is no longer clearly visible in the ERBL region. Since the jump discontinuities disappear at $Q^{2}$ above the model scale, the model scale discontinuities do not present a problem for factorization at scales where QCD factorization is expected to be relevant.

III.1 Polynomiality and form factors

An especially remarkable property of GPDs is polynomiality, which ensures that the $s$th Mellin moment of a GPD is a polynomial in $\xi$ of order $s$ or less444 This is true for the gluon GPD if the Ji convention Ji (1998) is used. If the Diehl convention Diehl (2003) is used, this is true instead of the $(s-1)$th Mellin moment. . For the helicity-independent GPDs of spin-half particles, the polynomials are even in $\xi$555 GPDs that are odd in $\xi$ exist, even at leading twist. $H_{4}(x,\xi,t)$ and $\tilde{H}_{3}(x,\xi,t)$ for spin-one targets are such examples Diehl (2003). The matrix element of the bilocal correlator must be time reversal invariant, and if it’s possible to construct T-odd Lorentz structures in its decomposition, the Lorentz-invariant function multiplying it must likewise be T-odd (i.e., odd in $\xi$) for the product of both to be T-even. . Polynomiality is a consequence of Lorentz covariance, with the reality of the coefficients following from hermiticity of the bilocal operator and the evenness of the polynomial from time reversal symmetry of its matrix element. Since we have observed complete Poincaré covariance throughout the calculation, our model GPDs exhibit polynomiality.

Of particular interest are the cases $s=1$, which reproduce partonic contributions (without charge weights) to electromagnetic form factors (EMFFs), and $s=2$, which give gravitational form factors (GFFs). These relationships are given in Eqs. (16,17).

The EMFFs of the proton have been previously calculated within the NJL model (see, e.g., Ref. Cloët et al. (2014)), but the GFFs have not been. In Fig. 8, we present numerical results for the EMFFs and GFFs extracted from the GPDs calculated in this work. Since pion loops have not been included in this calculation, our EMFFs should be compared to the “bse” results from Ref. Cloët et al. (2014).

There are several constraints that the form factors must obey due to symmetries and conservation laws. Charge and momentum conservation require $F_{1}(0)=1$ and $A(0)=1$, respectively, and our results satisfy these relations. Conservation of angular momentum gives the Ji sum rule Ji (1997a) \frac{1}{2}\big{(}A(0)+B(0)\big{)}=\frac{1}{2}, or equivalently (when combined with momentum conservation) $B(0)=0$—a statement otherwise known as the vanishing of the anomalous gravitomagnetic moment. We find that indeed $B(0)=0$ within our results, and remark that this is an inevitable consequence of observing Lorentz covariance throughout the calculation.

The remaining static observables $F_{2}(0)$ and $C(0)$ are not constrained by symmetries or conservation laws. The first of these gives the anomalous magnetic moment of the proton. Empirically, $F_{2}(0)=\kappa_{p}=1.793$, but we underestimate this, finding $F_{2}(0)=1.39$. It was observed in Ref. Cloët et al. (2014) that perturbatively introducing a pion cloud can significantly close the gap between the calculated and model values. The second of these, $C(0)$, is commonly known as the “D-term” Polyakov and Weiss (1999), and is neither constrained by symmetries nor by experiment.666 There does however exist a recent model-dependent phenomenological extraction of the quark contribution from JLab Hall B data Burkert et al. (2018); Kumerički (2019).

We find $C(0)=-1.08$.

The form factor $C(t)$ has been interpreted as describing the distribution of forces within hadrons Polyakov and Schweitzer (2018), and the fact that $C(0)<0$ is understood as an important stability criterion. Since $C(0)$ is not constrained by any symmetries, it is possible for operator dressing to change its value—in contrast to $A(0)$ or $F_{1}(0)$. In Ref. Freese and Cloët (2019), $C_{\pi}(t)$ was found to change by a factor of $\sim 3$ from dressing the quark-graviton vertex, and introducing the dressing was necessary to satisfy a low-energy pion theorem. In a similar vain, we find that dressing the nonlocal correlator via the BSE depicted in Fig. 3 is necessary for proton stability. If the bare nonlocal operator is used to calculate the proton GPDs, then we obtain $C(0)=0.85>0$, at stark odds with the apparent stability of the proton.

To understand the necessity of operator dressing for observing proton stability, we reiterate that dressing can be seen as accounting for the elementary current quark substructure of the three dressed quarks that comprise the proton. In the absence of dressing, we would be attempting to describe the proton as being made of three current quarks, which would not be accurate, and would account for only the leading Fock state of the proton. That we get $C(0)>0$ in this case tells us that a hypothetical hadron with the mass and quantum numbers of the proton that is made of only three elementary current quarks would not be mechanically stable in this framework.

Since the NJL model contains only quarks and $C(0)$ is not constrained by any symmetries, it is possible that $C(0)$ may change significantly with the introduction of gluons. Therefore, our finding of $C(0)=-1.08$ should at best be interpreted as a prediction for the quark contribution to the D-term, rather than for the overall D-term of the proton.

IV Summary and outlook

We have calculated the helicity-independent, leading-twist GPDs of the proton at finite skewness, in a confining version of the NJL model. Dressing of non-local operator defining the light cone correlator—which happens as a result of DCSB—was necessary for sensible results to be obtained, including a negative D-term compatible with the stability of the proton. The Lorentz covariance of the model and all approximations made ensured that polynomiality and sum rules relating to charge, energy-momentum, and angular momentum conservation were all obeyed, giving the vanishing of the anomalous gravitomagnetic moment as a corollary.

It will be possible in future work to apply the same formalism to the helicity-dependent and helicity-flip GPDs of the proton. Moreover, the Lagrangian can be generalized to include immersion in a finite-density medium, allowing predictions of GPDs in nuclear matter and predictions for a generalized EMC effect.

Acknowledgments

We would like to thank Chao Shi for illuminating discussions that helped contribute to our investigation. This work was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, contract no. DE-AC02-06CH11357, and an LDRD initiative at Argonne National Laboratory under Project No. 2017-058-N0.