Decay of Boson Stars with Application to Glueballs and Other Real Scalars

One of the most interesting candidates for dark matter are massive real scalar particles. A well-motivated example is from a pure Yang-Mills hidden sector, which locks up into glueballs in the early universe. The lightest glueball states are scalar particles and can act as a form of bosonic dark matter. If self-interactions are repulsive this can potentially lead to very massive boson stars, where the inward gravitational force is balanced by the repulsive self-interaction. This can also arise from elementary real scalars with a regular potential. In the literature it has been claimed that this allows for astrophysically significant boson stars with high compactness, which could undergo binary mergers and generate detectable gravitational waves. Here we show that previous analyses did not take into proper account $3 \to 2$ and $4 \to 2$ quantum mechanical annihilation processes in the core of the star, while other work miscalculated the $3 \to 1$ process. In this work, we compute the annihilation rates, finding that massive stars will rapidly decay from the $3 \to 2$ or $4 \to 2$ processes (while the $3 \to 1$ process is typically small). Using the Einstein-Klein-Gordon equations, we also estimate the binding energy of these stars, showing that even the densest stars do not have quite enough binding energy to prevent annihilations. For such boson stars to live for the current age of the universe and to be consistent with bounds on dark matter scattering in galaxies, we find the following upper bound on their mass for $O(1)$ self-interaction couplings: $M_*<10^{-18} M_{sun}$ when $3 \to 2$ processes are allowed and $M_*<10^{-11} M_{sun}$ when only $4 \to 2$ processes are allowed. We also estimate destabilization from parametric resonance which can considerably constrain the phase space further. Furthermore, such stars are required to have very small compactness to be long lived.

Perhaps the best motivation for physics beyond the Standard Model is the presence of dark matter which comprises most of the mass of the universe. The most exciting aspect of dark matter is that it may represent an entirely new sector of physics. Due to the current lack of discovery of any dark matter candidates that have direct couplings to the Standard Model, including WIMPs, it raises the possibility that dark matter may be part of some hidden sector and/or associated with new very heavy particles (e.g., see Refs. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]).
Hidden sectors may include pure Yang-Mills interactions, i.e., new collections of massless spin 1 particles with self interactions. Such constructions are entirely plausible from the point of view of fundamental physics. They are also considered to be entirely natural, since they do not appeal to any unnecessarily small parameters. In particular, they have no allowed elementary masses, and are fully described by only two quantities: the scale at which the theory becomes strong coupled Λ (which can be naturally small compared to any unification scale due to the logarithmically slow running of coupling) and the size of the gauge group, such as SU (N ). It is also possible of course that physics beyond the Standard Model includes various other kinds of particles. Importantly, this may include elementary spin 0 particles. Both the hidden spin 1 and spin 0 particles are of course bosons, and as such they allow for a rich phenomenology; under some conditions they can organize into interesting states of high occupancy, as we will discuss in this paper.
Let us address in more detail the case of hidden Yang-Mills. If there are indeed such sectors of interacting spin 1 particles, it leads to various questions: What happens below the confinement scale? If the particles are thermally produced in the early universe, what is their energy density today? Do the strong interactions lead to inconsistency with constraints on dark matter scattering in galaxies? And, very importantly, are there novel observational signatures of these sectors?
To address these issues, let us start in the very early universe. We know that at high temperatures of the early universe, the particles exhibit asymptotic freedom, and so we have a gas of almost free particles. When the hidden sector temperature is on order of the strong coupling scale ∼ Λ the particles are expected to lock up into color neutral states, known as "glueballs" (for a review, see Ref. [17]). The lightest glueball states are expected to be spin 0 particles, whose mass m φ is of the order the strong coupling scale Λ. As the temperature of the dark sector T d lowers further, these glueballs become non-relativistic, and can act as a form of dark matter in the late universe. If one just focusses on freeze-out, the relic abundance can be estimated to be roughly [7] Ω φ ∼ N 2 ξ 3 (Λ/(100 eV)) where ξ ≡ T d /T is the ratio of temperatures of the dark to visible sectors. If ξ ∼ 1 then one evidently needs rather small values of Λ, no larger than a few eV, to avoid overproduction. However, there are important regimes in which the abundance can be modified due to self-interactions leading to 3 → 2 processes (this will be relevant to the rates we compute in stars later). For aspects of glueball and self-interacting dark matter relic abundances, see Refs. [18][19][20][21][22][23][24][25].
On the other hand, if the temperature of the hidden sector is small (ξ 1) then one can have larger strong coupling scales without overclosing the universe. As some of us showed in Ref. [26], there is a very reasonable scenario in which the inflaton ϕ decays predominantly to the Standard Model through the dimension 3 coupling to the Higgs ϕ H † H, while its decays to hidden sector Yang-Mills is suppressed as it would occur through the dimension 5 coupling ϕ G µν G µν . This leads to the expectation ξ 1 (such as ξ ≈ 0.006 for reasonable parameters considered in Ref. [26]).
The scattering cross section to mass ratio is expected to be on the order σ φ /m φ ∼ few/(Λ 3 N 4 ) [27].
While constraints from bullet cluster implies σ φ /m φ cm 2 /gram ≈ 1/(60 MeV) 3 [28,29]. Thus to satisfy this bound, one needs Λ Λ bc where Λ bc ∼ 100 MeV/N 4/3 . This in turn implies that Ω φ is large unless ξ 1, which is compatible with the reasoning in Ref. [26] (although for extremely high N , one would need to compensate with extremely small ξ, which seems less plausible). Or alternatively, if 3 → 2 processes are significant, this can reduce the abundance further.
If we move beyond the glueball motivation, we can in fact consider any dark matter candidate which organizes into massive scalars with self-interactions. As an example, we could just study an elementary scalar with a renormalizable potential. Generally, these are all interesting candidates that are the focus of this study.

A. Novel Signature
An interesting possibility is that the dark matter scalars organize into gravitationally bound systems, socalled "boson stars" (for a review see Refs. [30][31][32][33]). For both glueballs and for elementary scalars with a regular potential, one anticipates self-interactions, including ∼ λ 4 φ 4 , etc. For repulsive self-interactions this can give rise to very dense boson stars. This was carefully studied originally in Ref. [34] for a complex scalar field with a global U (1) symmetry (recent work includes Ref. [35]). They solved the full Einstein-Klein-Gordon equations of motion for spherically symmetric time independent solutions, finding that the maximum star mass is (1) 19 GeV). At this maximum mass, the physical radius is only a factor of ∼ 2 larger than the corresponding Schwarzschild radius. For masses above M max , boson stars solutions do not exist; the density is so high that such configurations can collapse to a black hole. Interestingly, if λ 4 = O(1), then this scaling is similar to the Chandrasekhar mass for white dwarf stars M wd ∼ M 3 Pl /m 2 p ; one is effectively replacing a repulsion from Pauli exclusion of fermions by a self-interaction repulsion of bosons. Then if m φ GeV, this can be of the order of a solar mass M sun ∼ 10 57 GeV, or larger.
In the literature, this result has been applied to boson stars from real scalars, including glueballs and elementary scalars. It is not clear that glueballs have the required repulsive interaction (and in the case of the SU (2) gauge group, some recent lattice calculations suggest it may in fact be attractive [36,37]). Since we do not know the sign for a generic gauge group, we will simply work under the assumption that it may be repulsive for some cases, and this is of course certainly possible for scalars in other kinds of theories. To apply the result to glueballs, one estimates the quartic coupling as λ 4 ∼ (4π) 2 /N 2 and m φ ∼ few Λ. Then using the above bound on Λ Λ bc in order to satisfy bullet bluster constraints, one finds that the maximum mass can be as large as M max ∼ M 3 Pl N 5/3 /(100 MeV) 2 ∼ 100M sun N 5/3 . Hence one can readily have boson star masses that are larger than a solar mass by choosing the strong coupling scale accordingly. Alternatively, if one is simply studying elementary scalars, one can just impose the appropriate values of m φ and λ 4 to obtain such massive stars and satisfying bullet cluster constraints.
Very interestingly, Refs. [27,[38][39][40] pointed out that if one has M max O(10)M sun and if these boson stars undergo a merger, they can emit gravitational waves with a frequency and amplitude that could possibly be detectable at LIGO, and even heavier stars at LISA, with a signal that could be potentially distinguished from that of black hole mergers (for a focus on complex scalars, see Refs. [41,42]). This presents a novel signature of these hidden sectors. This is such an exciting possibility, it acted as a motivation for the present work.

B. Outline of this work
In this work, we critically examine whether such stars are in fact long lived. While a complex scalar with an internal U (1) symmetry, of the sort studied in Ref. [34], carries a conserved particle number, the same is not true for a real scalar. In this case there is no distinction between the particle and its antiparticle. Such real scalars can undergo quantum annihilation processes in the core of the star, including 3φ → 2φ and 4φ → 2φ processes. Such processes will cause the star to evaporate away. In this paper we compute these annihilation rates for the boson stars, finding that while the processes are negligible for very low mass stars, they are extremely important for high mass stars including those above a solar mass. We find that this leads to short lifetimes, unless the couplings are very small (associated with huge gauge groups for glueballs). However, even for small couplings, the stars whose mass would be relevant to LIGO/LISA typically have low compactness and so the gravitational wave emission would be suppressed (unless one takes extreme parameters). Furthermore, we estimate possible decays from parametric resonance, finding that this cuts down the available parameter space considerably further.
One might hope that the star carries enough binding energy to prevent these number changing processes. We find that even though the most massive stars have appreciable binding energy, it is not sufficient to prevent the radiation. Finally, for completeness, we also compute 3φ → 1φ processes, which had previously been claimed to be the most important process in the work of Ref. [43] (in the context of axions). We show that the estimates provided in that work are missing an important factor and we provide the correct scaling, including the recognition that such a process is actually provided by classical field theory. For the massive stars of interest, these 3φ → 1φ classical processes are found to typically be small, while the 4φ → 2φ or 3φ → 2φ quantum processes that we focus on can be very important.
Our paper is organized as follows: In Section II we present the basic effective field theory field theory, the corresponding classical equations of motion for a spher-ically symmetric star, and recap the properties of stars in a single harmonic approximation. In Section III we compute the quantum annihilation rates in the core of the stars. In Section IV we use these results to derive bounds on the mass and compactness. In Section V we also consider parametric resonance. In Section VI we estimate whether the star's binding energy can prevent decays. In Section VII we compute the classical radiation and compare to existing claims in the literature. In Section VIII we briefly discuss decays in boson stars supported by quantum pressure. Finally, in Section IX we conclude.

II. BOSON STARS
When glueballs form they organize into spin 0 scalar particles. In the effective field theory formalism, we can organize this into a scalar field φ. These scalars interact directly with one another in the effective theory, which is a consequence of the microscopic hidden gluon interactions. In this work, the most important form of the interactions will be the leading order operators as organized into a scalar potential V . We can expand this as where λ n are dimensionless couplings. There could be a tower of higher dimension operators, including higher order derivative operators in the effective theory, but these leading terms will be sufficient to describe the most relevant effects in this work. This formalism is obviously also applicable to other scalars, including elementary scalars with a renormalizable potential. So our analysis will be quite general. In the case of glueballs formed from the group SU (N ), the couplings are expected to scale as λ n ∼ (4π/N ) n−2 , with n = 3, 4, . . .. For a generic real scalar the same hierarchy of couplings may occur too But we do not suppose that the prefactors are necessarily close to 1. In the case of a real scalar, one can also imagine that it is endowed with some discrete Z 2 symmetry φ → −φ, forbidding all the odd powers. We shall consider this case in this work also. Of course the field is also coupled to gravity with action (we use units c = = 1 and signature − + ++) where R is the Ricci scalar, G is Newton's gravitational constant, and g µν is the metric of space-time.

A. Spherical Symmetry
Under some conditions, the scalars can condense into states of high occupancy, which is usually well described by classical field theory (we shall return to quantum corrections in the next section). A condensate for a system with gravity is a localized system; either a gravitationally bound "boson star", or a potentially bound (for attractive interactions) "oscillon". This work will focus on repulsive interactions mainly (though some of our reasoning will be relevant for attractive interactions too), so the only relevant solution is the gravitationally bound boson star. The lowest energy configurations are expected to be spherically symmetric.
In this case, it is useful to go to spherical co-ordinates (r, θ, ϕ), where quantities only depend on radius r (as well as time t). Any spherically symmetric space-time metric can be written in the Schwarzschild co-ordinates where g rr = A and g tt = B are functions that can depend on radius and time that we need to solve for; we also need to solve for the scalar field φ(r, t).

B. Einstein-Klein-Gordon Equations
The Einstein equations G µν = 8πG T µν follow from varying the action S above with respect to g µν and extremizing. Here G µν is the Einstein tensor and T µν is the energy-momentum tensor provided by the scalar field φ. Due to spherical symmetry, we only need to report on the time and radius components of the Einstein equations. The Einstein tensor can be shown to be Here we use a notation in which a prime means a derivative with respect to radius r and a dot means a derivative with respect to time t. The energy momentum tensor arises from varying the scalar field contribution of S with respect to the metric The time and radius components are readily found to be Combining eqs. (6,7) with eqs. (9,10) we have a pair of Einstein equations for A and B. Note that these equations only involve spatial derivatives for A and B, which reflects the fact that these are constrained variables; there are no dynamical (tensor) components of the metric due to spherical symmetry. Finally we need the equation of motion for the scalar field φ. This comes from varying the action S with respect to φ, giving Substituting in the above metric and carrying out the derivatives leads to There is also a 3rd Einstein equation from the time-space components, which arises from G tr =Ȧ/(rA) and T tr = φ φ . However, this is redundant with the above set of equations and so is not needed.

C. Single Harmonic Approximation
In general a solution for a real scalar field will have a complicated time dependence. A boson star φ * or oscillon is a periodic solution of the classical equations of motion. It can therefore be expanded in a harmonic expansion (this solution exists up to exponentially small corrections, which we will return to in Section VII). For small amplitude stars one can be sure that this expansion is dominated by the first harmonic with ω close to, but slightly less than, the mass m φ . For very large amplitude stars, especially those that are at masses approaching the maximum mass, higher order terms can become more important. By inserting this into the above equations of motion, along with a similar harmonic expansion for metric components A(r, t) and B(r, t), one obtains an infinite tower of coupled non-linear ordinary differential equations for the prefactors Φ n (r), A n (r), B n (r). In principle, one can numerically solve this complicated system. Interesting work on this was done in Refs. [44,45] for the φ 4 potential, where the tower of harmonics was constructed and numerical investigations were done. However, we leave a full treatment of all this for future work. A useful approximation is to assume that the system is dominated by only one harmonic of frequency ω. We write (the factor of √ 2 is for convenience). This is accurate at small amplitudes, as it is dominated by the fundamental in that case anyhow, and provides a rough estimate of the behavior at high amplitudes too. Of course by inserting this into the equations of motion it will not satisfy them precisely as it will generate higher harmonics. So the idea of the single harmonic approximation is to time average the equations of motion over a period of oscillation T = 2π/ω. Self-consistently then, we need to approximate the metric as time independent, writing A(r, t) ≈Ā(r) and B(r, t) ≈B(r). On the right hand side of the Einstein equations we need the time average of the energy momentum tensor, which is In addition, a useful way to handle the Klein-Gordon equation in eq. (12) is to multiply throughout by √ 2 cos(ωt) and then time average. This gives The time average of the potential V is simple if the potential has only even powers of φ. If there are odd powers of φ, like ∼ λ 3 m φ φ 3 , then the averaging is more subtle. Naively one just time averages those terms to zero. But in fact the presence of those terms means that the field undergoes asymmetric oscillations; it spends more time on one side of the potential than the other. A proper treatment of this would include a time independent term in φ, providing a non-zero mean φ = 0 (e.g., see the asymmetric terms in the study of oscillons in Ref. [46]).
However, we can instead pass to the low momentum effective theory; this is valid since we will only be studying boson star solutions whose size is much larger than the inverse mass of the scalar m φ . Here one integrates out intermediate φ exchange processes from ∼ λ 3 m φ φ 3 . This replaces interactions by only contact ∼ λφ 4 interactions, i.e., at large distances the interaction between scalars acts a kind of delta-function interaction. In the low momentum effective theory this leads simply to a shift in the effective quartic coupling, which one can readily show is (e.g., see Ref. [49]) In the effective theory one needs λ 4eff > 0 in order to have the repulsive interaction needed to support the stars of interest in this work (stars supported by quantum pressure are discussed briefly in Section VIII). The corresponding time averages are then (we truncate the potential to quartic order here for simplicity) Bound state solutions are found by numerically searching for configurations that obey the boundary conditions A → 1,B → 1, Φ → 0 as r → ∞ and Φ → 0 as r → 0. The ground state solution is identified by having no nodes.
One can readily integrate up the G tt = 8πG T tt equation to solve forĀ in terms of the enclosed mass M enc (r) asĀ The enclosed mass is just defined as the (weighted) integral of the energy density The corresponding mass of the star, or total energy, is Note thatĀ → 1 as r → 0. However, the value ofB or Φ as r → 0 is not specified uniquely. Let us call them B 0 ≡B(0) and Φ 0 ≡ Φ(0), respectively. Their values are related to the value of ω. It is useful to trade inB for another dimensionless function as By scanning different values of β 0 ≡ β(0), or equivalently by scanning different values for Φ 0 , and numerically solving the equations, one finds a range of boson star solutions. An example is given in Fig. 1

D. Large Coupling Regime
As shown in Ref. [34], the structure of the solutions is controlled by the following dimensionless parameter (Note: our g is what Ref. [34] calls Λ, but we have used Λ as strong coupling scale, and our λ 4eff /4 is effectively playing the role of λ in Ref. [34].) If g is small, then the self-interactions are negligible and the solution is just provided by gravity balanced by quantum pressure of the bosons (which is classical pressure within the field formalism); we shall study this in Section VIII. However, if g is large, then the behavior changes: there exist massive solutions in which gravity is balanced by the repulsive selfinteraction, which will be the main focus of this work (in addition there always exists very light solutions in which one can once again ignore self-interactions, but they are less interesting to us here). For glueballs, one anticipates m φ ≪ M Pl and λ 4eff a parameter that is not especially small, with λ 4eff ∼ (4π/N ) 2 . If we take N = O(2) and m φ ∼ 0.1 GeV then g ∼ 10 40 , i.e., it is extremely large.
In this large g regime, one can further simplify the equations of motion (as pointed out in Ref. [34]). To make it manifest which terms are important and which terms can be ignored at large g, it is useful to pass to the dimensionless variableŝ Exact Solution for g=100 Large g Approximation Boson star profile Φ(r) for the core value Φ0 ≈ 0.028 MPl and parameter g = 100. The orange is the exact numerical solution, while the purple is from scaling the large g approximation (this is similar in form to a plot in Ref. [34].) By re-writing the above Einstein-Klein-Gordon equations in terms of these variables and re-scaled metric component β and then neglecting all terms that are subdominant in the g → ∞ limit, one obtains the following simplified form of the equations One can readily integrate up the equation forĀ, as above, to obtain the enclosed mass in this limit as where In these dimensionless variables, one can anticipate that the maximum value of the integral here, when integrating over the whole star, is I max = O(1) (in addition to a family of lighter stars with I 1). In fact a precise calculation reveals that I max ≈ 0.2. Thus when multiplying by the prefactor in eq. (30) gives the maximum star mass Fig. 7), consistent with our earlier discussion in the introduction. The corresponding radius of the star is R * ∼ few × G M max , i.e., it is comparable to, though a little larger than, the Schwarzschild radius. Ratio of the mass of a boson star using the large g approximation scheme to the exact mass using the full equations (albeit always in the single harmonic approximation) with β0 = 1.585. Note that the ratio approaches 1 at large g.
The Klein-Gordon equation (29) can now be trivially solvedΦ(r) = β(r) − 1. This evidently only makes sense in the regime β > 1. For β < 1 the above approximations break down. This is because even though g is large the validity of this scaling rests upon the assumption thatΦ remains appreciable, say O(1). However, at large radius the exact solution has an exponentially small tail, as seen in Fig. 1 (orange curve), and so this assumption no longer holds. Instead this approximation is accurate in the bulk of the star. One can easily solve these approximate equations to compare. Fig. 1 (purple curve) shows a plot for g = 100, where the approximation is seen to be somewhat accurate in the bulk. However, as mentioned above, we are mainly interested in extremely large values of g, in which the approximation becomes extremely accurate (away from the exponentially suppressed tail).
Furthermore, we can compare the star's total mass M * in the exact versus approximate methods as a function of the dimensionless parameter g. We have fixed β 0 = 1.585 and numerically obtained the exact equations as well as the approximate result. The ratio of the masses is plotted as a function of g in Fig. 2. Note that the ratio M approx /M exact → 1 as g increases, as expected. For very large g it becomes prohibitively difficult to solve the exact equations. This is because at fixed Φ 0 , one has to search for the corresponding value of β 0 to obtain a zero node solution that asymptotes to flat space at large radii. This requires continually increasing precision. So instead by passing to the large g equations we can readily make progress in this regime.

III. QUANTUM ANNIHILATION IN CORES
The above analysis is all rigorous in the case of a complex scalar field ψ with a global U (1) symmetry. In that case the above solution is closely related to an exactly harmonic solution by writing ψ(r, t) = Φ(r) e −iωt . This corresponds to a boson star made out of a collection of particles (or anti-particles). The particles are then stable against annihilation into one another due to the symmetry. One may not expect global symmetries to be exact, but we put aside the discussion of this issue here. For the present purposes of a glueball or single elementary scalar φ, the boson star is made of of φ particles which are their own anti-particle. Their stability against annihilations is therefore not protected by any symmetry. Since the above very dense stars exist due to the selfinteractions λ n φ n , one should be concerned that these same interactions may lead to rapid annihilation in the star's core. Previous work on this topic has often ignored these details. A notable exception is the work of Ref. [43], which studied 3φ → 1φ annihilations (with a focus on axion stars); we shall return to a discussion of that work in Section VII.

A. Perturbative Rates
To discuss this issue, suppose we had the full classical field boson star solution, including its full time dependence. As we will return to in Section VII this can never be exact due to outgoing classical radiation, but that turns out be typically very small. On the other hand, at any order in the harmonic expansion, one can find a periodic solution that we can refer to as the classical star solution φ * (r, t); that will suffice in this section. The approximate classicality can be justified by various considerations, including averaging and decoherence at high occupancy [47,48]. In order to discuss the inevitable quantum fluctuations on top of this, we can work in the Heisenberg picture and write the field operator aŝ By treating the quantum fluctuations as small δφ φ * , we can write down the Heisenberg equations of motion and linearize them. At small couplings, one can then further obtain solutions for δφ working perturbatively. One assumes that the operator begins in the Minkowski vacuum and then evolves the system from there. Since the equations are linear this is doable in principle. If the background φ * were homogeneous, the perturbations would be easily diagonalized in k-space, leading to standard Floquet theory. However, since our background of interest φ * (along with the metric components A and B) depend on space, all the k-modes are coupled to one another. This makes the analysis non-trivial, even though the system is linear, but can be formulated as a generalized type of Floquet theory. This was all studied by some of us systematically in Refs. [49,50]; other work appears in Refs. [51][52][53][54][55]. The important result was that when one is in a small coupling regime, the final resonance matches the result from standard perturbation theory of particles annihilating in vacuum, i.e., the Bose enhancement is shut off for sufficiently small coupling. As shown in that work, the requirement for this is that the maximum parametric resonance Floquet rate µ from a homogeneous oscillating condensate is smaller than the inverse size of the star 1/R * . The reason for this is that in this regime the resonantly produced scalar particles escape the condensate before Bose-Einstein statistics can be effective. The regime of parametric resonance will be studied in next Section V.
In any case, the perturbative decay rates normally act as a lower bound on the true decay rate. The only reasonable ways they could be shut off is (i) if the decay products were fermions, then there is the issue of Pauli blocking [56]. However, this is irrelevant here since our decay products are the bosons φ themselves. Or (ii) if the star carried sufficient binding energy to make the process kinematically impossible. This is not an issue for dilute stars which carry very small binding energy, while highly compact stars will be addressed in Section VI.
These perturbative rates can be readily calculated. The boson star (at least away from very high compactness) can be viewed as a condensate of non-relativistic particles. Their large de Broglie wavelengths mean that they overlap with one another and can therefore annihilate through the contact interactions of the potential. This leads to the star emitting energy in the form of pairs of scalar particles. As shown in Ref. [49] the normalized perturbative rates for energy output Γ d = |dE * /dt|/E * (where E * = M * is the energy of the star) for 3φ → 2φ and 4φ → 2φ annihilation processes are given by where k 32 ≈ √ 5 m φ /2 and k 42 ≈ √ 3 m φ are the momenta of the 2 outgoing particles in each process, respectively. Here n * (x) is the number density of particles within the boson star, which can be approximated as n * (x) = ρ * (x)/m φ , where ρ * (x) is the local energy/mass density. Here the coefficients are given by where for completeness we have included the sextic term λ 6 φ 6 /(6! m 2 φ ) in eq. (34) (although we did not include it earlier when discussing the profile of the star). We note that for special choices of couplings, such as λ 2 4 = λ 6 , these rates can be shut off; but this is not the generic situation. Interestingly, this condition λ 2 4 = λ 6 occurs when one Taylor expands a cosine potential, which may describe an axion; however since it has attractive interactions, it is not of direct relevance here.

B. Application to Boson Stars
To evaluate these rates we need the star's density n * (x). It is difficult to provide exact analytical results as one needs to solve nonlinear differential equations to obtain the star's profile Φ(r); important work includes Refs. [57,58]. To proceed it is useful to use an approximate form for the shape of a star, which captures the following 2 basic ideas: (i) it is flat near its core, i.e., Φ (0) = 0, and (ii) it falls off exponentially at large distances. For example, as used in Ref. [59], a useful representation that has these properties and is somewhat accurate (at least for stars that are not too compact) is a sech ansatz where R is length scale that should be adjusted to minimize the energy to obtain the most accurate solution. The corresponding energy/mass density is ρ * (x) = m 2 φ |Φ| 2 , and the normalization ensures that the mass of the star is M * . This analysis is only accurate in the nonrelativistic regime of low compactness stars. For compact stars, they involve orbital boson speeds that are O(1/2) that of light. In that regime these estimates are still correct to within a factor of ∼ 2, which suffices for our main results.
Star solutions involve a relationship between the total mass/energy of the star M * and its physical radius R * . A precise definition of the physical radius of a star is the region that contains a fixed percentage of the mass of the star. For definiteness we will take this to be the radius that encloses 90% of the mass. This turns out to be related to the scale that appears in the argument of the above sech profile by R * ≈ d R, with d ≈ 2.8 for the sech profile. However, an O(1) change in this definition will not be important for our main results.
The radius-mass relationship for the condensate arises from minimizing the energy from kinetic energy (pressure), repulsive self-interaction, and gravitation. It can be shown this leads to [57][58][59] where a, b, c = O(1) numbers that depend on the choice of ansatz. For example, for the above sech ansatz they are: a = (12 + π 2 )/(6π 2 ), b = 6(12ζ(3) − π 2 )/π 4 , c = (π 2 −6)/(8π 5 ). Here we are using dimensionless variables Note that at largeM * 1 the radius approaches a constant R → 3c/b, giving a star size of This is as expected from the scalings of the full Einstein-Klein-Gordon equations in the previous section.
Let us focus on this asymptotic regime of the more massive stars, since they are of most interest astrophysically. We can readily carry out the above integrals to obtain the decay rates as

IV. ASTROPHYSICAL BOUNDS
If the above annihilation rates (38,39) are much bigger than the current Hubble rate H 0 of the universe, then the stars are unlikely to be cosmologically relevant. As an example, consider the case with λ = O(1) and m φ ∼ 0.1 GeV, giving rise to M max of the order of 10s of solar masses. Then, in order for the decay rates eqs. (38,39) to be smaller than today's Hubble rate H 0 ≈ 10 −33 eV, we obtain the following bound on the mass of the star: M * 10 −18 M sun if 3φ → 2φ processes are active and M * 10 −11 M sun if 4φ → 2φ processes are active. Such stars have very low compactness and are nowhere close to being relevant to the LIGO/LISA bands. This already undermines the claims of Refs. [27,[38][39][40], as well as other work on real scalars, including Ref. [60], since the heavy stars of these analyses would be too short lived.
A much more general exploration of the constrained parameter space is provided in the Figs. 3,4,5,6, where we have included detailed information in their captions. The decay rates Γ d are plotted versus star mass M * in Figs. 3,4, with the coupling λ = (2π) 2 fixed in Fig. 3, and the scattering cross-section fixed in Fig. 4 (see next subsection for explanation). Then in Figs. 5,6 we fix the decay rate to be today's Hubble rate Γ d = H 0 , with contours of fixed mass M * indicated; these provided upper bounds on the allowed mass for the star to live to present time. We also indicated the compactness (see upcoming subsection for explanation).
In the remainder of this section and the next section we would like to provide more details on the ingredients that have gone into these figures.

A. Scattering in Galaxies
Apart from the boson stars, one expects there to be a large collection of diffuse φ particles acting as a form of dark matter. Due to the above self-interactions they will undergo scattering in the galaxy. The 2φ → 2φ scattering cross section for non-relativistic particles is readily obtained as If the scalar particles φ make up a significant fraction, or all, of the dark matter then there are bounds on this scattering. Some of the best bounds come from observations of collisions of galaxies, such as the bullet cluster, which are essentially consistent with non-interacting dark matter. This imposes an observational upper bound on the scattering cross section of [28,29] σ 2→2 m φ ζ u cm 2 /g where ζ u is argued to be O(1), depending on the analysis. For concreteness, we will take the bound with ζ u = 1 in this work. In Figs. 3,5,6, we have imposed this bound, which rules out the green region. We add that if the scalar particles φ only make up a tiny fraction of the dark matter (yet there is still enough to provide some boson stars), then these bounds are significantly weakened; in this case the green region can largely be ignored. On the other hand, it has been argued that the cores of galaxies are explained precisely by the presence of such scattering [61]. To do so requires a lower bound on the cross section of σ 2→2 /m φ ζ l cm 2 /g, where ζ l has also argued to be O(1). Altogether this provides some motivation to consider the case of σ 2→2 /m φ ∼ cm 2 /g, although it is still controversial. (See Ref. [62] for a critical examination of using ultralight scalars to solve this problem.) We have fixed the scattering to be σ 2→2 /m φ ∼ cm 2 /g in Fig. 4. In that plot we are showing decay rate versus star mass M * for different choices of the coupling λ. In the upper left region one would have λ 1, which is forbidden by unitarity. For concreteness we take λ = (2π) 2 as the upper value (although one could argue for even smaller values to be safer from unitarity considerations).

B. Compactness
Of particular interest for the production of gravitational waves is the compactness C of a star. We shall define this as the ratio of the corresponding Schwarzschild radius R S = 2GM * and the physical radius of the star If the compactness is O(1), then it is undergoing strong gravity in its vicinity. When M * → M max , as defined in Section II D, then we are indeed in this regime. The compactness there is roughly C ∼ 1/2 and mergers can produce gravitational waves with significant amplitudes. Of course the compactness parameter cannot be larger than 1 or the system would have collapsed to a black hole. In Figs. 3,4,5 we have indicated this in the orange region (In Fig. 6 this occurs at much higher values of C than those displayed). On the other hand, the compactness can be small C 1, describing the dilute boson star. In this regime the gravitational wave signal from mergers is expected to be suppressed.
We can consider a family of solutions that exist at some compactness C. For example, one may imagine that the stars will continue to accrete until they have achieved their maximum compactness C ∼ 1/2. For any particle mass m φ we can consider this possibility. So we use eqs. (37,42), eliminate m φ in favor of C, and insert into the decay rate formulas to obtain where β 3 , β 4 are O(1) prefactors. The sech ansatz estimates their values to be β 3 ≈ 0.005, β 4 ≈ 0.0001. Contours of fixed compactness are provided in Fig. 3 for fixed coupling and in Fig. 5 for fixed decay rate.

C. Implications for Gravitational Waves
Importantly, by fixing the decay rate to be Hubble, we plot contours of maximum allowed star mass M * in the {C, λ}-plane in Fig. 6. This shows that for reasonable couplings, and to obey the bullet cluster bound, the mass and compactness need to be small for longevity. This has significant implications for gravitational waves from possible mergers of these stars. It implies that one is well outside of both the LIGO and LISA bands. Gravitational wave detection requires O(10 − 10 6 ) solar mass objects with high compactness. However these results indicate that the decay rates are too fast to achieve this. For completeness, we have indicated both the LIGO and LISA bands in Figs. 3,4,5. In Figs. 3,4 the signals are only appreciable in the upper right hand region, which involves rapid decay. While in Fig. 5, where the decay is fixed to Hubble, the signal is only appreciable in the lower left hand region, involving extremely tiny self couplings, and may be unlikely to persist due to the phenomenon of parametric resonance, which we turn to now.

V. PARAMETRIC RESONANCE
In addition to the quantum mechanical perturbative decays discussed above, one can also enter a regime in which the coherently oscillating boson star condensate drives parametric resonance of its own field fluctuations. Since the stars of interest are wide compared to the inverse mass of the particle (R * ∼ √ g/m φ 1/m φ ) they are susceptible to parametric resonance. This would represent a type of instability against linear perturbations; differing views on this have been expressed in the literature [63,64]. As mentioned earlier, the criteria for this to occur is µ R * > 1 (45) where µ is the maximum exponential growth rate ("Floquet exponent") within the homogeneous background approximation. The intuition behind this is that when this inequality is not satisfied, the produced particles escape the condensate before Bose-Einstein statistics are effective. In any case, it was earlier established in Ref. [49] (also see earlier work in the context of axion-photons in Refs. [51,52,55]), so we shall not repeat the derivation here. Hence we do not need to consider the full complications of expanding around the inhomogeneous star background, we can focus on the corresponding homogenous configuration with amplitude matching the star's core amplitude, obtain µ, and check on this inequality. Let us perturb around the background as Here one should, in principle, also allow for fluctuations in the metric. However, the metric fluctuations are primarily only important to describe long wavelength perturbations around the homogeneous background. This leads very importantly to the collapse of homogeneous structure and is the source of structure formation, and ultimately to the formation of boson stars, etc. However, what we are interested in is to imagine a star has formed and we are only interested in these annihilation processes inside its core. These are particle number changing processes and are mediated by self-interactions of the potential, and are not mediated by gravity (one could consider resonance into gravitons but this is normally highly suppressed). Hence we can focus on flat space fluctuations for the purpose of this discussion. The perturbed equation of motion is where V I is the interaction potential. The advantage of studying this homogenous condensate is that it can be readily diagonalized by passing to Fourier space δφ → δφ k and −∇ 2 δφ → k 2 δφ k . This makes eq. (47) a form of Hill's equation since it is a linear differential equation with a periodically changing prefactor V I (φ 0 (t)).

A. Floquet Exponents
Let us begin with the case when both odd and even powers in the potential are included. The most obvious version is the cubic term λ 3 m φ φ 3 /3! and quartic λ 4 φ 4 /4!. However to obtain the relevant resonance is slightly complicated (as can be appreciated by drawing all the corresponding Feynman diagrams). So for the sake of simplicity, let us focus on an interaction provided by V I = λ 5 φ 5 /(5! m), which provides 3φ → 2φ with no intermediate propagators. Although, one should also anticipate the presence of cubic terms, but the final result will have a similar scaling (assuming λ n ∼ λ (n−2)/2 ).
To first approximation, the oscillating condensate is mainly driven by the mass term. So the oscillations of the background are where ω 0 ≈ m φ and φ a is the amplitude of oscillations. By inserting this back into eq. (47) one obtains the interaction term as a pair of harmonics due to the factor (3 cos(ω 0 t) + cos(3ω 0 t)) (49) One can then expand δφ in harmonics too. One can have resonance from long wavelength perturbations, which can be driven by the leading harmonic term. However, as mentioned above this is not important for us. We know that this only leads to the destabilization of the homogeneous condensate towards a boson star etc. Instead we are interested in the possible resonance at the higher harmonic 3ω 0 ≈ 3m φ . This is potentially resonant for δφ k ∝ e i3ω0/2 , since when we insert this into the equation of motion, the driving term will have a frequency that matches the input frequency. In turn this matches the natural frequency for ω k ≡ √ k 2 + m 2 = 3ω 0 /2 ≈ 3m φ /2, which means k ≈ k 32 ≈ √ 5 m φ /2 the outgoing wave number mentioned earlier in our perturbative analysis in eq. (33).
Hence in the vicinity of this resonance of interest, we can write the equation of motion as where "n.r" refers to the non-resonant term from the cos(ω 0 t) in eq. (49). Ignoring the non-resonant piece, this is a form of the Mathieu equation which is known to possess exponential growth in some band of wavenumbers. Here we can identify For small amplitudes, the Floquet exponent is known to be [65] As we scan over different k-values this is clearly maximal when A k = 1, i.e., for ω k = 3ω 0 /2 as expected. This gives the maximum Floquet exponent of (using ω 0 ≈ m φ ) where we have indicated in the final step that when we include the cubic and quartic couplings, we can generalize the result to |λ 5 | → |λ 5 + 5λ 3 (λ 2 3 + 3λ 4 )/12| = λ 3/2 32 , since we know this arises from the scattering amplitudes.
If there are no odd powers of φ in the potential, we can still have parametric resonance from the quartic term V I = λ 4 φ 4 /4!. In this case the analysis is rather more complicated. This can be seen by the fact that there are now several Feynman diagrams associated with the 4φ → 2φ process. The full details of this analysis was carried out by some of us in Ref. [65]. A simpler calculation arises from a sextic term V I = λ 6 φ 6 /(6! m 2 φ ) because it occurs through a process without any intermediate processes.
For this we can again follow again the above analysis, this time expanding V I (φ 0 (t)) ∝ φ 4 0 (t) in harmonics and focussing on the cos(4ω 0 t) term. The result is where we have indicated that we can generalize this to λ 6 → |λ 6 − λ 2 4 | = λ 2 42 , since we know this is the relevant contribution from both processes.

B. Application to Boson Stars
We now would like to apply these results to the boson star. Firstly we need to relate φ a to the properties of the star. The idea of the criteria µ R * > 1 for resonance is that one replaces the the amplitude of the homogenous configuration φ a by the corresponding amplitude at the core of the star, i.e., where f is yet another O(1) prefactor; in the sech ansatz it is given by f = 3 d 3 /π 3 ≈ 1.46. Hence the product of interest is found to be where γ 3 = f 3 b 7/4 /(108 √ 2 3 3/4 c 7/4 d 7/2 ) and γ 4 = f 4 b 5/2 /(3456 √ 3 d 5 c 5/2 ). In the sech ansatz they are given by γ 3 ≈ 2.1 and γ 4 ≈ 1.8.
The region in which the inequality for parametric resonance is obeyed µ R * > 1 is indicated by the light blue region in Figs. 3,4,5,6. We see that it can constrain the allowed parameter space considerably. However, its scaling is different to the perturbative decays, so it is often complementary.
If we now eliminate the particle mass m φ to rewrite µ R * in terms of the compactness C, as we did in Section IV B, we obtain where δ 3 , δ 4 are O(1) prefactors. The sech ansatz estimates their values to be δ 3 ≈ 0.4, δ 4 ≈ 0.1. This result suggests something important: By imposing the bound on mass M * M Pl / √ λ 4eff , so that we are in the selfinteraction supported regime (as opposed to the quantum pressure supported regime) and using λ n ∼ λ (n−2)/2 , we have µ 3→2 R * C and µ 4→2 R * C 3/2 . This means that if one considers stars of large compactness, say C ∼ 1/2, then the inequality eq. (45) is satisfied and parametric resonance is expected to occur. This suggests very compact stars that could undergo mergers and be relevant to LIGO/LISA are likely to be quite unstable to parametric resonance. However, since we have not solved for the time dependence fully we can not be certain of this conclusion. This leads us to search for additional clues, as we do in the next section.

VI. BINDING ENERGY
The results of the previous sections indicate that any massive compact stars will decay rapidly. However, one might consider the possibility that the boson stars carry so much binding energy that these processes are shut off. The above analyses were done reliably in the weak gravitational field regime, where the binding energy is small and is unlikely to be able to shut off these processes. However, in the strong gravity regime this becomes at least conceivable.
To examine this fully for a real scalar, we would really need to solve the full set of time dependents equations by expanding in a tower of harmonics, as we described in Section II C. However, as we also discussed there, this is a rather difficult task and will be left for future work.
For now we will simply take our clues from the much simpler case of a complex scalar field theory. To define this, we can return to our starting action eqs. (4), replace the real scalar φ by a complex scalar ψ, and endow the theory with a global U (1) symmetry. To be concrete, the kinetic term is now ∆L = −|∂ψ| 2 /2 and we choose the potential now as V = m 2 φ |ψ| 2 /2 + λ 4 |ψ| 4 /16. The U (1) symmetry ensures there is a conserved current [66] with a corresponding conserved particle number One can now rigorously define the boson star as a state that minimizes the energy subject to the constraint of (or ψ(r, t) = Φ(r) e +iωt for a star of anti-particles). The conserved particle number is Note that for the complex scalar, the metric of a single boson star is exactly static, so we could writeĀ = A,B = B here. By solving the large g Einstein-Klein-Gordon equations, we can determine the star's energy/mass and number for a family of solutions; this is given in Fig. 7. The left hand solid curves are the usual solutions that we have focussed on in this work; they are stable for a complex scalar. The right hand dashed curves may have instabilities of a variety we have not discussed in this work; we leave their analysis for future work. Now we would like to infer any ramifications for the case of a real scalar. We can be concrete about this in the following way: Suppose we have the above complex scalar with a U (1) symmetry, and now we introduce small U (1) breaking terms, which can in principle mediate particle number changing processes. In order for gravity to prevent an annihilation process of the form n φ → 2φ, the change in energy with respect to particle number would need to be less than 2m φ /n. So to kinematically forbid 3 → 2 annihilations one would need dE * /dN * < 2m φ /3, while to kinematically forbid 4φ → 2φ annihilations one would need dE * /dN * < m φ /2. In Fig. 8 we show our results for dE * /dN * as a function of particle number. We see that for all stars up to the maximum value allowed (solid branch), we find (and dE * /dN * > 0.79 m φ including the dashed branch). So, while there can be an appreciable amount of binding energy for the most massive stars, it does not appear to be quite enough to prevent annihilations. So although we have not performed the full analysis for the real scalar, this provides some circumstantial evidence that binding energy, even for the most compact stars, will not be enough to prevent the decay processes computed in the previous sections.

VII. CLASSICAL "ANNIHILATION" IN CORES
So far we have considered quantum annihilation in the core of a star. What we mean by this is the expansion given in eq. (32), i.e., we take the star solution φ * (r, t), which we take to be exactly periodic and localized, and we add the inevitable fluctuations required by quantum mechanics by the operator δφ(x, t). This is very important in the perturbative regime, since it is this bath of quantum fluctuations that get driven through forced resonance to provide a stream of outgoing particles, with corresponding annihilation rates Γ d (3φ → 2φ) or Γ d (4φ → 2φ). As a matter of principle it is also important for the parametric resonance, even though the growth (Floquet) rates can be computed classically. That is because one also needs to perturb the solution in order to see the growth. Such perturbations can very easily happen from imperfect initial conditions around the star solutions, and so could be seen classically, but if one starts exactly on the star solution, then it is quantum fluctuations that are important.
However, in addition to this is the fact that the boson star solution φ * (r, t) is in general not expected to be an exact solution of the classical equations of motion. This is because it is generally very difficult to have any exact periodic solutions of non-linear partial differential equations. One of the only known counter examples is the sine-Gordon breather in 1 + 1-dimensions, but that is not our focus here. So if one expands the solution in harmonics at small amplitudes, one anticipates that the expansion is only an asymptotic series, missing exponentially small terms. These corrections are sometimes referred to as part of a "hyperasymptotic series" [67]. Such series have been addressed in the literature for quite sometime [68], so we will only report on some of its basic features here.
In order to address this, one can write the full classical boson star solution as where the sum is the exactly periodic piece, summed up to P terms in the harmonic expansion, and χ(r, t) is the inevitable correction that survives at that order. If we then insert this into the full equations of motion and expand to linear order in χ (since χ is expected to be very small), the structure of the resulting equation is where J(r, t) arises from the star solution which is never exact at any order P in the expansion. Due to the presence of nonlinear terms, the driving term J is will take the form J(r, t) = j(r) cos(ω 2 t) + . . .
where ω 2 is the frequency of the leading harmonic that is generated by the nonlinear equations of motion from the fundamental ω 1 . For a theory with odd powers of φ in the potential, ω 2 = 2ω 1 , while for a theory with only even powers of φ in the potential, ω 2 = 3ω 1 . Here j(r) is a function of radius that is determined by the star solutions Φ n (r). Eqs. (66,67) imply that the boson star acts as a coherent oscillating source that is generating its own classical scalar radiation χ. The relevant solution of eq. (66) is the particular solution, which is readily obtained in Fourier space as χ = d 4 k/(2π) 4 J(k, ω)e i(k·x−ωt) /(−ω 2 + k 2 ). In the far distance regime, well outside the star, the tail of this radiation takes the form where k 2 is the wavenumber of the outgoing on-shell radiation with ω 2 = k 2 2 + m 2 φ and j(k 2 ) is the Fourier transform of j(r) evaluated at k 2 .
As demonstrated in Ref. [49] the value j(k 2 ) is on the order of the star's solution itself j(k 2 ) ∼ m 2 φ φ * (k 2 ) (evaluated at, say, t = 0); this is reasonable since it is the star solution that provides the source for its own radiation. Then using eq. (68) the corresponding power output is therefore where in the last step we are using the single harmonic approximation, which was studied in detail in Section II C, and we are using ω ∼ m φ here. For the sake of concrete analytical results, let us use the sech ansatz of eq. (35), which is most trustworthy for dilute stars. We can readily obtain its Fourier transform as We need to evaluate this at the on-shell wavenumber k = k 2 = ω 2 2 − m 2 φ . For stars of low compactness, we have ω 1 ≈ m φ . So then can write k 2 = κ 2 m φ , where κ 2 = √ 3 when there are odd powers of φ in the potential, and κ 2 = √ 8 when there are only even powers of φ in the potential. The stars of interest are much wider than the inverse mass of the field, i.e., R m φ 1. So we are deep into the tails of both the tanh and sech functions above. In this regime, we can approximate tanh(πk 2 R/2) ≈ 1 and sech(πk 2 R/2) ≈ 2 exp(−πκ 2 m φ R/2). The corresponding decay rate due to this classical radiation output is Γ class = |dE * /dt|/E * (E * = M * is energy of the star). This gives (dropping O(1) factors) where c is an O(1) number; which has the value c = π/d ≈ 1.1 within the sech ansatz. However, we expect there to be an O(1) correction to c from the real solution, so its specific value here is not the focus. This leads to the expectation of exponential suppression. However, we note that the above Fourier transform is really only valid in a small amplitude expansion, ensuring that the star is wide and dilute and enters a scaling regime controlled by one scale R (we shall return to study this case in detail in Section VIII). In the more massive branch of solutions of most interest in this paper, the star actually has new features; including the somewhat more vertical shape at the edge of the star seen in Fig. 1 when g is large. Hence this simple estimate of the Fourier transform is only trustworthy for dilute stars, and may be less applicable for dense stars with another scale entering the analysis. Nevertheless, we only wish to comment on a very basic qualitative feature of the solution: In particular, for stars of large radius m φ R * 1, one can still anticipate the Fourier transform is somewhat small at k = k 2 = κ 2 m φ ∼ m φ . Indeed we know that for boson stars supported by repulsive self-interactions, the radius is on the order the coupling introduced earlier R * m φ ∼ √ g ∼ √ λ 4eff M Pl /m φ , which is assumed to be large so far in our work.
Hence, the classical radiation for M * M max boson stars is expected to be small. There may be a large change in this conclusion for somewhat compact stars, wherein more bosons are circulating in a relativistic fashion. In this case the Fourier transform could have appreciable support at the relevant wavenumber k 2 = κ 2 m φ ∼ m φ , and the classical radiation could be significant. But from our analyses in earlier sections, quite compact stars are already anticipated to have significant perturbative, and possible parametric, decays anyhow.

A. Comparison to Literature
The above analysis involves the star itself driving its own classical radiation. In a sense this can be viewed as a kind of 2 → 1 process, in the case in which there are odd powers of φ in the potential, or a a kind of 3 → 1 process, in the case in which there are only even powers of φ in the potential. While this does not immediately seem to be compatible with momentum conservation, if we recall that the outgoing waves are spherically symmetric, there is no real problem.
In Ref. [43], a similar type of 3 → 1 process was considered. There the focus was on axions in which there are expected to be only even powers of φ in the potential. In that work they considered the limit in which gravity is decoupled. In this regime, one is not in fact even studying boson stars. Instead in this regime these are oscillons ("axitons" in Ref. [69]); scalar field solutions in which the attractive self-interactions (λ 4 < 0) are balanced by the field's gradient pressure ("quantum pressure"). Even though this is a different regime to what we are mainly focussing on in this work, we can nevertheless compare the basic form of the radiation formulae.
In this case there is a well defined small amplitude expansion (although it suffers from a collapse instability in 3+1-dimensions), and one can compute the radiation output from oscillons (for example, see Refs. [46,49]) where the coefficient of the exponential can be reliably found to bec ≈ 1.2. Here ε 1 is a small expansion parameter; it is related to the shift in frequency of the fundamental by ω 2 1 = m 2 φ (1−ε 2 ). In this small amplitude regime, it also sets the size of the oscillon to be m φ R * ∼ 1/ε. Hence the scaling here is similar to the case above in eq. (71), albeit a difference is that R * is fixed by the gravitational mass-radius relation eq. (36) for the stars of interest, while it is not fixed to this value for the oscillon.
On the other hand, in Ref. [43] the corresponding 3 → 1 process was (i) claimed to be quantum mechanical and (ii) had a different scaling. In that work they claimed Γ ∼ m φ λ ε 2 exp(−c κ 2 /ε) (Ref. [43] result) (73) To compare notation to Ref. [43] replace ε → ∆ and λ → m 2 φ /f 2 a , where f a is the PQ scale of the axion. We note that while the exponent is the same as our results above, the prefactor of Ref. [43] is enhanced by a factor of 1/(λ ε). While this may not be too important if one focusses on glueballs with N ∼ 2 giving λ = O(1), it can be important for extremely small λ. In particular, for the axion, one does in fact anticipate such extremely tiny values. For example, for the QCD axion, one may conceivably have m φ ∼ 10 −5 eV and f a ∼ 10 12 GeV, giving 1/λ ∼ 10 52 . Combined with ε 1, this factor of 1/(λ ε) can be a huge enhancement.
However, we would like to point out that in Ref. [43] the result that they reported in eq. (73) is in fact simply the rate at which 3 particles are converted into 1 particle. This is not the decay rate of the condensate. To obtain the decay rate of the entire condensate, one must divide by the number of particles in the condensate (or say, half the number) in order to obtain the characteristic time for an appreciable fraction (say half) of the condensate to radiate away. The number of particles in the oscillon condensate, in this small amplitude expansion, can be readily shown to be N * ∼ 1/(λ ε). So by forming Γ/N * , we recover the correct scaling of eq. (72).
We note that if one were to persist with eq. (73) as a decay rate, one would find something very peculiar: if we were to reinstate factors of , starting with the classical field theory, one can readily show that it is proportional to 1/ . This means its classical field theory limit is badly behaved. Instead, by noting that this is not the physical rate of decay of the macroscopic condensate, but instead we need to divide by a factor on the order ∼ N * ∼ 1( λ ε), we then obtain the result of eq. (72) which is independent of and is indeed classical.

VIII. QUANTUM PRESSURE SUPPORTED STARS
For completeness, let us also briefly discuss more familiar boson stars: those that are not supported by repulsive λφ 4 interactions, but instead are supported by the field's gradient pressure (which is often referred to as "quantum pressure", as it is quantum mechanical from the particle point of view, since it originates from the de Broglie wavelength of the bosons). In this regime the couplings can in fact be taken to zero λ → 0 and the star's still persist. It can be readily shown that in the λ → 0 limit the maximum mass of such stars is on the order (e.g., see Ref. [70]) Which is a very different scaling to the maximum mass in the self-interaction supported regime of M max ∼ √ and long lived. We further investigated the possibility of parametric resonance of fluctuations in Section V, with the relevant dimensionless parameter being µ R * to indicate whether the Bose-Einstein statistics are effective or not, given in eqs. (56,57) or eqs. (58,59). We found that the resonance becomes the dominant mechanism for stability at smaller couplings and can considerably constrain the parameter space further. An important caveat is that all of our analysis was done with a simple starting point for the star, in which we describe it within the time averaged single harmonic approximation. To then ascertain decay rates, we then perturb around this to obtain quantum decay rates by the Heisenberg equation of motion due to forced resonance or exponential growth rates due to parametric resonance. This strategy should be valid for low compactness stars, which are essentially non-relativistic and possess a reasonable single harmonic approximation. However, this could be less accurate at high compactness, where many harmonics are expected to play a role, and this simple analysis could conceivably miss out on some of the physics. We did study the binding energy of the star in Section VI within this simplified treatment, finding that it does not appear to be enough to prevent decays, but a full treatment would be preferable.
Important future work is therefore to included the full time dependence, which can be written as a collection of harmonics and then to solve a system of ODEs for the coefficient functions of radius. Alternatively, one could run full simulations, which has the advantage of not only capturing all harmonics, but identifying instabilities readily. This is at least true for the instabilities associated with parametric resonance, while quantum radiation can be much harder to see in a classical simulation, as it requires seeding the fluctuations with a bath of zero point fluctuations, which can be difficult to track reliably in a finite simulation.
We also examined the classical decay rate in Section VII, which is ordinarily much smaller because the star's are wide. This means that, except perhaps for somewhat compact stars, their Fourier transform is typically small at the relevant resonant wavenumber, so they do not efficiently emit their own classical radiation. Conversely, the quantum radiation is set by and is not suppressed in most regimes of interest. We compared the basic structure of the classical rate to existing claims in the literature [43], regarding 3 → 1 processes for dilute axions, in which they studied the limit in which gravity is decoupled. Here the condensate is in fact a type of oscillon (or "axiton") and is held together by attractive selfinteractions. The work of Ref. [43] quite rightly claimed the process was exponentially suppressed at small amplitudes, but we found that they greatly over estimated the prefactor by at least 1/λ; this can be a huge enhancement for other dark matter candidates, such as axions. Instead, we explained that it is necessary to divide their rate by the number of particles in the condensate to obtain the physical decay rate, rather than merely the rate for any particle to meet others. This leads to a final result which is independent of , and is therefore classical, as expected.
For completeness, we then considered the more regular branch of solutions of boson stars in Section VIII, namely those that are supported against gravity by "quantum pressure", rather than repulsive self-interactions. In this regime, we found that both the perturbative and parametric resonance rates are suppressed. This is not too surprising, since the appearance of this branch only emerges in the limit in which self-interactions are small, but it is those same terms that are trying to drive the number changing processes. In any case, these stars are therefore the most robust against decays. On the other hand, they have the parametrically much smaller maximum mass of M max ∼ M 2 Pl /m φ , which can only be astrophysically significant for extremely small m φ , requiring extremely small λ m 2 φ /M 2 Pl to avoid the self-interaction corrections.
Other future directions are to return to consider the complex scalar ψ with a global U (1) symmetry [34]. This immediately prevents the decay processes studied here, since it has a conserved particle number. However, there are suggestions from quantum gravity that all global symmetries should be explicitly broken. This could introduce terms like ∆L ∼ ψ 4+n /M n Pl + c.c, and leads to the expectation that there may be particle number changing processes allowed in the star. If n = 1 this is roughly equivalent to replacing λ 3/2 32 → m φ /M Pl in our above decay rates, or if n = 2 this is roughly equivalent to replacing λ 2 42 → m 2 φ /M 2 Pl . This can be such a large suppression that it significantly opens up much more parameter space where the stars are stable; however, we leave a full investigation for future work.
Furthermore, one can consider other types of novel compact stars that may be allowed by other interesting dynamics in dark sectors. Alternatively, one may focus on other kinds of observational consequences that may arise from dilute stars, including fast radio bursts [53] or pulsar timing [71].