Dilute and dense axion stars

Axion stars are hypothetical objects formed of axions, obtained as localized and coherently oscillating solutions to their classical equation of motion. Depending on the value of the field amplitude at the core $|\theta_0| \equiv |\theta(r=0)|$, the equilibrium of the system arises from the balance of the kinetic pressure and either self-gravity or axion self-interactions. Starting from a general relativistic framework, we obtain the set of equations describing the configuration of the axion star, which we solve as a function of $|\theta_0|$. For small $|\theta_0| \lesssim 1$, we reproduce results previously obtained in the literature, and we provide arguments for the stability of such configurations in terms of first principles. We compare qualitative analytical results with a numerical calculation. For large amplitudes $|\theta_0| \gtrsim 1$, the axion field probes the full non-harmonic QCD chiral potential and the axion star enters the {\it dense} branch. Our numerical solutions show that in this latter regime the axions are relativistic, and that one should not use a single frequency approximation, as previously applied in the literature. We employ a multi-harmonic expansion to solve the relativistic equation for the axion field in the star, and demonstrate that higher modes cannot be neglected in the dense regime. We interpret the solutions in the dense regime as pseudo-breathers, and show that the life-time of such configurations is much smaller than any cosmological time scale.

In this article, we study the stability of axion stars as a function of the amplitude of the axion field at the core of the star |θ 0 | ≡ |θ(r = 0)|. Our results apply to the full range of axion masses for which QCD axions can comprise all of the dark matter. We identify three distinct branches of axion stars, distinguished by the field amplitude at the core, which in turn determines the density * Electronic address: luca.visinelli@fysik.su.se † Electronic address: sbaum@fysik.su.se ‡ Electronic address: jredondo@unizar.es § Electronic address: ktfreese@umich.edu ¶ Electronic address: wilczek@mit.edu of the star. We should keep in mind, that the axion is a periodic field with amplitude effectively restricted to the domain 0 ≤ |θ 0 | ≤ π. For small field values |θ 0 | 10 −6 10 −5 eV/m with m the axion mass, the axion field only probes the harmonic part of the potential, and it can be treated as a free field. In this regime, self-gravity is balanced by the kinetic pressure arising from the uncertainty principle. We call this the dilute axion star branch. We reproduce the previous findings in the literature for the mass-radius relationship, R ∝ M −1 , where R and M are the radius and mass of the star, respectively. In this regime, the configuration is stable against perturbation: For a given mass M , stars are pulled back to the equilibrium radius if they expand because then the (attractive) self-gravity is stronger than the (repulsive) kinetic pressure; conversely, if they are perturbed to smaller radii, they expand because kinetic pressure becomes stronger than self-gravity.
For configurations with M ∼ 10 −11 M 10 −5 eV/m 2 , self-interactions cannot be neglected anymore, although the amplitude is still comparatively small, |θ 0 | ∼ 10 −6 10 −5 eV/m . For QCD axions, the lowest order self-interaction is an attractive quartic term. For amplitudes |θ 0 | 10 −6 10 −5 eV/m , the attractive quartic self-interaction is stronger than gravity, which is negligible in this regime. In this critical branch, we find solutions when the quartic self-interaction balance the kinetic pressure with mass-radius relation R ∝ M . Note, that this relation implies that axion stars become lighter with growing density, such that they always have masses M 10 −11 M 10 −5 eV/m 2 in this branch. However, the solutions are unstable against perturbations: for a given mass M , stars expand when perturbed to radii larger than the equilibrium value since the quartic self-interactions are weaker than the repulsive pressure. Eventually the configuration relaxes to the dilute, interaction-free regime described in the previous paragraph. Conversely, if configurations are perturbed to radii smaller than the equilibrium value, the quartic interaction is too strong to be balanced by the pressure and the star collapses to even higher densities. It has recently been pointed out, that new stable configurations, called dense axion stars, are obtained when the amplitude of the axion field in the core reaches |θ 0 | = O(1) [27]. For such amplitudes, the axion field scans the full non-perturbative axion potential, and selfinteractions must be taken into account to all orders. Using the assumption that the axion field in the star is coherently oscillating at a single frequency, as commonly used in the literature, we obtain the mass-radius relation M ∝ R 3 , in agreement with Ref. [27]. However, we find that the single-harmonic approximation, which holds in the branches described above, is not accurate for the dense branch. Using a multi-harmonic expansion, we find that higher harmonics are generated with amplitudes comparable to the fundamental mode's amplitude. Heuristically, the presence of higher harmonics corresponds to the generation of (relativistic) axions by coalescence processes na → a. We find that configurations on the dense branch decay via emission of relativistic axions, with lifetimes of order τ life ∼ 10 3 /m, which are much shorter than any cosmological timescale.
When |θ 0 | O(1), axions stars are short lived solutions of the relativistic equation, elsewhere known as oscillons [37][38][39][40][41][42][43][44][45]. In the literature, similar objects have also been called pseudo-breathers [46], axitons [26], or oscillatons when driven by gravity [23,47,48]. Since gravity is negligible in the dense branch, the axion field is described by the Klein-Gordon equation with the QCD chiral potential (the χ-Gordon equation). There is a large but scattered literature on finding solutions to related equations. For example, in one dimension, assuming a cosine potential leads to the Sine-Gordon equation, which admits localized breather solutions that are not harmonic [49], i.e. which feature an infinite collection of higher harmonics. In three dimensions oscillons closely resemble the breather solutions of the one dimensional Sine-Gordon equation, but they differ in that they radiate energy and thus decay in a finite lifetime, though slowly relative to the "natural" timescale set by the inverse mass of the particles.
Justifying and expanding upon this concise summary, the remainder of this paper is as follows. In Sec. II we set out the basic equations. In Sec. III we find numerically stable solutions and provide quantitative results for the dilute and the critical axion star branches. In Sec. IV we discuss the dense branch, and analyze the equilibrium and metastability of dense configurations in a relativistic framework. In Sec. V we use the mass-radius diagram to sketch a qualitative storyline for axion stars, and in Sec. VI we summarize and conclude.

A. Axion Lagrangian
The axion results from promoting the flavor-neutral CP violating angle of the standard model, θ, to a dynamical field [1,2] in the Peccei-Quinn mechanism [10,11]. The canonical normalization of the dynamical angle θ(x) requires a new energy scale f , the axion decay constant, to define the axion field a(x) = θ(x)f . In the following we will refer to both θ and a as the axion field. The dynamics of the axion field under the influence of gravity are described by the action where the metric g µν is determined by the Einstein equation for the energy momentum tensor of the axion field T µν (a). We adopt the axion potential [50,51], where Λ 4 ≈ (75.5 MeV) 4 is the topological susceptibility [51][52][53] and c z ≈ z/(1 + z) 2 ≈ 0.22 with the ratio of the up and down quark masses z = m u /m d ≈ 0.48. Note, that the minimum of the potential is at V (0) = 0 and the maximum at V (π) = Λ 4 1 − √ 1 − 4c z /c z . The axion mass m and the quartic coupling constant λ are defined through Assuming spherical symmetry and expanding the metric to linear order about flat space yields the line element where φ is the gravitational potential, which satisfies the Poisson equation with energy density ρ = T 00 (a), and dΩ is the differential solid angle. In the following, we rescale time and radius as t → mt and r → mr, respectively, so that the Lagrangian in Eq. (1) reads where a dot indicates a derivative with respect to the rescaled time, a prime indicates a derivative with respect to the rescaled radius, andṼ (θ) ≡ V (θ)/Λ 4 . Coupling the Poisson equation with the equation of motion obtained from the Lagrangian density L gives where β ≡ Gf 2 = (f /m Pl ) 2 with the Planck mass m Pl = 1.221 × 10 19 GeV. The energy densityρ ≡ ρ/Λ 4 is dimensionless, and reachesρ ∼ 1 when |θ| ∼ π and the axion potential saturates. In Eq. (9), we denote the contributions to the energy density from the kinetic, gradient, and potential components separately. Note, that the gradient energy is due to the momentum of the axion arising from the uncertainty principle. So far, the only approximation used is that gravity is weak, φ 1. We anticipate one of the results of this paper, namely that the system can be studied in two different regimes depending on whether the axion field is |θ| 1 (the "dilute" and the "critical" axion star regimes) or |θ| 1 (the "dense" axion star regime). In the dilute and critical regimes, the axions comprising the star are nonrelativistic and the tools described in Sec. II B below apply. When |θ| ∼ 1, a full relativistic description is needed, as we sketch in Sec. III C.
B. Non-relativistic (single harmonic) limit When the non-relativistic limit applies, the axion mass is the largest energy scale in the problem, so that axion stars oscillate at a frequency very close to the axion mass m. Despite non-linear interactions arising from a cosine or a chiral potential precluding axion stars solutions from having one single frequency, for small field configurations |θ| 1, the one-frequency approximation suffices. Here, ω is the total energy of a constituent axion, in units of the axion mass. We write ω = 1+ , where accounts for the contribution from the binding, kinetic and self-interaction energies, while the one accounts for the rest mass energy. In the non-relativistic approximation, we have | | 1 and ω ≈ 1. We further assume that gravity is a weak effect, so that we can drop all terms containing φ in Eq. 7, except for the term 2φθ which is of the same order asθ + θ = 1 − ω 2 θ ≈ −2 θ. We split the potential into a mass term and the self interaction as Inserting the representation in Eq. (10) into Eqs. (7)- (9) and averaging over the period 2π/ω, we obtain ρ ρ kin +ρ grad +ρ pot .
In the last expressions, we have introduced the energy density terms and we have defined the effective self-interaction potential and its first derivative through For | | 1, Eq. (12) is a Schrödinger equation for the radial eigenfunction Θ with eigen-energy , while the energy density reduces toρ = Θ 2 /2 since the contributions from the gradient term and self-interactions are negligible.
We stress that our procedure, which involves the average over 2π/ω of the equation of motion leads to the same results as what was obtained in Ref. [54], where the authors neglect the rapidly oscillating terms proportional to powers of exp(iωt). As long as gravity is negligible and the single-harmonic approximation in Eq. (10) holds, Eqs. (12)- (14) are valid even for relativistic axions. We anticipate, that for (most of) the dense branch, gravity is indeed negligible but the single harmonic approximation no longer holds.

C. Axion potential
We expand the expression in Eq. (2) as In our numerical calculation, we truncate the sum in Eq. (18) to the first five terms h ≤ 5. This attains a precision below 1% with respect to the chiral potential in Eq. (2); this precision is better than the accuracy of the chiral perturbation theory itself. We slightly modify the coefficients v h so that the truncated potential shows: I) the same minimumṼ (0) = 0, II) the same mass V θθ = 1, and III) the same quartic couplingṼ θθθθ = λ φ as the full chiral potential in Eq. (2), where the (negative) quantity λ φ = −(1 − 3c z ) is related to the axion quartic self-interaction constant as λ = λ φ (m/f ) 2 . The numerical values of the corresponding corrected coefficients are given in Table I for z = 0.48. The effective non-relativistic potential in Eq. (16) is where J 0 (x) is the Bessel function of the first kind of order zero for the argument x. Notice that the cosine potential is recovered in the limit c z → 0, equivalent to setting v 0 = 1, v 1 = −1, and all other v h equal to zero in Eq. (18). The set of Eqs. (12)- (14) has been extensively applied to self-gravitating systems made of bosons. For the case of axions, the free case W 1 (Θ) = 0 has been studied in Refs. [23,55,56] following the seminal work in Refs. [17,18]. The potential expanded to the quartic interactions has been studied in Refs. [54,[57][58][59]. Ref. [27] considers the set of Eqs. (12)- (14) with the cosine potential, using the expression for the energy density (in our notation)ρ = Θ 2 /2, instead of our Eq. (15) obtained from the full energy-momentum tensor. This implicitly neglects contributions from self-interaction and kinetic energy to the energy density, which sources the gravitational potential. As we show below, those contributions to the energy density affect the results for the "dense" branch.

A. Axion star branches
We numerically solve for the radial profile Θ(r) appearing in the set of Eqs (12)- (14), as a function of the frequency ω. We impose the boundary conditions whereρ 0 is the rescaled energy density at r = 0 and the core amplitude Θ 0 is the amplitude of the axion field at r = 0. We obtain a radial profile Θ(r) via a shooting method, that is by varying the value of the core amplitude Θ 0 until we find a profile that decays as exp(−kr)/r at a sufficiently large r. The solution we seek shows no nodes, and corresponds to the lowest energy state for a given value of . See Ref. [60] for excited states of an axion star with a quartic potential. We find solutions for all values of ω within the range (0,1), although the numerics are particularly tricky as we approach ω = 0. For each value of ω, we obtain a unique value of the core amplitude and a unique profile. Given the radial profile, we obtain the total mass M = d 3 rρ and the radius R of the axion star, the latter defined as the radius containing 90 % of the energy [18]. In Fig. 1, we show the mass-radius relation for three values of f = {10 11 , 10 13 , 10 15 } GeV. 1 Each point on the line is characterized by a fixed value of ω and the core amplitude Θ 0 . For increasing value of Θ 0 , we identify three different regimes: the dilute branch (|Θ 0 | β 1/2 ), the unstable critical configurations (β 1/2 Θ 0 1), and the dense branch (Θ 0 1). For the critical line (red dashed line) and (most of) the dense branch (dashed black line), gravity is negligible. Then we find universal solutions when expressed in terms of the natural units of star mass, f 2 /m, and radius, 1/m. However, gravity is relevant in the dilute branch, where solutions depend on the value of f through β.

B. Non-relativistic solutions
In this section, we present heuristic arguments explaining the numerical results obtained in the previous Sec. III A for the dilute and critical branches where Θ 0 1; see also [57,58] for a similar approach. These branches can be understood in terms of the different contributions to the axion star energy U : the gravitational binding energy, the gradient energy, and the (quartic) self-interaction contribution, Here, α k and α 4 are dimensionless parameters which we insert to match the analytical results derived from Eq. (23) with the numerical solution. Estimating the mass of the axion star as  we can express the central amplitude as |Θ 0 | 2 ∼ M/(Λ 4 R 3 ), and the total energy U can be rewritten as In the last equality, we have used the scaling property of the Schrödinger-Poisson equation, writing the mass and the radius of the star in terms of dimensionless quantities, M = M (m/f 2 ) andR = mR. The natural scale for the mass and the radius of the axion star are then where M and R are respectively the mass and the radius of the Sun. The equilibrium configurations of the axion star can be qualitatively obtained by minimizing the energy density in Eq. (25) with respect toR, while fixing the axion star mass or, equivalently, the total number of axions N = M/m. This gives a quadratic equation whose solutions correspond to the radius of the star for either the dilute branch (R + ) or the critical branch (R − ), namelyR The stability of the solution is determined by the sign of ∂ 2 U/∂R 2 R=R± . Solutions in the dilute branch (ρ 0 β) are stable, while those in the critical branch (β ρ 0 1) are unstable. Matching onto our numerical results from section III A, we obtain α k = 9.9, independent of the value of β.
The dilute branch of the axion star corresponds to the equilibrium between the gradient energy and gravity. Depending on the value of the decay constant, equilibrium configurations of this type populate the line with negative slope in Fig. 1 with f = 10 11 GeV (blue), f = 10 13 GeV (green), or f = 10 15 GeV (orange), with the mass-radius relationR For configurations lying above this equilibrium line, the gravitational pull overcomes gradient pressure, so these configurations contract. On the contrary, configurations lying below the mass-radius line in Eq. (30) are restored to the equilibrium condition by the gradient pressure term. Hence, a restoring force acts to vanish any deviation from the stable equilibrium. The critical branch, the dashed red line in Fig. 1, corresponds to the balance of the gradient and the quartic selfinteraction energy contributions, with mass-radius rela-tionR Deviations from this configuration are pushed either further towards the dilute branch or to further contraction and are hence unstable. A solution for the radius of the axion star exists as long as the quantity below the square root in Eq. (28) is positive, that is when the mass of the star is smaller than the critical valuẽ which corresponds to the radiusR * = and to the core amplitudẽ The values ofM * andR * define the turning point in the top right corner of Fig. 1, corresponding to the transition from the dilute to the critical branch. In the critical branch, a denser solution corresponds to moving along the red dashed line in Fig. 1 towards the bottom left of the figure, with the star contracting and becoming lighter. Since in this branch the core amplitude increases as Θ 0 = Θ * 0 M * /M , non-perturbative dynamics becomes relevant when Θ 0 ≈ 1, or at a typical mass These values ofM (Θ 0 = 1) andR (Θ 0 = 1) mark the second turning point in the bottom-left region of Fig. 1.
For larger values of the core amplitude, the axion field explores the whole chiral potential and a different treatment is needed.

C. Non-perturbative solution
The axion star solutions found for Θ 0 1 correspond to a clump of axions whose total mass and radius are larger than the critical values in Eqs. (35) and (36). For such configurations, higher order terms in the attractive self-interacting potential cannot be neglected and a new regime is obtained, often referred to as the "dense" axion star regime in the recent literature [27,32]. We show the numerical results for the mass-radius relation obtained in the dense branch configuration with the solid black line in Fig. 1. Fitting the curve far from the turning point leads to the relationR = 0.6M 1/3 . This regime corresponds to classically stable configuration with an almost constant density ρ ∼ Λ 4 in the inner core. For the mass-radius relation, we have obtained the same power-law exponent (1/3) as in Ref. [27], because such dependence follows from the fact that the solution in the dense branch saturates the QCD potential and leads to a constant density of the star.
However, the structure of our solution differs greatly from what was obtained in Ref. [27]. We disagree on their interpretation of the equilibrium of the axion star in the dense branch for three main reasons. I) We have included the self-interactions and the gradient energy terms through Eq. (15). These terms cannot be neglected, as we show in Fig. 2. II) In Ref. [27] the set of equations is solved in the Thomas-Fermi approximation, that is neglecting the Laplacian of Θ appearing on the lefthand side of Eq. (12). III) Most importantly, the singleharmonic approximation in Eq. 10 does not hold in the non-perturbative regime.
In Fig. 2, we show the different contributions to the mass of the axion star, M = d 3 rρ, from the various components in Eq. (9), namely u α = d 3 rρ α /M , where α ∈ [kin, grad, pot], as a function of the core amplitude.
In the Θ 0 1 (Θ 0 1) regime shown, the star is in The frequency of the axion star ω (black solid line) as a function of the core amplitude Θ0 for our numerical solutions of the non-relativistic stability equations, (12)- (14). We also show the contributions to the total energy from the kinetic (blue dotted line), gradient (orange dashed line), and potential energy (red dot-dashed line). In the dense branch, i.e. Θ0 1, the solution is not consistent with the nonrelativistic approximation.
the critical (dense) branch. In the critical branch, the kinetic and potential energies both contribute a factor equal to 1/2. This result can be interpreted by the fact that the wave function of the coherent axion field undergoes harmonic oscillations, with the energy density equipartitioned between the kinetic and potential terms. However, as we approach the dense regime, the contribution from the gradient term increases, to the extent that for Θ 1 all three components contribute with a similar magnitude. Thus for dense axion stars the energy density must include all energy contributions. Also, the Thomas-Fermi approximation is not justified since the Laplacian term is crucial for solving Eq. (12) in the whole domain shown in Fig. 2 and 3.
In more detail, the structure of a dense axion star looks as follows. The stellar core is composed of relativistic axions since in that region ω 2 Θ ∼ ∇ 2 Θ, although selfinteractions are not entirely negligible. As we move out of the core, there is an intermediate region where the self-interactions balance the gradient term. Finally, in the outmost part self-interactions are again negligible .
To further illustrate that the axion field is relativistic in the dense regime, in Fig. 2 we show the axion energy per particle ω (black solid line), which drops to zero for Θ 1, due to the fact that self-interactions increase with Θ 0 . Then, the non-relativistic condition ω π/R, which expresses that the typical momentum of the axion is much smaller than its energy, no longer holds. Fig. 3 also shows this conclusion, since the quantity ωR decreases from being much larger than one to a constant value ∼ 3 for which the non-relativistic interpretation does no longer hold. The inequality mR 1, orR 1, which holds even in the dense branch, is not sufficient to justify a non-relativistic approach. In addition, our solution shows that gravity is negligible everywhere inside the star. The gravitational energy density at a distance r from the center of the star is ρ G = GρM (r)/r, where M (r) is the mass enclosed within the radius r, so we can write where in the last step we have used the parametriza-tionR = 0.6M 1/3 . Hence, gravity can be neglected for R 1/β. ForR = O(1), gravity can be safely neglected as long as f m Pl or β 1, which is the range of parameters considered in this work. However, for dense axion star solutions of larger mass, gravity could eventually become important again forR ≈ (4.6β) −1/2 . We do not consider this latter possibility here.
As we have previously discussed, the solutions obtained in the dense branch are not self-consistent because the single frequency approximation in Eq. (10) is not justified on the basis of the findings in Fig. 3. When the amplitude of the axion field becomes Θ = O(1), the axion fields probes the full chiral potential and all orders of self-interaction become relevant. Then, higher harmonic modes of the axion field whose frequency is a multiple of the fundamental mode ω = m are generated with amplitude comparable to that of the fundamental mode. In the next section we therefore start over from Eq. (7) and perform a multi-harmonic expansion.

A. Generalities on the relativistic equation
Based on the findings of the previous Section, axions in the dense regime Θ 0 O(1) can be studied using a relativistic approach and ignoring gravity. For simplicity, we derive results for the illustrative case of a cosine potential obtained from the chiral potential Eq. (2) for c z → 0.
In that case, the relativistic equation of motion is the Sine-Gordon equation We wish to identify the oscillon solutions of Eq. (39), namely the solutions that are spatially-localized and time-periodic. Such solutions circumvent Derrick's theorem [63], which states that the scalar field Lagrangian in Eq. (1) expressed in flat space-time does not admit time-independent, finite energy solutions because shrinking a non-zero field configuration effectively reduces the total energy of the system [64][65][66][67][68]. Although the ansatz we used previously, Eq. (10), is not a proper solution for the non-time averaged potential, we expect it to be a reasonable approximation at the transition from the non-relativistic to the relativistic domain when Θ 0 ∼ 1. There is a long history of searching for oscillons of the Sine-Gordon equation, with the most positive outcome being solutions that last O(100 − 1000) oscillations in units of 1/m [42,44,[69][70][71]. The general consensus is that absolutely stable solutions do not exist, although we know of no definite proof. In any case, is much that we can learn about unstable oscillons from the literature.
For axions in particular, Kolb and Tkachev [26] discovered the so called "axitons" when studying the cosmological evolution of the axion field in the dark matter context. They followed the evolution of the Sine-Gordon equation in an expanding Universe in which the axion mass strongly depends on the cosmic time, and identified an instability condition that leads to small clumps of the axion field with large values θ ∼ π to disappear in bursts of relativistic axions. This instability, which originates from the attractive quartic self-interaction term, is well known in the condensed matter community and has been recently revisited in Ref. [30]. In that paper, the authors follow the collapse of a dilute axion star with a mass slightly above the critical value M * . The axion star solution shows a self-similar collapse that ends when the central amplitude saturates the axion potential. Then, the axion field oscillates for a few times, radiating relativistic axions and relaxing to a small amplitude which is nevertheless larger than the starting value. Such instabilities are triggered for a few times until the central amplitude relaxes to the stability region described above. The simulations include gravity, so that the final state can still be a dilute axion star, but the dynamics of the collapse and the radiation of relativistic axions happens at very small radii where gravity is negligible compared with the self-interactions and gradients.
The simulations in Ref. [30] are of considerable phenomenological interest, since in principle the collapse of dilute stars is the most natural mechanism to produce dense axion stars. However, one can address the question of dense axion star stability separately from their possible cosmological origin. For such a task we need other means. A promising approach emerged in Ref. [72], where the authors convert the Sine-Gordon equation into a series of equations with different harmonics.

B. Beyond the 1st harmonic approximation
A general time-periodic solution can be written in terms of an infinite numerable set of harmonics. Thus we can write our oscillon ansatz as θ as [73] which, once plugged into the Sine-Gordon Eq. (39), yields a set of coupled equations for the different harmonics, Here, we have introduced the notation The set of Eq. (41) is a generalization of Eq. (10) when higher harmonics other than the fundamental mode ω are considered; when truncating the sum at n = 0 we obtain the single harmonic approximation Eq. (12) with Θ 1 ≡ Θ.
As an example, we consider the case where we also include the first term beyond the single-harmonic approximation besides the fundamental mode ω. This gives where we have approximated the computation of the coefficients I 1 and I 3 by expanding around Θ 3 = 0, with In fact, the solutions found in Sec. III C correspond to the zeroth-order approximation of the full non-linear solution, while solving the set of Eqs. (43)- (44) gives the next-to-leading order contribution. At r → ∞, solutions must approach zero, with Θ 1 , Θ 3 1. In this regime, the higher harmonic Θ 3 must satisfy which, in the range 1/3 < ω < 1, is an oscillatory solution in space with wavelength √ 9ω 2 − 1. Thus, if ω > 1/3, the axion star configuration radiates energy away through the third harmonic, with a contribution that increases with the total energy of the axion ω. Equation (47) can in principle be used to compute the lifetime of the axion star in the dilute branch, where the third harmonics is a very small perturbation of the exact solution. For values of ω < 1/3, there is no radiation solution at r → ∞ and higher harmonics have to be taken into account through Eq. (41). In Fig. 4, we solve the set in Eqs. (43)-(44) for a given frequency ω = 2π/T , with T = 7, using a shooting method to obtain the initial conditions for Θ 1 and Θ 3 at r = 0 that satisfy Θ 1 (+∞) = Θ 3 (+∞) = 0. We stress that the amplitudes of the 1st and 3rd harmonic are of the same order of magnitude everywhere in the star, demonstrating that the single harmonic approximation sufficient for the case of dilute axions stars does not suffice for the description of the dense regime.

V. DISCUSSION
We have shown that when Θ 0 O(1), axions are relativistic and axion stars enter the dense branch regime where the configuration behaves as a meta-stable oscillon of the χ-Gordon equation, with a characteristic lifetime. For a free field, bosons stream away from the oscillon core of a star of radius R, with a lifetime τ lin = 0.836 √ 2R 2 [74] and with a radiation spectrum peaking at ω lin with width Γ lin = (2τ lin ) −1 . Including a more realistic non-linear self-interaction potential modifies the spectrum by lowering the peak frequency at a lower value ω nl < ω lin , with a new width Γ nl < Γ lin . Following Ref. [74], an oscillon forms if the two spectra do not significantly overlap, that is when The computation of the oscillon lifetime for a quartic selfinteraction has been addressed in Refs. [74,75], where the relatively long lifetime (on the scale of the intrinsic timescale m −1 ) of oscillons is explained by the relatively small overlap between the oscillation frequencies. Following this method, we estimate of the lifetime of an oscillon for a cosine potential as where we have used the parameters α = 5 × 10 −5 , E osc = 402.1, and E ∞ = 372.8, following Refs. [74,75] with the axion Lagrangian in Eq. (1) and a Gaussian ansatz for the radial wave function. In short, the energy of an oscillon is described by its radius and amplitude, and damped oscillations in the oscillon develop along the line of constant minimum energy [42,43]. We performed an independent check of these results by using the solutions of the time-independent Eq. (12) in the dense branch as initial conditions which we time-evolve with the Sine-Gordon Eq. (39), as prescribed in Ref. [72]. Although this initial wave function is not a proper solution to the Sine-Gordon equation, our numerical solutions yield breather solutions. For a period T nl = 7.0 as considered above, we find that the solution decays after τ life ≈ 1200/m. This result is of the same order of magnitude as what we obtained using Eq. (49) The fact that pseudo-breathers exist has been shown in Ref. [46], where the existence of a finite life-time solution to the Sine-Gordon equation has been related to the singular behavior of the solution at zero, when an oscillating function has been imposed as the boundary condition of the solution at infinity. Pseudo-breathers are ultimately decaying states, as discussed in length in Ref. [26], where it is found numerically that such solutions are unstable and fragment into smaller clumps. The dynamics and the initial conditions considered in Ref. [26] are however different from ours, since the authors consider a cosmological evolution of the axion field with "white noise" initial conditions, and included the Hubble rate in the equation of motion.
Recently, Helfer et al. [31] have studied the stability of axion stars including gravity and non-linear effects, finding that stable dense profiles may be possible when f 0.1 M Pl , the exact value depending on the axion star mass. In any case, the energy scale f involved is well above the scales we consider here. For values of f below this critical value, the axion star either collapses to a black hole or dissolves by the emission of relativistic particles, consistently with the puffing out obtained in Ref. [30] and in this work.

VI. CONCLUSIONS
In this paper, we have discussed the properties of axion stars for all allowed values of the core amplitude of the axion field Θ 0 . In particular, we have discussed how classically stable solutions can arise from the interplay between self-gravity, axion self-interactions, the pressure due to the Heisenberg uncertainty principle, and the kinetic energy. Using assumptions commonly made in the literature, we have obtained a set of equations describing coherent axion field oscillations inside the axion star in the single-harmonic approximation. For small core amplitudes Θ 0 1, we confirmed known results for axion stars in the dilute and critical branches, and provided a heuristic interpretation of those results from first principles.
For Θ 1, the "dense" regime, we recover similar results to those in Ref. [27] when using the single-harmonic approximation, in particular, the mass radius relation R ∝ M 1/3 . However, we argue that the single-harmonic approximation does not hold for the dense regime and thus a different approach is needed, taking into account higher harmonics. In the end, we arrive a very different physical interpretation of the dense regime. We find gravity to be negligible for Θ = O(1). Dense axion stars should be solutions to the Sine-Gordon (or χ-Gordon) equation describing the axion field inside the star. We computed the lifetime of dense configurations using both the semi-analytical procedure described in [42,43,74,75] and by using our single-harmonic solutions as initial conditions, which we time-evolved numerically using the Sine-Gordon equation as prescribed in [72,73]. Both methods yield comparable lifetimes of order τ life ∼ 10 3 /m, much shorter than any cosmological time scale.
We conclude that if dense axion stars can be formed, they would immediately (on cosmological scales) radiate relativistic axions and decay. Since axion stars in the critical branch are unstable against perturbations and either expand to stable dilute configurations or contract to the dense branch and subsequently decay, stable axion stars with mass M >M * (f 2 /m) ∼ 10 −11 M 10 −5 eV/m 2 appear implausible.

ADDITIONAL NOTE
During the final preparation of the manuscript after completion of this work we received [76,77], partially overlapping with this work.