#  Superradiance in Stars:
Non-equilibrium approach to damping of fields in stellar media

F. Chadha-Day<EMAIL_ADDRESS>Institute for Particle
Physics Phenomenology, Department of Physics, Durham University, Durham DH1
3LE, United Kingdom  B. Garbrecht<EMAIL_ADDRESS>Technische Universität
München, Physik-Department, James-Franck-Straße, 85748 Garching, Germany  J.
I. McDonald<EMAIL_ADDRESS>Centre for Cosmology, Particle
Physics and Phenomenology, Université Catholique de Louvain, Chemin du
cyclotron 2, Louvain-la-Neuve B-1348, Belgium  Institute for Particle Physics
Phenomenology, Department of Physics, Durham University, Durham DH1 3LE,
United Kingdom

###### Abstract

Superradiance in black holes is well-understood but a general treatment for
superradiance in stars has until now been lacking. This is surprising given
the ease with which we can observe isolated neutron stars and the array of
signatures which would result from stellar superradiance. In this work, we
present the first systematic pipeline for computing superradiance rates in
rotating stars. Our method can be used with any Lagrangian describing the
interaction between the superradiant field and the constituents of the star.
Our scheme falls into two parts: firstly we show how field theory at finite
density can be used to express the absorption of long wavelength modes into
the star in terms of microphsyical scattering processes. This allows us to
derive a damped equation of motion for the bosonic field. We then feed this
into an effective theory for long wavelengths (the so-called worldline
formalism) to describe the amplification of superradiant modes of arbitrary
multipole moment for a rapidly rotating star. Our method places stellar
superradiance on a firm theoretical footing and allows the calculation of the
superradiance rate arising from any interaction between a bosonic field and
stellar matter.

Superradiance; neutron stars; astroparticle physics; axions

###### pacs:

95.35.+d; 14.80.Mz; 97.60.Jd

††preprint: IPPP/22/46

##  I Introduction

As we search for clues of new physics, astroparticle probes have in recent
years experienced an explosion of interest. It is becoming increasingly
apparent that physics beyond the standard model may lead to a vast array of
new and spectacular signatures arising from stars, black holes and galaxies
Berti  _et al._ (2022). This shift is driven both by the non-observation of
new physics at colliders (especially dark matter) and the advent of multi-
messenger astronomy, where the last decade has seen the emergence of
gravitational wave detectors and increasingly sophisticated X-ray and radio
telescopes.

In this regard superradiance Brito  _et al._ (2015); Bekenstein and Schiffer
(1998) sits perfectly at the intersection of many of these fields.
Superradiance is a general phenomenon which manifests itself in many systems.
In the present context, however, we are concerned with so-called rotational
superradiance. This variety of superradiance concerns the extraction of
rotational energy from fields scattering off spinning objects. This occurs
when low frequency waves are re-scattered with an amplitude larger than an
incident wave providing there is some mechanism for energy exchange with
object. This was originally illustrated by Zel’Dovich Zel’dovich (1971, 1972)
for the case of electromagnetic waves scattering off a rotating body with
finite conductivity (losses). This same mechanism manifests itself for
rotating Kerr black holes, where energy exchange is provided by dynamics
occurring close to the black hole111There is some discussion about whether or
not superradiant scattering requires a horizon, or simply the presence of an
ergo region Brito  _et al._ (2015); Vicente  _et al._ (2018); Eskin (2016);
Richartz  _et al._ (2009). For our purposes this is however a moot point,
since the ergoregion would be inside the star, where in any case one should
use an interior metric different from pure Kerr. We therefore investigate
other sources of direct dissipation coming from collisions of the matter in
the stellar interior.. Compact objects offer one additional ingredient, namely
that massive fields can become trapped around the black hole owing to the
existence of gravitational bound state solutions in black hole geometries.
This trapping and continual extraction of rotational energy leads to unstable
bound states for massive fields whose amplitude grows exponentially in time.
This is black hole superradiance. Given that superradiance will generically
occur for any rapidly rotating compact object, an immediate question is
whether the same mechanism will occur also in rapidly rotating neutron stars.

As with black holes, superradiance in stars would lead to a variety of
striking phenomena including gravitational waves Arvanitaki  _et al._ (2015);
Brito  _et al._ (2017a); Baryakhtar  _et al._ (2017); Brito  _et al._ (2017b),
polarisation signatures Chen  _et al._ (2020) and effects on binary inspirals
Baumann  _et al._ (2019). Superradiance can also be used to constrain the
existence of light fields Arvanitaki and Dubovsky (2011); Pani  _et al._
(2012); Cardoso  _et al._ (2018). It is also worth remarking on a recent
observation that the standard model photon itself may experience superradiance
growth around black holes owing to its finite plasma mass Blas and Witte
(2020); Wang  _et al._ (2022); Conlon and Herdeiro (2018)

Superradiance in stars has been studied for some specific cases Cardoso  _et
al._ (2015, 2017); Day and McDonald (2019) but a general treatment akin to
that for black holes has until now not been formulated. Here we present a
framework for computing the superradiance rates in stars given any fundamental
interaction.

Our key insight is to employ techniques from quantum field theory at finite
temperatures and densities. This allows us to describe, in a very general way,
an effective equation of motion for the bosonic field. This equation describes
the interaction between long-wavelengths of the superradiant field and the
constituents of the star, providing the missing link between microphysics and
damping terms that arise in the equations of motion for the boson.

This result is a satisfying illustration of the power of formal techniques in
field theory to answer questions in astrophysics. Such tools are more usually
brought to bear in the Early Universe to describe systems at finite
temperature and out-of-equilibrium (most usually in homogeneous but time-
varying environments). Applying these same procedures in an astrophysical
setting is quite novel. Here one must come to terms with the fact that the
star may provide a stationary, but not spatially invariant background.

Once an effective equation of motion is derived, the superradiance rate then
follows in a straightforward way by solving this equation for spherical modes
and mapping this onto superradiant bound states. Below we summarize our scheme
in more detail.

##  II Summary

The rate of stellar superradiance depends crucially on the damping rate of the
potentially superradiant field $\phi$ in the stellar background. In Section
III, we show how nonequilibrium field theory allows us to derive an effective
damped equation of motion for the scalar field

$\partial^{2}\phi+\mu^{2}\phi+\Gamma_{\phi}\dot{\phi}=0,$ (1)

where the damping is given by $\Gamma_{\phi}=-\lim_{p\rightarrow
0}\text{Im}\,\Pi^{R}(p)/p^{0}$ evaluated in the long-wavelength limit
$p\rightarrow 0$, where $\Pi^{R}$ is $\phi$’s in-medium self-energy

.

Despite the simple appearance of Eq. (1), it is worth emphasising that
demonstrating the validity of this result from first principles in non-
equilibrium quantum field theory and relating the damping rate to Green
functions and thence to collision integrals is a non-trivial and powerful
result, relying on a series of controlled approximations which are explained
step-by-step in later sections. This is testament to the power of fundamental
QFT and represents a novel application of finite density field theory to a
problem in astroparticle physics.

Below we summarise our scheme for computing superradiance rates in stars,
which consists of the following steps:

  * •

The first task is to calculate the imaginary part of the self energy that
yields $\Gamma_{\phi}$ per the above relation. This can be computed directly
within the framework of thermal field theory or, under approximations
discussed below, related via the optical theorem to a collision integral

$\displaystyle\text{Im}\,\Pi(q)=$
$\displaystyle\prod_{i}\int\frac{d^{3}\textbf{p}_{i}}{(2\pi)^{3}2E_{i}}f_{i}\prod_{j}\int\frac{d^{3}\textbf{p}_{j}}{(2\pi)^{3}2E_{j}}(1\pm
f_{j})\left|\mathcal{M}\right|^{2},$ (2)

where $f_{i}$ are distribution functions for particles in the medium and
$\mathcal{M}$ is a scattering amplitude.

From the classical equation of motion (1), one can derive an expression for
scattering of an arbitrary mode $(\ell,m,\omega)$ into a non-rotating star of
radius $R$

$Z_{\ell
m}\simeq-\frac{4R^{2}\omega(R\omega)^{2\ell+1}}{(1+2\ell)!!(2\ell+3)!!}\cdot\Gamma_{\phi}.$
(3)

  * •

The superradiant instability rate of a massive field gravitationally bound to
the star can then be calculated in the world line effective field theory
framework Endlich and Penco (2017), as described in Section IV. The scattering
described by equation (3) is essential in this process to match the effective
field theory to the low energy physics of the star. The superradiant growth
rate for harmonic bound state labeled by $(n,\ell,m)$

$\displaystyle\Gamma_{n\ell m}=C_{n\ell
m}\left(\frac{R}{r_{n\ell}}\right)^{(2\ell+3)}\frac{(m\Omega-\omega_{\ell
n})}{\omega_{\ell n}}\Gamma_{\phi},$ (4)

where $C_{n\ell m}$ is a combinatorial coefficient defined later on in Eq.
(72), $r_{n\ell}=(n+\ell+1)/(2GM\mu^{2})$ gives the characteristic radius of
the superradiant profile, and $\omega_{\ell n}$ gives the real part of the
discrete bound state frequencies Eq. (118). The rotational frequency of the
star is $\Omega$.

##  III Dissipation of Fields in QFT

Dissipation appears ubiquitously and is of importance for many physical
systems. This is also true for particle physics. While the phenomenon is
universal, there are many different calculational approaches, and the choice
of these can have various motivations. The following methods are of relevance
in the present discussion:

  * •

Given the utility of scattering theory in particle physics, inferring the
damping rate of a field from cross sections and the density of the medium is a
straightforward way that is often pursued. Limitations are met when other
medium effects such as corrections to the dispersion relations or screening
are of importance.

  * •

The Matsubara formalism for finite-temperature (equilibrium) field theory can
be applied to compute Green’s functions from which the damping rate and other
medium effects can be inferred. In the present context, it is of interest
whether such linear response coefficients also apply to the wave equation for
a classical field $\phi$.

  * •

The Schwinger-Keldysh closed-time-path (CTP) approach Keldysh (1964);
Schwinger (1961) allows one in principle to calculate the real time evolution
of the classical field and of correlation functions. In practice however, far-
reaching approximations and truncations must be applied to the equations that
directly descend from first principles. Such approximations result in
descriptions which can be solved numerically or analytically.

In the present work, we carry out some developments on the CTP approach in
order to apply these to the damping of scalar field waves. In Subsection
III.1, we review the CTP formalism, where the purpose is to declare the basic
conventions as well as to provide a very basic introduction to the matter.
Reviews and research articles that offer elaborations that are both, more
detailed and more synoptic, include Refs. Berges (2004); Calzetta and Hu
(2008, 1988); Prokopec  _et al._ (2004a). We then proceed to derive equations
of motion for a inhomogeneous scalar field $\phi$ in a general background to
linear approximation in $\phi$. These can be viewed as modified wave equations
that include dispersive and absorptive effects. Complementary to this are
studies where the evolution of a spatially homogeneous condensate is treated
to nonlinear level also using CTP methods Ai  _et al._ (2021a); Cheung  _et
al._ (2015)

In Subsection III.2, we carry out the transformation to Wigner space, which
allows to truncate the corrections to the wave equations so that these are
local, i.e. they do not feature memory integrals. In other words, the integro-
differential equations are approximated by the a wave equation in the usual
form of the partial differential equation (26). We discuss various additional
approximations, in particular the truncation of the gradient expansion, as
well.

###  III.1 Closed-time-path approach and equations of motion

When first introduced to the subject of quantum field theory one begins by
considering fluctuations about the quantum vacuum, where physical observables
are associated to the expectation values of some set of operators. The basic
objects in any field theory are the two-point functions, for instance, the
Wightman function defined by

$\braket{\varphi(x)\varphi(y)}_{\rm
vac.}=\bra{0,t\rightarrow\infty}\varphi(x)\varphi(y)\ket{0,t\rightarrow-\infty}$
(5)

where $\ket{0}$ is the ground state of the theory. Further, one assumes that
the one-point function of a generic field $\varphi$ will have a constant
expectation value $\phi(x)=\braket{\varphi(x)}=\text{const.}$, there are no
statistical fluctuations and the system is entirely described by a pure vacuum
state. Generally $\varphi$ may refer here to an arbitrary number of fields, so
that there may be flavour and Lorentz tensor or spinor structure that we
suppress for the time being.

By contrast, many physical systems are in a state that is quite different from
the vacuum $|0\rangle$. Of interest can be situations with a finite density of
particles, described in general by a mixed (possibly thermal) state. Fields
may attain non-trivial, spacetime-dependent, expectation values
$\phi=\braket{\varphi}\neq\text{const.}$ which vary throughout time and space.
The evolution of these expectation values can influence and be influenced by
the presence of other degrees of freedom. The latter may or may not be close
to thermal equilibrium. Clearly such systems require some self-consistent
treatment to address the interconnected dynamics of each field and its out-of-
equilibrium evolution.

The description of such systems rests on two fundamental points: that they
evolve unitarily according to a Hamiltonian $\mathcal{H}$ and that the state
of the system is described by a density matrix $\rho$. (Often, unitarity is
not respected by truncations. If the resulting equations describe e.g. the
state of a thermodynamic system they are typically irreversible.) In such
systems, generalized correlators are defined by

$\braket{\varphi(x)\varphi(y)}={\rm
tr}\left[\rho(t)\varphi(x)\varphi(y)\right]$ (6)

with similar definitions existing for the other Wightman function (with
exchanged coordinates), and the time-ordered (Feynman) and anti-time-ordered
(anti-Feynman) correlation functions. Such correlations ought to be
expressible via functional methods by taking derivatives of an appropriately
defined generating functional. Clearly the trace necessitates matrix elements
where the ket and bra states are the same - the so-called “in-in” scenario.
This is in contrast to Eq. (5) where the ingoing and outgoing states are
defined at different times. These are instead referred to as “in-out”
expectation values. Furthermore, we know that matrix elements can be re-
expressed in terms of path integrals. Matrix elements of the “in-in” type,
where both ket and bra states are evaluated at the same time $t$, can be
evaluated by inserting a complete basis of states at some intermediate time
$\tau$. The path integral can then be used to evolve forward from $t$ to
$\tau$ and then backwards again in the reverse direction Schwinger (1961);
Keldysh (1964), forming a closed loop in the time-domain, as illustrated in
FIG. 1. This corresponds to the following generating functional Calzetta and
Hu (1988, 2008)

$\displaystyle Z[J_{+},J_{-}]$ $\displaystyle=$
$\displaystyle\int\cal{D\varphi(\tau)}{\cal D}\varphi^{-}_{\rm in}{\cal
D}\varphi_{\rm in}^{+}\langle\varphi^{-}_{\rm
in}|\varphi(\tau)\rangle\langle\varphi(\tau)|\varphi_{\rm
in}^{+}\rangle\langle\varphi^{-}_{\rm in}|\varrho|\varphi_{\rm in}^{+}\rangle$
$\displaystyle=$ $\displaystyle\int{\cal D}\varphi^{-}{\cal D}\varphi^{+}{\rm
e}^{{\rm i}\int d^{4}x\\{{\cal L}(\varphi^{+})-{\cal
L}(\varphi^{-})+J_{+}\varphi^{+}-J_{-}\varphi^{-}\\}}\langle\varphi^{-}_{\rm
in}|\varrho|\varphi_{\rm in}^{+}\rangle\,.$ (7)

It generates all correlators, in particular the one-point function

$\phi=\braket{\varphi}=\pm\left.\frac{\delta Z[J_{+},J_{-}]}{\delta
J_{\pm}}\right|_{J_{\pm}=0}$ (8)

but also the set of four Green functions

$\displaystyle{\rm i}\Delta^{ab}(x,y)=$
$\displaystyle-\frac{\delta^{2}}{\delta J_{a}(x)\delta J_{b}(y)}\log
Z[J_{+},J_{-}]\Big{|}_{J_{\pm}=0}$ $\displaystyle=$ $\displaystyle{\rm
i}\langle{\cal C}[\varphi^{a}(x)\varphi^{b}(y)]\rangle\,,$ (9)

where we refer to $a,b=\pm$ as Schwiger–Keldysh indices and $\mathcal{C}$
defines a path ordering along the so-called Schwinger Keldysh contour in FIG.
1. We can then generate

$\displaystyle G^{>}(x,y)=$ $\displaystyle G^{-+}(x,y)\,,\quad
G^{<}(x,y)=G^{+-}(x,y)\,,$ $\displaystyle G^{T}(x,y)=$ $\displaystyle
G^{++}(x,y)\,,\quad G^{\bar{T}}(x,y)=G^{--}(x,y)\,,$ (10)

which are the two Wightman and time and anti-time ordered propagators,
respectively. In turn, linear combinations of these give the retarded and
advanced Green functions

$\displaystyle G^{\rm r}=$ $\displaystyle G^{T}-G^{<}=G^{>}-G^{T}\,,$
$\displaystyle\quad G^{\rm a}=$ $\displaystyle G^{T}-G^{>}=G^{<}-G^{\bar{T}}.$
(11)

We write here $G$ instead of $\Delta$ for generality because these definitions
apply as well for other two-point functions than correlators, in particular
for self energies.

Figure 1: The Schwinger Keldysh Contour. The closed path in the complex time
domain along which fields are evolved, corresponding to Eq. (7), where the
time $t$ corresponds to $x^{0}$.

Although we have presented a formalism for computing correlation functions for
fields out of equilibrium, we need to go one step further in actually
describing their real-time evolution through effective equations of motion.
This is especially true of the one-point function $\phi=\braket{\varphi}$,
whose equation of motion will eventually be identified with an effective
description for the dynamics of a superradiant field. What we need then is a
toolkit for describing the interconnected evolution of correlators. For this
we shall use the two-particle irreducible (2PI) effective action formalism,
which allows us to produce an intercoupled set of equations for one and two-
point functions.

The 2PI effective action Cornwall  _et al._ (1974); Calzetta and Hu (1988);
Berges (2004) has been established as a powerful formalism for providing a
description of the real time evolution of nonequilibrium systems, where the
relevant dynamics is typically described by two-point functions. For the
present work, the main target is the evolution of the one-point function
$\phi$, and backreaction to the finite density system is neglected. For this
reason, the relevant equations of motion could also be derived here from the
generating functional or the one-particle effective action. However, this
would not amount to a substantial simplification, and in view of applications
to nonequilibrium systems accounting for backreactions from $\phi$, we carry
out the formal developments here in the 2PI framework.

The 2PI effective action for a number of fields $\varphi$ and their
propagators $\Delta$ (where flavour, tensor, spinor and Schwinger–Keldysh
indices are suppressed) is given by

$\displaystyle\Gamma_{2\rm PI}[\phi,\Delta]$ $\displaystyle=$ $\displaystyle
S[\phi]+{\rm i}\,{\rm tr}[{\Delta^{(0)}}^{-1}\Delta]+{\rm i}\,{\rm
tr}\log\Delta^{-1}+\Gamma_{2}[\phi,\Delta].$ (12)

where $S$ is the classical action and ${\rm i}{\Delta^{(0)}}^{-1}=\delta
S/(\delta\phi\delta\phi)$ is the classical inverse propagator and $\Gamma_{2}$
is given by $(-{\rm i})$ times the sum of the 2PI vacuum graphs. The traces
are then to be taken over all indices, whether they pertain to flavour,
tensor, spinor structures or to the Schwinger–Keldysh contour. Note that
generally, neither ${\rm i}\Delta^{(0)}$ nor ${\rm i}\Delta$ are diagonal in
flavour but, of course, spacetime symmetry prohibits the mixing of fields from
different representations of the Lorentz algebra. Finally, the traces are
understood to act in complex field space, so that for a real scalar field, a
factor of $1/2$ is implied.

A closed set of equations of motion follows from the stationary points in the
functional dependence of $\Gamma_{2\rm PI}$ on $\phi$ and $\Delta$:

$\displaystyle\frac{\delta\Gamma_{2\rm
PI}[\phi,\Delta]}{\delta\phi(x)}=0,\qquad\frac{\delta\Gamma_{2\rm
PI}[\phi,\Delta]}{\delta\Delta(x,y)}=0.$ (13)

Again, all indices are suppressed here. Note that for the second of these
equations, there are four different Schwinger–Keldysh tensor components. These
can be decomposed into the celebrated Kadanoff–Baym equations as well as
equations for the retarded or advanced Green functions, see e.g. Ref. Prokopec
_et al._ (2004a).

In order to illustrate how dissipation of a classical field can result from
its interactions with another sector, we consider a simple model consisting of
a bosonic field $\varphi$ coupling to a massive fermion $\psi$ via

$\mathcal{L}=\frac{1}{2}\partial^{2}\varphi+V(\varphi)+i\bar{\psi}\not{\partial}\psi+m_{\psi}\bar{\psi}\psi+\lambda\varphi\bar{\psi}\psi.$
(14)

(To roll back to the above notation, refer to the scalar and the fermion field
collectively as $\varphi$, i.e. $\left\\{\varphi,\psi\right\\}\to\varphi$.)
This provides a toy proxy for the kind of interactions between a
(pseudo)scalar field and the constituents of the star which we will consider
in subsequent sections. It should also serve as an example to clarify how to
derive equations of motion for scalar fields subject to different
interactions.

We are principally interested in how the scalar dynamics is affected by the
presence of the second, fermionic field. This is described by the first
equation in (13), i.e. $\delta\Gamma/\delta\phi=0$, leading to

$\displaystyle\partial^{2}\phi(x)+V_{,\phi}(\phi(x))$
$\displaystyle+\frac{\delta\left(i{\rm tr}[{\Delta^{(0)}}^{-1}\Delta]+i{\rm
tr}[{{\cal S}^{(0)}}^{-1}{\cal S}]+\Gamma_{2}\right)}{\delta\phi}=0,$ (15)

where

$\displaystyle i\Delta^{(0)\,-1}$
$\displaystyle=\partial^{2}+V_{,\phi\phi}(\phi),$ (16) $\displaystyle i{\cal
S}^{(0)\,-1}$ $\displaystyle=i\not{\partial}+m_{\psi}+\lambda\phi$ (17)

and ${\cal S}$ and $\Delta$ are the full fermionic and bosonic propagators,
respectively. The term with the derivative $\delta/\delta\phi$ in Eq. (15)
introduces the effect of loop corrections on the dynamics of $\phi$. We take
the potential to be of the form

$V(\phi)=\frac{\mu^{2}}{2}\phi^{2}.$ (18)

From Eq. (16) we see then that $\Delta^{(0)\,-1}$ has no explicit $\phi$
dependence and is therefore annihilated by the action of the functional
derivative with respect to $\phi$. (Recall that the fields $\phi$ and the
propagators $\Delta$ are independent variables in the 2PI effective action.)
Furthermore, for the Yukawa interaction, insertions of $\phi$ into diagrams
are absorbed within the fermion propagator, cf. Eq. (17). Hence $\Gamma_{2}$
does not explicitly depend on $\phi$ and does not contribute to the derivative
term in Eq. (15). (Generally, when there are explicit insertions of $\phi$
into vertices of 2PI diagrams, the terms from $\Gamma_{2}$ enter the equations
of motion as well.) Therefore, in the present example, only the term ${\rm
tr}[{{\cal S}^{(0)}}^{-1}{\cal S}]$ survives the functional differentiation
with respect to $\phi$, which appears in ${{\cal S}^{(0)}}^{-1}$, leaving us
with

$\displaystyle\partial^{2}\phi(x)+\mu^{2}\phi(x)+\lambda i{\rm tr}[{\cal
S}(x,x)]=0\,.$ (19)

Now, while the fermion propagator $i{\cal S}$ does not explicitly depend on
$\phi$, it does so implicitly through its equation of motion in which $\phi$
appears. When $i{\rm tr}[{\cal S}(x,x)]_{\phi=0}\not=0$ (i.e. for constant,
vanishing field configurations), the fermion loop shifts the time-independent
solution for the expectation value of $\phi$. This contribution may be removed
by adding appropriate counter terms or simply by a redefinition of the field
$\phi$ by a constant value.

Assuming that the dynamics is adequately described by a linear
approximation,222At later stages of evolution, when field amplitudes become
large, non-linearities can enter, and bosenova collapse can occur for
superradiant fields Yoshino and Kodama (2012, 2014) in $\phi$, but for now we
are only interested in addressing superradiance with a linear stability
analysis. we functionally expand the implicit dependence of the solution for
$i{\cal S}$ to linear order in $\phi$. That is, we approximate

$\lambda\,{\rm tr}[{\cal S}(x,x)]=\int d^{4}y\Pi(x,y)\phi(y),$ (20)

where $\Pi$ is simply the scalar self energy evaluated in the background
$\phi(x)\equiv 0$. The resulting term in the wave equation is represented by
the Feynman diagram in Fig. 2. This leads to the following form of the wave
equation:

$\partial^{2}\phi^{a}(x)+\mu^{2}\phi^{a}(x)+\int
d^{4}y\,\Pi^{ab}(x,y)h^{bc}\phi^{c}(y)=0.$ (21)

where we have now made the Schwinger–Keldysh indices $a,b,c=\pm$ explicit,
with $h^{ab}$ defining a Schwinger-Keldysh metric given by $h={\rm
diag}(1,-1)$ which enforces the signs inherent in Eqs. (7)–(III.1). Taking
$a=+$, one can see that by using $\phi^{+}=\phi^{-}=\phi$ and the definitions
(III.1) and (11) generalized to self-energies, the final term can be expressed
as

$\partial^{2}\phi(x)+\mu^{2}\phi(x)+\int d^{4}y\,\Pi^{R}(x,y)\phi(y)=0,$ (22)

where $\Pi^{R}(x,y)=\Pi^{++}-\Pi^{+-}$ is the retarded self-energy. Equations
of this form are similar to those derived in Refs. Cheung  _et al._ (2015); Ai
_et al._ (2021a) for homogeneous field configurations. The fact that a memory
integral, i.e. a convolution with the retarded self energy, appears here
implies that this term in principle depends on values of $\phi(y)$ with
$y^{0}<x^{0}$. It may therefore be characterized as non-Markovian.

Figure 2: Tadpole diagram contributing to ${\rm tr}[S(x,x)]$ at linear order
in $\phi$. Double lines denote full fermion propagators in the event that
$\psi$ is coupled to some other sector different from $\phi$. Similarly, thick
dots may also represent vertex corrections from such a sector.

The fact that the convolution of $\phi$ with the retarded self energy appears
in Eq. (22) is a generic feature when the loop propagators lineraly depend on
$\phi$. Another important situation is the quadratic dependence, i.e. for a
quartic self interaction $\varphi^{4}$ or a coupling to another scalar $\chi$
via $\varphi^{2}\chi^{2}$. Then, tadpole corrections of the form

$\displaystyle\Pi^{\rm TP}(x)\phi(x)$

appear in the equations of motion. With the help of a $\delta$-function, these
can be of course brought to the form of the self energy term as in Eq. (22).

###  III.2 Wigner space and gradient expansion

It is apparent that Eq. (22) describes the effects of other degrees of freedom
on the evolution of the scalar field and that these are encoded in the self
energy. However, the equation in this form is still a little formal in several
aspects. Firstly, it is integro-differential which is impractical from the
point of view of finding analytic solutions. Secondly, the self-energy is
currently expressed in configuration space whilst we would intuitively expect
information about the effects of a fermionic medium on the scalar to be
expressible in terms of microphysical scattering processes, according to some
particle or quasiparticle approximation. Indeed, we would hope to make
quantitative predictions by calculating such processes and identifying them
with $S$-matrix elements. Finally, this equation in its present form does not
yet bear much resemblance to the usual form for a damped scalar field. It
turns out that all three of these issues can be solved by introducing the
notion of a local (for each space-time point) momentum space, as captured by
the Wigner transformation defined by

$\displaystyle\Pi^{R}(q,x)=\int d^{4}r\,e^{iq\cdot
r}\,\Pi^{R}\left(x+\frac{r}{2},x-\frac{r}{2}\right).$ (23)

The approximations will rely on the mild dependence of $\Pi^{R}(q,x)$ on the
average coordinate $x$. This is usually given by the macroscopic dependence on
the background that is characterized by small derivatives $\partial_{x}^{\mu}$
with respect to $x$. On the other hand, the dependence on the relative
coordinate corresponds to the quantum-physical correlations and is of order of
the intrinsic energy-momentum scales given by $q^{\mu}$. Provided these
generic assumptions hold, we can carry out a gradient expansion in terms of
$q\cdot\partial_{x}$.

Now, with this Wigner space description, we shall derive the damping term in
the classical wave equation from collision integrals in a medium, as
advertised in the Introduction. This connection appears as something that one
would very much expect intuitively, and certainly it has been applied in one
form or another in many instances. Nonetheless, we are not aware of a thorough
derivation, which is what we shall present now.

The Wigner transformation allows us to write

$\displaystyle\int d^{4}y\,\Pi^{R}(x,y)\phi(y)$ $\displaystyle=\int
d^{4}y\int\frac{d^{4}q}{(2\pi)^{4}}\Pi^{R}\left(q,\frac{x+y}{2}\right)\phi(y)e^{-iq\cdot(x-y)}.$

Our goal is to use the gradient expansion in order to reexpress this as a
series of local terms so that we obtain an ordinary differential equation
through its truncation. This is achieved through the Taylor expansion of
$\Pi^{R}(q,(x+y)/2)$ about $q=0$. After repeated use of integration by parts
(see Appendix A), one obtains the following expression for the convolution

$\displaystyle\int d^{4}y\,\Pi^{R}(x,y)\phi(y)$
$\displaystyle=\left[e^{i\partial_{q}\cdot\partial_{y}}\Pi^{R}\left(q,\frac{x+y}{2}\right)\phi(y)\right]_{q=0,y=x}\,.$
(25)

This allows us to express Eq. (22) terms of the Wigner-transform of the self
energy as

$\partial^{2}\phi+\mu^{2}\phi+\left[e^{i\partial_{q}\cdot\partial_{y}}\Pi^{R}\left(q,\frac{x+y}{2}\right)\phi(y)\right]_{q=0,x=y}\\!\\!\\!=0.$
(26)

Instead of an integro-differential equation, this is now in the form of an
ordinary differential equation of infinite order.

We note that for a real field $\phi$, ${\rm Re}[\Pi(q,x)]$ is an even function
under point reflections in the four momentum $p$, while the spectral self
energy

$\displaystyle\frac{i}{2}\left(\Pi^{>}(q,x)-\Pi^{<}(q,x)\right)=-{\rm
Im}[\Pi^{R}(q,x)]$ (27)

is odd. Therefore, Eq. (26) can be written as

$\displaystyle\partial^{2}\phi+\mu^{2}\phi$ $\displaystyle+$
$\displaystyle\left[\cos(\partial_{q}\cdot\partial_{y}){\rm
Re}\left[\Pi^{R}\left(q,\frac{x+y}{2}\right)\right]\phi(y)\right]_{q=0,x=y}$
$\displaystyle-$ $\displaystyle\left[\sin(\partial_{q}\cdot\partial_{y}){\rm
Im}\left[\Pi^{R}\left(q,\frac{x+y}{2}\right)\right]\phi(y)\right]_{q=0,x=y}$
$\displaystyle=$ $\displaystyle 0.$ (28)

In this form, it becomes explicit that the gradient expansion and Wigner space
representation of the self energy respect that the field $\phi$ is real.
Furthermore, the dispersive effects from ${\rm Re}[\Pi^{R}]$ and the
absorptive ones from ${\rm Im}[\Pi^{R}]$ are thus separated.

The wave equation for $\phi$ as in Eqs. (26,28) is still non-Markovian because
it is equivalent with Eq. (22). In general, the non-Markovian property of an
integro-differential equation implies a Wigner space description with an
infinite series of derivative operators. In many different situations, it is
however possible to truncate or to resum the derivatives so that one arrives
at an effectively Markovian description. This is because the boundary data for
a differential equation of finite order is fully specified by the field and an
finite number of its derivatives at a given time. A typical approximation is
to neglect the gradients $\partial_{y}\Pi^{R}$ when these correspond to
macroscopic scales and are thus small compared to the scales that determine
$\partial_{p}\Pi^{R}$, i.e. the density or temperature of the medium.

####  III.2.1 Damping in the long wavelength limit and superradiance

The gradients $\partial_{y}\phi$ may or may not be small. A generic situation
is however $(\partial_{q}\Pi^{R})\cdot(\partial_{y}\phi)\ll\Pi^{R}\phi$. This
typically occurs when the wave number $|\mathbf{k}|$ and frequency $|k^{0}|$
(corresponding to $\partial_{y}$) are small compared to the density or
temperature $T$ of the medium (corresponding to $\partial_{q}$)—which is the
case in the present application to stellar superradiance. Then, from Eq. (28)
it can be seen how damping enters the picture. When neglecting all derivatives
in $\partial_{q}\cdot\partial_{y}$ but the leading one, the contribution
involving ${\rm Im}[\Pi^{R}]$ in Eq. (28) takes the form of a classical
friction term in covariant form:

$\displaystyle\partial^{\mu}_{q}{\rm
Im}\left[\Pi^{R}(q,x)]\right]_{q=0}\,\partial_{\mu}\phi.$ (29)

We typically envisage a situation in which the damping arises as the result of
the absorption of a soft $\phi$ mode onto the external leg of an existing in-
medium scattering process—see Fig. 4. In this case, the leading four-momentum
transfer into the medium comes from $q^{0}$ component whilst three-momentum
shift q is sub-leading, so that the axion absorption is elastic with respect
to the three-momenta of the medium particles. This allows us to neglect
spatial momenta in $\Pi^{R}$ to leading order. This can be seen if one takes a
diagram in Fig. 4, where the axion momentum enters via the nucleon propagator
$1/[(p+q)^{2}-m^{2}]=1/[2p^{0}q^{0}-2\textbf{p}\cdot\textbf{q}+q^{2}]$, since
$p^{2}=m^{2}$ where $m$ and $p$ is the mass of an in-medium particle. Clearly
for non-relativistic particles, $p^{0}\simeq m\gg|\textbf{p}|$ so that the
process is dominated by $q^{0}$ with the q dependence being suppressed by the
heavy mass of the medium particles. As a result we can approximate

$\Pi^{R}(q^{0},\textbf{q},x)\simeq\Pi^{R}(q^{0},0,x).$ (30)

This means that the spatial derivative term in Eq. (29) is effectively
suppressed. Then, since ${\rm Im}[\Pi^{R}]$ is odd in $q^{0}$ it follows that
we can derive an effective damping rate

$\displaystyle\Gamma_{\phi}$
$\displaystyle=-\left.\partial_{q_{0}}\text{Im}\Pi^{R}(q^{0},0,x)\right|_{q^{0}\rightarrow
0}$ $\displaystyle=-\lim_{q^{0}\rightarrow
0}\frac{\text{Im}\,\Pi^{R}(q^{0},0,x)}{q^{0}}.$ (31)

leading to an equation

$\partial^{2}\phi+\mu^{2}\phi+\Gamma_{\phi}\dot{\phi}=0.$ (32)

####  III.2.2 Damping in the short wavelength limit (WKB)

For completeness of this discussion, we also consider the case when the
gradients $\partial_{y}\phi$ are not small, i.e. when the wave number
$|\mathbf{k}|$ or frequency $|k^{0}|$ of the field is of the same order or
large compared to the temperature $T$ or density. Another setup where these
gradients cannot be neglected is when the wave propagates in vacuum and the
invariant mass of $\phi$ is such that the decay into two or more particles is
possible. The corresponding momentum thresholds in ${\rm Im}[\Pi^{R}]$ can
then invalidate the truncation of the gradient expansion. For comparison, note
that for the convolution of two-point functions in Wigner space, large
derivatives because of large wave numbers or frequencies do not occur because
these are accounted for by the momentum variables in Wigner space. However, as
we shall see now, for approximate plane waves, the derivatives can be summed
to all orders so that once again an effectively Markovian description results.

In order to see that we recover in such situations from Eq. (28) the standard
corrections to the dispersion relations and the damping effects, assume for
now that the background is spatially homogeneous and time-independent, i.e.
$\Pi^{R}(q,x)\equiv\Pi^{R}(q)$. A relaxation of the homogeneity conditions may
be dealt with using the Wentzel–Kramers–Brillouin approximation (the
applicability of which will depend on the specific problem).

Then, the spatial part of the solutions can be expanded in a basis of momentum
eigenstates $\phi(x)\sim\exp(i\mathbf{k}\cdot\mathbf{x})$. When summing the
spatial derivatives to all orders, we recognize a Taylor expansion of the self
energy (note that this relies on the exponential form of the solution) so that
Eq. (26) becomes

$\displaystyle(\partial_{t}^{2}+\omega^{2})\phi(t,\mathbf{x})+\left[e^{i\partial_{q^{0}}\partial_{y}^{0}}\Pi^{R}(q^{0},\mathbf{k})\phi(y^{0})\right]_{q^{0}=0,y^{0}=t}=0,$
(33)

where $\omega=(\mathbf{k}^{2}+m^{2})^{1/2}$.

Neglecting a possible dependence of $\Pi^{R}$ on $\phi$, this equation is
linear so that the solution is of the form
$\phi(x)=\exp(-i\omega^{(1)}t+i\mathbf{k}\cdot\mathbf{x})$, where the notation
$\omega^{(1)}$ indicates that i corresponds to a correction to $\omega$ due to
$\Pi^{R}$. Then, also the temporal gradients can be summed, what leads to

$\displaystyle(\partial_{t}^{2}+\omega^{2})e^{-i\omega^{(1)}t}+\Pi^{R}(\omega^{(1)},\mathbf{k})e^{-i\omega^{(1)}}t=0$
(34)

and therefore

$\displaystyle{\omega^{(1)}}^{2}=\omega^{2}+\Pi^{R}(\omega,\mathbf{k}).$ (35)

In the argument of $\Pi^{R}$, we have replaced $\omega^{(1)}$ with $\omega$.
This leads to a difference at second order in the self energy, while we are
aiming here to consistently include the first order effects.

##  IV Superradiance with worldline effective field theory

To describe stellar superradiance, we are interested in the interaction of our
field with a rotating star. One strategy would be to attempt to include
rotation in the formalism described in Section III by trying to introduce a
bulk azimuthal fluid velocity into the distribution functions. It remains to
be seen how this would be realised in detail. Instead, in this work, we will
make use of the worldline effective field theory (EFT) approach developed in
Endlich and Penco (2017). This will simplify the calculations and also gives
insight into the physics of superradiance. This provides a relation between
the absorption of modes into a static non-rotating star and the amplification
of modes into a rotating star, with the transformation being inferred from an
appropriate transformation between a rotating and non-rotating frame.

As with all EFTs we must exploit a hierarchy of scales. In the present case,
the high energy scale is set by the radius of the star or black hole $\sim R$,
and the low energy scale corresponds to wavelengths $\sim 1/k$ of the
superradiant field which are large compared to this scale. The EFT is then an
effective expansion in the parameter $R/k^{-1}\ll 1$.

In this picture, the extended object is essentially integrated out, leaving a
power series of contact interactions between an effective point-like object
and the field around it. In this way, finite size effects are then described
by an infinite series of effective operators that act in the worldline of the
object. The symmetries of these operators are determined by the symmetries of
the object. As described in Delacrétaz  _et al._ (2014), a spherically
symmetric object is described by operators transforming under $SO(3)$.

For example, for an extended object interacting with a scalar field $\phi$,
the interaction Hamiltonian in the worldline EFT is

$\displaystyle H_{\rm int}(t,{\bf x})=$
$\displaystyle\partial^{I}\phi(x)\mathcal{O}^{(1)}_{I}(x)\delta^{(3)}({\bf
x}-{\bf y}(t))+$ (36)
$\displaystyle\partial^{I}\partial^{J}\phi(x)\mathcal{O}^{(2)}_{IJ}(x)\delta^{(3)}({\bf
x}-{\bf y}(t))+\dots,$

where ${\bf y}(t)$ is the worldline of the extended object and
$\mathcal{O}^{(n)}$ are effective operators describing the interaction between
the object and the scalar field. $I,J$ are spatial indices, so that the
operators $\mathcal{O}^{(n)}$ are not Lorentz invariant. This arises because a
macroscopic object breaks Lorentz symmetry through its rest frame, location
and orientation in space.

The worldline EFT has many applications, including in black hole superradiance
Endlich and Penco (2017). One powerful aspect is the ease with which rotation
of the extended object can be added. Firstly, we can integrate the Hamiltonian
over space to project the Hamiltonian onto the worldline of the star:

$\displaystyle H_{\rm int}(t)$ $\displaystyle=\int d^{3}{\bf x}H_{\rm
int}(t,{\bf x})$
$\displaystyle=\partial^{I}\phi(t)\mathcal{O}^{(1)}_{I}(t)+\partial^{I}\partial^{J}\phi(t)\mathcal{O}^{(2)}_{IJ}(t)+...,$
(37)

where the fields are now understood to be evaluated at the spatial position of
the object which we will take to be ${\bf x}=0$. The EFT associated to a
rotating object can be inferred through the action of a time-dependent
rotation matrix on the operators $\mathcal{O}^{n}$

$\mathcal{O}^{(n)}_{I_{1}\cdots
I_{n}}(t)\rightarrow\tensor{R}{{}^{J_{1}}_{I_{1}}}(t)\cdots\tensor{R}{{}^{J_{\ell}}_{I_{n}}}(t)\mathcal{O}^{(n)}_{J_{1}\cdots
J_{n}}(t),$ (38)

where, without loss of generality we have chosen coordinates in which the
object spins about the $z$-axis so that

$R(t)=\left(\begin{array}[]{ccc}\cos(\Omega t)&-\sin(\Omega t)&0\\\
\sin(\Omega t)&\cos(\Omega t)&0\\\ 0&0&1\end{array}\right).$ (39)

The EFT for the field interacting with a rotating object, as as seen by an
inertial (non-rotating), observer then reads

$\displaystyle H_{\rm int}(t)$
$\displaystyle=\partial^{I}\phi(t)\tensor{R}{{}_{I}^{J}}(t)\mathcal{O}^{(1)}_{J}(t)$
$\displaystyle+\partial^{I}\partial^{J}\phi(t)\tensor{R}{{}_{I}^{K}}(t)\tensor{R}{{}_{J}^{L}}(t)\mathcal{O}^{(2)}_{KL}(t)+...,$
(40)

In the following sections we will show how this EFT can be used to describe
superradiant scattering and superradiant instabilities of bound modes around a
rotating object.

###  IV.1 Superradiant scattering

First, let us imagine the interaction Hamiltonian is switched on over some
finite time $\mathcal{T}$. The evolution operator is then given by

$\displaystyle S$ $\displaystyle=T\exp\left\\{-i\int_{0}^{\mathcal{T}}dt\int
d^{3}{\bf x}H_{\mathrm{int}}(t,{\bf x})\right\\}$ (41)
$\displaystyle=T\exp\left\\{-i\int_{0}^{\mathcal{T}}dtH_{\mathrm{int}}(t)\right\\}.$
(42)

Let us now consider the effect of this interaction on the absorption of the
field $\phi$ by the rotating object, following the method introduced in
Endlich and Penco (2017). We compute the absorption of a spherical mode of the
field $\phi$ with frequency $\omega$ and angular momentum labels $\ell,m$. The
probability for this process is:

$P_{\mathrm{abs}}=\sum_{X_{f}}\frac{\left|\left\langle
X_{f};0|S|X_{i};\omega,\ell,m\right\rangle\right|^{2}}{\langle\omega,\ell,m\mid\omega,\ell,m\rangle},$
(43)

where $X_{i}$ and $X_{f}$ are the initial and final states of the star. Note
that the leading order contribution to scattering of a mode with angular
momentum $\ell$ comes from the operator

$H_{\ell}(t)=\partial^{I_{1}}\cdots\partial^{I_{\ell}}\phi\tensor{R}{{}_{I_{1}}^{J_{1}}}\cdots\tensor{R}{{}_{I_{\ell}}^{J_{\ell}}}O_{J_{1}...J_{\ell}},$
(44)

in Eq. (36). Upon inserting (44) into (41) and expanding to leading order in
interactions, one obtains

$\displaystyle\sum_{X_{f}}\left|\left\langle
X_{f};0|S|X_{i};\omega,\ell,m\right\rangle\right|^{2}=$
$\displaystyle\int_{0}^{\mathcal{T}}dt\int_{0}^{\mathcal{T}}dt^{\prime}$
$\displaystyle\bra{0}\partial^{I_{1}}\cdots\partial^{I_{\ell}}\phi(t^{\prime})\tensor{R}{{}_{I_{1}}^{J_{1}}}(t^{\prime})\cdots\tensor{R}{{}_{I_{\ell}}^{J_{\ell}}}(t^{\prime})\ket{\omega,\ell,m}$
$\displaystyle\bra{\omega,\ell,m}\partial^{K_{1}}\cdots\partial^{K_{\ell}}\phi(t)\tensor{R}{{}_{K_{1}}^{L_{1}}}(t)\cdots\tensor{R}{{}_{K_{\ell}}^{L_{\ell}}}(t)\ket{0}$
$\displaystyle\braket{O_{J_{1}...J_{\ell}}(t^{\prime})O_{L_{1}...L_{\ell}}(t)},$
(45)

where $\braket{O_{J_{1}...J_{\ell}}(t^{\prime})O_{L_{1}...L_{\ell}}(t)}$
denotes the expectation value for the initial state $|X_{i}\rangle$. All that
remains is to manipulate the matrix elements associated to the spherical modes
of $\phi$. This is a lengthy but straightforward computation, the details of
which can be found in Endlich and Penco (2017). The result is

$\displaystyle\sum_{X_{f}}\left|\left\langle
X_{f};0|S|X_{i};\omega,\ell,m\right\rangle\right|^{2}$ (46)
$\displaystyle=\frac{\ell!q^{2\ell+2}}{2\pi
v\omega(2\ell+1)!!}\int\frac{d\omega^{\prime}}{2\pi}\Delta_{\ell}(\omega^{\prime})\left|\int_{0}^{\mathcal{T}}dt{\rm
e}^{i\omega^{\prime}t}{\rm e}^{-i(\omega-m\Omega)t}\right|^{2},$ (47)

where $q=|\textbf{q}|$ and $v=\sqrt{1-\mu^{2}/\omega^{2}}$ is the group
velocity of a field mass $\mu$ asymptotically far from the star and
$\Delta_{\ell}(\omega)$ is the Fourier transform of the Wightman correlation
function defined by

$\braket{O_{J_{1}...J_{\ell}}(t^{\prime})O_{L_{1}...L_{\ell}}(t)}=\delta_{J_{1}...J_{\ell}}^{L_{1}...L_{\ell}}\int\frac{d\omega}{2\pi}\Delta_{\ell}(\omega){\rm
e}^{i\omega(t-t^{\prime})}.$ (48)

The above form is permitted by isotropy (spherical symmetry of the star),
leading to delta function prefactors, and the fact the background is
stationary so that time-translation symmetry means the correlators can only
dependent on the relative coordinate $t-t^{\prime}$. For sufficiently long
times the $t$ integral can be approximated by a delta function, multiplied by
$\mathcal{T}$ leading to

$\displaystyle\sum_{X_{f}}\left|\left\langle
X_{f};0|S|X_{i};\omega,\ell,m\right\rangle\right|^{2}$
$\displaystyle=\frac{\ell!q^{2\ell+2}}{2\pi
v\omega(2\ell+1)!!}\int\frac{d\omega^{\prime}}{2\pi}2\pi\mathcal{T}\delta_{\ell}(\omega^{\prime}-(\omega-m\Omega))\Delta_{\ell}(\omega^{\prime})$
$\displaystyle=\frac{\ell!q^{2\ell+2}}{2\pi
v\omega(2\ell+1)!!}\Delta_{\ell}(\omega-m\Omega)\mathcal{T}.$ (49)

We now turn to the normalisation of states in the denominator, which, as we
show in appendix B is normalised to 333Note that our treatment of limits in
deriving the absorption and emission rates differs slightly from that in
Endlich and Penco (2017). In Endlich and Penco (2017), the $\delta(0)$ terms
corresponding to the infinite time and infinite radius limits are incorrectly
canceled. This error does not affect their conclusions.

$\braket{\omega\ell m\left missing}{\omega\ell m}=\frac{2\mathcal{R}}{v},$
(50)

where $\mathcal{R}$ is a finite spherical radius used to regulate the inner
product of non-normalised states. This is the spherical equivalent of the
usual replacement $(2\pi)^{3}\delta^{(3)}({\bf k}=0)\rightarrow V$. Hence we
can write

$P_{\rm
abs}=\frac{\ell!q^{2\ell+2}}{4\pi\omega(2\ell+1)!!}\Delta_{\ell}(\omega-m\Omega)\frac{\mathcal{T}}{\mathcal{R}}.$
(51)

Dividing by $\mathcal{T}$, we then obtain an absorption rate

$R_{\rm
abs}=\frac{\ell!q^{2\ell+2}}{4\pi\omega(2\ell+1)!!}\Delta_{\ell}(\omega-m\Omega)n_{r},$
(52)

where $n_{r}=1/\mathcal{R}$ is the particle number per spherical shell of unit
radius. We can similarly calculate an emission probability

$R_{\rm
em}=\frac{\ell!q^{2\ell+2}}{4\pi\omega(2\ell+1)!!}\Delta_{\ell}(m\Omega-\omega)n_{r}.$
(53)

noting that the argument of $\Delta_{\ell}$ appears now with the opposite sign
to Eq. (51). We can then define an amplification factor for spherical waves as

$\displaystyle Z_{\ell m}$ $\displaystyle=\frac{\Phi_{\rm out}-\Phi_{\rm
in}}{\Phi_{\rm in}}$ $\displaystyle=\frac{R_{\rm em}-R_{\rm abs}}{n_{r}v}$
$\displaystyle=\frac{\ell!\,q^{2\ell+2}}{4\pi(2\ell+1)!!v\omega}\,\rho_{\ell}(m\Omega-\omega),\,$
(54)

where $\Phi_{\rm in}=n_{r}v$ is the ingoing flux, which is related to the
outgoing flux by $\Phi_{\rm out}=\Phi_{\rm in}[1+P_{\rm em}-P_{\rm ab}]$ and
where

$\rho_{\ell}(\omega)=\Delta_{\ell}(\omega)-\Delta_{\ell}(-\omega),$ (55)

is the spectral response function for the star. The key point is that
$\rho_{\ell}(\omega)$ is an odd function of $\omega$, hence for frequencies
$\omega>m\Omega$ energy is absorbed ($Z_{\ell m}<0$) whilst for
$\omega>m\Omega$ energy is extracted from the star ($Z_{\ell m}>0$). This is
superradiant scattering.

###  IV.2 Superradiant instabilities

We now turn to a similar question of bound states, extending the method of
Endlich and Penco (2017). It is precisely these bound states, which, due to
their gravitational confinement, exhibit instabilities through the
superradiant extraction of rotational energy. Consider bound states of a
spin-0 field with mass $\mu$. The scalar field operator can be rewritten in
terms of bound state creation and annihilation operators:

$\displaystyle\hat{\phi}=\sum_{n\ell
m}\frac{1}{\sqrt{2\omega_{n\ell}}}\Big{(}$ $\displaystyle\hat{a}_{n\ell
m}f_{n\ell m}(r,\theta,\phi){\rm e}^{-i\omega_{n\ell}t}$
$\displaystyle+\hat{a}^{\dagger}_{n\ell m}f^{\star}_{n\ell
m}(r,\theta,\phi){\rm e}^{i\omega_{n\ell}t}\Big{)},$ (56)

where $a_{n\ell m}$ are creation operators for a mode $(n,\ell,m)$. The mode
functions

$f_{n\ell m}(r,\theta,\phi)=R_{n\ell}(r)Y_{\ell m}(\theta,\phi),$ (57)

are the bound state solutions on a Schwarzschild black-hole background.
Expressions for $R_{n\ell}$ are given in appendix C. They are normalised as:

$\int d\Omega drr^{2}f_{n\ell
m}(r,\theta,\phi)f^{\star}_{n^{\prime}\ell^{\prime}m^{\prime}}(r,\theta,\phi)=\delta_{nn^{\prime}}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}.$
(58)

$n$ and $l$ take non-negative integer values and $m$ takes integer values from
$-l$ to $l$. We wish to calculate the absorption and emission probability of a
bound state $\ket{n,\ell,m}=\hat{a}^{\dagger}_{n\ell m}\ket{0}$ by the star.
We have:

$P_{\rm abs}=\sum_{X_{f}}|\bra{X_{f};0}S\ket{X_{i};n\ell m}|^{2},$ (59)

where $X_{i}$ and $X_{f}$ are the initial and final states of the star. Note
that the bound states are normalised to 1. As before, the dominant
contribution to this process comes from the interaction of Eq. (44), leading
to

$\displaystyle P_{\rm abs}\simeq\sum_{X_{f}}\int^{\mathcal{T}}_{0}$
$\displaystyle
dtdt^{\prime}\bra{X_{i};n,\ell,m}H_{\ell}(t^{\prime})\ket{X_{f};0}$ (60)
$\displaystyle\bra{X_{f};0}H_{\ell}(t)\ket{X_{i};n,\ell,m}.$

Whilst the general $\ell$ case for superradiant scattering has been covered in
Endlich and Penco (2017), to which we referred the reader earlier in
subsection IV.1, the general $\ell$ bound state case was not presented there.
We therefore collect the relevant details in Appendix C. There, we show that
for large $\mathcal{T}$, one finds

$\displaystyle P_{\rm abs}=$ $\displaystyle A_{n\ell
m}\left(\frac{1}{r_{n\ell}}\right)^{2\ell+3}\int\frac{d\omega}{2\pi}\frac{\Delta_{\ell}(\omega)}{2\omega_{\ell
n}}\mathcal{T}(2\pi)\delta(m\Omega+\omega-\omega_{n\ell}),$ (61)

where $r_{n\ell}=(n+\ell+1)/(2GM\mu^{2})$ and $\omega_{\ell n}$ are the
discrete bound state frequencies defined in Eq. (118) and

$A_{n\ell
m}=\frac{\ell!(2\ell+1)\text{!!}(2\ell+n+1)!}{8\pi(\ell+n+1)n!(2\ell+1)!^{2}}.$
(62)

We then obtain the absorption rate

$\Gamma_{\rm abs}=A_{n\ell
m}\left(\frac{1}{r_{n\ell}}\right)^{2\ell+3}\frac{1}{2\omega_{\ell
n}}\Delta_{\ell}(\omega_{\ell n}-m\Omega).$ (63)

Similarly, the emission rate is

$\Gamma_{\rm em}=A_{n\ell
m}\left(\frac{1}{r_{n\ell}}\right)^{2\ell+3}\frac{1}{2\omega_{\ell
n}}\Delta_{\ell}(m\Omega-\omega_{\ell n}),$ (64)

and so the superradiance rate is simply the difference of these two quantities

$\displaystyle\Gamma_{n\ell m}$ $\displaystyle=\frac{A_{n\ell
m}}{2\omega_{\ell
n}}\left(\frac{1}{r_{n\ell}}\right)^{2\ell+3}\rho_{\ell}(m\Omega-\omega_{\ell
n}).$ (65)

As we have seen in the previous section, and as described in Endlich and Penco
(2017), we can obtain $\rho_{\ell}$ by matching to the superradiant scattering
amplitude (IV.1), so that scattering off a static star leads immediately to
the superradiant growth rate of bound states. By definition, gravitational
bound states have frequencies $\omega<\mu$ below the mass of the field. Hence
superradiance requires knowledge of $\rho_{\ell}(\omega<\mu)$ and since a
massive propagating wave has $\omega>\mu$, it follows that carrying out
scattering of massless modes with $\omega<\mu$ is sufficient to extract the
spectral response at arbitrarily small arguments. This can then be used to
infer the superradiance rate for massive fields which requires $\rho(\omega)$
in the region $\omega<\mu$. We illustrate this in the next subsection for a
uniform star model.

The crucial point of this section is that correlators
$\braket{O_{J_{1}...J_{\ell}}(t^{\prime})O_{L_{1}...L_{\ell}}(t)}$ encode all
the information about the response of a point-like star to a harmonic stimulus
$\sim e^{-i\omega t}$ from the external field which is essentially applied
homogeneously across the star in the long wavelength limit. Hence knowledge of
these correlators, as embodied in $\rho$ is sufficient to calculate
superradiance rates and this can be extracted from any physical process to
which $\rho$ contributes.

###  IV.3 Superradiance rate from microphysics

As an example, let us consider superradiance for a scalar field around a
uniform star. To obtain the spectral function $\rho$, we will calculate the
scattering of a massless scalar field $\phi$ from the star when it is not
rotating. For a uniform star the self-energy $\Pi^{R}$ is constant inside
$r\leq R$ and zero elsewhere. The equation of motion for $\phi$ is given by
equation (32). Under these approximations, we can compute the absorption of a
mode $(\ell,m)$ by considering solutions of the form

$\phi=\psi(r)Y_{\ell m}(\theta,\varphi)e^{-i\omega t},$ (66)

which gives a radial equation

$r^{2}\psi^{\prime\prime}+2r\psi^{\prime}+\left[r^{2}(\omega^{2}+i\omega\Gamma_{\phi})-\ell(\ell+1)\right]\psi=0.$
(67)

We recognise this immediately as the equation for spherical Bessels functions
which admits two linearly-independent solutions. Inside the star, by
regularity at the origin, the solution must take the form
$j_{\ell}\left(r\sqrt{\omega^{2}+i\omega\Gamma_{\phi}}\right)$ whilst outside
it can be most conveniently parameterised in terms of spherical hankel
functions $h_{\ell}^{(1,2)}(\omega r)$. Hence at infinity, the solution takes
the form

$\psi(r\rightarrow\infty)=\phi_{\rm out}\frac{e^{i\omega r}}{r}+\phi_{\rm
in}\frac{e^{-i\omega r}}{r}.$ (68)

Imposing continuity of the solution and its first derivative at $r=R$ gives
expressions for $\phi_{\rm in}$ and $\rm\phi_{\rm out}$. One can then define
an amplification factor for a non-rotating star in flat space:

$\displaystyle Z_{\ell m}=\frac{|\phi_{\rm out}|^{2}}{|\phi_{\rm in}|^{2}}-1.$
(69)

Using the above continuity conditions, expanding in the limit $\omega R\ll 1$
and to linear order in $\Gamma_{\phi}$, we find (recalling we deal first with
a non-rotating star)

$Z_{\ell
m}\simeq-\frac{4R^{2}\omega(R\omega)^{2\ell+1}}{(1+2\ell)!!(2\ell+3)!!}\cdot\Gamma_{\phi}$
(70)

repeating the expression found in Cardoso  _et al._ (2015). We then equate the
amplification factors in equations (70) and (IV.1) for a massless scalar field
scattering from a non-rotating star. In that limit we have $v=1$, $k=\omega$,
which allows us to extract an expression for the spectral function
$\rho(\omega)$ which can then be substituted into the superradiance rate (65)
leading to

$\displaystyle\Gamma_{n\ell m}=C_{n\ell
m}\left(\frac{R}{r_{n\ell}}\right)^{(2\ell+3)}\frac{(m\Omega-\omega)}{\omega}\Gamma_{\phi},$
(71)

where

$C_{n\ell
m}=\frac{(2\ell+1)\text{!!}(2\ell+n+1)!}{(\ell+n+1)n!(2\ell+1)!^{2}(2\ell+3)!!}.$
(72)

It is interesting to compare the amplification factor for a rotating star

$Z^{\rm rot.}_{\ell
m}\simeq\frac{4R^{2}(m\Omega-\omega)(R\omega)^{2\ell+1}}{(1+2\ell)!!(2\ell+3)!!}\cdot\Gamma_{\phi}$
(73)

inferred from Eq. (IV.1), with that of a black hole Brito  _et al._ (2015)

$\displaystyle Z^{\rm BH}_{lm}$ $\displaystyle=-8GMr_{+}(\omega-m\Omega_{\rm
H})\omega^{2\ell+1}\left(r_{+}-r_{-}\right)^{2l}$
$\displaystyle\times\left[\frac{(\ell!)^{2}}{(2\ell)!(2\ell+1)!!}\right]^{2}\prod_{k=1}^{l}\left[1+\frac{(GM)^{2}}{k^{2}}\left(\frac{\omega-m\Omega_{\rm
H}}{\pi r_{+}T_{H}}\right)^{2}\right]$ (74)

where $r_{\pm}=GM(1\pm\sqrt{1-\tilde{a}^{2}})$ where $\tilde{a}=a/(GM)$ and
$T_{H}=(r_{+}-r_{-})/(4\pi r_{+}^{2})$ is the BH temperature. Expanding this
to leading order in spin and the splitting $(\omega-m\Omega_{H})$, we find

$\displaystyle Z^{\rm BH}_{lm}$ $\displaystyle\simeq 4r_{s}^{2}(m\Omega_{\rm
H}-\omega)\left(\omega
r_{s}\right)^{2\ell+1}\left[\frac{(\ell!)^{2}}{(2l)!(2\ell+1)!!}\right]^{2}\cdot\frac{1}{r_{s}}\
$ (75)

where $r_{s}=2GM$ is the Schawrzschild radius. Hence up to combinatorial
factors, we can see by comparing (71) and (75) that there is a one-to-one
correspondence between the superradiant scattering for stars and black holes
given by replacing the Schwarzshild radius, $r_{s}$ with the Neutron star
radius $\rightarrow R$ and identifying the damping $\Gamma_{\phi}$ with
$1/r_{s}$. Hence we see that superradiance in stars will become comparable to
that of stellar mass black holes when the mean free path $\lambda$, associated
to the microphysics $\lambda^{-1}\simeq\Gamma_{\phi}$, is less than or equal
to the gravitational radius of the black hole.

##  V Axions

We shall now consider the case where the bosonic field is an axion-like
particle (ALP) which couples to nuclear matter within the star. Ultra-light
bosonic fields occur in many well-motivated extensions of the Standard Model.
In particular, axions are an attractive dark matter candidate Kim (1979);
Shifman  _et al._ (1980); Dine  _et al._ (1981); Zhitnitsky (1980), as well
offering a solution to the strong CP problem444Note some recent debate over
whether the $\theta$ term actually appears in observables Ai  _et al._
(2021b). Peccei and Quinn (1977); Weinberg (1978); Wilczek (1978) and arising
generically in string compactifications Svrcek and Witten (2006); Arvanitaki
_et al._ (2010). Indeed in recent years, neutron stars have become an exciting
environment in which to study axions Pshirkov and Popov (2009); Garbrecht and
McDonald (2018); Hook  _et al._ (2018); Huang  _et al._ (2018); Battye  _et
al._ (2020); Leroy  _et al._ (2020); Darling (2020a, b); Witte  _et al._
(2021); Battye  _et al._ (2021); Millar  _et al._ (2021); Battye  _et al._
(2022); Foster  _et al._ (2022)

We shall model the nuclear interactions using chiral perturbation theory. In
this instance, the principle absorption channel is $a+N+N\rightarrow N+N$,
which proceeds via the interactions Brinkmann and Turner (1988)

$\displaystyle\mathcal{L}_{aNN}$
$\displaystyle=\frac{g_{an}}{2m}\partial_{\mu}a\bar{N}\gamma^{\mu}\gamma_{5}N,$
(76) $\displaystyle\mathcal{L}_{\pi NN}$
$\displaystyle=\frac{2m}{m_{\pi}}f\pi_{0}\bar{N}\gamma^{5}N,$ (77)

where $f\approx 1$ and $N=p,n$ are the proton or neutron.

Figure 3: Axion self energies whose optical cuts give rise to
$N+N+a\rightarrow N+N$ processes corresponding to ${\rm Im}\,\Pi^{R}$.

Note that we have chosen a basis in which the axion-nucleon interactions are
“pseudo-vector” (PV) (axion derivatives coupling to the axial nucleon current)
and the neutral pion coupling is “pseudo-scalar” (PS) (non-derivative coupling
of the pion to nucleons). Since the pion and axion are pseudo Nambu-Goldstone
bosons, we can always perform chiral transformations, leading to different
combinations of PV and PS couplings to nucleons. This was demonstrated
explicitly in Choi  _et al._ (1989) which illustrated the various combinations
of PS and PV couplings which can be achieved. In particular, our
parametrisation (76)-(77) allows us to make contact with the classic reference
Brinkmann and Turner (1988) in which axion-nucleon absorption processes are
considered.

We shall now see how the microphysical scattering processes which contribute
to $\Pi^{R}$ and thence to the damping term $\Gamma_{a}$ of the axion emerge
naturally from cuts through an appropriate self-energy, in accordance with the
optical theorem. For the interactions (76)-(77) there are four such diagrams
shown in figure 3 whose cuts lead to processes $a+N+N\leftrightarrow N+N$. We
do not consider diagrams at lower order in $\mathcal{L}_{\pi NN}$ whose cuts
would lead to $a+N\leftrightarrow N+\pi$ processes. These scattering do not
occur for the low energy axion modes we consider as the pion is too heavy.

As explained in Section III, the spectral self-energy (which encodes
information about dissipation) can be easily obtained once either of the
Wightman self-energies is known (see Eq. (27)). Let us therefore compute
$\Pi^{>}$, which is obtained by summing the internal vertices in the diagram
over Schwinger-Keldysh indices. The internal propagators connecting different
$\pm$ Schwinger-Keldysh therefore each carry $\pm$ labels. For Fermions, in
momentum space, these free thermal propagators are given by (following the
conventions of Prokopec  _et al._ (2004a, b))

$\displaystyle i{\cal S}^{++}(p)$
$\displaystyle=\phantom{-}\frac{i(\not{p}+m)}{p^{2}-m^{2}+{\rm
sign}(p^{0})\epsilon}$ $\displaystyle-2\pi(\not{p}+m)\delta(p^{2}-m^{2}){\rm
sign}(k^{0})\,F(p^{0}),$ $\displaystyle i{\cal S}^{--}(p)$
$\displaystyle=-\frac{i(\not{p}+m)}{p^{2}-m^{2}+i{\rm sign}(p^{0})\epsilon}$
$\displaystyle+2\pi(\not{p}+m){\rm
sign}(p^{0})\delta_{D}(p^{2}-m^{2})\,\left(1-F(p^{0})\right)$ $\displaystyle
i{\cal S}^{+-}(p)$ $\displaystyle=-2\pi\,{\rm
sgn}(p^{0})\,(\not{p}+m)\,\delta_{D}(p^{2}-m^{2})\,F(p^{0}),$ $\displaystyle
i{\cal S}^{-+}(p)$ $\displaystyle=2\pi\,{\rm
sgn}(p^{0})\,(\not{p}+m)\,\delta_{D}(p^{2}-m^{2})\left[1-F(p^{0})\right],$

For bosons the relevant propagators are

$\displaystyle i{\cal D}^{++}(p)$ $\displaystyle=\frac{i}{p^{2}-m^{2}+i{\rm
sign}(p^{0})\epsilon}$ $\displaystyle+2\pi\delta_{D}(p^{2}-m^{2})\,{\rm
sign}(p^{0})B(p^{0}),$ $\displaystyle i{\cal D}^{--}(p)$
$\displaystyle=-\frac{i}{p^{2}-m^{2}+i{\rm sign}(p^{0})\epsilon}$
$\displaystyle+2\pi\delta_{D}(p^{2}-m^{2})\,{\rm sign}(p^{0})(1+B(p^{0})),$

Here

$F(p^{0})=\frac{1}{e^{(p^{0}-\mu)/T}+1},\quad
B(p^{0})=\frac{1}{e^{(p^{0}-\mu)/T}-1}$ (80)

are the Fermi-Dirac and Bose-Einstein distributions respectively with chemical
potential $\mu$.

Let us start by computing diagram (i) of Fig. 3. As stated above, we are only
interested in processes $a+N+N\leftrightarrow N+N$ which correspond to the
explicit choice of internal Schwinger Keldysh indices shown. The other choices
of internal indices would correspond to cuts resulting in processes
$a\leftrightarrow N+N$ or $a+N\leftrightarrow N+\pi$ which are kinematically
forbidden for the low-energy axions considered here. The only choice of
indices giving rise to $a+N+N\leftrightarrow N+N$ processes corresponds to a
cut down the centre line. The corresponding expression for this contribution
to $\Pi^{>}$ is then given by

$\displaystyle i\,\Pi^{>}_{(i)}(q)=$
$\displaystyle\left(-\frac{ig_{an}}{2m}\right)^{2}\left(-\frac{i2mf}{m_{\pi}}\right)^{4}\int\frac{d^{4}p_{1}}{(2\pi)^{4}}\frac{d^{4}p_{2}}{(2\pi)^{4}}\frac{d^{4}p_{3}}{(2\pi)^{4}}\frac{d^{4}p_{4}}{(2\pi)^{4}}(2\pi)^{4}\delta^{(4)}(q+p_{1}+p_{2}-p_{3}-p_{4})i{\cal
D}^{--}(p_{4}-p_{2})i{\cal D}^{++}(p_{4}-p_{2}),$ $\displaystyle{\rm
Tr}\left[\gamma^{5}i\not{q}iS^{+-}(p_{1})\gamma^{5}i\not{q}i{\cal
S}^{--}(q+p_{1})\gamma^{5}i{\cal S}^{-+}(p_{3})\gamma^{5}i{\cal
S}^{++}(q+p_{1})\right]{\rm Tr}\left[\gamma^{5}i{\cal
S}^{+-}(p_{2})\gamma^{5}i{\cal S}^{-+}(p_{4})\right],$ (81)

where $q$ is the axion momentum. We can now insert the explicit expressions
(V)-(V) for the propagators. The Feynman and anti-Feynman propagators with
indices $++$ and $--$, respectively, have a vacuum (virtual) and thermal on-
shell part. The finite temperature piece can again be neglected on the grounds
it corresponds to thermal corrections to the pion and electron propagators,
which are again suppressed inside the neutron star since $T/m_{\pi},\,T/m\ll
1$. As a result, these propagators are dominated by their vacuum zero-
temperature piece, corresponding to the usual propagation of virtual
particles. This allows us to make connection with zero-temperature matrix
elements. With these approximations in mind, after taking the Dirac traces,
the remaining integral takes the form

$\displaystyle i\,\Pi^{>}_{(i)}(q)=$
$\displaystyle\int\frac{d^{4}p_{1}}{(2\pi)^{3}}\frac{d^{4}p_{2}}{(2\pi)^{3}}\frac{d^{4}p_{3}}{(2\pi)^{3}}\frac{d^{4}p_{4}}{(2\pi)^{3}}(2\pi)^{4}\delta^{(4)}(q+p_{1}+p_{2}-p_{3}-p_{4})$
$\displaystyle{\rm sgn}\,(p_{1}^{0})\delta(p_{1}^{2}-m^{2}){\rm
sgn}\,(p_{2}^{0})\delta(p_{2}^{2}-m^{2}){\rm
sgn}\,(p_{3}^{0})\delta(p_{3}^{2}-m^{2}){\rm
sgn}\,(p_{4}^{0})\delta(p_{4}^{2}-m^{2})$ $\displaystyle
F(p_{1}^{0})F(p_{2}^{0})(1-F(p_{3}^{0}))(1-F(p_{4}^{0}))\,C(q,p_{1},p_{2},p_{3},p_{4}),$
(82)

where $C$ is an effective matrix element which results from the trace parts
(to leading order in axion momentum $q_{\mu}$)

$\displaystyle C(q,p_{1},p_{2},p_{3},p_{4})$
$\displaystyle=\frac{4g_{an}^{2}f^{4}m^{2}}{m_{\pi}^{4}}\frac{1}{((p_{4}-p_{2})^{2}-m_{\pi}^{2})^{2}}\frac{16\left[4(q\cdot
p_{1})^{2}(m^{2}-p_{1}\cdot p_{3})+q^{2}m^{2}(6p_{1}\cdot
p_{3}-4m^{2})\right](m^{2}-p_{2}\cdot p_{4})}{((q+p_{1})^{2}-m^{2})^{2}},$
(83)

In what follows we will drop terms directly involving $q^{2}$ but retain those
arising from $(p_{i}\cdot q)^{2}$ in order to make contact with the estimates
made in e.g. Harris (2020); Brinkmann and Turner (1988). In those works, the
authors had in mind relativistic (i.e. lightlike) axions, and so they did not
retain terms $q^{2}$ in the matrix elements. In the traces, this is equivalent
to keeping terms such as $\textbf{p}_{i}.\textbf{p}_{j}(p.q)^{2}\sim
p_{F}^{2}m^{2}m_{a}^{2}$ whilst dropping terms $p^{4}q^{2}\sim
m^{4}m_{a}^{2}$. We repeat this step, allowing us to quote existing results
from the literature involving matrix elements and their interferences. For the
non-relativistic axions considered here, the latter terms would in fact
enhance the matrix element by $~{}(m/p_{F})^{2}$ so this approximation is
actually conservative. We leave these corrections for future work, with the
view that our main goal is to demonstrate a general procedure for computing
superradiant rates.

Next we must bring Eq. (V) into the form of a phase-space integral over on-
shell modes. This we do by expanding the delta functions using
$\delta(p^{2}-m^{2})=\frac{\delta(p^{0}-E_{p})}{2E_{p}}+\frac{\delta(p^{0}+E_{p})}{2E_{p}}$
where $E_{p}=\sqrt{\textbf{p}^{2}+m^{2}}$. We also use the relation
$F(-E)=1-{\bar{F}}(E)$, where ${\bar{F}}(E)$ is the phase space distribution
for the antiparticle given by replacing $\mu\rightarrow-\mu$ in Eq. (80).
Using these relations and performing the $p^{0}$ integrations we arrive at

$\displaystyle i\Pi^{>}_{(i)}(q)=\int$
$\displaystyle\frac{d^{3}p_{1}}{2E_{1}(2\pi)^{3}}\frac{d^{3}p_{2}}{2E_{2}(2\pi)^{3}}\frac{d^{3}p_{3}}{2E_{3}(2\pi)^{3}}\frac{d^{3}p_{4}}{2E_{4}(2\pi)^{3}}(2\pi)^{4}\delta^{(3)}({\bf
q}+{\bf p_{1}}+{\bf p_{2}}-{\bf p_{3}}-{\bf p_{4}})$ (84)
$\displaystyle[C(q,p_{1},p_{2},p_{3},p_{4})F(E_{1})F(E_{2})(1-F(E_{3}))(1-F(E_{4}))\delta(E_{q}+E_{1}+E_{2}-E_{3}-E_{4})+$
$\displaystyle
C(q,p_{1},\bar{p}_{2},\bar{p}_{3},p_{4})F(E_{1}){\bar{F}}(E_{3})(1-{\bar{F}}(E_{2}))(1-F(E_{4}))\delta(E_{q}+E_{1}+E_{3}-E_{2}-E_{4})+$
$\displaystyle
C(q,\bar{p}_{1},p_{2},\bar{p}_{3},p_{4}){\bar{F}}(E_{3})F(E_{2})(1-{\bar{F}}(E_{1}))(1-F(E_{4}))\delta(E_{q}+E_{3}+E_{2}-E_{1}-E_{4})+$
$\displaystyle
C(q,p_{1},\bar{p}_{2},p_{3},\bar{p}_{4})F(E_{1}){\bar{F}}(E_{4})(1-F(E_{2}))(1-{\bar{F}}(E_{3}))\delta(E_{q}+E_{1}+E_{4}-E_{2}-E_{3})+$
$\displaystyle
C(q,\bar{p}_{1},p_{2},p_{3},\bar{p}_{4})F(E_{2}){\bar{F}}(E_{4})(1-{\bar{F}}(E_{1}))(1-F(E_{3}))\delta(E_{q}+E_{2}+E_{4}-E_{1}-E_{3})+$
$\displaystyle
C(q,\bar{p}_{1},\bar{p}_{2},\bar{p}_{3},\bar{p}_{4}){\bar{F}}(E_{3}){\bar{F}}(E_{4})(1-{\bar{F}}(E_{1}))(1-{\bar{F}}(E_{2}))\delta(E_{q}+E_{3}+E_{4}-E_{1}-E_{2})],$

where we have used the fact that the energy conserving delta function can only
be satisfied if there are two neutrons each in the initial and final states
and $p,\bar{p}=(\pm E_{p},{\bf p})$ gives the on-shell 4-momenta of positive
and negative energies. Note we have also neglected the momentum transfer
associated to the axion in the $3$-momentum delta function.

Next we assume that we are in the regime $\mu/T\gg 1$, so that neutrons
dominate over anti-neutrons. Under these approximations, all those terms
containing a $\bar{F}$ in the initial state can be discarded, leading to

$\displaystyle i\Pi^{>}_{(i)}(q)=\int$
$\displaystyle\frac{d^{3}p_{1}}{2E_{1}(2\pi)^{3}}\frac{d^{3}p_{2}}{2E_{2}(2\pi)^{3}}\frac{d^{3}p_{3}}{2E_{3}(2\pi)^{3}}\frac{d^{3}p_{4}}{2E_{4}(2\pi)^{3}}$
$\displaystyle(2\pi)^{4}\delta^{(4)}(q+p_{1}+p_{2}-p_{3}-p_{4})F(E_{1})F(E_{2})(1-F(E_{3}))(1-F(E_{4}))C_{(i)}(q,p_{1},p_{2},p_{3},p_{4}).$
(85)

Thus we have brought the self-energy into the form of a standard collision
integral. This again represents the power of thermal and non-equilibrium field
theory not only to derive an effective equation of motion for the scalar
field, but also to express the damping coefficients in that equation in terms
of straightforward collision integrals.

To make contact with the calculations in Brinkmann and Turner (1988) we now go
to the non-relativistic limit and expand the matrix element to leading order
in $\textbf{p}_{i}/m$. In the limit of non-relativistic neutrons, we find, to
leading order in $q$:

$C(q,p_{1},p_{2},p_{3},p_{4})\simeq\frac{32f^{4}g_{an}^{2}m^{2}}{m_{\pi}^{4}}\frac{|{\bf
k}|^{4}}{(|{\bf k}|^{2}+m_{\pi}^{2})^{2}}.$ (87)

where $\textbf{k}=\textbf{p}_{2}-\textbf{p}_{4}$. Note this is precisely the
form of the matrix element associated to the modulus-squared of the top-left
diagram in Fig. 4, as can be seen from the appendix of ref. Brinkmann and
Turner (1988).

Figure 4: Matrix Elements for axion absorption corresponding to the cuts
through neutron lines in Fig. 3. The left and right columns correspond to
$t$\- and $u$-channel processes respectively.

Thus, as expected by the optical theorem, cuts through diagrams in fig. 3
correspond to $\Pi^{>}$ which feeds into the imaginary part of the retarded
self-energy. In turn these cuts can be interpreted in terms of S-matrix
elements for the process $a+N+N\rightarrow N+N$.

One can repeat this process for other diagrams in Fig. 3 which can be
similarly identified with matrix elements (and their interferences) in Fig. 4.
This is of course just the optical theorem, which demands that the sum over
diagrams $(i)$-$(iv)$ must give

$\displaystyle i\Pi^{>}(q)=\int$
$\displaystyle\frac{d^{3}p_{1}}{2E_{1}(2\pi)^{3}}\frac{d^{3}p_{2}}{2E_{2}(2\pi)^{3}}\frac{d^{3}p_{3}}{2E_{3}(2\pi)^{3}}\frac{d^{3}p_{4}}{2E_{4}(2\pi)^{3}}$
$\displaystyle(2\pi)^{4}\delta^{(4)}(q+p_{1}+p_{2}-p_{3}-p_{4})F(E_{1})F(E_{2})(1-F(E_{3}))(1-F(E_{4}))\sum_{\text{spins}}|\mathcal{M}|^{2}$
(88)

where $\mathcal{M}$ is the full matrix element associated to the diagrams in
Fig. 4. This is a classic result which has already been computed a long time
ago Brinkmann and Turner (1988); Iwamoto (1984)

$\displaystyle\sum_{\text{spins}}|\mathcal{M}|^{2}=\frac{256}{3}\frac{f^{4}m^{2}g_{an}^{2}}{m_{\pi}^{4}}\left[\frac{\mathbf{k}^{4}}{\left(\mathbf{k}^{2}+m_{\pi}^{2}\right)^{2}}+\frac{\mathbf{l}^{4}}{\left(\mathbf{l}^{2}+m_{\pi}^{2}\right)^{2}}+\frac{\mathbf{k}^{2}\mathbf{l}^{2}-3\left(\mathbf{k}\cdot\mathbf{l}\right)^{2}}{\left(\mathbf{k}^{2}+m_{\pi}^{2}\right)\left(\mathbf{l}^{2}+m_{\pi}^{2}\right)}\right].$
(89)

Note that strictly speaking, as can be seen in Brinkmann and Turner (1988),
this is actually the matrix element given by averaging over the direction of
the axion three-momentum, so that in the matrix element, one replaces e.g.
$(\hat{\textbf{k}}\cdot\hat{\textbf{q}})^{2}\rightarrow\braket{(\hat{\textbf{k}}\cdot\hat{\textbf{q}})^{2}}=1/3$.
Here hats denote unit vectors. However, by symmetry we know the final self
energy $\Pi^{>}(q)=\Pi^{>}(q^{0},\textbf{q})$ must be isotropic with respect
to q so that
$\Pi^{>}(q^{0},\textbf{q})=\Pi^{>}(q^{0},\left|\textbf{q}\right|)$. It
therefore follows that the self-energy is invariant under the averaging over
the direction of $\hat{\textbf{q}}$.

Formally this means $\Pi^{>}=\braket{\Pi^{>}}$. Furthermore, when applying the
operator $\braket{\cdots}$ to Eq. (V), the averaging can be taken past the
distribution functions and energy conserving delta functions (provided we
neglect the axion 3-momentum in the spatial 3-momentum conserving delta
function) and directly onto the matrix element in the collision integral. As a
result, to leading order in axion-momentum, the answer is unchanged by using a
matrix element averaged over the direction of axion momentum.

Having related the self-to a collision integral (V), whose matrix element (89)
is known, all that remains is to actually compute this integral. Fortunately,
this integral can be directly identified with the so-called axion mean free
path $\lambda^{-1}_{>}$ defined by Ishizuka and Yoshimura (1990); Burrows  _et
al._ (1990); Iwamoto (2001); Dietrich and Clough (2019); Harris (2020)

$\lambda^{-1}_{>}=\frac{i\Pi^{>}(q)}{2q^{0}}.$ (90)

The identification of a $>$ subscript with the mean free path will become
clearer shortly. For now though, we simply note that by looking up the pre-
existing calculations of the mean-free path $\lambda^{-1}_{>}$ in the
literature, we can immediately read off $\Pi^{>}$. The damping rate then
follows immediately from Eq. (27).

In the degenerate limit appropriate for neutron stars, this integral has in
fact already been computed in a Fermi Surface approximation Ishizuka and
Yoshimura (1990); Iwamoto (2001); Harris (2020), in that limit one finds

$\lambda^{-1}_{>}=\frac{1}{18\pi^{5}}\frac{f^{4}g_{an}^{2}m^{2}}{4m_{\pi}^{4}}p_{Fn}Q(y)\frac{\omega^{2}+4\pi^{2}T^{2}}{1-e^{-\omega/T}},$
(91)

where

$\displaystyle Q(y)=$ $\displaystyle
4-\frac{1}{1+y^{2}}+\frac{2y^{2}}{\sqrt{1+2y^{2}}}\arctan{\left(\frac{1}{\sqrt{1+2y^{2}}}\right)}$
$\displaystyle-5y\arcsin{\left(\frac{1}{\sqrt{1+y^{2}}}\right)},$ (92)

and $y=m_{\pi}/(2p_{F})$ with $p_{F}$ the neutron Fermi momentum. We can
combine Eq. (27) with the KMS relation $\Pi^{>}(q)=e^{\omega/T}\Pi^{<}(q)$ to
obtain

$\displaystyle{\rm Im}[\Pi^{R}(p,x)]$
$\displaystyle=-\frac{i}{2}\left(\Pi^{>}(p,x)-\Pi^{<}(p,x)\right)$
$\displaystyle=-\frac{i}{2}\left(1-e^{-\omega/T}\right)\Pi^{>}(p,x).$ (93)

Inserting (90) and (91) into this we can immediately compute the anti-
hermitian part of the self energy in the degenerate limit, which gives

$\text{Im}\,\Pi^{R}={\color[rgb]{.75,0,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.75,0,.25}-}\frac{\omega}{18\pi^{5}}\frac{f^{4}g_{an}^{2}m^{2}}{4m_{\pi}^{4}}p_{F}Q(y)\left(\omega^{2}+4\pi^{2}T^{2}\right)$
(94)

Putting this together, from Eq. (4), we are able to get a superradiance
expressed in terms of microphysical parameters:

$\Gamma^{a}_{n\ell
m}=\frac{1}{18\pi^{5}}\frac{f^{4}g_{an}^{2}m^{2}}{4m_{\pi}^{4}}p_{F}Q(y)4\pi^{2}T^{2}\left(\frac{R}{r_{n\ell}}\right)^{2\ell+3}\frac{(2\ell+1)\text{!!}(2\ell+n+1)!}{(\ell+n+1)n!(2\ell+1)!^{2}(2\ell+3)!!}\frac{(m\Omega-\omega_{\ell
n})}{\omega_{\ell n}},$ (95)

Let us make some estimates of this rate. We can take typical neutron star
values $\Omega=600{\rm Hz}$, $M=1.4M_{\odot}$ and $R=10{\rm km}$ and
$p_{F}\simeq 300{\rm MeV}$. For couplings constants we can take $f\simeq 1$
and the maximum limit set by Supernova constraints is $g_{an}\simeq 10^{-9}$.
Furthermore, for hot young stars, the rate is at its highest, here (although
briefly) temperatures could reach up to $T=10^{11}{\rm K}$ Raffelt (1996)
giving an upper bound on the superradiance rate. For these numbers, the
$\ell=m=1,n=0$ has a superradiance timescale $\Gamma_{011}^{-1}\simeq
10^{13}{\rm yr}$, somewhat longer than the age of the Universe. For older,
cooler neutron stars, the rate would of course be lower. Hence at least for
the simplest axion scenarios which couple axions to nuclear matter, the star
is protected from superradiant instabilities over timescales of the Universe.

At this point a few remarks are in order concerning the microphysical damping
rate as embodied in our Eq. (94) and the notion of the mean-free path found
elsewhere in the literature Ishizuka and Yoshimura (1990); Burrows  _et al._
(1990); Iwamoto (2001); Dietrich and Clough (2019); Harris (2020). Formally
the damping rate of the field in the equation of motion is defined as the
difference of two Wightman self energies

$\text{Im}\,\Pi^{R}=\frac{1}{2i}\left(\Pi^{>}-\Pi^{<}\right)$ (96)

where, by virtue of the KMS relation, we can write the self-energies more
transparently as

$i\Pi^{>}=-(1+B)\text{Im}\Pi^{R},\qquad i\Pi^{<}=-B\,\text{Im}\Pi^{R},$ (97)

where $B$ is the Bose-Einstein distribution function defined in (80). Notice
therefore that the true mean free path involves two opposite contributions
where $\Pi^{>}$ and $\Pi^{<}$ are to be interpreted as gain and loss
processes, respectively, as is apparent from the presence of the $B$ and $1+B$
prefactors. Notice however, that the collision integral in Harris  _et al._
(2020) which the authors claim corresponds to the mean free path, actually
only give the piece $\Pi^{>}$, explaining the appearance of a prefactor
$1/(1-e^{-\omega/T})=1+B$. In reality, the authors should have taken both
contributions, which would have given a mean free path
$\lambda^{-1}=(\text{Im}\Pi^{R}/\omega)$.

The fact that it is ${\rm Im}\Pi$ and not just one of $\Pi^{>}$ or $\Pi^{<}$
which should be identified with the mean free path is also apparent if one
write the spectral propagator of the axion $G_{s}(q)$ in Breit-Wigner Form
Blaizot and Iancu (2002)

$G_{s}(q)=\frac{2\text{Im}\Pi^{R}}{\left((q^{0})^{2}-E_{\textbf{q}}-\text{Re}\,\Pi^{R}\right)^{2}+(\text{Im}\Pi^{R})^{2}}.$
(98)

This spectral propagator describes the propagating axion modes in the medium.
In the lossless limit $\text{Im}\Pi^{R}\rightarrow 0$, we recognise the above
form as the Lorenz representation of the delta function so that the spectral
function $G_{s}(q)$ simply becomes a mass-shell condition
$G_{a}(q)\propto\delta\left((q^{0})^{2}-E_{\textbf{q}}-\text{Re}\,\Pi^{R}\right)$
of a zero width mode where $\text{Re}\,\Pi^{R}$ encodes the so-called
dispersive part of $\Pi^{R}$ which gives radiative corrections to the axion
dispersion relation. This provides yet another justification that
$\text{Im}\,\Pi^{R}$ sets the width of axion states associated to the lifetime
of the axion in-medium.

The question is then whether the authors of Harris  _et al._ (2020)
misidentifying the mean free path with the wrong type of self-energy has any
quantitative bearing on their calculation. Since there the authors typically
consider axions with energies $\omega\gtrsim T$, we have that $B\lesssim 1$,
so that proved one considers axion frequencies at or above the temperature,
$B$ factors become negligible and one finds
$\text{Im}\Pi^{R}\simeq\Pi^{>}/2i$. This misidentification would therefore
only change the mean free path by $\mathcal{O}(1)$ so long as one sticks to
$\omega\gtrsim T$.

In fact, the contribution to the damping rate associated to $\Pi^{>}$
corresponds to stimulated emission of axions. This means that for
superradiance, since we are in the regime $\omega\ll T$ we have $B>>1$ owing
to large bose enhancement. This means that the two Wightman contributions
almost cancel, with their difference being suppressed by $\omega/T$ so that to
leading order in this ratio we have
$\text{Im}\Pi^{R}\simeq\frac{\omega}{T}\Pi^{>}$. Hence in our case, since
$\omega/T\simeq 10^{-17}$ for $T=10^{8}$K and $\omega=10^{-13}$eV, this would
lead to a radically different estimate for the superradiance rate if we had
misidentified the wrong collision integral (associated to $\Pi^{>}$) with the
true damping rate of the field. This again emphasises the importance of
carrying out a rigorous first principles calculation of the damping rate
within non-equilibrium field theory and not simply guessing it from a
collision integral.

##  VI Discussion and Conclusions

In this paper we have presented the first general treatment of superradiance
in stars, which through straightforward adaptation allows the reader to
compute the superradiant instability for any fundamental interaction.

One key development was to employ the full power of field theory at finite
density. As shown in Section III, this allowed us to understand how the
damping of fields should be identified with the self energy as embodied in Eq.
(32) and what precise approximations are necessary to produce a simple damping
term. Crucially, we were then able to relate the self-energy (and thence this
damping coefficient) to microphysical scattering processes, as illustrated in
the case of axions in Section V.

By contrast, in previous work a damping term was added in a somewhat schematic
way Cardoso  _et al._ (2015) and it remained opaque as to how to relate this
to microphysics. It is worth noting that Cardoso  _et al._ (2017) made some
attempt to interpret the damping coefficients for gauge field in terms of QED-
like processes but it remained unclear how these damping terms emerge in a
more general theories and indeed how they can be justified from a first
principles approach in QFT.

We are now able to justify and explain the existence of damping terms for
fields, and explicitly calculate them in terms of microphysical absorption
processes. In the case of axions, this corresponded to inverse Bremsstrahlung
processes $a+N+N\rightarrow N+N$. The final expression for the damping
coefficient takes the form of a collision integral, which for the present
interactions had already been derived elsewhere in the literature when
computing the mean free path of axions in the context of energy transport in
neutron stars.

Having computed the microphysical damping coefficient for the scalar field, it
was then straightforward to compute the absorption rate of long-wavelength
modes into the whole star. Our next key insight was that this was sufficient
to compute the superradiant instability rates for massive fields in
gravitational bound states of a rotating star. This we achieved through
appealing to an EFT of long-wavelength modes Endlich and Penco (2017). This
approach allowed us to compute the absorption of arbitrary multiple modes for
a static, non-rotating star, and then through symmetry arguments, relate this
through a simple algebraic substitution for the rate of a rotating star. This
completely circumvents the need to introduce bulk currents for the rotating
matter in the star when computing the effective equation of motion (32) as was
done in Cardoso  _et al._ (2017). Thus rather elegantly, one simply needs to
compute the absorption of a massless mode into a non-rotating star, and the
superradiance rate can be derived in a quick and straightforward way.

In summary, we now have a fully working pipeline to compute superradiance
rates in stars for any interaction. Having presented this general treatment,
one can now of course generalise this to a broader class of couplings and
fields to see for which fundamental interactions the rate is highest, and
investigate their corresponding signatures.

Acknowledgements

JIM is supported by an FSR Incoming Postdoctoral Fellowship. FCD is supported
by Stephen Hawking Fellowship EP/T01668X/1 and STFC grant ST/T001011/1. BG
acknowledges support from the Collaborative Research Centre SFB 1258 of the
Deutsche Forschungsgemeinschaft. We thank M. Drewes for useful conversations.
We are indebted to Paolo Pani, Richard Britto and Vitor Cardoso for comments
on the manuscript and to Steven Harris for his input on axion mean free paths.

Data Access Statement

No data sets are used or created in this paper.

##  Appendix A Convolutions of Green Functions

The full derivation of the result (25) is as follows. We begin with

$\displaystyle[\Pi^{R}\odot\phi](x)=\int d^{4}y\Pi^{R}(x,y)\phi(y)$
$\displaystyle=$ $\displaystyle\int
d^{4}y\int\frac{d^{4}p}{(2\pi)^{4}}\Pi^{R}\left(p,\frac{x+y}{2}\right)\phi(y)e^{-ip\cdot(x-y)}.$
(99)

where we defined a convolution operator $\odot$ to compactify notation.
Expanding $\Pi^{R}$ about $p=0$ to all orders, we obtain

$\displaystyle[\Pi^{R}\odot\phi](x)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\int
d^{4}y\int\frac{d^{4}p}{(2\pi)^{4}}\frac{p^{n}}{n!}\left[\partial^{n}_{p}\Pi^{R}\left(p,\frac{x+y}{2}\right)\right]_{p=0}\phi(y)$
$\displaystyle\hskip 56.9055pt\times e^{-ip\cdot(x-y)}.$ (100)

This can be written as

$\displaystyle[\Pi^{R}\odot\phi](x)$ $\displaystyle=$ $\displaystyle\int
d^{4}y\int\frac{d^{4}p}{(2\pi)^{4}}\left[\Pi^{R}\left(p,\frac{x+y}{2}\right)\phi(y)e^{-i\overleftarrow{\partial}_{p}\cdot\overrightarrow{\partial}_{y}}\right]_{p=0}$
$\displaystyle\hskip 56.9055pt\times e^{-ip\cdot(x-y)}.$ (101)

Performing the $p$ integration gives a delta function $\delta(x-y)$. One can
then integrate by parts with respect to $y$ to arrive at

$\displaystyle[\Pi^{R}\odot\phi](x)$ $\displaystyle=$ $\displaystyle\int
d^{4}y\delta(x-y)\left[e^{i\partial_{p}\cdot\partial_{y}}\Pi^{R}\left(p,\frac{x+y}{2}\right)\phi(y)\right]_{p=0}$
$\displaystyle=$
$\displaystyle\left[e^{i\partial_{p}\cdot\partial_{y}}\Pi^{R}\left(p,\frac{x+y}{2}\right)\phi(y)\right]_{p=0,y=x},$
(102)

which is the result (25).

##  Appendix B Normalisation of angular momentum eigenstates

In this appendix we will derive the normalisation of angular momentum
eigenstates given in equation (50). As usual in quantum field theory, delta
functions arising from the normalisation of states correspond to the infinite
volume of spacetime. Our free angular momentum eigenstates are normalised as

$\braket{\omega\ell m\left missing}{\omega\ell
m}=2\pi\lim_{\omega^{\prime}\rightarrow\omega}\delta(\omega-\omega^{\prime})=\frac{2\pi}{v}\lim_{k^{\prime}\rightarrow
k}\delta(k-k^{\prime})$ (103)

This arises from the normalisation spherical bessel functions associated to
spherical waves:

$\int_{0}^{\infty}dr\,r^{2}\,j_{\ell}(kr)j_{\ell}(k^{\prime}r)=\frac{\pi}{2k^{2}}\delta(k-k^{\prime}),$
(104)

which can be regulated using a finite integration range up to radius
$\mathcal{R}$

$\displaystyle\int_{0}^{R}dr\,r^{2}\,j_{\ell}(kr)j_{\ell}(kr)$ (105)
$\displaystyle=\frac{\pi\mathcal{R}^{2}\left(J_{l+\frac{1}{2}}(k\mathcal{R}){}^{2}-J_{\ell-\frac{1}{2}}(k\mathcal{R})J_{\ell+\frac{3}{2}}(k\mathcal{R})\right)}{4k}.$

Hence we can regulate

$\displaystyle\lim_{k^{\prime}\rightarrow k}\delta(k-k^{\prime})$
$\displaystyle=\lim_{R\rightarrow\infty}\frac{2k^{2}}{\pi}\frac{\pi\mathcal{R}^{2}\left(J_{l+\frac{1}{2}}(k\mathcal{R}){}^{2}-J_{l-\frac{1}{2}}(k\mathcal{R})J_{l+\frac{3}{2}}(kR)\right)}{4k}.$
(106)

We now use the asymptotic form of the Bessel function

$J_{\alpha}(z)\rightarrow\sqrt{\frac{2}{\pi z}}{\rm
cos}\left(z-\frac{\alpha\pi}{2}-\frac{\pi}{4}\right),$ (107)

for large $z$ and half-integer $\alpha$. This gives

$\displaystyle\delta(k-k^{\prime})$
$\displaystyle=\lim_{\mathcal{R}\rightarrow\infty}\frac{2k^{2}}{\pi}\frac{\mathcal{R}}{2k^{2}}\left({\rm
sin}^{2}\left(k\mathcal{R}-\frac{\ell\pi}{2}\right)+{\rm
cos}^{2}\left(k\mathcal{R}-\frac{\ell\pi}{2}\right)\right)$
$\displaystyle=\lim_{\mathcal{R}\rightarrow\infty}\frac{\mathcal{R}}{\pi}$
(108)

Note that for a spherical wave with energy density $\rho\propto 1/r^{2}$, the
number of particles per unit radius (i.e. number of particles in a spherical
shell of unit width) is constant. We can therefore interpret the factor
$1/\mathcal{R}$ as the number of particles per unit radius. This is the
spherical equivalent of the usual replacement $(2\pi)^{3}\delta^{(3)}({\bf
k}=0)\rightarrow V$.

##  Appendix C Superradiant bound states

In this appendix we detail the spherical integrals used to derive equation
(65) in section IV.2. Equation (IV.2) gives explicitly

$\displaystyle P_{\rm abs}\simeq\int^{\mathcal{T}}_{0}$ $\displaystyle
dtdt^{\prime}\braket{O_{J_{1}...J_{\ell}}(t^{\prime})O_{L_{1}...L_{\ell}}(t)}$
$\displaystyle\bra{n,l,m}\partial^{I_{1}}...\partial^{I_{\ell}}\phi(t^{\prime})\ket{0}$
$\displaystyle\bra{0}\partial^{K_{1}}...\partial^{K_{\ell}}\phi(t)\ket{n,l,m}$
$\displaystyle
R_{I_{1}}^{J_{\ell}}(t^{\prime})...R_{I_{\ell}}^{J_{\ell}}(t^{\prime})R_{K_{1}}^{L_{\ell}}(t)...R_{K_{l}}^{L_{\ell}}(t).$
(109)

We identify the Wightman correlation function:

$\braket{O_{J_{1}...J_{\ell}}(t^{\prime})O_{L_{1}...L_{\ell}}(t)}=\delta_{J_{1}...J_{\ell}}^{L_{1}...L_{\ell}}\int\frac{d\omega}{2\pi}\Delta_{\ell}(\omega){\rm
e}^{i\omega(t-t^{\prime})}.$ (110)

Inserting this into $P_{\rm abs}$ gives:

$\displaystyle P_{\rm
abs}=\int\frac{d\omega^{\prime}}{2\pi}\Delta_{\ell}(\omega^{\prime})$ (111)
$\displaystyle\left|\int^{\mathcal{T}}_{0}dt{\rm
e}^{i\omega^{\prime}t}\bra{0}\partial^{K_{1}}...\partial^{K_{\ell}}\phi(t)\ket{n\ell
m}R^{L_{1}}_{K_{1}}(t)...R^{L_{l}}_{K_{\ell}}(t)\right|^{2},$

where we use a compact notation meaning that the free $K_{i}$ indices become
contracted with themselves upon taking the modulus squared. Using the mode
expansion Eq. (IV.2), we can write explicitly

$\displaystyle\bra{0}\partial^{K_{1}}...\partial^{K_{\ell}}\phi(t)\ket{n\ell
m}=\sum_{n^{\prime}l^{\prime}m^{\prime}}\frac{1}{\sqrt{2\omega_{n^{\prime}l^{\prime}}}}$
$\displaystyle\times\bra{0}\hat{a}_{n^{\prime}\ell^{\prime}m^{\prime}}\partial^{K_{1}}...\partial^{K_{l}}f_{n^{\prime}l^{\prime}m^{\prime}}(0){\rm
e}^{-i\omega_{n^{\prime}\ell^{\prime}}}\ket{n\ell m}.$

$f_{n\ell m}$ is evaluated at the spatial location of the star. Using the
orthogonality of the states $\ket{n\ell m}$, we have:

$\bra{0}\partial^{K_{1}}...\partial^{K_{\ell}}\phi(t)\ket{n\ell
m}=\frac{1}{\sqrt{2\omega_{n\ell}}}\partial^{K_{1}}...\partial^{K_{\ell}}f_{n\ell
m}(0){\rm e}^{-i\omega_{n\ell}t}.$ (112)

This gives

$\displaystyle P_{\rm
abs}=\int\frac{d\omega}{2\pi}\frac{\Delta_{\ell}(\omega)}{2\omega_{\ell n}}$
$\displaystyle\times\left|\int^{\mathcal{T}}_{0}dt{\rm
e}^{i(\omega-\omega_{n\ell})t}\partial^{I_{1}}...\partial^{I_{\ell}}f_{n\ell
m}(r=0)R^{J_{1}}_{I_{1}}(t)...R^{J_{\ell}}_{I_{\ell}}(t)\right|^{2}.$ (113)

We will now use the identities Thorne (1980)

$Y_{lm}(\theta,\phi)=\sqrt{\frac{(2\ell+1)!!}{4\pi\ell!}}V^{m}_{I_{1}...I_{\ell}}\hat{r}^{I_{1}}...\hat{r}^{I_{\ell}},$
(114)

where

$\displaystyle V_{I_{1}\cdots I_{\ell}}^{m}=\sum_{j=0}^{[(\ell-m)/2]}\nu^{\ell
mj}\left(\delta_{\left(I_{1}\right.}^{1}+i\delta_{\left(I_{1}\right.}^{2}\right)\cdots\Big{(}\delta_{I_{m}}^{1}+i\delta_{I_{m}}^{2}\Big{)}$
$\displaystyle\times\delta_{I_{m+1}}^{3}\cdots\delta_{I_{\ell-2j}}^{3}\left(\delta_{I_{\ell-2j+1}}^{a_{1}}\delta_{I_{\ell-2j+2}}^{a_{1}}\right)\cdots\left(\delta_{I_{\ell-1}}^{a_{j}}\delta_{\left.I_{\ell}\right)}^{a_{j}}\right),$
(115)

with

$\nu^{\ell
mj}=\sqrt{\frac{\ell!(\ell-m)!}{(2\ell-1)!!(\ell+m)!}}\frac{(-1)^{m+j}(2\ell-2j)!}{2^{\ell}j!(\ell-j)!(\ell-m-2j)!}.$
(116)

In order to progress, we must evaluate the multi-derivative operator acting on
$f_{nlm}$ and extarct its leading order contribution. For this we shall
require an expression for the bound states Day and McDonald (2019)

$\displaystyle R_{n\ell}=$ $\displaystyle\mu^{3/2}\alpha_{\ell
n}^{3/2}\sqrt{\frac{n!}{2(n+2\ell+1)!(n+\ell+1)}}e^{-x/2}x^{\ell}L_{n}^{2\ell+1}\left[x\right],$
(117)

where

$\displaystyle x=r\mu\alpha_{\ell n},$ $\displaystyle\alpha_{\ell n}=\frac{\mu
r_{s}}{\ell+n+1},\quad\omega^{2}_{\ell
n}=\mu^{2}\left(1-\frac{\alpha^{2}_{\ell n}}{4}\right)\simeq\mu^{2}$ (118)

and $L_{n}^{2\ell+1}(x)$ is a generalised Laguerre polynomial. Note that these
satisfy the correct normalisation $\int drr^{2}R_{n\ell}^{2}=1$. From this we
can combine Eqs. (114) and (117) to derive the leading order result in
$\alpha_{\ell n}$ of the form

$\displaystyle\partial^{I_{1}}...\partial^{I_{\ell}}f_{n\ell m}(r=0)\simeq$
$\displaystyle\sqrt{\frac{l!(2\ell+1)!!}{4\pi}}V^{I_{1}...I_{\ell}}_{m}\left.\frac{R_{n\ell}(r)}{r^{\ell}}\right|_{r=0},$
(119)

where we used

$\displaystyle
V^{m}_{J_{1}...J_{\ell}}\partial^{I_{1}}...\partial^{I_{\ell}}r^{J_{1}}...r^{J_{\ell}}$
$\displaystyle=l!V^{m}_{I_{1}...I_{\ell}}.$ (120)

Next, making use of the identities

$\displaystyle V_{I_{1}\cdots I_{\ell}}^{m}R^{I_{1}}{}_{J_{1}}\cdots
R^{I_{\ell}}{}_{J_{\ell}}=V_{J_{1}\cdots J_{\ell}}^{m}e^{im\Omega t},$ (121)
$\displaystyle(V^{m}_{J_{1}\cdots J_{\ell}})^{*}V_{m^{\prime}}^{J_{1}\cdots
J_{\ell}}=\delta^{m}_{m^{\prime}}$ (122)

which leads to

$\displaystyle P_{\rm abs}=$
$\displaystyle\frac{l!(2\ell+1)!!}{4\pi}\left|\frac{R_{n\ell}(r)}{r^{\ell}}\right|^{2}_{r=0}$
(123)
$\displaystyle\times\int\frac{d\omega}{2\pi}\frac{\Delta_{\ell}(\omega)}{2\omega_{\ell
n}}\left[\int^{\mathcal{T}}_{0}dt{\rm e}^{i(m\Omega+\omega-\mu)t}\right]^{2}.$

Explicitly we have

$\displaystyle\left.\frac{R_{n\ell}(r)}{r^{\ell}}\right|_{r=0}=$
$\displaystyle\mu^{\ell+3/2}\alpha_{\ell
n}^{3/2+\ell}\sqrt{\frac{n!}{2(n+2\ell+1)!(n+\ell+1)}}$
$\displaystyle\times\frac{(2\ell+n+1)!}{(1+2\ell)!\,n!}$ (124)

where we used $L_{n}^{2\ell+1}[0]=(1+2l+n)!/[(1+2l)!\,n!]$. Putting this
together, for sufficiently long times we can therefore write

$\displaystyle P_{\rm abs}$ $\displaystyle=A_{n\ell
m}\left(\frac{1}{r_{n\ell}}\right)^{2\ell+3}\int\frac{d\omega}{2\pi}\frac{\Delta_{\ell}(\omega)}{2\omega_{\ell
n}}\mathcal{T}(2\pi)\delta(m\Omega+\omega-\mu),$ $\displaystyle=A_{n\ell
m}\left(\frac{1}{r_{n\ell}}\right)^{2\ell+3}\frac{\Delta_{\ell}(\omega)}{2\omega_{\ell
n}}\mathcal{T},$ (125)

where

$r_{n\ell}=\frac{n+l+1}{\mu^{2}r_{s}}$ (126)

and

$\displaystyle A_{n\ell m}$
$\displaystyle=\frac{\ell!(2\ell+1)\text{!!}(2\ell+n+1)!}{8\pi(\ell+n+1)n!(2\ell+1)!^{2}}.$
(127)

## References

  * Berti  _et al._ (2022) E. Berti _et al._ , in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.07984 [hep-ph] .
  * Brito  _et al._ (2015) R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, pp.1 (2015), arXiv:1501.06570 [gr-qc] .
  * Bekenstein and Schiffer (1998) J. D. Bekenstein and M. Schiffer, Phys. Rev. D 58, 064014 (1998), arXiv:gr-qc/9803033 .
  * Zel’dovich (1971) Y. B. Zel’dovich, Pis’ma Zh. Eksp. Teor. Fiz. 14, 270 [JETP Lett. 14, 180 (1971)] (1971).
  * Zel’dovich (1972) Y. B. Zel’dovich, Zh. Eksp. Teor. Fiz 62, 2076 [Sov.Phys. JETP 35, 1085 (1972)] (1972).
  * Vicente  _et al._ (2018) R. Vicente, V. Cardoso, and J. C. Lopes, Phys. Rev. D 97, 084032 (2018), arXiv:1803.08060 [gr-qc] .
  * Eskin (2016) G. Eskin, Rev. Math. Phys. 28, 1650025 (2016), arXiv:1509.03197 [math-ph] .
  * Richartz  _et al._ (2009) M. Richartz, S. Weinfurtner, A. J. Penner, and W. G. Unruh, Phys. Rev. D 80, 124016 (2009), arXiv:0909.2317 [gr-qc] .
  * Arvanitaki  _et al._ (2015) A. Arvanitaki, M. Baryakhtar, and X. Huang, Phys. Rev. D 91, 084011 (2015), arXiv:1411.2263 [hep-ph] .
  * Brito  _et al._ (2017a) R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein, and P. Pani, Phys. Rev. Lett. 119, 131101 (2017a), arXiv:1706.05097 [gr-qc] .
  * Baryakhtar  _et al._ (2017) M. Baryakhtar, R. Lasenby, and M. Teo, Phys. Rev. D 96, 035019 (2017), arXiv:1704.05081 [hep-ph] .
  * Brito  _et al._ (2017b) R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein, and P. Pani, Phys. Rev. D 96, 064050 (2017b), arXiv:1706.06311 [gr-qc] .
  * Chen  _et al._ (2020) Y. Chen, J. Shu, X. Xue, Q. Yuan, and Y. Zhao, Phys. Rev. Lett. 124, 061102 (2020), arXiv:1905.02213 [hep-ph] .
  * Baumann  _et al._ (2019) D. Baumann, H. S. Chia, and R. A. Porto, Phys. Rev. D 99, 044001 (2019), arXiv:1804.03208 [gr-qc] .
  * Arvanitaki and Dubovsky (2011) A. Arvanitaki and S. Dubovsky, Phys. Rev. D 83, 044026 (2011), arXiv:1004.3558 [hep-th] .
  * Pani  _et al._ (2012) P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, Phys. Rev. Lett. 109, 131102 (2012), arXiv:1209.0465 [gr-qc] .
  * Cardoso  _et al._ (2018) V. Cardoso, O. J. C. Dias, G. S. Hartnett, M. Middleton, P. Pani, and J. E. Santos, JCAP 03, 043 (2018), arXiv:1801.01420 [gr-qc] .
  * Blas and Witte (2020) D. Blas and S. J. Witte, Phys. Rev. D 102, 123018 (2020), arXiv:2009.10075 [hep-ph] .
  * Wang  _et al._ (2022) Z. Wang, T. Helfer, K. Clough, and E. Berti, Phys. Rev. D 105, 104055 (2022), arXiv:2201.08305 [gr-qc] .
  * Conlon and Herdeiro (2018) J. P. Conlon and C. A. R. Herdeiro, Phys. Lett. B 780, 169 (2018), arXiv:1701.02034 [astro-ph.HE] .
  * Cardoso  _et al._ (2015) V. Cardoso, R. Brito, and J. L. Rosa, Phys. Rev. D 91, 124026 (2015), arXiv:1505.05509 [gr-qc] .
  * Cardoso  _et al._ (2017) V. Cardoso, P. Pani, and T.-T. Yu, Phys. Rev. D 95, 124056 (2017), arXiv:1704.06151 [gr-qc] .
  * Day and McDonald (2019) F. V. Day and J. I. McDonald, JCAP 10, 051 (2019), arXiv:1904.08341 [hep-ph] .
  * Endlich and Penco (2017) S. Endlich and R. Penco, JHEP 05, 052 (2017), arXiv:1609.06723 [hep-th] .
  * Keldysh (1964) L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964).
  * Schwinger (1961) J. S. Schwinger, J. Math. Phys. 2, 407 (1961).
  * Berges (2004) J. Berges, AIP Conf. Proc. 739, 3 (2004), arXiv:hep-ph/0409233 .
  * Calzetta and Hu (2008) E. A. Calzetta and B.-L. B. Hu,  _Nonequilibrium Quantum Field Theory_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2008).
  * Calzetta and Hu (1988) E. Calzetta and B. L. Hu, Phys. Rev. D 37, 2878 (1988).
  * Prokopec  _et al._ (2004a) T. Prokopec, M. G. Schmidt, and S. Weinstock, Annals Phys. 314, 208 (2004a), arXiv:hep-ph/0312110 [hep-ph] .
  * Ai  _et al._ (2021a) W.-Y. Ai, M. Drewes, D. Glavan, and J. Hajer, JHEP 11, 160 (2021a), arXiv:2108.00254 [hep-ph] .
  * Cheung  _et al._ (2015) Y.-K. E. Cheung, M. Drewes, J. U. Kang, and J. C. Kim, JHEP 08, 059 (2015), arXiv:1504.04444 [hep-ph] .
  * Cornwall  _et al._ (1974) J. M. Cornwall, R. Jackiw, and E. Tomboulis, Phys. Rev. D 10, 2428 (1974).
  * Yoshino and Kodama (2012) H. Yoshino and H. Kodama, Prog. Theor. Phys. 128, 153 (2012), arXiv:1203.5070 [gr-qc] .
  * Yoshino and Kodama (2014) H. Yoshino and H. Kodama, PTEP 2014, 043E02 (2014), arXiv:1312.2326 [gr-qc] .
  * Delacrétaz  _et al._ (2014) L. V. Delacrétaz, S. Endlich, A. Monin, R. Penco, and F. Riva, JHEP 11, 008 (2014), arXiv:1405.7384 [hep-th] .
  * Kim (1979) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ (1980) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 166, 493 (1980).
  * Dine  _et al._ (1981) M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
  * Zhitnitsky (1980) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980).
  * Ai  _et al._ (2021b) W.-Y. Ai, J. S. Cruz, B. Garbrecht, and C. Tamarit, Phys. Lett. B 822, 136616 (2021b), arXiv:2001.07152 [hep-th] .
  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Svrcek and Witten (2006) P. Svrcek and E. Witten, JHEP 06, 051 (2006), arXiv:hep-th/0605206 .
  * Arvanitaki  _et al._ (2010) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Pshirkov and Popov (2009) M. S. Pshirkov and S. B. Popov, J. Exp. Theor. Phys. 108, 384 (2009), arXiv:0711.1264 [astro-ph] .
  * Garbrecht and McDonald (2018) B. Garbrecht and J. I. McDonald, JCAP 07, 044 (2018), arXiv:1804.04224 [astro-ph.CO] .
  * Hook  _et al._ (2018) A. Hook, Y. Kahn, B. R. Safdi, and Z. Sun, Phys. Rev. Lett. 121, 241102 (2018), arXiv:1804.03145 [hep-ph] .
  * Huang  _et al._ (2018) F. P. Huang, K. Kadota, T. Sekiguchi, and H. Tashiro, Phys. Rev. D 97, 123001 (2018), arXiv:1803.08230 [hep-ph] .
  * Battye  _et al._ (2020) R. A. Battye, B. Garbrecht, J. I. McDonald, F. Pace, and S. Srinivasan, Phys. Rev. D 102, 023504 (2020), arXiv:1910.11907 [astro-ph.CO] .
  * Leroy  _et al._ (2020) M. Leroy, M. Chianese, T. D. P. Edwards, and C. Weniger, Phys. Rev. D 101, 123003 (2020), arXiv:1912.08815 [hep-ph] .
  * Darling (2020a) J. Darling, Astrophys. J. Lett. 900, L28 (2020a), arXiv:2008.11188 [astro-ph.CO] .
  * Darling (2020b) J. Darling, Phys. Rev. Lett. 125, 121103 (2020b), arXiv:2008.01877 [astro-ph.CO] .
  * Witte  _et al._ (2021) S. J. Witte, D. Noordhuis, T. D. P. Edwards, and C. Weniger, Phys. Rev. D 104, 103030 (2021), arXiv:2104.07670 [hep-ph] .
  * Battye  _et al._ (2021) R. A. Battye, B. Garbrecht, J. I. McDonald, and S. Srinivasan, JHEP 09, 105 (2021), arXiv:2104.08290 [hep-ph] .
  * Millar  _et al._ (2021) A. J. Millar, S. Baum, M. Lawson, and M. C. D. Marsh, JCAP 11, 013 (2021), arXiv:2107.07399 [hep-ph] .
  * Battye  _et al._ (2022) R. A. Battye, J. Darling, J. I. McDonald, and S. Srinivasan, Phys. Rev. D 105, L021305 (2022), arXiv:2107.01225 [astro-ph.CO] .
  * Foster  _et al._ (2022) J. W. Foster, S. J. Witte, M. Lawson, T. Linden, V. Gajjar, C. Weniger, and B. R. Safdi, (2022), arXiv:2202.08274 [astro-ph.CO] .
  * Brinkmann and Turner (1988) R. P. Brinkmann and M. S. Turner, Phys. Rev. D 38, 2338 (1988).
  * Choi  _et al._ (1989) K. Choi, K. Kang, and J. E. Kim, Phys. Rev. Lett. 62, 849 (1989).
  * Prokopec  _et al._ (2004b) T. Prokopec, M. G. Schmidt, and S. Weinstock, Annals Phys. 314, 267 (2004b), arXiv:hep-ph/0406140 .
  * Iwamoto (1984) N. Iwamoto, Phys. Rev. Lett. 53, 1198 (1984).
  * Ishizuka and Yoshimura (1990) N. Ishizuka and M. Yoshimura, Prog. Theor. Phys. 84, 233 (1990).
  * Burrows  _et al._ (1990) A. Burrows, M. T. Ressell, and M. S. Turner, Phys. Rev. D 42, 3297 (1990).
  * Iwamoto (2001) N. Iwamoto, Phys. Rev. D 64, 043002 (2001).
  * Dietrich and Clough (2019) T. Dietrich and K. Clough, Phys. Rev. D 100, 083005 (2019), arXiv:1909.01278 [gr-qc] .
  * Harris (2020) S. P. Harris,  _Transport in Neutron Star Mergers_ , Ph.D. thesis, Washington U., St. Louis (2020), arXiv:2005.09618 [nucl-th] .
  * Raffelt (1996) G. G. Raffelt,  _Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles_ (Chicago Press, 1996).
  * Harris  _et al._ (2020) S. P. Harris, J.-F. Fortin, K. Sinha, and M. G. Alford, JCAP 07, 023 (2020), arXiv:2003.09768 [hep-ph] .
  * Blaizot and Iancu (2002) J.-P. Blaizot and E. Iancu, Phys. Rept. 359, 355 (2002), arXiv:hep-ph/0101103 .
  * Thorne (1980) K. S. Thorne, Rev. Mod. Phys. 52, 299 (1980).