Abstract

We study the energy-momentum tensor of a baryon in a top-down holographic QCD. In holographic QCD, the baryons are represented as solitons in a 5-dimensional gauge theory. We obtain the soliton solution by solving the equations of motion numerically. Using this result, the energy-momentum tensor and related quantities such as the mass, mean square radii, and the D-term (druck term) are computed. The evaluated D-term is about $-2.05$, whose absolute value is significantly larger than that in the previous work [Fujita:2022jus].

1 Introduction

The gravitational form factors are important quantities to probe the internal structures of hadrons. They are defined by the matrix elements of the energy-momentum tensor (EMT) and carry the information of the distributions of mass, spin, pressure, and shear force inside the hadron. Recently, these quantities have been extracted from experimental data [Kumano:2017lhr, Burkert:2018bqq, Kumericki:2019ddg, Burkert:2021ith, Kou:2021qdc, Wang:2022ndz] and triggered a lot of theoretical developments. (See [Polyakov:2019lbq, Burkert:2023wzr] and the references therein.) In general, the evaluation of hadronic matrix elements in quantum chromodynamics (QCD), especially at low energies, is complicated because it requires serious analysis in a strongly coupled system involving the physics of bound states. Therefore, it would be worthwhile to investigate various methods to calculate these quantities and compare the predictions with each other and with experimental data.

The main goal of this paper is to report an improved numerical analysis of the gravitational form factors, in particular the so-called D-term for the nucleon, using a top down holographic QCD proposed in [Sakai:2004cn, Sakai:2005yt]. The holographic description of QCD is obtained from a configuration of D-branes that realizes QCD in string theory. Since string theory contains gravity, it provides a natural framework to read off the gravitational form factors through the gravitational interactions of the hadrons.

The explicit procedure to obtain the gravitational form factors for baryons in this model was investigated in [Fujita:2022jus] in detail and the value of the D-term was estimated. It was pointed out that the gravitational form factors (in the leading order of the $1/N_{c}$ expansion) can be expressed as a sum over the contributions of glueball exchange diagrams, which was referred to as the glueball dominance. Furthermore, it was shown that in the zero momentum transfer limit, the gravitational form factors can be estimated by the EMT obtained in the meson effective theory. Based on these results, the value of the D-term was estimated to be around $-0.14$. However, the analysis in [Fujita:2022jus] cannot be considered complete. One of the drawbacks was that the soliton configuration was obtained not by solving the equations of motion (EOMs), but by smoothly connecting the solutions near the center and the boundary. In this paper, we improve this point by numerically solving the EOMs. It turns out that the value of the D-term is sensitive to the configuration of the gauge field in the intermediate region and it is crucial to accurately obtain the solution there. Our prediction of the value of the D-term is around $-2.05$. We also compute mean square radii, energy density, pressure, and shear force.

The organization of the paper is as follows. We start by reviewing the gravitational form factor in section 2 and summarizing useful formulas. In section 3, we investigate the EMT obtained from the meson effective action in holographic QCD for the soliton solution under Witten’s ansatz [Witten:1976ck]. Our numerical results are summarized in section 4. The estimated values of the D-term and mean square radii are listed in Table 1. Section 5 is devoted to the conclusion and discussion of possible future directions.

2 Gravitational form factor

2.1 Definition of the gravitational form factors and the D-term

In this subsection, we briefly review the definition of the gravitational form factors (GFFs) and the D-term based on [Polyakov:2018zvc]. The GFFs are given by the matrix element of the EMT operator $\hat{T}_{\mu\nu}$. According to Lorentz invariance, it is written as

$$\begin{split}&\matrixelement{p^{\prime},s^{\prime}}{\hat{T}_{\mu\nu}(x)}{p,s}\\ &=\bar{u}\quantity[A(t)\frac{P_{\mu}P_{\nu}}{m}+J(t)\frac{i(P_{\mu}\sigma_{\nu\rho}+P_{\nu}\sigma_{\mu\rho})\Delta^{\rho}}{2m}+D(t)\frac{\Delta_{\mu}\Delta_{\nu}-\eta_{\mu\nu}\Delta^{2}}{4m}+\bar{c}(t)\eta_{\mu\nu}]ue^{-i\Delta x}.\end{split}$$ (2.1)

for a spin-1/2 particle with mass $m$. Here $\ket{p,s}$ is a one-particle state with momentum $p$ and helicity $s$, which satisfies the normalization condition $\innerproduct{p^{\prime},s^{\prime}}{p,s}=2p^{0}(2\pi)^{3}\delta({\bf\it p}^{\prime}-{\bf\it p})\delta_{s^{\prime},s}$. $P$, $\Delta$, and $t$ are defined as $P:=(p+p^{\prime})/2$, $\Delta:=p^{\prime}-p$, and $t:=-\Delta^{2}$, respectively. $u(p,s)$ is the four-component spinor normalized as $\bar{u}(p,s)u(p,s)=2m$. $A(t)$, $J(t)$, and $D(t)$ are the GFFs which are renormalization scale invariant scalar functions representing the inner structure of the particle. $\bar{c}(t)$ corresponds to the breaking of EMT conservation and it vanishes for the total EMT. It appears when the EMT is separated into some parts (in this paper, we will separate the EMT into the $SU(N_{f})$ part and the $U(1)$ part). The number of form factors depends on the spin of the particle. For example, the $J$ term is absent for scalar particles.

It is convenient to introduce the Breit frame $P=(E,{\bf\it 0})$, $\Delta=(0,{\bf\it\Delta})$. In the Breit frame, the static EMT matrix elements are expressed as

$$\matrixelement{p^{\prime},s^{\prime}}{\hat{T}^{00}(0)}{p,s} =2mE\quantity[A(t)-\frac{t}{4m^{2}}\quantity(A(t)-2J(t)+D(t))]\delta_{s^{\prime}s},$$ (2.2)

$$\matrixelement{p^{\prime},s^{\prime}}{\hat{T}^{0i}(0)}{p,s} =2mE\quantity[J(t)\frac{(-i{\bf\it\Delta}\times{\bf\it\sigma}_{s^{\prime}s})_{i}}{2m}],$$ (2.3)

$$\matrixelement{p^{\prime},s^{\prime}}{\hat{T}^{ij}(0)}{p,s} =2mE\quantity[D(t)\frac{{\bf\it\Delta}_{i}{\bf\it\Delta}_{j}-\delta_{ij}{\bf\it\Delta}^{2}}{4m^{2}}]\delta_{s^{\prime}s}.$$ (2.4)

Here, $\sigma_{s^{\prime}s}^{i}=\chi^{\dagger}_{s^{\prime}}\sigma^{i}\chi_{s}$ with the Pauli spinors $\chi_{s^{\prime}}$ and $\chi_{s}$ in the respective rest frames normalized as $\chi^{\dagger}_{s^{\prime}}\chi_{s}=\delta_{s^{\prime}s}$. Using the Fourier transformation

$$T_{\mu\nu}({\bf\it x}):=\int\frac{\differential^{3}{\bf\it\Delta}}{(2\pi)^{3}2E}\matrixelement{p^{\prime},s}{\hat{T}_{\mu\nu}(0)}{p,s}e^{-i{\bf\it\Delta}\cdot{\bf\it x}},$$ (2.5)

$A(t)$, $J(t)$, and $D(t)$ satisfy

$$A(t)-\frac{t}{4m}\quantity(A(t)-2J(t)+D(t)) =\frac{1}{m}\int\differential^{3}{\bf\it x}T^{00}({\bf\it x})e^{i{\bf\it\Delta}\cdot{\bf\it x}},$$ (2.6)

$$\quantity(\delta_{ij}+\quantity(t\delta_{ij}+{\bf\it\Delta}_{i}{\bf\it\Delta}_{j}){\derivative{t}})J(t)s^{j} =\int\differential^{3}{\bf\it x}\epsilon_{ijk}x^{j}T^{0k}({\bf\it x})e^{i{\bf\it\Delta}\cdot{\bf\it x}},$$ (2.7)

$$\quantity(1+\frac{4}{3}t{\derivative{t}}+\frac{4}{15}t^{2}\derivative[2]{t})D(t) =-\frac{2}{5}m\int\differential^{3}{\bf\it x}\quantity({\bf\it x}^{i}{\bf\it x}^{j}-\frac{1}{3}\delta^{ij}{\bf\it x}^{2})T_{ij}({\bf\it x})e^{i{\bf\it\Delta}\cdot{\bf\it x}}.$$ (2.8)

Here, $s$ dependence of $T^{\mu\nu}({\bf\it x})$ is implicit and $s^{i}=(\sigma_{ss})^{i}$ in (2.7) is the unit vector along the spin direction.

In the forward limit $t\to 0$, we have

$$A(0) =\frac{1}{m}\int\differential^{3}{\bf\it x}T^{00}({\bf\it x})=1,$$ (2.9)

$$J(0)s_{i} =\int\differential^{3}{\bf\it x}\epsilon_{ijk}x^{j}T^{0k}({\bf\it x})=\frac{1}{2}s_{i},$$ (2.10)

$$D(0) =-\frac{2}{5}m\int\differential^{3}{\bf\it x}\quantity({\bf\it x}^{i}{\bf\it x}^{j}-\frac{1}{3}\delta^{ij}{\bf\it x}^{2})T_{ij}({\bf\it x}).$$ (2.11)

The value of $A(0)$ is fixed to 1. This corresponds to the mass of the particle being $m$. In the same way, $J(0)$ is fixed to 1/2, which reflects the fact that the particle has spin 1/2. On the other hand, the value of $D=D(0)$ is not determined in general. This $D$ is called the D term and reflects the internal force distribution of the particle. The D-term is as fundamental as the mass and spin of a particle and is often referred to as “the last unknown global property” [Polyakov:2018zvc]. One of our goals is to estimate the value of the D-term of the baryon.

The value of the D-term for the nucleon has been evaluated using experimental data [Burkert:2018bqq, Kumericki:2019ddg, Burkert:2021ith, Kou:2021qdc, Wang:2022ndz, Goharipour:2025lep], bag model [Neubelt:2019sou], lattice QCD [Shanahan:2018nnv, Shanahan:2018pib, Hackett:2023rif], light-cone QCD sum rules [Anikin:2019kwi, Azizi:2019ytx, Dehghan:2025ncw], Skyrme model and its generalization [Cebulla:2007ei, Jung:2013bya, Kim:2012ts, GarciaMartin-Caro:2023klo, GarciaMartin-Caro:2023toa], chiral quark soliton model [Goeke:2007fp, Wakamatsu:2007uc], and bottom up holographic QCD [Mamo:2019mka, Mamo:2021krl, Mamo:2022eui]. It is expected to be in the range $-5<D<-1$. For details, see [Polyakov:2018zvc] and [Cao:2024zlf].

2.2 Spherically symmetric case

In a spherically symmetric system, the stress tensor $T_{ij}$ can be decomposed into a trace part and a traceless part

$$T_{ij}({\bf\it x})=p(r)\delta_{ij}+s(r)\quantity(n_{i}n_{j}-\frac{1}{3}\delta_{ij}),$$ (2.12)

where $p(r)$ is the pressure and $s(r)$ is the shear force inside the particle. Substituting the spherically symmetric stress tensor (2.12) into (2.8) and solving the differential equation, the form factor $D(t)$ can be expressed as

$$D(t)=16\pi m\int\differential rr^{4}s(r)\frac{j_{2}\quantity(\sqrt{-t}r)}{tr^{2}},$$ (2.13)

where $j_{2}(x)$ is the spherical Bessel function. Taking the forward limit $t\to 0$, we have

$$D=-\frac{16\pi}{15}m\int\differential rr^{4}s(r).$$ (2.14)

The conservation law $\partial_{\mu}T^{\mu\nu}=0$ yields a constraint between $p$ and $s$,

$$p^{\prime}(r)+\frac{2}{3}s^{\prime}(r)+\frac{2}{r}s(r)=0.$$ (2.15)

From this equation, it can be shown that the pressure $p$ satisfies the von-Laue condition [laue1911dynamik]

$$\int\differential rr^{2}p(r)=0.$$ (2.16)

This implies that the pressure has positive and negative regions. Using (2.15) and integration by parts, we can express the D-term in terms of $p(r)$ instead of $s(r)$,

$$D=4\pi m\int\differential rr^{4}p(r).$$ (2.17)

Note that (2.14) and (2.17) are equivalent only for the total EMT. For the separated EMT ($SU(N_{f})$ part and $U(1)$ part, for example), the non-conserving term $\bar{c}$ also appears on the right-hand side of (2.17).

3 The energy-momentum tensor in holographic QCD

3.1 The model

The model we consider is the top-down holographic QCD model introduced in [Sakai:2004cn, Sakai:2005yt]. It is a holographic dual of $SU(N_{c})$ QCD with $N_{f}$ massless quarks realized in a system of $N_{c}$ D4-branes compactified on $S^{1}$ of radius $M_{\mathrm{KK}}^{-1}$ with $N_{f}$ pairs of D8-$\mathrm{\overline{D8}}$ branes. The holographic dual description is obtained by replacing the D4-brane with the corresponding supergravity background [Witten:1998zw]. This gives a system of $N_{f}$ D8-branes embedded in this background, and closed and open strings correspond to glueballs and mesons, respectively.

The low energy effective action of open strings on the D8 branes is the $U(N_{f})$ Yang-Mills theory with the Chern-Simons term on a 1+4 dimensional curved background spacetime,

$$S=-\frac{\kappa}{2}\int\differential^{4}x\differential z\tr\quantity(h(z)\mathcal{F}_{\mu\nu}^{2}+2k(z)\mathcal{F}_{z\mu}^{2})+\frac{N_{c}}{24\pi^{2}}\int\omega_{5}(\mathcal{A}).$$ (3.1)

We take the Kaluza-Klein mass $M_{\mathrm{KK}}=1$ for simplicity. Here $\mu=0,1,2,3$ are 1+3 dimensional Lorentz indices, and $z$ is the coordinate of the fifth space-like dimension. $\mathcal{A}=\mathcal{A}_{\mu}\differential x^{\mu}+\mathcal{A}_{z}\differential z$ is the $U(N_{f})$ gauge field and $\mathcal{F}:=\differential\mathcal{A}+i\mathcal{A}\wedge\mathcal{A}$ is its field strength. We can decompose $\mathcal{A}$ into $SU(N_{f})$ part and $U(1)$ part,

$$\mathcal{A}=\mathcal{A}^{SU(N_{f})}+\frac{1}{\sqrt{2N_{f}}}1_{N_{f}}\mathcal{A}^{U(1)}.$$ (3.2)

$\omega_{5}(\mathcal{A})$ is the Chern-Simons 5-form of $\mathcal{A}$. The constant $\kappa$ is related to the ’t Hooft coupling $\lambda$ (at the scale of $M_{\rm KK}$) as

$$\kappa=\frac{N_{c}\lambda}{216\pi^{3}}.$$ (3.3)

$h(z)$ and $k(z)$ are given by

$$h(z)=(1+z^{2})^{-\frac{1}{3}},\,k(z)=1+z^{2}.$$ (3.4)

The EOM for the action (3.1) is

$$\begin{split}-2\kappa\quantity(\mathcal{D}_{\nu}(h(z){\mathcal{F}}^{\nu\mu})+\mathcal{D}_{z}\quantity(k(z)\mathcal{F}^{z\mu}))&=\frac{N_{c}}{16\pi^{2}}\epsilon^{\mu\nu\rho\sigma z}({\mathcal{F}}_{\rho\sigma}\mathcal{F}_{\nu z}+\mathcal{F}_{\nu z}\mathcal{F}_{\rho\sigma})\\ -2\kappa\mathcal{D}_{\mu}(k(z)\mathcal{F}^{\mu z})&=\frac{N_{c}}{32\pi^{2}}\epsilon^{\mu\nu\rho\sigma z}\mathcal{F}_{\mu\nu}\mathcal{F}_{\rho\sigma}\end{split}$$ (3.5)

where $\mathcal{D}_{\mu}=\partial_{\mu}+i\commutator{\mathcal{A}_{\mu}}{\cdot}$ and $\mathcal{D}_{z}=\partial_{z}+i\commutator{\mathcal{A}_{z}}{\cdot}$ are the covariant derivatives, and $\epsilon_{\mu\nu\rho\sigma}$ is a 4-dimensional Levi-Civita symbol with $\epsilon_{0123z}=+1$. In this model, the fluctuations of the $U(N_{f})$ gauge field represent massless pions and massive (axial) vector mesons, while baryons are obtained as topological solitons in 1+4 dimensional spacetime [Sakai:2004cn, Sakai:2005yt, Hata:2007mb]. The soliton corresponding to a baryon is an instanton on the 4-dimensional space which carries a non-trivial topological number given by the second Chern number

$$B=\frac{1}{8\pi^{2}}\int\differential^{3}{\bf\it x}\differential z\epsilon^{0ijkz}\tr(\mathcal{F}_{ij}\mathcal{F}_{kz})=-\frac{1}{8\pi^{2}}\int\differential^{3}{\bf\it x}\differential z\epsilon_{0ijkz}\tr(\mathcal{F}_{ij}\mathcal{F}_{kz}).$$ (3.6)

Here ${\bf\it x}=(x^{1},x^{2},x^{3})$ and $i,j,k=1,2,3$. This is related to the baryon number in the Skyrme model.

3.2 Soliton solution

We take $N_{f}=2$ in this paper for simplicity. In order to construct a soliton solution with baryon number $B=1$, we set Witten’s ansatz [Witten:1976ck]

$$\mathcal{A}_{z}^{SU(2)} =A_{z}(r,z)n^{a}\tau_{a}$$ (3.7)

$$\mathcal{A}_{i}^{SU(2)} =A_{r}(r,z)n_{i}n^{a}\tau_{a}+\frac{\Phi_{1}(r,z)}{r}\quantity(\tau_{i}-n_{i}n^{a}\tau_{a})-\frac{\Phi_{2}(r,z)+1}{r}\epsilon_{ija}n_{j}\tau_{a}$$ (3.8)

$$\mathcal{A}_{0}^{U(1)} =a_{0}(r,z),$$ (3.9)

$$\mathcal{A}_{0}^{SU(2)} =\mathcal{A}_{z}^{U(1)}=\mathcal{A}_{i}^{U(1)}=0,$$ (3.10)

where $r=|{\bf\it x}|$ is the radial coordinate and ${\bf\it n}:={\bf\it x}/r$ is the unit vector of the 3-dimensional space. $\tau_{a}:=\sigma_{a}/2,\,(a=1,2,3)$ are the generators of $SU(2)$ Lie algebra. This ansatz is invariant under $SO(3)$ “spatial rotation” combined with the global $SU(2)$ transformation

$$\begin{split}\mathcal{A}_{z}^{SU(2)}({\bf\it x},z)&=U(R)\mathcal{A}_{z}^{SU(2)}(R^{-1}{\bf\it x},z)U(R)^{\dagger},\\ \mathcal{A}_{i}^{SU(2)}({\bf\it x},z)&=U(R){R_{i}}^{j}\mathcal{A}_{j}^{SU(2)}(R^{-1}{\bf\it x},z)U(R)^{\dagger},\end{split}\quad R\in SO(3)$$ (3.11)

where $U(R)$ is the $SU(2)$ representation of $R$. Notice that this ansatz has a residual $U(1)$ gauge transformation

$$\mathcal{A}_{M}(x)^{SU(2)}\to e^{i\lambda(r,z){\bf\it n}\cdot{\bf\it\tau}}\quantity(\mathcal{A}_{M}^{SU(2)}(x)-i\partial_{M})e^{-i\lambda(r,z){\bf\it n}\cdot{\bf\it\tau}}\quad(M=0,1,2,3,z).$$ (3.12)

Under this transformation, $\Phi:=\Phi_{1}+i\Phi_{2}$, $A_{r}$, and $A_{z}$ transform as

$$\begin{split}\Phi\to e^{-i\lambda}\Phi,\,\,A_{\alpha}\to A_{\alpha}-\partial_{\alpha}\lambda\quad(\alpha=r,z).\end{split}$$ (3.13)

Therefore, $(A_{r},A_{z})$ and $\Phi$ can be considered as a $U(1)$ gauge field and a $U(1)$ charged scalar, respectively, in the $r$ - $z$ plane. It is convenient to introduce the $U(1)$ field strength and covariant derivative on the $r$ - $z$ plane

$$F_{\alpha\beta}:=\partial_{\alpha}A_{\beta}-\partial_{\alpha}A_{\beta},\,D_{\alpha}\Phi:=(\partial_{\alpha}-iA_{\alpha})\Phi.$$ (3.14)

and the 1+4 dimensional Yang-Mills theory is reduced to 2 dimensional Abelian-Higgs theory. The Chern number (3.6) is written as the winding number of the vortex

$$B=-\frac{1}{2\pi}\int\differential r\differential z\quantity(F_{rz}(1-\absolutevalue{\Phi}^{2})+2\imaginary(D_{r}\Phi^{\dagger}D_{z}\Phi))=:\int\differential r\differential z\mathrm{ch}_{2}(r,z),$$ (3.15)

which only depends on the boundary values of $\Phi$ and $A_{\alpha}$ as

$$B=\int\differential r\quantity[A_{r}+\imaginary(\Phi^{\dagger}D_{r}\Phi)]_{z=-\infty}^{z=\infty}-\int_{-\infty}^{\infty}\differential z\quantity[A_{z}+\imaginary(\Phi^{\dagger}D_{z}\Phi)]_{r\to 0}^{r\to\infty}.$$ (3.16)

With this ansatz, the action (3.1) is written in terms of $S=-\int\differential tM_{\mathrm{sol}}$, where $M_{\mathrm{sol}}$ is the classical mass of the soliton

$$\begin{split}M_{\mathrm{sol}}&=4\pi\kappa\int\differential r\differential z\quantity(\frac{r^{2}}{2}kF_{rz}^{2}+h|D_{r}\Phi|^{2}+k|D_{z}\Phi|^{2}+\frac{h}{2r^{2}}(1-\absolutevalue{\Phi}^{2})^{2})\\ &-4\pi\kappa\int\differential r\differential z\frac{r^{2}}{2}\quantity(h(\partial_{r}a_{0})^{2}+k(\partial_{z}a_{0})^{2})\\ &+4\pi\kappa\gamma\int\differential r\differential za_{0}(F_{rz}(1-\absolutevalue{\Phi}^{2})+2\imaginary(D_{r}\Phi^{\dagger}D_{z}\Phi)),\end{split}$$ (3.17)

where $\gamma=\frac{N_{c}}{16\pi^{2}\kappa}$. The Euler-Lagrange equations of (3.17) are

$$(hD_{r}^{2}+D_{z}(kD_{z}))\Phi+\frac{h}{r^{2}}(1-\absolutevalue{\Phi}^{2})\Phi =i\gamma\quantity(\partial_{r}a_{0}D_{z}\Phi-\partial_{z}a_{0}D_{r}\Phi),$$ (3.18)

$$\partial_{z}(r^{2}kF_{zr})+2h\imaginary(\Phi^{\dagger}D_{r}\Phi) =\gamma\partial_{z}a_{0}(1-\absolutevalue{\Phi}^{2}),$$ (3.19)

$$\partial_{r}(r^{2}F_{rz})+2k\imaginary(\Phi^{\dagger}D_{z}\Phi) =-\gamma\partial_{r}a_{0}(1-\absolutevalue{\Phi}^{2}),$$ (3.20)

$$(h\partial_{r}(r^{2}\partial_{r})+r^{2}\partial_{z}(k\partial_{z}))a_{0} =-\gamma\quantity(F_{rz}(1-\absolutevalue{\Phi}^{2})+2\imaginary(D_{r}\Phi^{\dagger}D_{z}\Phi)).$$ (3.21)

These equations can also be derived from substituting the ansatz into the original EOM (3.5).

3.3 The energy-momentum tensor in Witten’s ansatz

In AdS/CFT correspondence, the proper method to calculate the expectation value of the EMT has been established [Balasubramanian:1999re, deHaro:2000vlm]. The application of this method to a baryon in holographic QCD was investigated in [Fujita:2022jus], and it was shown that the value of the D-term can be calculated from the classical 1+3 dimensional EMT defined as

$$\begin{split}T^{\mu\nu}&:=\evaluated{\frac{2}{\sqrt{-g}}\functionalderivative{S}{g_{\mu\nu}}}_{g=\eta}=2\kappa\int\differential z\tr(k\mathcal{F}^{z\mu}\mathcal{F}^{z\nu}+h\mathcal{F}^{\mu\rho}{\mathcal{F}^{\nu}}_{\rho}-\frac{1}{4}\eta^{\mu\nu}\quantity(h\mathcal{F}_{\rho\sigma}^{2}+2k\mathcal{F}_{z\rho}^{2})).\end{split}$$ (3.22)

Substituting the ansatz (3.7) to (3.10), we have the spherically symmetric EMT

$$T^{00}({\bf\it x}) =\epsilon(r),$$ (3.23)

$$T^{0i}({\bf\it x}) =0,$$ (3.24)

$$T^{ij}({\bf\it x}) =p(r)\delta_{ij}+s(r)\quantity(n_{i}n_{j}-\frac{1}{3}\delta_{ij}),$$ (3.25)

where the energy density $\epsilon(r)$, the pressure $p(r)$, and the shear force $s(r)$ are

$$\begin{split}\epsilon(r)&=\epsilon^{SU(2)}(r)+\epsilon^{U(1)}(r)=\int\differential z\epsilon_{z}^{SU(2)}(r,z)+\int\differential z\epsilon_{z}^{U(1)}(r,z),\\ \epsilon_{z}^{SU(2)}(r,z)&=\kappa\quantity(\frac{1}{2}kF_{rz}^{2}+\frac{1}{r^{2}}(h|D_{r}\Phi|^{2}+k|D_{z}\Phi|^{2})+\frac{h}{2r^{4}}(1-\absolutevalue{\Phi}^{2})^{2}),\\ \epsilon_{z}^{U(1)}(r,z)&=\kappa\frac{1}{2}(h(\partial_{r}a_{0})^{2})+k(\partial_{z}a_{0})^{2}),\end{split}$$ (3.26)

$$\begin{split}p(r)&=p^{SU(2)}(r)+p^{U(1)}(r),\\ p^{SU(2)}(r)&=\frac{\kappa}{3}\int\differential z\quantity(-\frac{1}{2}kF_{rz}^{2}+\frac{1}{r^{2}}(h|D_{r}\Phi|^{2}-k|D_{z}\Phi|^{2})+\frac{h}{2r^{4}}(1-\absolutevalue{\Phi}^{2})^{2}),\\ p^{U(1)}(r)&=\frac{\kappa}{3}\int\differential z\quantity(\frac{1}{2}h(\partial_{r}a_{0})^{2}+\frac{3}{2}k(\partial_{z}a_{0})^{2}),\end{split}$$ (3.27)

$$\begin{split}s(r)&=s^{SU(2)}(r)+s^{U(1)}(r),\\ s^{SU(2)}(r)&=\kappa\int\differential z\quantity(kF_{rz}^{2}+\frac{1}{r^{2}}(h|D_{r}\Phi|^{2}-k|D_{z}\Phi|^{2})-\frac{h}{r^{4}}(1-\absolutevalue{\Phi}^{2})^{2}),\\ s^{U(1)}(r)&=-\kappa\int\differential zh(\partial_{r}a_{0})^{2}.\end{split}$$ (3.28)

The classical mass of the soliton is calculated as

$$M_{\mathrm{EMT}}=M^{SU(2)}+M^{U(1)}=4\pi\int\differential rr^{2}\quantity(\epsilon^{SU(2)}(r)+\epsilon^{U(1)}(r)).$$ (3.29)

It is easy to show that $M_{\mathrm{sol}}=M_{\mathrm{EMT}}$ using the $U(1)$ part Gauss-law equation (3.21).

As introduced in section 2.1, we can derive the D-term in two ways, from the pressure $p$ (2.17) and the shear force $s$ (2.14)

$$D_{p} =D_{p}^{SU(2)}+D_{p}^{U(1)}=4\pi M_{\mathrm{sol}}\int\differential rr^{4}(p^{SU(2)}(r)+p^{U(1)}(r))$$ (3.30)

$$D_{s} =D_{s}^{SU(2)}+D_{s}^{U(1)}=-\frac{16\pi}{15}M_{\mathrm{sol}}\int\differential rr^{4}(s^{SU(2)}(r)+s^{U(1)}(r)).$$ (3.31)

$D_{p}$ and $D_{s}$ must be equal due to the conservation law (2.15).

3.4 The baryon mean square radii

There are several ways to estimate the size of the soliton. One way is the mean square radius of the energy density

$$\expectationvalue{r^{2}}_{\epsilon}:=\frac{1}{\int\differential^{3}{\bf\it x}\epsilon(r)}\int\differential^{3}{\bf\it x}r^{2}\epsilon(r)=\frac{1}{M_{\mathrm{sol}}}\int\differential^{3}{\bf\it x}r^{2}\epsilon(r).$$ (3.32)

In the spherically symmetric case, it has been argued that the radial force $F_{r}$

$$F_{r}:=T_{ij}n^{i}n^{j}=p(r)+\frac{2}{3}s(r)$$ (3.33)

is always positive in a stable system [Perevalova:2016dln] (see also [Polyakov:2018zvc]). The mean square radius of $F_{r}$

$$\expectationvalue{r^{2}}_{\mathrm{mech}}:=\frac{\int\differential^{3}{\bf\it x}r^{2}F_{r}}{\int\differential^{3}{\bf\it x}F_{r}}.$$ (3.34)

is called the mechanical radius.

4 Analysis and Results

In this section, we show our results for a soliton solution with $B=1$ obtained by solving the EOMs (3.18)-(3.21) numerically. The parameters $M_{\mathrm{KK}}$ and $\lambda$ are taken to be $949\,$\mathrm{M}\mathrm{e}\mathrm{V}$$ and $16.65$, respectively, which are the values used in [Sakai:2004cn, Sakai:2005yt] to reproduce the pion decay constant $f_{\pi}$ and the mass of the $\rho$ meson. We use the coordinate $w:=\arctan{z}$ instead of $z$ for the calculation. In this coordinate, we have $k(w)=1$, $h(w)=(\cos{w})^{-\frac{4}{3}}$, and the boundary corresponding to $z\to\pm\infty$ is $w=\pm w_{\mathrm{max}}=\pm\frac{\pi}{2}$. The cutoff of $r$ is set to $R_{\mathrm{max}}=5\pi$ in the $M_{\mathrm{KK}}=1$ unit, and this is physically $3.3\,$\mathrm{f}\mathrm{m}$$. The $r$ - $w$ plane is first reduced to $w\geq 0$ by adopting the $\mathbb{Z}_{2}$ symmetry associated with $w\to-w$, and discretized to $30\times 300$ square lattice with spacing $\Delta=\pi/60$ or $0.01\,$\mathrm{f}\mathrm{m}$$. EOMs (3.18)-(3.21) are solved by the Gauss-Seidel method. For more details, see the Appendix A.

Figure 1 are the plots of $4\pi r^{2}\epsilon_{w}(r,w):=4\pi r^{2}{\derivative{z}{w}}(\epsilon_{z}^{SU(2)}+\epsilon_{z}^{U(1)})$ and $\mathrm{ch}_{2}(r,w)$ (3.15) on the $r$ - $w$ plane. The peak of $\epsilon_{w}$ is located at $w\neq 0$. This is because $\epsilon_{w}$ includes the factor $h(w)$, which is an increasing function for $w>0$. In contrast, the center of $\mathrm{ch}_{2}$, which does not depend on $h(w)$, stays on the $w=0$ axis.

Integrating Figure 1 along the $w$ axis, we obtain Figure 2. In the outer region $r\gtrsim 1.0\,$\mathrm{f}\mathrm{m}$$, $r^{2}\epsilon(r)$ behaves as $r^{-4}$, while $\mathrm{ch}_{2}(r):=\int\differential w\mathrm{ch}_{2}(r,w)$ decreases as $r^{-7}$. This is consistent with the asymptotic behavior in large $r$ derived analytically in [Cherman:2011ve] and the EMT of the Skyrme model [Cebulla:2007ei]. The pressure $p(r)$ and the shear force $s(r)$ defined in (2.12) are shown in Figure 4. $p(r)$ is positive in $r\leq 0.64\,$\mathrm{f}\mathrm{m}$$ and negative outside. This behavior is consistent with experimental results [Burkert:2018bqq, Wang:2022ndz, Goharipour:2025lep] and calculations in other models [Neubelt:2019sou, Shanahan:2018nnv, Anikin:2019kwi, Dehghan:2025ncw, Cebulla:2007ei, Jung:2013bya, GarciaMartin-Caro:2023toa, Goeke:2007fp, Kim:2012ts, Mamo:2019mka, Mamo:2021krl, Mamo:2022eui]. To check the conservation law (2.15), we evaluate

$$\int_{0}^{\infty}\differential r\absolutevalue{p^{\prime}+2s^{\prime}+\frac{2}{r}s}\simeq 1.4\,$\mathrm{M}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{f}\mathrm{m}^{3}$.$$ (4.1)

This is about $1\,\%$ of $p(0)\simeq 114\,$\mathrm{M}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{f}\mathrm{m}^{3}$$, and it is small enough within the numerical precision.

Table 1 is the list of global properties and the sizes of the baryon evaluated in our numerical calculation along with the value of D-term in [Fujita:2022jus]. The classical mass of the baryon (3.29) is $M_{SU(2)}=950\,$\mathrm{M}\mathrm{e}\mathrm{V}$$, $M_{U(1)}=232\,$\mathrm{M}\mathrm{e}\mathrm{V}$$, and $M_{\mathrm{sol}}=1.18\,$\mathrm{G}\mathrm{e}\mathrm{V}$$ in total. This value agrees well with the value obtained in [Hori:2023fxq], although this is $416\,$\mathrm{M}\mathrm{e}\mathrm{V}$$ smaller than the value evaluated in [Hata:2007mb]. The values of D-term are $D_{p}=-2.06$ and $D_{s}=-2.05$. The difference between $D_{p}$ and $D_{s}$ is less than $1\,\%$, and agrees with good precision. As mentioned in section 2.2, the individual values of the $SU(2)$ part and $U(1)$ part of $D_{p}$ and $D_{s}$ do not agree because $D_{p}^{SU(2)}$ and $D_{p}^{U(1)}$ involve the contribution of a non-conserving term $\bar{c}$. The $t$-dependence of $D(t)$ can be obtained from (2.13), and the result is shown in Figure 4. $D(t)$ is negative anywhere, and converges to $0$ at large $t$. The slope at $t=0$ is infinitely steep, since $\derivative{t}D(t)$ is proportional to $\int\differential rr^{6}s(r)$ at small $t$, which diverges in the chiral limit.

The two mean radii defined in section 3.4 are $\expectationvalue{r^{2}}_{\epsilon}=(0.662\,$\mathrm{f}\mathrm{m}$)^{2}$ and $\expectationvalue{r^{2}}_{\mathrm{mech}}=(0.938\,$\mathrm{f}\mathrm{m}$)^{2}$. These values are comparable to other recent results .

To compare our numerical results with the previous calculation in [Fujita:2022jus], we use $D_{s}$ because the D-term was calculated from the traceless part of the stress tensor in [Fujita:2022jus]. We find qualitatively similar behavior such that the $SU(2)$ part is negative and $U(1)$ is positive, being negative in total in both cases. However, the absolute value of the $SU(2)$ part from our numerical result is about 4 times larger than that obtained in [Fujita:2022jus], making the total D-term more negative.

	$M_{\mathrm{sol}}$	$D_{p}$	$D_{s}$	$D$ in [Fujita:2022jus]	$\expectationvalue{r^{2}}_{\epsilon}$	$\expectationvalue{r^{2}}_{\mathrm{mech}}$
results	$1.18\,[$\mathrm{G}\mathrm{e}\mathrm{V}$]$	$-2.06$	$-2.05$	$-0.140$	$(0.662)^{2}\,[$\mathrm{f}\mathrm{m}$^{2}]$	$(0.938)^{2}\,[$\mathrm{f}\mathrm{m}$^{2}]$
$SU(2)$ part	$950\,[$\mathrm{M}\mathrm{e}\mathrm{V}$]$	$-3.42$	$-2.54$	$-0.685$	$(0.692)^{2}\,[$\mathrm{f}\mathrm{m}$^{2}]$	$(1.38)^{2}\,[$\mathrm{f}\mathrm{m}$^{2}]$
$U(1)$ part	$232\,[$\mathrm{M}\mathrm{e}\mathrm{V}$]$	$1.36$	$0.489$	$0.543$	$(0.529)^{2}\,[$\mathrm{f}\mathrm{m}$^{2}]$	$(0.267)^{2}\,[$\mathrm{f}\mathrm{m}$^{2}]$

Table 1: Mass, D-term and mean square radii.

Figure 1: Energy density and Chern number density of the baryon on $r$ - $w$ plane. The plot range is restricted to $0\leq w\leq w_{\mathrm{max}}=\pi/2$, which corresponds to $0\leq z<+\infty$.

Figure 2: Energy density and Chern number density along the $r$ direction. The black dashed line is $\mathrm{ch}_{2}(r,w)$ integrated along the $w$ axis. The blue solid line is $4\pi r^{2}\epsilon(r)$ normalized with $M_{\mathrm{sol}}$. The right figure is plotted in the log-scale. The red solid and dashed lines are $1/r^{7}$ and $1/r^{4}$, respectively.

Figure 3: $p(r)$ and $s(r)$ multiplied with $4\pi r^{2}$. The blue solid line is $p(r)$ and the red dashed line is $s(r)$.

Figure 4: $t$ dependence of $D(t)$. The graph is derived from (2.13) for $t\neq 0$, and $D(0)$ is fixed to $D_{s}=-2.05$.

5 Conclusions

In this paper, we investigated the gravitational form factors of a baryon using a top-down holographic description of QCD proposed in [Sakai:2004cn, Sakai:2005yt, Hata:2007mb]. In this framework, a baryon is described as a topologically stable solitonic gauge configuration in a 5-dimensional gauge theory. In order to obtain a classical solution corresponding to a baryon, we used Witten’s ansatz and solved the equations of motion numerically. Using this numerical solution, the energy-momentum tensor was evaluated and the quantities such as energy density, pressure, and shear force were extracted. The results are summarized in Table 1 and Figures 1–4. In particular, the D-term turned out to be significantly larger in absolute value compared to the previous result in [Fujita:2022jus].

Our analysis is not complete and there are several directions to improve the calculation. First of all, our analysis is based on a static classical solution of the soliton corresponding to the baryon. It would be interesting to combine our numerical analysis with the quantization of the fluctuations around the soliton solution investigated in [Hata:2007mb, Panico:2008it, Hashimoto:2008zw, Cherman:2011ve]. An important ingredient is the spin of the baryon, which can be introduced by quantizing the collective coordinates related to the rotational degrees of freedom. Once it is included, it will be possible to analyze the $J(t)$ form factor for nucleon (see (2.1) and (2.3)) and the gravitational form factors for higher spin baryons (e.g. $\Delta$ resonance). (See [Kim:2020lrs] for a work in this direction.) The analysis of the electromagnetic form factors done in [Hashimoto:2008zw] can also be improved by using our numerical analysis. Secondly, the model used in this paper describes QCD with massless quarks and the effect of the current quark mass is neglected. There are some proposals to add the quark masses in this system [Aharony:2008an, Hashimoto:2008zw]. Since the effect of the quark mass drastically changes the IR behavior of some quantities, such as $\frac{d}{dt}D(t)|_{t\rightarrow 0}$ mentioned in section 4, it would be worthwhile to include it in the analysis of the form factors. For the study of hyperons, it will be crucial to include the strange quark mass to obtain realistic predictions. (See [Hashimoto:2009st]) Thirdly, our analysis is based on the expression of the EMT in (3.22). As it was shown in [Fujita:2022jus], (3.22) is obtained as an approximation that is valid when the momentum transfer is smaller than the mass scale of the lightest glueball, which is of order 1 GeV. Therefore, for $t\sim{\cal O}(1)\,\mbox{GeV}^{2}$ or higher, the glueball propagators (or, in other words, the bulk to boundary propagator) should be taken into account properly. In addition, the analysis in the holographic description is done in the supergravity approximation, in which the $\alpha^{\prime}$ corrections (corresponding to the $1/\lambda$ corrections) as well as quantum gravity effects (corresponding to the $1/N_{c}$ corrections) are neglected. It would be nice if these corrections could be incorporated.

Appendix A Boundary conditions

In this appendix, we summarize the boundary conditions we used to solve the EOM (3.18) - (3.21). Under the ansatz (3.7) - (3.10), the 5-dimensional gauge field $\mathcal{A}$ is reduced to real and complex scalars $a_{0}$ and $\Phi$, and a $U(1)$ gauge field $(A_{z},A_{r})$ on a 2-dimensional $z$-$r$ plane. Replacing $z$ with $w=\arctan(z)$, the boundaries of this plane are $r=0$, $w=\pm\frac{\pi}{2}$, and $r\to\infty$.

From the regularity of $\mathcal{A}$, the boundary condition at $r=0$ are

$$\Phi=-i,\,A_{w}=0,\,\partial_{r}A_{r}=0,\,\partial_{r}a_{0}=0,\quad(r=0).$$ (A.1)

Because of the divergent factor $h(w)=\cos(w)^{-\frac{4}{3}}$, $\mathcal{F}_{\mu\nu}$ must vanish at $w=\pm\frac{\pi}{2}$ to make energy finite. Therefore, $\mathcal{A}_{\mu}$ is pure gauge at $w=\pm\frac{\pi}{2}$. We choose a gauge such that $\mathcal{A}_{\mu}\to 0$ at this boundary in this paper. Therefore we have

$$\Phi=-i,\,A_{r}=0,\,a_{0}=0,\quad\quantity(w=\pm\frac{\pi}{2}).$$ (A.2)

Similarly, $\mathcal{F}_{\mu\nu}$ and $\mathcal{F}_{w\mu}$ should be $0$ at $r\to\infty$. However, we cannot let $\mathcal{A}\to 0$ as the boundary condition at $r\to\infty$. This is because we have set $\mathcal{A}_{\mu}\to 0$ at $w=\pm\frac{\pi}{2}$, and the baryon number (3.16) is now

$$B=\evaluated{-\frac{1}{2\pi}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}A_{w}\differential w}_{r\to\infty}=\evaluated{-\frac{1}{2\pi}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\imaginary(\Phi^{\dagger}\partial_{w}\Phi)\differential w}_{r\to\infty}.$$ (A.3)

This leads to the non-trivial behavior of $A_{w}$ and $\Phi$ at $r\to\infty$. We take a simple configuration

$$\Phi\to ie^{-2iw},\,A_{w}\to-2,\,A_{r}\to 0,\,a_{0}\to 0,\quad(r\to\infty).$$ (A.4)

In the chiral limit, it is known that the EMT of a baryon decays with a power law as $r\to\infty$. Thus, merely using (A.4) with the finite cutoff is not appropriate. Because of this, we set the asymptotic behavior of the fields following the discussion in [Cherman:2011ve, Bolognesi:2013nja] as

$$A_{w}+2=\order{1/r^{2}},\,|\Phi|^{2}-1=\order{1/r^{4}},$$ (A.5)

at $r=R_{\mathrm{cutoff}}$. For other fields, we take the value at $r\to\infty$ as the boundary condition on $r=R_{\mathrm{cutoff}}$, since they converge faster than $1/r^{5}$.

[title = References]