Decay of ultralight axion condensates

Axion particles can form macroscopic condensates, whose size can be galactic in scale for models with very small axion masses $m\sim10^{-22}$ eV, and which are sometimes referred to under the name of Fuzzy Dark Matter. Many analyses of these condensates are done in the non-interacting limit, due to the weakness of the self-interaction coupling of axions. We investigate here how certain results change upon inclusion of these interactions, finding a decreased maximum mass and a modified mass-radius relationship. Further, these condensates are, in general, unstable to decay through number-changing interactions. We analyze the stability of galaxy-sized condensates of axion-like particles, and sketch the parameter space of stable configurations as a function of a binding energy parameter. We find a strong lower bound on the size of Fuzzy Dark Matter condensates which are stable to decay, with lifetimes longer than the age of the universe.


Introduction
Axions are hypothetical particles introduced to solve the Strong CP problem in QCD [1,2,3,4,5,6,7,8]. Axions are self-adjoint bosons, with no conserved discrete quantum numbers to guarantee particle number conservation. The axion potential can be written in terms of an angular variable with a 2π shift symmetry. Axion-like scalar particles also appear in a variety of models beyond QCD, especially in low-energy limits of string theories [9,10,11]. Those axions have properties similar to QCD axions, but their mass scales and decay constants are vastly different.
Axion-like particles are among the prime candidates for the composition of dark matter [12,13,14,15,16,17]. Axions, being scalar bosons, can condense. Axion condensates have been discussed by numerous groups, with condensate sizes ranging from galaxy or galaxy cluster scale [18,19], down to stellar size and smaller (termed "axion stars") [20,21,22,23,24,25], all the way to radii of a few meters [26,27]. Also of interest is the cosmological evolution of the axion field, which has been worked on extensively in [28,29,30,31,32], but this will not be disussed further here.
Uncondensed QCD axions are not stable, as they can decay to photons through a process a → 2 γ, but the decay rate is slow enough such that most axions created after the Big Bang would survive many Hubble times [33]. Axion condensates, however, may also decay slowly due to the self-interaction of axions [34,35]. The self-interaction term of their Lagrangian (for both the instanton [36,37] and chiral [37,38] cases) have only terms containing an even number of axion fields. Thus, disregarding the rare decay into photons, the axion number is conserved only modulo 2.
In a recent paper [34] we have investigated the decay of weakly bound axion stars due to the self-interaction of axions. The decay proceeds mostly through a sequence of processes, where A k is an axion star, which is a condensate containing k axions, and a p denotes an axion in a scattering state with the magnitude of the momentum p. The process (1.1) is the simplest of many possible decay modes responsible for the decay of axion stars. 1 This process is allowed by energy-momentum conservation, provided the binding energy of a bound axion is small enough that a relativistic particle can be produced: δE < 2 m / 3, where m is the mass of a free axion.
In [34] we used an axion field operator, which was the generalization of the field proposed by Ruffini and Bonazzola [39]. To facilitate the decay process, terms creating and annihilating axions in the continuum of scattering states were included in the quantum field of axions, in addition to terms creating and annihilating bound axions [39]. The Ruffini-Bonazzola method is based on taking the expectation value of the quantum Einstein and Klein-Gordon equations in the condensate to derive equations of motion for the metric components and the scalar field. We solved the equations of motions numerically in the weak gravity and weak binding (δE m) limits, to find the wave function of axions in the condensate [24].
The bound axions are not in momentum eigenstates. They have definite energies, but their wave functions extend over the size of the axion star. Accordingly, the bound axions have an extended momentum distribution as well. In the recent publications [35,40], the authors questioned the validity of the decay mechanism proposed in [34], arguing that momentum is not conserved in (1.1), and that the decay rate through this process is exactly zero by the Optical Theorem. These authors have further suggested that one can show the rate to be zero by the classical equation of motion for the condensate. We will address these issues and explain our response in Appendix A.
In the present paper we will apply our method of discussing the decay of dilute axion stars [34] to condensates of cosmological size; such models have been referred to previously as Fuzzy Dark Matter (FDM) [41]. Condensates of galactic sizes have been considered by a number of authors, and typically correspond to a scalar particle mass of m ∼ 10 −22 − 10 −21 eV [42,43,44,41,45,46,47,48,49,50,19]. 2 When such condensates are formed from real scalars, a version of the decay analysis of [34] applies, and we will investigate whether interesting bounds can be placed on these models by taking decays into account. As we will explain in the next section, we will utilize the axion potential with a cosine dependence on the field; other proposals, for example a cosh potential [52,53], have been investigated in the context of ultralight scalars as well.
It is also an aim of this paper to emphasize the inclusion of axion self-interactions in investigations of axion condensates. Although the self coupling λ ∼ m 2 /f 2 ∼ 10 −95 ≪ 1 for typical FDM models (f is the axion decay constant), the astronomical number of axions in a condensate participating in these interactions could lead to large corrections to certain important physical quanties. We investigate the macroscopic properties of these condensates using the fully self-interacting analysis and emphasize the differences from the non-interacting limit.
In the next section, we give a more detailed explanation of how axion star decay through the process (1.1) can be calculated. In Section 3, we will outline the calculation of the wavefunction, following largely [24]. In Section 4, we apply the formulas for the macroscopic properties and decay rates to condensates formed from ultralight axion-like particles. We conclude in Section 5.
We will use natural units throughout, where = c = 1.
2 For a brief but recent review of ultralight scalar field dark matter models, see [51] and references therein.
There is a variety of methods for the quantitative investigation of axion condensates, as classical [33,21,20,54,55], quantum mechanical [56,57,25,58,59,60,61], and field theoretic [26,24,27]. Our field theoretic discussion of the decay of an axion condensates into relativistic axions [34] was based on an extension of the Ruffini-Bonazzola operator [39], by the addition of scattering state contributions. 3 Thus, we proposed to extend the expansion of the boson field using the form [34] where E 0 is the energy eigenvalue of a single bound axion, and where R(r) and a 0 are the wave function and annihilation operator of the axions in the condensate (respectively). ψ f (r, t) is the annihilation part of a complete system of free axion operators expanded in scattering states, where ω p and a lm (p) are the energy eigenvalue and the annihilation operator of the scattering state axion, with quantum numbers l and m, respectively. The functions j l (x) and Y m l (x) are spherical Bessel functions and spherical harmonics, respectively.
The annihilation operator in the spherical wave basis is defined by where a( p) is the annihilation operator for a particle which is the eigenstate of the momentum operator with eigenvalue p (i.e. a plane wave). This annihilation operator and its adjoint, the creation operator, satisfy the commutation relation Note that (2.2) is exactly equal to the negative frequency part of a complete system of free axion states given by 1 which was used in [34] to investigate the decay of QCD axion stars.
We will see later that the bosons produced by the decay of a weakly bound boson condensate are relativistic. Therefore, we have chosen to use free particle states to approximate the scattering states. For the purposes of this calculation, this choice is admissible, because the energy level of produced axions is sufficiently high compared to the effective depth of the potential created by gravitation and self-interactions. This is explained in greater detail in Appendix B.
In a future work we will take into account corrections to this approximation.
Let us consider now (1.1) in the Born approximation. The axion self-interaction potential can be approximated by the so-called instanton potential [36,37] where m and f are the axion mass and decay constant (respectively). In the Ruffini-Bonazzola paradigm, one finds the expectation value of eq. (2.6) transforms the cosine into a Bessel function J 0 [24]. 4 One also finds that the transition matrix element for the process (1.1) is where X(r) = 2 √ N R(r) / f is the rescaled wave function of the condensate which, as we will see in Section 3, can be obtained by solving the equations of motion [24]. We are considering transitions of the form (1.1), where N | is the initial N particle condensate (the left hand side of (1.1)), and |N − 3, φ(p) is the direct product of the final state N − 3 particle condensate and a scattering state φ(p) of momentum magnitude p (the right hand side of (1.1)).
We restrict this work to non-rotating axion condensates; the reason for this is twofold. First, for the potential in eq. (2.6) and the parameters we use here, only the inner cores of galaxies composed of axion particles can be described as a condensate; outside of this inner region, the dark matter halo is described by a virialized gas of particles [19] and cannot participate in the decay processes we describe here. 5 Because the condensed core is small compared to the full radius of the halo, it carries at most a tiny fraction of the angular momentum of the galaxy, so as a first approximation we believe restricting to = 0 angular momentum states is appropriate.
The second reason is that a full treatment of rotating axion condensates has not yet been done, though slowly rotating condensates were analyzed in a particular limit by [14]. This is a topic we hope to return to in the near future.
Because we work in the limit of zero angular momentum, annihilation processes of the form a a → G, where G is a spin-2 graviton [62,63], have a rate of zero. Such an interaction would require the participating axions to have at least of angular momentum each; and even if we accounted for the nonzero rotation speed of the galaxy, by our estimation the vast majority 4 More generally, the annihilation process for k bound axions generates an effective potential proportional to the Bessel function J k , as explained in the Appendix of [34]. 5 For ultralight bosons with repulsive self-interaction, like those presented in e.g. [45,48], the condensates can be much larger, and can even constitute the entire dark matter halo.
of axions in a galactic condensate would have far less angular momentum than what would be necessary for this process to occur.
For static condensates, note that the integration over Ω r in eq. (2.7) vanishes for all but swave axions. 6 In that case, the scattering state axion is described by the zero angular momentum contribution only, |φ(p) = a † 00 (p)|0 . Then the wave function of the emitted axion is The integration over t also fixes the energy of the outgoing axion to ω p = 3 E 0 . The matrix element takes the form [34] where the dimensionless integral is The symmetry of the integrand for the substitution r → −r has been used to extend the integration region to all real values of r, and to switch from sin(p r) to e i p r in the integrand. In the second equality, we expanded the Bessel function J 3 to leading order, an appropriate limit for dilute axion stars. Now observe that for dilute axion stars the radius of the star R is very large. In other words, [24,34]. As a result, the range of p, as represented by the momentum uncertainty δp ∼ m ∆ m, is very small. Then due to the delta function in eq. (2.9), enforcing energy conservation, the emitted axion has a momentum peaked at a very large value, p √ 8 m; as a result, the matrix element (2.9) is very small for weak binding. However, as the binding energy δE increases, M 3 will take larger values and the decay rate Γ ∼ |M 3 | 2 also increases. Now to bring out issues related to momentum conservation, we can define the momentum representation wave function as (2.11) Then we can rewrite eq. (2.9) as (2.12) 6 Should we consider rotating axion stars, higher angular momentum scattering states would also contribute.
Since for weakly bound condensates p √ 8 m, the magnitude of the transition matrix depends crucially on the large q tail of momentum distribution Ξ(q). However, for calculational purposes, it is still advantageous to use (2.9) rather than (2.12). One can rely on the numerical solution of the equations of motion, as explained below, using a simple approach to estimate approximate behavior of Ξ(q) at large q [34].
Note that the process (1.1) is not the only possible channel through which decay can proceed; however, it is by far the dominant process. First, we have shown previously [34] that processes of the form A N → A N −(2 j+1) + a p , are suppressed by higher powers of ∆ for each higher j > 1.
In the weak-binding limit, where ∆ 1, these corrections are completely negligible. On the other hand, this argument breaks down for dense axion stars [27,64], where ∆ = O(1); we will return to this case in a future publication.
Second, there are processes of the form A N → A N −k + µ a p , where µ > 1 axions are emitted at once. Unlike the process (1.1), the emission of µ > 1 axions from a condensate can proceed onshell, implying that the corresponding decay rate has a weak dependence on ∆. Nonetheless, as shown in [34], these processes are suppressed by the very small factor m 2 /f 2 for each additional axion in the final state. Since in FDM m 2 /f 2 ∼ 10 −95 ≪ 1, we can safely neglect these corrections as well. We conclude that (1.1) is by far the dominant contribution to the decay of axion condensates.

The calculation of the wave function of the condensate
We review here the calculation of the condensate wavefunction X(r) [24]. The matrix elements of the rr and tt components of the Einstein equation, along with the Klein-Gordon equation, form a closed set of equations for the metric and the axion field, X(r): where the metric is with δ = f 2 /M P 2 and M P = 1 / √ G = 1.22 × 10 19 GeV (the Planck mass). As above, we have taken only the leading contribution to the Bessel function which represents the self-interaction potential; doing so preserves the leading attractive X(r) 4 interaction term.
Assuming that δ 1, a condition satisfied in applications where gravity is weak (Newtonian limit), we can write A = 1 + δ a and B = 1 + δ b, where a, b = O(1). Furthermore using the large radius approximation and the definition of the scale parameter ∆ = can introduce dimensionless radial coordinate as x = m r ∆. In that case the axion field also scales with its engineering dimension, such that X(r) = ∆ Y (x), leading to the following system of equations for a, b, and Y in leading order of ∆ and δ [24]: is the effective coupling constant of the field Y (x) to gravity. Further details, and a more comprehensive justification of this double expansion of the equations, can be found in [24].
Solutions of the equations of motion (3.5) correspond to ground-state configurations of axions, which can be stable or metastable. In [24], we solved these equations and found a spectrum of solutions which were parameterized by κ (or, equivalently, by ∆). When applied to QCD axion parameters m = 10 −5 eV and f = 6 × 10 11 GeV, we found a maximum gravitationally stable mass of M c ∼ 10 19 kg. 8 By rescaling these solutions to values of m and f appropriate for FDM, we can analyze the properties of galactic-size condensates in a way that includes the self-interaction term in the potential. The physical interpretation of these condensates is that they form the cores of FDM halos; they are surrounded by a virialized distribution of axions which extend to the outer edge of the dark matter halo.
To analyze the decay of these condensates through processes like (1.1), we investigate the singularity structure of solutions of (3.5). Now, (3.5) is a system of equations with two singular points, x = 0 and x = ∞. Using boundary conditions we require that the solutions are regular at x = 0 and decrease exponentially at x → ∞. The solutions are even functions of x, so they also approach zero at x → −∞. Then they are real analytic functions at −∞ < x < ∞ and can be continued into the complex plane of x. As they fast vanish at infinity, the contour of integration can be moved up along the imaginary axis in the rescaled version of the integrals in eqs. (2.9) and (2.10), until we reach a singularity in the complex plane. The contribution of that singularity dominates the decay rate integrals at large momentum.
It is easy to show that the Klein-Gordon equation (3.5), in which the the leading order singular terms are retained, is 9 This definition of κ appears to differ by a factor of 8π compared with [34], because in that work we wrote δ in terms of the reduced Planck mass. In fact, the two definitions of κ are equivalent. 8 The maximum masses for attractive interactions were discussed by Stoof [68] in the context of condensed matter BECs, and in the context of boson stars by Chavanis and Delfini [56,57]. 9 In fact, this expression contains the next to leading order term proportional to Y (x).
Near the singular point x = i ρ, this has a solution of the form The parameter ρ is an integration constant, having a one-to-one relationship with the rescaled central density of the axion field Y (0) 2 , and in turn, with mass and the radius of the condensate.
In fact, high order Taylor series expansion of equations around x = 0 show that the singularity closest to the origin is indeed of the form (3.7), connecting the value of Y (0) with ρ [34]. In principle, gravitational interactions have an effect on the solutions Y (x). In practice, however, the term in the equations of motion (3.5) which couple Y (x) to gravity give a subleading contribution to the singularity. We can therefore solve the equations in the limit κ 1, i.e. where gravity decouples. In that limit, the nontrivial solution has Y (0) = 12.268, which implies a fixed value ρ = .603156.
Finally, we can rewrite the Fourier transform of eq. (2.11) as In this form, it is clear that at small ∆ the singular term of (3.7) term dominates the integral.
To calculate the decay rate, we follow the procedure of [34]: we take the leading order solution for Y (x) near the singularity x = i ρ, given by eq. (3.7), and evaluate I 3 (p) in the matrix element of eq. (2.9). The result is The decay rate for the process (1.1) is then where T is the duration of the decay process. Then the lifetime of the condensate through this Further details regarding the evaluation of eq. (3.11) can be found in [34]; the result for the where y M 75.4 is determined by the relationship between M and ∆ in the large ∆ region [24]. The lifetime is a monotonically decreasing function of ∆ in the relevant range; in the case of QCD axions, we found in [34] that above a value ∆ .05 − .06, axion condensates become very unstable to decay to relativistic axions, their lifetimes becoming shorter than the age of the universe. We will examine the consequences of this fact in the context of ultralight axions in the next section.  [41], have been written about extensively [42,43,44,41,45,46,47,48,49,50,19]. While there are significant constraints on these models 10 , they remain a viable alternative to WIMP or QCD axion models of dark matter. 10 While this work was being finalized, a paper appeared suggesting a strong tension between the preferred mass scale for FDM, m ∼ 10 −22 − 10 −21 eV, and data from Lyman-α forest simulations [65]. We will not comment here about whether such constraints could rule out FDM as a viable paradigm.  We will consider an ultralight axion in this context, using the potential of eq. (2.6). The mass of the ultralight axion in question will be taken to be m ∼ 10 −22 eV, which gives the right approximate scale for the size of dark matter halos [19], where v is the velocity in the halo. This choice is also consistent with the known epoch of matter-radiation equality. Further, a decay constant of f ∼ 5 × 10 16 GeV naturally leads to the correct relic density, and can thus account for the observed dark matter abundance [19]; however, to remain as general as possible, we allow f to deviate from this value by a few orders of magnitude. At the upper limit of what we consider, f = 10 18 GeV is still below the Planck scale, implying that the parameter δ = f 2 /M 2 P ≈ .007 1; thus, the weak-gravity approximation holds reasonably well over our entire range.
In [24], we found numerically the solutions to the system (3.5) over a wide range of κ. In the FDM picture, these solutions correspond to the cores of FDM halos discussed (most recently) in [19]. We found that there exists a maximum mass at κ ≈ . 34 We observe in Figure 1  We can also analyze the relationship between M and R 99 , which were investigated for both attractive and repulsive self-interactions in [56,57]. In a recent paper [19], the authors present a bound on the product where R 1/2 is the radius inside which .5 of the mass of the condensate is contained; the inequality is saturated for stationary, ground state configurations, i.e. for the condensates considered here.
In our calculation, on the stable branch of solutions (where κ > .34), we find the product interactions.
It is worth noting also that in the limit f ∼ M P , the mass-radius relationship for axion condensates approaches that of a black hole. This is easy to see using eq. (4.2): implying that R 99 ∼ R S , the Schwarzschild radius.
We must also ensure that the weak-binding approximation, on which our analysis [34] and the classical one of [19] depends, is also valid. 12 We observe in Figure 2 that cores of radius R 1 pc have ∆ .3, and become relatively strongly bound. Such cores would not be well-described by our weak-binding analysis.
An estimate of the decay rate, obtained from the expression derived in [34], is given in eq.

Conclusions
In a previous publication [24] we established scaling relations for the mass and radius of weakly bound condensates of interacting axions, as functions of the mass of the axion, its decay constant, and the particle energy (or alternatively the central density). We also found the maximum mass and size of the bound states as functions of those parameters. In this paper we have applied those results to condensates of axions forming FDM, providing corrections to similar calculations which neglect the self-interaction of axions [44,41,45,47,49,50,19].
In another publication [34] we developed methods to calculate the lifetime of axion condensates due to their self-interaction, through the four-axion interaction term in which three bound axions produce a single free relativistic axion. Here we have applied those results to estimate the lifetime of condensates formed from FDM. We have found that, provided the decay constant of FDM axions satisfies f 0.05 M P , all condensates having binding energy smaller than that of those of maximal mass have lifetimes greater than the age of the universe making them viable candidates for forming central regions of galactic halos. We have also explained in details the decay mechanism described in [34] and further clarified the justification of its validity.
The methods we have described here, based on previous work in [24,34], rely on a double expansion to leading order in δ and ∆. This is appropriate for so-called dilute axion stars, which are weakly bound. However, it has been pointed out that an effective short-distance repulsive interaction in the axion potential also gives rise to dense axion stars [27,64,73], which are at least energetically stable. We plan to extend our methods to this regime to analyze the properties of these states in the near future.
Collapsing boson condensates have been investigated by a number of groups [67,68,69,70,64,71,72,73]. Recently, we found that supercritical QCD axion condensates, having masses larger than the maximal allowed stable mass, collapse towards the global minimum of the effective axion potential [64,73]. Similar arguments indicate that FDM axion condensates which exceed the maximum mass M c will also collapse in this way. Such supercritical condensates can form during galactic collisions, in a manner similar to the mechanism outlined in [74]; such events would lead to collapse, causing the condensate to emit a large number of relativistic particles. Consequences of such events will be studied in a future publication.
The results of our paper, "The Lifetime of Axion Stars" [34], have been called into question in recent publications [35,40]. Before discussing this issue in detail, we review the premises of our work. We proposed a way to discuss the decay of axion stars through the repeated elementary decay mode where A k denotes an axion star, which is a condensate of k axions and a p denotes a (relativistic) free axion labeled by its momentum p. This calculation was performed using an extended axion field operator, eq. (2.1), which included both bound and scattering state contributions.
First of all, we need to establish the fact that there is no conservation law that would forbid (A.1). Axions are real bosons, and consequently the axion number is not conserved. Axions, being coupled to the electromagnetic field as Φ E · B, can decay to photons through the slow process a → 2 γ, though this decay process does not significantly affect the lifetime of axion stars [33]. Disregarding axion decay to photons, the axion number is conserved modulo 2, as the self-interaction terms of the axion potential contains only even powers of the axion field, through a dependence of cos(Φ / f ) in eq. (2.6). Momentum and energy conservation would allow the decay process to proceed even if the condensates were in momentum eigenstates: a decay process is always allowed if the sum of the masses of the decay products is smaller than the mass of the decaying object, unless the conservation of discrete quantum numbers prevents it.
At any rate, the condensates represented in (A.1) are not in momentum eigenstates. They are quantum objects, which have extended wave functions localized on a large radius R. Consequently, though they have mean momentum of zero, their momentum distributions are smeared.
In fact, the momentum distributions extend from zero to infinity, albeit with fast decreasing amplitudes.
We have to emphasize that a Gross-Pitäevskii approach, being the non-relativistic limit of the relativistic quantum field theory, cannot be used to discuss particle number violating processes.
In the non-relativistic limit the axion number is conserved and process (A.1) is forbidden. Note however that a method to adapt the Gross-Pitäevskii approach to these processes was formulated by the authors of [75]. We also point out Gell-Mann's totalitarian principle [76], borrowed from T. H. White's "The Once and Future King", which when it is applied to physics is: "Everything not forbidden is compulsory." This principle implies that the process (A.1) should exist. In the present context, this principle also implies that taking the nonrelativistic limit discards important contributions that make such a transition possible.
Now the questions raised in [35,40]  The argument of [35,40] pertaining to the microscopic side (a) is that in the overwhelmingly dominating Born approximation, the process around which the decay process (1.1) is built is where a c represents a bound axion and a f represents the final state axion, is inadmissible due to energy-momentum conservation. The total energy of the three bound axions is E tot 3 m, which is certainly sufficient to produce a free relativistic axion. However, if the three axions were in momentum representation, having zero momentum in the rest system of the axion star, then the axion in the final state would not have the required momentum of p √ 8m.
To resolve this problem, observe that the three axions contributing to the decay process are not at rest. The notion of "rest" is tied to particles in momentum representation, and for the axions in the axion star, this is not so. They have momentum space wave functions, which, though peaked near zero, extend to large momenta (albeit with very small probability).
As we discussed above, the uncertainty of the momentum of each axion is δp ∼ R −1 . The probability that three of the bound axions have sufficient total momentum to create a free axion, p = √ 9 E 0 2 − m 2 , decreases rapidly with the size of a condensate, though it is never exactly zero. Consequently, local momentum conservation always allows the decay of condensate. The question of whether this decay process affects cosmology, due to the survival or non-survival of the condensate, is a question of numerical calculations and depends on the parameters of the axion theory, as well as the size and mass of the condensate.
Another argument presented in [35,40]  Consider, however, that the condensate itself is not in momentum representation either, unlike condensates in a container, for which the boundary conditions are set by the wall of the container. Still, one could evoke the valid argument that just like the expectation value of the coordinate, the expectation value of the momentum of the condensate must be zero.
Consequently, the only constraint we can impose on the decay process is the conservation of the expectation value of the momentum. Now, the condensate of the final state of (A.1) is not in momentum representation either, though its average momentum is zero as well. Consider now the created scattering state axion. It is produced as a zero angular momentum spherical wave, going with the same probability into every direction. Though the magnitude of the momentum is sharply peaked at a particular value, the spherical wave also has a vanishing average momentum.
Thus, the average momentum is conserved in (A.1).
One should not confuse the emission of the axion with its detection. Suppose we detect an emitted axion, the decay product of (A.1). By performing a measurement we alter the system.
Just like in the case of the collapse of a simple wave packet by performing a measurement, the conservation of the average momentum is valid only if we include the measuring device, which absorbs the appropriate amount of momentum, to make the average momentum of the complete system zero.
We turn now to the last critique presented in [35,40], namely that the classical equation of motion for the axion star precludes a linear coupling to a scattering state axion. In [35], the authors write about our previous work [34] as follows: They expanded the scalar axion field φ around the classical field φ 0 : φ = φ 0 +φ, wherẽ φ is the quantum fluctuation field. The Hamiltonian includes a term proportional to φ 0 3φ from the axion interaction potential. They claimed that this term produces We can see this more quantitatively by analyzing the equation of motion. In an "exact" field theory, any matrix element of the form ψ|KG[Φ]|ψ = 0, including However, because this exact theory is not known, one is restricted to using some ansatz for the field Φ. We have chosen the Ruffini-Bonazzola ansatz, expanded to include the scattering state solutions, which was given in eq. (2.1). With this choice, the wavefunction R(r) (which is equivalent to the Gross-Pitaevskii wavefunction), does not satisfy eq. (A.5), and in this same parameterization, the matrix element for the process (A.1) is nonzero as well.
We are working to extend the Ruffini-Bonazzola paradigm so that eq. (A.5), as well as higherorder expressions in the operator Klein-Gordon equation, can be simultaneously satisfied. In this extension, we find that the rate for the process (A.1) is still nonzero, and is equal at leading order to the result we obtain here. We will present these and related findings in a future publication.

Appendix B: The continuous spectrum in the Ruffini-Bonazzola equations
The most general bound solution, discussed in [39], for a real scalar field in spherically symmetric metric of eq. (3.4) and satisfying the non-interacting Klein-Gordon equation is where c nlm are arbitrary amplitudes of the nth bounds state solution with angular momentum l. The wave functions R nl (r) satisfy wave equations To be able to describe condensates of bosons, Ruffini and Bonazzola [39] introduced second quantization by promoting coefficients c nlm to creation and annihilation operators c nlm → a nlm and c * nlm → a † nlm ; these operators satisfy [a nlm , a † n l m ] = δ nn δ ll δ mm . (B.8) Now notice that Φ b is not the most general solution of the wave equation. Because the attractive gravitational potential, which makes bound states solutions possible, vanishes at large r, (B.7) has scattering state solutions as well. After quantization these states can be written as Φ s = 1 2 π 2 d 3 k 2 ω k ml f l (k)Y m l (θ, φ) e −i ω k t a lm (k) + h.c. , (B.9) where k = ω 2 k − m 2 is the momentum and [a lm (k), a † l m (k )] = (2 π) 3 2 ω k δ ll δ mm δ(k − k ). (B.10) The complete set of solutions to the Klein-Gordon equation is Φ = Φ b + Φ s .
The Ruffini-Bonazzola method was used by Barranco and Bernal [26] and by us [24] to describe axion stars. It is sufficient to use of quantized field Φ b in leading order to describe static axion stars. However, the term Φ s becomes significant if we notice that in first order of the expansion of the axion potential using Φ results in an operator This has been described in detail in [34] and also applied in the present paper.
The second subject we discuss in this Appendix is the reason why free spherical waves, f l (k) → j l (k r), can be used to calculate the decay rate of axion stars. This is the consequence of the fact that dilute axion stars, those we consider in this work, are weakly bound. In fact, only axion stars whose particle energy satisfies |E 0 − m| 0.002 m, corresponding to a value of ∆ ≈ .05, could survive from the big bang until the present epoch. where we defined an effective potential Now considering that b(r) < 0, b (r) > 0 and X (r) < 0 over the whole range of r, as shown by the numerical calculations [24], the effective potential V eff satisfies 3 (E 0 − m) V eff (r) ≤ 0 (B.14) The quantity 3 (E 0 − m) −0.006 m at ∆ .05. This shows that the depth of the effective potential is of O(10 −3 m) or smaller, implying that produced relativistic axions, which have energy E 3 m, can well be regarded as free. The situation would be different in a discussion of dense axion star states, where the binding energy is much more significant; we will return to this topic in a future work.