Gravitational form factors of light mesons

TL;DR

This paper computes the gravitational form factors of light mesons using the NJL model, revealing insights into their internal structure and mass distribution, and compares these with charge radii and experimental data.

Contribution

It introduces a method to calculate meson gravitational form factors using the canonical EMT and demonstrates the importance of dressing the graviton vertex for accurate results.

Findings

Light cone mass radii: pion 0.27 fm, sigma 0.32 fm, rho 0.25 fm.

Pion and rho mass radii are smaller than their charge radii.

The pion mass radius agrees with KEKB experimental data.

Abstract

We calculate the gravitational form factors of the pion, sigma meson, and rho meson in the Nambu-Jona-Lasinio (NJL) model of quantum chromodynamics. The canonical energy-momentum tensor (EMT) is used in their derivation, allowing the possibility of an antisymmetric contribution when the hadron has intrinsic spin. We show that the asymmetric graviton vertex arising from the canonical EMT satisfies a simpler Ward-Takahashi identity (WTI) than the symmetric graviton vertex of the Belinfante EMT. The necessity of fully dressing the graviton vertex through the relevant Bethe-Salpeter equation is demonstrated for observing both the WTI and a low-energy pion theorem. Lastly, we calculate static moments of the meson EMT decompositions, obtaining predictions for the meson mass radii. We find light cone mass radii of 0.27 fm for the pion, 0.32 fm for the sigma, and 0.25 fm for the rho. For the pion and rho, these are smaller than the light cone charge radii, respectively 0.51 fm and 0.45 fm, while we have a sigma charge radius of zero. Our light cone pion mass radius agrees with a phenomenological extraction from KEKB data.

Tables

Table 1. Table 1: NJL model parameters used in this work. All but G f subscript 𝐺 𝑓 G_{f} are adopted from Ref. Cloët  et al. ( 2014 ) . G f subscript 𝐺 𝑓 G_{f} is determined by requiring an f 1 subscript 𝑓 1 f_{1} pole at m f 1 = 1.28 subscript 𝑚 subscript 𝑓 1 1.28 m_{f_{1}}=1.28  GeV in the NJL T-matrix. Couplings are in units GeV -2 , while the dressed quark mass M 𝑀 M and regulators are in units GeV.

${\Lambda }_{\mathrm{IR}}$	${\Lambda }_{\mathrm{UV}}$	$M$	$G_{\pi }=G_{\eta }$	$G_{\rho }$	$G_{\omega }$	$G_f$	$G_T$
0.240	0.645	0.4	19.0	11.0	10.4	0.82	0

Table 2. Table 2: Mean squared mass radius of spin-zero mesons in the NJL model, both using the Breit frame and light cone prescriptions. For the pion, an empirical value for the light cone mass radius extracted from KEKB data is included for comparison. All values are in fm.

	Breit frame	Light cone	Empirical light cone Kumano et al. (2018) [see text]
Pion	1.28	0.27	$0.26\sim 0.32$
Sigma	0.56	0.32

Table 3. Table 3: Static properties of the rho meson in the NJL model. Electromagnetic properties taken from the BSE calculation (without pion cloud effects) for ρ + superscript 𝜌 \rho^{+} from Cloët  et al. ( 2014 ) . All radii are in fm, the mass quadrupole moment is in units of m ρ subscript 𝑚 𝜌 m_{\rho} -fm 2 , and the electric quadrupole moment is in e 𝑒 e -fm 2 .

	$\sqrt{{⟨r^2⟩}_{\mathrm{mass}}}$	$\sqrt{{⟨r^2⟩}_{\mathrm{elec}.}}$	${𝒬}_{\mathrm{mass}}$	${𝒬}_{\mathrm{elec}.}$
Breit frame	0.45	0.67	-0.0224	-0.0200
Light cone	0.25	0.45

Equations

$$\mathcal{L}$$

$$+\frac{1}{2}{G}_{\eta }[(\overline{\psi }\mathbf{\tau }\psi {)}^2-(\overline{\psi }{\gamma }_5\psi {)}^2]-\frac{1}{2}{G}_f(\overline{\psi }{\gamma }_{\mu }{\gamma }_5\psi {)}^2-\frac{1}{2}{G}_T[(\overline{\psi }i{\sigma }^{\mu \nu }\psi {)}^2-(\overline{\psi }i{\sigma }^{\mu \nu }\mathbf{\tau }\psi {)}^2],$$

$$M=m+2i{G}_{\pi }(2N_c)\int \frac{{\mathrm{d}}^{4}k}{(2\pi {)}^4}{\mathrm{Tr}}_D\left[S(k)\right]=m+8i{G}_{\pi }(2N_c)\int \frac{{\mathrm{d}}^{4}k}{(2\pi {)}^4}\frac{M}{k^2-M^2+i0},$$

$$i{{\mathcal{T}}_{\pi }(q)}_{\alpha \beta ,\gamma \delta }$$

$$i{{\mathcal{T}}_{\sigma }(q)}_{\alpha \beta ,\gamma \delta }$$

$$i{{\mathcal{T}}_{\rho }(q)}_{\alpha \beta ,\gamma \delta }$$

$$i{{\mathcal{T}}_{f_1}(q)}_{\alpha \beta ,\gamma \delta }$$

$$1+2G_{\pi }{\Pi }_{PP}(m_{\pi }^2)$$

$$1-2G_{\pi }{\Pi }_{SS}(m_{\pi }^2)$$

$$1+2G_{\rho }{\Pi }_{VV}(m_{\rho }^2)$$

$$1+2G_f{\Pi }_{AA}^{(T)}(m_{f_1}^2)$$

$${\Gamma }_{\pi }^i$$

$${\Gamma }_{\sigma }^i$$

$${\Gamma }_{\rho }^i$$

$$\displaystyle Z_{\pi}^{-1}=-\frac{\partial}{\partial q^{2}}\Pi_{PP}(q^{2})\bigg{|}_{q^{2}=m_{\pi}^{2}}$$

$$\displaystyle Z_{\sigma}^{-1}=+\frac{\partial}{\partial q^{2}}\Pi_{SS}(q^{2})\bigg{|}_{q^{2}=m_{\sigma}^{2}}$$

$$\displaystyle Z_{\rho}^{-1}=-\frac{\partial}{\partial q^{2}}\Pi_{VV}(q^{2})\bigg{|}_{q^{2}=m_{\rho}^{2}}\,.$$

$$\frac{1}{X^n}=\frac{1}{(n-1)!}{\int }_0^{\infty }\mathrm{d}\tau {\tau }^{n-1}{e}^{-\tau X}\to \frac{1}{(n-1)!}{\int }_{1/{\Lambda }_{\mathrm{UV}}^2}^{1/{\Lambda }_{\mathrm{IR}}^2}\mathrm{d}\tau {\tau }^{n-1}{e}^{-\tau X}.$$

$${\Delta }_{\mu }{\Gamma }_G^{\mu \nu }(p^{\prime },p)=p^{\nu }{S}^{-1}(p^{\prime })-p^{\prime \nu }{S}^{-1}(p).$$

$$M^{\mu \alpha \beta }=x^{\alpha }{T}^{\mu \beta }-x^{\beta }{T}^{\mu \alpha }+S^{\mu \alpha \beta }$$

$$\displaystyle\vbox{\hbox{\includegraphics[scale={0.5}]{graviton_qq.pdf}}}=\gamma_{G}^{\mu\nu}\left(k+\frac{\Delta}{2},k-\frac{\Delta}{2}\right)=\gamma^{\mu}k^{\nu}-g^{\mu\nu}(\not{k}-m)\,.$$

$${\Delta }_{\mu }{\Gamma }_G^{\{\mu \nu \}}(p^{\prime },p)=p^{\nu }{S}^{-1}(p^{\prime })-p^{\prime \nu }{S}^{-1}(p)+\frac{1}{2i}{\Delta }_{\mu }\left[S^{-1}(p^{\prime }){\Sigma }^{\mu \nu }-{\Sigma }^{\mu \nu }{S}^{-1}(p)\right],$$

$$\displaystyle T^{\mu\nu}(x)=\bar{\psi}(x)(i\gamma^{\mu}\overleftrightarrow{\partial}^{\nu})\psi(x)-g^{\mu\nu}\bar{\psi}(x)(i\overleftrightarrow{\not{\partial}}-m)\psi(x)-\frac{1}{2}g^{\mu\nu}\sum_{\Omega}G_{\Omega}\big{(}\bar{\psi}(x)\Omega\psi(x)\big{)}\big{(}\bar{\psi}(x)\Omega\psi(x)\big{)}\,,$$

$$\displaystyle\vbox{\hbox{\includegraphics[scale={0.5}]{graviton_qqqq.pdf}}}=\gamma_{Gqqqq}^{\mu\nu}=-g^{\mu\nu}\sum_{\Omega}2G_{\Omega}\Omega\otimes\Omega\,.$$

$$\displaystyle\vbox{\hbox{\includegraphics[scale={0.5}]{graviton_qq_condensate.pdf}}}=2iG_{\pi}(2N_{c})g^{\mu\nu}\int\frac{\mathrm{d}^{4}k}{(2\pi)^{4}}\mathrm{Tr}_{D}\left[S(k)\right]=g^{\mu\nu}(M-m)\,,$$

$$\displaystyle\vbox{\hbox{\includegraphics[scale={0.5}]{graviton_qq_dressed.pdf}}}=\Gamma_{Gqq}^{\mu\nu}\left(k+\frac{\Delta}{2},k-\frac{\Delta}{2}\right)=\gamma^{\mu}k^{\nu}-g^{\mu\nu}(\not{k}-M)+\frac{\Delta^{\mu}\Delta^{\nu}-\Delta^{2}g^{\mu\nu}}{4M}C_{Q}(t)+\frac{i\epsilon^{\mu\nu\Delta\sigma}\gamma_{\sigma}\gamma_{5}}{4}D^{\prime}_{Q}(t)\,.$$

$$C_Q(t)$$

$$D_Q^{\prime }(t)$$

$$⟨p^{\prime }\mid T_{\mu \nu }^a(0)\mid p⟩$$

$$\sum _{a=q,g}{\bar{c}}_a(t)=0$$

$$\sum _{a=q,g}{A}_a(0)=1.$$

$$\underset{m_{\pi }\to 0}{\lim }\sum _{a=q,g}{C}_a(0)=-1.$$

$$⟨p^{\prime },{\lambda }^{\prime }\mid T_{\mu \nu }^a(0)\mid p,\lambda ⟩$$

$$-\frac{1}{2}({\Delta }_{\mu }{\Delta }_{\nu }-{\Delta }^2{g}_{\mu \nu })\left[({ϵ}^{\prime \ast }ϵ){\mathcal{G}}_3^a(t)-\frac{(\Delta {ϵ}^{\prime \ast })(\Delta ϵ)}{2m_{\rho }^2}{\mathcal{G}}_4^a(t)\right]+P_{\{\mu }\left({ϵ}_{\nu \}}^{\prime \ast }(\Delta ϵ)-{ϵ}_{\nu \}}(\Delta {ϵ}^{\prime \ast })\right){\mathcal{G}}_5^a(t)$$

$$+\frac{1}{2}\left[{\Delta }_{\{\mu }\left({ϵ}_{\nu \}}^{\prime \ast }(\Delta ϵ)+{ϵ}_{\nu \}}(\Delta {ϵ}^{\prime \ast })\right)-{ϵ}_{\{\mu }^{\prime \ast }{ϵ}_{\nu \}}{\Delta }^2-g_{\mu \nu }(\Delta {ϵ}^{\prime \ast })(\Delta ϵ)\right]{\mathcal{G}}_6^a(t)$$

$$+{ϵ}_{\{\mu }^{\prime \ast }{ϵ}_{\nu \}}{m}_{\rho }^2{\mathcal{G}}_7^a(t)+g_{\mu \nu }{m}_{\rho }^2({ϵ}^{\prime \ast }ϵ){\mathcal{G}}_8^a(t)+\frac{1}{2}{g}_{\mu \nu }(\Delta {ϵ}^{\prime \ast })(\Delta ϵ){\mathcal{G}}_9^a(t)$$

$$+P_{[\mu }\left({ϵ}_{\nu ]}^{\prime \ast }(\Delta ϵ)-{ϵ}_{\nu ]}(\Delta {ϵ}^{\prime \ast })\right){\mathcal{G}}_{10}^a(t)+{\Delta }_{[\mu }\left({ϵ}_{\nu ]}^{\prime \ast }(\Delta ϵ)+{ϵ}_{\nu ]}(\Delta {ϵ}^{\prime \ast })\right){\mathcal{G}}_{11}^a(t),$$

$$\sum _{a=q,g}{\mathcal{G}}_7^a(t)=\sum _{a=q,g}{\mathcal{G}}_8^a(t)=\sum _{a=q,g}{\mathcal{G}}_9^a(t)=\sum _{a=q,g}{\mathcal{G}}_{11}^a(t)=0,$$

$$\sum _{a=q,g}{\mathcal{G}}_1^a(0)=1$$

Full text

I Introduction

In recent years, the energy-momentum tensor (EMT) of hadrons has become an increasingly popular object of study in our quest to better understand hadronic structure—and, in turn, quantum chromodynamics (QCD) itself. (See Refs. Leader and Lorcé (2014); Polyakov and Schweitzer (2018) and references therein.) Understanding the EMT can help address such fundamental questions as where does the mass of the proton come from, and where is all of its spin? It also opens new avenues for exploration, including not just spatial distributions of energy, momentum, and angular momentum Polyakov (2003); Leader and Lorcé (2014); Lorcé (2018a); Lorcé et al. (2018), but also the distribution of pressure and shear forces Polyakov and Shuvaev (2002); Polyakov (2003); Polyakov and Schweitzer (2018); Lorcé et al. (2019).

It has been remarked that the EMT introduces a new intrinsic global quantity in addition to mass and spin, called the “D-term” Polyakov and Shuvaev (2002); Polyakov and Schweitzer (2018). The D-term, which quantifies the strength of the forces binding the hadron together, is not constrained by conservation laws nor by representation theory of the Lorentz group. However, in certain cases it may be constrained by other considerations—as it is with Nambu-Goldstone bosons Novikov and Shifman (1981); Voloshin and Zakharov (1980); Polyakov and Weiss (1999). Since the pion plays a vital role in QCD as the Nambu-Goldstone boson of dynamical chiral symmetry breaking (DCSB) Nambu and Jona-Lasinio (1961a, b), understanding how its D-term comes to be constrained may play as important a role in fully grasping QCD as understanding the origin of proton mass and spin.

DCSB has come to be understood as one of the central features of QCD, and is intimately involved in the observed mass of hadrons Bashir et al. (2012). Accordingly, understanding the pion is vital to understanding QCD, and pion structure has long been a major topic of study in hadron physics. Its EMT in particular has recently been computed on the lattice Brommel (2007); Shanahan and Detmold (2019). Additionally, it has become possible through dispersive analysis of KEKB data for the reaction $\gamma^{*}\gamma\rightarrow\pi^{0}\pi^{0}$ to extract an empirical parametrization of the pion EMT in the timelike region Kumano et al. (2018). The prospect of comparison to empirical values makes model calculations of the pion EMT especially important now.

In this work, we investigate the EMT of the pion and of other mesons in the Nambu–Jona-Lasinio (NJL) model Vogl and Weise (1991); Klevansky (1992); Hatsuda and Kunihiro (1994) of QCD. The NJL model is a Poincaré covariant quantum field theory that successfully reproduces low-energy properties of QCD such as DCSB. Additionally, since the model does not contain gluons, we can defer issues regarding gauge invariance of the EMT. These properties make the NJL model an ideal framework in which to investigate the EMT of the pion, especially aspects of the EMT that arise from DCSB. The sigma and rho mesons are also investigated as a point of contrast, with the latter serving to illustrate the ways in which spin manifests in the EMT.

This work is organized into the following sections. We first give an overview of the NJL model in Sec. II. We then develop the formalism needed to calculate the meson EMTs in Sec. III, where the formalism is cast as a study of quark-graviton and meson-graviton interactions. In Sec. IV, we develop the formalism needed to calculate the axial form factors of the rho meson, which are needed to numerically demonstrate a correspondence between certain gravitational and axial form factors. Results for the meson EMTs are given in Sec. V, where we also explore implications of our results. Finally, we conclude in Sec. VI and give a brief outlook, including the prospects of performing similar calculations in QCD.

II The NJL model

The Nambu–Jona-Lasinio (NJL) model, originally proposed as a theory of elementary nucleons Nambu and Jona-Lasinio (1961a, b), is used as a low-energy effective field theory that models the dynamical chiral symmetry breaking of quantum chromodynamics through a four-fermi contact interaction Vogl and Weise (1991); Klevansky (1992); Hatsuda and Kunihiro (1994). It has successfully been applied to modeling the physical properties of both mesons Vogl and Weise (1991); Klevansky (1992); Cloët et al. (2014); Ninomiya et al. (2017) and baryons Ishii et al. (1993a, b, 1995); Cloët et al. (2014). Accordingly, we employ the NJL model in this work to calculate the energy-momentum tensor of light mesons.

The two-flavor NJL model Lagrangian we use is Cloët et al. (2014):

[TABLE]

where $\hat{m}=\mathrm{diag}[m_{u},m_{d}]$ is the current quark mass matrix (where we take $m_{u}=m_{d}\equiv m$ in this work), $\tau_{i}$ are the isospin matrices, and the $G_{i}$ are four-fermi coupling constants. The expression in Eq. (1) is invariant under $\mathrm{SU}(2)_{V}$ and $\mathrm{SU}(2)_{A}$ transformations, but $\mathrm{U}(1)_{A}$ symmetry requires the additional constraints $G_{\eta}=G_{\pi}$ and $G_{T}=0$, which we will assume in this work. We take the Lagrangian in Eq. (1) to symbolically represent a Lagrangian that has already been Fierz symmetrized, so that only direct terms need to be calculated in the interaction kernel (See Ref. Ishii et al. (1995) for a detailed description of the Fierz symmetrization procedure.)

The NJL model dynamically generates a dressed quark mass $M$ when $G_{\pi}>G_{\mathrm{critical}}$. This dynamical mass generation is described by the gap equation:

[TABLE]

where the trace is over the Clifford matrix structure (with $2N_{c}$ having come from the color and isospin traces already).

Mesons appear in the NJL model as bound state poles in the quark-antiquark T-matrix, which itself can be found from solving an inhomogeneous Bethe-Salpeter equation, depicted diagrammatically in Fig. 1. Poles corresponding to mesons with the quantum numbers of various mesons are present in the NJL T-matrix. Our interest in this work lies primarily with the pion, sigma, and rho, but we also consider the $f_{1}$ meson as a means of determining the coupling constant $G_{f}$. The T-matrices for these four mesons can be written111 The T-matrix for the pion can have additional structure generated by $\pi$-$a_{1}$ mixing, but we neglect this as a higher-order effect in this work. :

[TABLE]

where the bubble diagrams are defined in App. A. Since the mesons appear as poles in these T-matrices, their masses are given by the conditions:

[TABLE]

The residues at the poles can be interpreted as effective quark-meson coupling constants, which allow us to find the properly normalized homogeneous Bethe-Salpeter vertex functions:

[TABLE]

where

[TABLE]

The NJL model is non-renormalizable, owing to the four-fermi contact interaction. To fully define the model, it is necessary to introduce a regularization scheme. We follow Refs. Ebert et al. (1996); Hellstern et al. (1997); Cloët et al. (2014) in using proper time regularization, with both an infrared and an ultraviolet regulator. The regularization proceeds formally through the substitution:

[TABLE]

Only the UV regulator is necessary to make the model finite, but the presence of an IR regulator ensures the two-body T-matrix is always real, and prevents the decay of mesons into two quarks. As is customary with non-remormalizable theories, these regulators are kept finite as additional model parameters.

We adopt the NJL model parameters used in previous work Cloët et al. (2014), along with a value for $G_{f}$ which produces an $f_{1}$ pole in the T-matrix with mass 1.28 GeV. The model parameters we use are given in Tab. 1.

III Gravitation and the EMT

The energy-momentum tensor (EMT) has long had strong ties to gravitation. The equivalence between inertial and gravitational mass has long suggested that energy is the charge upon which gravitation acts, and this equivalence principle is currently canonized as one of the central premises of general relativity (and its extensions, such as Einstein-Cartan theory Cartan (1923, 1924)). Accordingly, it is helpful in theoretical investigations of the EMT—even if gravitation is not the intended topic of study—to consider the physics of graviton coupling. If the EMT is the current upon which gravitons act, then one can constrain the EMT of a system of interest using traditional field theoretic methods that have been applied to other currents such as the electromagnetic current. Just as with coupling to a photon, one can obtain a Ward-Takahashi identity for graviton coupling and solve Dyson-Schwinger equations for the interaction between a graviton and a fully dressed field in a strongly-interacting, non-perturbative regime.

III.1 Gravitational Ward-Takahashi identity

The gravitational Ward-Takahashi identity (GWTI) is given by Brout and Englert (1966):

[TABLE]

This relation has applicability to fields of general spin if the canonical EMT is used as the source of the gravitational field. A proof of the general validity of Eq. (18) can be found in App. B.

The canonical EMT for which Eq. (18) has general applicability naturally arises as the Noether current associated with spacetime translation symmetry. This current is not symmetric in its indices for particles with intrinsic spin Leader and Lorcé (2014); Lorcé (2018b); Florkowski and Ryblewski (2018), as a consequence of the generalized angular momentum

[TABLE]

being a conserved quantity—in particular, the Noether current associated with Lorentz transformations. Here $S^{\mu\alpha\beta}$ is the intrinsic generalized angular momentum, and as a consequence of the conservation laws $\partial_{\mu}M^{\mu\alpha\beta}=0$ and $\partial_{\mu}T^{\mu\nu}=0$, one has $T^{\alpha\beta}-T^{\beta\alpha}=-\partial_{\mu}S^{\mu\alpha\beta}$. For particles with spin, this is not expected to be zero, so the EMT is not expected to be symmetric.

It is possible to obtain a symmetric, “Belinfante-improved” EMT by adding the divergence of a superpotential to the EMT derived through Noether’s theorem Leader and Lorcé (2014). The resulting tensor is also a conserved current, but there is no guarantee that it encodes exactly the same physical properties as the Noether current associated with spacetime translations. It has been observed that adding a total derivative to the EMT can in fact change the values of measurable quantities such as the D-term Hudson and Schweitzer (2017). We thus choose not to add any total derivatives to the EMT in this work.

One reason for preferring the Belinfante-improved EMT over the canonical EMT is that general relativity assumes the EMT to be symmetric, and spacetime to accordingly be torsion-free. However, we do not a priori know whether spacetime is really torsion-free, and whether gravitons can couple to the antisymmetric component of the EMT222 Einstein-Cartan theory Cartan (1923) is a straightforward extension of general relativity that incorporates spacetime torsion and uses an asymmetric EMT as the source of gravitation. . Since in this work we consider graviton coupling only as a theoretical means of developing the formalism for calculating the EMT—which is empirically accessed through other means, such as DVCS Ji (1997) and hadron pair production in diphoton collisions Kumano et al. (2018)—we find it most prudent to consider the graviton as capable of containing torsion and proceed to use the canonical, asymmetric EMT.

As remarked above, one major consequence of using the canonical rather than the Belinfante EMT is that Eq. (18) holds in general. One can readily confirm this for simple examples of particles with spin, such as an elementary free fermion, for which the gravitational vertex is:

[TABLE]

On the other hand, the symmetrization of the vertex $\gamma_{G}^{\{\mu\nu\}}$—which arises from using the Belinfante EMT as the gravitating current in the spin-half case—does not obey Eq. (18). In fact, several authors DeWitt (1967); Bessler et al. (1969) have found that the symmetrized graviton vertex obeys a different Ward-Takahashi identity with a more complicated structure:

[TABLE]

where $\Sigma^{\mu\nu}$ is the generator of Lorentz transforms. We consider the fact that it obeys a simpler WTI one of the virtues of the canonical EMT.

III.2 Graviton-quark vertices in NJL model

Application of Noether’s theorem to the full NJL model Lagrangian of Eq. (1) gives the following EMT:

[TABLE]

which contains the non-linear contact interaction terms. The presence of a four-fermi interaction term in the EMT entails the existence of a five-point vertex function (4 quarks and 1 graviton line) in addition to the usual three-point vertex (2 quarks and 1 graviton).

At the level of truncation we are considering here, it is not necessary to dress the five-point vertex. We do not dress the four-fermi contact interaction, after all, and it would be inconsistent to dress the graviton’s interaction with the contact interaction while not dressing the interaction itself.

The effective Feynman rule for the five-point vertex can be read off from the EMT in Eq. (22) directly:

[TABLE]

We shall see presently that the existence of this five-point vertex has non-trivial consequences, including both a vacuum condensate contribution to the Bethe-Salpeter equation for the three point vertex (see second-from-the-right diagram of Fig. 2), and the existence of a “bicycle diagram” in the calculation of bound state EMTs (see rightmost diagram of Fig. 3). Crucially, these diagrams are both necessary for energy-momentum conservation to be observed.

For the fully-dressed three-point vertex, we must solve an inhomogeneous Bethe-Salpeter equation, depicted in Fig. 2. In contrast to the photon vertex BSE, there are two driving terms. The second diagram from the right arises because it is possible for the five-point vertex to contribute to the dressing of the three-point vertex. Most interestingly, this new driving term is directly proportional to the vacuum condensate, which—as a consequence of the gap equation, Eq. (2)—is proportional to the mass dressing $(M-m)$. In particular, we have:

[TABLE]

meaning that the two driving terms together contribute $\gamma^{\mu}k^{\nu}-g^{\mu\nu}(\not{k}-M)$ to the dressed graviton vertex.

When we take into consideration that the NJL model propagator takes the same functional form as the bare propagator, but with the current quark mass replaced by the dressed mass, we observe the remarkable property that the driving terms in the graviton vertex BSE already satisfy the gravitational WTI on their own. This means the rightmost diagram in Fig. 2 must contribute a quantity transverse to the momentum transfer $\Delta$. Additionally, since the NJL interaction kernel contains no explicit momentum dependence and the $k$ dependence in the dressed kernel is integrated over in the loop, the rightmost diagram can only depend on $\Delta$. These considerations, and the fact that the EMT is $\mathsf{P}$ and $\mathsf{T}$ even, gives us the following most general form for the dressed three-point vertex:

[TABLE]

The functions $C_{Q}(t)$ and $D^{\prime}_{Q}(t)$ can be determined by algebraically solving the BSE. We find:

[TABLE]

with the bubbles given in App. A. Remarkably, $C_{Q}(t)$ has a pole in the timeline region when $t=m_{\sigma}^{2}$, and $D_{Q}^{\prime}(t)$ has a timelike pole when $t=m_{f_{1}}^{2}$, by comparison to the pole conditions in Eqs. (8,10).

III.3 EMT decomposition and gravitational form factors

Now that we have constructed the necessary formalism, we proceed to consider matrix elements of the NJL model EMT between momentum and spin eigenstates of a particular hadron, viz., $\langle p^{\prime},\lambda^{\prime}|T_{\mu\nu}(0)|p,\lambda\rangle$. For conciseness, we refer to this matrix element as the EMT of the hadron.

To begin, we must consider three diagrams when calculating the EMT of a meson in the NJL model. These diagrams are given in Fig. 3. The first two (triangle) diagrams are completely analogous to the diagrams contributing to the electromagnetic or axial vector current. The third (bicycle) diagram is new to graviton interactions, and is a consequence of the equivalence principle. The four-fermi contact interaction is a sort of potential energy in the NJL model, and gravitation couples to all energy in exactly the same way.

The EMT of any hadron can be decomposed into a finite number of independent algebraic structures, each multiplied by a Lorentz-invariant function only of the invariant momentum transfer $t=(p^{\prime}-p)^{2}$. These functions are called gravitational form factors (GFFs), and are analogous to the electromagnetic form factors appearing in decompositions of the electromagnetic current $\langle p^{\prime},\lambda^{\prime}|j_{\mu}(0)|p,\lambda\rangle$. Like with electromagnetic form factors, the number of GFFs depends on the spin of the hadron, and this number increases with increasing spin.

The most general form the spin-zero EMT can take for a particular parton flavor $a$ is Polyakov and Schweitzer (2018):

[TABLE]

where $a=q,g$, $P=(p+p^{\prime})/2$ is the average momentum between initial and final states, and $\Delta=p^{\prime}-p$ is the momentum transfer. The following sum rules follow from conservation of energy and momentum:

[TABLE]

Since the NJL model does not contain gluons, the former sum rule translates for us to $\bar{c}_{q}(t)=0$. Another rule, following from a low-energy pion theorem Novikov and Shifman (1981); Voloshin and Zakharov (1980); Polyakov and Schweitzer (2018), is:

[TABLE]

This is a rule that holds exactly only in the chiral limit, and is a consequence of the pion being the Nambu-Goldstone boson of chiral symmetry breaking. Since the NJL model is a model of dynamical chiral symmetry breaking, we should observe Eq. (31) when taking the pion mass to zero. Moreover, one expects a quantity close to (although not exactly) $-1$ even at the physical pion mass.

The most general form the spin-one EMT can take is Cosyn et al. (2019); Polyakov and Sun (2019):

[TABLE]

where $A^{\{\mu\nu\}}=\frac{1}{2}(A^{\mu\nu}+A^{\nu\mu})$ and $A^{[\mu\nu]}=\frac{1}{2}(A^{\mu\nu}-A^{\nu\mu})$. The following sum rules follow from energy-momentum conservation, and the assumption that the symmetric and antisymmetric components of the EMT are separately conserved333 This assumption is true of the gauge-invariant kinetic EMT used in Cosyn et al. (2019), and is also of the canonical EMT of the quark field in the NJL model. However, the symmetric and antisymmetric components of the canonical EMT of the gluon field are not separately conserved in QCD. One should note that the gauge-invariant kinetic form of the gluon EMT does not give rise to a graviton vertex that satisfies the gravitational WTI of Eq. (18) that we employ in this work, while the canonical gluon EMT does. :

[TABLE]

but since there are no gluons in the NJL model, each of these form factors should be identically zero for the quarks. We have two additional sum rules, the first also following from energy-momentum conservation, and the latter from angular momentum conservation:

[TABLE]

Lastly, there is a correspondence between the GFF $\mathcal{G}_{10}(t)$ for a spin-one hadron is related to its isoscalar axial form factors:

[TABLE]

This correspondence is a consequence of the QCD equation of motion Leader and Lorcé (2014):

[TABLE]

However, Eq. (37) also holds in the NJL model (see App. C for a proof), and therefore Eq. (36) is also expected to hold in the NJL model.

III.4 Energy-momentum tensor of free quarks

Although quarks are not asymptotic states in the NJL model owing to the introduction of an infrared regulator, we may ask what the EMT and gravitational form factors of a free dressed quark would look like if we were to remove this regulator. The most general from that the EMT of a spin-half particle can take is Lorcé et al. (2019):

[TABLE]

Using the quark-graviton vertex in Eq. (25) and placing the quark on-shell gives, once the quark flavors have been summed over, $A_{Q}(t)=1$, $B_{Q}(t)=0$, $C_{Q}(t)$ as given in Eq. (26), $\bar{c}_{Q}(t)=0$, and $D_{Q}(t)=-1+D_{Q}^{\prime}(t)$, with $D_{Q}^{\prime}(t)$ given in Eq. (27). This allows the functions $C_{Q}(t)$ and $-1+D_{Q}^{\prime}(t)$ we found by solving the BSE to be interpreted as gravitational form factors of dressed quarks.

We have plotted in Fig. 4 the two non-trivial GFFs of a dressed quark, leaving out the trivial $A_{Q}(t)=1$ and $B_{Q}(t)=0$, using the NJL model parameters given in Tab. 1. In Fig. 4, both $C_{Q}(t)$ and $D_{Q}(t)$ can be seen to approach their elementary values of [math] and $-1$, respectively. It should be noted that $C(t)=0$ for an elementary free Dirac particle Hudson and Schweitzer (2018), and that a finite outcome for this form factor can only be produced by interaction dynamics.

IV Axial form factors

In order to check Eq. (36) within the NJL model, we must also calculate the axial form factors of the rho meson. Since gravitation is an isoscalar interaction, we look specifically at the isoscalar axial vector current. The matrix element of the axial vector current takes the following most general form for a spin-one particle Frederico et al. (1992); Berger et al. (2001):

[TABLE]

In order to calculate the matrix element in this equation, we must evaluate the triangle diagrams in Fig. 5 with the axial vector vertex for a dressed quark, and this vertex itself must be found through an inhomogeneous Bethe-Salpeter equation. In this section, we will find the axial form factors of the rho by first finding the requisite axial vector vertex, and in the process explore its properties and demonstrate that it satisfies the axial Ward-Takahashi identity.

IV.1 Axial vector and pseudoscalar vertices

We begin by finding the dressed axial vector and pseudoscalar vertices in the NJL model. The relevant currents are isoscalar, and the driving terms for the axial vector and pseudoscalar vertices are respectively $\gamma^{\mu}\gamma_{5}$ and $\gamma_{5}$. We denote the dressed vertices $\Gamma_{5}^{\mu}(p^{\prime},p)$ for the axial vector and $\Gamma_{5}(p^{\prime},p)$ for the pseudoscalar current.

The most general forms that the dressed vertices can take are:

[TABLE]

where $a_{i}(t)$ and $g_{i}(t)$ are off-shell form factors to be solved for algebraically. The resulting solutions to the two BSEs are:

[TABLE]

where the bubbles are given in App. A, and

[TABLE]

IV.2 Axial vector Ward-Takahashi identity

The axial vector and pseudoscalar vertices we have found should satisfy an axial vector Ward-Takahashi identity, which can be stated as Adler (1969):

[TABLE]

Here, $m$ is the current quark mass rather than the dressed quark mass. One can immediately observe that:

[TABLE]

meaning the axial vector WTI requires that:

[TABLE]

This can be proved, but requires a little work to show. The equality between the second and third expressions in Eq. (51) can be shown using a little algebra using $\Pi_{PA}(t)=-\Pi_{AP}(t)$ and $\Pi_{AA}^{(L)}(t)=-2\Pi_{PA}(t)$. The equality between the first and second expressions requires:

[TABLE]

The remaining bubbles can be found to evaluate to:

[TABLE]

where we leave the regularization scheme unspecified and implicit, except to assume that several basic operations are allowed, namely: (1) differentiation and integration with respect to variables other than $k$ commutes with the $k$ integral, and (2) one can cancel factors of $k^{2}-B(M,x,t)$ between the numerator and denominator and still obtain the same result. These restrictions, it should be noted, are met by the proper time regularization scheme we employ. With these rules, we can find that:

[TABLE]

Comparison with the gap equation, Eq. (2), gives us:

[TABLE]

and therefore the axial vector WTI requires that:

[TABLE]

which holds true if and only if $G_{\eta}=G_{\pi}$. This requirement is not surprising, since $\mathrm{U}(1)_{A}$ invariance is the symmetry responsible for the partial conservation of isoscalar axial symmetry, and therefore for the isoscalar axial vector WTI. It was remarked in Cloët et al. (2014) that $G_{\eta}=G_{\pi}$ is a necessary condition for the NJL model Lagrangian (1) to satisfy $\mathrm{U}(1)_{A}$ symmetry, and the necessity of this equality in the final step of our proof is simply a manifestation of that fact.

V Results

We now present results for the pion and rho meson EMTs as calculated in the two-flavor NJL model. Since the EMT is an isoscalar quantity, and since we take $m_{u}=m_{d}$ in this work, the up and down quark contributions to the EMT are equal, and we present the meson EMT obtained after summing over both quark flavors.

V.1 Pion and sigma

The GFFs of the pion and sigma meson have been plotted in Fig. 7. One can see immediately that the expected constraints appear to hold. In particular, we have $A_{\pi,\sigma}(0)=1$ and $C_{\pi}(0)\approx-1$. Not included in the plot are the results $\bar{c}_{\pi,\sigma}(t)=0$, which are found to hold exactly. In addition, we find that in the chiral limit, $C_{\pi}(0)=-1$ exactly, satisfying the low-energy pion theorem of Eq. (31).

It’s worth stressing the importance of the “quark D term” [$C_{Q}(t)$ in Eq. (25), not $D_{Q}(t)$] in satisfying the low-energy pion theorem. This has implications for the calculation of generalized parton distributions (GPDs). The leading-twist GPDs $H_{\pi}^{q,g}(x,\xi,t)$ are related to GFFs through the second Mellin moment444 This is true for the gluon GPD is the Ji convention Ji (1998) is used. If the Diehl convention Diehl (2003) is used, the gluon GFFs come from half the first Mellin moment of the GPD. , viz.,

[TABLE]

which is a specific case of the more general property that Mellin moments of GPDs produce matrix elements of local operators.

If a bare non-local operator is used to calculate the GPDs, then its moments can produce only matrix elements of bare local operators. For the second moment in particular, this will give the “$++$” component of the bare three-point graviton vertex. As a consequence, if one does not dress the non-local operator defining the GPDs, an incorrect value will be inferred for $C_{a}(t)$.

To illustrate this, we define “partially dressed” gravitational form factors as those that incorporate the vacuum condensate diagram, but for which $C_{Q}(t)$ is set to zero. We note that we need to include the vacuum condensate diagram for consistency, as otherwise the gravitational WTI is violated. Additionally, this diagram has no bearing on the Mellin moments of leading-twist GPDs, since its contribution contains a factor $g^{\mu\nu}$, and $g^{++}=0$. On the other hand, $C_{Q}(t)$ should appear in the Mellin moment with a factor $\xi^{2}$, so we expect $C_{\pi}(t)$ to be incorrect with this term missing.

In Fig. 8, we show the difference between $C_{\pi}(t)$ as calculated in the NJL model with and without the quark D-term $C_{Q}(t)$ folded in. The contrast is actually stark, with about two-thirds of the pion D-term coming from the quark dressing. We repeat and reemphasize that the “partially dressed” curve in this plot is the $C_{\pi}(t)$ that one would find from taking the second Mellin moment of the leading-twist pion GPD if it were calculated using a bare non-local operator, and that this speaks to the necessity of dressing operators in GPD calculations.

We illustrate this point further by finding analytic expressions for $C_{Q}(0)$ and $C_{\pi}(0)$ within proper time regularization. We find:

[TABLE]

where the $\tau$ integration is from $\Lambda_{\mathrm{UV}}^{-2}$ to $\Lambda_{\mathrm{IR}}^{-2}$. One can indeed see that $C_{\pi}(0)=-1$ exactly in the chiral limit. Additionally, we instead get exactly $-\frac{1}{3}$ if the contribution from $C_{Q}(0)$ is neglected—meaning a full two-thirds of the expected value come from the quark dressing.

The necessity for fulling dressing the relevant operator to observe the low-energy pion theorem—either the three-point graviton vertex or the bilocal light cone operator—draws an analogy to other results of dynamical chiral symmetry breaking, such as the dressing of quark mass. These non-perturbative phenomena both require a self-consistent solution of the relevant Dyson-Schwinger equations to manifest, and both result in significant changes to physical, observable quantities that are not constrained by conservation laws. This can be contrasted with electric charge, which remains unaltered between the bare electromagnetic vertex and the dressed vertex found by solving the inhomogeneous BSE.

In contrast to the pion, the sigma meson is not a Nambu-Goldstone boson and its “D-term” is not constrained in any a priori fashion. One can observe that $C_{\sigma}(0)=-2.27\not\approx-1$.

On another point, it is worth mentioning that the vanishing of $\bar{c}_{\pi,\sigma}(t)$ is possible only when all three diagrams have been included. Each of the diagrams on its own makes a non-zero contribution to the spin-zero meson EMT. The triangle diagrams contribute $\frac{1}{2}Z_{\pi}\Pi_{PP}(m_{\pi}^{2})$ or $\frac{1}{2}Z_{\sigma}\Pi_{SS}(m_{\sigma}^{2})$ each to $\bar{c}_{\pi}(t)$ or $\bar{c}_{\sigma}(t)$, respectively (with direct application of the WTI of Eq. (18) being the most straightforward way to find this result), while the bicycle diagram contributes 2G_{\pi}Z_{\pi}\Big{(}\Pi_{PP}(m_{\pi}^{2})\Big{)}^{2} or -2G_{\sigma}Z_{\sigma}\Big{(}\Pi_{SS}(m_{\sigma}^{2})\Big{)}^{2} to the same form factor. The three diagrams sum to

[TABLE]

which both vanish because of the pion and sigma mass-shell conditions in Eqs. (7,8).

V.1.1 Spin-zero meson mass radius

With the EMT of the spin-zero mesons in hand, it is also possible to study their static mechanical properties. It is conventional in much of the literature Polyakov (2003); Lorcé et al. (2018) to define a Breit frame EMT through:

[TABLE]

where the requirement that $\mathbf{P}=\mathbf{p}^{\prime}-\mathbf{p}=0$ gives us $P^{0}=\sqrt{m_{\pi}^{2}+\bm{\Delta}^{2}/4}$ and $\Delta^{0}=0$. One can then evaluate moments of $\langle T_{\mu\nu}\rangle_{\mathrm{Breit}}(\mathbf{r})$ to obtain static properties such as the mass radius.

Interpretation of Eq. (63) as an actual spatial distribution has been called into question Miller (2019). Despite this, one could still formally define a mass radius in terms of the Breit frame “density,” as is often done with the proton charge radius. Nonetheless, we run into problems for the pion. The limit $\bm{\Delta}\rightarrow 0$ must be taken when evaluating multipole moments, which requires the existence of a rest frame for the particle. Thus, Eq. (63) cannot be used in the chiral limit.

One may alternatively define a two-dimensional spatial distribution using light cone quantization Dirac (1949), where the spatial dimensions are the transverse light cone coordinates. The transverse spatial distribution is then:

[TABLE]

where the requirement $\Delta^{+}=0$ implies $P^{+}=p^{+}=p^{\prime+}$. Since in light cone quantization $P^{+}$ is kinematic and $P^{-}$ is dynamical Dirac (1949); Brodsky et al. (1998)555 The roles of $P^{+}$ and $P^{-}$ in Dirac (1949) are the opposite as in this work. , $P^{+}$ has no dependence on $\bm{\Delta}$, and thus no dependence on $t$, by contrast to $P^{0}$ in the Breit frame. We shall see the significance of this presently.

The mass radius can be found as the mean value of either $\mathbf{r}^{2}$ weighted by $\frac{1}{2P^{0}}\langle T^{00}\rangle_{\mathrm{Breit}}(\mathbf{r})$, or $\mathbf{r}_{\perp}^{2}$ weighted by $\frac{1}{2P^{+}}\langle T^{++}\rangle_{\mathrm{LC}}(\mathbf{r}_{\perp})$. We find:

[TABLE]

These radii differing (beyond a factor of the number of spatial dimensions) ultimately amounts to the mass form factor being different in the Breit frame and on the light cone, respectively $P^{0}A(t)$ and $P^{+}A(t)$. In the relevant frames, $P^{0}$ has $t$ dependence while $P^{+}$ does not. This is analogous to the electric charge distribution of the nucleon Miller (2007) or deuteron Carlson and Vanderhaeghen (2009) differing between the Breit frame and the light cone. In the case of the nucleon, the Sachs form factor $G_{E}(t)$ describes the electric charge distribution in the Breit frame, while the Dirac form factor $F_{1}(t)$ does instead on the light cone.

A remarkable property of the Breit frame radius (65) is that it remains finite even for point particles, for which $\frac{\mathrm{d}A(t)}{\mathrm{d}t}=0$. By contrast, the light cone radius (66) is zero for point particles. The former can be understood as owing to spatial distributions not being invariant under Lorentz boosts, while the latter is due to the fact that transverse boosts in light cone coordinates are Galilean Brodsky et al. (1998). In more detail, spatial distributions in the Breit frame are found by integrating over a collection of reference frames where the system of interest is in motion by different amounts, but without applying any corrections to counteract Lorentz contractions. These relativistic corrections are intrinsically accounted for by using light cone coordinates, however.

Because of the $\mathcal{O}(m_{\pi}^{-2})$ term in Eq. (65), the Breit frame radius of the pion blows up in the chiral limit. On the other hand, Eq. (66) remains finite at zero pion mass. Both equations can be used at physical pion mass, but produce staggeringly different values for the mass radius. Numerical values can be found in Tab. 2. For the empirical value of the pion mass radius, we look to the extraction in Kumano et al. (2018), where a dispersive analysis of KEKB data for $\gamma^{*}\gamma\rightarrow\pi^{0}\pi^{0}$ was done to extract gravitational form factors. The formula used in Kumano et al. (2018) was the same as our Eq. (66) for the squared light cone pion mass radius, but with a factor 6 instead of 4. We have thus scaled down the range of $0.32\sim 0.39$ fm reported in Kumano et al. (2018) by a factor $\sqrt{2/3}$. The NJL model result for this radius agrees with the empirical range, but falls on the low end.

In Ref. Cloët et al. (2014), the Breit frame charge radius of the pion was found to be 0.62 fm, which is significantly smaller than the Breit frame mass radius we have found. However, as we have discussed, the Breit frame mass radius made artificially large by $\mathcal{O}(m_{\pi}^{-1})$ terms that are not present in the light cone mass radius. The pion charge radius of Ref. Cloët et al. (2014), when scaled by $\sqrt{2/3}$, gives a light cone charge radius of 0.51 fm, which is instead larger than the light cone mass radius. The Breit frame and light cone prescriptions for radii thus suggest strongly divergent pictures of the relative distribution of mass and charge in the pion. Since doubt has been cast on the interpretation of the Breit frame radius as the moment of an actual density Miller (2019), the picture painted by light cone coordinates seems more plausible.

Because the form factor $A_{\sigma}(t)$ falls off faster than $A_{\pi}(t)$ (see Fig. 7), the light cone mass radius of the sigma is larger than that of the pion. On the other hand, because of its greater mass, the sigma meson has a smaller Breit frame radius than the pion. There is a stark difference between the results in these two frames, again demonstrating the magnitude and importance of correctly accounting for relativistic effects—by, for instance, using light cone coordinates in defining spatial distributions.

Finally, we remark again that getting the correct value of $C_{\pi}(0)$ by fully dressing the quark-graviton vertex is vital here. If one neglects the “quark D-term,” or equivalently obtains the GFFs through Mellin moments of bare GPDs, so that $C_{\pi}(0)\approx-\frac{1}{3}$, then one finds $\langle r^{2}\rangle_{\mathrm{Breit}}<0$ at the physical pion mass—an obvious absurdity that violates the weak energy condition Hawking and Ellis (2011).

V.2 Rho meson

The rho meson, as a spin-one hadron, has many more GFFs than the pion. There are 11 GFFs total, with 4 of these being non-conserved Cosyn et al. (2019). As required from the lack of gluons in the NJL model, the four non-conserved GFFs vanish: $\mathcal{G}_{7}=\mathcal{G}_{8}=\mathcal{G}_{9}=\mathcal{G}_{11}=0$. There are thus seven non-zero GFFs to consider.

Three of the rho GFFs have direct analogues to the pion GFFs, in particular, $\mathcal{G}_{1}(t)\sim A_{\pi}(t)$, $\mathcal{G}_{3}(t)\sim C_{\pi}(t)$, and $\mathcal{G}_{8}(t)\sim-2\bar{c}_{\pi}(t)$. The behavior of these form factors is remarkably similar to those of the pion. Firstly, we find $\mathcal{G}_{1}(0)=1$, as is required by momentum conservation. Additionally, we curiously find $\mathcal{G}_{3}(0)\approx-1$, and even more curiously find that this becomes exactly $-1$ in the chiral limit. Since the rho meson is not a Nambu-Goldstone boson, we do not know of any theorems requiring that this be the case, as we did for $C_{\pi}(0)$. Lastly, as with $\bar{c}_{\pi}(t)$ for the pion, we find that including the contribution of the bicycle diagram (rightmost diagram in Fig. 3) is necessary for $\mathcal{G}_{8}(t)$ to fully vanish.

Among the new conserved form factors, $\mathcal{G}_{5}(t)$ describes the spatial distribution of total angular momentum. An angular momentum sum rule Abidin and Carlson (2008) requires that $\mathcal{G}_{5}(0)=2$, and we satisfy this sum rule in the NJL model. The non-conserved form factor $\mathcal{G}_{7}(t)$ also contributes to this spatial distribution, but vanishes in the NJL model. The remaining form factors contribute only to higher multipole moments of the energy-momentum tensor.

Of the seven non-zero GFFs, the six $\mathcal{G}_{1-6}(t)$ appear in the symmetric component of the energy-momentum tensor. These six form factors are of special phenomenological interest, since they can be found from second Mellin moments of leading-twist generalized parton distributions Cosyn et al. (2018, 2019), which can be measured in hard exclusive reactions such as DVCS and DVMP. These six GFFs have been plotted in Fig. 9.

In Cosyn et al. (2019), the Breit frame multipole moments of the spin-one EMT were found. The mean squared mass radius and gravitational quadrupole moment are:

[TABLE]

where we have neglected the non-conserved form factors, since they are zero in the NJL model. Since the Breit frame density cannot be literally interpreted as an actual spatial density Miller (2019), we also determine the light cone transverse mass radius as a point of contrast, which we find to be:

[TABLE]

Since the light cone density is a two-dimensional quantity, we will defer exploration of the quadrupole moment—both quantitative and conceptual—to a future work on the light cone interpretation of the EMT.

The numerical values have been computed and tabulated in Tab. 3, along with the equivalent electric (Coulomb) quantities Cloët et al. (2014). The quadrupole moments are remarkably close when comparable units are used, but the electric charge radius is larger than the mass radius, for both the Breit frame and light cone radii. This suggests a highly inhomogeneous distribution of electric charge in the rho meson, as could occur (for instance) in a configuration with a small negatively-charged core surrounded by a shell of positive charge. Curiously, this occurs even for a positive rho meson—consisting of positively charged up and anti-down quarks—but we remark that the dressed quarks themselves have spatially extended electric charge distributions.

The conserved, non-zero GFF $\mathcal{G}_{10}(t)$ appears in the antisymmetric component of the EMT. In both QCD and the NJL model, this GFF has a special relationship with the axial form factors, given in Eq. (36). To demonstrate this correspondence, we plot the relevant form factors in Fig. 10. Because of this relation, $\mathcal{G}_{10}(t)$ encodes information about the distribution of quark spin, with the quantity $-\frac{1}{2}\mathcal{G}_{10}(0)=0.526$ giving the proportion of the rho meson’s total angular momentum carried by quark spin. If we take $G_{f}=0$, we instead get $0.556$, a number consistent with previous work on the rho meson in the NJL model Ninomiya et al. (2017) that implicitly assumed $G_{f}=0$.

VI Conclusions

In this work, we have found the gravitational form factors appearing in the decompositions of the canonical energy-momentum tensor for the pion, sigma, and tho mesons in the NJL model. In the process of obtaining these, we proved a gravitational Ward-Takahashi identity for the canonical EMT of fields with arbitrary spin, which takes a simpler form than the equivalent WTI for the Belinfante EMT. This simpler WTI makes consistency cross-checks between the dressed graviton vertex and dressed propagator within model calculations easier, and will also extend to models with spin-one constituents.

We found that the non-linear four-fermi interaction of the NJL model entails a five-point gravitational vertex (with four quark lines and one graviton line) in addition to the usual three-point vertex (with two quark lines and one graviton). This finding is a result of the equivalence principle, and is necessary for conservation of energy and momentum to be observed. In solving the Bethe-Salpeter equation for the dressed three-point vertex, we found it was necessary to include contributions from the five-point vertex for the gravitational WTI to be satisfied. Additionally, the five-point vertex contributed to meson EMT calculations directly in the form of a bicycle diagram, which was needed in order for the non-conserved GFFs $\bar{c}(t)$ and $\mathcal{G}_{8}(t)$ to vanish.

The necessity of the five-point vertex suggests difficulties in the prospect of calculating gravitational form factors in QCD. The QCD Lagrangian contains quark-gluon, three-gluon, four-gluon, and (in a covariant gauge) ghost-gluon vertices in its Lagrangian, and the equivalence principle—encoded by the presence of $-g^{\mu\nu}\mathcal{L}$ in the EMT—requires that the graviton be able to couple directly to every one of these vertices. Inclusion of all these graviton interactions is likely necessary for energy-momentum conservation to be observed, just as inclusion of the five-point vertex was in the NJL model. Moreover, consistent solution of Dyson-Schwinger equations for the dressings of these vertices is likely necessary for the gravitational WTI to be satisfied.

The meson GFF results demonstrate several fascinating properties of not just the mesons themselves, but the field theoretical framework which gives rise to them. A low-energy pion theorem requires $C_{\pi}(0)\approx-1$, with the approximation becoming exact in the chiral limit. The NJL model correctly reproduces this (including the exactness in the chiral limit), and we found fully dressing the three-point graviton vertex by solving its inhomogeneous Bethe-Salpeter equation to be necessary to reproduce this behavior—an observation that has implications for model calculations of generalized parton distributions of the pion.

We also observed the importance of correctly accounting for relativistic effects when describing spatial properties of light mesons, especially the pion—and more especially when considering the chiral limit. The Breit frame mass radius of the pion is significantly larger than its transverse light cone radius, only the latter of which is even finite in the chiral limit. We concur with previous literature that light cone coordinates are necessary to meaningfully define spatial distributions. We find that the light cone mass for the pion predicted by the NJL model agrees with a phenomenological extraction from KEKB data.

In future work, we plan to extend the methods developed here to baryons. The results found here will be directly applicable within a quark-diquark model, where scalar and axial vector diquarks mimic closely the structure of pions and rho mesons. Additionally, the gravitational WTI we derived can be used as a consistency cross-check for the off-shell diquark-graviton vertex.

Acknowledgements

We would like to thank Rafael Badui, Wim Cosyn, Sabrina Cotogno, and Cédric Lorcé for illuminating discussions that helped contribute to our investigation. This work was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, contract no. DE-AC02-06CH11357. AF was supported by an LDRD initiative at Argonne National Laboratory under Project No. 2017-058-N0.

Appendix A Bubbles in the NJL model

The bubbles are defined using the following convention:

[TABLE]

where $\Omega_{X}$ and $\Omega_{Y}$ represent particular vertices that appear in the bubble diagrams. A diagrammatic depiction of the bubble is given in Fig. 11, where the factor $i^{2}$ from the propagators cancels the factor $(-1)$ appearing because the diagram has a closed fermion loop. An overall factor $i$ is included in the definition (70) since this factor typically appears in diagrammatic equations involving bubbles as sub-diagrams, and this factor additionally makes the bubbles purely real below the two-particle production threshold. Specific bubbles that appear in this work are defined as follows:

[TABLE]

Appendix B Proof of the gravitational Ward-Takahashi identity

We here prove the gravitational Ward-Takahashi identity given in Eq. (18), first proved for spin-zero fields in Brout and Englert (1966), is true for fields with arbitrary spin, provided the graviton couples to the canonical energy-momentum tensor.

We must start with a definition of the fully dressed 3-point gravitational vertex function $\Gamma^{\mu\nu}(p^{\prime},p)$. The definition we use is:

[TABLE]

where $p$ and $p^{\prime}$ are respectively the initial and final momentum of a quantum of a field $\phi$, $S$ is the fully-dressed propagator of said field, $T^{\mu\nu}$ is the canonical energy-momentum tensor (EMT), and we explicitly notate the internal degrees of freedom of $\phi$ with indices $r,s$.

We begin by contracting the left-hand side of Eq. (76) with $q_{\mu}$. Let us use $L^{\mu\nu}_{rs}(p^{\prime},p,q)$ as shorthand for the entire LHS, for better use of space. We have:

[TABLE]

Now, in differentiating this time-ordered product, we note that $\partial_{\mu}T^{\mu\nu}(0)=0$ as a consequence of Noether’s theorem. (We note that differentiation with respect to the first index is crucial here, as conservation with respect to the second index is not guaranteed.) Thus, any $x$ dependence in the time-ordered product can only come from the time-ordering itself. We have, explicitly:

[TABLE]

where we have used the canonical commutation relation

[TABLE]

from which $P^{\mu}=\int\mathrm{d}^{3}x\,T^{0\mu}(x)$ follows. (Note that such a commutation relation does not hold for $T^{\mu 0}(x)$, giving us a second point at which privileging the first index is important.) From the definition of the propagator, we obtain:

[TABLE]

This enables us to find:

[TABLE]

At this point, we compare the RHS of Eq. (76) and integrate over $q$. We have:

[TABLE]

Now, by cancelling out the factors $i$, multiplying by the inverse propagators, and dropping the $r$ and $s$ indices, we get Eq. (18), as required.

Appendix C Proof of axial-gravitational correspondence in NJL model

To prove Eq. (37) holds in the NJL model, we follow the derivation in Sec. 3.2 of Leader and Lorcé (2014) for QCD. Using the NJL model Lagrangian (1), we obtain the following equations of motion:

[TABLE]

From here, we use Eqs. (129) of Leader and Lorcé (2014) with $A_{\mu}=0$, which can be stated:

[TABLE]

With this, we find that:

[TABLE]

What remains to be shown is that $\overline{\psi}[\hat{M},\sigma^{\mu\nu}]\psi=0$. Using the identities:

[TABLE]

we find:

[TABLE]

Finally, sandwiching this between $\psi$ and $\overline{\psi}$, we get

[TABLE]

an expression which is clearly both symmetric and antisymmetric under the swap $(\mu\leftrightarrow\nu)$, and which is therefore zero. We have thus proved Eq. (37) holds in the NJL model.