Revisiting the mechanical properties of the nucleon

We discuss in detail the distributions of energy, radial pressure and tangential pressure inside the nucleon. In particular, this discussion is carried on in both the instant form and the front form of dynamics. Moreover we show for the first time how these mechanical concepts can be defined when the average nucleon momentum does not vanish. We express the conditions of hydrostatic equilibrium and stability in terms of these two and three-dimensional energy and pressure distributions. We briefly discuss the phenomenological relevance of our findings with a simple yet realistic model. In the light of this exhaustive mechanical description of the nucleon, we also present several possible connections between hadronic physics and compact stars, like e.g. the study of the equation of state for matter under extreme conditions and stability constraints.

The paper is organized as follows. In Section II we explain how to describe a relativistic quantum system localized in phase space and present the matrix elements of the general asymmetric EMT. Mechanical properties of the nucleon defined in the instant form of dynamics are discussed in Section III. After a quick review of the proper decomposition of the nucleon EMT into quark and gluon contributions, we present the three-dimensional spatial distributions defined in the Breit frame where the system is in average at rest. We extend the study to the more general class of elastic frames introduced in [29], where two-dimensional spatial distributions can be defined in the plane transverse to the average motion of the system. In Section IV, we discuss for the first time the mechanical structure of the nucleon using the front form of dynamics. In such a framework, we define two-dimensional spatial distributions free of relativistic corrections and compare them with the corresponding distributions in instant form. The questions of hydrostatic equilibrium and stability conditions are discussed in Section V. Finally, we conclude the paper with Section VI where we summarize our results. In order to define the distribution of a physical quantity inside a system, one should first localize the system in both position and momentum space. In a quantum theory, this can be achieved in the Wigner sense [30]. A system with average position R and average momentum P is described by the covariant phase-space density operator ρ R,P = dP 2 2π where the delta functions 1 ensure that initial and final states have the same mass M . It follows from the relativistic normalization for momentum eigenstates p |p = 2p 0 (2π) 3 δ (3) (p − p) that Tr[ρ R,P ] = 1. Defining the covariant "position" states as with normalization x |x = , the covariant phase-space density operator can alternatively be expressed as If we integrate the covariant phase-space density operator over the position R, we recover the density operator in momentum space and if we integrate over the momentum P , we recover the density operator in position space d 3 P (2π) 3 2P 0 ρ R,P = |R R| .
and translation invariance implies that where x = X − R is the relative average position.
In the present work, we are interested in matrix elements of the (renormalized) nucleon EMT T µν (X). The latter can be expressed in terms of gravitational (or energy-momentum) form factors (GFFs), first introduced 2 by Kobzarev and Okun [32] and later by Pagels [33]. In the general case of a local, gauge-invariant asymmetric EMT for a spin-1 2 target, the standard parametrization reads [28,29,34,35] p , s |T µν a (0)|p, s =ū(p , s ) where p (p ) and s (s ) are the four-momentum and canonical polarization of the initial (final) nucleon of a mass M , η µν = diag(+, −, −, −) is the Minkowski metric, and t = ∆ 2 with ∆ = p − p and P = (p + p)/2. For convenience, we introduced the symmetrizer a {µ b ν} = a µ b ν + a ν b µ , the antisymmetrizer a [µ b ν] = a µ b ν − a ν b µ , and the notation iσ µ∆ = iσ µλ ∆ λ . The label a in Eq. (8) refers to the contribution from a particular type of constituents, typically a = q for quarks and a = G for gluons. The total EMT is then simply obtained by summing over all the constituent types T µν = a T µν a . The generic matrix element (8) is parametrized in terms of five GFFs, namely A a (t), B a (t), C a (t),C a (t) and D a (t). Beside their t-dependence, the GFFs are also usually scale and scheme-dependent. Except for D a (t), these additional dependences drop out when summing over all quark and gluon contributions. In particular, three of the GFFs associated with the symmetric part of the EMT satisfy the sum rules which arise from the Poincaré invariance of the theory [28,36,37]. In particular, the vanishing of the total anomalous gravitomagnetic moment a B a (0) = 0 is related to the equivalence principle in General Relativity [32,36,38] and holds separately for each Fock component of the state [39]. The GFFC a (t) accounts for the non-conservation of the partial EMT p , s |∂ µ T µν a (0)|p, s = i∆ ν Mū(p , s )u(p, s)C a (t) and should naturally vanish when summed over all the constituents. The GFFs have been studied in a variety of theoretical approaches, see e.g. [40] and references therein, and also using Lattice QCD simulations [4][5][6][7][41][42][43].
Although a direct measurement of nucleon scattering by a gravitational field is not realistic with the current technology, it is remarkable that the GFFs can in principle be extracted from experimental data. Ji [1,44] showed that the three GFFs A a (t), B a (t) and C a (t) can be obtained from the second Mellin moment of leading-twist generalized parton distributions (GPDs) [2,45,46], which are accessible in several exclusive processes, like deeply virtual Compton scattering [47] and meson production [48]. Recently, the corresponding GFFs for a pion target have been extracted from Belle data on γ γ → π 0 π 0 [49]. The GFFC a (t), which can formally be obtained from higher-twist GPDs [28,50,51], is related to the σ πN term extracted from πN scattering amplitudes [52,53], and to the trace anomaly which can be studied through the production of heavy quarkonia at threshold [54][55][56][57][58]. Finally, the GFF associated with the antisymmetric part of the EMT is directly related to the axial-vector form factor [29] and hence can be obtained from quasi-elastic neutrino scattering and pion electroproduction processes [59]. For illustrative purposes, we will adopt in this work a simple multipole Ansatz for the GFFs which is supported by model calculations for |t| < 1 GeV 2 [26], together with the parameters given in Table I. We adopt a standard dipole Ansatz (i.e. n F = 2) for A a ,C a and D a , but for B a and C a we choose a tripole Ansatz (i.e. n F = 3) in order for the energy and pressure distributions to be realistic, see discussion in Section V. The normalization A q (0) ≈ 0 .55 is taken from the recent MMHT2014 analysis [60]. We set B q (0) ≈ −0 .07 as suggested by the AdS/QCD correspondence [61,62] and which agrees in magnitude with recent estimates from Lattice QCD [7,63]. We also use the values C q (0) = d q 1 (0)/5 with d q 1 (0) ≈ −1.59 obtained in a dispersive analysis of deeply virtual Compton scattering [64] which is close to a recent experimental extraction reported in [65],C q (0) ≈ −0 .11 given by a phenomenological estimate [66] supported by a recent Lattice calculation [6], and D q (0) = −G q A (0) ≈ −0 .33 obtained from a leading twist NNLO analysis by the HERMES collaboration [67]. The quark multipole masses Λ Fq with F = A, B, C,C are motivated 3 by results obtained in the chiral quark-soliton model [26,68], and Λ Dq = Λ G q A is taken from a recent Lattice estimate [69]. In the gluon sector, the normalizations A G (0), B G (0) andC G (0) are determined by the sum rules (9). As suggested in [58], we use the simple relation C G (0) = 16 3n f C q (0) with n f = 3. Since we lack information about the gluon GFFs, we simply set Λ F G = Λ Fq for F = A, B, C,C.
This simple parametrization should be considered only as a naive model with the sole aim of allowing us to illustrate the various distributions discussed in the rest of the paper. In particular, one of its advantages is to permit an evaluation of the two and three-dimensional Fourier transforms of the multipole distributions in closed forms where b = |b ⊥ | with b ⊥ a two-dimensional vector, and r = |r| with r a three-dimensional vector.

III. DISTRIBUTIONS IN INSTANT FORM
Performing the integrals over P 2 and ∆ 0 in the covariant phase-space density operator (1) leads to where the initial and final target energies are given by The non-explicitly covariant form (15) coincides with that of Appendix B in Ref. [70] if we replace the normalization factor 2P 0 by 2 p 0 p 0 . The difference in the normalization comes from the fact that "position" states in Ref. [70] were defined with the non-relativistic normalization x |x = δ (3) (x − x) at equal times. This difference will however not concern us since we will essentially be interested in the case ∆ 0 = 0, when p 0 = p 0 . Setting the origin at the average position of the system R = 0 and denoting x µ = (0, r), the static EMT encoding the distribution of energy and momentum inside the system with canonical polarization s and average momentum P , is defined as the following Fourier transform of the off-forward amplitude [29] T µν Note that x 0 = 0 because the positions X and R are considered at the same time X 0 = R 0 . If we allow X 0 and R 0 to be different, the above static EMT T µν a (r; P ) can be recovered by considering the time average dX 0 /(2π δ(0)) in a frame where the energy transfer vanishes ∆ 0 = 0 [10].
The GFFs in Eq. (8) are multiplied by two types of Dirac bilinears, namelyū u andū iσ µ∆ u. Using the same canonical polarization for both initial and final states, these bilinears can be expressed in instant form as [71] where S µ = (0 , s) and N = p 0 + M p 0 + M . In the present work, we will restrict ourselves to the case of an unpolarized target, which amounts to setting S µ = 0 in the above expressions. The unpolarized off-forward amplitude then reads T µν a (0) = Except the global normalization factor, the first line is identical to the standard parametrization for a spin-0 state. The second and third lines can be interpreted as the polarization-independent distortion arising due to the spin of the target. Indeed, it has been shown in [10,29] that the combinations [A a (t) + B a (t)]/2 and −D a (t)/2 provide the information about the spatial distribution of total angular momentum and spin associated with constituent type a. The static EMT can receive a (quasi-)probabilistic interpretation only when no energy is transferred to the system ∆ 0 = 0 [29]. Since the onshell conditions impose that ∆ 0 = P · ∆/P 0 , we will consider the following cases: (a) ∆ = 0 -forward limit (FL); (c) P · ∆ = 0 -elastic frame (EF); (d) P · ∆ = 0 and |P | → ∞ -infinite-momentum frame (IMF).
Note that when ∆ 0 = 0 the normalization factor appearing in Eq. (21) reduces to N = P 0 + M .

A. Forward limit
Since ∆ is the Fourier conjugate variable to relative position r, the forward limit (FL) ∆ = 0 is obtained by integrating the static EMT over r d 3 r T µν a (r; P ) = P, s|T µν a (0)|P, s 2E P where E P = P 2 + M 2 . Note that the dependence on the nucleon spin disappears in the FL. This is expected because the EMT is the Noether current associated with invariance under translations and not Lorentz transformations. Focusing on the T 00 a component in the rest frame P = 0, Ji [72,73] proposed a decomposition of the nucleon mass based on Eq. (22). Recently, a covariant treatment of T µν a revealed that the gravitational charges A a (0) andC a (0) can be interpreted in terms of partial internal energy density and isotropic pressure [66]. Indeed, denoting the proper volume by V and the boost factor by γ = E P /M , one finds that the average density has the same structure as the EMT of an element of perfect fluid [74] θ µν (r) = (ε +p) u µ u ν −p η µν .
The four-velocity of the nucleon being given in the FL by u µ = P µ /M , this suggests that the following combinations can be interpreted 4 as the spatial average of partial energy density and isotropic pressure associated with constituent type a. The contributions to proper internal energy and pressure-volume work are then given by The nucleon being a stable system with mass M , one obtains a mass sum rule and a stability constraint consistent with Eq. (9). While A q (0) is well determined [60],C q (0) is poorly known. Based on the phenomenological estimates in [75],C q (0) seems to be negative and sizeable [66,76], in agreement with the MIT Bag Model prediction [77] and recent Lattice estimates [6]. In contrast, it has been suggested based on the instanton picture of the QCD vacuum thatC q (0) at low scale may be small and positive [68].

B. Breit frame
Since the work of Sachs on the electromagnetic form factors [78], the Breit frame (BF) defined by P = 0 became a popular frame for the physical analysis of form factors in instant form. The phase-space perspective adopted in the present work shows that working in the Breit frame amounts to looking at the system which is in average at rest and sitting in average at the origin.
The unpolarized off-forward amplitude (21) reduces in the BF to where we used P µ = η µ0 P 0 . Interestingly, the GFF D a (t) does not contribute in the BF. We therefore recover the case of a symmetric EMT studied by Polyakov et al. [10,11,26]. After Fourier transform, we find the following unpolarized static EMT where x µ = (0, r) and the three-dimensional Fourier transforms of GFFs are denoted by with t = −∆ 2 . We observe that the unpolarized static EMT in the BF (29) has the same structure as the EMT of an anisotropic spherically symmetric compact star [79] Θ µν (r) = [ε(r) where u µ and χ µ = x µ /r are unit timelike and spacelike four-vectors orthogonal to each other. The functions ε(r), p r (r) and p t (r) represent the energy density, radial pressure and tangential pressure, respectively. As noticed by The radial pressure pr(r) and the tangential pressure pt(r) at a distance r from the center of the system (31). Spherical symmetry imposes only the equality of the two tangential pressures.
Einstein and developed first by Lemaitre in 1933 [80], spherical symmetry requires only the equality of the two tangential pressures, see Fig 1. The tensor (31) can alternatively be written as with h µν = u µ u ν − η µν . Isotropic pressure p(r) and pressure anisotropy s(r) are related to radial and tangential pressures as follows The comparison of the unpolarized static EMT in the BF (29) with the EMT of an anisotropic spherically symmetric compact star (31) or (32) with u µ = η µ0 suggests that the following combinations ε a (r) = M A a (r) +C a (r) can be interpreted as the partial energy density, radial pressure, tangential pressure, isotropic pressure, and pressure anisotropy associated with constituent type a, respectively. They can alternatively be written as  Table I, see Eq. (45) for the definition.
As indicated by the presence of B a (t), the non-zero spin of the target affects only the energy distribution in the BF. Classically, we indeed expect angular momentum to push matter away from the center. The above distributions are illustrated in Figs. 2-7 in units of GeV/fm 3 = 1.7827 × 10 15 g/cm 3 using the multipole model (11) with parameters given in Table I. The energy density in Fig. 2 is always positive and is approximately shared equally between quark and gluon contributions. One defines the corresponding average squared mass radius as In our simple model, we find R M = 0 .905 fm which is a bit larger than the charge radius R Q = 0 .841 fm extracted from muonic hydrogen spectroscopy [81,82] and R Q = 0 .879 fm extracted from electron-proton scattering [83].
Knowing the distribution of energy density, it is also easy to derive the standard mass function widely used in General Relativity which represents the mass contained within a sphere of radius r, see Fig. 3. While the total radial pressure in Fig. 4 is always positive and largely dominated by the quark contribution, the total tangential and isotropic pressures in Figs. 5 and 6 switch from positive sign at the center of the nucleon (where it is dominated by the quark contribution) to negative sign at the periphery (where it is dominated by the gluon contribution). The pressure anisotropy in Fig. 7 vanishes at the center of the nucleon, as required by spherical symmetry, and is positive anywhere else, indicating that the radial pressure is always larger than the tangential one. Looking at the separate contributions, we see that the quark and gluon radial forces are both repulsive and of similar range. For the tangential forces, the quark contribution appears to be mostly repulsive and short range whereas the gluon contribution appears to be mostly attractive and long range. If we integrate the energy density and the isotropic pressure over the whole volume, we naturally recover the FL (26) discussed in the former section One can also relate the value of the GFF C a (t) at t = 0 to a weighted integral of the pressure anisotropy (43) [10,11,26] Summing over the constituents, one obtains the following additional relations [26,84] Table I, see Eq. (38) or Eq. (43) for definitions in terms of GFFs.
Several hints suggest that C(0) is likely negative [11,26,85,86]. Hudson and Schweitzer [85] observed that while adding total derivatives to the EMT leaves the Poincaré generators unaffected (and hence the particle mass and spin), it does change C(0). Since this term can be extracted from experimental data, one should not be allowed to add these divergence terms without changing some scheme prescriptions, contrary to the common belief. The same conclusion is reached when one consistently treats intrinsic angular momentum at the level of spatial distributions [29,70]. Interestingly, in view of the energy density and pressure conditions encountered in a nucleon, one may conjecture that studies of the nucleon EMT will shed some light on the EoS inside compact stars [66], which so far remains largely unknown [87], and will therefore complement efforts based on heavy-ion collisions [88] and gravitational wave observations [89][90][91][92][93][94]. Eliminating the radial variable r in Eqs. (34)-(37), we obtain the nucleon EoS for radial pressure p r (ε), tangential pressure p t (ε), and isotropic pressure p(ε). The results plotted in Figs. 8 and 9 show a pretty stiff behavior compatible with the observation of supermassive (∼ 2 M ) compact stars [90,95,96] well above the Chandrasekhar mass limit 1.44 M [97]. Although our multipole model is very naive, it supports the idea of an exciting crosstalk between hadronic physics and compact stars. An example of such a connection is given by the use of quark models to study the EoS of potential strange quark matter inside compact stars [14,15,[98][99][100][101][102][103].

C. Elastic frame
Spatial distributions with quasi-probabilistic interpretation can also be introduced when P = 0 [29]. In order to maintain the condition ∆ 0 = 0, we have to restrict ourselves to the set of elastic frames (EF) defined by P · ∆ = 0. They can be obtained by integrating the static EMT over the longitudinal coordinate r = r · P /|P | where ∆ ⊥ and b ⊥ = r ⊥ are vectors lying in the two-dimensional plane orthogonal to P , see Fig. 10. The unpolarized off-forward amplitude (18) in the EF takes the form If we further integrate the static EMT over the impact-parameter b ⊥ , which amounts to setting ∆ ⊥ = 0 ⊥ in Eq. (50), we recover the FL expression (22). If we set P = 0 in Eq. (50), we recover the BF expression (28) integrated over r , i.e. with ∆ = 0.
Since P 0 depends on ∆ ⊥ in the EF, we could not find simple expressions for the spatial distributions T µν a (b ⊥ ; P ) in terms of Fourier transforms of GFFs. Moreover, these spatial distributions will be |P |-dependent. Let us choose for convenience the z-axis along P . If we restrict ourselves to the (1+2)-dimensional (or transverse) static EMT 5 with components α, β ∈ {0, 1, 2}, its structure looks the same 6 as that of an anisotropic axially symmetric compact star in two dimensions where v α = (1, 0 ⊥ ) and χ α = (0 , b ⊥ /b). The functions ρ(b, P z ), σ r (b, P z ) and σ t (b, P z ) represent the twodimensional version of energy density, radial pressure and tangential pressure, respectively. Note that we have included explicitly a factor of γ 2 in the energy density component that comes from the longitudinal Lorentz boost allowing us to compare directly ρ(b, P z ) for different values of P z . In particular, for the total EMT we have d 2 b ⊥ ρ(b, P z ) = M . The tensor t αβ can alternatively be written as The two-dimensional isotropic pressure σ(b, P z ) and pressure anisotropy Π(b, P z ) are related to radial and tangential pressures as follows The particular case P z = 0 is simple and corresponds to the BF. We find 6 Note that the proper identification t αβ ∼ γ T αβ involves a boost factor accounting for the Lorentz contraction of the volume just like in the FL, see Eq. (23).
where x α = (0, b ⊥ ) and the two-dimensional Fourier transforms of GFFs are denoted by with t = −∆ 2 ⊥ . Comparing Eq. (54) with the EMT of an anisotropic axially symmetric compact star in two dimensions (51) and (52) suggests that the following combinations can be interpreted as the two-dimensional partial energy density, radial pressure, tangential pressure, isotropic pressure, and pressure anisotropy associated with constituent type a. They can alternatively be written as Integrating Eq. (29) over z allows us to alternatively express these quantities in terms of the three-dimensional energy density and pressures as follows with r = √ b 2 + z 2 . These distributions are illustrated in Figs. 11-15 in units of GeV/fm 2 = 178.27 g/cm 2 using the multipole model (11) with parameters given in Table I. Their behavior turns out to be similar to the corresponding three-dimensional distributions, see Figs. 2-7, except for the gluon contribution to the radial pressure, compere Figs. 4 and 12. This is an effect of the projection onto the transverse plane which mixes three-dimensional radial and tangential pressures together, as indicated by Eq. (67). While the latter have the same sign for quarks, they are of opposite sign for gluons.  ρ(b, 0), using the multipole model (11) with parameters given in Table I σt(b, 0), using the multipole model (11) with parameters given in Table I, see Eq. (58) or (63) for the definition in terms of GFFs.  σ(b, 0), using the multipole model (11) with parameters given in Table I, see Eq. (59) or (64) for the definition in terms of GFFs.  Π(b, 0), using the multipole model (11) with parameters given in Table I Similarly to the three-dimensional case, we can relate the value of the GFF C a (t) at t = 0 to a weighted integral of the pressure anisotropy Summing over the constituents, we also find the additional relations which are simply the two-dimensional version of those appearing in Eq. (48).
The last term in Eq. (50) makes the EMT asymmetric. In particular, the density of longitudinal momentum T 03 a is not equal to the longitudinal flux of energy T 30 a when the GFF D a (t) is not identically zero. Using Poincaré invariance, the antisymmetric part of the EMT T [µν] can be expressed in terms of the intrinsic angular momentum current S λµν as follows [28,29] In the case of QCD, we get for the corresponding off-forward matrix elements where the matrix elements of the axial-vector current are parametrized as follows in terms of the axial-vector and induced pseudoscalar form factors. According to Ref. [71], the corresponding bilinears can be expressed in instant form as For an unpolarized target S µ = 0, we can write and hence we find in the EF Using now the expression in Eq. (50) for the left-hand side, we recover the relation D q (t) = −G q A (t). We are now ready to discuss the distribution of spin inside an unpolarized target. The quark spin operator being given by 1 2 ψ(0)γ i γ 5 ψ(0), the corresponding distribution in the EF is given by We see that even in an unpolarized target, a nonzero quark spin distribution appears when the target is moving, see Fig. 16. The spin direction is orthogonal to both the target momentum and the impact parameter. This is reminiscent of transverse shifts observed for transversely polarized moving target, see [70] and references therein. In the latter case, the transverse shifts appear because of a nonzero net orbital angular momentum in the system. In the present case, the appearance of transverse spin distribution is a result of spin-orbit coupling [104][105][106].

D. Infinite-momentum frame
The infinite-momentum frame (IMF) is a special case of the EF obtained by considering the limit P z → ∞ [107]. We find that the (1+2)-dimensional unpolarized off-forward amplitude (50) reduces in the IMF to After two-dimensional Fourier transform, we obtain for where x α = (0, b ⊥ ) and therefore, by comparison with Eqs. (51) and (52), we can write keeping only the leading terms in P 2 z P 2 and taking the boost factor γ = E P /M into account in the identification t αβ ∼ γ T αβ .
While the two-dimensional pressures are the same in both the BF (P z = 0) and IMF (P z → ∞), the energy densities (56) and (83) differ. This may be interpreted by the expectation that kinetic energy grows with P z whereas binding energy associated with pressure forces remains constant. In the IMF, kinetic energy becomes by far the dominant contribution and we recover the parton picture where quarks and gluons behave as almost free massless particles. Since the two-dimensional pressures in the IMF coincides with those shown in Figs. 12-15, we simply illustrate in Fig. 17 the energy density in the IMF, using the multipole model (11) with parameters given in Table I. The comparison with the energy density in the BF shows that the kinetic energy is more concentrated around the center of the nucleon than the binding energy. Naturally, once summed over all the constituents and integrated over the impact parameter space, the difference disappears.

IV. DISTRIBUTIONS IN FRONT FORM
The interpretation of form factors in the Breit frame are known to be plagued by relativistic corrections [10,108,109]. In particular, one could multiply the off-forward amplitude (18) by some function f (P 0 /E P ) normalized as f (1) = 1 to correct for Lorentz contractions effects. While such correction factor does not change the integrated quantities, it does change the spatial distributions since it introduces an additional ∆-dependence. Moreover, the correction factor cannot be determined in practice in a model-independent way because Lorentz boosts depend on the dynamics of the interacting system.
An interpretation free of such relativistic corrections can however be obtained using the light-front (LF) formalism [110,111], which amounts to adopting the point of view of a massless observer [70]. This remarkable feature is explained by the fact that, in the LF formalism, the subgroup of Lorentz transformations associated with the transverse plane is Galilean [27,112]. Burkardt used this formalism and introduced the boost-invariant impact-parameter distributions (IPDs) of quarks and gluons [27,108]. Including parton transverse momentum to the picture led then to the concept of relativistic phase-space (or Wigner) distributions [104,[113][114][115]. Longitudinal LF momentum IPDs [49,61,62,[116][117][118][119][120][121][122] and longitudinal angular momentum IPDs [29,62,123] have also recently been discussed in the literature. Our aim here is to introduce the IPDs associated with the EMT.  Table I, see Eqs. (56) or (61), and (83) for the definition in terms of GFFs.

A. Light-front components
The LF components of a four-vector are given by where a ± = (a 0 ± a 3 )/ √ 2. In terms of these, the scalar product of two four-vectors reads One conventionally chooses x + to the represent the LF time coordinate. It then follows from Eq. (89) that p − represents the LF energy. The other conjugate components x − and p + represent the longitudinal LF coordinate and momentum, respectively. Without loss of generality, we choose the z-axis along P , so that we simply have We will denote the spatial LF three-vectors asã and use a dotted notation to keep the tensor expressions compact. A LF three-vector with a dotted index will indicate a minus first component followed by the transverse ones, while an undotted index will indicate a plus first component followed by the transverse ones, i.e. aα = (a − , a ⊥ ) and a α = (a + , a ⊥ ). Thus, in position space, the spatial LF coordinates read xα = (x − , x ⊥ ), and in momentum space they read p α = (p + , p ⊥ ). We can then rewrite the scalar product of two four-vectors (89) as withx ·p = xα p α = xα p α = x − p + − x ⊥ · p ⊥ . Using LF components is convenient because they behave in a simple way under LF boosts [110,111]. In particular, performing a longitudinal LF boost amounts to a mere rescaling of the LF components [a + , a − , a ⊥ ] → [e −ω a + , e ω a − , a ⊥ ], where γ = cosh ω.
Because of the Galilean symmetry in the transverse plane, it is interesting to organize the LF components of the EMT as follows [124] The upper left corner corresponds to a (1+2)-dimensional Galilean EMT T αβ . The upper right corner corresponds to a "mass" current J α = T α+ , where the role of "mass" in the transverse plane is played by the longitudinal LF momentum. Since we are considering only the static EMT, the above currents are conserved in the (1+2)-dimensional subspace Together, they form the so-called covariant non-relativistic stress-energy tensor T αµ appearing in the context of Newton-Cartan geometries, which find important applications in condensed matter and in the study of non-relativistic holographic systems, see e.g. [125] and references therein. The last line in Eq. (92), which describes the flux of energy and momentum along the spatial LF direction x − , does not have any known simple interpretation within the Galilean picture.

B. Light-front amplitudes
We can repeat the same procedure as in Section III in the LF formalism. Integrating the covariant phase-space density operator (1) over P 2 and ∆ − leads to with ∆ − = −∆ + P − /P + , since P ⊥ = 0 ⊥ (90), and The non-explicitly covariant form (94) coincides with that in Ref. [29] if we replace the normalization factor 2P + by 2 p + p + . Once again, the difference in the normalization comes from the fact that "position" states in Ref. [29] were defined with the non-relativistic normalization x |x = δ (3) (x −x) at equal LF times. This difference will however not concern us since we will essentially be interested in the case ∆ − = p − − p − = 0. Setting the origin at the average position of the system, the LF distribution of energy and momentum inside the system with canonical polarization s and average LF momentumP = (P + , 0 ⊥ ) is given by the Fourier transform of the LF off-forward amplitude [29] T µν a (0) = p , s|T µν a (0)|p, s 2P + .
Note that x + = 0 because the LF positions X − and R − are considered at the same LF time X + = R + . If we allow X + and R + to be different, the above static EMT T µν a (x;P ) can be recovered by considering the LF time average The Dirac bilinears appearing in the parametrization of the EMT (8) can be expressed in front form as [71] where S µ ⊥ = [0 , 0 , s ⊥ ] and N = p + p + . The unpolarized off-forward LF amplitude then reads Note that the structure is slightly simpler than the corresponding expression in instant form (21). The static LF EMT can receive a (quasi-)probabilistic interpretation only when no LF energy is transferred to the system ∆ − = 0 [29]. Since the onshell conditions impose that ∆ − = −∆ + P − /P + and since P ± > 0 for a massive target, we will consider the following cases: (a) ∆ + = 0 -Drell-Yan frame (DYF); Note that in both cases the normalization factor appearing in Eq. (100) reduces to N ≈ P + .

C. Drell-Yan frame
The Drell-Yan frame defined by ∆ + = 0 can be seen as the LF version of the elastic frame defined byP ·∆ = 0. Distributions in the DYF can be obtained by integrating the static LF EMT over the longitudinal LF coordinate x − T µν a (b ⊥ ;P ) = dx − T µν a (x;P ) = where b ⊥ = x ⊥ is the same impact parameter as in instant form. The unpolarized off-forward LF amplitude (97) in the DYF takes the form If we further integrate the static LF EMT over the impactparameter b ⊥ , which amounts to setting ∆ ⊥ = 0 ⊥ in Eq. (102), we recover the FL expression (22) with E P replaced by P + .
However for ∆ ⊥ = 0 ⊥ , unlike P 0 in instant form (50), P + in front form is an independent variable which does not depend on ∆ ⊥ . One can therefore write a relatively simple expression for the static LF EMT in terms of Fourier transforms of GFFs 7 . The (1+2)-dimensional Galilean "mass" current take the simple form Its LF time component J + a = T ++ a can be interpreted as the density of longitudinal LF momentum carried by constituent type a and has been discussed in [49,61,62,[116][117][118][119][120][121][122]. As expected, the distribution of longitudinal LF momentum coincides with the properly rescaled instant-form energy density in IMF µ a (b) = ρ IMF (b)P + /M . For the (1+2)-dimensional Galilean EMT, we find Note that the Galilean EMT T αβ a has a global boost factor M/P + . We can get rid of this factor by integrating T µν a (x;P ) over the coordinate r L = r − P + /M invariant under longitudinal LF boosts instead of r − . The structure of the Galilean EMT then looks like where n α = η α− ,nβ = ηβ + , and This tensor can alternatively be written as with l αβ = n αnβ − η αβ . We can therefore define P + -independent two-dimensional Galilean energy density, radial pressure, tangential pressure, isotropic pressure, and pressure anisotropy associated with constituent type a as The two-dimensional pressures are the same as the ones obtained in the BF (57)- (60). For this reason, we used the same notation as in instant form. It is also not surprising that the energy densities ρ a and µ a defined, respectively, in instant form and front form differ because they are simply related to different components of the EMT. The latter is illustrated in Fig. 18 using the multipole model (11) with parameters given in Table I. The difference in the behavior around b = 0 between the quark and gluon contributions comes from the fact D g (t) = 0 whereas D q (t) = 0. One might also be at first sight puzzled by the fact that total LF energy is given by a d 2 b ⊥ µ a (b) = M/2 instead of M . This simply comes from our definition of the LF components which implies that d 3x T +− (x;P ) = M 2 /2P + .

D. Infinite-momentum frame
In Ref. [129,130] three-dimensional LF distributions for finite P + have been defined by means of a Fourier transform with respect to the longitudinal LF boost-invariant skewness variable ξ = −∆ + /2P + instead of ∆ + . The problem with these distributions is that the (quasi-)probabilistic interpretation is lost owing to the non-vanishing LF energy transfer ∆ − = 0. Moreover, the center of the target with respect to which the transverse coordinates are defined differs between the initial and final states when ξ = 0 [131].
The above issues disappear when ξ = 0. In the literature, one usually considers finite P + and hence ∆ + = 0, reducing the LF distributions to two spatial dimensions. This option was discussed in the previous section. The other possibility is to consider P + ∆ + , √ P 2 , i.e. the IMF within the LF formalism. In this case, we can formally define another type of three-dimensional LF distributions free of the aforementioned problems. To the best of our knowledge, this is the first time that such an option is explored. The reason can likely be attributed to the fact that one usually has in mind reaching the IMF through an infinite longitudinal boost. In that case, both P + and ∆ + will get large with their ratio ξ fixed. Actually, one should just consider P + → ∞ with ∆ + fixed. The information about the longitudinal spatial structure is then encoded around ξ ≈ 0.
Expanding the unpolarized off-forward LF amplitude (100) in powers of 1/P + , we find Its Fourier transform T µν a can be expressed in terms of 3-dimensional Fourier transforms of GFFs. The expression we obtain is however so complicated that we were not able to recognize the EMT structure of any known continuous medium discussed in the literature. Note also that since ∆ − ∝ 1/(P + ) 2 , we can write t ≈ −∆ 2 ⊥ . The GFFs therefore do not contribute to the ∆ + -dependence in the LF IMF and the longitudinal structure is essentially determined by Lorentz symmetry. As a final remark, we naturally recover the DYF results in the limit ∆ + → 0. This means that the projection of the distributions defined in the LF IMF onto the transverse plane coincides with the two-dimensional distributions defined in the DYF.

V. DISCUSSION
Having defined the notions of energy density and pressure inside the nucleon, we can go on and discuss the questions of hydrostatic equilibrium and stability constraints. While the former are automatically satisfied once GFFs are determined, the latter provide new constraints particularly useful for the phenomenology of high-energy scatterings involving nucleons.
A. Hydrostatic equilibrium

Three-dimensional case
Conservation of the total EMT ∂ µ T µν (x) = 0 implies that the static total EMT T µν = a T µν a satisfies in the BF or equivalently dp r (r) dr = − 2s(r) r , (114) which is the equation of hydrostatic equilibrium in the presence of pressure anisotropy 8 . The RHS of Eq. (114) represents the effects of a force arising from the anisotropic nature of the medium. When p r (r) < p t (r), we get s(r) < 0 and the force is repulsive 9 . When p r (r) > p t (r), we get s(r) > 0 and the force is attractive. This force is the mechanical origin of the surface tension between a liquid and its vapor [133]. In the bulk of ordinary fluids, the density is essentially constant and the pressure is isotropic. The density changes however drastically across the interface, and generates an asymmetry in the stress tensor which can be understood as arising from a difference of ranges between attractive and repulsive interactions. The interface being usually extremely thin for ordinary fluids, the asymmetry can usually be modelled by a simple surface tension. This picture is supported by a molecular dynamics simulation of molecules interacting via a Lennard-Jones potential [133,134], which shows that the transition in the time-averaged density profile from the high-density liquid to the low-density gas takes place in a very narrow region that is a few molecules wide. At the same time, the profile of the stress anisotropy indicates that there is a force localized in the same narrow region acting in the direction parallel to the interface. If r is the coordinate normal to the interface, the surface tension γ is obtained from the Bakker equation [135] (also known as the Kirkwood and Buff method [136]) where D is the domain where s(r) is significant. When the domain is very narrow, we can make the approximation s(r) ≈ γ δ(r − R). Integrating then the equation of hydrostatic equilibrium (114) over the radial coordinate yields the well-known Young-Laplace relation p r (0) = 2γ/R for a spherical drop with radius R [137]. For compact stars and hadrons, the anisotropic stress spreads over a significant fraction of the volume of the system and cannot be realistically approximated by a delta function.
Many relations can be derived from Eq. (114), as discussed in [11,26]. Let f (r) be some radial function, one finds using integration by parts For f (r) = 1, one obtains the generalized Young-Laplace relation If p r (r) is finite at r = 0 and decays faster than 1/r N with N > 0 for r → ∞, the choice f (r) = r N leads to The case N = 3 is known as the von Laue condition [138] ∞ 0 dr r 2 p(r) = 0 , 8 For spherically symmetric compact stars, the line element in Schwarzschild coordinates is given by ds 2 = e ν(r,t) dt 2 − e −λ(r,t) dr 2 − r 2 dθ 2 + sin 2 θ dϕ 2 , and the equation of hydrostatic equilibrium, which directly derives from the vanishing of the covariant divergence of the EMT in the static limit, reads dpr(r) dr = − 1 2 dν(r) dr ε(r) + pr(r) − 2s(r) r (115) in the presence of anisotropic matter [19]. If we switch off gravitational effects by sending Newton's constant to zero, both functions ν(r, t) and λ(r, t) vanish and we recover Eq. (114). 9 This seems to be the favored case for neutrons stars since it allows the construction of more compact and more stable objects than with ordinary isotropic matter [23,132], relaxing therefore the tension between observations and theoretical bounds. and indicates that the isotropic pressure has to change sign. We illustrate this in Fig. 19 using the multipole model (11) with parameters given in Table I indicates that the tangential pressure also changes sign. The net tangential force in the inner region is repulsive whereas it is negative in the outer region, leading once more to the conclusion that C(0) < 0 based on Eq. (48). The other remarkable values are N = 1, 6 5 , 12 All the above relations are automatically satisfied by our expressions (35) and (38) when summed over the constituents. The reason for this is that the parametrization (8) together with the sum rules (9) already encode the conservation of the total EMT, implying that the equation of hydrostatic equilibrium (114) is identically satisfied. These relations can however be particularly useful to test model predictions where full Lorentz covariance are often absent.

Two-dimensional case
The discussion in the EF proceeds analogously to the BF. The main difference is that the spatial distributions are now two-dimensional instead of three-dimensional. Note also that since the two-dimensional pressure distributions are the same in both EF (with P z = 0 or P z → ∞) and DYF, the following results apply to both instant and front forms.
Using the conservation of the static total EMT in two dimensions we find that the equation of hydrodynamic equilibrium in the EF takes the form where we omitted the P z -dependence for convenience. If f (b) is some radial function, we obtain using integration by parts For f (b) = 1, we get When the pressure anisotropy is concentrated within a thin region, it can be approximated by Π(b) ≈ τ δ(b − R). Equation (130) reduces then to the two-dimensional version of the Young-Laplace relation σ r (0) = τ /R, where τ can be thought of as some sort of effective string tension. More generally, the effective string tension can be defined from the two-dimensional version of the Bakker equation If σ r (b) is finite at b = 0 and decays faster than The case N = 2 corresponds to the two-dimensional version of the von Laue condition and indicates that the isotropic pressure has to change sign. We illustrate this in Fig. 20 using the multipole model (11) with parameters given in Table I indicates that the tangential pressure also changes sign. Once again, the net tangential force in the inner region is repulsive whereas it is negative in the outer region, in agreement with C(0) < 0 according to Eq. (72). The other remarkable values are For the same reason as in the three-dimensional case, all the above relations are automatically satisfied by our expressions (59) and (60) when summed over the constituents.

B. Stability
We have seen earlier that the study of the nucleon EMT might provide some clues about the EoS for the matter lying in the heart of compact stars. The stability of compact stars made of anisotropic matter has been extensively discussed in [20,139,140]. We suggest that applying these results to the case of the nucleon can in turn provide new constraints on the nucleon EMT, and hence on the GPDs. For a stable system, it is expected that 10 (i) ε(0) < ∞, p(0) < ∞ and s(0) = 0 ; (ii) ε(r) > 0 and p r (r) > 0 ; (iii) dε(r) dr < 0 and dpr(r) dr < 0 .
All these constraints are satisfied by our simple multipole model (11) with parameters given in Table I. We observe in particular that the constraint (i) rules out the dipole Ansatz for the GFFs B a (t) and C a (t) sometimes used in the literature, because it generates a 1/r pole in ε(r) and p(r), and leads to s(0) = 0. We therefore used in our simple multipole model the tripole Ansatz which does not have the same problem and which agrees with the asymptotic behaviour expected from the quark counting rules [141][142][143][144].
Violations of these additional constraints are observed over some range in r within our simple multipole model. 10 Note that for compact stars it is also expected that pt(r) > 0, and hence p(r) > 0. This does not contradict Eqs. (120) and (121) since the gravitational force is attractive and long range, leading to a significant stress anisotropy. In some sense, we could interpret the gravitational contribution in Eq. (115)  Finally, there exist also energy conditions which reflect the principles of relativity and play an important role in General Relativity. They constitute an essential ingredient for establishing general results like the no-hair theorem, the laws of black hole thermodynamics or the singularity theorems of Penrose and Hawking [126,[145][146][147]. The most popular ones are the Null Energy Condition (NEC), Weak Energy Condition (WEC), Strong Energy Condition (SEC), and Dominant Energy Condition (DEC) NEC ε(r) + p i (r) ≥ 0 , (139) WEC ε(r) + p i (r) ≥ 0 and ε(r) ≥ 0 , SEC ε(r) + p i (r) ≥ 0 and ε(r) + 3 p(r) ≥ 0 , where i = r, t. Some of these energy conditions have a simple physical interpretation. Namely, the weak energy condition arises from the requirement that the energy density is non-negative for any observer, and the dominant energy condition ensures that the energy flow cannot exceed the speed of light for any observer. All known forms of matter so far satisfy these energy conditions. Our simple multipole model (11) with parameters given in Table I also satisfy these energy conditions except the SEC ε(r) + 3 p(r) ≥ 0 for some range in r.
All the above conditions on the distributions of energy density and pressure are extremely interesting since, when transposed to the nucleon case using our expressions (34)- (38), they can provide new phenomenological constraints on the GFFs and hence on the GPDs of the nucleon 11 . Recently, criterion (ii) has been considered in [84] and led to the conclusion that C(0) should be negative owing to Eq. (48). We note that this agrees with criterion (iii) which implies that s(r) > 0 (as assumed in [26]) using the equation of hydrostatic equilibrium (114), and in turn C(0) < 0 according to Eq. (47). Note also that the inequalities can in principle be easily transposed to the two-dimensional case, which in conjunction with our expressions (56)-(60), may lead to yet further new constraints on the GPDs.
Satisfying all the constraints at the same time is not at all a trivial task. Some may even perhaps be unapplicable to the nucleon case. For these reasons, we refrain from developing here a more realistic model since its main purpose in the present study was only to illustrate the various distributions. A detailed analysis of the stability constraints in the hadronic context goes beyond the scope of the present and is left for future investigations.

VI. CONCLUSIONS
We revisited the interpretation of spin-1/2 gravitational form factors in terms of the mechanical properties of hadrons. We rederived and significantly extended the existing literature, which so far has been exploiting only the Breit frame and the instant form of dynamics. In particular we discussed here both the instant and front forms of dynamics, and the separate quark and gluon contributions to the energy and pressure distributions in the Breit, elastic, infinite-momentum and Drell-Yan frames. Our key results are contained in Eqs. (34)- (38) for the Breit frame, Eqs. (56)-(60) for the elastic frame, and Eqs. (107)-(111) for the Drell-Yan frame. We illustrated our argument with a simple phenomenological model, and highlighted the constraints coming from mechanical properties that should generically be satisfied in model-building.
We put a special emphasis on the pressure anisotropy, which offers an innovative perspective on the nucleon structure from our rapidly growing knowledge on compact stellar objects. The direct observation of neutron star mergers in terrestrial gravitational wave observatories indeed already brought constraints on the equation of state of nuclear matter at high density and low temperature. Orders of magnitude and model studies suggest that the nucleon itself may be described with the same concepts and pictures. In particular, the stability and hydrostatic equilibrium conditions are presumably the best theoretical ingredients to elaborate on this picture since most of the gravitational form factors are within contemporary experimental reach through hard exclusive experiments and generalized parton distributions. However, determining how much we can learn about the physics of compact stars from the nucleon structure (or conversely) is still an exciting open question.