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

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.


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 [1]. 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 [2,3] 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 [4,5] 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 hole 1 . 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 [9][10][11][12], polarisation signatures [13] and effects on binary inspirals [14]. Superradiance can also be used to constrain the existence of light fields [15][16][17]. 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 [18][19][20] Superradiance in stars has been studied for some spe-cific cases [21][22][23] 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 φ 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 where the damping is given by Γ φ = − lim p→0 Im Π R (p)/p 0 evaluated in the long-wavelength limit p → 0, where Π R is φ'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 Γ φ 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 where f i are distribution functions for particles in the medium and M is a scattering amplitude.
From the classical equation of motion (1), one can derive an expression for scattering of an arbitrary mode ( , m, ω) into a non-rotating star of radius R • 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 [24], 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, , m) where C n m is a combinatorial coefficient defined later on in Eq. (72), r n = (n + + 1)/(2GM µ 2 ) gives the characteristic radius of the superradiant profile, and ω n gives the real part of the discrete bound state frequencies Eq. (C10). The rotational frequency of the star is Ω.

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 φ.
• The Schwinger-Keldysh closed-time-path (CTP) approach [25,26] 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 A, 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. [27][28][29][30]. We then proceed to derive equations of motion for a inhomogeneous scalar field φ in a general background to linear approximation in φ. 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 [31,32] In Subsection III B, 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.
A. 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 where |0 is the ground state of the theory. Further, one assumes that the one-point function of a generic field ϕ will have a constant expectation value φ(x) = ϕ(x) = const., there are no statistical fluctuations and the system is entirely described by a pure vacuum state. Generally ϕ 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 . 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 φ = ϕ = 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 H and that the state of the system is described by a density matrix ρ. (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 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 τ . The path integral can then be used to evolve forward from t to τ and then backwards again in the reverse direction [25,26], forming a closed loop in the time-domain, as illustrated in FIG. 1. This corresponds to the following generating functional [28,29] It generates all correlators, in particular the one-point function but also the set of four Green functions where we refer to a, b = ± as Schwiger-Keldysh indices and C defines a path ordering along the so-called Schwinger Keldysh contour in FIG. 1. We can then generate 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 We write here G instead of ∆ for generality because these definitions apply as well for other two-point functions than correlators, in particular for self energies.
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 φ = ϕ , 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 [27,29,33] 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 φ, 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 oneparticle effective action. However, this would not amount to a substantial simplification, and in view of applications to nonequilibrium systems accounting for backreactions from φ, we carry out the formal developments here in the 2PI framework.
The 2PI effective action for a number of fields ϕ and their propagators ∆ (where flavour, tensor, spinor and Schwinger-Keldysh indices are suppressed) is given by (12) where S is the classical action and i∆ (0) −1 = δS/(δφδφ) is the classical inverse propagator and Γ 2 is given by (−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 i∆ (0) nor i∆ 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 Γ 2PI on φ and ∆: 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. [30]. 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 ϕ coupling to a massive fermion ψ via (To roll back to the above notation, refer to the scalar and the fermion field collectively as ϕ, i.e. {ϕ, ψ} → ϕ.) 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. δΓ/δφ = 0, leading to where and S and ∆ are the full fermionic and bosonic propagators, respectively. The term with the derivative δ/δφ in Eq. (15) introduces the effect of loop corrections on the dynamics of φ. We take the potential to be of the form From Eq. (16) we see then that ∆ (0) −1 has no explicit φ dependence and is therefore annihilated by the action of the functional derivative with respect to φ. (Recall that the fields φ and the propagators ∆ are independent variables in the 2PI effective action.) Furthermore, for the Yukawa interaction, insertions of φ into diagrams are absorbed within the fermion propagator, cf. Eq. (17). Hence Γ 2 does not explicitly depend on φ and does not contribute to the derivative term in Eq. (15). (Generally, when there are explicit insertions of φ into vertices of 2PI diagrams, the terms from Γ 2 enter the equations of motion as well.) Therefore, in the present example, only the term tr[S (0) −1 S] survives the functional differentiation with respect to φ, which appears in S (0) −1 , leaving us with Tadpole diagram contributing to tr[S(x, x)] at linear order in φ. Double lines denote full fermion propagators in the event that ψ is coupled to some other sector different from φ. Similarly, thick dots may also represent vertex corrections from such a sector. Now, while the fermion propagator iS does not explicitly depend on φ, it does so implicitly through its equation of motion in which φ appears. When itr[S(x, x)] φ=0 = 0 (i.e. for constant, vanishing field configurations), the fermion loop shifts the time-independent solution for the expectation value of φ. This contribution may be removed by adding appropriate counter terms or simply by a redefinition of the field φ by a constant value.
Assuming that the dynamics is adequately described by a linear approximation, 2 we functionally expand the implicit dependence of the solution for iS to linear order in φ. That is, we approximate where Π is simply the scalar self energy evaluated in the background φ(x) ≡ 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: where we have now made the Schwinger-Keldysh indices a, b, c = ± explicit, with h ab defining a Schwinger-Keldysh metric given by h = diag(1, −1) which enforces the signs inherent in Eqs. (7)- (9). Taking a = +, one can see that by using φ + = φ − = φ and the definitions (10) and (11) generalized to self-energies, the final term can be expressed as where Π R (x, y) = Π ++ − Π +− is the retarded self-energy. Equations of this form are similar to those derived in Refs. [31,32] 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 φ(y) with y 0 < x 0 . It may therefore be characterized as non-Markovian. The fact that the convolution of φ with the retarded self energy appears in Eq. (22) is a generic feature when the loop propagators lineraly depend on φ. Another important situation is the quadratic dependence, i.e. for a quartic self interaction ϕ 4 or a coupling to another scalar χ via ϕ 2 χ 2 . Then, tadpole corrections of the form appear in the equations of motion. With the help of a δ-function, these can be of course brought to the form of the self energy term as in Eq. (22).

B. 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 The approximations will rely on the mild dependence of Π 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 ∂ µ x 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 µ . Provided these generic assumptions hold, we can carry out a gradient expansion in terms of q · ∂ 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 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 Π 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 This allows us to express Eq. (22) terms of the Wignertransform of the self energy as (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 φ, Re[Π(q, x)] is an even function under point reflections in the four momentum p, while the spectral self energy is odd. Therefore, Eq. (26) can be written as In this form, it becomes explicit that the gradient expansion and Wigner space representation of the self energy respect that the field φ is real. Furthermore, the dispersive effects from Re[Π R ] and the absorptive ones from Im[Π R ] are thus separated.
The wave equation for φ as in Eqs. (26,28) is still non-Markovian because it is equivalent with Eq. (22).
In general, the non-Markovian property of an integrodifferential 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 ∂ y Π R when these correspond to macroscopic scales and are thus small compared to the scales that determine ∂ p Π R , i.e. the density or temperature of the medium.

Damping in the long wavelength limit and superradiance
The gradients ∂ y φ may or may not be small. A generic situation is however ( This typically occurs when the wave number |k| and frequency |k 0 | (corresponding to ∂ y ) are small compared to the density or temperature T of the medium (corresponding to ∂ 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 ∂ q · ∂ y but the leading one, the contribution involving Im[Π R ] in Eq. (28) takes the form of a classical friction term in covariant form: We typically envisage a situation in which the damping arises as the result of the absorption of a soft φ mode onto the external leg of an existing in-medium scattering process-see Fig. 4. In this case, the leading fourmomentum 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 Π 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 where m and p is the mass of an in-medium particle. Clearly for non-relativistic particles, p 0 m |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 This means that the spatial derivative term in Eq. (29) is effectively suppressed. Then, since Im[Π R ] is odd in q 0 it follows that we can derive an effective damping rate leading to an equation 2. Damping in the short wavelength limit (WKB) For completeness of this discussion, we also consider the case when the gradients ∂ y φ are not small, i.e. when the wave number |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 φ is such that the decay into two or more particles is possible. The corresponding momentum thresholds in Im[Π 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 timeindependent, i.e. Π R (q, x) ≡ Π 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 φ(x) ∼ exp(ik · 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 where ω = (k 2 + m 2 ) 1/2 . Neglecting a possible dependence of Π R on φ, this equation is linear so that the solution is of the form φ(x) = exp(−iω (1) t + ik · x), where the notation ω (1) indicates that i corresponds to a correction to ω due to Π R . Then, also the temporal gradients can be summed, what leads to and therefore In the argument of Π R , we have replaced ω (1) with ω. 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 [24]. 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 ∼ R, and the low energy scale corresponds to wavelengths ∼ 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 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 [36], a spherically symmetric object is described by operators transforming under SO(3).
For example, for an extended object interacting with a scalar field φ, the interaction Hamiltonian in the worldline EFT is where y(t) is the worldline of the extended object and O (n) are effective operators describing the interaction between the object and the scalar field. I, J are spatial indices, so that the operators 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 [24]. 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: (37) where the fields are now understood to be evaluated at the spatial position of the object which we will take to be x = 0. The EFT associated to a rotating object can be inferred through the action of a time-dependent rotation matrix on the operators O n O (n) where, without loss of generality we have chosen coordinates in which the object spins about the z-axis so that The EFT for the field interacting with a rotating object, as as seen by an inertial (non-rotating), observer then reads 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.

A. Superradiant scattering
First, let us imagine the interaction Hamiltonian is switched on over some finite time T . The evolution operator is then given by Let us now consider the effect of this interaction on the absorption of the field φ by the rotating object, following the method introduced in [24]. We compute the absorption of a spherical mode of the field φ with frequency ω and angular momentum labels , m. The probability for this process is: 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 comes from the operator in Eq. (36). Upon inserting (44) into (41) and expanding to leading order in interactions, one obtains where O J1...J (t )O L1...L (t) denotes the expectation value for the initial state |X i . All that remains is to manipulate the matrix elements associated to the spherical modes of φ. This is a lengthy but straightforward computation, the details of which can be found in [24].
The result is where q = |q| and v = 1 − µ 2 /ω 2 is the group velocity of a field mass µ asymptotically far from the star and ∆ (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 timetranslation symmetry means the correlators can only dependent on the relative coordinate t − t . For sufficiently long times the t integral can be approximated by a delta function, multiplied by T leading to We now turn to the normalisation of states in the denominator, which, as we show in appendix B is normalised to where 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π) 3 δ (3) (k = 0) → V . Hence we can write Dividing by T , we then obtain an absorption rate where n r = 1/R is the particle number per spherical shell of unit radius. We can similarly calculate an emission probability noting that the argument of ∆ appears now with the opposite sign to Eq. (51). We can then define an amplification factor for spherical waves as where Φ in = n r v is the ingoing flux, which is related to the outgoing flux by Φ out = Φ in [1+P em −P ab ] and where is the spectral response function for the star. The key point is that ρ (ω) is an odd function of ω, hence for frequencies ω > mΩ energy is absorbed (Z m < 0) whilst for ω > mΩ energy is extracted from the star (Z m > 0). This is superradiant scattering.

B. Superradiant instabilities
We now turn to a similar question of bound states, extending the method of [24]. 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 µ. The scalar field operator can be rewritten in terms of bound state creation and annihilation opera- where a n m are creation operators for a mode (n, , m).
The mode functions are the bound state solutions on a Schwarzschild blackhole background. Expressions for R n are given in appendix C. They are normalised as: (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 |n, , m =â † n m |0 by the star. We have: 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 P abs Whilst the general case for superradiant scattering has been covered in [24], to which we referred the reader earlier in subsection IV A, the general bound state case was not presented there. We therefore collect the relevant details in Appendix C. There, we show that for large T , one finds where r n = (n+ +1)/(2GM µ 2 ) and ω n are the discrete bound state frequencies defined in Eq. (C10) and A n m = !(2 + 1)!!(2 + n + 1)! 8π( + n + 1)n!(2 + 1)! 2 .
We then obtain the absorption rate Similarly, the emission rate is and so the superradiance rate is simply the difference of these two quantities As we have seen in the previous section, and as described in [24], we can obtain ρ by matching to the superradiant scattering amplitude (54), so that scattering off a static star leads immediately to the superradiant growth rate of bound states. By definition, gravitational bound states have frequencies ω < µ below the mass of the field. Hence superradiance requires knowledge of ρ (ω < µ) and since a massive propagating wave has ω > µ, it follows that carrying out scattering of massless modes with ω < µ 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 ρ(ω) in the region ω < µ. We illustrate this in the next subsection for a uniform star model.
The crucial point of this section is that correlators O J1...J (t )O L1...L (t) encode all the information about the response of a point-like star to a harmonic stimulus ∼ e −iω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 ρ is sufficient to calculate superradiance rates and this can be extracted from any physical process to which ρ contributes.

C. Superradiance rate from microphysics
As an example, let us consider superradiance for a scalar field around a uniform star. To obtain the spectral function ρ, we will calculate the scattering of a massless scalar field φ from the star when it is not rotating. For a uniform star the self-energy Π R is constant inside r ≤ R and zero elsewhere. The equation of motion for φ is given by equation (32). Under these approximations, we can compute the absorption of a mode ( , m) by considering solutions of the form which gives a radial equation We recognise this immediately as the equation for spherical Bessels functions which admits two linearlyindependent solutions. Inside the star, by regularity at the origin, the solution must take the form j r ω 2 + iωΓ φ whilst outside it can be most conveniently parameterised in terms of spherical hankel functions h (1,2) (ωr). Hence at infinity, the solution takes the form Imposing continuity of the solution and its first derivative at r = R gives expressions for φ in and φ out . One can then define an amplification factor for a non-rotating star in flat space: Using the above continuity conditions, expanding in the limit ωR 1 and to linear order in Γ φ , we find (recalling we deal first with a non-rotating star) repeating the expression found in [21]. We then equate the amplification factors in equations (70) and (54) for a massless scalar field scattering from a non-rotating star.
It is interesting to compare the amplification factor for a rotating star inferred from Eq. (54), with that of a black hole [2] Z BH where r ± = GM (1 ± √ 1 −ã 2 ) whereã = a/(GM ) and T H = (r + − r − )/(4πr 2 + ) is the BH temperature. Expanding this to leading order in spin and the splitting (ω − mΩ H ), we find 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 → R and identifying the damping Γ φ 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 λ, associated to the microphysics λ −1 Γ φ , 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 [37][38][39][40], as well offering a solution to the strong CP problem 4 [42][43][44] and arising generically in string compactifications [45,46]. Indeed in recent years, neutron stars have become an exciting environment in which to study axions [47][48][49][50][51][52][53][54][55][56][57][58][59] We shall model the nuclear interactions using chiral perturbation theory. In this instance, the principle absorption channel is a + N + N → N + N , which proceeds via the interactions [60] where f ≈ 1 and N = p, n are the proton or neutron. Note that we have chosen a basis in which the axionnucleon interactions are "pseudo-vector" (PV) (axion derivatives coupling to the axial nucleon current) and the neutral pion coupling is "pseudo-scalar" (PS) (nonderivative 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 [61] 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 [60] in which axion-nucleon absorption processes are considered.
We shall now see how the microphysical scattering processes which contribute to Π R and thence to the damping term Γ 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 ↔ N + N . We do not consider diagrams at lower order in L πN N whose cuts would lead to a + N ↔ N + π 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 Π > , which is obtained by summing the internal vertices in the diagram over Schwinger-Keldysh indices. The internal propagators connecting different ± Schwinger-Keldysh therefore each carry ± labels. For Fermions, in momentum space, these free thermal propagators are given by (following the conventions of [30,62]) For bosons the relevant propagators are Here are the Fermi-Dirac and Bose-Einstein distributions respectively with chemical potential µ.
Let us start by computing diagram (i) of Fig. 3. As stated above, we are only interested in processes a + N + N ↔ 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 ↔ N + N or a + N ↔ N + π which are kinematically forbidden for the low-energy axions considered here. The only choice of indices giving rise to a + N + N ↔ N + N processes corresponds to a cut down the centre line. The corresponding expression for this contribution to Π > is then given by Tr where q is the axion momentum. We can now insert the explicit expressions (78)-(79) 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 π , T /m 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 where C is an effective matrix element which results from the trace parts (to leading order in axion momentum q µ ) In what follows we will drop terms directly involving q 2 but retain those arising from (p i · q) 2 in order to make contact with the estimates made in e.g. [60,68]. 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 p i .p j (p.q) 2 ∼ p 2 F m 2 m 2 a whilst dropping terms p 4 q 2 ∼ m 4 m 2 a . 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. (82) into the form of a phase-space integral over on-shell modes. This we do by expanding the delta functions using δ( where E p = p 2 + m 2 . We also use the relation F (−E) = 1 −F (E), whereF (E) is the phase space distribution for the antiparticle given by replacing µ → −µ in Eq. (80). Using these relations and performing the p 0 integrations we arrive at 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,p = (±E p , 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 µ/T 1, so that neutrons dominate over anti-neutrons. Under these approximations, all those terms containing aF in the initial state can be discarded, leading to  Fig. 3. The left and right columns correspond to t-and u-channel processes respectively. 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 [60] we now go to the non-relativistic limit and expand the matrix element to leading order in p i /m. In the limit of nonrelativistic neutrons, we find, to leading order in q: where k = p 2 − 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. [60]. Thus, as expected by the optical theorem, cuts through diagrams in fig. 3 correspond to Π > 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 → 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 where 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 [60,63] spins |M| 2 = 256 3 Note that strictly speaking, as can be seen in [60], 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. (k ·q) 2 → (k ·q) 2 = 1/3. Here hats denote unit vectors. However, by sym-metry we know the final self energy Π > (q) = Π > (q 0 , q) must be isotropic with respect to q so that Π > (q 0 , q) = Π > (q 0 , |q|). It therefore follows that the self-energy is invariant under the averaging over the direction ofq. Formally this means Π > = Π > . Furthermore, when applying the operator · · · to Eq. (85), 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 (88), 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 λ −1 > defined by [64][65][66][67][68] 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 λ −1 > in the literature, we can immediately read off Π > . 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 [64,66,68], in that limit one finds where Q(y) =4 − 1 1 + y 2 + 2y 2 1 + 2y 2 arctan 1 1 + 2y 2 − 5y arcsin 1 and y = m π /(2p F ) with p F the neutron Fermi momentum. We can combine Eq. (27) with the KMS relation Π > (q) = e ω/T Π < (q) to obtain Inserting (90) and (91) into this we can immediately compute the anti-hermitian part of the self energy in the degenerate limit, which gives Putting this together, from Eq. (4), we are able to get a superradiance expressed in terms of microphysical parameters: Γ a n m = 1 18π 5 (2 + 1)!!(2 + n + 1)! ( + n + 1)n!(2 + 1) Let us make some estimates of this rate. We can take typical neutron star values Ω = 600Hz, M = 1.4M and R = 10km and p F 300MeV. For couplings constants we can take f 1 and the maximum limit set by Supernova constraints is g an 10 −9 . Furthermore, for hot young stars, the rate is at its highest, here (although briefly) temperatures could reach up to T = 10 11 K [69] giving an upper bound on the superradiance rate. For these numbers, the = m = 1, n = 0 has a superradiance timescale Γ −1 011 10 13 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 [64][65][66][67][68]. Formally the damping rate of the field in the equation of motion is defined as the difference of two Wightman self energies where, by virtue of the KMS relation, we can write the self-energies more transparently as where B is the Bose-Einstein distribution function defined in (80). Notice therefore that the true mean free path involves two opposite contributions where Π > and Π < 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 [70] which the authors claim corresponds to the mean free path, actually only give the piece Π > , explaining the appearance of a prefactor 1/(1 − e −ω/T ) = 1 + B.
In reality, the authors should have taken both contributions, which would have given a mean free path λ −1 = (ImΠ R /ω).
The fact that it is ImΠ and not just one of Π > or Π < 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 [71] This spectral propagator describes the propagating axion modes in the medium. In the lossless limit ImΠ R → 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) ∝ δ (q 0 ) 2 − E q − Re Π R of a zero width mode where Re Π R encodes the so-called dispersive part of Π R which gives radiative corrections to the axion dispersion relation. This provides yet another justification that Im Π R sets the width of axion states associated to the lifetime of the axion in-medium. The question is then whether the authors of [70] 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 ω T , we have that B 1, so that proved one considers axion frequencies at or above the temperature, B factors become negligible and one finds ImΠ R Π > /2i. This misidentification would therefore only change the mean free path by O(1) so long as one sticks to ω T .
In fact, the contribution to the damping rate associated to Π > corresponds to stimulated emission of axions. This means that for superradiance, since we are in the regime ω 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 ω/T so that to leading order in this ratio we have ImΠ R ω T Π > . Hence in our case, since ω/T 10 −17 for T = 10 8 K and ω = 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 Π > ) 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 [21] and it remained opaque as to how to relate this to microphysics. It is worth noting that [22] made some attempt to interpret the damping coefficients for gauge field in terms of QEDlike 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 → 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 [24]. 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 [22]. 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. 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.