Within the framework of string theory, a number of new fields are possible correcting the Einstein-Hilbert action, including a Kalb-Ramond two-form field. In this work we derive explicitly first order relativistic corrections to conservative dynamics with a Kalb-Ramond field, using the effective field theory approach. The resulting additional terms in the Lagrangian governing conservative binary dynamics are presented explicitly.

§ INTRODUCTION

String theory postulates a theory of quantum gravity, which generically includes extra dimensions. This leads to modifications of Einstein's theory of general relativity. However, it is not known at which scale these modifications will manifest themselves, and this depends heavily on which “string compactification" is employed.  At lower energies gravitational effects in the extra dimensions will give rise to new fields called moduli fields, generically present in most models derived from string theory. These fields, and in particular their axionic components and the Peccei–Quinn axion <cit.>, have been suggested as natural candidates for light dark matter, where a plethora of physical models and searches have been suggested for their observations <cit.>.
The effects of light particles and QCD axions in gravitational mergers have also been studied in the literature before, see for example <cit.>.

In addition to moduli fields, string theory and more generally supergravity theories, are equipped with a Kalb-Ramond two-form field <cit.>. This field was originally introduced as a classical generalisation of the electromagnetic interaction of point particles to one-dimensional strings in the Feynman-Wheeler picture <cit.>.  The Kalb-Ramond field is now understood as a generic feature of supergravity theories, forming part of the gravitational supermultiplet.[Here, we are explicitly focussing only on contributions from the Kalb-Ramond field, ignoring the dilaton and other supergravity effects. Generic four-dimensional solutions with such fields turned on were found in <cit.>.] Classically upon dualising, the Kalb-Ramond field may be thought of as a “fundamental string axion” <cit.>, although the physical and quantum behaviours of this fundamental axion are slightly different to the QCD axion. For example, the fundamental Kalb-Ramond axion is required to be massless at the level of the classical supergravity action, as its coupling is dictated by classical gauge symmetry.

In string theory the Kalb-Ramond two-form is responsible for mediating the force between fundamental strings, whose natural energy scale is the string scale. Being a massless mediator, the classical force is expected to be long range. However, as the couplings are string scale suppressed, any direct measurement of such interactions seem a long way off. The recent discoveries by NANOGrav are however interesting in this regard, as cosmic superstrings provide a potential explanation of the data <cit.>. As fundamental strings, these would then also provide large sources of potentially detectable Kalb-Ramond charge <cit.>.

Another possibility is that a large number of strings can somehow amalgamate to produce a compound effect which can be measured. One such scenario where this can potentially happen is in the case of fuzzballs, see <cit.> and references therein. This string-inspired alternative to black holes has a large number of strings condense at the horizon of the black hole, where space time effectively ends. The accumulated Kalb-Ramond charge may then have a chance of being observed, and it is this coupling we wish to model in the present paper. We'll do this by classically modelling the Kalb-Ramond field as an axionic scalar, where astrophysical objects can be charged under this field. At the quantum level, the Kalb-Ramond field should be treated as a two-form, whose interactions are rather more subtle and involved. We will leave the quantum analysis for future work.

In this work, we derive classical corrections to the gravitational potential for a binary system. Recent years have seen substantial improvements in constraining binary motion for relativistic systems. One of these is the pulsar timing with radio telescopes of the double pulsar <cit.>. Another is the constraints from gravitational wave signals from binary inspiral <cit.>. These are expected to improve by orders of magnitude over the coming decades <cit.>.

To derive the corrections due to the Kalb-Ramond field, we employ here the effective field theory formalism <cit.>. This formalism allows an explicit separation of scales and has been used in a number of other studies, see e.g. <cit.> and others. Our work here serves to illustrate the utility of this formalism for studying additional fields from particle physics and including them in gravitational interactions. For completeness we also explicitly include pure graviton interactions at the same order. As noted, we restrict ourselves to tree-level computations and do not include quantum effects in our analysis. Our main result is equation (<ref>). This differs qualitatively from previous analysis involving axionic scalars <cit.>. In particular, the Kalb-Ramond and classical General Relativity corrections have very similar $\tfrac{1}{r}$-behaviours in the $\tfrac{v}{c}$ expansion. We speculate how this relates to understanding the Kalb-Ramond field in the context of generalised geometry,[For an introduction to the mathematical framework of generalised geometry, see <cit.>, and for applications in supergravity, see e.g. <cit.> with references therein.] where the metric and Kalb-Ramond field are put on the same footing in the generalised metric.

§ EFFECTIVE ACTION FOR THE KALB-RAMOND FIELD

Our aim is to calculate the first-order relativistic correction to the interaction of two bodies charged under the Kalb-Ramond field. We therefore study a model in which the Kalb-Ramond field is minimally coupled to gravity alone.

The Kalb-Ramond field $B_{\mu\nu}$ is a two-form on spacetime with a field strength $H\ind{_{\mu\nu\rho}}$ given by
H_μνρ ≡B_μνρ = B_νρ,μ + B_ρμ,ν + B_μν,ρ

The action and equations of motion (EoM) for the metric including the Kalb-Ramond field are[We here define the action as $S = \int L \dd{t} = \int \mathcal{L} \dd[3]{x} \dd{t}$, thus we may write $\int \mathcal{L} \dd[4]{x} = cS $ for notational convenience.]

\begin{align}
& cS = \int \dd[4]{x} \sqrt{-g} \left(\frac{1}{2\kappa} R + \frac{1}{12} H\indices{_{\mu\nu\rho}} H\indices{^{\mu\nu\rho}} + \mathcal{L}_\text{matter} \right), \label{eq:Action:full} \\
& \nabla\ind{_\sigma} H\indices{^\sigma_{\mu\nu}} = 0, \label{eq:EoM:H}
\end{align}

§.§ A solution for the Kalb-Ramond field

The solution for the $H$-field, equation (<ref>), has been known for a long time, and can be found in <cit.>. Since $H\ind{_{\mu\nu\rho}}$ is exact by equation (<ref>), we make the ansatz that $H\ind{_{\mu\nu\rho}}$ can be described by the derivative of a scalar function with an antisymmetric tensor structure, such that

\begin{align}
H\indices{_{\alpha\beta\gamma}} = \varepsilon\indices{_{\alpha\beta\gamma\delta}} K\indices{^{;\delta}}\:, \label{eq:sol:H}
\end{align}

where $\varepsilon\indices{_{\alpha\beta\gamma\delta}} = \sqrt{-\abs{g}} \epsilon\indices{_{\alpha\beta\gamma\delta}}$ is the totally anti-symmetric volume form of four space-time dimensions, i.e. the Levi-Civita tensor.
This is a solution of equation (<ref>) based on the antisymmetry of the indices $\alpha \beta \gamma \delta$, and the symmetry of the derivative indices $\alpha \delta$: $\nabla_\alpha H\indices{^{\alpha}_{\beta\gamma}}=\nabla_\alpha \nabla_\delta \varepsilon\indices{^\alpha_{\beta\gamma}^\delta} K$.

The anzats (<ref>) can be further constrained by making sure $H\indices{_{\alpha\beta\gamma}}$ is exact by obeying equation (<ref>).

\begin{align}
\dd H\indices{_{\alpha\beta\gamma\delta}} = \varepsilon\indices{_{\alpha\beta\gamma\sigma}} K\indices{^{,\sigma}_\delta} -  \varepsilon\indices{_{\delta\alpha\beta\sigma}} K\indices{^{,\sigma}_\gamma} + \varepsilon\indices{_{\gamma\delta\alpha\sigma}} K\indices{^{,\sigma}_\beta} - \varepsilon\indices{_{\beta\gamma\delta\sigma}} K\indices{^{,\sigma}_\alpha} = 0.
\end{align}

Since $\varepsilon\indices{_{\alpha\beta\gamma\delta}}$ is fully antisymmetric, it follows that $\varepsilon\indices{_{\alpha\beta\gamma\delta}} \neq 0$ iff $\alpha\neq\beta\neq\gamma\neq\delta$. Hence, expressions of the form $\varepsilon\indices{_{\alpha\beta\gamma\sigma}} K\indices{^{,\sigma}_\delta}$ only get contributions when $\sigma=\delta$, all other terms cancelling, and giving all the terms a common factor of $K\indices{^{,\sigma}_\sigma}$.[This factor is not summed over $\sigma$.]

\begin{align}
\dd H\indices{_{\alpha\beta\gamma\delta}} & = \varepsilon\indices{_{\alpha\beta\gamma\delta}} K\indices{^{,\delta}_\delta} - \varepsilon\indices{_{\delta\alpha\beta\gamma}} K\indices{^{,\gamma}_\gamma} + \varepsilon\indices{_{\gamma\delta\alpha\beta}} K\indices{^{,\beta}_\beta} - \varepsilon\indices{_{\beta\gamma\delta\alpha}} K\indices{^{,\alpha}_\alpha} \nonumber \\
& = \varepsilon\indices{_{\alpha\beta\gamma\delta}} \left( K\indices{^{,\delta}_\delta} + K\indices{^{,\gamma}_\gamma} + K\indices{^{,\beta}_\beta} + K\indices{^{,\alpha}_\alpha} \right) = \varepsilon\indices{_{\alpha\beta\gamma\delta}} \dalembertian K = 0.
\end{align}

Thus the solution for the $H$-field is described by the gradient of a massless scalar field

\begin{align}
\begin{split}
H\indices{_{\alpha\beta\gamma}} & = \varepsilon\indices{_{\alpha\beta\gamma\delta}} K\indices{^{,\delta}},
\\ \quad \qq*{where} K\indices{_{,\mu}^{\mu}} & = 0 .\end{split} \label{eq:EoM:Lambda}
\end{align}

It also follows that the action can be rewritten to that of a massless scalar field[To obtain this result one needs the relation $\varepsilon\indices{^{\mu\nu\rho}_{\sigma}} \varepsilon\indices{_{\mu\nu\rho\delta}} = -6 g\ind{_{\sigma\delta}}$. This relation can be argued for from $\varepsilon\indices{^{\mu\nu\rho\sigma}} \varepsilon\indices{_{\mu\nu\rho\sigma}}=-4!$ in 4 space-time dimensions. Then, it follows straightforwardly that $\varepsilon\indices{^{\mu\nu\rho\sigma}} \varepsilon\indices{_{\mu\nu\rho\sigma}}=g\indices{^{\sigma\delta}} \left(\varepsilon\indices{^{\mu\nu\rho}_{\sigma}} \varepsilon\indices{_{\mu\nu\rho\delta}}\right) = -24 = g\indices{^{\sigma\delta}} \left(-6 g\ind{_{\sigma\delta}} \right)$. ]

\begin{align}
\begin{split}
H\indices{_{\mu\nu\rho}} H\indices{^{\mu\nu\rho}} & = g\ind{^{\sigma\alpha}} \varepsilon\indices{_{\mu\nu\rho\sigma}} \varepsilon\indices{^{\mu\nu\rho}_{\delta}} g\ind{^{\delta\beta}} K\indices{_{,\alpha}} K\indices{_{,\beta}} = -6 g\ind{^{\sigma\alpha}} g\indices{_{\sigma\delta}} g\ind{^{\delta\beta}} K\indices{_{,\alpha}} K\indices{_{,\beta}} = -6 K\indices{_{,\sigma}} K\indices{^{,\sigma}}, \\
\rightarrow \quad cS_B & = \frac{1}{12} \int \dd[4]{x} \sqrt{-\abs{g}} H\indices{_{\mu\nu\rho}} H\indices{^{\mu\nu\rho}} = -\frac{1}{2}\int \dd[4]{x} \sqrt{-\abs{g}} K\indices{_{,\sigma}} K\indices{^{,\sigma}}.
\end{split} \label{eq:action:K}
\end{align}

§.§ Relativistic expansion of the action

Expressing the Kalb-Ramond action in terms of the scalar field $K$, like in equation (<ref>), one can relativistically expand the action in the so-called post-Newtonian (PN) expansion. In post-Newtonian calculations, the metric is split into the usual weak-field approximation $g\ind{_{\mu\nu}} = \eta\ind{_{\mu\nu}} + \lambda h\ind{_{\mu\nu}}$, where $ \eta\ind{_{\mu\nu}} $ is the Minkowski metric, choosing suitable Cartesian coordinates such that $\eta\ind{_{\mu\nu}} = \text{diag}(-1,1,1,1)$. $h\ind{_{\mu\nu}}$ denotes the metric perturbations, and will also be referred to here as the graviton field. $\lambda$ is a scaling factor that governs the relativistic expansion, which will be shown to be $\propto \sqrt{G}$, the Newtonian gravitational constant.

Using (<ref>), the metric splitting gives an expansion of the Kalb-Ramond action in terms of the weak-(gravitational)field as

\begin{align}
cS_K & = -\frac{1}{2} \int \dd[4]{x} \sqrt{-\abs{g}} \left( \eta\indices{^{\mu\nu}} + \lambda h\indices{^{\mu\nu}} \right) K\indices{_{,\mu}} K\indices{_{,\nu}} \\
& = -\frac{1}{2} \int \dd[4]{x} \left( 1 + \frac{1}{2} \lambda h +\dots \right) \left( \eta\indices{^{\mu\nu}} + \lambda h\indices{^{\mu\nu}} \right) K\indices{_{,\mu}} K\indices{_{,\nu}} \\
& \approx -\frac{1}{2} \int \dd[4]{x} \left[ \left(1 + \frac{1}{2}\lambda h\right) K\indices{_{,\mu}} K\indices{^{,\mu}} + \lambda h\indices{^{\mu\nu}} K\indices{_{,\mu}} K\indices{_{,\nu}} \right]
\end{align}

The determinant of the metric was expanded using the following steps, assuming $\lambda h \ll 1$ and using $\eta\indices{_{\mu\nu}}$ for raising and lowering indices, which will be the case for the rest of this paper.

\begin{align}
\begin{split}
&\det(g\indices{_{\rho\nu}}) = \det( \eta\indices{_{\rho\mu}} g\indices{^\mu_{\nu}}) =  \det(\eta\indices{_{\rho\mu}}) \det(g\indices{^\mu_{\nu}}) = -\det(g\indices{^\mu_{\nu}}) = -\det(\delta\indices{^\mu_\nu} + \lambda h\indices{^\mu_\nu} )\\
& = -\exp(\ln( \det(\delta\indices{^\mu_\nu} + \lambda h\indices{^\mu_\nu} ) )) = -\exp(\tr( \ln(\delta\indices{^\mu_\nu} + \lambda h\indices{^\mu_\nu} ) )) \approx -\exp(\tr( \lambda h\indices{^\mu_\nu} )) \approx -1 - \lambda h.
\end{split} \label{eq:det(g)}
\end{align}

Then by also Taylor expanding the square root the following relation is also established, which appears in the action integral.

\begin{align}
\sqrt{-\det(g)} \approx \sqrt{1+\lambda h} \approx 1 + \frac{\lambda}{2} h + \order{\lambda^2 h^2}. \label{eq:sqrt(det(g))}
\end{align}

The form of the action (<ref>) is fairly similar to that which is investigated in <cit.> and <cit.> for massive axions coupling to gravity. The main difference is the absence of a mass term for the Kalb-Ramond field.

For the coupling term to sources like fuzzballs, we introduce effective Kalb-Ramond charges $q_a$, $p_a$, etc and generic interaction terms:

\begin{align}
S_{\text{int}K} = -\sum_{a} \int \frac{\dd{c\tau_a}}{c^2} \dd[3]{x} \left[ q_a K + p_a K^2 +\dots \right] \sqrt{ -g\ind{_{\mu\nu}} \dot{x}^\mu \dot{x}^\nu } \dirac{3}{\vb*{x}-\vb*{x}_a(\tau_a)}. \label{eq:action:interaction:LambdaPP}
\end{align}

The term in the square root is needed to make the point particle action reparametrization invariant.[Which is required by relativity since there are no privileged time coordinates, and to have classical particles have the correct number of degrees of freedom. This requirement is also equivalent to the condition $\dot{x}\ind{^\alpha} \pdv{L}{\dot{x}\ind{^\alpha}} = L$ <cit.>.]

§.§ The graviton and point-particle action

The linearized action of gravity is <cit.>

\begin{align}
cS_h = \int \dd[4]{x} \left( -\frac{1}{4} h\indices{_{\mu\nu,\rho}}h\indices{^{\mu\nu,\rho}} + \frac{1}{8} h\indices{_{,\mu}} h\indices{^{,\mu}} + \order{h^3}\right), \label{eq:Graviton:Action}
\end{align}

with the effective coupling to matter via its energy-momentum tensor $T\indices{^{\mu\nu}}$:

\begin{align}
cS_{\text{int}h} = \int \dd[4]{x} \frac{\lambda}{2} h\indices{_{\mu\nu}} T\indices{^{\mu\nu}}.
\end{align}

For objects influenced only by gravity, the action can be approximated at large scales by the free point-particle action from GR:

\begin{align}
& S_{fpp} = -\sum_a m_a c \int \dd{s_a} = -\sum_{a} \int m_a c \sqrt{-g\indices{_{\mu\nu}} \dd{x}^\mu_a \dd{x}^\nu_a } \nonumber \\ & \hspace{0.7cm} = -\sum_{a} m_a c \int \sqrt{1-\lambda h\indices{_{\mu\nu}} \frac{ \dot{x}_a^\mu(x^0) }{c} \frac{ \dot{x}_a^\nu(x^0) }{c}}  \gamma_a^{-1} \dd{x^0} \nonumber \\
& = \sum_{a} \int \frac{\dd[4]{x}}{c} \gamma_a^{-1} \left[ -m_a c^2 + \frac{\lambda}{2} h\indices{_{\mu\nu}} m_a \dot{x}\indices{^\mu} \dot{x}\indices{^\nu} + \frac{m_a\lambda^2}{8c^2} \left( h\indices{_{\mu\nu}} \dot{x}\indices{^\mu} \dot{x}\indices{^\nu} \right)^2 +\dots \right] \dirac{3}{ \vb*{x} - \vb*{x}_a(x^0) } .
\end{align}

The point-particle action can be further extended to include the effective coupling to the $K$-field most simply by adding the interaction term (<ref>)

\begin{align}
\begin{split}
S_{pp} = S_{fpp} + S_{\text{int}K} = \sum_a & \int \frac{\dd[4]{x}}{c} \gamma_a^{-1} \left(m_a c^2 + q_a K + p_a K^2 +\dots \right) \\
& \cdot \left[ -1 + \frac{1}{2} \lambda h\indices{_{\mu\nu}} \frac{\dot{x}\indices{^\mu}}{c} \frac{\dot{x}\indices{^\nu}}{c} + \frac{1}{8} \left( \lambda h\indices{_{\mu\nu}} \frac{\dot{x}\indices{^\mu}}{c} \frac{\dot{x}\indices{^\nu}}{c} \right)^2 +\dots \right] \dirac{3}{ \vb*{x} - \vb*{x}_a(x^0) }.
\end{split} \label{eq:ppAction full}
\end{align}

The different cross-terms in this action will here be represented by Feynman diagrams in figure <ref>-<ref>. Note that for multiple massive point particles the energy-momentum tensor is $T\ind{^{\mu\nu}}(x) = \sum_a \gamma_a^{-1} m_a \dot{x}\ind{^\mu} \dot{x}\ind{^\nu} \dirac{3}{\vb*{x} -\vb*{x}_a } $, which appears in the cross-term between the mass term and the linear term in $\lambda h\ind{_{\mu\nu}}$, hence the interaction term of gravity is to linear order $\frac{\lambda}{2} h\ind{_{\mu\nu}} T\ind{^{\mu\nu}}$.

§.§ Equations of motion for the gravity-Kalb-Ramond sector

Let $S_\text{total} = S_K + S_h + S_{pp} = \int \mathcal{L} \dd[4]{x}/c$. Variation of this action produces the equations of motion,

\begin{align}
& \partial_\mu \fdv{\mathcal{L}}{h\indices{_{\alpha\beta,\mu}}} = -\frac{1}{2} \dalembertian h\indices{^{\alpha\beta}} + \frac{1}{4}\eta\ind{^{\alpha\beta}} \dalembertian h = \fdv{\mathcal{L}}{h\indices{_{\alpha\beta}}} = \frac{\lambda}{2}\left(T\indices{^{\alpha\beta}} + t^{\alpha\beta}\right) - \frac{\lambda}{2} \left( \frac{1}{2} \eta\indices{^{\alpha\beta}} K_{,\sigma} K^{,\sigma} + K^{,\alpha} K^{,\beta} \right),\nonumber \\
& \Rightarrow \quad \dalembertian h\indices{_{\mu\nu}} = -\lambda P\indices{_{\mu\nu\alpha\beta}} \left( T\indices{^{\alpha\beta}} + t\indices{^{\alpha\beta}} - \frac{1}{2} \eta\indices{^{\alpha\beta}} K_{,\sigma} K^{,\sigma} - K^{,\alpha} K^{,\beta} \right) \label{eq:EoM:h} \\
\mathrm{and} \nonumber \\
&\partial_\mu \fdv{\mathcal{L}}{K\indices{_{,\mu}} } = -\partial_\mu  \left[ \left(1+\frac{\lambda}{2}h\right) K\indices{^{,\mu}} + \lambda h\indices{^{\mu\sigma}} K\indices{_{,\sigma}} \right] = \fdv{\mathcal{L}}{K} = -\overbrace{\sum_{a} \gamma_a^{-1} \left[ q_a + 2p_a K \right] \dirac{3}{\vb*{x} - \vb*{x}_a(x^0)} }^{=J_K} \nonumber \\
&\Rightarrow \quad \dalembertian K = J_K - \lambda \left( h\indices{^{\mu\nu}} K\indices{_{,\mu\nu}} + h\indices{^{\mu\nu}_{,\mu}} K\indices{_{,\nu}} \right) - \frac{\lambda}{2} \left(h\indices{_{,\mu}} K^{,\mu} + h\dalembertian K \right) \label{eq:EoM:K}
\end{align}

The $P\indices{_{\mu\nu\alpha\beta}}$ is the combination of $\eta\ind{_{\mu\nu}}$'s which when acted upon the first equality of equation (<ref>) simplifies it to just $-\frac{1}{2}\dalembertian h\ind{_{\mu\nu}}$.[For those familiar with the bar operator introduced in <cit.>, $P\indices{_{\mu\nu}^{\alpha\beta}}$ is the bar operator, which is also its own inverse.] The exact expression for $P\indices{_{\mu\nu\alpha\beta}}$ can be found in equation (<ref>). $t\indices{^{\alpha\beta}}$ represents the energy-momentum tensor of the graviton field itself, arising from cubic terms in the graviton action (<ref>). To leading order $t\indices{_{\alpha}^{\beta}} = \frac{1}{2} h\ind{_{\mu\nu,\alpha}}h\ind{^{\mu\nu,\beta}}$.

Both of the equations (<ref>) and (<ref>) can be solved by the method of Green's functions.

\begin{align}
& h\indices{_{\mu\nu}}(x) = -\lambda \int \dd[4]{y} \Delta(x^\sigma - y^\sigma) P\indices{_{\mu\nu\alpha\beta}} \left(T\indices{^{\alpha\beta}} + t\indices{^{\alpha\beta}} - \frac{1}{2} \eta\indices{^{\alpha\beta}} K_{,\sigma} K^{,\sigma} - K^{,\alpha} K^{,\beta} \right), \label{eq:FieldExpansion:h:general} \\
& K(x) = \int \dd[4]{y} \Delta(x^\sigma-y^\sigma) \left(  J_K - \lambda \left( h\indices{^{\mu\nu}} K\indices{_{,\mu\nu}} + h\indices{^{\mu\nu}_{,\mu}} K\indices{_{,\nu}} \right) - \frac{\lambda}{2} \left(h\indices{_{,\mu}} K^{,\mu} + h\dalembertian K \right) \right), \label{eq:FieldExpansion:K:general}
\end{align}

where $\Delta(x-y)$ is the Green's function of the d'Alembertian:

\begin{align}
& \Delta(x^\sigma-y^\sigma) = \int \frac{\dd[4]{k}}{(2\pi)^4} \frac{-1}{ k\indices{_\sigma} k\indices{^\sigma} } e^{ik_\mu \left( x^\mu - y^\mu \right)}, \label{eq:GreensFunction:d'Alembert} \\
& P\indices{_{\mu\nu\alpha\beta}} = \frac{1}{2}\left( \eta_{\mu\alpha} \eta_{\nu\beta} + \eta_{\mu\beta}\eta_{\nu\alpha} - \eta_{\mu\nu} \eta_{\alpha\beta} \right). \label{eq:Projection tensor}
\end{align}

To compute the leading order terms of the action in the post-Newtonian formalism interactions between different and similar fields can be neglected. Thus only the leading term of (<ref>) and (<ref>) contribute. Also to get the static, non-retarded, Newtonian potential the Green's function (<ref>) must be expanded in a power series and truncated at leading order. This effectively approximates it as the Green's function of the Laplacian operator $\nabla^2=\partial\ind{_i}\partial\ind{^i}$ (with corrections scaling as $(v/c)^{2n}$):

\begin{align}
& \Delta(x^\mu-y^\mu) = \int \frac{\dd[4]{k}}{(2\pi)^4} \frac{ e^{ik_\mu \left( x^\mu - y^\mu \right)} }{ -k_\sigma k^\sigma} = \int_k \frac{ e^{ik_\mu \left( x^\mu - y^\mu \right)} }{ k_0^2 - \vb*{k}^2 } = \int_k \frac{ e^{ik_\mu \left( x^\mu - y^\mu \right)} }{ -\vb*{k}^2 \left( 1 - \frac{k_0^2}{\vb*{k}^2} \right) } \nonumber \\
&\approx \int_k \frac{ e^{ik_\mu \left( x^\mu - y^\mu \right)} }{ -\vb*{k}^2 } \left( 1 + \left(\frac{k_0}{\vb*{k}}\right)^2 +  \left(\frac{k_0}{\vb*{k}}\right)^4 + \dots \right) = \Delta^{(0)}_\text{inst}(x-y) + \Delta^{(2)}_\text{inst}(x-y) + \dots  \label{eq:Propagator expansion}\\
& = \frac{-\delta (x^0 - y^0)}{4\pi \abs{\vb*{x} -\vb*{y}} } - \frac{ \vb*{\dot{x}} \vb*{\cdot} \vb*{\dot{y}} - (\vb*{\dot{x}}\vb*{\cdot} \vb*{\hat{r}})(\vb*{\dot{y}}\vb*{\cdot} \vb*{\hat{r}}) }{8\pi c^2 \abs{\vb*{x} -\vb*{y}} }\delta (x^0 - y^0) +\dots = \frac{-\delta(x^0 - y^0)}{4\pi r} \left( 1 - \frac{\eta}{2} \frac{\vb*{v}^2 - (\vb*{v\cdot \hat{r}})^2 }{ c^2 } + \dots \right). \nonumber
\end{align}

Diagrammatically the velocity expansion (<ref>) of the Green's function is represented by the $\otimes$-symbol on top of the propagator <cit.>. The number of $\otimes$ on the line indicates which order in the expansion is being considered, the leading order Newtonian propagator having none, the 1PN having one, etc.

Here the superscript (0) in $\Delta^{(0)}_\text{inst}(x-y)$ is to indicate this term scale as $(v/c)^0$ and thus belong in the 0PN order, while $\Delta^{(2)}_\text{inst}(x-y)$ scale as $(v/c)^2$ and belongs to the 1PN. It is given the subscript `inst' as this approximation makes the interaction instantaneous (via the factor of $\delta (x^0 - y^0)$).

With this all in mind leading-order terms of the action governed by graviton and Kalb-Ramond fields are

\begin{align}
&S_{pp}^{(0)} = \int \frac{\dd[4]{x}}{c} \sum_{a} \gamma_a^{-1} \left[ -m_ac^2 + \frac{m_a \lambda}{2} h_{00}(x) \dot{x}^0 \dot{x}^0 - q_a K(x) \right] \dirac{3}{ \vb*{x} - \vb*{x}_a(x^0) } \\
&= \int \frac{\dd[4]{x}}{c} \sum_{a} \gamma_a^{-1} \Biggl[ \sum_{b>a} \left\{ \frac{m_a \lambda}{2} \frac{\lambda m_b \gamma_b c^2}{8\pi \abs{\vb*{x} - \vb*{x}_b }} (\gamma_a c)^2 + q_a \frac{\gamma_b^{-1} q_b}{4\pi\abs{\vb*{x} - \vb*{x}_b }} \right\} -m_ac^2 \Biggr] \dirac{3}{ \vb*{x} - \vb*{x}_a(x^0) } \\
&= \int \left[ \frac{1}{2}m_1 v_1^2 + \frac{1}{2}m_2 v_2^2 + \frac{\lambda^2 m_1c^2 m_2c^2}{16\pi r} + \frac{q_1 q_2 }{4\pi r} \right] \dd{t}.
\end{align}

Thus we recognise the action of point particles in a Newtonian gravitational potential, also interacting with a Coulomb-like potential. The Newtonian limit implies

\begin{align}
\lambda^2 = \frac{16\pi G}{c^4} = \frac{16\pi}{ F_\text{Pl}}, \qq{$F_\text{Pl}$ being ``the Planck force''.}
\end{align}

§ EXPANDING THE ACTION INTO THE 1PN ORDER

§.§ Velocity corrections to the 0PN diagrams

Feynman diagrams of the 0PN potential of GR (fig:Feynman:H-type:0PN) and Kalb-Ramond (fig:Feynman:H-type:0PN:K), with 1PN velocity corrections (fig:Feynman:H-type:VelocityExp)-(fig:Feynman:H-type:ox) and (fig:Feynman:H-type:VelocityExp:K)-(fig:Feynman:H-type:ox:K). The wavy lines represent graviton propagators, the dashed lines Kalb-Ramond propagators, and the solid lines represent non-propagating, point-particle, world-line sources <cit.>.

In order to obtain the point particle action contribution at next to leading order in $(v/c)^{2n}$, the so-called 1PN order, both velocity modifications to the previous potentials as well as higher order terms of the potential must be accounted for. Diagrammatically the velocity corrections to the leading order term can be expressed using the Feynman diagrams of figure <ref>. Figure <ref> and <ref> represent the 0PN contributions already calculated

\begin{align}
V_{\ref{fig:Feynman:H-type:0PN}} &= -\frac{Gm_1 m_2}{r}, \\
V_{\ref{fig:Feynman:H-type:0PN:K}} &= -\frac{q_1 q_2}{4\pi r}.
\end{align}

Diagrams <ref> and <ref> belongs to the 1PN and represent corrections to the Green's function as described in equation (<ref>):

\begin{align}
V_{\ref{fig:Feynman:H-type:ox}} &= -\frac{Gm_1 m_2}{r}\left( \frac{\vb*{v}_1\vb*{\cdot} \vb*{v}_2 - \left( \vb*{v}_1 \vb*{\cdot} \vb*{\hat{r}} \right) \left( \vb*{v}_2 \vb*{\cdot} \vb*{\hat{r}} \right) }{2c^2} \right), \\
V_{\ref{fig:Feynman:H-type:ox:K}} &= -\frac{q_1 q_2}{4\pi r}\left( \frac{\vb*{v}_1\vb*{\cdot} \vb*{v}_2 - \left( \vb*{v}_1 \vb*{\cdot} \vb*{\hat{r}} \right) \left( \vb*{v}_2 \vb*{\cdot} \vb*{\hat{r}} \right) }{2c^2} \right).
\end{align}

The last two diagrams, <ref> and <ref>, represent general velocity expansion of the interaction term, mainly the expansion of Lorentz factors $\gamma_a^{-1} = 1 - \frac{v_a^2}{2c^2} - \frac{v^4_a}{8c^4} - \dots$, but for the graviton interaction one must also include spatial contributions of the source four-velocities, like $h\ind{_{ij}} \dot{x}\ind{^i} \dot{x}\ind{^j}$ and  $h\ind{_{0i}} \dot{x}\ind{^0} \dot{x}\ind{^i}$.

\begin{align}
V_{\ref{fig:Feynman:H-type:VelocityExp}} &= \begin{cases}
\textcolor{blue}{-\frac{Gm_1 m_2}{r} \frac{3\vb*{v}_1^2}{2c^2} } \\ \textcolor{red}{-\frac{Gm_1 m_2}{r} \frac{3\vb*{v}_2^2}{2c^2} } \\ \textcolor{teal}{\phantom{-}\frac{Gm_1 m_2}{r} \frac{4\vb*{v}_1 \vb*{\cdot} \vb*{v}_2}{c^2} }
\end{cases}, \\
V_{\ref{fig:Feynman:H-type:VelocityExp:K}} &= \begin{cases}
\textcolor{blue}{\phantom{-}\frac{q_1 q_2}{4\pi r} \frac{\vb*{v}_1^2}{2c^2} } \\ \textcolor{red}{\phantom{-}\frac{q_1 q_2}{4\pi r} \frac{\vb*{v}_2^2}{2c^2} }
\end{cases}.
\end{align}

For higher PN orders these expansions are continued and mixed, for example including a mix of diagram <ref> and <ref>, and a diagrams with two $\otimes$, and interaction points of the type $v^2 - v^2$ and $v^4 - v^0$ in analogy to diagrams <ref> and <ref>. But for the 1PN, figure <ref> contains all the contributing single propagator diagrams.

§.§ Higher order diagrams: seagull

Seagull-type diagrams contributing at 1PN. These diagrams originate from higher-order terms in the interaction expansion of the point particle action, see equation (<ref>).

Moving on to higher order diagrams, that is diagrams of higher order powers of the fields $h\ind{_{\mu\nu}}$ and $K$. Figure <ref> depicts the next-to-leading order interaction terms of equation (<ref>), namely $\text{fig. } \ref{fig:Feynman:Seagull:HH} \simeq m_a c^2 \times \frac{1}{8} \left( \lambda h\indices{_{\mu\nu}} \frac{\dot{x}\indices{^\mu}}{c} \frac{\dot{x}\indices{^\nu}}{c} \right)^2$, $\text{fig. } \ref{fig:Feynman:Seagull:KK} \simeq p_a K^2 \times \left(-1\right)$, and $\text{fig. } \ref{fig:Feynman:Seagull:HK} \simeq q_a K \times \frac{1}{2}\lambda h\indices{_{\mu\nu}} \frac{\dot{x}\indices{^\mu}}{c} \frac{\dot{x}\indices{^\nu}}{c}$. Inserting the leading order solution for $\lambda h\indices{_{00}}(x) = \sum_b \frac{\lambda^2 m_b \gamma_b c^2 \delta(x^0 - x^0_b)}{8\pi \abs{\vb*{x}-\vb*{x}_b }}$ and $K(x) = \sum_b -\frac{\gamma_b^{-1} q_b \delta(x^0 - x^0_b)}{4\pi \abs{\vb*{x} - \vb*{x}_b}}$, the resulting terms of the action originating from these diagrams and subsequently their potential is easily obtained.

\begin{align}
V_{\ref{fig:Feynman:Seagull:HH}} &= -\frac{ G^2 (m_1 m_2^2 + m_2 m_1^2) }{ 2 c^2 r^2 }, \\
V_{\ref{fig:Feynman:Seagull:KK}} &= \phantom{-} \frac{ (p_1 q_2^2 + p_2 q_1^2) }{ 16 \pi^2 r^2 }, \\
V_{\ref{fig:Feynman:Seagull:HK}} &= \phantom{-} \frac{ Gq_1 q_2 (m_1+m_2) }{ 4\pi c^2 r^2}.
\end{align}

§.§ Higher order diagrams: field-interactions

The diagrams with field-interactions contributing at 1PN. Wavy lines denote graviton propagators and dashed lines denote propagators of the $K$-scalar Kalb-Ramond field.

The last set of diagrams contributing to the 1PN action originates from corrections to the solution of the fields themselves, see equation (<ref>)-(<ref>), and are depicted diagrammatically in figure <ref> as field-interactions. The explicit form of the action term is here

\begin{align}
S_{\text{fig.}\ref{fig:3Point:diagrams}} = \int \frac{\dd[4]{x}}{c} \sum_{a} \gamma_a^{-1} \left[ \frac{m_a \lambda}{2}h_{00}^{(2)}(x) (\gamma_a c)^2 - q_a K^{(2)}(x) \right] \dirac{3}{\vb*{x} - \vb*{x}_a(x^0)}, \label{eq:ThreePoint:Action:intro}
\end{align}

where the corrections in the fields $h_{00}^{(2)}(x)$ and $K^{(2)}(x)$ are the items of interest.

§.§.§ The corresponding potential of diagram <ref>

Diagram <ref> originates from the graviton energy-momentum (pseudo-)tensor $t\ind{^{\alpha\beta}}$ contribution to the equation of motion (<ref>) with the corresponding solution:

\begin{align}
\begin{split}
h^{(2)}_{00}(x) \ni -\lambda \int \dd[4]{y} \Delta(x^\sigma - y^\sigma) P\indices{_{00\alpha\beta}} \frac{1}{2} h\ind{_{\mu\nu}^{,\alpha}} h\ind{^{\mu\nu,\beta}}.
\end{split}
\end{align}

To solve this, one can write out the $h\ind{_{\mu\nu}}$ fields as their leading order solution in Fourier space

\begin{align}
h\ind{_{\mu\nu}^{,\alpha}}(y) = & -\lambda \iint \frac{\dd[4]{p}}{(2\pi)^4} \frac{ ip\ind{^\alpha} e^{i k\ind{_{\rho}} \left( y\ind{^{\rho}} - z\ind{^{\rho}} \right) }}{ -p\ind{_{\sigma}} p\ind{^{\sigma}} } P\ind{_{\mu\nu00}} T\ind{^{00}}(z) \dd[4]{z}, \\
& P\ind{_{\mu\nu00}} = P\ind{_{00\mu\nu}} = \frac{1}{2} \delta\ind{_{\mu\nu}}.
\end{align}

For algebraic convenience, one might also utilise the following relation valid for any field $\Phi$ and $\Theta$:

\begin{align}
\Phi\ind{_{,\{\alpha}} \Theta\ind{_{,\beta\}}} = \frac{1}{2}\left( \left(
\Phi \Theta \right)\ind{_{,\alpha\beta}} - \Phi\ind{_{,\alpha\beta}} \Theta - \Phi \Theta\ind{_{,\alpha\beta}} \right). \label{eq:derivative trick}
\end{align}

Using the notation $\int \frac{\dd[4]{k}}{(2\pi)^4} = \int_k$, and taking advantage of rewriting the derivatives in accordance with (<ref>) one obtains

\begin{align}
\begin{split}
h^{(2)}_{00}(x) \ni -\lambda^3 \int \int_k \int_{p_1, p_2} \frac{ e^{ik_\sigma(x^\sigma - y^\sigma) } }{ k\ind{_\sigma} k\ind{^\sigma} } \frac{ e^{ip_{1\sigma}(y^\sigma - z_1^\sigma) } }{ p\ind{_{1\sigma}} p^{\sigma}_1 } \frac{ e^{ip_{2\sigma}(y^\sigma - z_2^\sigma) } }{ p\ind{_{2\sigma}} p^{\sigma}_2 } \frac{1}{8} \delta\ind{_{\alpha\beta}}\times \\
\Biggl[ (p_1^\alpha+p_2^\alpha)(p_1^\beta+p_2^\beta)
- p_1^\alpha p_1^\beta - p_2^\alpha p_2^\beta \Biggr] T\ind{^{00}}(z_1) T\ind{^{00}}(z_1) \dd[4]{y}
\end{split}
\end{align}

Separating out all the factors of, and integrating over $y$ results in a conservation of momentum factor

\begin{align*}
\int e^{-iy\ind{^{\mu}}\left( k\ind{_\mu} - p_{1\mu} - p_{2\mu} \right) } \dd[4]{y} = (2\pi)^4 \delta\hspace{-2pt}\left(k\ind{_\mu} - p_{1\mu} - p_{2\mu} \right).
\end{align*}

Then performing the $k$ integral makes $k\ind{_\mu} = p_{1\mu} + p_{2\mu} \equiv q\ind{_\mu}$, and tidying up the expression gives

\begin{align}
\begin{split}
h^{(2)}_{00}(x) \ni -\lambda^3 \int_{p_1, p_2} \frac{ e^{iq_\sigma x^\sigma } }{ q\ind{_\sigma} q\ind{^\sigma} } \frac{ e^{-ip_{1\sigma} z_1^\sigma } }{ p\ind{_{1\sigma}} p^{\sigma}_1 } \frac{ e^{-ip_{2\sigma} z_2^\sigma } }{ p\ind{_{2\sigma}} p^{\sigma}_2 } \frac{1}{8} \delta\ind{_{\alpha\beta}} \Biggl[ q^\alpha q^\beta - p_1^\alpha p_1^\beta - p_2^\alpha p_2^\beta \Biggr] T\ind{^{00}}(z_1) T\ind{^{00}}(z_2).
\end{split}
\end{align}

The final step is to again let $\frac{1}{k_\mu k^\mu} \to \frac{1}{\vb*{k}^2}$, and likewise ignore factors of $p^0$ appearing in the numerator as these also will result in velocity corrections and thus belong in higher PN order terms, see equation (<ref>). Then we finally have

\begin{align}
h^{(2)}_{00}(x) \ni -\frac{\lambda^3}{8} \int_{p_1, p_2} &\frac{ e^{iq_\sigma x^\sigma } }{ \vb*{q}^2 } \frac{ e^{-ip_{1\sigma} z_1^\sigma } }{ \vb*{p}_1^2 } \frac{ e^{-ip_{2\sigma} z_2^\sigma } }{ \vb*{p}_2^2 } \Biggl[ \vb*{q}^2 - \vb*{p}_1^2 - \vb*{p}_2^2 \Biggr] T\ind{^{00}}(z_1) T\ind{^{00}}(z_2) \nonumber \\
= -\frac{\lambda^3}{8} \int_{p_1,p_2} &\Biggl[ \frac{ e^{ip_{1\sigma} (x^\sigma - z_1^\sigma) } }{ \vb*{p}_1^2 } \frac{ e^{ip_{2\sigma} (x^\sigma - z_2^\sigma) } }{ \vb*{p}_2^2 } - \frac{ e^{iq_{\sigma} (x^\sigma - z_1^\sigma) } }{ \vb*{q}^2 } \frac{ e^{ip_{2\sigma} (z_1^\sigma - z_2^\sigma) } }{ \vb*{p}_2^2 } \label{eq:htripp:intermediate}  \\
&\quad \quad - \frac{ e^{iq_{\sigma} (x^\sigma - z_2^\sigma) } }{ \vb*{q}^2 } \frac{ e^{ip_{2\sigma} (z_2^\sigma - z_1^\sigma) } }{ \vb*{p}_1^2 } \Biggr] T\ind{^{00}}(z_1) T\ind{^{00}}(z_2). \nonumber
\end{align}

In the first term $\vb*{q}^2$, the first denominator is cancelled, and remembering that $q\ind{^\mu} = p_1^\mu + p_2^\mu$ the remaining exponential factor can be split up and distributed into the remaining exponential factors. Thus the first term looks like

\begin{align}
\begin{split}
& -\frac{\lambda^3}{8} \int_{p_1, p_2} \frac{ e^{ip_{1\sigma} \left( x\ind{^\sigma} - z_1^\sigma \right) } }{ \vb*{p}_1^2 } \frac{ e^{ip_{2\sigma}\left( x\ind{^\sigma} - z_2^\sigma \right) } }{ \vb*{p}_2^2 }T\ind{^{00}}(z_1) T\ind{^{00}}(z_2) \\
&= -\frac{\lambda^3}{8} \Delta_\text{inst}^{(0)}(x\ind{^\mu}-z_1^\mu) \Delta_\text{inst}^{(0)}(x\ind{^\mu}-z_2^\mu) T\ind{^{00}}(z_1) T\ind{^{00}}(z_2).
\end{split}
\end{align}

The procedure is the same for the two last terms of (<ref>), with the addition of shifting one of the integration variables to $q^\mu$.

Returning to the action (<ref>) for diagram <ref> it is

\begin{align}
S_\text{fig.\ref{fig:Feynman:3p:HHH}} = \int \sum_{a\neq b,c} \frac{\lambda^4 m_a m_b m_c c^6}{16 (4\pi)^2} & \left( \frac{-1}{\abs{\vb*{x}_a - \vb*{z}_b} \cdot \abs{\vb*{x}_a - \vb*{z}_c}} + \frac{1}{\abs{\vb*{x}_a - \vb*{z}_b} \cdot \abs{\vb*{z}_b - \vb*{z}_c}} + \frac{1}{\abs{\vb*{x}_a - \vb*{z}_c} \cdot \abs{\vb*{z}_c - \vb*{z}_b}} \right) \dd{t} \nonumber \\
&= \int -\frac{G^2 m_1 m_2(m_1+m_2)}{c^2 r^2} \dd{t}, \\
V_{\ref{fig:Feynman:3p:HHH}} &= \frac{G^2 m_1 m_2(m_1 + m_2)}{c^2 r^2}. \label{eq:potential:3p:HHH}
\end{align}

The requirement $b\neq a \neq c$ of the sum is to make sure the gravitational potential acting on particle number $a$ is not sourced by itself. Discarding then all divergent terms[For more rigour impose particle coordinates before performing the momentum integrals. Then one obtains terms like $\int_{p_2} e^{ip_{2\sigma} ( 0^\sigma )}/\vb*{p}_2^2 = \int_{p_2} 1/\vb*{p}_2^2 $ which also diverge to infinity.] $\propto \frac{1}{\abs{\vb*{z}_i - \vb*{z}_i }}$ one is left with the potential (<ref>).

§.§.§ The corresponding potential of diagram <ref>

Diagram <ref> is quite similar to diagram <ref>, but rather accounts for the Kalb-Ramond scalar contribution to $h_{00}^{(2)}(x)$ from the equation of motion (<ref>):

\begin{align}
h^{(2)}_{00}(x) \ni & \lambda \int \dd[4]{y} \Delta(x^\sigma - y^\sigma) P\indices{_{00\alpha\beta}} \left\{ \frac{1}{2} \eta\indices{^{\alpha\beta}} K_{,\sigma}(y) K^{,\sigma}(y) + K^{,\alpha}(y) K^{,\beta}(y) \right\} \\
\begin{split}
& = \frac{i^2\lambda}{2} \int_{p_1, p_2} \frac{e^{iq_\mu x^\mu } }{-\vb*{q}^2 } \frac{e^{-ip_{1\mu} z_1^\mu } }{ -\vb*{p}_{1}^2 } \frac{e^{-ip_{2\mu} z_2^\mu } }{-\vb*{p}_{2}^2 } \Biggl[ \vb*{q}^2 - \vb*{p}_{1}^2 - \vb*{p}_2^2 \Biggr] J(z_1) J(z_2), % red 1 blir \frac{1}{2} hvis red 3 egentlig er 1
\end{split}
\end{align}

where again $q^\mu = p_1^\mu + p_2^\mu$. These integrals are identical to the ones found in (<ref>), and are solved the same way.

\begin{align}
h_{00}^{(2)} \ni & \frac{\lambda}{2} \Bigl[ \Delta_\text{inst}(x-z_1) \Delta_\text{inst}(x-z_2) - \Delta_\text{inst}(x-z_1) \Delta_\text{inst}(z_1-z_2) - \Delta_\text{inst}(x-z_2) \Delta_\text{inst}(z_2-z_1)\Bigr] J(z_1) J(z_2) \nonumber \\
& = \sum_{b, c} \frac{\lambda q_b q_c}{32 \pi^2}\left[ \frac{ 1 }{ \abs{ \vb*{x} - \vb*{z}_b} \cdot \abs{ \vb*{x} - \vb*{z}_c } } - \frac{ 1 }{ \abs{ \vb*{x} - \vb*{z}_b } \cdot \abs{ \vb*{z}_b - \vb*{z}_c } } - \frac{ 1 }{ \abs{ \vb*{x} - \vb*{z}_c } \cdot \abs{ \vb*{z}_c - \vb*{z}_b } } \right].
\end{align}

By substituting this result back into the action (<ref>), the associated potential can be read off.

\begin{align}
S_\text{fig.\ref{fig:Feynman:3p:KKH}} = \int \sum_{a\neq b,c} \frac{ \lambda^2 m_a q_b q_c c^2}{64\pi^2} \Bigl( & \frac{ 1 }{ \abs{ \vb*{x} - \vb*{z}_b} \cdot \abs{ \vb*{x} - \vb*{z}_c } } - \frac{ 1 }{ \abs{ \vb*{x} - \vb*{z}_b } \cdot \abs{ \vb*{z}_b - \vb*{z}_c } } - \frac{ 1 }{ \abs{ \vb*{x} - \vb*{z}_c } \cdot \abs{ \vb*{z}_c - \vb*{z}_b } }\Bigr) \dd{t} \nonumber \\
& = \int \frac{G(m_1q_2^2+m_2 q_1^2) }{4\pi r^2 c^2} \dd{t}, \\
& V_{\ref{fig:Feynman:3p:KKH}} = - \frac{G(m_1q_2^2+m_2 q_1^2) }{4\pi r^2 c^2}.
\end{align}

§.§.§ The corresponding potential of diagram <ref>

For diagram <ref> it is the next-to-leading-order terms in the equation of motion for the Kalb-Ramond scalar (<ref>) which is included.

\begin{align}
K^{(2)}(x) = -\lambda \int \dd[4]{y} \Delta(x^\sigma-y^\sigma) \left( h\indices{^{\mu\nu}} K\indices{_{,\mu\nu}} + h\indices{^{\mu\nu}_{,\mu}} K\indices{_{,\nu}}  + \frac{1}{2} \left(h\indices{_{,\mu}} K^{,\mu} + h\dalembertian K \right) \right). \label{eq:Lambda(2):correction:intro}
\end{align}

Since the derivative indices are either contracted with themselves or with the graviton field  $h^{\mu\nu}\partial_\mu \partial_\nu$ they are symmetric under $\mu \leftrightarrow \nu$, and one may therefore use (<ref>) also in this context. Using this and the other tricks used on diagrams <ref> and fig:Feynman:3p:KKH one can demonstrate

\begin{align}
\begin{split}
K^{(2)}(x) = & -\frac{i^2\lambda^2}{2} \int_{p,k}  \frac{e^{iq_\sigma x^\sigma }}{-\vb*{q}^2} \frac{e^{-ik_\sigma z^\sigma_b }}{\vb*{k}^2} \frac{e^{-ip_\sigma z^\sigma_c }}{-\vb*{p}^2} \left( \vb*{q}^2 + \vb*{p}^2 - \vb*{k}^2 \right) T\ind{_{00}}(z_b) J(z_c) \label{eq:Lambda(2):correction} \\
= & \sum_{b,c} \frac{\lambda^2 m_b q_c c^2}{32 \pi^2} \left( \frac{1}{\abs{ \vb*{x} - \vb*{z}_b } \cdot \abs{ \vb*{x} - \vb*{z}_c }} + \frac{1}{\abs{ \vb*{x} - \vb*{z}_c } \cdot \abs{ \vb*{z}_c - \vb*{z}_b }} - \frac{1}{\abs{ \vb*{x} - \vb*{z}_b } \cdot \abs{ \vb*{z}_b - \vb*{z}_c }}\right)
\end{split}
\end{align}

Again $q^\mu = k^\mu + p^\mu$.

\begin{align}
S_\text{fig.\ref{fig:Feynman:3p:HKK}} = \int \sum_{a\neq b,c} \frac{ G q_a m_b q_c}{2\pi c^2} &\left( \frac{1}{\abs{ \vb*{x} - \vb*{z}_b } \cdot \abs{ \vb*{x} - \vb*{z}_c }} + \frac{1}{\abs{ \vb*{x} - \vb*{z}_c } \cdot \abs{ \vb*{z}_c - \vb*{z}_b }} - \frac{1}{\abs{ \vb*{x} - \vb*{z}_b } \cdot \abs{ \vb*{z}_b - \vb*{z}_c }}\right) \dd{t} \nonumber \\
& V_{\ref{fig:Feynman:3p:HKK}} = -\frac{ G q_1 q_2 (m_1+m_2)}{2\pi r^2 c^2} \label{eq:3p:Potential:HKK}
\end{align}

Notice how $V_{\ref{fig:Feynman:3p:HKK}} = -2\cdot V_{\ref{fig:Feynman:Seagull:HK}}$ and will thus `destructively interfere' with one another, analogous to how $V_{\ref{fig:Feynman:3p:HHH}} = -2 \cdot V_{\ref{fig:Feynman:Seagull:HH}}$ in GR.

§ THE RESULTING BINARY LAGRANGIAN

The kinetic term of the binary Lagrangian is obtained from equation (<ref>) by the cross-term $\sum_a \gamma_a^{-1} m_a c^2 \times (-1)$ and expanding the Lorentz factor in powers of velocity. To 1PN this is

\begin{align}
L_{pp}^\text{kin} & = \sum_a -m_a c^2 + \frac{1}{2} m_a v^2_a + \frac{1}{8} m_a \frac{v^4_a}{c^2} + \order{\frac{v^6}{c^4}} = -Mc^2 + \frac{1}{2}\mu v^2 + \frac{1-3\eta}{8} \mu \frac{v^4}{c^2} + \order{\frac{v^6}{c^4}}.
\end{align}

The constant term can be ignored for all dynamical purposes. In the last equality effective one-body coordinates have been implemented, where $M$ is the total mass $M=m_1+m_2$, $\mu$ is the reduced mass $\mu=\frac{m_1 m_2}{m_1+m_2}$, $\eta$ is the symmetric mass ratio $\eta=\frac{\mu}{M}$, and $v$ is the relative velocity $\vb*{v} = \dot{\vb*{r}} = \vb*{v}_1 - \vb*{v}_2$ when $\vb*{r}=\vb*{r}_1 - \vb*{r}_2$.

The contributions from standard general relativity are derived from all the diagrams only containing graviton propagators

\begin{align}
L_{pp}^\text{GR} = & \frac{Gm_1m_2}{r} + \frac{Gm_1m_2}{r}\left( \frac{3}{2}\frac{v_1^2+v_2^2}{c^2} - \frac{4\vb*{v}_1 \vb*{\cdot} \vb*{v}_2}{c^2} + \frac{\vb*{v}_1 \vb*{\cdot} \vb*{v}_2 - (\vb*{v}_1 \vb*{\cdot} \vb*{\hat{r}})(\vb*{v}_2 \vb*{\cdot} \vb*{\hat{r}}) }{2c^2} - \frac{G(m_1+m_2)}{2 r c^2} \right) \nonumber \\
= & \frac{GM\mu}{r} + \frac{GM\mu}{r}\left( \frac{3}{2} \left(1-2\eta\right) \frac{ v^2 }{c^2} + \frac{4\eta v^2}{c^2} - \eta \frac{ v^2 - (\vb*{v \cdot \hat{r}})^2}{2c^2} - \frac{GM}{2rc^2} \right) +\order{\frac{v^3}{c^3}}.
\end{align}

Next the additional potentials from the inclusion of the Kalb-Ramond (effectively scalar) field

\begin{align}
L_{pp}^\text{KR} &= \frac{q_1 q_2}{4\pi r} + \frac{q_1 q_2}{4\pi r}\left( -\frac{v_1^2+v_2^2}{2c^2} + \frac{\vb*{v}_1 \vb*{\cdot} \vb*{v}_2 - (\vb*{v}_1 \vb*{\cdot} \vb*{\hat{r}})(\vb*{v}_2 \vb*{\cdot} \vb*{\hat{r}}) }{2c^2} -\frac{\left(\frac{p_1 q_2}{q_1} + \frac{p_2 q_1}{q_2}\right)}{4\pi r} + \frac{G(m_1+m_2)}{r c^2} + \frac{G \left( \frac{m_1 q_2}{q_1} + \frac{m_2 q_1}{q_2} \right) }{r c^2}\right) \nonumber \\
&= \frac{q_1 q_2}{4\pi r} +  \frac{q_1 q_2}{4\pi r} \left( - \frac{1-2\eta}{2}\frac{v^2}{c^2} - \eta \frac{v^2 - \left(\vb*{v \cdot \hat{r}}\right)^2}{2c^2} - \frac{p}{4 \pi r} + \frac{GM}{rc^2} + \frac{G\left( \frac{m_1 q_2}{q_1} + \frac{m_2 q_1}{q_2} \right)}{ rc^2 } \right). \label{eq:Lagragian:KR}
\end{align}

In the last line $p=\frac{p_1 q_2}{q_1} + \frac{p_2 q_1}{q_2}$. Notice that if $p_a \propto q_a$ it follows that $p\propto q_1 + q_2$, i.e. the total charge, making it analogous to how diagram <ref> is a factor of total mass correction to Newton's potential in GR ($V_\text{\ref{fig:Feynman:Seagull:HH}}=-\frac{GM\mu}{r} \cdot \frac{-GM}{2 r c^2}$).

For comparison here is the Lagrangian obtained by  for axions in <cit.>, rewritten to have consistent notation with this paper:

\begin{align}
L_{pp}^\phi = & \frac{ q_1 q_2 }{4\pi r} e^{-m_s r}  + \frac{ q_1 q_2 }{4\pi r} e^{-m_s r} \left( \frac{\vb*{v}_1 \vb*{\cdot} \vb*{v}_2 - (\vb*{v}_1 \vb*{\cdot} \vb*{\hat{r}})(\vb*{v}_2 \vb*{\cdot} \vb*{\hat{r}})(1+m_s r) }{2c^2} -\frac{\left(\frac{p_1 q_2}{q_1} + \frac{p_2 q_1}{q_2}\right)}{2\pi r} e^{-m_s r} - \frac{G(m_1+m_2)}{r c^2} \right) \nonumber \\
& - \frac{G (m_1 q_2^2 + m_2 q_1^2)}{16\pi r} m_s \left[-1 + e^{-2m_s r} + 2m_s r \text{Ei}(-2m_s r) \right] + \frac{G M q_1 q_2}{2\pi rc^2} m_s \mathcal{I}(m_s r).
\end{align}

§ CONCLUSION AND OUTLOOK

We have derived the leading order relativistic Lagrangian corrections for a Kalb-Ramond field coupled to gravity, Eqn.(<ref>). Our computation would apply in the relativistic regime, but not the quantum regime. Astrophysical compact objects may or may not have measurable Kalb-Ramond charge. In the stringy fuzzball proposal, black holes are composed of macroscopic strings <cit.> and these may well have non-zero Kalb-Ramond charge. Fuzzballs should potentially also have non-trivial dilation profile, as happens in other black-hole solutions in dilaton gravity and supergravity <cit.>. But unlike the Kalb-Ramond (and photon) field, the dilaton is not protected by a symmetry and so can dynamically acquire a mass, leading to a short range force even at the classical level. If macroscopic strings exist at cosmological scales, then they may also interact via a Kalb-Ramond mediated force <cit.>.

In our derived Lagrangian, note the similar $\tfrac{1}{r}$-behaviours in the relativistic expansions for both the Kalb-Ramond and classical General Relativity corrections. This is perhaps not too surprising in view of generalised geometry, where the metric and Kalb-Ramond fields are put on the same footing in the generalised metric <cit.>. This point of view also leads us to speculate that this similar behaviour will persist to higher orders. In order to investigate this further, it would be very useful to extend the perturbative Feynman diagram formalism of <cit.> to the generalised framework. This framework will also be relevant when we come to understand the theory in higher dimensions or with quantum corrections, where the Kalb-Ramond field cannot be treated as a scalar anymore.

The approach presented here also demonstrates the power of the effective field theory approach to classifying and calculating possible deviations from Einstein's theory of gravity. This is relevant in an era of precision pulsar and gravitational wave observations that are testing Einstein's theory in the strongly relativistic regime.

§ ACKNOWLEDGEMENTS

We thank Jahed Abedi and Edward Hardy for useful conversations.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

journaltitlePhys. Rev. Lett.
titleCP Conservation in the Presence of Instantons

journaltitlePhys. Rev. Lett.
titleProblem of Strong $P$ and $T$ Invariance in the Presence of Instantons

journaltitlePhys. Rev. Lett.
titleA New Light Boson?

journaltitlePhys. Rev. Lett.
titleWeak Interaction Singlet and Strong CP Invariance

journaltitlePhys. Rev. Lett.
titleParity Conservation in QCD

titleThe QCD axion, precisely

titleAxions In String Theory

journaltitlePhys. Rev. D
titleString Axiverse

journaltitlePhys. Rept.
titleAxion Cosmology

American Physical Society (APS)

journaltitlePhysical Review D
titleProspects for axion searches with Advanced LIGO through binary mergers

journaltitlePhys. Rev. D
titleClassical direct interstring action

journaltitleNucl. Phys. B
titleOn spherically symmetric string solutions in four-dimensions

journaltitleEur. Phys. J. C
titleCausal structures and dynamics of black-hole-like solutions in string theory

journaltitleAstrophys. J. Lett.
titleThe NANOGrav 15 yr Data Set: Search for Signals from New Physics

titleCosmic F and D strings

journaltitleFortsch. Phys.
titleThe Fuzzball proposal for black holes: An Elementary review

journaltitlePhys. Rev. X
titleStrong-Field Gravity Tests with the Double Pulsar

journaltitlePhys. Rev. D
titleTests of General Relativity with the Binary Black Hole Signals from the LIGO-Virgo Catalog GWTC-1

journaltitleLiving Rev. Rel.
titleProspects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA

titleScience Case for the Einstein Telescope

journaltitleGen. Rel. Grav.
titleProspects for Fundamental Physics with LISA

American Physical Society (APS)

journaltitlePhysical Review D
titleEffective field theory of gravity for extended objects

titleLes Houches Lectures on Effective Field Theories and Gravitational Radiation

journaltitlearXiv preprint arXiv:1008.0973
titleLectures on generalized geometry

titleSupergravity as Generalised Geometry I: Type II Theories

titleSupergravity as Generalised Geometry II: $E_{d(d)} \times \mathbb{R}^+$ and M theory

journaltitleProQuest Dissertations and Thesis
titleModeling Axion-Graviton Interactions from an Effective Action Approach

Berlin Heidelberg

seriesGraduate Texts in Physics
titleSpecial Relativity in General Frames

titleFeynman lectures on gravitation

Physical Review D

journaltitlePhysical Review D
titleFactorization and polarization in linearized gravity

Elsevier BV

journaltitlePhysics Reports
titleThe effective field theorist’s approach to gravitational dynamics

journaltitlePhys. Rev. D
titleFirewall black holes and echoes from an action principle


# The first quenched galaxies, when and how?

Lizhi Xie (谢利智<EMAIL_ADDRESS>Tianjin Normal University,
Binshuixidao 393, Xiqing, 300387, Tianjin, China
INAF – Astronomical Observatory of Trieste, via G.B. Tiepolo 11, I-34143
Trieste, Italy
Gabriella De Lucia INAF – Astronomical Observatory of Trieste, via G.B.
Tiepolo 11, I-34143 Trieste, Italy
IFPU - Institute for Fundamental Physics of the Universe, via Beirut 2, 34151,
Trieste, Italy
Fabio Fontanot INAF – Astronomical Observatory of Trieste, via G.B. Tiepolo
11, I-34143 Trieste, Italy
IFPU - Institute for Fundamental Physics of the Universe, via Beirut 2, 34151,
Trieste, Italy
Michaela Hirschmann Institute for Physics, Laboratory for Galaxy Evolution and
Spectral Modelling,
Ecole Polytechnique Federale de Lausanne, Observatoire de Sauverny, Chemin
Pegasi 51, CH-1290 Versoix, Switzerland  INAF – Astronomical Observatory of
Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy
Yannick M. Bahé Laboratory for Astrophysics, École Polytechnique Fédérale de
Lausanne (EPFL), Observatoire de Sauverny, Chemin Pegasi 51,
CH-1290 Versoix, Switzerland  Michael L. Balogh Department of Physics and
Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada  Adam
Muzzin Department of Physics and Astronomy, York University, 4700 Keele St.,
Toronto, Ontario MJ3 1P3, Canada  Benedetta Vulcani INAF- Osservatorio
Astronomico di Padova, Vicolo Osservatorio 5, 35122 Padova, Italy  Devontae C.
Baxter Department of Astronomy & Astrophysics, University of California, San
Diego, 9500 Gilman Dr, La Jolla, CA 92093, USA
Ben Forrest Department of Physics and Astronomy, University of California
Davis, One Shields Avenue, Davis, CA 95616 USA  Gillian Wilson Department of
Physics, University of California Merced, 5200 Lake Road, Merced, CA 95353,
USA  Gregory H. Rudnick Department of Physics and Astronomy, University of
Kansas, Lawrence, KS 66045, USA  M. C. Cooper Department of Physics and
Astronomy, University of California, Irvine, 4129 Frederick Reines Hall,
Irvine, CA 92697, USA  Umberto Rescigno Instituto de Astronomía y Ciencias
Planetarias de Atacama (INCT), Universidad de Atacama, Copayapu 485, Copiapó,
Chile
Instituto de Astrofisica, Universidad Andres Bello, Fernandez Concha 700, Las
Condes, Santiago RM, Chile

###### Abstract

Many quenched galaxies discovered in the early Universe by JWST raise
fundemental question s on when and how these galaxies became quiescent. Making
use of the latest version of the semi-analytic model GAEA that provides good
agreement with the observed quenched fractions up to $z\sim 3$, we make
predictions for the expected fractions of quiescent galaxies up to $z\sim 7$
and analyze the main quenching mechanism. We find that in a simulated box of
$685~{}{\rm Mpc}$ on a side, the first quenched massive ($M_{\star}\sim
10^{11}{\rm M}_{\odot}$), Milky Way mass, and low mass ($M_{\star}\sim
10^{9.5}{\rm M}_{\odot}$ ) galaxies appear at $z\sim 4.5$, $z\sim 6.2$, and
before $z=7$. Most quenched galaxies identified at early redshifts remain
quenched for more than 1 Gyr. Independently of galaxy stellar mass, the
dominant quenching mechanism at high redshift is accretion disk feedback
(quasar winds) from a central massive black hole, which is triggered by
mergers in massive and Milky Way mass galaxies, and by disk instabilities in
low-mass galaxies. Environmental stripping become increasingly more important
at lower redshift.

Galaxy — Quenching — Simulation

##  1 Introduction

The cessation of star formation in galaxies has drawn considerable attention
in recent years, especially given the large numbers of quiescent massive
galaxies that have been found in the early Universe (Schreiber et al., 2018;
Merlin et al., 2019; Girelli et al., 2019; Glazebrook et al., 2017;
Nanayakkara et al., 2022), where the timescale available to assemble and
quench these systems is short. Spectroscopic confirmation suggests that these
galaxies experience short periods of intense star formation, grow up to
$10^{11}{\rm M}_{\odot}$ in the first one or two billion years of the
universe, and then stop forming stars within the next few million years
(Forrest et al., 2020; Valentino et al., 2020; Kakimoto et al., 2023; Carnall
et al., 2023a). This rapid assembly and quenching process might challenge our
current understanding of galaxy formation (Finkelstein et al., 2023).

The number densities of quenched massive galaxies $M_{\star}>10^{10.5}{\rm
M}_{\odot}$ increase rapidly from $\sim 10^{-6}~{}{\rm Mpc}^{-3}$ at $z\sim 5$
to as much as a factor of $10$ times higher at $z\sim 3$ (Marsan et al.,
2022), although the measured number densities have a relatively large scatter
due to different selection criteria and cosmic variance (Valentino et al.,
2023). The classical UVJ color selection of quenched galaxies (Wuyts et al.,
2008) is found to be incomplete and underestimates the number of quenched
galaxies at $z>3$ (Schreiber et al., 2018). Some studies favour a NUVrJ
(Ilbert et al., 2013) or ugi color selection (Antwi-Danso et al., 2023) to
identify galaxies that have been quenched recently, which is important for
galaxies at $z>3$ (Gould et al., 2023; Kubo et al., 2024).

Despite large scatter in measurements, it is a solid conclusion that most
theoretical models under-predict the number densities of quenched galaxies
(Cecchi et al., 2019; Girelli et al., 2019; Gould et al., 2023) at $z>4$ by
about an order of magnitude. Weaver et al. (2023) found that quenched galaxies
in the SHARK model (Lagos et al., 2023) and IllustrisTNG simulation (Pillepich
et al., 2018) at $3.5<z<4.5$ are not as massive as the observed ones. Either
creating enough massive galaxies at early cosmic epochs or quenching them on a
short time scale remains a challenge for current galaxy formation models.

Various physical mechanisms have been proposed to explain the rapid assembly
of mass in the early Universe, including weaker feedback (Dekel et al., 2023),
enhanced star formation efficiencies (Wang et al., 2023), and a top-heavy IMF
(Trinca et al., 2023). The physical mechanisms driving quenching also remain
unclear. Quenching could be caused by internal feedback from active galactic
nuclei (AGN) and supernovae (SN) feedback, or by external physical processes
including galaxy-galaxy interactions and environmental stripping. With a
minimal physical model, Gelli et al. (2023) suggests that SN feedback is not
powerful enough to quench galaxies of $\sim 10^{8}{\rm M}_{\odot}$ at high
redshift. The fact that many high-$z$ quenched galaxies are found to host
luminous AGN (Ito et al., 2022; Shimakawa et al., 2023) suggests an important
contribution to quenching from feedback from their central supermassive black
holes (SMBH). This appears to be confirmed in recent theoretical works:
Kurinchi-Vendhan et al. (2023) and Kimmig et al. (2023) analyze the massive
quenched galaxies in TNG and Magneticum at $z\sim 3$ and show that these
galaxies are indeed quenched by AGN feedback. Lovell et al. (2023) found the
AGN feedback is the dominant quenching mechanism for galaxies above
$10^{9}{\rm M}_{\odot}$ at $z\sim 5$. Qin et al. (2017) use semi-analytic
models to identify analogs of quenched galaxies observed at $z\sim 5$ and show
that these have grown through significant mergers and host the most massive
black holes at their redshifts. Some quenched galaxies are found in over-dense
environments (Kubo et al., 2021; McConachie et al., 2022; Tanaka et al., 2023;
Alberts et al., 2023), suggesting environmental quenching may also contribute
as early as $z\sim 4$.

In our recent work (De Lucia et al., 2024, here after GAEA2023), we present
the latest version of our GAlaxy Evolution and Assembly (GAEA) model and show
that it can correctly reproduce the observed quenched fractions up to redshift
$\sim 3$ as well as the number densities of quenched galaxies up to redshift
$\sim 5$, better than many state-of-the-art models and simulations. The good
agreement with observations makes it a perfect tool for studying the physical
origin of quenched galaxies at high redshift.

In this work, we extend the analysis presented in De Lucia et al. (2024) to
the quenched fractions and their quenching mechanisms since $z\sim 7$. In
Section 2, we introduce the semi-analytic model and our methodology. In
Section 3 and Section 4, we present the results and give our conclusions.

##  2 Model and Simulation

GAEA2023 (De Lucia et al., 2014) now combines independent versions of the
model including an improved treatment for the supernovae feedback (Hirschmann
et al., 2016), of the multi-phase cold gas (Xie et al., 2017), of
environmental effects (Xie et al., 2020), and of AGN accretion and feedback
(Fontanot et al., 2020). In particular, GAEA2023 implements a treatment for
tidal stripping and ram pressure stripping that gradually removes hot gas, as
well as ram-pressure stripping of cold gas, for satellite galaxies. These
implementations help to solve the over-quenching of low-mass satellite
galaxies and to improve the model predictions in terms of gas fractions (Xie
et al., 2020). GAEA2023 also implements an updated modeling for cold gas
accretion onto supermassive black holes. Mergers and disk instability cause a
fraction of cold gas to lose angular momentum and flow towards the center,
where it forms an accretion disk around the super massive black hole. This
material is then accreted onto the SMBH on a viscous timescale: accretion can
heat the surrounding gas and cause an outflow (for details about the model, we
refer to Fontanot et al. 2020. Below, we refer to this process as accretion
disk feedback. GAEA2023 has been calibrated to reproduce the galaxy stellar
mass function up to $z\sim 3$, HI mass function and quenched fractions at
$z\sim 0$, as well as AGN luminosity function up to $z\sim 4$. In De Lucia et
al. (2024) we show that this model version reproduces well the observed
quenched fraction, as well as the stellar mass function of the quenched
population up to $z\sim 3$. GAEA2023 is in good agreement with the observed
quenched fraction, as well as the stellar mass function of the quenched
population up to $z\sim 4$.

The model is run on the Millennium Simulation (Springel et al., 2005) with a
box size of $685$ Mpc based on a WMAP 1-yr cosmology (Spergel et al., 2003)
with $\Omega_{m}=0.25$, $\sigma_{b}=0.045$, $\sigma_{8}=0.9$, and $h=0.73$.

In the following, we will compare GAEA2023 results with predictions from
TNG100 and TNG300 (Springel et al., 2018; Nelson et al., 2018; Naiman et al.,
2018; Marinacci et al., 2018; Pillepich et al., 2018). The TNG project is a
suite of cosmological magneto-hydro-dynamical simulations, adopting the Planck
cosmology (Planck Collaboration et al., 2016) with $\Omega_{m}=0.3089$,
$\Omega_{b}=0.0486$, $\sigma_{8}=0.8159$, and $h=0.6774$. The TNG100 and
TNG300 have model cubic boxes of side length approximately 100 and 300 Mpc.
TNG considers two modes of AGN feedback: for high accretion rates, the
surrounding gas is heated by thermal feedback from AGN. When the accretion
rates are low, gas instead receive a kinetic ‘kick’ that causes gas outflows.
In this work, we use the publicly available database to retrieve the simulated
stellar mass, star formation rate, and black hole mass within twice the
stellar half-mass radius.

##  3 Results

###  3.1 Quenched fraction

Figure 1 shows the evolution of quenched fractions as predicted by GAEA2023
for low-mass ($2-4\times 10^{9}{\rm M}_{\odot}$), Milky-Way-mass ($2-4\times
10^{10}{\rm M}_{\odot}$), and massive galaxies ($0.8-1.5\times 10^{11}{\rm
M}_{\odot}$). We consider different defintions for quenched galaxies: first,
we select a sample imposing specific star formation rate
$sSFR=SFR/M_{\star}<0.3/t_{H}$, where $t_{H}$ is the Hubble time at given
redshift (Franx et al., 2008). Then in order to get a fair comparison with
observational data, we also present quenched fractions based on a UVJ
(synthetic) color selection. We tried different cuts commonly adopted in the
literature (Williams et al., 2009; Whitaker et al., 2011; Muzzin et al., 2013;
Martis et al., 2016) and plot the scatter obtained as shaded regions in Fig.
1. The error boxes show observational measurements of quenched fractions for
galaxies from the UltraVISTA DR1 and 3D-HST surveys (Martis et al., 2016). In
the redshift range $0<z<3$, there is a good agreement between GAEA2023 and
data at all mass scales. The quenched fraction defined using the synthetic UVJ
photometry are similar to those defined by sSFR at $z<2$, but decrease more
rapidly at $z>2$. As for GAEA2023, the UVJ color selection fails to find
significant numbers of quiescent galaxies at $z>2$. For TNG the quenched
fractions of galaxies with $M_{\star}\sim 10^{11}{\rm M}_{\odot}$ are larger
than observational measurement.

Moving to higher redshifts, the predicted quenched fractions decrease. In the
framework of GAEA2023, the quenched fraction of low-mass galaxies is 0.2 per
cent at $z\sim 7.3$, i.e. 6 out of 2880 galaxies are quenched. The quenched
fraction remains below 1 per cent until $z\sim 3$. We traced 42 quenched low-
mass galaxies at $z\sim 6.2$ forward in time and found that 11 of them move
back to the main sequence within 0.5 Gyr. Most of the high-redshift quenched
low-mass galaxies, however, remain quiescent for more than 1 Gyr, indicating
they are not just temporarily quenched between two episodes of bursty star
formation.

The first quenched galaxies with mass similar to the Milky Way appear at
$z\sim 6.2$, i.e. 2 out of 116 Milky Way-like galaxies are quenched. One of
these two returns to the main sequence after $\sim 0.5$ Gyr. The other is a
satellite that remains quenched until it merges with another galaxy. The
quenched fraction, for galaxies of this mass, grows quickly to $10$ per cent
between $4<z<5$.

The first massive quenched galaxies (with $M_{\star}\sim 10^{11}{\rm
M}_{\odot}$) are found at $z\sim 4.5$. 2 out of 35 galaxies in this mass bin
are quenched at this redshift. In the subsequent evolution, these two galaxy
continue to have a low star formation rate for most of the time. We have
tested low cut $sSFR<10^{-11}~{}{\rm yr}^{-1}$ to define quiescent galaxies
and find our results do not depend significantly on the SFR cut adopted.

Predictions from TNG are quite different from those based on GAEA2023, with
systematically lower quenched fractions at high redshift, and no quenched
galaxy at $z>4$. It is in clear tension with the observed fact of
spectroscopically confirmed quenched massive galaxies at $z>4$(e.g. Carnall et
al., 2023b). A similar result is reported in Merlin et al. (2019), where a
lower sSFR cut $10^{-11}~{}{\rm yr}^{-1}$ was used. Though most massive
galaxies at $z\sim 4$ have a central massive black hole of $10^{8}{\rm
M}_{\odot}$, only a minor fraction of them are quenched. In TNG, kinetic
feedback from AGN represents the most efficient quenching mechanism for
massive galaxies (Kurinchi-Vendhan et al., 2023). However, the accretion rates
at high redshift are so large that the assumed mode for AGN feedback in TNG
comes in the form of thermal feedback, which is not efficient enough to quench
galaxies at $z>3$.

Figure 1: Quenched fractions as a function of redshift. Different colors
represent galaxies of different stellar mass. Vertical error bars show the
standard deviations obtained for 100 randomly selected sub-samples. Shaded
regions show the uncertainties in quiescent fractions from slightly different
cuts in the UVJ diagram (more details in text). The dash-dotted line shows
results from a GAEA run where accretion disk feedback is switched off. Dotted
and dashed lines show the quenched fractions measured from TNG300 and TNG100.
Error boxes are observed quenched fractions for UltraVISTA DR1 and 3D-HST
survey from Martis et al. (2016).

###  3.2 Quenching Mechanisms

Figure 2: Evolution of three representative galaxies. The upper and lower
panels show the evolution of sSFR and the SMBH accretion rate associated with
the disk accretion mode normalized by stellar mass. The sSFR is indicated by
solid and dotted lines when a galaxy is classified as central or satellite,
respectively. The gray straight line in the upper panel is the separation
between quenched and star-forming galaxies. Large and small triangles in the
lower panel indicate merger events with mass ratios above 0.3 and 0.01. Color
code is the same as in Figure 1.

Broadly speaking, we can consider two quenching scenarios: internal quenching,
i.e. AGN feedback and SN feedback, and external quenching, i.e. environmental
stripping and galaxy-galaxy interactions. In this section, we analyze the
relative importance of these quenching mechanisms at different redshifts in
the GAEA framework.

First of all, we traced the history of high-redshift quenched galaxies in the
three stellar mass ranges considered. Figure 2 shows the evolution histories
of three representative galaxies selected from each stellar mass bin. We find
that all quenched model galaxies have experienced high-rate black hole
accretion and suffered accretion disk feedback right before quenching,
suggesting that this mechanism is the main quenching channel at this redshift.
This is confirmed by the dot-dashed lines in Figure 1, showing predictions for
Milky Way mass galaxies from the a version of GAEA (Xie et al., 2020) that
does not include disk accretion feedback and that significantly under-predicts
the fractions of quenched galaxies at high redshift (see also in De Lucia et
al., 2024).

Since the large accretion rates give rise to luminous quasars, in the top
panel of Figure 3 we compare the fraction of quasar-host galaxies in recently-
quenched galaxies (dashed) and the entire population (solid). About 70 and 40
per cent of massive and Milky Way mass galaxies host AGN with bolometric
luminosity brighter than $10^{44}$ erg/s at $z>2$. These fractions rise to 100
per cent for recently quenched galaxies at $z>2$. Surprisingly, more than
$75\%$ of low-mass quenched galaxies at $z>2$ also host luminous AGN, whereas
the fraction is only 20 per cent for all galaxies in this mass range. The
elevated fraction of luminous AGN for quenched galaxies confirms that the disk
accretion feedback from SMBH is the dominant quenching mechanism for high
redshift galaxies. The fraction of luminous AGN decreases at lower redshift,
and the differences between all galaxies and quenched samples also reduce: at
low redshift, accretion disk feedback is less important for quenching.

The disk accretion feedback is triggered by both disk instabilities and
mergers. While tracing evolution of individual galaxies, we find the large
accretion rates are associatd with merger events for Milky Way mass and
massive galaxies. In the middle panel of Figure 3 we show the fraction of
galaxies that just experienced mergers between the recently quenched and in
the entire population. All mergers with a mass ratio larger than $1:100$ are
considered, motivated by the rapid growth of black holes driven by multiple
minor mergers or even very minor mergers, which we find to be common for
high-z quiescent galaxies. Compared to the entire sample of Milky Way mass and
massive galaxies, newly quenched galaxies have much higher merger rates,
especially at high redshift. Therefore, mergers represent the main channel for
black hole accretion at these mass ranges in the GAEA framework.

For low-mass galaxies, more than half of the recently quenched galaxies
haven’t experienced any mergers around the time of quenching. We thus conclude
that their SMBH accretion are not primarily driven by mergers. We find that
low-mass galaxies are more likely to have unstable disks at high redshift, and
it is this disk instability that triggers efficient black hole accretion.

The connection between the quenching process and black hole accretion becomes
weaker at lower redshift, where environmental effects become increasingly
important. The bottom panel of Figure 3 shows the fraction of galaxies that
are satellites. In GAEA2023, satellite galaxies lose hot gas and cold gas
gradually by tidal stripping and ram-pressure stripping, whereas central
galaxies are unaffected. For low-mass galaxies, a larger fraction of satellite
galaxies are quenched at as early as $z\sim 6$, so the dependence of quenching
on the environment starts at very early epochs in the framework of our model.
The phenomenon of environmental quenching starts since $z\sim 3$ for Milky-
Way-mass galaxies. The difference in satellite fraction between quenched and
all massive galaxies is negligible, which is consistent with previous results
that massive galaxies are mainly quenched by AGN feedback (Xie et al., 2020).

Figure 3: The top, middle, and bottom panels show the fraction of galaxies
hosting luminous AGN ($\log L_{bol}/({\rm erg/s})>44$), the fraction of
galaxies that have experienced recent mergers with a mass ratio larger than
1/100, and the fraction of satellites, respectively. Solid, dotted, and dashed
lines correspond to the total, quenched, and newly-quenched galaxy
populations, respectively. Color code is the same as in Figure 1.

##  4 Conclusion

In this work we use the semi-analytic model GAEA2023 to study the quenched
fractions predicted for massive ($M_{\star}\sim 10^{11}{\rm M}_{\odot}$),
Milky Way mass ($M_{\star}\sim 10^{10.5}{\rm M}_{\odot}$), and low mass
galaxies ($M_{\star}\sim 10^{9.5}{\rm M}_{\odot}$) since $z=7$. GAEA2023
predictions are in good agreement with the observed quenched fraction measured
from UltraVista and 3D-HST in the redshift range $0<z<3$.

The quenched fractions defined by UVJ color are consistent with those defined
by sSFR up to $z\sim 3$. At higher redshift, the quenched fractions are under-
estimated by a UVJ color selection. When adopting a sSFR selection, about
$5\%$ of massive galaxies are found to be quenched at $z\sim 4.5$, about $2\%$
of Milky Way mass galaxies are firstly found to be quenched at $z\sim 6.2$,
and the quenched fraction of low-mass galaxies is $0.2\%$ at $z\sim 7$. More
than half of galaxies are not temporarily quenched due to bursty star
formation episodes, as they maintain a low star formation rate for over $\sim
1$ Gyr.

Taking advantage of our model, we analyse the quenching mechanism at different
redshifts. We find that all recently quenched Milky Way and massive galaxies,
and more than 75 per cent of low mass newly quenched galaxies at $z>2$ host
luminous quasars (with bolometric luminosity brighter than $10^{44}$ erg/s).
This suggests that accretion disk feedback from SMBHs is the main reason for
quenching at high redshift. This is also confirmed by analysing predictions
from an alternative model where this physical process is switched off. We find
that disk accretion feedback responsible for quenching is driven by galaxy
mergers for massive and Milky Way mass galaxies, and by disk instabilities for
lower mass galaxies. Environmental effects become increasingly important for
low-mass galaxies at $z<6$, and for Milky Way mass galaxies at $z<2$. Massive
galaxies are not quenched by environmental effects(Hirschmann et al., 2016; De
Lucia et al., 2019).

The earliest quenched galaxy so far is at $z\sim 7$ (Looser et al., 2023), and
it has a stellar mass of $M_{\star}\sim 5\times 10^{8}{\rm M}_{\odot}$. At
$z\sim 5$, most quenched galaxies discovered are massive galaxies. Based on
predictions from our model, we expect to find non-negligible numbers of
quenched galaxies with stellar mass $\sim 10^{9.5}{\rm M}_{\odot}$ at $z\sim
7$ or even higher redshift. A large fraction of these galaxies are expected to
host liminous quasars.

##  5 Acknowledgments

LZX acknowledges support from the National Natural Science Foundation of China
(No. 12041302, No. 11903023). We also acknowledge the support of Lorentz
Center, ISSI and the use of their facilities, which facilitated discussions on
this project. BV acknowledges the support from the INAF Mini Grant 2022
“Tracing filaments through cosmic time” (PI Vulcani). MLB acknowledges support
from an NSERC Discovery Grant. GHR acknowledges support from NSF Astronomy and
Astrophysics grant 2206473, HST grant GO-16300.004-A, and NASA ADAP grant
80NSSC19K0592. GW gratefully acknowledges support from the National Science
Foundation through grant AST-2205189 and from HST program number GO-16300. YMB
gratefully acknowledges funding from the Netherlands Organization for
Scientific Research (NWO) under Veni grant number 639.041.751 and financial
support from the Swiss National Science Foundation (SNSF) under funding
reference 200021_213076.

## References

  * Alberts et al. (2023) Alberts, S., Williams, C. C., Helton, J. M., et al. 2023, arXiv e-prints, arXiv:2312.12207, doi: 10.48550/arXiv.2312.12207
  * Antwi-Danso et al. (2023) Antwi-Danso, J., Papovich, C., Esdaile, J., et al. 2023, arXiv e-prints, arXiv:2307.09590, doi: 10.48550/arXiv.2307.09590
  * Carnall et al. (2023a) Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2023a, Nature, 619, 716, doi: 10.1038/s41586-023-06158-6
  * Carnall et al. (2023b) Carnall, A. C., McLeod, D. J., McLure, R. J., et al. 2023b, MNRAS, 520, 3974, doi: 10.1093/mnras/stad369
  * Cecchi et al. (2019) Cecchi, R., Bolzonella, M., Cimatti, A., & Girelli, G. 2019, ApJ, 880, L14, doi: 10.3847/2041-8213/ab2c80
  * De Lucia et al. (2024) De Lucia, G., Fontanot, F., Xie, L., & Hirschmann, M. 2024, arXiv e-prints, arXiv:2401.06211, doi: 10.48550/arXiv.2401.06211
  * De Lucia et al. (2019) De Lucia, G., Hirschmann, M., & Fontanot, F. 2019, MNRAS, 482, 5041, doi: 10.1093/mnras/sty3059
  * De Lucia et al. (2014) De Lucia, G., Tornatore, L., Frenk, C. S., et al. 2014, MNRAS, 445, 970, doi: 10.1093/mnras/stu1752
  * Dekel et al. (2023) Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201, doi: 10.1093/mnras/stad1557
  * Finkelstein et al. (2023) Finkelstein, S. L., Leung, G. C. K., Bagley, M. B., et al. 2023, arXiv e-prints, arXiv:2311.04279, doi: 10.48550/arXiv.2311.04279
  * Fontanot et al. (2020) Fontanot, F., De Lucia, G., Hirschmann, M., et al. 2020, MNRAS, 496, 3943, doi: 10.1093/mnras/staa1716
  * Forrest et al. (2020) Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020, ApJ, 903, 47, doi: 10.3847/1538-4357/abb819
  * Franx et al. (2008) Franx, M., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2008, ApJ, 688, 770, doi: 10.1086/592431
  * Gelli et al. (2023) Gelli, V., Salvadori, S., Ferrara, A., & Pallottini, A. 2023, arXiv e-prints, arXiv:2310.03065, doi: 10.48550/arXiv.2310.03065
  * Girelli et al. (2019) Girelli, G., Bolzonella, M., & Cimatti, A. 2019, A&A, 632, A80, doi: 10.1051/0004-6361/201834547
  * Glazebrook et al. (2017) Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Nature, 544, 71, doi: 10.1038/nature21680
  * Gould et al. (2023) Gould, K. M. L., Brammer, G., Valentino, F., et al. 2023, AJ, 165, 248, doi: 10.3847/1538-3881/accadc
  * Hirschmann et al. (2016) Hirschmann, M., De Lucia, G., & Fontanot, F. 2016, MNRAS, 461, 1760, doi: 10.1093/mnras/stw1318
  * Ilbert et al. (2013) Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55, doi: 10.1051/0004-6361/201321100
  * Ito et al. (2022) Ito, K., Tanaka, M., Miyaji, T., et al. 2022, ApJ, 929, 53, doi: 10.3847/1538-4357/ac5aaf
  * Kakimoto et al. (2023) Kakimoto, T., Tanaka, M., Onodera, M., et al. 2023, arXiv e-prints, arXiv:2308.15011, doi: 10.48550/arXiv.2308.15011
  * Kimmig et al. (2023) Kimmig, L. C., Remus, R.-S., Seidel, B., et al. 2023, arXiv e-prints, arXiv:2310.16085, doi: 10.48550/arXiv.2310.16085
  * Kubo et al. (2024) Kubo, M., Nagao, T., Uchiyama, H., et al. 2024, MNRAS, 527, 403, doi: 10.1093/mnras/stad3210
  * Kubo et al. (2021) Kubo, M., Umehata, H., Matsuda, Y., et al. 2021, ApJ, 919, 6, doi: 10.3847/1538-4357/ac0cf8
  * Kurinchi-Vendhan et al. (2023) Kurinchi-Vendhan, S., Farcy, M., Hirschmann, M., & Valentino, F. 2023, arXiv e-prints, arXiv:2310.03083, doi: 10.48550/arXiv.2310.03083
  * Lagos et al. (2023) Lagos, C. D. P., Bravo, M., Tobar, R., et al. 2023, arXiv e-prints, arXiv:2309.02310, doi: 10.48550/arXiv.2309.02310
  * Looser et al. (2023) Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2023, arXiv e-prints, arXiv:2302.14155, doi: 10.48550/arXiv.2302.14155
  * Lovell et al. (2023) Lovell, C. C., Roper, W., Vijayan, A. P., et al. 2023, MNRAS, 525, 5520, doi: 10.1093/mnras/stad2550
  * Marinacci et al. (2018) Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113, doi: 10.1093/mnras/sty2206
  * Marsan et al. (2022) Marsan, Z. C., Muzzin, A., Marchesini, D., et al. 2022, ApJ, 924, 25, doi: 10.3847/1538-4357/ac312a
  * Martis et al. (2016) Martis, N. S., Marchesini, D., Brammer, G. B., et al. 2016, ApJ, 827, L25, doi: 10.3847/2041-8205/827/2/L25
  * McConachie et al. (2022) McConachie, I., Wilson, G., Forrest, B., et al. 2022, ApJ, 926, 37, doi: 10.3847/1538-4357/ac2b9f
  * Merlin et al. (2019) Merlin, E., Fortuni, F., Torelli, M., et al. 2019, MNRAS, 490, 3309, doi: 10.1093/mnras/stz2615
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18, doi: 10.1088/0004-637X/777/1/18
  * Naiman et al. (2018) Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206, doi: 10.1093/mnras/sty618
  * Nanayakkara et al. (2022) Nanayakkara, T., Glazebrook, K., Jacobs, C., et al. 2022, arXiv e-prints, arXiv:2212.11638, doi: 10.48550/arXiv.2212.11638
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624, doi: 10.1093/mnras/stx3040
  * Pillepich et al. (2018) Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077, doi: 10.1093/mnras/stx2656
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Qin et al. (2017) Qin, Y., Mutch, S. J., Duffy, A. R., et al. 2017, MNRAS, 471, 4345, doi: 10.1093/mnras/stx1852
  * Schreiber et al. (2018) Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018, A&A, 618, A85, doi: 10.1051/0004-6361/201833070
  * Shimakawa et al. (2023) Shimakawa, R., Pérez-Martínez, J. M., Koyama, Y., et al. 2023, arXiv e-prints, arXiv:2306.06392, doi: 10.48550/arXiv.2306.06392
  * Spergel et al. (2003) Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175, doi: 10.1086/377226
  * Springel et al. (2005) Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629, doi: 10.1038/nature03597
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676, doi: 10.1093/mnras/stx3304
  * Tanaka et al. (2023) Tanaka, M., Onodera, M., Shimakawa, R., et al. 2023, arXiv e-prints, arXiv:2311.11569, doi: 10.48550/arXiv.2311.11569
  * Trinca et al. (2023) Trinca, A., Schneider, R., Valiante, R., et al. 2023, arXiv e-prints, arXiv:2305.04944, doi: 10.48550/arXiv.2305.04944
  * Valentino et al. (2020) Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93, doi: 10.3847/1538-4357/ab64dc
  * Valentino et al. (2023) Valentino, F., Brammer, G., Gould, K. M. L., et al. 2023, ApJ, 947, 20, doi: 10.3847/1538-4357/acbefa
  * Wang et al. (2023) Wang, Y.-Y., Lei, L., Yuan, G.-W., & Fan, Y.-Z. 2023, ApJ, 954, L48, doi: 10.3847/2041-8213/acf46c
  * Weaver et al. (2023) Weaver, J. R., Davidzon, I., Toft, S., et al. 2023, A&A, 677, A184, doi: 10.1051/0004-6361/202245581
  * Whitaker et al. (2011) Whitaker, K. E., Labbé, I., van Dokkum, P. G., et al. 2011, ApJ, 735, 86, doi: 10.1088/0004-637X/735/2/86
  * Williams et al. (2009) Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879, doi: 10.1088/0004-637X/691/2/1879
  * Wuyts et al. (2008) Wuyts, S., Labbé, I., Förster Schreiber, N. M., et al. 2008, ApJ, 682, 985, doi: 10.1086/588749
  * Xie et al. (2020) Xie, L., De Lucia, G., Hirschmann, M., & Fontanot, F. 2020, MNRAS, 498, 4327, doi: 10.1093/mnras/staa2370
  * Xie et al. (2017) Xie, L., De Lucia, G., Hirschmann, M., Fontanot, F., & Zoldan, A. 2017, MNRAS, 469, 968, doi: 10.1093/mnras/stx889