Global View of Axion Stars with (Nearly) Planck-Scale Decay Constants

We show that axion stars formed from axions with nearly Planck-scale decay constants $f$ are unstable to decay, and are unlikely to have phenomenological consequences. More generally, we show how results at smaller $f$ cannot be naively extrapolated to $f=\mathcal{O}(M_P)$ as, contrary to conventional wisdom, gravity and special relativity can both become relevant in the same regime. We clarify the rate of decay by reviewing and extending previous work on oscillons and axion stars, which imply a fast decay rate even for so-called dilute states at large $f$.

Axion stars on the dilute branch are generally stable both under perturbations [35,36] and to decay [40,41]. Recently both simulations [42][43][44] and analytic arguments [45] have suggested that they can efficiently form * Electronic address:joshaeby@gmail.com † Electronic address:streetlg@mail.uc.edu ‡ Electronic address:peter.suranyi@gmail.com § Electronic address:rohana.wijewardhana@gmail.com in the early universe; as such, dilute axion stars have been investigated for possible phenomenological effects, including recent analyses of radio photon emission [46][47][48][49] and gravitational lensing [50][51][52]. For the attractive φ 4 potential, there is a maximum stable mass of order M P / |λ 4 | = M P f /m (for λ 4 = −m 2 /f 2 < 0, and with M P = 1.22 × 10 19 GeV) which signals a crossover to the structurally unstable transition branch [35,36,53,54]. Because transition states are unstable to perturbations, they are unlikely to have observable consequences. On the third branch, so-called dense axion stars have received considerable interest in recent years. The "dense" moniker was coined in [55] when the configurations were investigated using the nonrelativistic equation of motion; note however that they are fundamentally no different from axitons, discussed decades earlier in a cosmological context by Kolb and Tkachev [56]. Today, dense axion stars / axitons are understood to be strongly bound with momentum dispersion of order m, and as a consequence, they decay to relativistic axions with a short lifetime and are unlikely to have phenomenological consequences [20,21].
In this work, we investigate the nature of axion stars at large f , approaching the Planck scale. The reader may wonder whether the case of large f is well-motivated enough to warrant a full investigation; to allay this critique, we review several contexts in which this parameter space is commonly invoked: 1. The QCD axion [1][2][3][4][5][6][7][8] is motivated by its ability to solve the strong CP problem while simultaneously producing a good DM candidate. This model imposes a relation between the two relevant scales such that m f 6 × 10 −3 GeV 2 [57,58]. The abundance of dark matter QCD axions depends on cosmological assumptions, with masses m µeV required if the global symmetry of the axion is broken after inflation, corresponding to relatively low f 6 × 10 12 GeV (see [59] for a recent, more stringent bound on the post-inflationary scenario). On the other hand, if symmetry breaking occurs before inflation, then masses all the way down to m 5 × 10 −13 eV are allowed (at the expense of tuning the initial misalignment angle, see e.g. [58], or non-minimal model-building, e.g. [60]), in which case f M P is obtained. Even in such a model, where there are no large primordial density fluctuations, it may be possible to form axion stars through large misalignment [61]. 2. In other particle physics models, the scales m and f characterizing the axion are independent. ALPs may, for example, originate from string theory at low-energy from string compactification or as a consequence of anomaly cancellation [62][63][64]. In such cases the decay constant f is naturally very large, of order the Grand Unified Theory (GUT) scale Λ GU T 10 16 GeV or higher. A widely-known example is Ultralight DM (ULDM) [9-12, 14, 15, 17-19], where m 10 −22 − 10 −20 eV implies galaxyscale axion wavelengths, and axion stars forming in the central cores of galaxies [42,[65][66][67][68]]. Another well-motivated particle physics model is coherent relaxion dark matter, which may also require f Λ GU T in order to be consistent with constraints from fifth-force experiments [69].

A more phenomenologically-driven parameter
choice is f ∼ M P and m ∼ 10 −10 eV, as the dilute bound states become very compact, having (at their maximum stable mass) roughly the size and mass of neutron stars: R = O(few) km and M = O(M ). Such states may collide with other astrophysical bodies, giving rise to gravitational wave signals [70,71].
4. Axion stars have also been investigated as the origin of certain black holes in the universe. At small f , it is known that axion stars which grow unstable and collapse do not form black holes, as they are stabilized by a core repulsion; they instead explode in a burst of scalar radiation in their approach to the dense configuration, a process known as a Bosenova [20,[72][73][74]. But at large f 0.1M P (whereM P = 2.4 × 10 18 GeV is the reduced Planck mass), it has been suggested that collapse can indeed lead to formation of black holes [74][75][76][77]. If true, this process could explain the tentative recent observation of a black hole of M ∼ 2 − 3M by the LIGO/VIRGO collaborations [78], which lies in a mass gap not easily explained by standard theories of black hole formation.
On the basis of this and other work, we conclude that the parameter range Λ GUT f M P is motivated, and warrants the study we put forward here.
In this work, we show that the naive picture of three branches of axion stars outlined above breaks down as the decay constant f approaches M P , due to the usual assumptions about the relevance of gravity and special relativity breaking down. This fact implies that previous naive estimations of axion star parameters in this regime have neglected important contributions. In addition, we will point out that axion stars with large f become unstable to decay to relativistic particles, even on the dilute branch; as a result, such axion stars are unlikely to be phenomenologically relevant in the manner described above. We will make this point by reviewing calculations for the decay rate in previous literature, extending them to include gravity and higher-order self-interactions, and finally determining the full range of stable axion star solutions.
We will use natural units throughout, where = c = 1.

A. Non-Relativistic Bound States
We first review a few relevant facts about axion stars, which are derivable using a number of viable methods [39]; for our purposes, the formalism of Ruffini and Bonazzola (RB) [27] is the most useful. RB used an expansion of the axion field operator which was linear in creation and annihilation operators to describe axionic bound states when self-interactions were absent; their work was first extended to the case of an attractive φ 4 self-interaction in [79]. Some of the present authors further extended the analysis to fully characterize the attractive φ 4 case [53], and also to include contributions from scattering states that give rise to decay processes [40,41] and relativistic corrections to bound states [80].
The generic, spherically-symmetric RB field operator can be written in the form φ(t, r) = R(r) e −i m t a 0 + h.c. (RB) where R(r), R k (r), and ψ f (t, r) are the N -particle bound ground state, higher-harmonic states, and single-particle scattering state wavefunctions (respectively), a 0 and its conjugate are the annihilation and creation operators for the bound state, and m < m is the bound state eigenenergy. The scattering and higher-order harmonic modes were not included by RB, and will be discussed in the next sections. In the weak binding limit, where | − 1| 1, it is appropriate to rescale the wavefunction as Y (x) = 2 √ N R(r)/(f ∆) with the coordinate rescaled as x = ∆ m r [53], where ∆ ≡ √ 1 − 2 . Bound-state configurations can be determined by solving the coupled Einstein+Klein-Gordon (EKG) equations for the wavefunction Y (x), along with the functions determining the gravitational metric. In the Newtonian and nonrelativistic limits, the equations of motion are [53] where the spherically-symmetric metric function 1 g tt = 1 − δ b(x) with δ = 8π f 2 /M 2 P ; the other relevant metric function g rr = 1 − δ a(x) was eliminated using the rr Einstein equation. Note that we have taken the potential to be that of Eq. (1); it has been shown that in this case the above equations are equivalent to the Schrödinger+Poisson equations often used in the study of boson stars [81].
The set of equations (3)(4) represents the leading order in a double-expansion of the relativistic EKG equations in the two small parameters δ 1 (representing weak, Newtonian gravity) and ∆ 1 (representing weak binding, or the nonrelativistic limit); as written, they are correct to O(δ, ∆ 2 ) [53], and in the next section we extend to the next order in special-relativistic corrections, O(∆ 4 ). Previous perturbative studies have been restricted to either the fully nonrelativistic limit (e.g. [35] and many others), or they include leading-order special relativity but neglect Newtonian gravity [82][83][84][85]; in the next section, we attempt to bridge this gap by including both effects. For a detailed discussion of the relativistic expansion parameters, see [39].
There exist higher-order corrections proportional to δ ∆ 2 or δ 2 , which we require to be subleading by restricting the maximum parameter values δ max = 0.1 ∆ 2 max = 1/2 throughout; the effect of these subleading terms was partially investigated in [86], but for simplicity we leave a full study for future work.  The one free parameter in Eqs. (3)(4) is κ = δ/∆ 2 , which is the effective coupling to gravity. For a given value of κ, these coupled equations have a unique solution that is easily determined by the standard shooting method: one varies the central values This can be done to arbitrarily high precision (see [53] for details). For clarity we translate between values of δ, f , and κ in Table I for the range we investigate here.
Eqs. (3)(4) are, in most applications, sufficient to describe the dilute and transition branches of axion star solutions. On the stable dilute branch, the total mass M ∝ ∆ and the radius R ∝ ∆ −1 , so that M ∝ R −1 . Stable solutions have κ ≥ κ c ≡ 0.34; the inequality is saturated at a maximum stable mass M c = 10.2 M P f /m, which corresponds to a critical value of Beyond this critical point, the central density Y (0) grows while the gravitational coupling κ continues to decrease below unity, leading to a decoupling of gravity which characterizes the transition branch. From here we can already see that naive extrapolation to f = O(M P ) will fail, as Eq. (5) predicts ∆ c = O(1) or larger. In reality, we will find in the next section that above some value of f this naive maximum mass point is no longer attainable. Indeed, when ∆ becomes large, the eigenenergy in the axion star approaches m, implying a breakdown of the nonrelativistic approximation [39]. This also affects the bound state and leads to relativistic decay processes that render the star unstable; we discuss both effects below.

B. Bound States with Relativistic Corrections
There have recently been several independent efforts to quantify the effect of relativistic corrections on axion stars, which are relevant at large densities [82][83][84][85]. In one such work, it was shown how such corrections can be organized as a power series in the parameter ∆ by extending the RB framework, thereby defining a Generalized Ruffini-Bonazzola (GRB) procedure [80]. In that work gravity was negligible, because (a) the focus was on the transition / dense branches of solutions, and (b) the input parameters were assumed to be in the standard range for QCD axion experiments, where f ∼ 10 10−12 GeV. However at large f , the parameters δ and ∆ can be of the same order, and therefore both Newtonian gravity and special-relativistic corrections must be taken into account.
At leading order in higher-harmonic (GRB) corrections of Eq. (2), Eq. (3) is modified as [80] Table I; the blue curve corresponds to the limit of gravity decoupling, κ → 0. Configurations on the solid curves are stable to decay processes with long lifetimes τ0 > τU ; those along the dashed curves are unstable to decay. The filled circles represent the endpoint of each curve at ∆ = 1/ √ 2.
where we continue to use the potential of Eq. (1). Once again, note that we ignore corrections proportional to δ ∆ 2 , which appeared in [86], or proportional to δ 2 , which are post-Newtonian, because such corrections are small for the parameter space we consider. The mass M of the axion star at leading-order in GRB We define the radius as R 99 , inside which 0.99 of the mass is contained, To the extent that the number of particles is conserveda topic we will return to in the next section-we can treat axion stars as N -particle states with N = M/m. The system of equations (4 , 6) requires two input parameters, ∆ and κ, though at fixed ∆ we can trade κ for the more intuitive parameter δ. Then the solutions 2 For more information, see [80] as well as Appendix B of [87]. have masses and radii computable using Eqs. (7)(8); the results are depicted in Fig. 1, where the curves are labeled by values of δ and color-coded as in Table I. The blue line corresponds to κ → 0 or full gravitational decoupling. The endpoints of each curve (marked with a filled circle) represent ∆ = ∆ max = 1/ √ 2, essentially an arbitrary cutoff; at ∆ max , not only does the GRB expansion break down, but the binding energy is so large that the momentum uncertainty becomes comparable to the mass m, and the usual methods of analysis fail [39]. Note for future reference that on the red curve, f 2.4 Λ GUT , and on the purple curve f 0.3M P .
For low f 10 17 GeV, we see in Fig. 1 the standard separation into three distinct branches of solutions. The dilute and transition branches are joined at the (local) maximum mass M max , while the transition and dense branches are joined at a (local) minimum mass M min ∼ few × 100f 2 /m. As f grows, we find at this order in ∆ that the transition branch shrinks, while simultaneously the cutoff at high ∆ limits the extent of the dense branch, leaving only the dilute branch intact. It is possible that between our cutoff ∆ max = 1/ √ 2 and the maximum possible value ∆ = 1 that some portion of the dense branch persists; however, the transition branch is almost certainly lost, as we see the dilute and dense branches joining to a single point between δ = 10 −1 and δ = 10 −2 .
In previous literature, it has always been assumed that the region where relativistic corrections are important is nonoverlapping with the region where gravity is important; usually the dense and dilute branches are separated by a large transition parameter space. However, we see in Fig. 1 that this assumption is badly violated as f grows above Λ GUT . At δ = 10 −1 (purple curve), the high-∆ cutoff occurs very close to M max ; we were unable to push to larger δ in this analysis due to the appearance of higher-order corrections, but it is plausible that for δ = O(1) the naive maximum mass is never reached. Other aspects of the solutions in Fig. 1 will be discussed in the next section.
An important consequence of Eq. (6) for large f is that due to shifts in the parameter ∆ as f increases, decay processes can be relevant even on the dilute and transition branches. We address this in the next section.

III. DECAY RATE
Axion stars decay through emission of relativistic particles. There are two classes of decay modes, which proceed through tree-level [40,41,82] and loop-level [88] diagrams 3 . Importantly, the momentum distribution of bound axions has a finite width, which allows for tree-level decay processes that would be forbidden by energy-momentum conservation for particles in momentum eigenstates [40,41]. For a φ n potential, the leading process of this type is (n − 1) a c → a f , where n − 1 bound ("condensed") axions a c annihilate to a single relativistic axion a f emitted from the star. We focus below on this tree-level decay process, though of course our conclusions would only be made stronger if higher-order processes were included as well.
In the specific case of the axion potential in Eq. (1), the leading self-interaction potential at O(φ 4 ) gives rise to a 3 a c → a f decay process. One can approximate the decay rate Γ 3 for this process by taking the matrix element of the self-interaction potential, between the initial state N a c | and final state |(N − 3) a c + 1a f . The outgoing particle is emitted as a spherical wave, conserving average momentum for the [remaining condensate + free particle] system (see [41] for a detailed discussion).
Previous work on relativistic decay processes in axion stars [21,40,41,73] have assumed the decoupling of gravity in the region where decay becomes relevant. However, we find in this work that gravity does not decouple in this way when f becomes large, and thus we repeat the calculation of the lifetime in the Appendix A, including both effects.
A second assumption made in previous work is that the axion star can track its equilibrium configuration as it decays, which we refer to as the adiabatic approximation. However, simulation results [21,73], previous semianalytic estimations [40,41], as well as the present work (see Appendix A), all suggest that when decay becomes relevant, the resulting explosion of relativistic particles is extremely rapid, likely to greatly outpace the relaxation of the star back to equilibrium. Therefore in this work, we argue that an instantaneous approximation for the lifetime of axion stars is more appropriate. We thus define the lifetime as following Eq. (4.1) of [40] evaluated in the instantaneous limit 4 . This lifetime can be uniquely determined for each configuration (on any branch of solutions), which we compare to the age of the universe τ U ≈ 13.8 × 10 9 years. As described in the Appendix A, we find that there is a sharp transition from τ 0 τ U (stability) to τ 0 τ U (instability) over a narrow range of M , which we define M τ (δ). In Fig. 1, we mark M τ by the change from solid (stable, τ > τ U ) to dashed (unstable, τ < τ U ) lines on each curve. Importantly, at large enough f , the critical point M τ occurs on the dilute branch, which has until now been regarded as perfectly stable.
The precise value of M τ depends exponentially on ∆ but only polynomially on m. Owing to this functional dependence, the transition from τ < τ U to τ > τ U occurs at a very sharply defined M τ , varying only by at most a factor of 2 as the axion mass m varies between 10 −20 − 1 eV.
We further illustrate in Fig. 2 how the decay instability (blue shaded region) sets in on the dilute branch, at lower mass than the onset of gravitational instability (black shaded region), when δ 10 −3 (f 8 × 10 16 GeV). The limit of very large binding energy, where ∆ = 1/ √ 2, where our approximation is no longer valid, is illustrated by the green shaded region. The blue band represents M τ (δ), and its width represents its variation upon varying m in the range 10 −20 − 1 eV. The red band similarly represents the onset of decay instability on the transition branch of solutions for smaller δ.

IV. DISCUSSION
In this work, we have analyzed the structure and stability of axion stars when the decay constant f is large, with a focus on the approach to the Planck scale M P . Such scenarios are motivated by a significant body of literature, and previous work has typically assumed that results at smaller f can be extrapolated to large f . We have shown that this extrapolation is erroneous; at the order of ∆ that we consider, the dense and transition branches of solutions no longer exist above δ 10 −1 (f 8 × 10 17 GeV) up to very large ∆ = 1/ √ 2, and the decay instability point M τ occurs on the dilute branch for δ 10 −3 (f 8 × 10 16 GeV), as shown in Fig. 1.
This result has applications in any scenario where axion stars form with f 10 17 GeV, including those outlined in the introduction. For example, for ULDM with particle mass m ∼ 10 −22 eV, the correct relic abundance is obtained for f ∼ 10 17 GeV [18]. ULDM simulations generally find axion stars forming the cores of galaxies in this scenario [42,[65][66][67][68], with masses that are safely below the instability points we find here; however, if the axion star masses had been a factor of ∼ 80 larger at formation, or if they can accrete mass efficiently, then ULDM axion stars would not merely collapse but also decay on the dilute branch, strengthening the argument of previous studies [92] (see Appendix B for details). Such an effect would not be seen in standard ULDM simulations, as they typically neglect both the self-interaction potential as well as relativistic effects. As simulations of axion star formation and accretion become more precise, or as nonminimal axion models are investigated, this must be taken into account in the final analysis.
Further, phenomenological studies of axion stars (including those looking for gravitational wave signals [70,71] or new black hole formation mechanisms [74][75][76][77]) must acknowledge that decay may make their proposed configurations unstable. For f 7 × 10 17 GeV, the mass of truly stable axion stars is already a factor of ∼ 3 lower than the usual boundary of gravitational stability, as shown in Fig. 2. This in fact precludes axion stars with neutron star-like masses and radii by nearly an order of magnitude.
We should point out, of course, that our own results should not be naively extrapolated either. For example, in the limit f → ∞ the axion self-interactions will decouple, and in that limit neither nonrelativistic collapse nor relativistic decay processes discussed here destabilize the star. Taking f > M P may, however, be in tension with theoretical considerations like the Weak Gravity Conjecture (see e.g. [93]), though there may be ways to reconcile them [94,95].
Previous works have also considered decay modes other than 3 → 1, including so-called quantum decay (e.g. 4 → 2) [89], other classical decay processes (e.g. 5 → 1) for alternate axion potentials [96], or for repulsive selfinteractions [91]. In some contexts it is claimed that the resulting decay modes can actually dominate the decay rate. Because we have not included such contributions in this work, we emphasize that the lifetime we calculate is merely an upper bound on the true axion star lifetime.
We found in our study that higher-order corrections to the EKG equations appear not only at higher orders in ∆ 2 (as found in [80]), but also as powers of the product δ ∆ 2 and as δ 2 ; as a result we were not able to robustly analyze configurations with δ 0.1, corresponding to f 8×10 17 GeV. We leave the effect of these corrections, and the (in)stability of axion stars with δ = O(1), for a future study, It would also be worthwhile to make a thorough comparison of our results, which apply to static axionic configurations, to the dynamical works of [74,77] which use very different initial conditions; see Appendix B for a more heuristic comparison.

Note Added
After completion of this project, a related work [97] appeared, which illustrates a different method but comes to many of the same conclusions. We believe the two papers complement one another. To analyze decay processes, in previous work some of us showed how to extend the RB field operator to include a scattering state contribution in Eq. (2), allowing axion quanta with energy ω p > m [40,41]. The scattering state wavefunction ψ f , at leading order in spherical harmonics, is given by where a 00 (p) is the annihilation operator for a scattering state of momentum p labeled by its angular momentum quantum numbers = z = 0, and j 0 is the zeroth spherical Bessel function. The decay rate was then analyzed for an attractive φ 4 potential, where the rate is always nonzero due to an essential singularity in the equation of motion; however, the matrix element is exponentially suppressed for weakly-bound axion stars. In the range of parameters considered there (for example, f ∼ 10 10−12 GeV), the decay process was irrelevant on the dilute branch and became important on the transition branch. The lifetime of an axion star on the transition branch was found to be where the constant y t was determined by fitting the curve M (∆) on the transition branch, and the constant y I characterized the position of the singularity in the complex plane at x = i y I The lifetime in Eq. (A2) is approximately correct in its region of applicability, but rests on assumptions about the input parameters. For example, in [40,41] it was assumed that the effect of gravity could be neglected (using the κ → 0 limit of Eq. (3)), which becomes appropriate near M ≈ M max and remains so on the transition branch. Then the Klein-Gordon equation for Y (x) has a unique solution corresponding to Y (0) = 12.268, and the constant value for y I ≈ 0.602 is uniquely determined. Because gravity does not decouple in this way at large f (see Main Text), we determine the position of the singularity for nonzero κ below. We also show how the singularity is shifted to larger values of y I at leading-order in GRB, leading to slower decay rates than one would have obtained using the constant y I ≈ 0.6. The calculation is detailed below, and the result is given in Fig. A 1.
The decay rate depends upon the integral [40] ("Bessel"), (A3) which, at small enough ∆, can be approximated by 8 is the momentum of the outgoing relativistic axion for a 3 a c → a f annihilation, in units of the axion mass. These integrals, labeled for future convenience, can be calculated directly by numerical integration only at relatively large ∆, as the integrand is highly oscillatory at small ∆.
A simpler approach is available, as the integration is dominated by the leading singularity of the wavefunction Y (x) in the complex plane. To determine the position of this singularity, we follow the prescription of [40] and use the following ansatz for the wavefunction in the vicinity of the singularity: When κ > 0, the gravitational potential b(x) shifts the position of the singularity, and must be determined selfconsistently with the wavefunction.
To determine the effect of gravity on the singularity, observe that under spherical symmetry, we can rewrite Eq. (4) as Then substituting the ansatz of Eq. (A5) for Y (x) we can perform the integrals explicitly; we find Of course, the combination Y s and b s satisfies the Poisson equation [81] ∇ 2 Thus we can see that the solution for b(x) near the singularity at x = i y I is regular. We now Taylor expand the wavefunction and gravitational potential around x = 0, as Evaluating the equations of motion (3-4), we obtain recursion relations among the expansion parameters η n and β n ; for example, the coefficient multiplying x 2n in Eq. (3) must evaluate to zero, which implies (2n + 2)(2n + 1)η n+1 + 4(n + 1)η n+1 and similarly for β n and Eq. (4). We therefore leave η 0 and β 0 as free parameters and iteratively solve for η n>0 and β n>0 using these recursion relations, truncating at some large maximum integer in the expansion. The firstorder terms in the expansions imply the coefficients whereas the second order implies and so on. We can then match Eq. (A5) and Eq. (A9) to obtain relations for A and y I in terms of the parameters η n ; at nth order in the expansion, we obtain [40] A = (−1) n η n y 2n+2 Noting that using the recursion relations above, we have η n = η n (η 0 , β 0 , κ), we can take a given solution and evaluate the position of the singularity y I (η 0 , β 0 , κ). In the limit κ → 0 (gravitational decoupling), we recover the result of [40] that A = 8 y I and y I = 0.602. More generally at non-zero κ, we find A = 8 y I still holds, and illustrate the position of the singularity y I in the left panel of Figure A 1.

Position of the Singularity in GRB
We turn now to the structure of the essential singularity including relativistic corrections, using the Generalized Ruffini-Bonazzola (GRB) formalism. The equation of motion for Y (x) at next-to-leading order in the GRB formalism is Eq. (6). Using the recursion procedure above, we can determine the singularity structure of this equation as well. There are small modifications to the coefficient relations, e.g. η 1 in Eq. (A11) is modified at but the basic procedure is the same as above.
The GRB equation of motion depends on an additional parameter ∆ in addition to κ, though at fixed ∆ we trade the latter for δ = ∆ 2 κ. Therefore we characterize our solutions by the input parameters {η 0 , β 0 , ∆, δ}. We solved the GRB equation over a large range inside the bounds 10 −4 ≤ ∆ ≤ 1/ √ 2 and 10 −6 ≤ δ ≤ 10 −1 (this corresponds to 10 −4 κ 10 3 ); at the largest values of ∆ we expect the higher-order contributions in GRB to be very relevant, and worse yet, such solutions are unphysical due to extremely high binding energies [39]. Still, we can analyze the structure of solutions as an academic exercise, and we will see that solutions with very large values of ∆ are not phenomenologically relevant anyway due to fast decay rates.
In the right panel of Figure A 1, we illustrate the shift in the singularity position for different choices of δ, with the goal of approaching f = O(M P ). We see that, as expected, when f is sufficiently small (given by the blue curve, or f 2.4 × 10 15 GeV), the singularity is unaffected by the GRB correction over the full range of κ we analyzed. As f increases, the deviation in the singularity point appears at larger values of κ, though at the same time the cutoff at ∆ = 1/ √ 2 (given by the circles in the Figure) reduces the relevant physical range.
We can now proceed to evaluate the lifetime of axion stars at large f and moderate ∆, at leading-order in GRB. With the position of the singularity y I in hand, we can use the residue theorem to compute I 3 in Eq. (A4) and find (A14) which depends exponentially both on the inverse of the binding energy parameter 1/∆ and the position of the singularity y I .
However, note that at large ∆ the result in Eq (A14) will not be correct, as the application of the residue theorem assumed the validity of Eq. (A4). When the product ∆ Y (x) is large, we must not expand the Bessel function J 3 (∆ Y (x)) in Eq. (A3); fortunately, it is in this range that the integrand does not oscillate fast and we can integrate I 3 directly. To briefly summarize: • The analytic "Residue" result of Eq. (A14) is applicable at low ∆ 0.1; • The numerical result "Bessel" from integrating J 3 in Eq. (A3) is applicable at the largest ∆ ∼ ∆ max = 1/ √ 2; • Both should roughly agree, with each other and with Eq. (A4) ("Y 3 "), in an intermediate range We illustrate the results of these estimations (blue for residues, green for Bessel, red for Y 3 ) in Figure A 2. The black line, which we use to calculate the decay rate in the next section, interpolates between Eq. (A14) at low ∆ and Eq. (A3) at high ∆; when in doubt we used the smaller estimation of |I 3 | to make a conservative estimation of the decay rate.

Decay Rate and Lifetime
The decay rate for a single annihilation process is given by [40,41] We show the resulting rate of annihilations/sec in We observe a very rapid increase of the decay rate, from extremely small values 1 annihilation/sec to greater than 10 40 annihilations/sec around (for example) ∆ 0.03 for δ = 10 −4 . This increase occurs at slightly higher values of ∆ for larger m or for larger f (i.e. larger δ).
The lifetime is proportional to |I 3 | −2 ∝ exp(y I /∆), implying that dilute axion stars (corresponding to the smallest allowed ∆) may be stable and survive longer than the age of the universe. However, this has only been explicitly investigated at small values of f , and we have seen that naive extrapolation of such results to large f may not be appropriate. Indeed, we find that contrary to conventional wisdom about axion stars, the lifetime for f 10 17 GeV is shorter than the age of the universe even on the dilute branch, as explained below.
We estimate the lifetime numerically as follows. In a region of parameter space where the decay rate is large, we may approximate it as constant, as the axion star explodes in a rapid Bosenova of relativistic particles; we call this the instantaneous approximation. In that case (assuming a constant decay rate), the timescale for the evaporation of the star would be where the prefactor N/3 is the number of annihilations necessary to completely deplete the star (for each annihilation, 3 bound axions are lost). Then, because M has a one-to-one relationship with ∆, we can compute the decay timescale uniquely as a function of ∆. In practice we use the full solution for M (∆) in the numerical results, which is illustrated in Fig. 6. Note that the lifetime in Eq. (A16) is identical to the one used in [40,41], but here it is evaluated in the instantaneous limit; we clarify the difference at the end of this section. Using fits on each branch of axion stars, we obtain an analytic form for the lifetime which may help guide the reader's intuition. Firstly, in [40,41], we noted that (at small f 10 12 GeV) the decay rate was exponentially small on the dilute branch and that the lifetime only became short on the transition branch; in that case, the relationship between M and ∆ is given by with y t 75.4 is determined by fitting the mass function shown in Fig. 6 [41]. Then the timescale for instantaneous decay is For large f , we find in this work that the decay timescale can, on the dilute branch, already be lower than the age of the universe, as ∆ becomes large before the crossover point. In that case the mass function is 6), and the lifetime takes the form (A20) Due to the rapid turn-on of the decay rate, aside from a very narrow region near τ 0 τ U , the lifetime of axion stars is either (i) so long that the star can be treated as stable with a conserved particle number N , or (ii) so short that the star decays almost instantly. In the latter case, the instantaneous approximation above seems very well-justified, as the relaxation time for the axion star to remain in its equilibrium configuration at each instant as it evolves will likely be much longer than the decay timescale τ 0 . Therefore, we use the instantaneous case to define the crossover from stability to instability under decay which we describe in the Main Text. The result for the instantaneous lifetime as a function of M and δ is given in Fig. 7  other hand, the instantaneous case of Eq. (A18) gives τ ∼ a I (∆/m) exp(c/∆), where a I is another constant. The latter is precisely the scaling found in previous investigations using classical field theory; see e.g. Eq. (24) of [89] or Eq. (50) of [90]. In all cases, the prefactor retains no explicit dependence on f .
so in the largest halos simulated, where M h 10 12 M , a typical ULDM candidate with f = 10 17 GeV will form an axion star core with mass ∼ 80 times below its maximum mass. If axion stars can accrete enough mass after formation to cover this gap, they run the risk of not only collapsing but also decaying, because as illustrated in Fig.  1 the decay instability sets in also in the same parameter range.
Collapse Simulations: Other simulations have probed the fate of axion stars at large f = O(M P ) by evolving the classical equations of motion [74,77]. One of the purposes of these simulations was to determine under which conditions an axion star might collapse directly to a black hole; their analysis suggests a triple point in the axion star phase diagram at f T P 0.3M P and M T P 2.4M 2 P /m, which is δ T P = 0.09 and M T P = (2.4/δ)(f 2 /m) = 26.6 f 2 /m in the notation of this paper. We confirm some of these results but not all, as explained below.
The "dispersal region" in the phase diagram of [74], which we understand to be the transition branch of solutions in the equations of motion, disappears above the triple point at f 0.3M p . We observe this as the disappearance of the transition branch, which occurs for δ 10 −1 (the purple curve in Fig. 1), roughly confirming their result for f T P ; we do not, however, confirm the mass at the triple point, which is a factor of few larger in their result compared to ours. It is possible this difference originates in the use of the "classical" equations of motion in [74] (i.e. a cosine rather than a Bessel function in the self-interaction potential); we have pointed out previous that this change, when evaluated on the transition branch of axion stars, leads to differences compared to the semi-classical analysis we have used here [39].
Secondly, in spite of analyzing states at large binding energy, we find General Relativistic effects to be negligible everywhere, and therefore do not confirm the formation of black holes observed in [74,77]. For each solution at large δ, we confirm first that the Schwarzschild radius R S = 2 G M is always much smaller than the radius of the star R 99 . We further checked that the condition R S (r) = 2 G M(r) < r is satisfied inside the star at every r < R 99 , as most of the density is in the inner region. For δ = 0.1 (near the triple point of [74]), we see in Fig. 8 that the ratio R S (r)/r 1 for all points in our solution space, and therefore these states are not black holes. Note however, that the analyses of [74,77] were dynamical, focusing on the collapse of axionic objects, whereas ours is static by construction; the results of this work apply to structurally stable (or metastable) configurations only. In those works the initial profiles were also The ratio of the Schwarzschild radius RS(r) to the radius in the axion star r, for radii less than or equal to the total radius R99, at δ = 10 −1 and at different (labeled) values of ∆. When the ratio is 1, the star is far from the General Relativistic limit.
far from the physical axion star profiles we analyze here. These factors may account for any discrepancy between the two results, though a more thorough investigation may be warranted.