Nucleon D-term in holographic QCD

The D-term is one of the conserved charges of hadrons defined as the forward limit of the gravitational form factor $D(t)$. We calculate the nucleon's D-term in a holographic QCD model in which the nucleon is described as a soliton in five dimensions. We show that the form factor $D(t)$ is saturated by the exchanges of infinitely many $0^{++}$ and $2^{++}$ glueballs dual to transverse-traceless metric fluctuations on the Wick rotated AdS$_7$ black hole geometry. We refer to this phenomenon as `glueball dominance', in perfect analogy to the vector meson dominance of the electromagnetic form factors. However, the value at vanishing momentum transfer $D(t=0)$ can be interpreted as due to the exchange of pairs of pions and infinitely many vector and axial-vector mesons without any reference to glueballs. We find that the D-term is slightly negative as a result of a cancellation between the isovector and isoscalar meson contributions.


Introduction
The nucleons (protons and neutrons) have a number of conserved charges. They have mass of about 1 GeV associated with translational invariance. They have spin 1/2 associated with rotational invariance. They have electric charges and magnetic moments from electromagnetic gauge invariance. They have a baryon number +1 associated with the global U(1) symmetry. There are also approximately conserved charges related to isospin (flavor) symmetry. All these charges are well known and very accurately measured, and their values are well documented in textbooks and literature [1]. There is, however, one exactly conserved charge, the so-called D-term, whose value is presently unknown. The D-term is the forward limit of one of the gravitational form factors D(t = −k 2 ) defined through the off-forward hadronic matrix element of the QCD energy momentum tensor p |T µν |p (k = p − p), see [2] for a recent review. It was first discovered in the 60's [3,4] but has kept a remarkably low-profile status known only to a subset of theorists in the nuclear physics community. The situation has changed drastically in the past several years, mainly fueled by the anticipation of the future Electron-Ion Collider (EIC) experiments [5][6][7] dedicated to the study of the nucleon structure. In particular, the first attempt to extract the D-term from experimental data has been made [8].
The value of the nucleon D-term has long remained unknown because there is no known way to directly measure it. This would require a controlled experimental setup to scatter a nucleon from gravitons, which is possible only in thought experiments. Yet, there are indirect ways to measure the D-term in lepton-nucleon scattering. The price to pay, however, is that one can only access the quark and gluon components separately from different experiments where D u is from the up-quark part of the energy momentum tensor, D g is from the gluon part, etc. Each component depends on the renormalization scale, and only the sum is renormalizationgroup invariant. The light-quark contribution D u,d can be in principle accessed in Deeply Virtual Compton Scattering (DVCS) by separately measuring the real and imaginary parts of the Compton form factors [9,10]. The gluon contribution D g , on the other hand, can be accessed in nearthreshold photo-and lepto-production of heavy quarkonia such as J/ψ and Υ [11][12][13][14][15][16][17][18] (see however, [19][20][21][22][23]). Similarly, the strangeness contribution D s can be probed in near-threshold φ-meson leptoproduction [14]. Alternatively, these components can also be calculated in lattice QCD simulations, see for example [24] and references therein.
Although the magnitude and even the sign of the D-term are unknown, a general expectation is that it is negative. This is based on an analogy to the mechanical stability analysis of classical systems. Being the spatial component of the energy momentum tensor, after Fourier transforming to the coordinate space the D-term form factor may be associated with the 'shear' and 'pressure' of a classical spherical system [2,25]. A positive D-term would imply overall positive outward force, making the system unstable. We however note that there is no field-theoretical proof of the connection between the negativity of the D-term and the stability of hadronic bound states. Besides, the term 'pressure' should not be taken literally in its usual sense in thermodynamics.
In this paper, we calculate the (total) nucleon D-term in the chiral limit using gauge/string duality based on a holographic QCD model proposed in [26,27]. The model is 'top-down', meaning that it has been directly derived from string theory, and is realized in a D4/D8 brane configuration in type IIA superstring theory. As such, it does not rely on ad hoc assumptions and allows for a systematic truncation and/or inclusion of various higher order corrections (stringy effects, 1/N c corrections, etc.). In this model, the baryons are realized as solitons in a five-dimensional gauge theory [26,[28][29][30][31], and various properties including electromagnetic form factors [32][33][34] have been investigated using this description. As a matter of fact, holographic approaches are particularly suited for the study of the gravitational form factors because one can literally exchange gravitons, albeit in extra dimensions. Indeed, there have been previous attempts to compute the D-term in 'bottom-up' holographic models [19,35,36]. However, the outcomes of these studies are rather mixed: Ref. [35] found that the D-term was simply zero. Ref. [19] argued that the Dform factor was proportional to the A-form factor (defined in (2.2) below), but the proportionality constant remained undetermined. The latter work has been recently revisited in [18,36] where the proportionality constant was found to be subleading in the 1/N c expansion. There is also a work based on an AdS/QCD-inspired quark-diquark model [37], but holography was not used in actual calculations. In view of this, it is worthwhile to see what top-down holographic models have to say about the D-term. Our model is sophisticated enough to accommodate infinite towers of meson, baryon and glueball resonances. We shall be particularly interested in how these degrees of freedom contribute to the D-term.
Our calculation bears some resemblance to those in the chiral soliton model and the Skyrme model [38][39][40][41]. This is so because a baryon in our model is described as a soliton (an instanton) in five-dimensions, similarly to the Skyrmion in four dimensions. In fact, it is known that one can derive the Skyrme model from our model [26][27][28].
In the next section, we give a brief review of the gravitational form factors in QCD. In Section 3, we introduce our model and take a first look at the energy momentum tensor in this model. In Section 4, we calculate the D-term in the 'classical' approximation by Fourier transforming the soliton energy momentum tensor. Then in Section 5, we discuss the form factor D(t) using holographic renormalization and establish its connection to scalar and tensor glueballs in QCD. Finally, in Section 6 we conclude with physics interpretations and future perspectives.

Gravitational form factors
In this section, we quickly introduce the nucleon gravitational form factors and set up our notations. More details can be found in a recent review [2]. In accordance with the string theory literature, we use the 'mostly plus' metric η µν = (−1, 1, 1, 1). The QCD energy momentum tensor then takes the from where a = 1, 2, .., N 2 denotes symmetrization.ψ = ψ † β = iψ † γ 0 and the gamma matrices satisfy the Dirac algebra γ µ γ ν + γ ν γ µ = 2η µν . The nucleon gravitational form factors are defined by the off-forward matrix element of the energy momentum tensor [3,4] . M is the nucleon mass. The nucleon spinors are normalized asū(p)u(p) = 2M . (2.2) is the most general parameterization given that T µν is symmetric and conserved ∂ µ T µν ∼ k µ T µν = 0. Energy conservation also implies that the three form factors A, B, D are renormalization group invariant. Their values at t = 0 are of particular interest. It is known that A(t = 0) = 1 from momentum conservation and B(t = 0) = 0 from angular momentum conservation. However, the value D(t = 0) is not constrained by any symmetry, and is currently unknown. Our main goal of this paper is to study D(t) in the holographic QCD model proposed in [26,27].
We shall be working in the Breit frame where p = k/2 = − p, so that P = 0 and k 0 = p 0 −p 0 = 0. In this frame,ū(p )u(p) = 2p 0 δ s s = 2 M 2 + k 2 /4 δ s s , where s , s = ± 1 2 denote the spin states of the nucleon, and Therefore, the D-term can be obtained by reading off the coefficient of For a later discussion, following [11], let us also introduce the 'transverse-traceless' (TT) part T µν TT of the energy momentum tensor. It is defined as a part of T µν that satisfies the conditions ∂ µ T µν TT = (T TT ) µ µ = 0. Explicitly, Eq. (2.2) can then be rewritten as The first line on the right hand side is the matrix element of T µν TT and the second line is from the trace part. Note that the structure k µ k ν − k 2 η µν characteristic of the D-term is now present in both parts.
Before leaving this section, for the convenience of the reader, we note the definition of the gravitational form factors in the 'mostly minus' metric η µν = (1, −1, −1, −1) which is used in most QCD literature. The gamma matrices in the two conventions are related as γ µ mostly minus = iγ µ mostly plus . (2.6) With the QCD energy momentum tensor in this metric

Nucleon in holographic QCD
The two-flavor (N f = 2) meson-baryon sector of our model [26,27] is defined by a U(2) gauge theory in a curved five-dimensional spacetime (x µ=0,1,2,3 , z) supplemented with the Chern-Simons (CS) term are the warp factors along the fifth dimension. F µν = F a µν τ a 2 , F µz = F a µz τ a 2 are the SU(2) field strength tensors (τ a=1,2,3 being the Pauli matrices normalized as tr[τ a τ b ] = 2δ ab ), and F µν is the U(1) field strength tensor. The Lorentz indices of these gauge fields are raised and lowered by the flat five-dimensional metric (−1, 1, 1, 1, 1). The Chern-Simons term S CS prescribes the interaction between the SU(2) and U(1) fields, but its explicit form is not needed for the present discussion. The parameter is proportional to the number of colors N c and the 't Hooft coupling λ = g 2 N c . All dimensionful scales have been made dimensionless by appropriately rescaling by the model's only mass parameter M KK (e.g., x µ → x µ = x µ M KK ). These parameters was determined in [26,27]  Mesons with isospin quantum numbers π, ρ, a 1 , · · · are described by the fluctuations of the SU(2) gauge field F µν . We consider the chiral limit, so the pions are massless. Iso-singlet mesons ω, η , · · · are described by the U(1) field. A baryon is realized by a static (independent of t = x 0 ), soliton-like configuration of F µν which satisfies the equation of motion of the five-dimensional gauge theory (3.1) [26,30]. This is charged under the U(1) gauge field A µ through the CS term, and the charge is identified with the baryon number. At strong coupling λ 1 and in the small-z region where the metric is approximately flat h(z) ≈ k(z) ≈ 1, the solution is simply given by the BPST instanton [42] in four-dimensional Euclidean space (x i=1,2,3 , z) [30]

5)
ρ is the instanton 'size' and ( X, Z) denotes the 'center' of the instanton. When the soliton is quantized, these parameters, together with the 'orientation' of the instanton in the flavor SU(2) space F µν → V F µν V −1 , are promoted to time-dependent operators (for example, ρ →ρ(t)), a procedure known as the collective coordinate quantization. While this has been done in previous applications of the model [30,32], in this work we shall treat them as c-numbers, leaving their quantum treatment for future work. This means that we eventually set X = 0 (without loss of generality) and employ the value Z = 0 which minimizes the soliton potential in the Z-direction [30]. It should be kept in mind that, by neglecting quantization, we are effectively treating the nucleon as a scalar particle because the quantization of the SU(2) orientation is what makes the soliton a spin-1/2 particle. The D-form factor exists also for scalar hadrons, but the B-form factor (see (2.2)) does not. In order to compute the latter, soliton quantization is crucial.
When |z| 1, the flat space approximation h ≈ k ≈ 1 breaks down. Exact analytical solutions in this region are no longer available. However, when |z| 1, the following approximate solution has been constructed [32] in the so-called singular gauge In [32], the non-Abelian commutator terms in F g iz and F g ij have been dropped since they are negligible when |z| 1. Here we have restored them for a later purpose. The Green functions In order to compute the D-term, or more generally the gravitational form factors, it is desirable to have an approximate solution which smoothly interpolates the above solutions in the two limits |z| 1 and |z| 1. Such a solution can be readily found for the U(1) part. We start with the small-z region and consider the following equation of motion for A 0 in the SU(2) instanton background 9) or in the Fourier space, In the flat space approximation h(z) ≈ k(z) ≈ 1, the solution regular at ξ = 0 is given by (3.5).
The boundary conditions are such that A 0 (z → ∞) = 0 and the solution is smooth at z = 0, namely, ∂ z A 0 (z)| z=0 = 0. When |z| 1, the right hand side of (3.9) is negligible, and the equation becomes identical to that for G in (3.8). Therefore, the solution of (3.9) with the said boundary conditions smoothly interpolates the solutions at |z| 1 and |z| 1. We will use it as an approximate solution in the whole region 0 < |z| < ∞.
The situation is more complicated for the SU(2) part. The asymptotic solution (3.6) with (3.8) has been obtained by neglecting the nonlinear terms in the Yang-Mills equation. In the smallinstanton regime ρ 1, or equivalently, the strong coupling regime λ 1 (note the correspondence ρ ∼ 1/ √ λ [30]), the large-z and small-z solutions have an overlapping region of validity ρ z 1 where they can be smoothly matched [32]. We extrapolate G and H to small-z by solving (3.6) with the instanton configuration (3.5) rotated to the singular gauge A µ → A g µ and substituted into the left hand side. This gives (3.12) In the ρ → 0 limit, G, H reduce to the flat space Green function −1 4π 2 ξ 2 in four dimensions. Regarding the right hand side of (3.12) as a regularized form of δ( x − X)δ(z − Z), we are led to consider the following equations instead of (3.8). 5 Eq. (3.13) is to be solved with boundary conditions G, H → 0 as |z| → ∞ and ∂ z G| z=0 = ∂ z H| z=0 = 0. At large-z, the right hand side is negligible, and the equation reduces to (3.8). At small-z, G and H smoothly connect to the instanton solution by construction. Therefore, we can use (3.6) with (3.13) as an approximate solution in the whole range 0 < |z| < ∞ even when ρ is order unity, provided the non-Abelian commutator terms in F ij and F iz in (3.6) are kept. To go beyond this approximation, one has to numerically solve the Yang-Mills equation in the curved background as was done in [43][44][45].
Let us now take a first look at the energy momentum tensor in this model. We can write down the following 'classical' energy momentum tensor in four-dimensions by varying the action (3.1) with respect to the flat metric η µν [46]. The Chern-Simons term does not contribute since it is independent of the metric. When evaluated on-shell, (3.14) is conserved We have explicitly verified (3.15) using the equation of motion derived in [32].
We expect that (3.14) is a reasonable approximation to the full result, when the baryon is treated as heavy (at least parametrically) as in large-N c QCD. (In Section 5, we shall discuss the corrections to this formula due to glueballs.) In particular, the nucleon mass can be calculated as where we have already substituted the classical configuration A 0 = A i = A z = 0. Eq. (3.16) is consistent with Eq. (3.18) of [30] after using the equation of motion. The latter expression is simply the minus of the on-shell action dtM = −S (this time including the Chern-Simons term), which is appropriate for a classical, heavy particle at rest. For the instanton solution, the integrals in (3.16) can be explicitly evaluated and lead to the structure [30] (Remember the mass is measured in units of M KK .) The leading term 8π 2 κ ∼ O(λN c ) comes from the SU(2) part. The subleading terms ρ 2 and 1/ρ 2 come from the SU(2) and U(1) fields, respectively. The SU(2) fields tend to shrink the instanton size ρ → 0, while the U(1) field tends to expand the size ρ → ∞, making the system unstable. As a compromise, the minimum energy is achieved when The situation is entirely analogous to the Skyrme model where the baryon (realized as a solitonic configuration of pions) is stabilized by introducing the ω-meson [47]. The numerical value of ρ has been fixed in this way in the instanton approximation [30] ρ = 27π λ Together with (3.4), this gives M ∼ 1.5 GeV which is larger than the observed value. However, the value of M is subject to changes after the collective coordinate quantization [30].
Returning to (3.14), our main interest in this paper is the spatial i, j = 1, 2, 3 components Classically, one might expect that the form factor D(| k|) could be calculated by simply inserting the above solutions into (3.19) and Fourier transforming to momentum space x → k. (Below we shall often use the notation D(| k|) instead of D(t).) Since T cl ij is conserved, it must have the structure However, it turns out that this naive approach is valid only at k = 0. In the next section, we shall numerically evaluate the D-term D(0) in this way. A more general analysis valid for arbitrary values of | k| will be presented in Section 5.

'Classical' calculation of the D-term
In this section, we calculate the D-term 'classically', by Fourier transforming the naive energy momentum tensor (3.19) with the approximate solutions A 0 , F ij and F iz constructed in the previous section. This is analogous to what has been done in the Skyrme model or chiral soliton models. For a reason to be clarified in the next section, the calculation in the present section is valid only at vanishing momentum transfer | k| = 0. When | k| = 0, one must carry out a fully holographic calculation which will be discussed in the next section. We shall calculate the contributions from the U(1) and SU(2) fields separately In the remainder of this section, only three-momenta k, q appear. We thus write | k| ≡ k, | q| ≡ q below to simplify the notation.

U(1) part
In Section 3, we have constructed an approximate solution for the U(1) gauge potential A 0 (x, z) which can be used in the entire range −∞ < z < ∞. This is obtained by numerically solving (3.10) with the boundary conditions A 0 (k, z = ∞) = 0 and ∂ z A(k, z)| z=0 = 0. The next step is to substitute this solution into (3.19).
In the second equality we integrated by parts in z. The coefficients P 1 and Q 1 can be calculated as follows Despite the singular prefactor 1/k 2 , P 1 (k = 0, z) is finite as one can see by expanding the integrand in k and performing the angular integral. On the other hand, Q 1 (k = 0, z) is actually divergent. If the bulk equation of motion is solved exactly, this divergence should be canceled by the contribution from the SU(2) field. Moreover, after the cancellation the coefficients of k i k j and −δ ij k 2 must agree exactly due to the conservation law (cf. (3.20)). Indeed, after some manipulations (including the addition of total derivative terms), one can show that The expression inside the brackets is connected to the SU(2) fields via the equation of motion (3.10). In practice, since our solution is approximate and obtained only numerically, it may be difficult to achieve a precise cancellation. We thus focus on the coefficient of k i k j which is safely calculable in the present approach, and leave a more complete analysis for future work.

SU(2) part
We now turn to the SU(2) fields. Let us first point out that the instanton solution, valid in the small-z region, does not give rise to the structure where the z-integral should be cut off around z O(1). Note that the coefficient of δ ij is finite at k = 0, meaning that (4.6) gives a divergent contribution to D(k = 0) since T ij ∼ D(k)δ ij k 2 . This is expected in view of our discussion in the previous subsection. If the equation of motion is solved exactly, the divergence must be canceled by the contributions from the U(1) field as well as that from the SU(2) field in the large-z region (where the instanton approximation breaks down). However, since we decided to focus on the coefficient of k i k j , we do not dwell on (4.6) further.
To calculate the D-term, we need to extrapolate the above solution to the small-z region. The main complication is the nonlinearity of the SU(2) Yang-Mills equation. We have argued in Section 3 that, at least in the strong coupling limit where 1 ρ, the linear approximation is valid up to z ρ, and an approximate solution in this regime can be obtained by replacing (3.8) with (3.13).
The remaining parametrically small region ρ > z > 0 does not contribute to the coefficient of k i k j as we have seen in (4.6), so we may further extrapolate this solution down to z = 0. However, as the coupling is lowered and ρ exceeds unity, the non-Abelian commutator terms in F ij and F iz become important. This is why we have kept them in (3.6).
Since the full SU(2) energy momentum tensor is local and involves up to four powers of the gauge fields, it is much more advantageous to work in the original coordinate space. In (3.19), there are two terms that can give rise to the structure k i k j . Noting that G and H depend only on the magnitude x = | x| (and z), we write After a rather tedious calculation we find ∂ Z G can be eliminated by using the formula (see Eq. (2.79) of [32] and (3.13)) after which we may set Z = 0. The first relation was originally derived in the large-z, linear regime, but it is also valid in the small-z regime where Given X, we can evaluate the SU(2) contribution to the D-term as where j 2 is the spherical Bessel function.

Numerical result
We have solved (3.10) and (3.13) numerically with the Neumann boundary conditions at x = 0 and z = 0 and the Dirichlet boundary conditions A 0 = G = H = 0 at infinity. Great care is needed in order to obtain a stable solution for H (and especially its x, z-derivatives) because of the massless pion pole, the n = 0 term in (4.10), as well as the pole 1/ξ 2 on the right hand side of (3.13). To cope with these, we first make the following shift (cf. (3.11)) to avoid the singular behavior near the origin. 6 We then solve the resultant differential equations forG(x, z) andH(x, z) by employing the pseudospectral method [48]. In this method, the functions G(x, z),H(x, z) are expanded by tensor products u m (x; L x )u n (z; L z ) (1 ≤ m, n ≤ N ) where the basis functions u n (y; L) = y 2 /(y + L) 2 TL n (y; L) − 1, which individually satisfy the boundary conditions, are related to the rational Chebyshev functions TL n (y; L) on the semi-infinite interval 0 < y < ∞ via the "basis recombination" [48]. The adjustable parameters L x and L z are set to 1 for solvingG(x, z) while we use L x = 16 and L z = 1 forH(x, z) to better stabilize the x-dependence of H. The number of basis functions N is varied in the range 60 ≤ N ≤ 120, and the variation in the results is used as an estimate of systematic errors for each value of k. Libraries SciPy [49], Eigen [50] and Boost.Multiprecision [51] have been helpful for these numerical analyses.
The solutions just described depend on the input value of ρ. For the instanton solution, ρ is given by (3.18). Since we go beyond the instanton approximation, ρ needs to be recalculated accordingly. For our new solution, the nucleon mass (3.16) takes the form We find M (ρ) has a minimum at (compare to (3.18)) We have thus evaluated the integrals (4.3) and (4.16) using the solutions A 0 , G and H with ρ = ρ * , and extrapolated them to k → 0 to obtain where the errors for the U(1) part are negligibly small. After a cancellation between the positive U(1) and negative SU (2) contributions, the total D-term turns out to be slightly negative. That the U(1) contribution is positive is intuitively easy to understand. The U(1) field is analogous to the static electric field E ∼ r/r 3 of a point charge. The energy momentum tensor (Maxwell's stress tensor) of a point charge takes the form The 'pressure' p(r) ∼ 1/(6r 4 ) (see (1.2)) is everywhere positive and hence the D-term D ∼ d 3 rr 2 p(r) is also positive (actually divergent). On the other hand, the SU(2) fields may be thought of as the 'pion cloud' in traditional hadron physics. In the chiral quark soliton model, it has been argued that the pion cloud is responsible for making the D-term negative [38]. Our result is consistent with this argument, although in the present approach the 'cloud' is made up of not only pions but also infinitely many vector and axial-vector mesons. Incidentally, if we neglect the non-Abelian commutator terms (terms proportional to ρ 2 in (4.14)), the SU(2) contribution also becomes positive.
Our result provides a new perspective on the stability of nucleons in holographic QCD. If we interpret the D-term as a measure of outward radial force, the U(1) and SU(2) fields generate positive and negative forces, respectively, and they tend to expand and shrink the system. This is entirely analogous to the calculation of the nucleons mass as briefly mentioned below (3.17) and discussed in detail in [30]. Namely, the U(1) field (iso-singlet mesons, in particular the ω meson) tends to expand the nucleon by preferring large instanton sizes ρ → ∞, and this is counterbalanced by the SU(2) fields (iso-vector mesons π, ρ, a 1 , · · · ) which prefer small sizes ρ → 0. Neither one of them alone can stabilize the nucleon. Thus, there seems to be a direct link between the stability arguments in terms of the nucleon's mass and D-term when they are decomposed into contributions from different subsystems. On the other hand, the present discussion does not indicate whether the sign of the total D = D U(1) + D SU (2) , which happens to be slightly negative, is of particular significance regarding stability (cf., [52,53]).

Coupling to gravity
In gauge/string duality, the proper method to calculate the field theory expectation value of the energy momentum tensor has been established [54,55]. In this section, we apply the framework of holographic renormalization developed in [55,56] to the calculation of the form factor D(| k|) and elucidate its connection to the glueball spectrum. We shall also explain why the classical calculation in the previous section is valid only at | k| = 0. For previous attempts in bottom-up holographic models, see [19,36].

Setup
The basic idea of our holographic calculation is that matter fields in the 'bulk' perturb the metric, and this 'wake' is propagated to the boundary and recorded as the field theory expectation value T µν . In bottom-up holographic QCD models, gravity is confined in the same five-dimensional (deformed) anti-de Sitter spaces where the matter fields live. However, the situation is different in our top-down model. The action (3.1) is a low-energy effective theory on the 'flavor' D8 branes embedded in a ten-dimensional curved space-time in type IIA supergravity. The latter is further derived from eleven-dimensional supergravity (M-theory) on doubly Wick rotated AdS 7 black hole ×S 4 [57] after compactifying the eleventh dimension x 11 on a circle of radius where g s is the string coupling and l s = √ α is the string length. The following relation is useful The radial coordinate r (∞ > r > R) is related to z in (3.1) as Below we shall use r and z interchangeably keeping in mind that they are related as (5.5).
The AdS 7 black hole geometry is sourced by the N c D4 branes spanning the coordinates (x µ=0,1,2,3 , τ ). In this background, N f D8 branes are placed at τ = 0 and τ = π/M KK corresponding to the z > 0 and z < 0 branches, respectively, which are smoothly connected at r = R. We adopt the probe approximation N c N f and the backreaction of the geometry due to the D8 branes are neglected. The S 4 part of the D8 branes has been already integrated out in (3.1). The dual field theory on the boundary of AdS 7 at r → ∞ is a six-dimensional conformal field theory coupled with four-dimensional fermions at the location of the D8/D8 branes. 7 After compactifying the x 11 direction (5.2) and furthermore the τ -direction on a circle τ ∼ τ +2π/M KK with supersymmetry breaking boundary conditions, the theory becomes a four-dimensional, confining SU (N c ) Yang-Mills theory coupled with N f Dirac fermions, that is SU (N c ) QCD with N f flavors, at low energies. Our task is to calculate the induced energy momentum tensor T ij in this boundary theory sourced by the soliton, which corresponds to the nucleon, living on the D8 branes.
While the method of holographic renormalization has been extended to non-conformal theories [56], for our purposes, it is more convenient to work in the eleven-dimensional (or sevendimensional, after reducing on S 4 ) setting (5.1) instead of ten-dimensional type-IIA supergravity with the dilaton. In this setting, we calculate the metric fluctuation caused by the bulk soliton (x M = (x 0,1,2,3 , τ, r, x 11 )) δg M N (τ, x 11 , x, r) = κ 2 7 dτ dx 11 d 4 x dr −G (7) ×G M N AB (τ − τ , x 11 − x 11 , x − x , r, r )T AB (τ , x 11 , x , r ), (5.6) in the axial gauge and study the behavior of the four-dimensional components δg µν (µ, ν = 0, 1, 2, 3) near the boundary r → ∞. The expectation value of the energy momentum tensor is then proportional to the coefficient of the 1/r 4 term integrated over the extra dimensions (See [55] for the precise prescription.) In (5.6), G M N AB is the graviton propagator and G (7) = −(r /L) 10 is the determinant of the AdS part of the metric. T AB (A, B = 0, 1, 2, 3, τ, r, 11) is the seven-dimensional energy momentum tensor (already integrated over S 4 ) of the soliton. κ 2 7 is related to the eleven-dimensional gravitational constant as For a static soliton, T AB does not depend on x 0 nor x 11 . Going to the Fourier space, we get where G( k, r, r ) is the momentum space graviton propagator with k τ = k 11 = k 0 = 0, and we defined T AB ( k, r ) ≡ 2πR 11 dτ T AB (τ , k, r ). (5.11) The τ -integral is trivial because the configuration of D8 branes described above instructs us to write, in the z-coordinates, The bulk energy momentum tensor T AB can be obtained by varying the D8 brane action with respect to the eleven dimensional metric (5.1) Vol(S 4 ), (5.13) where we only keep the quadratic terms in the field strength tensor in the D8 brane action with C = (64π 6 l 5 s ) −1 .g is the induced metric on the D8-branẽ with a nontrivial dilaton field e −Φ = 1 gs (L/r) 3/2 . (In the present case, ∂τ /∂r = 0.) In the elevendimensional notation, (5.14) takes the form This leads to with T µr = ∂z ∂r T µz = 3r 5 zR 6 T µz , T rr = 9r 10 z 2 R 12 T zz . (5.21) [For simplicity only the SU(2) part is shown.] It is important to notice that (5.16) is independent of G τ τ , hence T τ τ = 0 for the soliton configuration. This is simply because the D8/D8 branes do not extend in the τ -direction. On the other hand, T 11,11 is nonvanishing even though the D8/D8 branes do not extend in the x 11 direction. This can be easily understood as the coupling between the gauge fields and the dilaton field in the type IIA description in ten dimensions.

Glueballs in the AdS 7 black hole
In this paper, we do not attempt to compute (5.10) in its full glory because we do not have the exact (numerical) solution of the soliton configuration. Instead, assuming that the latter is known, we will analyze the bulk Einstein equation in detail and establish the connection between the gravitational form factors and the known glueball spectrum of this theory [58,59] (see, also, [60][61][62]).
The relevance of glueballs to the present problem can be understood as follows. Since the spatial part of the field theory energy momentum tensor is transverse k i T ij = 0, near the boundary the metric fluctuation must also be transverse k i δg ij = 0 and can be parameterized as 8 δg TT is the so-called transverse-traceless (TT) part where it is understood that 'trace' is taken in six-dimensions x a = (x µ=0,1,2,3 , x 11 , τ ). As for the trace part, since we neglect the backreaction of the D8/D8 branes in the probe approximation, 9 the geometry is asymptotically AdS with no dilaton. In the axial gauge, the Einstein equation then automatically requires that δg a a ≈ 0 in the asymptotically AdS regime. 10 Therefore, the trace term can be neglected and the D-term entirely originates from the TT modes of the AdS 7 black hole. Glueballs are nothing but the normalizable TT modes of the linearized Einstein equation. They can be classified according to spin under the rotation group SO(3) in physical space and interpreted as the actual spectrum of glueballs on the boundary field theory. Among the 14 independent TT modes, those relevant to the present discussion are referred to as 'T 4 ' and 'S 4 ' in the classification of [59]. 11 The T 4 mode consists of 2 ++ , 1 ++ and 0 ++ glueballs. They are transverse and traceless in five dimensions x 0,1,2,3,11 and have vanishing components in the other dimensions. Despite the difference in spin, all these glueballs obey the same bulk equation of motion, namely the massless Klein-Gordon equation on the AdS 7 black hole. Therefore, their masses are degenerate. The T 4 (2 ++ ) glueballs are transverse-traceless already in four dimensions x 0,1,2,3 . Among the five independent 2 ++ polarizations, the following mode stands out as it features the same spatial tensor structure as in (2.3) where i satisfies k i i = 0. The metric fluctuation corresponding to the T 4 (0 ++ ) glueballs are traceless in five dimensions [58,59] δg T 4 ab (0 ++ ) ∼ δg µν δg µ,11 δg 11,µ δg 11,11 ∼ η µν − kµkν with k µ = (0, k). The spatial part again contains the structure k i k j − δ ij k 2 we look for. Note that the largest component is in the eleventh dimension. Therefore, the T 4 (0 ++ ) mode may be thought of as the counterpart of the dilaton in type IIA supergravity. 8 Near the boundary r → ∞, the metric fluctuation in the d = 6 dimensional subspace x a = (x 0,1,2,3,11 , τ ) can be generically written as in momentum space. Taking k a = δ a i k i and imposing the condition k i δg i j = 0 we arrive at (5.23). 9 The backreacted geometry for the present problem has been calculated in [63]. 10 If one works in the type IIA setup or in AdS 5 /QCD 4 models, one must include this term and solve the Einstein equation coupled to the dilaton field [19,56]. 11 In this paper, we do not consider metric fluctuations on S 4 . Actually, if one does that, there exists another 0 ++ glueball mode called 'L 4 ' [59,60] which can be sourced by the S 4 part of the energy momentum tensor (obtained similarly to T 11,11 in (5.18)). However, this is not a TT mode, and one can show that it does not contribute to the boundary energy momentum tensor.
On the other hand, the S 4 mode consists only of 0 ++ glueballs. Their masses are not degenerate with the T 4 glueballs, and actually the lightest scalar glueball belongs to this class. In [58] this solution was dubbed 'exotic' because it has a nonzero component in the τ -direction Far away from the boundary, the metric fluctuation cannot be written in this simple form, and it is no longer possible to literally maintain the 'transverse-traceless' condition. Besides, the explicit solution constructed in [58] features nonvanishing components in the r-direction (see (5.42) below).
We thus see that there are three different types of glueballs that can potentially contribute to the D-term. Accordingly, the bulk energy momentum tensor (5.17)-(5.18) can be decomposed as 12 The boundary energy momentum tensor T ab of the six-dimensional field theory is therefore given by a linear combination Comparing the four-dimensional part of this to (2.3) or (2.5), we can identify 13 Note that t 2 and t 0 + s 0 both have a pole 1/ k 2 but they cancel in the sum. The components of (5.29) are then given by in agreement with (2.3). Moreover, the traceless condition in six dimensions leads to the QCD trace anomaly relation in four dimensions In order to determine the parameters t 0 , t 2 , s 0 , one needs to fully solve the linearized Einstein equation with the (numerically obtained) bulk energy momentum tensor. This will be discussed in the next subsection. A naive expectation, however, is that the S 4 mode decouples or is suppressed compared to the T 4 modes, namely, |s 0 | |t 0,2 |. Indeed, there have been arguments in the literature [58,62] that the S 4 glueball states may not survive in the 'continuum limit' M KK → ∞ keeping meson masses fixed and including the stringy 1/λ corrections to all orders. The corresponding metric fluctuation (5.26) has the largest component in the τ -direction whose compactification radius 1/M KK shrinks to zero in this limit. Such arguments are also practically motivated since there is an excess of 0 ++ glueballs from holography compared with lattice QCD results [59,62]. We however note that there is no symmetry argument which prevents the S 4 mode from surviving the continuum limit. In this paper, we will not try to include the stringy corrections, but work in the supergravity approximation without taking the continuum limit. Within this approximation, we will see that the S 4 mode actually contributes to the gravitational form factors.

Solving the Einstein equation
Let us investigate in detail whether and how the various glueball states described in the previous subsection couple to the soliton. For this purpose, we turn to the linearized Einstein equation where the covariant derivative and raising/lowering of indices are calculated with respect to the background metric g AB in (5.1). We do not impose δg C C ≈ 0 at this point, since this condition holds only near the boundary r → ∞. After eliminating T C C by taking a trace, (5.35) can be cast into an equivalent form Let us first substitute the following parametrization relevant to the T 4 (2 ++ ) mode (see (5.24)), we find, for k µ = (0, k), with all the other components vanishing. Here, ∇ 2 is the massless Klein-Gordon operator 39) and the propagation is diagonal in Lorentz indices. Naturally, this mode can be sourced by a part of the bulk T µν that has the same tensorial structure as in (5.38). As for the T 4 (0 ++ ) mode, we use (see (5.25)) 40) and find that Again there is no mixing of Lorentz indices. This explains the above-mentioned degeneracy between the T 4 (2 ++ ) and T 4 (0 ++ ) glueballs. Turning to T AB on the right hand side, we see from (5.18) that the component T 11,11 is nonvanishing. Therefore, the soliton can act as a source for the T 4 (0 ++ ) mode (or the dilaton in the type IIA language).
(5. 45) We see that the spatial components H ij have the expected tensor structure. On the other hand, unlike for the T 4 modes above, now the τ τ component H τ τ ∝ T τ τ is nonvanishing. Since T τ τ = 0 for the soliton configuration as we pointed out below (5.21), naively this implies that the soliton does not excite the S 4 mode. However, the situation is more complex because of the presence of T other in (5.27) and its induced metric fluctuation δg other . The primary role of this term is to account for the T µr and T rr components of the bulk energy momentum tensor. These are nonvanishing for the soliton solution (see (5.19)-(5.21)), but they do not directly couple to any of the above TT modes. In the axial gauge δg M r = 0, the components of the Einstein equation having an r-index serve as first-order constraints [64] that can be used to eliminate non-TT modes via the conservation law ∇ M T M n = 0 (n = 0, 1, 2, 3, 11), or explicitly, Let us employ the following ansatz δg other µν ( k, r) = r 2 L 2 k µ k ν k 2 a(r, k), δg other 11,11 ( k, r) = (r, k)). This solves the linearized Einstein equation with the bulk energy momentum tensor of the form:  where the right hand side of these equations is the bulk energy momentum tensor (5.21). Because the gauge configuration under consideration is spherically symmetric, one can show that T µr is proportional to k µ and hence the first equation of (5.54) can be solved immediately. The second equation of (5.54) contains at most first order derivatives in r and it can be uniquely solved by imposing the boundary condition a → 0 at r → ∞. With the conditions (5.54), we find that the metric fluctuation (5.48) satisfies the µr and rr components of the linearized Einstein equation (5.36). Note that the above parameterization (5.48) is not unique. Different parameterizations lead to different decompositions of the Einstein equation (5.27) and (5.28) without changing the total induced metric δg ab and hence the boundary energy momentum tensor T ab . We are led to the choice (5.48) because then the reduced 5D energy momentum tensor is a total derivative in z, 15 a feature that will turn out to be convenient shortly. As r goes to infinity, T other rr = T rr ∼ 1/r 11 and T other µr = T µr ∼ 1/r 10 as can be seen by substituting (3.6) into (5.19) and (5.20) and noticing that G ∼ 1/z ∼ 1/r 3 and H ∼ 1/z 2 ∼ 1/r 6 in this limit. This implies a ∼ 1/r 9 and b ∼ 1/r 3 and hence δg other µν ∼ 1/r 7 and δg other 11,11 ∼ δg other τ τ ∼ 1/r. We see that δg other µν decays too fast to contribute to T µν , but g other 11,11 and δg other τ τ decay slowly and potentially contribute to T τ τ and T 11,11 16 (but not to T µν ). On the other hand, near z = 0 or r = R where the instanton approximation is good, we can use (3.5) to deduce the singular behavior as z → 0 The leading singularity comes from the U(1) part.) Comparing with (5.51) and (5.52), we find that Since the τ τ component of the total bulk energy momentum tensor T τ τ , as well as the T 4 part T T(2) τ τ and T T(0) τ τ , is zero, (5.53) means that an effective source term for the S 4 mode is induced. To compensate (5.53) by the S 4 mode, we impose (see (5.43)), which implies with H S ∼ 1/r 3 and h S ∼ 1/r 3 as r → ∞. Comparing the large-r behavior of (5.58) with that of (5.44), we can identify We thus conclude that, contrary to the naive expectation, the S 4 glueballs in general contribute to the boundary energy momentum tensor, albeit in a somewhat indirect way. Consequently, the source term for the T 4 glueballs is not the total T µν (5.17), but instead 17 Note that the leading term δg S µν ∼ r 2 h S ∼ 1/r as r → ∞ is canceled against a like term in δg T µν ∝ T T µν ∼ −T S µν . Indeed, from (5.43) we find Similarly, from (5.37) and (5.40) Therefore, near the boundary r → ∞, Since T T µν + T S µν = T µν − T other µν ∼ 1/r 7 , most of the terms in δg T/S µν which decay slower than 1/r 7 cancel in the sum δg T µν + δg S µν . The leading surviving contribution is the 1/r 4 term This is precisely what contributes to the boundary energy momentum tensor via holographic renormalization (5.8). Remarkably, when k = 0, C µν is directly proportional to the classical energy momentum tensor (3.14) as can be seen by writing at k = 0. In the last equality, we have used the identity (5.45) in order to recast the integrand (cf. (5.59)) as a total derivative. The first term on the left hand side is proportional to (3.14) and the second term vanishes (cf. (5.57)) This means that the naive expression (3.14) can be used to calculate gravitational form factors at k = 0. Below we shall reconfirm this important observation in a different way.

Glueball dominance of the gravitational form factors
Armed with the general discussion in the previous subsection, we are now ready to compute the form factor D( k) for generic values of k. From (5.33), we can write where each component can be determined by solving the corresponding wave equations. It should be clear by now that the split D T ∼ t 2 + t 0 is actually not necessary. As we have seen above, the propagators of the T 4 (2 ++ ) and T 4 (0 ++ ) modes are identical and diagonal in Lorentz indices whereG T is the Green function for the massless Klein-Gordon equation We can thus write, in (5.10), where T T µν is given by (5.63). This leads to the formula where we used (5.17) and in the second line we switched to the variable z using (5.5). The omitted terms include the subtraction −T S ij −T other ij in (5.63). Similarly, with the help of the identity (5.45), we introduce the Green function for the S 4 mode up to terms irrelevant (suppressed by powers of 1/r) to holographic renormalization.
There is, however, a subtlety in the above discussion. As we have pointed out in the previous subsection, both δg T µν and δg S µν have a component which decays slowly ∼ 1/r. While they cancel eventually in the sum δg T µν + δg S µν , it is convenient to extract the O(1/r 4 ) components in δg T/S µν by subtracting the slowly decaying terms, because each term in our formal expression of the Green functionsG T/S (r, r ) in (5.83) behaves as 1/r 6 as r → ∞ and the 1/r terms in δg T/S µν cannot be captured without taking the infinite sum. For this purpose, we define and the corresponding subtracted bulk energy momentum tensor With this in mind, we now solve (5.72) and (5.75) by introducing the complete set of eigenfunctions − 1 r 3 L 4 ∂ r (r 7 − rR 6 )∂ r Ψ T n (r) = (m T n ) 2 Ψ T n (r), (n = 1, 2, 3, · · · ) (5.79) (5.80) which are regular at r = R and vanish in the limit r → ∞. They satisfy the completeness relations The solution of (5.72) and (5.75) can then be formally written as (5.83) m T/S n are nothing but the masses of the T 4 and S 4 0 ++ glueballs [59]. Numerically (see Table 1 of [59]), It is easy to see that the normalizable modes Ψ We can thus write the leading term as The relation between T µν and C T/S µν is fixed by holographic renormalization [55,56]. This is conveniently done in a different coordinate system = L 4 in which the metric (5.1) takes the Fefferman-Graham form where KK = L 4 /R 2 . From this expression we can read off the boundary energy momentum tensor To get the four-dimensional energy momentum tensor, we have to integrate over the extra dimensions x 11 and τ , but this has been already done in (5.10).
We have thus seen that the D-term gravitational form factor is saturated by the exchange of an infinite tower of 2 ++ and 0 ++ glueballs. Schematically, This should be contrasted with the tripole form The summation over n cannot be performed in a closed analytic form. However, the value at the special point k = 0 can be exactly determined with the help of (5.89) lim r→∞G T/S (r, r , k = 0) = 1 Remarkably, this is independent of r , so the convolution with the graviton propagator disappears. Thus, the distinction between the S 4 and T 4 propagators become irrelevant at this point. It then follows that This is nothing but the naive classical formula (3.19) evaluated in the previous section except for the subtraction term −T other in (5.63) which, however, can be ignored because it is a total derivative (5.69). We have thus arrived at the same conclusion as in the previous subsection that the value D(0) can be calculated by simply Fourier transforming the classical energy momentum tensor (3.19) without any reference to glueballs. At k = 0, exchanging infinitely many glueballs is tantamount to exchanging no glueball at all.
Another regime where an analytical expression of the propagator is available is the large momentum region | k| M KK . 18 In this regime, we expect thatG T/S (r → ∞, r , k) is exponentially small in | k| for | k|L 2 > r . 19 In the opposite region r > | k|L 2 = 3R | k| M KK R, (5.101) 18 The following discussion should be taken with caution because the present model starts to deviate from QCD as the momentum transfer gets much larger than M KK . 19 When | k|L 2 r R, one can see from (5.102) that the propagator is exponentially suppressed, and the suppression gets stronger as r approaches R from above. When r ∼ R, the AdS approximation breaks down. Nevertheless, by formally Taylor expanding (5.83) the geometry is approximately AdS 7 (instead of AdS 7 black hole), and the propagator is explicitly known in terms of the modified Bessel functions [68] G T/S (r, r , k) ≈ L 7 r 3 r 3 K 3 (| k|L 2 /r )I 3 (| k|L 2 /r) ≈ | k| 3 L 13 48r 6 r 3 K 3 (| k|L 2 /r ), (5.102) for r → ∞ > r . (5.102) reduces to (5.98) by formally setting k = 0. Noting that r ≈ |z| 1/3 R in this regime, we see that the D-term can be calculated from the convolution Since the z-integral is effectively limited to |z| 27| k| 3 /M 3 KK , we see that the gravity effect strongly suppresses the D-term at large-| k|.
To summarize, we arrive at the following strategy to calculate the energy momentum tensor and the gravitational form factors. The value at the special point k = 0 can be obtained 'classically' as demonstrated in the previous section. When 0 < | k| M KK ∼ 1 GeV, they can be calculated, in principle, by (5.93) and (5.94) where one has to explicitly evaluate the Fourier modes of the bulk energy momentum tensor T µν (r, k) and perform the infinite sum. Eq. (5.94) exhibits the phenomenon of 'glueball dominance', namely the form factor is saturated by the exchange of the T 4 and S 4 glueballs. When | k| M KK , we may continue to use holographic renormalization. The bulk spacetime effectively reduces to AdS 7 (as opposed to the AdS 7 black hole) and the connection to glueballs becomes less obvious. Note that the first region | k| 1 GeV practically covers the phenomenologically important region [11,12,15]. The region | k| 1 GeV is also interesting for certain applications [13,14,21], but in this regime one should use perturbative QCD techniques [66,67] rather than holography.
Finally, let us comment on the relation to previous approaches in the literature [11,19,36]. In bottom-up holographic models, one works in asymptotically AdS 5 spaces. Needless to say, the S 4 mode is nonexistent in these spaces. There are five transverse-traceless graviton polarizations in AdS 5 , in one-to-one correspondence with the T 4 (2 ++ ) modes in AdS 7 . They are interpreted as the 2 ++ glueballs in the boundary field theory. Via holographic renormalization, the AdS 5 counterpart of (5.24) induces the transverse-traceless part of the QCD energy momentum tensor, namely, the first line of (2.5) [11]. Separately to this, one has to exchange the dilaton dual to the 0 ++ glueballs which leads to the trace part of the energy momentum tensor, the second line of (2.5). The A-form factor in the coefficient of k µ k ν − η µν k 2 cancels between the two, and only the D-form factor remains. However, the split (2.5) is somewhat artificial, and the cancellation may become a delicate issue [36]. In our AdS 7 /CFT 6 /QCD 4 setup, the 2 ++ and 0 ++ glueballs (the latter would correspond to the dilaton in a five-dimensional setup) are treated on an equal footing. They are both transverse-traceless from seven-dimensional point of view, and their distinction is immaterial for the purpose of calculating the gravitational form factors.

Discussions and conclusions
Let us try to interpret our results diagrammatically. In our model, the form factor D(t = − k 2 ) is given by the convolution of the soliton energy momentum tensor T AB and the graviton propagator in the AdS 7 black hole geometry which contains information about the 2 ++ and 0 ++ glueballs. When t = 0, the graviton propagator trivializes. The D-term D(0) is then related to the Fourier transform of the 'classical' energy momentum tensor T cl (3.14) evaluated at on-shell. Since T cl ∼ G 2 , H 2 (see (4.13)) can be expressed by the square of vector meson propagators near the boundary, one may say that a graviton couples to the nucleon by exchanging pairs of pions and infinitely many vector and axial-vector mesons. 20 This is depicted in the left diagram of Fig. 1. Roughly the D-term may be understood as the coupling between a graviton and a 'meson cloud', similarly to the pion cloud in the chiral quark soliton model (see, e.g., [38]). We however emphasize that in our calculation the pion contribution is just one term n = 0 in the infinite sum (4.10). When 0 < | k| M KK ∼ 1 GeV, one sees the individual contributions of glueballs (5.83). The emerging picture in this regime is that the meson pairs from T AB in the bulk couple with glueball intermediate states which eventually interact with an external graviton (Fig. 1, middle). One may interpret this situation as the 'glueball dominance' of the gravitational form factors, in perfect analogy to the vector meson dominance of the electromagnetic form factors previously observed in our model [32]. In both cases, one has to exchange infinitely many resonances. Curiously, as we pointed out below (4.11), the masses of the T 4 glueballs are strikingly close to those of two vector mesons. To our knowledge, this coincidence does not seem to have been pointed out in the literature. It would be very interesting to explore the reason and implications of this observation. Finally, in the large-k region, the geometry becomes effectively AdS 7 (as opposed to the AdS 7 black hole) and the information about the glueball spectrum is lost. While the meson exchange picture is still possible, in this region one should rather use perturbative QCD techniques to calculate form factors. In particular, the D-term is given by the two gluon exchange to leading order (Fig. 1, right) [66].
The connection between the D-term and glueballs has been recently emphasized by Mamo and Zahed using a bottom-up AdS 5 /QCD 4 model [18,36]. As we commented at the end of Section 5, in their calculation the D-term results from a cancellation between the transverse-traceless (TT) graviton and dilaton exchanges in AdS 5 (in other words, between the first and second lines in (2.5)). They argue that, due to the degeneracy of the 2 ++ and 0 ++ glueball masses, the cancellation is complete (meaning D = 0) in the large-N c limit. They further argue that the D-term becomes nonzero only after including 1/N c corrections, suggesting that the large-N c counting argument [2] D ∼ O(N 2 c ), (6.1) does not hold. In contrast, in our AdS 7 /CFT 6 /QCD 4 calculation, the contributions from the 2 ++ and 0 ++ T 4 glueballs add up (see (5.71), (5.73)), after which their distinction becomes immaterial. Moreover, the value D(k = 0) can be calculated 'classically' as in Section 4 without any reference to glueballs. Our result is consistent with the large-N c counting since (M KK ∼ O(N 0 c ) since it is related to the ρ-meson mass.) Interestingly, however, the actual numerical value calculated in Section 4 is small due to a cancellation between the iso-singlet and iso-vector contributions.
In conclusion, we have demonstrated that our model of QCD offers a novel holographic framework to calculate the nucleon D-term and possibly other gravitational form factors, with a vivid physical interpretation in terms of meson and glueball exchanges. There are a number of avenues for improvement and generalization. First, the quantization of the collective coordinates should be implemented. It is known that the quantization brings in uncertainties in baryon masses associated with the zero-point (vacuum) energy of harmonic oscillators [30]. We expect that the D-term, being a nonforward matrix element, is less affected by this problem. However, without quantization, the spin-1/2 nature of the nucleon is lost, and one cannot talk about the B-form factor (2.2). Second, in order to calculate D(k) for nonvanishing values of k, it is desired to convolute a more accurate soliton configuration in the entire space 0 < |z| < ∞ obtained numerically (see e.g., [43][44][45]) with the graviton propagator to incorporate the contributions from the glueballs. This is crucial, in particular, for the calculation of the 'mechanical radius' [2] 21 .

(6.3)
Since the entire region in k is involved, a proper evaluation of the mechanical radius must include the glueball degrees of freedom. The calculation is significantly more complex than that of the electromagnetic form factors and the associated 'charge radius' in this model [32] which only requires the knowledge of the asymptotic behavior z → ∞ of the bulk gauge fields. There are other interesting directions such as the effect of finite quark masses (finite pion mass) and generalizations to mesons and excited baryons [70], and perhaps even to atomic nuclei [71]. We hope to address these issues in future.