Effects of long-range forces on the D-term and the energy-momentum structure

The hadronic form factors of the energy-momentum tensor (EMT) have attracted considerable interest in recent literature. This concerns especially the D-term form factor D(t) with its appealing interpretation in terms of internal forces. With their focus on hadron structure, theoretical studies so far have concentrated on strongly interacting systems with short-range forces. Effects on the EMT due to long-range forces like the electromagnetic interaction have not yet been studied. Electromagnetic forces play a small role in the balance of forces inside the proton, but their long-range nature introduces new features which are not present in systems with short-range forces. We use a simple but consistent classical field theoretical model of the proton to show how the presence of long-range forces alters some notions taken for granted in short-range systems. Our results imply that a more careful definition of the D-term is required when long-range forces are present.


I. INTRODUCTION
The matrix elements the EMT, T µν , [1] can be explored through studies of generalized parton distribution functions in hard exclusive reactions [2,3] and contain information on the basic properties of a particle: mass, spin, and the equally important but far less known D-term [4]. The information content of EMT form factors is visualized in terms of EMT densities [5] which allow us to learn about properties like energy density, angular momentum distribution, or internal forces in hadrons [5][6][7][8][9][10][11]. EMT properties were studied in hadronic models, chiral perturbation theory, lattice QCD and other strongly interacting systems  which had one common feature: these systems were governed by short-range forces. The goal of this work is to investigate the impact of long-range forces on the EMT properties.
For our study, we employ a classical model of the proton which is of interest for its own sake. Classical models of an extended electric charge have a long history dating back to the works of Abraham and Lorentz [52,53]. It was recognized by Poincaré that in order to compensate the electrostatic repulsion one must introduce cohesive forces, known as Poincaré stresses [54], which were introduced in an ad hoc manner [52][53][54][55][56][57][58]. The model of Bia lynicki-Birula [59] used in this work is to the best of our knowledge the first fully consistent classical model of an extended charged particle where the Poincaré stresses are generated dynamically in a local, relativistic, classical field theory.
In this model, "dust particles" carry an electric charge e, and strong charges g S and g V and interact with the electromagnetic 4-potential A µ , and strong scalar and vector fields, φ and V µ . The attractive (due to φ) and repulsive (due to V µ and A µ ) forces on the dust particles exactly compensate each other, such that the dust particles are in stable, static equilibrium and occupy a finite spherically symmetric region of radius R. Using nuclear phenomenology to fix model parameters, the model can describe a particle with the charge, mass, and size of the proton [59].
This classical system is well-suited for our purposes. It exhibits strong short-range forces which play an overwhelmingly important role in the internal structure of the proton, and at the same time consistently includes the effects of the long-range electromagnetic field. Two aspects are of importance for our study, namely (i) an internally consistent theoretical description of a stable particle, and (ii) the correct description of the long-range electromagnetic effects. The model of Ref. [59] satisfies both requirements. The classical aspect of the model is not a hindrance. Rather it is a virtue allowing us to investigate the effects of long-range forces undistracted by technical difficulties associated with computations in more realistic strongly interacting quantum systems.
The outline of this work is as follows: In Sec. II we briefly introduce the model, and apply it to the description of EMT densities in Sec. III, showing that the model is consistent and, in the region of r 2 fm, in good agreement with results from studies of systems with strong short-range forces. In Sec. IV we focus on distances beyond r 2 fm where new features appear which were not encountered before in systems with short-range forces. In Sec. V we show that our results regarding the long-distance properties of the EMT densities are model independent, compute the form factor D(t), and discuss the implications for experimental measurements and theoretical calculations of the D-term. The Sec. VI contains the conclusions.

II. THE CLASSICAL MODEL
In this section, we briefly introduce the model of Ref. [59] which consists of "dust particles" bound in a spherically symmetric region of radius R by the interplay of three types of fields: a massive scalar field φ, a massive vector field V µ and an electromagnetic field described by the 4-potential A µ . The particles couple to these fields respectively through the coupling constants g S , g V , e. The motion of the particles is described by a scalar phase-space distribution Γ( r, p, t). The system is defined by the following classical field equations (we use = c = 1 unless otherwise stated) The force F = f /u 0 is expressed in terms of the components of the 4-force The 4-velocity u α is defined by p α = m u α = (E p , p) and v = u/u 0 = p/E p . ∇ r and ∇ p denote derivatives with respect to positions r and momenta p of the particles. The 4-current j α and scalar density ρ are defined in terms of the phase-space distribution as Despite the non-covariant appearance of Eq. (1) the theory described by Eqs. (1-4) is relativistically invariant [59], and is a generalization of the Vlasov-Maxwell equations used in plasma physics [61]. The solution of Eqs. (1-4) is most conveniently expressed in the static case, where the particles are at rest with u α = (1, 0, 0, 0) and described by the phase space distribution Γ( r, p, t) = δ (3) ( p) ρ(r) with r = | r|. In this frame, the scalar density ρ and zeroth component of j α coincide, A α and V α only have zero components, and Eqs. (1-4) become Notice that the condition (7) provides a constraint on the fields only for r ≤ R where matter is present (i.e. ρ = 0), and is trivially satisfied in the region r > R with no matter (where ρ = 0). The density is normalized as d 3 r ρ(r) = 1.
In the region r > R the solutions of Eqs. (1)(2)(3)(4) are given by ρ(r) = 0, (16) eA 0 (r) = e 2 4πr (17) Eq. (16) means there are no particles outside the radius R, and Eq. (17) is the Coulomb potential. The six parameters b V , b S , d + , d − , 2E B , R are fixed by requiring the fields A 0 (r), V 0 (r), φ(r) to be continuous and differentiable at r = R. The reason why the constant 2E B has been named in this peculiar way will become clear shortly.
In order to apply the model to the description of the proton the following parameters were used in Ref. [59] where for convenience the constants c = 197 MeV fm are restored such that in all expressions, r is in units of fm, energies in units of MeV, ρ(r) in units of fm −3 , etc. The parameters m S and m V correspond respectively to the masses of a σ-meson and ω-meson as used in nuclear matter models. The coupling constants g S , g V are taken from the model QHD-I of the mean field theory of nuclear matter [62]. This means that in this model the proton is bound by nuclear forces [59]. For completeness, we remark that with these parameters, the requirements of continuity and differentiability of φ(r), The parameter E B = −15.71 MeV is to be confronted with the value of the bulk binding energy per nucleon in nuclear matter of −15.75 MeV [62]. Since the electric charge density is given by e ρ(r), the electric mean square radius is given by r 2 ch = d 3 r r 2 ρ(r) (recall that d 3 r ρ(r) = 1). The model yields r 2 ch 1/2 = 0.714 fm, which underestimates the experimental value by 20 % but has the right order of magnitude. This model could be elaborated to give a more realistic description. However, the modest goal of Ref. [59] was to show that the model of the proton with the parameters (20) is "not completely out of touch with reality." This is sufficient for our purposes.
In Fig. 1 we show the density ρ(r), and the potentials φ(r), V 0 (r), A 0 (r) as functions of r. The jump in the matter distribution ρ(r) shows that the dust particles are inside the radius R = 1.05 fm and there is no matter outside. The potentials in Figs. 1b-d are scaled with their respective coupling constants such that they all have the same unit MeV and can be compared. The scalar potential g S φ(r) and vector potential g V V 0 (r) are associated with strong forces and are 2 orders of magnitude larger than the Coulomb potential eA 0 (r). In the inner region, g S φ(r) is somewhat larger than g V V 0 (r). At a larger r it is the opposite: both potentials decay at a rate of ∼ exp(−m i r)/r with m i = m S and m V respectively, but m S < m V so the more massive V 0 (r) has a shorter range compared to φ(r). The Coulomb potential eA 0 (r) is small in the interior region, but becomes the dominant field in the outer region thanks to its long range A 0 (r) ∼ 1/r. These results have already been discussed in Ref. [59]. In the next section we will discuss how the different fields contribute to the mass of the proton, and how the internal forces inside the proton balance each other.

III. THE ENERGY MOMENTUM TENSOR OF THE CLASSICAL MODEL
The energy momentum tensor in the classical model of the proton was derived in Ref. [59] and is given by This expression needs to be evaluated for the static solution where u α = (1, 0, 0, 0), φ = φ(r), V α = (V 0 (r), 0, 0, 0), A α = (A 0 (r), 0, 0, 0). In the following, we discuss the different components of the EMT, and focus initially on the region r ≤ 2 fm. The long-distance properties of the EMT at r > 2 fm will be discussed in the next section.
A. Energy density T00(r) The 00-component of the EMT describes the energy density. Evaluating the 00-component in Eq. (21) yields The mass of the solution is defined as M = d 3 r T 00 (r). Integrating the energy density (22) over space, exploring the normalization d 3 r ρ(r) = 1, performing partial integrations and using Eqs. (8)(9)(10), one obtains [59] The potentials g S φ(r), g V V 0 (r), eA 0 (r) are positive, see Fig. 1. The intermediate step in (23) shows that the scalar field makes a negative contribution to the total energy, lowers the binding energy, and makes the system more strongly bound. In contrast to this, g V V 0 (r) and eA 0 (r) enter with positive signs, i.e. make the system less strongly bound. These observations are not surprising and reflect the well-known facts that scalar forces are attractive, and vector forces (for equal sign charges) repulsive. We will come back to this point below when discussing the stress tensor. As mentioned in Sec. II, numerically, E B = −15.71 MeV. This value can be compared to the bulk binding energy per nucleon in nuclear matter [59]. We remark that alternatively one could define the parameter m in (20) to be a "bare nucleon mass" such that m + E B would be the physical mass of the free proton. Here we use the original version of the model as formulated in Ref. [59] and refrain from such a redefinition of model parameters.
It is instructive to discuss the individual contributions to the energy density which we define as follows T dust 00 (r) = m ρ(r) , such that they add up to the total T 00 (r) in Eq. (22). Recalling that d 3 r ρ(r) = 1, we see that the dust particles contribute m and by far the most to M . The relatively small value E B = −15.71 MeV may give the incorrect impression that the contributions from the fields to M are small. However, these contributions are given by Thus, the relatively small value of the binding energy is the result of large cancellations between different contributions. T scal 00 (r) is the only contribution which exhibits a discontinuity at r = R and is negative in the inner region. The energy density and its individual contributions are plotted in Fig. 2a.

B. The T0k components
For the 0k-components of the EMT we obtain T 0k = 0 which is not surprising. This is because we deal with a static solution. Since there is no rotation in the system, the classical angular momentum J i = d 3 r ijk x j T 0k of the system is zero. At this point, we disregard the fact that the proton has spin 1 2 . In principle, if interested, one could treat our classical solution as a soliton and use standard quantization techniques to assign a definite spin [63].

C. The stress tensor Tij
Finally, for the ij-components of the EMT we obtain the result In general, the stress tensor can be expressed in terms of a traceless part associated with shear forces s(r) and a trace associated with the pressure p(r) as follows where e i r is the unit vector in the radial direction. The model results for s(r) and p(r) are given by Due to the EMT conservation, ∂ µ T µν = 0, the shear forces and pressure are not independent of each other but connected by the differential equation [8] (we leave here the space dimension n = 3 general for later purposes) Another consequence of the EMT conservation is the von Laue condition This is a necessary (but not sufficient) condition for stability, and requires that internal forces inside a system must exactly balance each other. If the integral in (31) was positive (negative), then the system would explode (implode).
Notice that the dust particles do not contribute to the stress tensor densities s(r) and p(r). This is naturally explained in the hydrodynamic interpretation of the model, where the dust particles can be viewed as an ideal pressureless fluid of density ρ(r), which flows (without dissipation) with the 4-velocity u µ [59]. Therefore, the densities s(r) and p(r) receive contributions only from the fields.
The scalar field makes the largest contribution to s(r), see Fig. 2b, which is positive. The contributions from vector and Coulomb fields are both negative. Not surprisingly, the Coulomb field contribution is rather small. The shear forces behave like s(r) ∝ r 2 at small r < 0.1 fm, and exhibit a global maximum at r = 0.711 fm, which is numerically close (within 0.1 %) but not the same value as the charge radius r 2 ch 1/2 = 0.714 fm. For a large nucleus in the liquid drop model, s(r) would be a δ-function centered at the edge of the nucleus (with the coefficient in front of the δ-function given by the surface tension) [8]. The result for s(r) remotely resembles a strongly smeared out δ-function. This reflects the fact that the proton has no sharp edge, and is a much more diffuse object than a nucleus.
The pressure is positive in the inner region, changes sign at r = 0.788 fm and is negative thereafter, see Fig. 2c. The sign convention is such that p(r) > 0 means repulsive forces are directed towards the outside, while p(r) < 0 means attractive forces are directed towards the inside. The shape of the total pressure distribution is largely due to the cancellation between the large contributions from scalar and vector fields. In the inner region, the repulsive vector forces are stronger than the attractive scalar forces. In the outer region, it is vice versa since, due to m V > m S , the range of the vector forces is shorter. Throughout the region plotted in Fig. 2c, the Coulomb contribution plays a minor role (but is not negligible, see below) and contributes to p(r) with the same sign as the vector field.
In order to attest the consistency of our calculation, we notice that inserting the expressions (28, 29) for s(r) and p(r) into Eq. (30) yields 2 3 s (r) + 2 r s(r) + p (r) = e r [ρ F ] = 0 due to Eq. (7). Also, the von Laue condition (31) holds which was proven in [59], and is illustrated in Fig. 2c. If we define the individual contributions to the pressure as These results mean that with scalar forces alone, the system would implode, while with vector (or Coulomb) forces alone it would explode. Clearly, in Eq. (33) the contribution of the Coulomb force is minuscule compared to that of the strong forces, but not negligible, for this system would implode without the Coulomb force. Some comments regarding the size of the forces are in order. In our model, the pressure in the center of the proton is about 20 MeV/fm 3 . This is about an order of magnitude less than in the chiral quark soliton model [16]. This result is expected and plausible for the following reason. The forces in the model of Ref. [16] are the strong forces acting between quarks. In contrast to this here, the strong forces (the massive scalar and vector fields) are modeled using nuclear physics phenomenology. Such "residual nuclear forces" are about an order of magnitude weaker than the strong forces among quarks inside the proton, and this is what we observe.
To make an intermediate summary, the description of the EMT in the classical proton model is internally consistent since the relations (30,31) hold. The size of the internal forces is what one would expect from a model with forces which have the strength of residual nuclear forces. The results discussed so far reflect the same features as encountered in other EMT studies in strongly interacting systems, including the sign patterns for the EMT densities in Fig. 2.

IV. EMT DENSITIES AND LONG-RANGE FORCES
The particular feature of the model used in this work is that it explicitly includes long-range (Coulomb) forces. Before we study this aspect in detail, let us briefly review several common features observed in prior EMT studies of strongly interacting systems governed by short-range forces [8]. For ground state solutions, in systems governed by strong short-range forces, the following common features were observed so far: 1. the shear forces s(r) are positive at all r, 2. the pressure p(r) exhibits one node at r 0 with p(r) > 0 for r < r 0 and p(r) < 0 for r > r 0 , 3. the combination 2 3 s(r) + p(r), which is referred to as normal force (per unit area), is always positive, i.e.
Some comments are in order. We are not aware of a rigorous proof of the property (i), though it is plausible given the connection of s(r) to surface tension and surface energy, which are positive in stable hydrostatic systems [8]. The positivity of s(r) was observed in all studies so far. The property (ii) arises because p(r) must have at least one node to comply with the von Laue condition (31), and ground states exhibit a single node. The pattern in Fig. 2d follows from mechanical stability arguments: repulsive forces are required in the inner region to prevent collapse and attractive forces in the outer region to bind the system [8]. For excited states, the pressure can exhibit several nodes, but the pattern with p(r) > 0 in the center and p(r) < 0 at large distances remains [24]. The property (iii) is a mechanical stability criterion and means that the radial forces T ij dA j r , where dA j r = e j r r 2 dΩ, are directed towards the outside, and the point where they vanish (if we deal with a finite size system) marks the "edge" of the system [8,49].
As long as we consider distances r 2 fm, the EMT densities in the classical proton model exhibit the properties (i-iii) as observed in prior studies. But the situation changes when we consider distances r 2 fm.

A. Long-distance behavior of the EMT densities
The dust distribution is confined to the region r < R = 1.05 fm and anyway does not contribute to s(r) and p(r). The behavior of the EMT densities at long distances is therefore determined by the fields. The contributions of the fields representing the strong forces decay exponentially at large distances r 2 fm. Despite being very small in the inner region, the Coulomb contribution becomes comparable to the contributions of the strong fields at r ∼ 2-3 fm. The Coulomb contribution is the dominating field at long-distances r 3 fm due to the slow, power-like 1 r -decay of the Coulomb potential. This is illustrated in Fig. 3a for the pressure; the situation is very similar for T 00 (r) and s(r). Using the fine-structure constant in (20) we read off from Eqs. (22,28,29) the long-distance behavior of the densities where the dots indicate subleading, exponentially suppressed contributions from the strong interaction fields. The long-distance behavior (35) of T 00 (r) does not show anything unusual: T 00 (r) > 0 for all 0 ≤ r < ∞ which was also observed in all prior studies. It is noteworthy that the 1/r 4 -decay of T 00 (r) at large r guarantees the convergence of the total energy M = d 3 r T 00 (r). But the mean square radius of the energy density r 2 E = d 3 r r 2 T 00 (r)/M diverges and cannot be defined in this model. New features emerge for the stress tensor densities s(r) and p(r). From Eq. (35) we see that s(r) is negative at large r. The classical proton model is in agreement with the property (i) and exhibits a positive s(r), see Fig. 2b, throughout the region 0 < r < r 0,s . But s(r) changes sign at the point r 0,s = 2.144 fm, and remains negative for r > r 0,s . The asymptotic expression for s(r) in Eq. (35) works with an accuracy of 2 % or better for distances r 3.1 fm.
Another new feature is that p(r) is positive at large r, see Eq. (35). Throughout the region 0 < r < r 0,p , the model conforms to the property (ii) and p(r) exhibits the characteristic pattern: positive p(r) in the inner region, a single node, and negative pressure in the outer region, see Figs. 2c and 2d. But then at r 0,p = 2.394 fm, the pressure exhibits an additional (second) change of sign after which it remains positive. The asymptotic expression (35) for p(r) works with an accuracy of 1 % or better beyond r 3.4 fm. We stress that the classical model describes a ground state (and, in fact, no excited solutions exist in this model) [59]. But nevertheless p(r) exhibits two nodes.
Finally, the normal force 2 3 s(r) + p(r) is positive for 0 ≤ r < r 0,n in agreement with condition (iii), until exhibiting a node at r 0,n = 1.881 fm, after which it is negative, another new feature. Notice that 2 3 s(r) + p(r) ∝ 1 r 4 + . . . , implying that the mechanical mean square radius r 2 mech = d 3 r r 2 ( 2 3 s(r) + p(r))/ d 3 r ( 2 3 s(r) + p(r)) diverges. The 3 new features consist of (i) a node in s(r), (ii) a second node in p(r), and (iii) a node in the normal force. After the appearance of the nodes, the respective densities exhibit opposite signs as compared to prior studies. It is worth remarking that in Eq. (35), the asymptotics of T 00 (r) and p(r) are such that T 00 (r) − 3 p(r) = 0, which reflects the tracelessness of the EMT tensor T µ µ = 0 (in classical electrodynamics). The asymptotic expressions (35) for s(r) and p(r) satisfy the differential equation (30) which is dictated by the conservation of the EMT.
The size of s(r) and p(r) is very small in the regions where the new features occur. For instance, the second node of p(r) at r 0,p = 2.394 fm is beyond the range of Figs. 2c and 2d. However, had we tried to show it there, then the second node would be hardly visible on the scales of the Figs. 2c and 2d. In order to visualize the new features, we multiply the respective densities by r 4 such that the Coulomb contributions proportional to 1/r 4 appear as constant lines at large r, see Figs. 3b-3d. Despite the factor r 4 which enhances the densities at large r, the Coulomb contribution is small even in these plots. In particular, it is so small in the case of r 4 p(r) that an insert is necessary (with the scale on the y-axis enlarged by a factor of 10) to clearly show the second zero of p(r) in Fig. 3c.

B. The divergence of the D-term
The presence of the long-range electromagnetic forces also affects the D-term, which is an important particle property and on the same footing as the mass, spin or electric charge [8]. The D-term has two equivalent definitions in terms of shear force and pressure where n = 3 is the space dimension, which we leave here general for later purposes. These expressions are equivalent due to the EMT conservation, i.e. D = D s = D p gives the same result. This can be proven exploring the differential relation (30), cf. Ref. [16]. In the classical model of the proton, the D-term is undefined because the integrals in Eqs. (36,37) diverge linearly due to the asymptotic behavior (35) of p(r) and s(r) at large distances. Notice that the Figs. 3b and 3c basically show the integrands in Eqs. (36,37). Even though the Coulomb contribution is very small, it clearly spoils the convergence of the integrals in Eqs. (36) and (37). This is a new feature not encountered in previous studies. Typically in strongly interacting systems, the EMT densities decay at long distances fast enough such that the integrals defining the D-term converge. In quantum field theoretical models of the nucleon, the EMT densities exhibit an exponential fall-off at large r for finite pion masses. In the chiral limit, when the Goldstone boson (pion) becomes strictly massless, s(r) and p(r) behave like 1/r 6 at large r, which is still sufficient to guarantee a finite, well-defined D-term [16].
In the remainder of this work we will address the following questions. Does the model constitute a mechanically stable solution? Is it possible to obtain a prediction for the D-term in this model? And, are our observations model-dependent or of general character?

C. Mechanical stability in the model
Even though constructed as a consistent classical mechanical model of the proton [59], we find that the mechanical stability criterion (34) is not satisfied at large r. This issue needs to be resolved. Let us stress that positivity of the normal force (34) is a stability criterion of mechanical continuum systems [64]. Care may be needed when carrying over such criteria to quantum systems [8]. But here we deal with a classical continuum system and the criterion (34) must hold.
In this context, it is interesting to recall how the criterion (34) is used to determine the radius of a neutron star: density and radial pressure in the neutron star interior are governed by the Tolman-Oppenheimer-Volkoff equations which include general relativity effects, and are connected to each other by an equation of state of nuclear forces. The equation of state contains information on the compressibility of nuclear matter. The solution of these equations yields the normal force (which in neutron star literature is often referred to as "radial pressure" or simply "pressure", not to be confused with the pressure p(r) in this work.) The normal force is positive, but at some point it turns negative. This point marks the radius of the neutron star. Could we apply the same procedure to our case?
The answer is no. In our case, this procedure would mean to declare the radius r 0,n = 1.881 fm where the normal force exhibits a node to be the "edge" of the system. This "works" in the following sense. If we multiply (30) by r 3 , integrate over a finite integral 0 ≤ r ≤ R n , and perform integrations by parts, we obtain [8] r 0 dr r 2 p(r ) = r 3 3 This means that the von Laue condition (31) could also be satisfied by integrating over a finite interval from zero up to the node of 2 3 s(r) + p(r) at r 0,n = 1.881 fm. From a mechanical stability point of view, we could be happy about such a solution. From physical point of view, we are not. While the effects of the short-range strong fields are practically negligible beyond r 0,n = 1.881 fm one cannot ignore the effect of the long-range Coulomb field, which "communicates" the presence of an electric charge. We recall that the model [59] was designed with the specific purpose to have a mechanical model of an electric charge -which inevitably includes a correct description (within Maxwell's equations) of its long-range Coulomb potential. The "truncation" of the system at a finite value of r is therefore unacceptable (in our case; for the macroscopic neutron stars it surely works). We must seek a solution along different paths.
One possible resolution lies in exploring the force concept in a classical system. Notice that in our system, the matter ("continuous medium") is described by the scalar density ρ(r) localized within the radius R = 1.05 fm. Thus, in the volume where the dust particles are present, the criterion (34) is satisfied. The violation of (34) occurs where no matter is present. Hence, it does not affect the mechanical stability of the medium. This argument cannot be applied to quantum field theoretical systems where the distributions of "matter" and "field energy" cannot in general be distinguished [16]. Only in classical systems, such as our model, is such a distinction unambiguous.
Thus, a possible resolution of the issue is that no violation of (34) occurs in our model, because no matter is present at the point where the normal forces become negative. Hence, the node in the normal force causes no mechanical instability, and the classical proton solution is consistent.
Finally, we comment on the mean square radius of energy density r 2 E and mechanical mean square radius r 2 mech . These radii "measure" the extent of the spatial distributions of energy density and normal forces, which include the contributions of the fields, and diverge due to the long-range Coulomb field. In this classical model, the "proton size" is associated with the localized distribution of matter. The matter particles carry electric charge, and we already saw that system has a well-defined finite charge radius which numerically has the right order of magnitude, see Sec. II.

D. Regularized result Dreg for the D-term
In this section, we address the question of whether the model can make a prediction for the D-term. Taken literally, the expressions (36) and (37) for the D-term diverge, but one may try to regularize them. In general, regularization is not unique, and one must define a "regularization prescription". Let us stress we are talking here about the small contribution of the Coulomb potential in Figs. 3b and 3c, which spoil the convergence of the integrals.
In our case, one can use the following procedure to obtain a finite result for the D-term. In order to motivate this procedure, we notice that in a theory, where (a) the EMT is conserved and (b) the integrals defining D p and D s in Eqs. (36,37) exist, one can compute the D-term in terms of an arbitrary linear combination of D p and D s as follows Notice that under these conditions, the same result, D(ζ) = D, is obtained for any value ζ which follows simply from the equivalence of the expressions for D s and D p . In our case, the condition (a) is of course satisfied, but (b) is not. As a consequence, the expression D(ζ) is divergent for all ζ except for one value ζ = ζ reg , which can be chosen such that D reg = D(ζ reg ) is finite. This value of ζ reg depends on the number of dimensions n = 3 and the power N = 4 in the long-distance asymptotic s(r) = a s /r N and p(r) = a p /r N , and can be determined as follows. 1 The coefficients a s and a p are not independent, but related to each other, as a p /a s = −(n−1)(N −n)/(nN ) due to Eq. (30). Thus, in the linear combination, s(r)/a p + p(r)/a s , the long-range Coulomb tail cancels out. This is the only linear combination which can give a convergent result in Eq. (39). Considering the prefactors in the definitions (36,37) of D s and D p , the required value of ζ reg is 2 This is ζ reg = 8 3 in our case for n = 3 dimensions and N = 4. The regularized D-term obtained in this way is finite, negative, expressed in terms of p(r) and s(r) as follows, and numerically given by Numerically, this is about an order of magnitude smaller than the D-term in the chiral quark soliton model [16]. This is understandable considering the D-term encodes information on internal forces. In the classical model, we deal with "residual nuclear forces" which are about an order of magnitude smaller than the forces among quarks in the chiral quark soliton model of Ref. [16]. Our regularization method is distinguished because it removes the divergences from D p and D s in an efficient and "minimalistic" way: the long-range QED contribution exactly cancels out in the linear combination [6p(r) + s(r)] in (41) and we introduce no "regulator dependence." It it is also a suitable regularization because it preserves the negative sign of the D-term observed so far in all theoretical studies.
Another interesting feature of this regularization prescription is that it removes the Coulomb contribution not only in the outer region r > R, which is required to make the D-term finite. Interestingly, in the linear combination of p(r) and s(r) in the integrand in Eq. (41), the contribution of electrostatic forces also cancel out exactly in the inner region r < R. In other words, D reg receives no contribution from the electric forces at all. Notice that this concerns only the Coulomb contribution to D reg in this particular regularization scheme. The electromagnetic contribution to the budget of the internal forces is well-illustrated by the von Laue condition in Eq. (33), where the numerical contribution of the Coulomb field is small, but indispensable to prevent the collapse of the proton in the classical model.
However, our result (41) is not unique, because in principle one could use other methods to regularize the divergences. It would be interesting to see whether other suitable regularization methods can be defined, and investigate the effects of the regularization scheme dependence on the D-term. This will be left to future studies.

V. MODEL-INDEPENDENT INSIGHTS, AND THE FORM FACTOR D(t)
In this section, we put our findings in a wider context and show that the observed long-distance properties of EMT densities are model independent. We discuss possible implications for theoretical and experimental studies of the EMT form factors, and especially the D-term form factor.

A. EMT densities in QED at long-distances
Our results for EMT densities are certainly model-dependent in the region up to r 2-3 fm, where strong forces dominate. The strong forces are modelled in a specific way in our approach. One could use a different model for strong forces and would obtain different results in the region r 2-3 fm. Independently, of the chosen model, the strong forces are short-ranged and their effects are faded out in the region r 3 fm. This region is governed by the long-range electromagnetic force. 3 In other words, if one were to solve QCD and QED exactly, e.g. in a lattice calculation [82][83][84][85][86][87], one would recover the same results for the EMT densities at long distances r 3 fm. The reason for that is obvious. The 1 r -behavior of the classical Coulomb potential is a consequence of the masslessness of the photon in QED. Indeed, the long-distance part of the EMT densities in QED was determined from 1-loop calculations in [60] and is given by where the parenthesis denotes more strongly suppressed terms. To this order, the results in Eq. (42) are the same for spin-0 and spin-1 2 particles (while the components T 0k QED ( r) obviously depend on the spin) [60]. The result (42) can be traced back to non-analytic terms in the EMT form factors at small t which arise from the masslessness of the photon [60]. The QED long-distance contributions to the EMT densities (42) coincide exactly with our results (35). This is not a coincidence, but due to the fact that QED must reproduce the classical Maxwell theory at long distances [60].
As we have seen, the long-distance behavior (35,42) implies a divergent D-term as well as the divergence of mean square radii associated with EMT densities. These divergences are also model-independent results. This can be seen without involving the notion of EMT densities [60]. For that we have to investigate the EMT form factors.
B. The D-term form factor EMT densities can be computed directly in classical models but not in quantum field theory 4 where all one can do is to evaluate matrix elements of the EMT operatorT µν . The information content of these matrix elements is described in terms of form factors which are defined in the case of the nucleon as p , s |T µν |p, s =ū(p , s ) A(t) where P = 1 2 (p + p), ∆ = p − p, t = ∆ 2 , and the spinors are normalized asū(p, s) u(p, s) = 2M . The EMT densities can be inferred indirectly from the form factors through an interpretation of the 3-dimensional Fourier transforms of the form factors. This interpretation is justified if the size of the particle is much larger than its Compton wavelength (which is the case for the proton to a good approximation) and is applicable for r λ c where λ c = /(M c) ≈ 0.2 fm denotes the Compton wave-length of the proton [8]. 3 The point where electromagnetic force becomes dominant, in our work r 3 fm, is model dependent. In QCD, the long distance-behavior of nucleon EMT densities is dictated by spontaneous chiral symmetry breaking and the emergence of Goldstone bosons, pions in SU(2) flavor case, whose contributions to s(r) and p(r) decay like 1 r 4 exp(−2mπr) for pion mass mπ = 0, and are proportional to 1 r 6 in the chiral limit [16]. At large enough distances, the contributions of electromagentic forces dominate the EMT densities in any case. 4 One exception are quantum field theoretical models based on the limit of a large number of colors Nc where the nucleon structure is described in terms of a mean field [65], like in chiral quark soliton or Skyrme model [16]. In the large-Nc limit, the 3-dimensional EMT density formalism for baryons is exact [8]. Form factors also have an interpretation in terms of 2-dimensional EMT densities, which is exact for all hadrons and valid for any Nc [66,67], cf. Ref. [10] specifically for the case of EMT densities.
The interpretation of the form factors in terms of densities is performed in a frame where t = − ∆ 2 . If the form factor D(t) is known, one way to determine the stress tensor densities s(r) and p(r) is as follows [8] Here we can proceed "backwards" and integrate Eqs. (45,46) to obtain D(r). This is most conveniently done by integrating the expressions (45,46) twice over the radial distances from r to infinity, where all densities vanish. Starting from Eq. (45) or (46) respectively yields the same result for D(r) which serves as a test for the calculation.
With the result for D(r), one can invert the Fourier transform (44) to compute D(t). The result is shown in Fig. 4. The Fig. 4a gives an overview of D(t). In the region 0.1 (−t) 0.5 GeV 2 , we observe a shape typical for hadronic form factors, which can be approximated by a quadrupole form D(t) approx ≈ D reg /(1−t/m 2 D ) 3 , with m D = 0.985 GeV. At smaller (−t), the QED long-distance effects become noticeable. At significantly larger (−t), the correct description of form factors requires short-distance properties of QCD which are not present in the model. As the t = 0 intercept of D(t) approx , we have chosen the regularized value D reg from Eq. (39), which is indicated in Fig. 4a (one could also choose a slightly different value). The quadrupole shape of D(t) approx is suggested by large-t QCD counting rules [8] (but one could also choose other shapes). For t > 0.1 GeV 2 , the classical model result for D(t) is in good qualitative agreement with more realistic models [16] but about an order of magnitude smaller which is expected, cf. discussion below Eq. (41). For (−t) < 0.1 GeV 2 , the long-distance QED effects start to become important, and cause the form factor D(t) to change sign at t = −2.8 × 10 −4 GeV 2 . This "transition region" is shown in Fig. 4b. In this region, we can start to compare our results to the predictions for D(t) due to the QED leading non-analytic terms [60].

C. Comparison to the leading non-analytic QED contributions to D(t)
The leading non-analytic QED terms in the small-t behavior of EMT form factors were derived in [60] (where the notation q 2 = t, F 1 (q 2 ) = A(t), F 2 (q 2 ) = 2J(t), F 3 (q 2 ) = 1 4 D(t) was used). The derivation of the long-distance QED contribution to the stress tensor quoted earlier in Eq. effective field theory techniques. For a charged spin-1 2 fermion, the small-t behavior of D(t) due to QED-effects is [60] where the dots indicate terms which are finite as t → 0. Notice that D(t) is multiplied by (∆ µ ∆ ν − g µν ∆ 2 ) in (43). Therefore, the matrix element p , s |T µν |p, s has a well-defined limit t → 0, but the form factor D(t) does not. The D-term given by D = D(0) = lim t→0 D(t) is divergent. The result for the EMT densities (42) was obtained from Eq. (47) by means of a Fourier transform. The EMT determines the metric through the Einstein equation, and from the long-distance QED contribution to the EMT densities (42), it is possible to reproduce the classical non-linear terms of the Reissner-Nordström metric for a non-spinning charge or Kerr-Newman metric for a spinning charge [60].
It is important to stress that the metric in general relativity is an inherently classical concept. The deeper reason why it is possible to determine quantum corrections to the metric lies in the massless nature of the photon, which causes long-distance effects much stronger than the gravitational effects, provided α G M 2 /( c) = M 2 /M 2 Planck , where G denotes the gravitational constant and M Planck is the Planck mass. Under this condition, quantum gravity corrections can be neglected, and gravity can be treated as a classical theory described in terms of a metric (in our case quantum gravity effects can safely be neglected: the proton mass M is 19 orders of magnitude smaller than M Planck = 1.2 × 10 19 GeV). Besides the photon [60], graviton effects can also be studied in this way [68][69][70][71].
The QED result (47) is shown in Fig. 4b with the label "QED leading." After the sign change, the model result for D(t) starts to slowly approach the QED result (47). For "asymptotically small" t in the region below (−t) ≤ 10 −6 GeV 2 , the classical model result for D(t) practically coincides with the QED result (47), see Fig. 4c. In particular, the model result for D(t) diverges like 1/ √ −t for small (−t). For completeness, we remark that the form factors A(t) and J(t) must satisfy the constraints A(0) = 1 and J(0) = 1 2 [72][73][74] and do so of course despite the presence of long-range QED corrections. In the case of these form factors, the leading non-analytic terms are of the type √ −t and (−t) log(−t) and well-behaving for t → 0. But the derivatives of these form factors with respect to t diverge in the limit t → 0. For instance, in the case of A(t), this implies the divergence of the "gravitational mean square radius" [60] is defined as r 2 grav = −6A (0). The mean square radius of the energy density is also related to A (0) [8] and divergent due to the long-range QED effects, see above.

D. Consequences for calculations and measurements of D(t)
In theoretical studies of the hadron structure, electromagnetic effects can often be neglected to a good approximation. But in the experiment, one certainly cannot neglect the electric charge of the proton and other electromagnetic effects. Even though not straightforward [8], the form factor D(t) of the proton can be extracted from analyses of hard exclusive reactions like deeply virtual Compton scattering [2,3] using dispersion relation methods [75][76][77][78] and first attempts were reported [79,80]. Can the divergent behavior of D(t) due to QED effects (47) be experimentally observed?
It is important to stress that the QED contribution to D(t) in Eq. (47) starts to become noticeable in our model only for (−t) 0.1 GeV 2 . In more realistic models, the contribution of strong forces to D(t) is an order of magnitude larger, implying that the transition region where D(t) changes sign, cf. Fig. 4b, sets in at even lower values of (−t). In addition, in the case of deeply virtual Compton scattering, it is necessary to consider higher order QCD corrections. When extracting information from electromagnetic processes, one must also consider QED radiative corrections, which are generically of the same order of magnitude as the effect (47). It is by no means obvious whether the result (47) can be disentangled from radiative corrections. Even if it can, it remains to be seen whether such corrections can be determined with sufficient precision to observe (47). Thus, from a practical point of view, one may never be able to reach the region of small enough (−t), cf. Figs. 4b and 4c, and sufficient precision to observe QED effects like (47).
But from a theoretical point of view, it is legitimate to ask how to calculate the D-term in a system with long-range forces. We do not know a definite answer to this question, though our work indicates a possible solution, namely to apply a regularization scheme. At first glance, it may appear unusual to invoke regularization in classical calculations. But the deeper reason why the D-term diverges is rooted in the masslessness of the photon, and hence related to infrared divergences in QED. The regularization prescription proposed in Sec. IV D is acceptable because: (i) it gives a finite result, (ii) preserves the negative value of the D-term in accordance with other theoretical studies, and (iii) the obtained numerical value (39) is useful to practically approximate D(t) at finite (−t) > 0.1 GeV 2 , see Fig. 4a. Thus, working with a such a regularization scheme is one practical way out. The regularization method of Sec. IV D works in our classical calculation. In perturbative QCD and QED calculations of the deeply virtual Compton scattering process, of course other "schemes" are invoked, and in nonperturbative lattice QCD calculations with included QED effects, one can use yet other regularization methods [82][83][84]. More theoretical work is required to compare results obtained in different schemes. These aspects go beyond the scope of this work.

VI. CONCLUSIONS
Prior EMT studies focused on applications to hadronic physics, and considered mainly strongly interacting systems with short-range forces. Long-range forces were not included. In systems governed by short-range forces, the D-term was always found to be well-defined, finite and negative. In this work, we have presented a study in a system where in addition to (strong) short-range forces, (electromagnetic) long-range forces are also present. We have encountered several interesting features not observed in prior studies of systems with short-range forces. The most interesting observation is that, when long-range forces are included, the D-term is no longer well-defined, and diverges.
For our study, we employed the classical proton model of Bia lynicki-Birula [59] which is of interest for its own sake. To the best of our knowledge, it is the first fully consistent classical model of an extended charged particle where the Poincaré stresses are generated dynamically in a local, relativistic, classical field theory. The model exhibits short-range strong forces, which are modelled after nuclear forces, and the electromagnetic long-range interaction [59]. The two crucial aspects for our study are the consistent description of a stable particle, and correct description of the long-range electromagnetic effects. The model of Ref. [59] satisfies both. The classical aspect of the model is an advantage in the sense that it allows us to concentrate on the effects of long-range forces undistracted by technical difficulties, which are inevitable in studies of more realistic, strongly interacting quantum systems.
In the region below r 2 fm, the classical proton model yields results for the energy density T 00 (r), shear force s(r), and pressure p(r) in good qualitative agreement with more realistic models like the chiral quark soliton or Skyrme model, except that s(r) and p(r) are about an order of magnitude smaller. This is of course to be expected in a model where the internal forces are modelled after the "residual nuclear forces." Otherwise, in this r-region, the classical proton model is in line with results from short-range systems.
The situation is different for r 2-3 fm, where the strong forces are faded out, and T 00 (r), s(r), p(r) exhibit tails proportional to 1 r 4 due to the long-range Coulomb field. This introduces new features, e.g. s(r) and p(r) show (additional) nodes and opposite signs at large-r as compared to systems with short-range forces. Another consequence is the divergence of the mean square radius of the energy density and the mechanical radius. These radii measure the spatial extensions of the energy density T 00 (r) and normal force 2 3 s(r) + p(r), which include contributions from fields. Due to the long-range of the Coulomb field, the size of the system, as measured by these radii, is consequently infinite. In contrast, the electric mean square radius is finite (and numerically of the right size) [59], as it measures the spatial extension of the electric charge distribution tight to the localized matter distribution in the model.
The most interesting new feature is that the D-term of a charged particle is divergent. Technically, this happens because the D-term is given in terms of integrals over s(r) or p(r) which diverge due to the 1 r 4 -tails of these densities. We proposed a regularization scheme which yields a finite, negative, and numerically reasonable value for a system with internal interactions of the strength of "residual nuclear forces." We computed the form factor D(t) which, for (−t) > 0.1 GeV 2 , is negative and shows a shape typical for hadronic form factors. But in the region (−t) 0.1 GeV 2 , the form factor changes sign, becomes positive and diverges as t → 0.
The observed long-distance properties of EMT densities and the related small-t divergence of D(t) are modelindepedent results. Both have been derived in Ref. [60] from QED diagrams. In a recent study of Q-balls carrying an electric charge, the same long-range tails were found as in our work [88]. The deeper reason for the emergence of these EMT properties can be traced back to the masslessness of the photon [60], which is reflected in the classical Maxwell's equations. Consequently, the EMT long-distance properties must be correctly reproduced in every system (classical, quantum mechanical, quantum field theoretical) where the electromagnetic interaction is correctly described.
Other EMT form factors are also affected by QED long-distance effects, but less strongly than D(t) for two reasons. First, the proton EMT form factors A(t) and J(t) are constraint to satisfy A(0) = 1 and J(0) = 1 2 [72][73][74]. QED longdistance effects must preserve these constraints, though they can (and do) affect the derivatives of these form factors (making them infinitely steep at t = 0 which in turn is connected to the divergence of the related mean square radii). The value of D(t) at t = 0 is in general not constrained by any principle, and can therefore show more variation than other form factors. Second, the D-term is intrinsically related to the internal forces and the dynamics in a system. As such, it is the particle property which exhibits by far the strongest sensitivity to variations in the system. It is therefore not unexpected that D(t) shows the most pronounced effects when long-range forces are included.
The long-distance QED effects on D(t) become noticeable at such small (−t) 0.1 GeV 2 that it is not clear whether they are, even in principle, measurable. More theoretical work is required to clarify this important point. Our results are very interesting from a theoretical point of view, and raise the question of how to define the D-term in a system with long-range forces. Another well-known long-range force is gravity. It would be very interesting to perform a fully consistent computation of EMT properties including general relativity. Solutions of the Einstein equation for a perfect charged fluids exist, see e.g. [81], but they typically make assumptions about the density or the equation of state. A consistent treatment of the D-term requires to exactly solve the dynamics of all involved fields, including the gravitational field. As noted in [46], due to its sensitivity to the details of the involved interactions, the correct definition of the D-term may require the consideration of all forces, including QED and perhaps even gravity.