Phenomenological assessment of proton mechanical properties from deeply virtual Compton scattering

Abstract

A unique feature of generalised parton distributions is their relation to the QCD energy–momentum tensor. In particular, they provide access to the mechanical properties of the proton i.e. the distributions of pressure and shear stress induced by its quark and gluon structure. In principle the pressure distribution can be experimentally determined in a model-independent way from a dispersive analysis of deeply virtual Compton scattering data through the measurement of the subtraction constant. In practice the kinematic coverage and accuracy of existing experimental data make this endeavour a challenge. Elaborating on recent global fits of deeply virtual Compton scattering measurements using artificial neural networks, our analysis presents the current knowledge on this subtraction constant and assesses the impact of the most frequent systematic assumptions made in this field of research. This study will pave the way for future works when more precise data will become available, e.g. obtained in the foreseen electron-ion colliders EIC and EIcC.

1 Introduction

Quantum chromodynamics (QCD) provides today’s standard description of hadrons. In this theory, hadrons are depicted as complex objects made out of fundamental building blocks referred to as partons: quarks and gluons. This picture has been extensively studied over the last fifty years, but many of the basic questions about partonic media still await a precise quantitative answer to this day. In particular, in spite of important theoretical advances since the beginning of the 2000s, our knowledge of distributions of certain properties like pressure and shear stress inside hadrons – the so-called “mechanical” properties – is still poor. Understanding those properties therefore remains among the main challenges faced by modern nuclear and high energy physics.

QCD factorisation theorems provide us with tools to study the properties of partonic media up to an arbitrarily high accuracy, only limited by the order of the perturbative expansion describing the short distance part of scattering processes. In particular, generalised parton distributions (GPDs) offer a rigorous theoretical framework that can be used to study the 3D structure of hadrons [1,2,3,4,5]. GPDs offer a consistent unified framework encompassing the usual parton distribution functions (PDFs) and elastic form factors (EFFs), and providing even more information. Such a non-trivial connection is used to obtain tomographic pictures of the nucleon, where spatial distributions of partons carrying a fraction of the nucleon momentum are projected onto the plane perpendicular to the direction of the nucleon motion [6,7,8]. Another exciting feature of GPDs is their relation to the QCD energy–momentum tensor (EMT), which would otherwise probably only be accessible via graviton scattering. This unique relation allows to evaluate the total angular momentum carried by gluons or quarks of a given flavour [2, 3], which is essential to solve the long-standing puzzle of the nucleon spin decomposition that emerged 30 years ago with the EMC measurements [9]. On top of this, the relation between GPDs and the EMT can also be used to access information about the mechanical properties of partonic systems [10, 11], like distributions of pressure inside the nucleon – which is the subject of this article.

The potential of studying the mechanical properties of the proton in the formalism of GPDs was first recognised in Ref. [12]. This subject has recently gained a lot of interest and seen vital theoretical progress. In addition, precise data have become available, recently allowing for the first data-driven estimates of a nucleon pressure distribution [13, 14]. These analyses aim at extracting the so-called subtraction constant from deeply virtual Compton scattering (DVCS) data. The constant naturally appears in the dispersive studies of DVCS amplitudes and it carries information about one of the EMT form factors that is crucial for the understanding of the mechanical properties. We point out that the subject in general has been recognised as a key element of current or future experiments, in particular those to be conducted in the electron-ion collider (EIC) [15, 16], Chinese electron-ion collider (EIcC) [17, 18], large hadron-electron collider (LHeC) [19].

In this article we benchmark the extraction of the subtraction constant from DVCS data. We present results of our extraction of the subtraction constant from the existing world DVCS proton data. To deliver a self-contained exposition, we remind the leading-order (LO) QCD evolution formulas, which in general are known, but appear in the literature in a scattered form. Compared to the first data-driven estimates, our analysis is characterised by a reduced model dependence, which is achieved by the use of artificial neural network techniques to describe the contribution of the four chiral-even leading-twist GPDs. The results of this analysis summarise our current knowledge of the DVCS subtraction constant and assess the impact of the most common systematic assumptions made in that research topic. This will pave the way for future studies when new precise data will be released.

This article is organised as follows. A brief introduction to the theoretical frameworks of EMT and GPDs is given in Sects. 2 and 3 respectively. Section 4 is dedicated to the LO evolution of the subtraction constant, and Sect. 5 to its modelling. The extraction based on the artificial neural network technique is presented in Sect. 6 and the corresponding results are discussed in Sect. 7. The summary is given in Sect. 8.

Notations Throughout the text we denote spatial vectors by boldface symbols. $\eta _{\mu \nu } = \text {diag}(+,-,-,-)$ is the Minkowski metric. For convenience, we use the compact notations: $a^{\{\mu }b^{\nu \}} = a^\mu b^\nu + a^\nu b^\mu$ and $a^{[\mu }b^{\nu ]} = a^\mu b^\nu - a^\nu b^\mu$.

2 Energy–momentum tensor

In the most general case, the proton matrix elements of the local gauge-invariant EMT operator can be parameterised in terms of five gravitational form factors (GFFs): $A_{a}(t)$, $B_{a}(t)$, $C_{a}(t)$, ${{\bar{C}}}_{a}(t)$ and $D_{a}(t)$; as follows [20, 21]:

$$\begin{aligned} \langle p',\varvec{s'}\vert T^{\mu \nu }_a(0)\vert p,\varvec{s}\rangle&= {\bar{u}}(p',\varvec{s'}) \Bigg \{ \frac{P^\mu P^\nu }{M}\,A_a(t) \nonumber \\&\quad + \frac{\varDelta ^\mu \varDelta ^\nu - \eta ^{\mu \nu }\varDelta ^2}{M}\, C_a(t) + M \eta ^{\mu \nu }{\bar{C}}_a(t)\nonumber \\&\quad + \frac{P^{\{\mu } i\sigma ^{\nu \}\rho }\varDelta _\rho }{4M}\left[ A_a(t)+B_a(t)\right] \nonumber \\&\quad + \frac{P^{[\mu } i\sigma ^{\nu ]\rho }\varDelta _\rho }{4M}\,D_a(t)\Bigg \} u(p,\varvec{s}). \end{aligned}$$

(1)

Here, $t=\varDelta ^2$ with $\varDelta =p'-p$ the four-momentum transfer to the proton, $P=(p'+p)/2$ is the average four-momentum, M is the proton mass, $u(p,\varvec{s})$, ${\bar{u}}(p',\varvec{s'})$ are Dirac spinors with the covariant normalisation ${\bar{u}}(p,\varvec{s}) u(p,\varvec{s}) = 2M$, and $\varvec{s}$, $\varvec{s'}$ are the rest-frame polarisation vectors. The label a denotes either the quark flavour ($a=q$) or the gluon ($a=g$) contribution to the EMT. We refer to Ref. [22] for a detailed discussion of the specific expressions of the EMT, and their decompositions into quark and gluon contributions.

In the following, a summation over parton types a is meant each time a generic GFF $G \in \{A, B, C, {\bar{C}}, D\}$ is written without any parton type index:

$$\begin{aligned} G(t) = \sum _{a=q,g} G_a(t) , \end{aligned}$$

(2)

and similarly a summation over quark flavours q (flavour-singlet combination) is meant each time this generic GFF G is written with the superscript S:

$$\begin{aligned} G^S(t) = \sum _{a=q} G_a(t) . \end{aligned}$$

(3)

When needed, we indicate the quark content with:

$$\begin{aligned} G^{u+d}(t)= & {} G^u(t) + G^d(t) , \end{aligned}$$

(4)

$$\begin{aligned} G^{u+d+s}(t)= & {} G^u(t) + G^d(t) + G^s(t) . \end{aligned}$$

(5)

Poincaré symmetry implies the following constraints on the GFFs [2, 23,24,25]:

$$\begin{aligned} A(0)&= 1 , \end{aligned}$$

(6)

$$\begin{aligned} B(0)&= 0 , \end{aligned}$$

(7)

$$\begin{aligned} {\bar{C}}(t)&= 0 , \end{aligned}$$

(8)

and [26]:

$$\begin{aligned} D^S(t)&= -G_A(t) , \end{aligned}$$

(9)

$$\begin{aligned} D_g(t)&= 0 , \end{aligned}$$

(10)

where $G_A(t)$ is the flavour-singlet axial-vector form factor. The GFFs $A_a(t)$, $B_a(t)$ and $C_a(t)$ can be related to leading-twist GPDs and thus can be probed at present experimental facilities, see Sect. 3. Accessing ${\bar{C}}_a(t)$ is more involved since it is related to higher-twist distributions [21, 27, 28].

The static EMT is defined as the following Fourier transform in the Breit frame, where ${\varvec{P}}={\varvec{0}}$ and $t=-\varvec{\varDelta }^2$ [11, 12, 29]:

$$\begin{aligned} {\mathcal {T}}^{\mu \nu }_a(\varvec{r})=\int \frac{\mathrm{d}^3\varvec{\varDelta }}{(2\pi )^3}\,e^{-i\varvec{\varDelta }\cdot \varvec{r}}\, \frac{\langle p',\varvec{s}\vert T^{\mu \nu }_a(0)\vert p,\varvec{s}\rangle }{2P^0} . \end{aligned}$$

(11)

It indicates how energy and momentum are distributed inside the proton with the canonical polarisation $\varvec{s}$. Writing $r = |{\varvec{r}}|$ the radial coordinate, the energy $\varepsilon _a(r)$, radial pressure $p_{r,a}(r)$ and tangential pressure $p_{t,a}(r)$ distributions can be expressed in the Breit frame as Fourier transforms of the GFFs $A_a(t)$, $B_a(t)$, $C_a(t)$ and ${\bar{C}}_a(t)$:

$$\begin{aligned} \varepsilon _a(r)&=M\int \frac{\mathrm{d}^3\varvec{\varDelta }}{(2\pi )^3}\,e^{-i\varvec{\varDelta }\cdot \varvec{r}} \nonumber \\&\quad \times \left\{ A_a(t)+\bar{C}_a(t) + \frac{t}{4M^2} [B_a(t)- 4C_a(t) ]\right\} , \end{aligned}$$

(12)

$$\begin{aligned} p_{r,a}(r)&= M\int \frac{\mathrm{d}^3\varvec{\varDelta }}{(2\pi )^3}\,e^{-i\varvec{\varDelta }\cdot \varvec{r}} \nonumber \\&\quad \times \left\{ -\bar{C}_a(t)-\frac{4}{r^2}\frac{t^{-1/2}}{M^2}\frac{\mathrm{d}}{\mathrm{d}t}\! [t^{3/2}\,C_a(t) ]\right\} , \end{aligned}$$

(13)

$$\begin{aligned} p_{t,a}(r)&= M\int \frac{\mathrm{d}^3\varvec{\varDelta }}{(2\pi )^3}\,e^{-i\varvec{\varDelta }\cdot \varvec{r}} \nonumber \\&\quad \times \left\{ -\bar{C}_a(t)+\frac{4}{r^2}\frac{t^{-1/2}}{M^2}\frac{\mathrm{d}}{\mathrm{d}t}\!\left( t\frac{\mathrm{d}}{\mathrm{d}t}\!\left[ t^{3/2}\,C_a(t)\right] \right) \right\} . \end{aligned}$$

(14)

In a relativistic system like the proton, pressure forces are usually not isotropic. The isotropic pressure $p_a(r)$ and pressure anisotropy $s_a(r)$ are then defined in terms of the radial and tangential pressures:

$$\begin{aligned} p_a(r)&= [p_{r,a}(r)+2p_{t,a}(r) ]/3 , \end{aligned}$$

(15)

$$\begin{aligned} s_a(r)&= p_{r,a}(r)-p_{t,a}(r) , \end{aligned}$$

(16)

or, in terms of GFFs:

$$\begin{aligned} p_a(r)&= M\int \frac{\mathrm{d}^3\varvec{\varDelta }}{(2\pi )^3}\,e^{-i\varvec{\varDelta }\cdot \varvec{r}} \left\{ -\bar{C}_a(t) +\frac{2}{3} \frac{t}{M^2}\,C_a(t)\right\} , \end{aligned}$$

(17)

$$\begin{aligned} s_a(r)&= -\frac{4 M}{r^2}\int \frac{\mathrm{d}^3\varvec{\varDelta }}{(2\pi )^3}\, e^{-i\varvec{\varDelta }\cdot \varvec{r}}\,\frac{t^{-1/2}}{M^2}\frac{\mathrm{d}^2}{\mathrm{d}t^2}\! [t^{5/2}\,C_a(t) ] . \end{aligned}$$

(18)

Among the five distributions: $\varepsilon _a(r)$, $p_{r,a}(r)$, $p_{t,a}(r)$, $p_a(r)$ and $s_a(r)$; the pressure anisotropy $s_a(r)$ is the only one that does not depend on the GFF $\bar{C}_a(t)$. This feature makes its current experimental access the least challenging.

These distributions permit the introduction of additional definitions of proton radius beyond the usual electromagnetic charge radii. Indeed the distributions $\varepsilon (r)=\sum _a\varepsilon _a(r)$ and $p_r(r)=\sum _ap_{r,a}(r)$ are expected to be positive and can be used to define the energy (or mass) and mechanical radii of the nucleon [11, 29] in terms of the GFFs A and C (summed over all constituents):

$$\begin{aligned} \langle r^2\rangle _E&=\frac{1}{M}\int \mathrm{d}^3{\varvec{r}}\,r^2\,\varepsilon (r)=6\left[ A'(0)-\frac{C(0)}{M^2}\right] , \end{aligned}$$

(19)

$$\begin{aligned} \langle r^2\rangle _\text {mech}&=\frac{1}{\mathcal {P}_r}\int \mathrm{d}^3\varvec{r}\,r^2\,p_r(r)=\frac{6C(0)}{\int _{-\infty }^0\mathrm{d}t\,C(t)} , \end{aligned}$$

(20)

with ${\mathcal {P}}_r=\int \mathrm{d}^3{\varvec{r}}\,p_r(r)$ and $A'(t)=\mathrm{d}A(t)/\mathrm{d}t$.

3 Generalised parton distributions

GPDs have been introduced by the end of the 90s [1,2,3,4,5] and now constitute a mature field of research [30,31,32,33,34,35] with experimental programmes carried out at CERN (COMPASS), DESY (HERMES, H1 and ZEUS), Jefferson Lab (Halls A, B and C) and in the future at the electron-ion collider EIC. References [36,37,38,39] provide a recent account of the experimental and phenomenological status and prospects of three exclusive channels which attract most of current experimental interest, namely DVCS, deeply virtual meson production (DVMP) and timelike Compton scattering (TCS).

The formal definition of GPDs in terms of matrix elements can be found in Ref. [31] and we will use the same notations. From the point of view of the current EMT studies, we are interested in the four leading-twist chiral-even GPDs: $H^{a}(x, \xi , t)$, $E^{a}(x, \xi , t)$, ${\widetilde{H}}^{a}(x, \xi , t)$ and ${\widetilde{E}}^{a}(x, \xi , t)$ where $a=q, g$. Here, x is the average longitudinal light-front momentum fraction of the active parton and $\xi$ is the skewness variable describing the transfer of longitudinal light-front momentum to the system. GPDs also depend on both factorisation, $\mu _{\mathrm{F}}$, and renormalisation, $\mu _{\mathrm{R}}$, scales, which will be implicit (unless when specifically needed) to keep concise notations.

The connection between GPDs and GFFs $A_{a}(t)$, $B_{a}(t)$ and $C_{a}(t)$ is the following [31]:

$$\begin{aligned} \int _{-1}^{1}\mathrm{d}x\,x\,H^q(x,\xi ,t)&=A_q(t)+4\xi ^2C_q(t) , \end{aligned}$$

(21)

$$\begin{aligned} \int _{-1}^{1}\mathrm{d}x\,x\,E^q(x,\xi ,t)&=B_q(t)-4\xi ^2C_q(t) , \end{aligned}$$

(22)

$$\begin{aligned} \sum _q\int _{-1}^{1}\mathrm{d}x\,{{\tilde{H}}}^q(x,\xi ,t)&=G_A(t) , \end{aligned}$$

(23)

for quarks, and:

$$\begin{aligned} \int _{-1}^{1}\mathrm{d}x\,H^g(x,\xi ,t)= & {} A_g(t)+4\xi ^2C_g(t) , \end{aligned}$$

(24)

$$\begin{aligned} \int _{-1}^{1}\mathrm{d}x\,E^g(x,\xi ,t)= & {} B_g(t)-4\xi ^2C_g(t) , \end{aligned}$$

(25)

for gluons. The GPD ${\widetilde{E}}^{a}(x, \xi , t)$ does not relate to these GFFs,¹ but is a significant ingredient in the description of DVCS amplitudes and therefore important in the phenomenological applications.

We focus now on the extraction of the GFF $C_{a}(t)$ from DVCS data. The DVCS amplitude can be parameterised in terms of structure functions known as Compton form factors (CFFs). For illustration we only give the formulas for the CFF ${\mathcal {H}}(\xi , t)$ related to the GPD $H^{a}(x, \xi , t)$:

$$\begin{aligned} \mathcal {H}(\xi ,t,Q^2) = \sum _q \mathcal {H}^q(\xi ,t,Q^2) + \mathcal {H}^g(\xi ,t,Q^2) , \end{aligned}$$

(26)

where the quark $\mathcal {H}^q$ and gluon $\mathcal {H}^g$ contributions read:

$$\begin{aligned} \mathcal {H}^q(\xi ,t,Q^2)&= \int _{-1}^{~1} \frac{\mathrm {d}x}{\xi } \bigg \{ T^q\!\left( \frac{x}{\xi }, \mu _{\mathrm{F}}^2, Q^2\right) H^q(x, \xi , t,\mu _{\mathrm{F}}^2) \bigg \} , \end{aligned}$$

(27)

$$\begin{aligned} \mathcal {H}^g(\xi ,t,Q^2)&= \int _{-1}^{~1} \frac{\mathrm {d}x}{\xi } \bigg \{T^g\!\left( \frac{x}{\xi }, \mu _{\mathrm{F}}^2, Q^2\right) \frac{H^g(x, \xi , t,\mu _{\mathrm{F}}^2)}{x} \bigg \} . \end{aligned}$$

(28)

Here, $Q^2$ is the virtuality of the photon mediating the interaction between the lepton beam and the proton target in DVCS. The functions $T^q$ and $T^g$ are the renormalised coefficient functions (see Ref. [41] and refs. therein, in slightly different form – here we added the explicit factors of x and $\xi$ to emphasize the $\frac{x}{\xi }$ dependence of $T^{q,g}$) which for the LO description of the DVCS hard scattering kernel read:

$$\begin{aligned} T^q= & {} -e_q^2\,\xi /(x+\xi -i\epsilon ) - (x \rightarrow -x) , \end{aligned}$$

(29)

$$\begin{aligned} T^g= & {} 0 , \end{aligned}$$

(30)

where $e_q$ is the quark fractional electric charge in units of the proton charge |e|.

In Eqs. (27) and (28) we see that CFFs are convolutions of GPDs. The so-called deconvolution problem, which consists in extracting GPDs from CFFs, is far from trivial, and is still an open question today. In the absence of non-parametric extractions of GPDs, an experimental determination of most GFFs thus comes with a significant model uncertainty. However it is possible to access the GFF $C_a(t)$ at the level of amplitudes, i.e. without the deconvolution of neither $H^{a}(x, \xi , t)$ nor $E^{a}(x, \xi , t)$. This allows one in principle to avoid any model dependence in the extraction of $C_a(t)$, which makes this situation unique in the context of GPD studies.

Thanks to the analytic properties of CFFs [42,43,44], a dispersion relation connects the real and imaginary parts of the CFF $\mathcal {H}$, which reads:

(31)

where denotes the principal value integral and ${\mathcal {C}}_H(t,Q^2)$ is the so-called subtraction constant (here constant over $\xi$). Similar equations hold for the CFFs ${\mathcal {E}}$, $\mathcal {{\widetilde{H}}}$ and $\mathcal {{\widetilde{E}}}$, but with the opposite subtraction constant for ${\mathcal {E}}$ and null subtraction constants for $\mathcal {{\widetilde{H}}}$ and $\mathcal {{\widetilde{E}}}$.

The D-term is defined by the coefficients proportional to the highest powers of $\xi$ in the polynomiality relations of a given GPD [31]. It therefore allows one to access the GFF $C_a(t)$ via Eqs. (21) and (24). Remarkably the subtraction constant ${\mathcal {C}}_H$ is related (to any order in QCD perturbation theory) to the quark and gluon D-terms in the following way:

$$\begin{aligned} {\mathcal {C}}_H(t,Q^2) = \sum _q {\mathcal {C}}^q_H(t,Q^2) + {\mathcal {C}}^g_H(t,Q^2) , \end{aligned}$$

(32)

with (omitting the dependence on $\mu _{\mathrm{F}}^2$ and $Q^2$ for the sake of conciseness):

$$\begin{aligned} {\mathcal {C}}^q_H(t,Q^2)= & {} \frac{2}{\pi } \int _1^\infty \mathrm {d}\omega \, \mathrm {Im} T^q(\omega ) \int _{-1}^{~1}\mathrm {d}z \, \frac{D^q_{\mathrm {term}}(z,t)}{\omega -z} , \nonumber \\\end{aligned}$$

(33)

$$\begin{aligned} {\mathcal {C}}^g_H(t,Q^2)= & {} \frac{2}{\pi } \int _1^\infty \frac{\mathrm {d}\omega }{\omega } \, \mathrm {Im} T^g(\omega ) \int _{-1}^{~1}\mathrm {d}z \, \frac{D^g_{\mathrm {term}}(z,t)}{\omega -z}.\nonumber \\ \end{aligned}$$

(34)

At LO the explicit gluon dependence of the subtraction constant vanishes and the imaginary part of the quark coefficient function merely selects the value $\omega = 1$:

$$\begin{aligned} {\mathcal {C}}_H(t,Q^2) = 2\sum _q e_q^2\,\int _{-1}^{~1}\mathrm{d}z\,\frac{D_{\mathrm {term}}^q(z,t,\mu _{\mathrm{F}}^2 \equiv Q^2)}{1-z} . \end{aligned}$$

(35)

Here, $Q^{2}$ is identified with $\mu _{\mathrm{F}}^2$. Because Gegenbauer polynomials diagonalise the ERBL evolution equations at LO, the D-term is often expanded into a series of these polynomials [30]:

$$\begin{aligned} D_{\mathrm {term}}^q(z,t,\mu _{\mathrm{F}}^2)&= (1-z^2) \sum _{\text {odd } n}d_n^q(t,\mu _{\mathrm{F}}^2)\, C_n^{3/2}(z) , \end{aligned}$$

(36)

$$\begin{aligned} D_{\mathrm {term}}^g(z,t,\mu _{\mathrm{F}}^2)&= \frac{3}{2} (1-z^2)^2 \sum _{\text {odd } n}d_n^g(t,\mu _{\mathrm{F}}^2)\, C_{n-1}^{5/2}(z) , \end{aligned}$$

(37)

where $d_n^q(t,\mu _{\mathrm{F}}^2)$, $d_n^g(t,\mu _{\mathrm{F}}^2)$ are the coefficients of the expansion. The values of the first three flavour-singlet coefficients:

$$\begin{aligned} d_1^{u+d}(0 , \mu _0^2)\simeq & {} - 4.0, \end{aligned}$$

(38)

$$\begin{aligned} d_3^{u+d}(0, \mu _0^2)\simeq & {} - 1.2, \end{aligned}$$

(39)

$$\begin{aligned} d_5^{u+d}(0, \mu _0^2)\simeq & {} - 0.4. \end{aligned}$$

(40)

from the chiral quark-soliton model ($\chi$QSM) at a low scale $\mu _0 \simeq 600~\mathrm {MeV}$, are reported in Ref. [30]. We observe that each coefficient is roughly three times bigger than its successor.

Gegenbauer polynomials $C_n^{\alpha }(z)$ are classically defined as a family of orthogonal polynomials with respect to the weight function $(1-z^2)^{\alpha -1/2}$, but they can also be regarded as coefficients of the formal series (in the variable t) of the following generating function:

$$\begin{aligned} \frac{1}{(1 - 2 z t + t^2)^\alpha } = \sum _{n=0}^\infty C_n^{\alpha }(z) t^n . \end{aligned}$$

(41)

For $|t|<1$ we observe that:

$$\begin{aligned}&\sum _{n=0}^\infty \int _{-1}^1 \mathrm{d}z\,\frac{(1-z^2)C_n^{3/2}(z)}{1-z}\,t^n\nonumber \\&\quad =\int _{-1}^{+1} \mathrm {d}z \, \frac{1+z}{(1 - 2 z t + t^2)^{3/2}} = \frac{2}{1-t} = \sum _{n=0}^\infty 2t^n . \end{aligned}$$

(42)

Therefore the DVCS subtraction constant at LO can be conveniently written as:

$$\begin{aligned} {\mathcal {C}}_H(t,Q^2) = 4\sum _q e_q^2\,\sum _{\text {odd } n} d_n^q(t,\mu _{\mathrm{F}}^2 \equiv Q^2) . \end{aligned}$$

(43)

The first element of this expansion carries information about the GFF $C_{a}(t)$. By definition of the D-term, using (21) and (24):

$$\begin{aligned} \int _{-1}^1 \mathrm {d}z\, z D_{\mathrm {term}}^q(z,t)&= 4 C_q(t) , \end{aligned}$$

(44)

$$\begin{aligned} \int _{-1}^1 \mathrm {d}z\, D_{\mathrm {term}}^g(z,t)&= 4 C_g(t) . \end{aligned}$$

(45)

Since $C_1^{3/2}(z) = 3z$ and $C_0^{5/2}(z) = 1$, the definition of Gegenbauer polynomials $C_n^{3/2}(z)$ and $C_n^{5/2}(z)$ as orthogonal polynomials associated to the respective weights $1-z^2$ and $(1-z^2)^2$, complemented by the expansions (36) and (37), yields:

$$\begin{aligned} d^q_1(t,\mu _{\mathrm{F}}^2)= & {} 5C_q(t,\mu _{\mathrm{F}}^2) , \end{aligned}$$

(46)

$$\begin{aligned} d^g_1(t,\mu _{\mathrm{F}}^2)= & {} 5C_g(t,\mu _{\mathrm{F}}^2) . \end{aligned}$$

(47)

Note that:

$$\begin{aligned} d(t)= & {} \sum _{q}d_1^{q}(t, \mu _{\mathrm{F}}^2) + d_1^{g}(t, \mu _{\mathrm{F}}^2) \nonumber \\= & {} 5 \Bigg ( \sum _{q}C_q(t, \mu _{\mathrm{F}}^2) + C_g(t, \mu _{\mathrm{F}}^2) \Bigg ) , \end{aligned}$$

(48)

becomes scale independent as it is one of the gravitational form factors of the full EMT.

Through Eq. (18), the two relations (46) and (47) provide a direct handle on the distribution of pressure anisotropy (for quarks and gluons) in the proton:

$$\begin{aligned} s_q(r)= & {} -\frac{4 M}{5 r^2}\int \frac{\mathrm{d}^3\varvec{\varDelta }}{(2\pi )^3}\, e^{-i\varvec{\varDelta }\cdot {\varvec{r}}}\,\frac{t^{-1/2}}{M^2}\frac{\mathrm{d}^2}{\mathrm{d}t^2}\!\left[ t^{5/2}\,d^q_1(t)\right] , \end{aligned}$$

(49)

$$\begin{aligned} s_g(r)= & {} -\frac{4 M}{5 r^2}\int \frac{\mathrm{d}^3\varvec{\varDelta }}{(2\pi )^3}\, e^{-i\varvec{\varDelta }\cdot {\varvec{r}}}\,\frac{t^{-1/2}}{M^2}\frac{\mathrm{d}^2}{\mathrm{d}t^2}\!\left[ t^{5/2}\,d^g_1(t)\right] , \end{aligned}$$

(50)

where the dependence on the factorisation scale $\mu _{\mathrm{F}}^2$ is implicit. The phenomenological challenge of the experimental determination of pressure forces in the proton thus boils down to the extraction of the first term in the expansion (43) of the (measurable) DVCS subtraction constant. In the next section we will show that each term in this expansion behaves differently under LO evolution, allowing one to separate $d^q_1(t,\mu _{\mathrm{F}}^2)$ and $d^g_1(t,\mu _{\mathrm{F}}^2)$ from the rest of the series (43).

4 Evolution of D-term and DVCS subtraction constant

The evolution of the D-term (and therefore of the DVCS subtraction constant) is governed by the ERBL evolution equations [45, 46]. These equations can be explicitly solved at LO once the D-term is expressed in terms of Gegenbauer polynomials, just like in Eqs. (36) and (37). We follow the presentation of Ref. [31].

Let $n_f$ denote the number of active quark flavours. The anomalous dimensions $\gamma _{n}$, $\gamma _n^{\pm }$ driving the LO evolution read:

$$\begin{aligned} \gamma _{n}&= \gamma _{QQ}(n) , \end{aligned}$$

(51)

$$\begin{aligned} \gamma _n^{\pm }&= \frac{1}{2}\bigg ( \gamma _{QQ}(n) + \gamma _{GG}(n)\nonumber \\&\quad \pm \sqrt{[\gamma _{QQ}(n) - \gamma _{GG}(n)]^2 + 4\gamma _{QG}(n)\gamma _{GQ}(n)} \,\bigg ) , \end{aligned}$$

(52)

where [47]:

$$\begin{aligned} \gamma _{QQ}(n)&= C_{F}\left( \frac{1}{2} - \frac{1}{(n+1)(n+2)} + 2\sum _{k=2}^{n+1} \frac{1}{k}\right) , \end{aligned}$$

(53)

$$\begin{aligned} \gamma _{QG}(n)&= -n_{f}\, T_{F}\frac{n^2+3n+4}{n(n+1)(n+2)} ,\end{aligned}$$

(54)

$$\begin{aligned} \gamma _{GQ}(n)&= -2\, C_{F} \frac{n^2+3n+4}{(n+1)(n+2)(n+3)} ,\end{aligned}$$

(55)

$$\begin{aligned} \gamma _{GG}(n)&= \frac{2}{3} n_{f} T_{F} +C_{A} \nonumber \\&\quad \times \left( \frac{1}{6} - \frac{2}{n(n+1)} - \frac{2}{(n+2)(n+3)} + 2 \sum _{k=2}^{n+1} \frac{1}{k} \right) . \end{aligned}$$

(56)

Here $C_{F} = \nicefrac {4}{3}$, $T_{F} = \nicefrac {1}{2}$ and $C_{A} = 3$.

At the leading logarithmic accuracy and for a non-singlet combination of $q_{1}$ and $q_{2}$ quark flavours, the evolution from the initial scale $\mu _{\mathrm{F,0}}^2$ to the final scale $\mu _{\mathrm{F}}^2$ is given by:

$$\begin{aligned}&d_n^{q_1}(t, \mu _{\mathrm{F}}^2) - d_n^{q_2}(t, \mu _{\mathrm{F}}^2) \nonumber \\&\quad = \left[ d_n^{q_1}(t, \mu _{\mathrm{F,0}}^2) - d_n^{q_2}(t, \mu _{\mathrm{F,0}}^2) \right] \left( \frac{ \alpha _{\mathrm{s}}(\mu _{\mathrm{F}}^2)}{\alpha _{\mathrm{s}}(\mu _{\mathrm{F,0}}^2)} \right) ^{\textstyle \frac{2\gamma _n}{\beta _0}} . \end{aligned}$$

(57)

Introducing the mixing coefficients, $a_n^{\pm }$:

$$\begin{aligned} a_n^\pm = 2\,\frac{n_f}{n}\,\frac{\gamma _n^{\pm }-\gamma _{n}}{\gamma _{QG}(n)} , \end{aligned}$$

(58)

the singlet combination of coefficients of the Gegenbauer expansion (36), and the gluon coefficients of the analogous expansion (37), read:

$$\begin{aligned}&\frac{1}{n_f} \sum _q d_n^q (t, \mu _{\mathrm{F}}^2) = d_n^{+}(t,\mu _{\mathrm{F}}^2) + d_n^{-}(t,\mu _{\mathrm{F}}^2) , \end{aligned}$$

(59)

$$\begin{aligned}&d_n^g(t,\mu _{\mathrm{F}}^2) = a_n^{+} d_n^{+}(t,\mu _{\mathrm{F}}^2) + a_n^{-} d_n^{-}(t,\mu _{\mathrm{F}}^2) . \end{aligned}$$

(60)

This decomposition is convenient since the coefficients $d_n^\pm$ evolve multiplicatively:

$$\begin{aligned} d_n^{\pm }(t,\mu _{\mathrm{F}}^2) = d_n^{\pm }(t,\mu _{\mathrm{F,0}}^2) \left( \frac{\alpha _{\mathrm{s}}(\mu _{\mathrm{F}}^2)}{\alpha _s(\mu _{\mathrm{F,0}}^2)}\right) ^{\textstyle \frac{2\gamma _n^{\pm }}{\beta _0}} . \end{aligned}$$

(61)

The inversion of the system Eqs. (59), (60) is straightforward:

$$\begin{aligned} d_n^\pm (t, \mu _{\mathrm{F,0}}^2) = \pm \frac{d_n^g(t, \mu _{\mathrm{F,0}}^2) - \displaystyle \frac{a_n^\mp }{n_f} \displaystyle \sum _q d_n^q(t,\mu _{\mathrm{F,0}}^2)}{a_n^+ - a_n^-} . \end{aligned}$$

(62)

It is useful to look at the behavior of the D-term in the $\mu _{\mathrm{F}}^2 \rightarrow \infty$ limit. Since all the anomalous dimensions $\gamma _n^{\pm }$ are positive, except $\gamma _{1}^{-} = 0$, all $d_{n}^{q}$ and $d_{n}^{g}$ coefficients vanish in this limit, except maybe the following two [30]:

$$\begin{aligned} \frac{1}{n_{f}}\sum _{q}d_1^q(t, \mu _{\mathrm{F}}^2)&\xrightarrow {\mu _{\mathrm{F}}^2 \rightarrow \infty } \, d_1^{-}(t, \mu _{\mathrm{F,0}}^2) , \end{aligned}$$

(63)

$$\begin{aligned} d_1^g(t, \mu _{\mathrm{F}}^2)&\xrightarrow {\mu _{\mathrm{F}}^2 \rightarrow \infty } \, a_1^{-} d_1^{-}(t, \mu _{\mathrm{F,0}}^2) , \end{aligned}$$

(64)

with:

$$\begin{aligned} d_1^{-}(t, \mu _{\mathrm{F,0}}^2)&= d(t, \mu _{\mathrm{F,0}}^2) \frac{T_F}{n_{f}\,T_F + 2 C_{F}} , \end{aligned}$$

(65)

$$\begin{aligned} a_1^{-} d_1^{-}(t, \mu _{\mathrm{F,0}}^2)&= d(t, \mu _{\mathrm{F,0}}^2) \frac{2 C_F}{n_{f}\,T_F + 2 C_{F}} , \end{aligned}$$

(66)

where $d(t, \mu _{\mathrm{F,0}}^2)$ has been defined in Eq. (48).

We conclude with the following remark. If the light flavours contribute equally, i.e. if $d_n^{q_1}(t, \mu _{\mathrm{F,0}}^2)$ = $d_n^{q_2}(t, \mu _{\mathrm{F,0}}^2)$ for the two light flavours $q_1$ and $q_2$, then the non-singlet combination (57) identically vanishes and the left-hand side of the singlet combination (59) reduces to the common value for the considered light flavours.

5 Modelling of DVCS subtraction constant

We have seen that an unbiased extraction at LO is theoretically possible through the lever arm in $Q^2$ of the subtraction constant, but the weak t and $Q^2$ dependence of current experimental data motivates us to introduce a simplified modelling of the D-term. Those simplifications are typical of analyses like the present one. We list below these assumptions and discuss their importance on the extraction of D-term information.

• 1. We restrict the analysis only to the first element of the Gegenbauer expansions (36) and (37), i.e. only $d_{1}^{q}(t, \mu _{\mathrm{F}}^2)$ is directly fitted.

  We may hope that these Gegenbauer expansions converge, and consequently that their general terms asymptotically vanish. However this does not justify a strict dominance of the first terms of the series; extrapolating the successive ratios ($\simeq 1/3$) of chiral quark-soliton model estimates (38)–(40) to infinite orders, we may think that retaining only $d_1$ generates a 50% systematic uncertainty. In order to fit coefficients proportional to other terms of these expansions, in particular $d_{3}^{q}(t, \mu _{\mathrm{F}}^2)$, one needs extra constraints separating various contributions, in this case either the t or $\mu _{\mathrm{F}}^2$ unique dependence. Since not much is known about the t-dependence from first principles (in particular there is no reason to think that all coefficients obey the same t-dependence), GPD evolution provides a natural tool to extract $d_{3}^{q}(t, \mu _{\mathrm{F}}^2)$ and even higher order terms. To achieve this, one needs both precise data and a large arm in $Q^2$ to study the evolution carefully. In the following we show that the current data does not allow for such kind of analysis. EIC and EIcC are natural candidates for changing this conclusion, but the quantitative assessment of their impact goes beyond the scope of the present study.

• 2. We assume an equal contribution $d_{n}^{uds}(t, \mu _{\mathrm{F}}^2)$ of light quarks to each coefficient $d_{n}^{q}$ of the D-term expansion (36), i.e.:

  $$\begin{aligned} d_{n}^{u} = d_{n}^{d} = d_{n}^{s} \equiv d_{n}^{uds} . \end{aligned}$$

  (67)

  In particular, the previous assumption focuses our study on $d_{1}^{uds}(t, \mu _{\mathrm{F}}^2)$ satisfying:

  $$\begin{aligned} d_{1}^{u} = d_{1}^{d} = d_{1}^{s} \equiv d_{1}^{uds} . \end{aligned}$$

  (68)

  This is a common assumption, met e.g. in the chiral quark–soliton model [30, 48]. We stress that disentangling the separate quark flavour contributions is not possible at LO with DVCS alone, but should be feasible for instance by studying hard exclusive meson production data. Early studies of this type [49] demonstrate the potential of such kind of analysis.

  In such a case, the relation (43) between the subtraction constant and the coefficients $d_1^q$ becomes:

  $$\begin{aligned} {\mathcal {C}}_H(t,Q^2) = \frac{8}{3} d_1^{uds}(t, Q^2) = \frac{8}{9} \sum _{q} d_{1}^{q}(t, Q^2) . \end{aligned}$$

  (69)

  As we will see in Sect. 7, some studies consider only two active quark flavours (i.e. no strange quark), which still equally contribute to the D-term. With an obvious adaptation of notations, the relation (43) now reads:

  $$\begin{aligned} {\mathcal {C}}_H(t,Q^2) = \frac{20}{9} d_1^{ud}(t, Q^2) = \frac{10}{9} \sum _{q} d_{1}^{q}(t, Q^2) . \end{aligned}$$

  (70)

• 3. Because of the absence of a direct sensitivity to the gluon D-term in a LO analysis of DVCS, $d_{1}^{g}$ is not fitted to the data.

  Instead, $d_{1}^{g}$ is radiatively generated starting from a low factorisation scale where a valence quark picture of the proton should hold. Although it contradicts the conclusion of Ref. [50], this assumption is frequently met in the computation of various parton distribution functions from quark models, and does not prevent an analysis of the existing DVCS data. Here we set:

  $$\begin{aligned} d_{1}^{g}(t, \mu _{\mathrm{F,0}}^{2}) = 0 \text { at } \mu _{\mathrm{F,0}}^{2} = 0.1~\mathrm {GeV}^2 . \end{aligned}$$

  (71)

  The sensitivity of this choice on the extraction of $d_{1}^{uds}(t, \mu _{\mathrm{F}}^2)$ is discussed in the following. The radiative generation of gluons will manifest itself with a non-vanishing $d_{1}^{g}(t, \mu _{\mathrm{F}}^{2})$ at scales $\mu _{\mathrm{F}}^2 > \mu _{\mathrm{F,0}}^{2}$.

• 4. Since a significant amount of current DVCS data possess $Q^{2}$ values larger than the squared mass of the charm quark, $m_{c}^{2} = (1.28~\pm ~0.03)^2~\mathrm {GeV}^{2}$, $d_{1}^{c}(t, \mu _{\mathrm{F}}^2)$ may contribute to the subtraction constant.

  This contribution is expected to be negligible in the considered kinematic range. However, it can be conveniently generated by the evolution equations, similarly to gluons. In this analysis we explore this possibility with the following boundary condition:

  $$\begin{aligned} d_{1}^{c}(t, m_{c}^{2}) = 0 . \end{aligned}$$

  (72)

  The radiative generation of charm quarks will manifest itself with a non-vanishing $d_{1}^{c}(t, \mu _{\mathrm{F}}^{2})$ at scales $\mu _{\mathrm{F}}^2 > m_c^{2}$.

• 5. We adopt a multipole form for the t-dependence of $d(t, \mu _{\mathrm{F}}^2)$ for $d \in \{d_{1}^{uds}, d_{3}^{uds}, d_1^{c}, d_1^{g}\}$:

  $$\begin{aligned} d(t, \mu _{\mathrm{F}}^2) = d(\mu _{\mathrm{F}}^2) \left( 1 - \displaystyle \frac{t}{\varLambda ^{2}} \right) ^{-\alpha } , \end{aligned}$$

  (73)

  where $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$, $d_{3}^{uds}(\mu _{\mathrm{F}}^2)$, $d_1^{c}(\mu _{\mathrm{F}}^2)$ and $d_1^{g}(\mu _{\mathrm{F}}^2)$ are parameters to be determined. Unless explicitly stated otherwise, the parameters $\varLambda = 0.8~\mathrm {GeV}$ and $\alpha = 3$ are kept fixed in our fits. In the following we will refer to it as the tripole Ansatz. The value of $\varLambda$ is motivated by the chiral quark-soliton model [10, 51], while that for $\alpha$ ensures a realistic shape of the pressure distribution at large t [29].

We will also explore the possibility of extracting one of the parameters $d_1^{g}(\mu _{\mathrm{F}}^2)$, $d_{3}^{uds}(\mu _{\mathrm{F}}^2)$, $\varLambda$ or $\alpha$ together with $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ in the following.

6 Input from global fits to DVCS data

The procedure described in this section is used to find replicas representing the subtraction constant from a given set of replicas representing CFFs. In particular we detail the implementation of the multipole Ansatz (73) driving the kinematic extrapolation to vanishing t.

In Ref. [52], 30 distinct DVCS observables spread over 2500 kinematic configurations and collected over 17 years were jointly analysed in terms of CFFs relying on a neural network approach. In this study the real and imaginary parts of each CFF were independently described and simultaneously fitted to experimental data. All four leading-twist and chiral-even GPDs were considered, and no assumption was made beyond the validity of a leading-twist analysis. The uncertainties of experimental data are reflected through a set of 101 replicas for each of the eight extracted functions (i.e. for real and imaginary parts of each CFF).

Let us denote the replicas associated to the real and imaginary parts of the CFF ${\mathcal {H}}$ by:

$$\begin{aligned} \mathrm {Re}\,{\mathcal {H}}_{i}^{\mathrm {NN}}(\xi , t, Q^2) \quad \text {and}\quad \mathrm {Im}\,{\mathcal {H}}_{i}^{\mathrm {NN}}(\xi , t, Q^2) , \end{aligned}$$

(74)

where $i=0, \ldots , 100$, and where the superscript $\mathrm {NN}$ indicates that a given quantity is obtained in the neural network analysis. Here, each replica is a function of the three variables $\xi$, t and $Q^2$, and represents a single neural network built of (i) three input neurons receiving the values of $\xi$, t and $Q^{2}$, (ii) one hidden layer with six neurons, and (iii) one output neuron returning either the real or imaginary part of the CFF ${\mathcal {H}}$. The parameters of this function, which are weights and biases “trained” to the experimental data, are real numbers (for details see Ref. [52]). Those numbers are not intuitive in the sense that they do not carry any physical meaning – they are not “human readable” or interpretable. The replica for $i = 0$ is obtained from a direct fit to the experimental data, while those for $i \ne 0$ are obtained from data where the central values are randomly smeared according to uncertainties. In this case i is used to distinguish between the unique values of the random seed used in the smearing procedure.

The subtraction constant is evaluated from a pair of i-th replicas ($0 \le i \le 100$) using the dispersion relation:

(75)

where the $\xi '$-integration ranges between $\epsilon = 10^{-6}$ and 1 instead of 0 and 1 as in Eq. (31). This corresponds to the kinematic domain probed by the replicas in the global CFF fit to DVCS data. This replacement of the integration range introduces a negligible ${\mathcal {O}}(1\%)$ bias, see Ref. [52] for further details. We remind that dispersion relations derive from first principles, which makes this approach model-independent.

In this analysis, the expectation value and the variance of the subtraction constant distribution at a kinematic configuration $(\xi , t, Q^2)$ are estimated with the empirical mean and the standard deviation of a replica subset of $\big (\mathrm {Re}{\mathcal {H}}_{i}^{\mathrm {NN}}, \mathrm {Im}{\mathcal {H}}_{i}^{\mathrm {NN}}\big )_{0 \le i \le 100}$ at the same kinematic configuration:

$$\begin{aligned}&\mu _{{\mathcal {C}}_{H}}(\xi , t, Q^2) = \frac{1}{|I|}\sum _{i \in I} {\mathcal {C}}_{H,i}(\xi , t, Q^2) , \end{aligned}$$

(76)

$$\begin{aligned}&\sigma _{{\mathcal {C}}_{H}}(\xi , t, Q^2) \nonumber \\&\quad = \sqrt{\frac{1}{|I|}\sum _{i \in I} \left[ {\mathcal {C}}_{H,i}(\xi , t, Q^2) - \mu _{{\mathcal {C}}_{H}}(\xi , t, Q^2)\right] ^2} , \end{aligned}$$

(77)

where $|I| \le 101$ is the number of replicas i in the subset I of $\{0, 1, \ldots , 100\}$ used in the estimation. Using most of, but not all, the available replicas is motivated by the presence of outliers: it happens sometimes that some replicas give exceptional or “exotic” values widely separated from those returned by the rest of the replica population. This typically signals problems with the supervised training of neural networks, either because some kinematic regions are not sufficiently covered by data, or e.g. when the involved training algorithm exceptionally converges to a local minimum [52]. The procedure to remove the outliers is described in the next paragraph.

It is assumed that at a given $(\xi , t, Q^2)$-point the replicas return normally distributed values, which in particular allow for a straightforward assignment of $68\%$ confidence level to the range of $\mu _{{\mathcal {C}}_{H}}(\xi , t, Q^2) \pm \sigma _{{\mathcal {C}}_{H}}(\xi , t, Q^2)$, and a simple propagation of uncertainties since the probability distribution of replicas becomes unambiguously characterised by the two parameters $\mu _{{\mathcal {C}}_{H}}(\xi , t, Q^2)$ and $\sigma _{{\mathcal {C}}_{H}}(\xi , t, Q^2)$. This assumption is in general fulfilled, however not in all cases, because of the presence of outliers. Since these issues are in general local in the $3\mathrm {D}$ phase-space of kinematic configurations $(\xi , t, Q^2)$, it is difficult to identify problematic replicas and entirely remove them from the analysis. Such procedure could also introduce a bias. Instead, the outliers are locally detected and removed with the classical three-sigma rule. This rule is known from big data analytics [53], where it is successfully used to improve the quality of data under consideration. The procedure is iterative. Starting from $I = \{0, 1, \ldots , 100\}$ (hence $|I| = 101$), it consists of the following steps:

• 1. evaluate the mean $\mu$, and standard deviation $\sigma$, of a given sample $\big ({\mathcal {C}}_{H,i}(\xi ,t,Q^2)\big )_{0 \le i \le 100}$ using Eqs. (76) and (77),

• 2. remove from that sample all replicas i for which ${\mathcal {C}}_{H,i}(\xi ,t,Q^2)$ does not lay in the $(\mu - 3\sigma , \mu + 3\sigma )$ range,

• 3a. if no elements have been removed, stop the procedure – the sample is free of outliers – and retrieve the subset I, the mean $\mu _{{\mathcal {C}}_{H}}(\xi , t, Q^2)$ and the standard deviation $\sigma _{{\mathcal {C}}_{H}}(\xi , t, Q^2)$,

• 3b. if any element has been removed, go to 1. and proceed with the next steps.

To some extent, the procedure is equivalent to a fit to a Gaussian Ansatz and a removal of all elements populating the exterior of the $\mu \pm 3\sigma$ range.

We will now fit the multipole Ansatz, see Eq. (73), to the subtraction constant replicas obtained in the neural network analysis. We work with the kinematic configuration $(\xi _{j}, t_{j}, Q^{2}_{j})_{1 \le j \le N_{\mathrm {pts}}}$ where DVCS measurements have been published, with $\xi = x_{\mathrm {Bj}}/(2-x_{\mathrm {Bj}})$ in the Bjorken regime and $N_{\mathrm {pts}} = 277$. We will explore various fitting scenarios where the extracted parameters constitute a subset Y of $\{d_{1}^{uds}(\mu _{\mathrm{F}}^2), d_{3}^{uds}(\mu _{\mathrm{F}}^2), d_{1}^{g}(\mu _{\mathrm{F}}^2), \varLambda , \alpha \}$ for some hadronic scale $\mu _{\mathrm{F}}^2 > m_c^{2}$. For each subtraction constant replica ${\mathcal {C}}_{H,i}$ with $i \in \{0, 1, \ldots , 100\}$, we define the $\chi ^{2}$-function $\chi _{i}^{2}(Y)$ for the generated points in the following way:

$$\begin{aligned}&\chi _{i}^{2}(Y) \nonumber \\&\quad = \sum _{j = 1}^{N_{\mathrm {pts}}} \left( \frac{ {\mathcal {C}}_{H,i}(\xi _{j}, t_{j}, Q^{2}_{j}) - {\mathcal {C}}_{H, i}^{\mathrm {\alpha P}}(\xi _{j}, t_{j}, Q^{2}_{j}, Y) }{\sigma _{{\mathcal {C}}_{H}}(\xi _{j}, t_{j}, Q^{2}_{j})} \right) ^{2} , \end{aligned}$$

(78)

where $(\xi _{j}, t_{j}, Q^{2}_{j})$ is the kinematic configuration generated for the point j, ${\mathcal {C}}_{H,i}(\xi _{j}, t_{j}, Q^{2}_{j})$ and $\sigma _{{\mathcal {C}}_{H}}(\xi _{j}, t_{j}, Q^{2}_{j})$ are, respectively, the neural network analysis values of the subtraction constant returned by replica i and its uncertainty at point j, and ${\mathcal {C}}_{H, i}^{\mathrm {\alpha P}}(\xi _{j}, t_{j}, Q^{2}_{j}, Y)$ is the value of the subtraction constant from the multipole Ansatz also at point j. The $\chi ^{2}$-function $\chi _{i}^{2}(Y)$ allows one to constrain the free parameters Y of the tripole or multipole Ansatz. In contrast to the neural network analysis, those parameters possess a physical interpretation, and in particular straightforwardly provide the value of the subtraction constant at $t=0$.

7 Results

7.1 Nominal fitting scenario

In this analysis, we first consider a nominal fitting scenario where only one parameter is directly fitted to data parameterised with neural networks, namely $d_{1}^{uds}(\mu _{\mathrm{F}}^{2})$. The values for gluons, $d_{1}^{g}(\mu _{\mathrm{F}}^{2})$, and charm quarks, $d_{1}^{c}(\mu _{\mathrm{F}}^{2})$, are indirectly fitted, i.e. they are generated “radiatively” by the LO evolution using the following boundary conditions:

$$\begin{aligned}&d_{1}^{g}(\mu _{\mathrm {F}, 0}^{2}) = 0 , \end{aligned}$$

(79)

$$\begin{aligned}&d_{1}^{c}(m_{c}^{2}) = 0 , \end{aligned}$$

(80)

where $\mu _{\mathrm {F}, 0}^{2} = 0.1~\mathrm {GeV}^{2}$ is the initial factorisation scale squared and where $m_{c}^{2} \simeq 1.64~\mathrm {GeV}^{2}$ is the mass of charm quark squared. For details see Sect. 5 above. The obtained numerical values are summarised in Table 1, where the values of $d_{1}^{uds}$, $d_{1}^{g}$ and $d_{1}^{c}$ are quoted at the typical hadronic scale $\mu _{\mathrm{F}}^2 = 2~\mathrm {GeV}^2$. The scale independent sum over partons, $d(t=0)$, is given as well.

Table 1 Values of $d_{1}^{uds}$, $d_{1}^c$ and $d_{1}^g$ for light, charm quarks and gluons quoted at $\mu _{\mathrm{F}}^2 = 2~\mathrm {GeV}^2$. Those values are obtained with a tripole Ansatz, where only $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ is directly fitted to experimental data parameterised with neural networks. In addition, given is the scale independent sum over partons defined in Eq. (48)

Means and uncertainties of the extracted parameters are estimated using Eqs. (76) and (77), respectively. The dependence of the subtraction constant on $\xi$, t and $Q^{2}$ is shown in Fig. 1 for both the tripole Ansatz and neural network replicas. One may notice that the uncertainties obtained with the tripole Ansatz are typically smaller than those obtained in the neural network analysis, which is expected and reflects the internal consistency of our approach. This is particularly visible in domains where the subtraction constant is strongly driven by the assumed functional form of the Ansatz, like in particular for large |t|. In Fig. 1 the difference between the confidence levels indicates the model uncertainty of the tripole Ansatz. The dependence of $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$, $d_{1}^{g}(\mu _{\mathrm{F}}^2)$ and $d_{1}^{c}(\mu _{\mathrm{F}}^2)$ at $\mu _{\mathrm{F}}^2 = 2~\mathrm {GeV}^2$ on the choice of the initial factorisation scale, $\mu _{\mathrm {F}, 0}^{2}$ is shown in Fig. 2. As one may notice, the impact of this choice on the extracted value of $d_{1}^{uds}(\mu _{\mathrm{F}}^2 = 2~\mathrm {GeV}^{2})$ is negligible compared to the estimated statistical uncertainties.

Fig. 1

Subtraction constant for the CFF ${\mathcal {H}}$ as a function of $\xi$ at $t = -\,0.3~\mathrm {GeV}^{2}$ and $Q^{2} = 2~\mathrm {GeV}^{2}$ (top), as a function of $-t$ at $\xi = 0.2$ and $Q^{2} = 2~\mathrm {GeV}^{2}$ (middle) and as a function of $Q^{2}$ at $\xi = 0.2$ and $t = -\,0.3~\mathrm {GeV}^{2}$ (bottom). The solid gray bands represent results of the neural network analysis of Ref. [52] while the green punctured bands represent the analysis based on a tripole Ansatz described in the text. All bands correspond to a $68 \%$ confidence level

Fig. 2

Values of $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ for light quarks (blue backward-hatched band – $\diagdown \diagdown$), $d_{1}^{c}(\mu _{\mathrm{F}}^2)$ for charm quarks (brown solid band) and $d_{1}^{g}(\mu _{\mathrm{F}}^2)$ for gluons (red forward-hatched band – $\diagup \diagup$) quoted at $\mu _{\mathrm{F}}^2 = 2~\mathrm {GeV}^2$ as a function of the initial factorisation scale squared $\mu _{\mathrm {F}, 0}^{2}$. The results for charm quarks are multiplied by 100. All bands correspond to $68 \%$ confidence level

Although error bars are too large to make a firm claim, we observe that a negative mean value of $d_{1}^{uds}(\mu _{\mathrm{F}}^2 = 2~\mathrm {GeV}^2)$ seems to be favoured by experimental data, as expected and observed in all stable systems [11]. The results of other phenomenological and theoretical analyses, including lattice-QCD predictions, are collected in Table 2. Their comparison to our results is summarised Fig. 3.

Table 2 Compilation of results for $\sum _{q} d_{1}^{q}(\mu _{\mathrm{F}}^{2})$ over a given number of quark flavours. Results (1) and (2) are originally only given for the DVCS subtraction constant considering three light flavours, and here are scaled by a factor 9/8 following Eq. (69). The dispersive evaluations (3) differ by the input pion PDF. The lattice results (9) are originally only given for the EMT form factor $C_{u}(0, \mu _{\mathrm{F}}^{2})+C_{d}(0, \mu _{\mathrm{F}}^{2})$, and here are scaled by a factor 5 following Eq. (46). Both differ by the extrapolation of lattice data to the chiral limit. The scale associated to all results coming from $\chi$QSM is assumed to be $0.6~\mathrm {GeV}$, which is the natural scale for this type of models as argued in Ref. [54]. The same scale is associated to the Skyrme model

Fig. 3

The sum $\sum _{q} d_{1}^{q}(\mu _{\mathrm{F}}^{2})$ as a function of $\mu _{\mathrm{F}}^{2}$ for this study (green band) and other phenomenological and theoretical analyses. See Table 2 for the description of each data point, including the marker legend. For the sake of legibility some markers are artificially shifted by a small distance in the horizontal direction

7.2 Relaxing constraints on the fitting Ansatz

We also tried to fit more than one parameter in the multipole Ansatz, namely $\{d_{3}^{uds}(\mu _{\mathrm{F}}^2), d_{1}^{g}(\mu _{\mathrm{F}}^2), \varLambda , \alpha \}$. The attempts at extracting two parameters at the same time are summarised in Table 3. We first observe that all five extracted values of $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ (from Tables 1, 3) are consistent and compatible with zero. The ratio of the extracted mean to standard deviation vary over an order of magnitude, between 0.04 and 0.5 in all five fitting scenarios, confirming the difficulty of extracting a statistically significant value of $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ from existing DVCS data.

Table 3 Results for four scenarios when $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ and one other parameter in $\{d_{3}^{uds}(\mu _{\mathrm{F}}^2), d_{1}^{g}(\mu _{\mathrm{F}}^2), \varLambda , \alpha \}$ are extracted from experimental data. This pair of parameters is noted Y. The values of are quoted at $\mu _{\mathrm{F}}^2 = 2~\mathrm {GeV}^2$. The last column indicates ranges in which parameters are allowed to vary in the fit procedure. If the range is not specified, the corresponding parameter is allowed to vary between $-\infty$ and $\infty$

The extraction introducing $d_{3}^{uds}(\mu _{\mathrm{F}}^2)$ as a free parameter has a somewhat different status from those proposing $\varLambda$ or $\alpha$ as their free parameters. While fixing $\varLambda$ and $\alpha$ imposes a parametric form to capture the t-dependence of experimental data, fitting only $d_1^{uds}(\mu _{\mathrm{F}}^2)$ represents a strong assumption on the behavior of the $(d_n^q)_{n \, \text {odd}}$ series. As mentioned in Sect. 5, this series has no expected reason to converge quickly and even less to stop after its first term. If a reasonable model can guide the choice of a specific t-dependence, a contrario neglecting $d_{3}^{uds}(\mu _{\mathrm{F}}^2), d_{5}^{uds}(\mu _{\mathrm{F}}^2), \ldots$ in the extraction may leave a major uncontrolled source of systematic uncertainties.

In this respect the four scenarios restricting the series (36) to its first term provide us values of $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ in good agreement. The multipole parameters $\varLambda$ and $\alpha$ are still subject to large uncertainties. However, we observe a dramatic change in the estimated mean values and standard deviations between the fits extracting $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ only on the one hand, and $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ and $d_{3}^{uds}(\mu _{\mathrm{F}}^2)$ simultaneously on the other hand. This signals a possible influence of higher-order terms in the series (36) and raises the question of determining the terms of the series (36) from the values of the subtraction constant and the specific $Q^2$-dependence of each term.

Separating $d_1$ from the higher-order terms $d_3, d_5, \ldots$ consists in carefully identifying in the data their own logarithmic $Q^2$-dependence, driven by distinct anomalous dimensions. This requires a large lever arm in $Q^2$ and very accurate data, making this study best adapted to collider settings like EIC or EIcC. We observe that the joint fit of $d_{1}^{uds}(\mu _{\mathrm{F}}^2)$ and $d_{3}^{uds}(\mu _{\mathrm{F}}^2)$ to existing measurements exhibits a large correlation between these two parameters with a Pearson’s correlation coefficient of $-\,0.997$, as shown in Fig. 4. This plot summarises the difference between fitting $d_{1}^{uds}$ only, or $d_{1}^{uds}$, $d_{3}^{uds}$ simultaneously and depicts the aforementioned change between these two fitting scenarios. It also shows that the accuracy of existing data cannot exclude the case where higher-order terms in the series (36) decrease slowly.

Fig. 4

Correlation between $d_1^{uds}(\mu _{\mathrm{F,0}}^2)$ and $d_3^{uds}(\mu _{\mathrm{F,0}}^2)$ obtained in the fitting scenario where both parameters are released at the same time. Each point corresponds to a single replica used in the extraction. Confidence levels of $68\%$ and $95\%$ are denoted by dashed and dotted curves, respectively. Pearson’s correlation coefficient is estimated to be $-\,0.997$. The result of the extraction where only $d_1^{uds}(\mu _{\mathrm{F,0}}^2)$ is fitted (and where $d_3^{uds}(\mu _{\mathrm{F,0}}^2) = 0$) is denoted by the green triangle

7.3 Discussion

Even if our results tend to favour a negative mean value for the first coefficient of the Gegenbauer expansion of the D-term $\sum _q d_{1}^{q}$ at $t=0$, we stress again that the estimated standard deviation makes it statistically compatible with 0.

While we used a non-parametric global fit of almost all existing DVCS experimental data to obtain this result, the authors of Ref. [13] reported a value of the same coefficient statistically incompatible with 0 (see Table 2) when using only a subset of the existing DVCS measurements (beam-spin asymmetries and cross sections published by CLAS in 2008 [60] and 2015 [61]) and a parametric fit of CFFs inspired by the KM model, see Ref. [37] and references therein. Extracting an apparently more accurate estimate of $d_1$ from less experimental data may come as a surprise, but it merely reflects the impact of systematic assumptions made in this type of studies and discussed in previous sections.

The conclusions of our nominal fitting scenario are in both quantitative and qualitative agreement with those of Ref. [14], which also rely on a neural network extraction of CFFs: $\sum _q d_{1}^{q}$ is found compatible with 0, with a mean to standard deviation ratio of about 0.5 (see Table 2). The difference between the results reported here and in Ref. [14] that is observed in Fig. 3 may be explained by a different extrapolation to $t=0$. We conducted here a careful study to evaluate the systematic uncertainties introduced by some classical assumptions, and we found no options allowing for a more accurate determination of $d_{1}^{uds}$ from the DVCS experimental data that have been collected so far.

Our result being consistent with 0 within error bars, we will not perform the extraction of a pressure anisotropy profile from experimental data along the lines of Sect. 2. In the same spirit, we do not attempt to evaluate the proton mechanical radius, which, from Eq. (20), is essentially the second moment of this pressure profile, and would come compatible with 0 as well.

We nevertheless make a simple remark: if we had performed such an extraction, the distribution of pressure forces would be encoded in the t-dependence of $d_1$, which in current studies is described by the classical multipole Ansatz (73). Stated differently, the qualitative shape of the pressure profile is selected before any fitting is realised, since it is merely the Fourier transform of a multipole. In this context the shape of the pressure profile is a model assumption, and this shape may vary with the choice of the multipole parameters $\varLambda$ and $\alpha$, and the overall normalisation $d_1(t=0)$.

We use Eq. (49) to point out the dependence of this profile on the assumed values of the tripole Ansatz parameters considering the results of Table 3. Since we did not estimate $\varLambda$ and $\alpha$ simultaneously from the data, we do not know the correlation matrix of these parameters. Thus we performed an impact study by varying only one parameter and keeping the other fixed. This may allow us to explore a larger parametric space than otherwise permitted by DVCS data, but the exercise reported in Fig. 5 demonstrates variations in the pressure profile over orders of magnitude. As a consequence, all phenomenological information about the distribution of pressure forces extracted to this day from current experimental data should be taken with great care. As emphasised above, future data from EIC or EIcC should improve these phenomenological results and clarify the situation.

Fig. 5

Profiles of pressure anisotropies for a single parton, see Eqs. (49) and (50), using the multipole Ansatz of $d^{q,g}_{1}(t)$, see Eq. (73), for various values of the parameters $\alpha$ and $\varLambda$. The upper plot corresponds to $\varLambda = 0.8~\mathrm {GeV}$ and various values of $\alpha$, while the lower one assumes $\alpha = 3$ and shows various values of $\varLambda$. In both cases the same normalisation factor is used: $d_{1}(0) = -3$

The improvement of the theoretical description of DVCS through the inclusion of NLO corrections [41, 62] may be needed for a reliable extraction of the EMT form factors. These corrections will change the relation between the subtraction constant and the D-term, which for LO is expressed by Eq. (35). A full NLO analysis is beyond the scope of this paper, but should be considered when more precise data becomes available. In addition, the higher twist corrections [63] may influence the extraction of the subtraction constant from experimental data. In the present analysis, the kinematic domain where these corrections are expected to be significant is rejected at the level of data selection [52]. The inclusion of higher-twist corrections may allow us to incorporate more data and extend the available range of $t = -\varDelta ^2$, which would be valuable since $\varDelta$ is Fourier-conjugated to the distance to the center of the nucleon. It is also beyond the scope of this paper, but should be considered as a possible improvement of this type of analyses.

8 Summary

This article summarises our analysis of the $d_{1}(\mu _{\mathrm{F}}^{2})$ coefficient extracted from the world experimental data for the DVCS process. This coefficient is of utmost importance for the understanding of the QCD energy–momentum tensor. It carries information on mechanical properties of partonic media, like shear stress and pressure.

We carefully detailed the extraction procedure, and plainly exhibited a set of assumptions that we made and that have been used in other analyses of similar type until now. The underlying motivation is the assessment of what can quantitatively be said to this day about the proton mechanical properties using DVCS measurements. We hope it will be a useful starting point for future studies when new data become available. We tried to make this document as self-contained as possible, including all equations with consistent notations, covering the whole reasoning leading to the extraction of pressure distributions, from proton structure to DVCS measurements. In particular, the evolution equations of the $d_{n}^{q,g}(\mu _{\mathrm{F}}^{2})$ coefficients are presented in a complete form and are used to estimate the contributions coming from gluons and charm quarks.

The comparison between analyses where either artificial neural networks or parametric functional forms are used to describe CFFs clearly demonstrates the model dependence of the latter approach, which results in considerably smaller uncertainty for the extracted $d_{1}^{uds}(\mu _{\mathrm{F}}^{2})$ coefficient, making it no longer compatible with 0. In such a case the systematic uncertainty associated to the choice of a specific functional form is difficult to evaluate and is left as an unknown. The extraction of $d_{1}^{uds}(\mu _{\mathrm{F}}^{2})$ independently from elements proportional to higher terms in the D-term Gegenbauer expansion, like $d_{3}^{uds}(\mu _{\mathrm{F}}^{2})$, cannot be performed at this moment. Removing the assumption regarding the symmetry of the quark coefficients $d_{n}^{uds}(\mu _{\mathrm{F}}^{2})$ most probably requires a multi-channel analysis.

Even if the link between the distribution of pressure forces in the proton and the DVCS subtraction constant is well-defined, and even if this link is subject to approximations that can be checked systematically, DVCS data do not allow yet a statistically significant extraction of these pressure forces. Moreover the present need of a prior knowledge of the t-dependence of the D-term drives the pressure profile. All phenomenological information about the distribution of pressure forces extracted to this day from current experimental data should be taken with great care. Our analysis establishes the need for more precise data and for an extension of the covered kinematic domain. This can be achieved by future experiments to be conducted at JLab, CERN, EIC and EIcC facilities.

Appendix A: Open source code

Replicas used in this analysis representing the subtraction constant parameterised with the tripole Ansatz are available in PARTONS framework [64] as the module GPDSubtractionConstantDLMSTW21. The code of this framework is open source and can be found online at https://drf-gitlab.cea.fr/partons/core/partons on version 3 of the GPL (GPLv3).