Collapse of Axion Stars

Axion stars, gravitationally bound states of low-energy axion particles, have a maximum mass allowed by gravitational stability. Weakly bound states obtaining this maximum mass have sufficiently large radii such that they are dilute, and as a result, they are well described by a leading-order expansion of the axion potential. Heavier states are susceptible to gravitational collapse. Inclusion of higher-order interactions, present in the full potential, can give qualitatively different results in the analysis of collapsing heavy states, as compared to the leading-order expansion. In this work, we find that collapsing axion stars are stabilized by repulsive interactions present in the full potential, providing evidence that such objects do not form black holes. In the last moments of collapse, the binding energy of the axion star grows rapidly, and we provide evidence that a large amount of its energy is lost through rapid emission of relativistic axions.


Introduction
The axion, a pseudoscalar particle originally associated with a solution of the strong CP problem in QCD [1][2][3][4][5][6][7][8], has been analyzed in a variety of astrophysical contexts, particularly in cosmological evolution [9][10][11][12][13] and as a candidate for dark matter [14][15][16][17][18][19]. Axions can condense into gravitationally bound objects, either in the early universe through largescale overdensities in a coherent axion field (called "miniclusters"), or through gravitational cooling and collapse (called "axion stars") [20][21][22]. The masses of weakly bound axion stars have been computed previously [23,24], and they are bounded above by gravitational stability [24][25][26]. Axion stars which exceed this maximum mass M c have a fate which remains an open question. Some authors [27,28] suggest that such configurations collapse to a compact, very dense state. Recently, the author of [29] examined the collapse of a boson star with an attractive self-interaction which has M > M c , using a dynamical equation derived from a Gaussian ansatz for its wavefunction. The author found that, as its potential is unbounded from below, a star of this kind collapses all the way to its Schwarzschild radius and forms a black hole. Indeed, the leading axion self-interaction is attractive; axionic or other bosonic objects with repulsive interactions have been considered by [30,31]. However, the axion potential contains additional terms which become increasingly important as the axion density becomes large. It is thus plausible to ask whether these higher-order terms, some of which give rise to repulsive self-interactions, can stabilize the collapsing axion star prior to its formation of a black hole state. In this note, we will consider the consequences of including the full self-interacting axion potential in the collapse of heavy, weakly bound axion stars.
In Section 2, we review the nonrelativistic limit of axion field theory in the description of axion stars; then in Section 3, we outline the variational method used to find energetically stable bound states, and the computation of the total collapse time for large mass solutions. We estimate the binding energy in Section 4, in the initial and final states, but also dynamically in time during collapse. As the binding energy increases, it is known [32] that the rate of decay for axion stars, through an annihilation process which ejects relativistic axions, rises quickly. We thus investigate whether collapsing axion stars emit a large fraction of their energy and decay due to quantum mechanical effects. Finally, we outline our conclusions in Section 5.

The Non-Relativistic Expansion for Axion Stars
Axions are real scalar fields, but in the nonrelativistic limit can be described by a complex wavefunction ψ, using the expansion [25,33] which preserves the Hermiticity of the axion field φ. At low temperatures, the wavefunction ψ describes collectively a condensed state of N axions, termed an axion star, and is normalized as d 3 r|ψ * ψ| = N . The Klein-Gordon equation for φ, expanded using eq. (2.1) in the non-relativistic limit, yields the Lagrangian density is the low-energy axion potential, with m and f the mass and decay constant of the axion, respectively. The gravitational potential representing the self-gravity of the condensate, can be added by hand [25,34]. Then the quantity is conserved. Here W (φ) describes the quantum self-interactions of the axion field, Note that the mass term in the first equality of eq. (2.5) is included to account for a cancellation in the non-relativistic limit between the potential and kinetic terms in L. In the nonrelativistic limit, the total energy per axion is E tot /N m; that is, the binding energy E tot /N − m in the axion star is small. In that case, we can expand eq. (2.5) using eq. (2.1) and drop the rapidly oscillating terms containing extra factors of e ±imt . The resulting equation of motion for ψ is the nonlinear Schrödinger equation.
We derive the total energy from eq. (2.4) in the following way. The nth term in eq. (2.5) contains the factor where 2n C n are binomial coefficients. Dropping the rapidly oscillating pieces, we obtain where J 0 (x) is a Bessel function of the first kind. 1 Including also the kinetic and gravitational terms, the total energy functional has the form A minimum of the energy correponds to a stable bound state, an axion star. Typically, one expands the Bessel function in eq. (2.7) to obtain the leading self-interaction term, which is proportional to (ψ * ψ) 2 . This leading self-interaction is attractive, and as we will explain below, this implies that the potential appears unbounded from below as the axion star size decreases, R → 0. There can nonetheless exist local energy minima, corresponding to metastable states which are dilute and weakly bound. However, there exists a critical particle number N c above which no energy minimum exists, local or global. As a result, it is often assumed that an axion star with M > m N c , being gravitationally unstable, will collapse all the way to a black hole state. A full description of this process can be found in [29], who used a Gaussian ansatz for the wavefunction and calculated the time for collapse to a black hole, which was on the order of an hour.
The full axion self-interaction potential, given by eq. (2.6), contains additional terms beyond the attractive (ψ * ψ) 2 , which depend on increasing powers of the field ψ. Indeed, these higher-order terms, beginning with a repulsive (ψ * ψ) 3 term, become increasingly relevant as the system increases in density, and we wish to investigate whether these terms have the effect of stabilizing the potential against complete collapse. To this end, we will examine the energy functional, including higher-order interactions, and determine whether the endpoint of collapse can lie at a radius greater than the Schwarzschild radius of the axion star. Such a result would be evidence that axion stars stabilize before they collapse to black holes.
We will use a variational ansatz for the wavefunction to calculate the energy in eq. (2.7) as a function of the condensate size, in order to estimate the positions of any energy minima. Using the result of [24], we know how the macroscopic parameters of a weakly bound axion star, the radius R and the axion number N , scale with the dimensionful parameters of the theory; we thus define the dimensionless quantities ρ and n by where δ ≡ f 2 /M 2 P and M P = G −1/2 is the Planck mass. For QCD axions, typical values are m = 10 −5 eV and f = 6 × 10 11 GeV, implying δ = O(10 −14 ) [35].
We will use a single variational parameter, the rescaled radius ρ, at fixed rescaled axion number n. Then the general form of a variational ansatz will be where at fixed ξ the function F (ξ) is independent of ρ, n, and δ. Then substituting the ansatz into the normalization condition of ψ gives the normalization constant as where we introduced the notation Using (2.7) we obtain for the energy functional and where we defined the functions Note that just like C k , also B 4 and D 2 are independent of the physical parameters n, ρ, and δ. The leading-order approximation of (3.5) is obtained when we take the small δ limit, at which (3.10) The minimization of (3.5) with respect to ρ locates the radii of metastable minima and maxima of the binding energy. The condition for the existence of metastable states constrains the reduced particle number n to a finite constant. Restricting ourselves to leading-order of δ, which is δ = O(10 −14 ) in QCD (for f = 6 × 10 11 GeV and m = 10 −5 eV), we obtain For n < n c , there exists a metastable minimum of the energy at a reduced radius which, at n = n c , has a value of

Gaussian Ansatz
Following [26,29], we use a Gaussian ansatz to approximate the axion star wavefunction: (3.14) which corresponds to eq. (3.2) with w = √ N /(π 3/4 σ 3/2 ) and F (ξ) = e −ξ 2 /2 . Note that when we talk about the "size" of such a condensate (whose wavefunction extends to r → ∞), we refer to the conventional R 99 , inside which .99 of the mass is contained. For the Gaussian ansatz, this occurs not at σ, but at a value closer to 3σ. Note also that we define ρ below using eq. (3.1) with R = σ, not R = R 99 .
The energy functional, given by eq. (3.5), depends on the coefficients computed using the Gaussian function in eq. (3.14). Written out explicitly, we have where in the second equality we have expanded J 0 and integrated term by term. Though no closed form exists for the integral in eq. (3.16), we show in the Appendix that it is finite as ρ → 0, and that the kinetic energy term (which is proportional to 1/ρ 2 ) dominates in this region. Consequently, the total energy is bounded from below in this formalism, and always has a global minimum. 2 In Figure 1, we show the position ρ GM of this global minimum of the energy functional in eq. (3.16) as a function of n. The global minimum always lies at a very small radius: ρ GM = O(10 −7 − 10 −6 ) depicted in the plot correspond to R GM ∼ 5 − 50 cm. Given that the Schwarzschild radius R S = 2 M/M 2 P , the ratio R/R S = ρ/(2 n δ). Thus ρ S = 2 n δ so at n = n c , ρ S = 10 −13 10 −7 < ρ GM , and this possible endpoint of collapse is not a black hole.  .16), as a function of the reduced particle number n. The global minimum always lies at a radius ρ GM > ρ eq , the position at which the kinetic energy ∼ 1/ρ 2 becomes dominant, which is depicted by the horizontal black line. The normalized energy per particle coming from the self-interaction is shown to be a constant −1/2 in the small ρ limit, so we can estimate the value of ρ 1 at which the kinetic and self-interaction energies are of the same order; we find comparable magnitudes at a radius of ρ eq ∼ 10 −7 , corresonding to roughly R eq ∼ 5 cm. ρ eq is shown as a horizontal black line in Figure 1. This radius R eq is of the same order as the axion reduced Compton wavelength, λ c = /m c ∼ 2 cm. It should be noted that on length scales of O(λ c ), neglecting higher powers of e ±imt in the expansion of eq. (2.1) would fail, as special relativistic corrections to the kinetic energy could be large. Nonetheless, weakly bound stars have radii much larger than this, and as we describe below, even collapsing stars are well described by the non-relativistic approximation until the last moments of collapse. We have estimated the leading correction to the kinetic energy, which is ∼ p 4 , and in the range of ρ considered here, its expectation value is down by a factor proportional to δ/ρ 2 1 compared to the leading-order term. We will thus postpone any further consideration of these relativistic corrections to the energy, which will be addressed in a future publication.
In this work we analyze the low-energy axion potential in eq. (2.3), sometimes called the instanton potential. But it is well-known (see e.g. [35,36]) that an improved approximation is the chiral potential where m π and f π are the mass and decay constant of the QCD pion. This expression takes into account the non-perturbative effects of up and down quark masses m u and m d . We find that substituting eq. (2.3) with this chiral potential does not qualitatively change the conclusions of this work: the global minimum of the energy in Figure 1 shifts down by at most a few percent, still significantly larger than the corresponding black hole state. We put off any further discussion of the chiral potential to a future publication. We also consider the effect of including a finite but increasing number of terms in the series of eq. (3.17). Because it has no closed form resummation, what is typically done is to truncate the series at some maximum k = K. We denote the truncated energy by E K (ρ), so that lim K→∞ E K (ρ) = E(ρ). We also define a dimensionless truncated energy The minima of e K (ρ) should, at sufficiently large K, approximate well the stable bound states of the full energy function.
The existence of a global minimum of the full energy functional in eq. (3.17) has important consequences. In particular, we have pointed out above that this minimum lies at a radius many orders of magnitude larger than the Schwarzschild radius of the axion star, providing evidence that such objects do not collapse to black holes. Further, we note that the terms contained in the series of eq. (3.18) alternate between attractive and repulsive interactions, even and odd k respectively. But as a result, a truncated energy e K (ρ) in eq. (3.18) for any even K has no global minimum, and thus such a truncation removes the possiblity of approximating the stable radius of the full energy functional. We thus submit that when considering dense configurations of axions or collapse of axion stars, it is important to truncate the series on a repulsive term to preserve the global minimum. The leading-order interaction term is contained in e 0 (ρ), and has been considered in great detail previously [24][25][26]29]. It has been pointed out that there exists a maximum particle number N = N c above which no stable energy minimum exists. This critical value corresponds to a radius of R 99 ∼ 500 km for QCD axions [24], and is approximated to the correct order of magnitude by the Gaussian ansatz, which gives a radius R 99 ∼ 200 km. In our notation, this critical particle number occurs at n c = 2π √ 3 and at a radius ρ * = 3/32π, for the Gaussian ansatz. We will use these as benchmark parameter values as we analyze the consequences of additional interaction terms in the axion potential. The energy functional in the vicinity of this minimum is shown in Figure 2. It is also worth noting that the inclusion of additional terms in the self-interaction potential introduces negligible differences in this range of ρ = O(1); the leading expansion is an extremely good approximation in this region. But as noted above, any e K (ρ) for even K (e.g. e 0 (ρ)) is unbounded from below and will not be applicable in approximating the global energy minimum of the full potential, which is at ρ 1. We turn now to e 1 (ρ), including the leading repulsive interaction which originates from a (ψ * ψ) 3 term in the potential: The energy e K (ρ) multiplied by the small parameter δ for N = .9N c at increasing odd orders in K: K = 1, K = 3, K = 5, and K = 7. The existence of a dense global energy minimum is preserved at each order, but shifts to smaller radii as K increases. The repulsive kinetic term ∼ 1/ρ 2 dominates the total energy at ρ = ρ eq ∼ 10 −7 .
In this case, the energy is bounded from below and has a minimum at a very small radius ρ = ρ D (in contrast to the result using only e 0 ). At these small values of ρ, the energy is well approximated by the self-interaction terms only (gravity and kinetic energy are negligible); thus we can use the analytic expression ρ D ≈ 2 π n δ 3 7/2 1/3 to approximate the position of the global minimum. At n = n c , ρ D ≈ 7 × 10 −6 , corresponding to R 99 ∼ 7 meters. Comparing with the global minimum of the full energy in Figure 1, we find a difference of only about a factor of 3 − 4 near this value of n ∼ n c , a reasonable order of magnitude agreement. This justifies our truncation of the energy at the leading repulsive term, i.e. e 1 (ρ), in this analysis. The difference between ρ D and ρ GM does become large if n increases far above n c . We find that the existence of a dense global energy minimum is preserved at any odd K in the approximation of eq. (3.18), and at increasing order, shifts to smaller radii (see Figure 3). Nonetheless, the kinetic energy term dominates the full potential below ρ eq ∼ 10 −7 , and the global minimum of the full energy is at ρ GM > ρ eq , for any n.
The collapse of dark matter halos consisting of condensed scalar particles was examined by [37], using a time-dependent formalism that originated in [38], and utilized by [26,39].
The application of this method to an axion star, at leading-order in the self-interaction potential, was recently performed by [29]. This collapse process is described by the dynamical equation where α = 3/4 for the Gaussian ansatz, E(R) is given by eq. (2.7) and E tot is a constant. R(t) is the size of the condensate, which varies with time during collapse. For a condensate with size R 0 at t = 0, the time required to reach some other size R(t) is given by , (3.20) where in the second equality we have rescaled the dimensionful quantities.
In the analysis of [29], E(R) was approximated by the leading-order expression E 0 (R), and the collapse from R 0 = R * to R → 0 was shown to last for a time which was on the order of an hour. We wish to investigate the effect of additional self-interactions in the axion potential on the collapse process. Including the first non-leading interaction piece, i.e. using e 1 (ρ), we have found that a global energy minimum exists at ρ D ; thus, we integrate eq. (3.20) not from ρ = 0 but rather from ρ = ρ D .
If the axion star begins its collapse at ρ 0 = ρ * , then of course at n = n c the collapse time is formally infinite, because the potential is flat at ρ * . We consider values of n which are slightly larger than n c and see how the collapse time changes. We also investigate the change in collapse time as the starting size ρ 0 deviates from ρ * . This latter case could be of interest, say, if axion star collapse can be catalyzed by collisions with other astrophysical sources. In that case, even condensates with N < N c can collapse, provided some catalyzing interaction which reduces its initial radius to R 0 < R * . These considerations are represented together in Figure 4.
We can also track the radius of the axion star as a function of time, throughout the collapse process; see Figure 5. For a large portion of the total collapse time, the radius changes little, as the star rolls slowly down a shallow potential, but later collapses fast to the dense minimum of radius ρ D .

Cosine Ansatz
The Gaussian ansatz is believed to be a reasonable approximation to the axion star wavefunction. However, in order to verify that our results are not an artifact of the wavefunction one chooses, we present a second ansatz for the variational analysis:  A comparison of the two ansätze we use is shown in Figure 6 for the same total size. 3 The energy functional, rescaled and truncated as above, depends on the coefficients This implies that, for the cosine ansatz, As before, we minimize the approximated energy e 1 (ρ) with respect to ρ, and find both a dilute and a dense minimum. The dilute minimum disappears above a critical particle number, corresponding to n c ≈ 12.6, where the radius is ρ * ≈ .44 (around 200 km). The dense minimum is at approximately ρ D ≈ 7.4 × 10 −6 n 1/3 , within a factor of 2 of the result from the Gaussian case, ρ D ≈ 4.5 × 10 −6 n 1/3 .

Decay of Collapsing Solutions
In a previous work [32], some of us found that axion stars can decay through repeated occurrences of the a process which ejects relativistic axions from the star. Such a process is not forbidden by any symmetry because axions, being Hermitian fields, do not have a conserved number, and because bound axions, along with the axion star itself, are not in momentum eigenstates. To describe this interaction, the spectrum of bound states describing the axion star was extended by a collection of scattering states, labeled by momentum p. The leading contribution to this process was an interaction of the form A N → A N −3 + a p , where A N denotes an axion star with N axions and a p denotes a relativistic axion with momentum p. Without the addition of these scattering states, the matrix element for this and many other interactions are identically zero. Our analysis assumed a small binding energy in the axion star. A contrarian point of view was expressed in [40]. We found in [32] that the lifetime of an axion star through emission of relativistic axions depends on a reduced binding energy parameter ∆ ≡ 1 − (E tot /N m) 2 . The leading-order expansion in ∆ 1 is equivalent to the infrared limit of the theory, where only the marginal φ 4 term appears in the interaction potential [24]. For weakly bound stars, the leading process A N → A N −3 + a p has a rate which, as a function of ∆, is dominated , and are thus very stable in this sense. More generally, we found that if a star has ∆ .05 − .06, then it is stable on timescales as long as the age of the universe, because the lifetime is a monotonically decreasing function of ∆ in the relevant range. The constant in eq. (4.2) has the value y M = 25.46. The dense energy minimum ρ D has a large binding energy, corresponding (in the Gaussian case) to ∆ = .56. A naïve application of eq. (4.2) at this large value of ∆ gives τ = 10 −9 sec; however it is not known whether this estimate is reliable, since the analysis of [32] applies only in the weak binding limit. Further, eq. (4.2) takes into account only the attractive φ 4 interaction, but this is a valid approximation throughout most of the decay process. Nonetheless, if valid, such a short lifetime would imply that these dense states, as the endpoint of collapse, would decay very quickly. However, our calculational method is not applicable to strongly bound systems, so we cannot make a definite statement about it. We hope to investigate the decay of strongly bound states in greater detail in the future. Recent investigations of collapse using a classical collapse analysis have concluded that collapsing axion stars lose a significant fraction of their mass through emission of relativistic axions [43,44].
In the weak binding region, where eq. (4.2) holds, we know that ∆ is a one-to-one function of ρ, and thus also of the collapse time t as defined in eq. (3.20). We find that the binding energy obtains ∆ ∼ .05 at ρ ∼ 10 −4 (compared with ρ D ∼ 10 −5 ). As a function of time, the binding energy only changes appreciably in the last fraction of a second of the collapse, but rises quickly to a strongly bound final state (see Figure 7). In these last moments, the decay rate in eq. (4.1) becomes astronomically large; Γ ∼ 1 emitted axion/sec at ∆ ∼ .0223, and rises to Γ ∼ 10 50 emitted axions/sec at ∆ ∼ .1. We therefore are led to the conclusion that axion stars, as they collapse, emit many highly energetic free axions. 4 Such an explosion, referred to as a Bosenova, has been observed experimentally by condensed matter physicists using cold atoms [42]. While this work was under review, a different group performing a numerical simulation also suggested that a large fraction of axion star energy is expelled during the collapse process through relativistic axion emission [43].
We emphasize again that the analysis of the decay process in [32] applies only at weak binding, when ∆ 1. This condition holds for the dilute state as well as throughout a large portion of the collapse process, but it is possible that some new dynamics take hold at truly strong binding ∆ = O(1), as for the dense global minimum of the energy where ∆ ∼ .56. We are led to the conclusion that relativistic axion emission becomes important during collapse, but it is possible that a stable, strongly bound remnant remains.

Conclusions
The contribution of the axion self-interaction potential to the total energy in the variational method can be computed to arbitrary order using an expansion in powers of the axion field. This expansion is equivalent to an expansion in the small parameter δ = f 2 /M 2 P 1. Because of the smallness of this parameter (δ ∼ 10 −14 for QCD axions), the potential is typically truncated at leading-order, including only the attractive (ψ * ψ) 2 term. This truncation works extraordinarily well at large radii, and the dilute radius R * found by multiple authors previously [24,25] is preserved. In the regime of larger δ (e.g. axions with f ∼ .1 M P ), some of these conclusions could be changed. While this work was being reviewed, an analysis performed in the classical limit [44] suggested that axion theories with large δ indeed allow collapse to black holes in some regions of parameter space. We are working out the mass spectrum of axion stars in such theories, which will be the topic of future work.
Going beyond the leading-order approximation, without truncation we have found that a global minimum of the full energy exists, which is not present in the leading-order expansion; we calculated its position, and it corresponds to a radius R GM many orders of magnitude larger than the corresponding Schwarzschild radius. We approximate this global minimum using a next-to-leading-order expansion, using the truncated energy of eq. (3.19), which has a global minimum at a radius R D R * . For m = 10 −5 eV QCD axions and using the Gaussian ansatz, the dilute radius R * ∼ 200 km, while R D ∼ 7 meters. R D is a good order of magnitude estimate of R GM .
Previous analyses of collapsing boson stars with an attractive self-interaction have concluded (correctly) that, with nothing to stabilize the potential as R → 0, the endpoint of collapse is a black hole state. For the axion potential, we have found higher-order self-interactions, some of which are repulsive, stabilize axion stars as they collapse and there exist energetically stable configurations at very small radii. These configurations correspond to dense axion star states which are nonetheless not black holes, and resemble closely the type of the dense states found by the authors of [27] using a different method. Dense configurations of this kind can exceed the maximum mass normally allowed for weakly bound axion stars, which is roughly M c ∼ 10 19 kg for m = 10 −5 eV axions [24].
We have examined the collapse dynamically in time, and find that masses M just above M c collapse from R = R * in a time on the order of hours, or tens of minutes. The radius changes slowly at first, then drops rapidly as the slope of the potential becomes increasingly steep. Stars which begin collapse at a radius R 0 < R * were also considered, a case which is interesting if, for example, axion star collapse is catalyzed by collisions of two lighter axion stars. This could occur even if these stars do not become gravitationally bound to each other. This topic will be pursued in a future work.
If stable, then heavy axion star states could be detectable via gravitational lensing experiments. Such states have large binding energies, and thus non-relativistic and nonperturbative corrections may become important in that regime. During collapse, however, when binding energies increase but are still sufficiently small, previous calculations [32] suggest that the rate of emission of relativistic axions from an axion star will rise very rapidly. The rate of decay through the leading number-changing interaction A N → A N −3 + a p rises to Γ 10 50 emitted axions/sec in the final moments of collapse, leading to rapid emission of axions in what is often called a Bosenova [42]. It is not clear in our analysis precisely what fraction of the energy of the star would be expelled through this process, or whether a stable dense state could remain. It would be interesting to investigate the energy spectrum of these collapses in detail, to determine if there are detectable consequences of such an explosion.