Recurrent Axinovae and their Cosmological Constraints

Axion-like dark matter whose symmetry breaking occurs after the end of inflation predicts enhanced primordial density fluctuations at small scales. This leads to dense axion minihalos (or miniclusters) forming early in the history of the Universe. Condensation of axions in the minihalos leads to the formation and subsequent growth of axion stars at the cores of these halos. If, like the QCD axion, the axion-like particle has attractive self-interactions there is a maximal mass for these stars, above which the star rapidly shrinks and converts an $\mathcal{O}(1)$ fraction of its mass into unbound relativistic axions. This process would leave a similar (although in principle distinct) signature in cosmological observables as a decaying dark matter fraction, and thus is strongly constrained. We place new limits on the properties of axion-like particles that are independent of their non-gravitational couplings to the standard model.


I. INTRODUCTION
The stars in our galaxy are made of ∼ 10 57 fermions bound together by gravity and protected from collapse by thermal pressure or fermion degeneracy pressure.In the presence of a light, long-lived boson similar gravitationally bound states of that boson may exist, but in the absence of nuclear burning they are instead supported by gradient pressure, which is a result of the uncertainty principle.Axion stars are one such example of these bosonic objects.
In the Standard Model, stars convert approximately 0.1% of their mass energy into radiation over their lifetime.The small energy released (compared to rest mass) in the pp-chain, for instance, is due to the relatively small binding energy inside the star.Only very compact objects like neutron stars are relativistic in nature.Moreover, because of baryon number conservation, there is a limitation on overall energy release given the (approximately degenerate) neutrons and protons which must remain in the final state.
However, in the dark sector, there are reasons to expect the overall energy conversion could be much higher if a similar process were to occur.Complete conversion of rest mass from e.g. a 3 → 1 process is possible because there is no "baryon-number" conservation for bosonic dark matter.For example, axion stars will collapse and emit relativistic axions when they reach a critical mass.We call such processes that drastically convert dark matter to dark radiation as Axinovae.There is no mechanism to quench the axionovae if axion stars form ubiquitously in the Universe, as expected in the post-inflationary scenario where axion miniclusters form after matter-radiation equality.Therefore, a large formation rate of axion stars that lead to axinovae is very constraining.We take the formation of axion stars as a concrete example to study but the result can apply to generic scalars whose self-interaction is attractive since the properties of axion stars do not depend on any inter-actions other than gravity and the axion self-couplings.
A natural cosmic history that can occur generically for these models is, after matter-radiation equality, these axions stars form, grow, and finally explode as an axinova, converting a significant fraction of energy into semirelativistic axions.After this, the remnant can continue to grow, until it explodes again.This process of recurrent axinovae can convert a significant fraction of the dark matter into relativistic energy, which is then constrained by cosmological observations.This paper is organized as follows: in Sec.II, we discuss the formation of enhanced structures at small scales due to the axion perturbations and study the formation history of axion stars inside those structures.In Sec.III, we study the constraints on axion parameter space by requiring the decay fraction of axion dark matter should not exceed an upper bound.In Sec.IV, we present our conclusions.

II. AXIONS, AXION MINIHALOS, AND AXION STARS
The axion is a well-motivated dark matter candidate, which can also leave unique fingerprint on the matter power spectrum at small scales if the PQ symmetry breaking occurs after inflation.In such scenarios, different horizon patches have different matter densities when the axion acquires its mass, leading to the formation of axion miniclusters or axion minihalos at matter-radiation equality [1,2].More interestingly, coherent objects called axion stars can form in the center of axion minihalos due to Bose-Einstein condensation [3], which may eventually accrete into a critical object and emit relativistic axions.We call such phenomenon axinovae, which can occur with an attractive axion self-coupling and the formation of axion minihalos at matter-radiation equality.
Originally proposed to solve the strong CP problem [4][5][6], the present-day landscape of axions and axion-like particles (ALPs) is broad.One common feature across this landscape is that the axion, φ, is a pseudo-Goldstone boson of a global U (1) P Q symmetry broken at a scale f a .The U (1) P Q is anomalous under a confining gauge group which means that the axion's potential is generated through instanton effects occurring at the compositeness scale of the gauge group, Λ, and takes the form In the case of the QCD axion Λ ≈ 200 MeV and c ud ≈ m u m d /(m u + m d ) 2 ≈ 0.2.In addition to the self couplings, coupling to gravity, and the anomaly-induced coupling to QCD (or QCD-like group), the axion may have model-dependent couplings to other SM gauge bosons and fermions.We focus here on the self and gravitational couplings only, which can already lead to interesting dynamics such as axinovae.

A. Axion Minihalos
In the post-inflationary scenario, the present-day Universe contains a large number of patches which were causally disconnected at the time of QCD phase transition.In each causally disconnected patch of the Universe, axion field values are uncorrelated.Once the axion acquires a mass, and Hubble friction is small enough, the axion behaves as cold dark matter and isocurvature fluctuations are present in the matter density.When the Universe becomes matter dominated this small-scale structure will start to collapse under gravity, leading to axion minihalos.Furthermore, there may be large overdensities of axions at even smaller scales arising from the evolution of the network of axion strings and domain walls [7] set up when the PQ symmetry breaks.Even for the much studied case of the QCD axion, there is controversy [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] as to what fraction of the relic dark matter axions arise from misalignment or from the decay of topological defects.Along with those topological defects, objects called oscillons or axitons that can contribute to the small scale overdensities will form after the axion accquires its mass [14,15].Those objects can form when the axion self-interaction dominates over the Hubble expansion term, which is easily satisfied in the early Universe when the self-interaction is strong due to the high density.As the axion density drops, the formation of oscillons will be turned off and oscillons themselves will dissipate via emitting relativistic axions.
It is worth noting that the post-inflationary scenario is not essential for the axinovae.Any matter power spectrum which is enhanced at small scales can lead to the formation of axion minihalos around matter-radiation equality, but the post-inflationary scenario is a minimal realisation.We take a simple ansatz for the spectrum of initial fluctuations in the axion field, namely that the spectrum of isocurvature fluctuations in the axion field follow a white-noise spectrum, cut off at small scales i.e.
Here k 0 ≈ a osc H osc is the (comoving) wavenumber determined by the horizon size at the time the axion starts to oscillate, i.e. m a (T osc ) ∼ 3H osc .While here we consider a pure white noise spectrum we extend this analysis to a more general power law spectrum in Appendix C. In reality one would expect a softening of the cutoff in the white noise power spectrum at small scales.The exact details of how this occurs is related to the dynamics of string network and axitons, and is unknown.It will not affect our conclusions, see Appendix C for details.As mentioned above, the contribution of strings and domain walls to the abundance of non-relativistic axions is uncertain and will impact the size of the power spectrum.Simulations typically show the density perturbations have A 0 ∼ 0.1 at k = a osc H osc , but they also show larger sub-horizon (larger k) fluctuations.These sub-horizon fluctuations can collapse earlier than those at the horizon scale, leading to high concentration minihalos.These halos are at smaller scales, k ≥ a osc H osc , and have larger δρ a /ρ a and lower axion speeds in the mini-halos, resulting in a faster star growth rate.We take the conservative limit of holding k 0 = H osc as the scale at which A 0 = 1.With this white noise power spectrum (2) the first structures, of mass M 0 , form at redshift z c ≈ A 1/2 0 z eq and the characteristic structure mass, defined as the peak in the distribution M dn/d log M , occurs at Where is the co-moving mass in the horizon at the time the axion starts to roll and ρ 0 is the present-day cosmological axion density.The minihalos have a distribution of masses but for simplicity we use the characteristic mass M peak to provide a measure of the overall behavior.The growth continues till around z ∼ 10 − 20 when the minihalos merge into standard CDM halos and their growth stalls [23].
We take the minihalos to have an NFW [24] density profile, defined by a scale radius r s and density ρ s , At the scale radius the circular speed is given by This speed will be relevant for the calculation of axion star formation rate, and in the minihalos that will form axion stars this speed is much smaller than typical speeds in the Milky Way.Numerical studies have shown that the dark matter halos at the characteristic mass will first reach a concentration factor [25] c ≈ 4 at the time of halo collapse [26].For lighter halos that grow more through accretion than mergers, the halo concentration will grow linearly with the scale factor due to the decreasing background density.In the next subsections we discuss further structures that can develop in the core of these minihalos.In addition to the low speeds in minihalos, the scale density in the mini-halos that form early is large, where ∆ 200 ≈ 200 in the spherical collapse model.For mini-halos that collapse at z eq this density is ρ s ≈ 10 14 ρ 0 .

B. Lifecycle of an Axion Star
We now turn to the question of formation of axion stars at the core of the minihalos discussed above.The subsequent growth and explosion of axion stars (axinovae) will also be studied.There are a few timescales we will discuss that are relevant in the lifecycle of axion stars: • The condensation timescale from gravitational interactions • The condensation timescale from axion selfinteractions.
• The evaporation timescale of light axion stars.
• The Hubble time when the axion star formation is active.
We will discuss those timescales later in this subsection.
For the parameter space that axinovae can place meaningful constraints on, the axion self-interaction always dominates over gravity.As discussed in Appendix A, there are two branches of axion star configurations: the dilute branch, which, below a certain mass, is stable and the dense branch which will explode and emit relativistic axions.There is a critical star mass (A5) that separates the two branches, which we denote as M max * .Therefore, if they continue to accrete mass, the lighter dilute axion stars will eventually become unstable in a minihalo environment.Axion minihalos are ideal environments for the axion star formation because they are dense and cold, owing to the high collapse redshifts and small virial masses.When the star formation rate is sufficiently large, stars will form in the minihalo center and grow to a critical mass star if the minihalo is massive enough.The critical star will contract under self-interaction and gravity, converting a large fraction of its mass to relativistic axions.Until the axinovae consume most of the minihalo mass, axion minihalos remain ideal environments for the axion star formation and axinovae shall occur again within the same timescale.Thus, we naturally expect the axinovae phenomenon to be recurrent, when the growth timescale is fast enough.The crucial calculation to determine the fate of axion stars is the formation rate in the minihalo center and the corresponding star mass.
Once minihalos exist, gravitational interactions or selfinteractions can subsequently lead to the formation of Bose-Einstein condensed axion stars at their center.The timescale for this formation, and subsequent growth, in an environment where the axions being captured have typical number density n and speed v is determined [27][28][29][30] by With σ the total scattering cross section.This formation rate is Bose-enhanced from the naive expectation due to the large phase space density, The gravitational Rutherford transport cross section is , where the Coulomb logarithm has been cutoff at a characteristic length scale of the minihalo, R. Attractive self-couplings can also lead to formation and the scattering cross section is . The total condensation time, considering both gravity and self-interaction, is With each individual process having a timescale of for gravity, and The parameters b, d ∼ O(1) are numerical coefficients that are extracted from numerical simulations [31].Comparing these two timescales, ( 10) and ( 11), we see that the self interactions will determine the axion star formation rate if f a < ∼ M pl v. Furthermore, if the relevant speed is determined by gravitational collapse of a minihalo (6) then self interactions dominate in the limit f 2 a < ∼ ρ s r 2 s .When determining the gravitational relaxation timescale for formation of axion stars in minihalos we take, as typical, the densities and speeds at the scale radius, see Eqs. ( 5) and (6).
In addition to the timescale for axion star growth there is also a rate for evaporation of the star.Axions in the halo that are not part of the star can collide with bound axions causing them to be ejected.The rate for this process shrinks with axion star mass and is approximately [32] Γ evap ∼ (m a vR * ) 2 τ −1 .The competition between growth and evaporation means only axion stars above a certain mass will gain mass by gathering axions from the halo.As observed in numerical simulations [27,32] such stars first appear after time τ and then proceed to grow.The growth is initially fast (d log M/dt is constant) but once the virial velocity of the minicluster falls below the speed of the axions in the axion star the rate of growth slows, d log M/dt becomes inversely proportional to (a power of) the star mass [32] which results in the mass growing with time as a power law.The characteristic axion star mass where this change in behavior occurs is obtained by equating the virial velocity of the minicluster with that of the axion star [27,28,[31][32][33][34] is where M h is the halo mass.The behavior of the growth rate once the axions in the star are moving faster than those in the halo is not definitively known, and there is evidence that it may continue to evolve with star mass [32].This would result in the mass growing as a power law with a running index.However, to simplify our analysis and to partially account for the numerical uncertainties, we will use a single power law but consider a range of possible powers.In particular, we parametrize the powerlaw mass growth as M * = M * (t/τ ) 1/α and vary α in the range of 1 to 5. With initial exponential growth followed by constant power law growth, the timescale to form an axion star at critical mass M max * depends in which regime the critical mass falls.Thus, (13) The numerical simulations discussed above have mostly been carried out assuming a homogenous gas of axions as the initial background upon which an axion star forms.For stars that form in minihalos the gas has a density and velocity profile.In Appendix B we argue that for an NFW profile the exponential growth is replaced with a power law, and the whole growth becomes a single power law, with α = 3/2 when self interactions dominate.
Given that the majority of the dark matter has collapsed into axion minihalos with a characteristic mass M peak (z), the total fraction of dark matter rest mass that has been converted to kinetic energy per unit time can be calculated as where κ is the fraction of the axion star's mass that is converted to relativistic axions during axinovae.From simulations of these processes [35], it is seen that approximately 50% of the star's mass is lost during the nova and of this about 20% is in the form of relativistic axions, so κ ≈ 0.1.The time to reach a critical star given in (13) assumes the star grows from an undistorted minihalo.After the first axinova there is a remanent of mass ∼ 0.5M max * already present and the time for this to grow to M max * is slightly shorter than for the first star.For the power law considered here this correction is small and we ignore it, assuming all subsequent stars take time t crit to explode.

III. COSMOLOGICAL CONSTRAINTS
A. The decay rate of axion stars The process of forming axion stars which subsequently become nova converts non-relativistic dark matter axions into boosted (γ ∼ O(few)) axions.The kinetic energy of the outgoing axions will red-shift away after the scale factor has grown by ∼ √ γ and thus the dark matter's contribution to the matter-energy budget is depleted.Here we study the impact of the cumulative loss of mass in the dark sector but it is possible that the temporary existence of a new relativistic species may lead to a measurable effect on large-scale structure and is worthy of future study.This process is closely related to the scenario of decaying dark matter, which is well constrained by recent cosmological data [36][37][38].For dark matter which decays after recombination, the decrease of the dark matter fraction will increase the angular diameter distance to the last scattering surface over time.Furthermore, the amount of CMB lensing is reduced due to a smaller gravitational potential than expected.This scenario is constrained by a combination of CMB [39] and, for very long lived dark matter, SDSS [40] data.If the decay of dark matter occurs well before recombination or even before matter-radiation equality, the primary effect of the decaying dark matter is to enhance N eff since the decay products behave as dark radiation.In the short-lived situation the constraints are primarily from CMB measurements.We will be interested in the long-lived case, and in particular decays which occur after matter-radiation equality but are no longer ongoing.The equivalent bound [38] for decaying dark matter on the fraction of the initial amount of dark matter that can decay is Although the cosmological evolution of the dark sectors for decaying dark matter and axinovae are not identical they are similar and since the above constraint is independent of decaying dark matter lifetime over a wide range of lifetimes we will use it to constrain axions.We leave a more detailed numerical analysis, and an investigation of other possible signals, for future work.Converting (15) to the case of axinova leads to the requirement that In the scenario of axinovae, the decay of dark matter occurs when axion miniclusters start to form, which is always after matter-radiation equality.To avoid the constraint of ( 16) requires either that the formation rate of axion stars is too small to be cosmologically relevant or that the formed axion star mass is smaller than the critical mass so there are no axinovae.Note that this bound does not rely upon there being a coupling to any SM particles e.g.photons, gluons, or SM fermions.However, our constraints do rely on the assumption that axions do make up the dark matter relic abundance and that the fluctuations in the axion field are isocurvature in nature and approximately power law.If axions make up a fraction of the dark matter, this can be encoded as a reduction of κ, see (14), and a corresponding weakening of the bounds.
For normal misalignment production of axions, where θ 2 ≈ 4 [41], the typical initial halos that form have a mass that depends upon the horizon size when the axion starts to oscillate m a (T osc ) = 3H(T osc )/2.For the QCD axion, where the temperature dependence of the axion mass is known, this oscillation time is uniquely determined.However, in more general axion scenarios the oscillation temperature, and therefore M 0 , is a free parameter.In the radiation dominated era H(T ) = π (8πg * (T )/90) 1/2 T 2 /M P l and the halos form with mass, The existence of DM structure down to small scales requires that the axions behave as dark matter by the time the temperature of the Universe is ∼ keV, i.e.T osc > ∼ 1 keV.Thus, there is an upper bound on the initial halo mass.More sophisticated analysis of the constraints on the axion isocurvature power spectrum at small scales can be found in Ref. [42].
Going forward we will assume that the axion makes up a sizable fraction of the dark matter abundance and place a bound on its self-coupling, equivalently f a , through recurrent axinova.There are four parameters that determine the amount of axion dark matter that is converted to dark radiation: the axion mass m a , the axion self coupling λ which in simple models is determined by the decay constant f a , the structure mass M 0 (or equivalently M peak (z c )), and the red-shift at which minihalos first form z c .Numerical simulations [14,15] indicate that the white noise spectrum has large amplitude at small scales A 0 ∼ O(1) and thus minihalos form as early as possible z c ∼ z eq , with mass given by (17).
As times evolves, the characteristic mass grows as M h ∼ (1 + z) −2 as minihalos merge with each other.Since a characteristic mass halo has concentration c ≈ 4 its scale radius and density vary with redshift as r s ∼ (1 + z) −5/3 , ρ s ∼ (1 + z) 3 and consequently the speed at the scale radius depends on redshift as v s ∼ (1 + z) −1/6 .From Eqs. (10) and (11) this implies that the time scales for collapse scale as This rapid lengthening of the axion star formation time as the Universe ages means that the dominant DM mass loss occurs as soon as the minihalo mass is larger than the critical star mass, and the earlier that occurs the greater the fraction lost.More precisely, assuming t crit is in the power law regime, the decay rate for halos of mass M 0 which initially form at redshift z c is, where we have suppressed the logarithmic corrections to the Rutherford cross section in (10), taken b = d = 1, and ρ col = (1 + z c ) 3 ρ 0 is the background density at the time of initial collapse.The from of ( 18) makes clear that the rate is peaked to early redshift and this rate is enhanced by decreasing both m a and f a .If the timescale for scattering is set by self interactions, i.e. f a < ∼ M the decay rate is constant.Furthermore, for any choice of parameters there is a maximal f a above which there is not enough time to form a critical mass star in a minihalo.This leads to a region, bounded from below (above), in m a − f a (m a − f −1 a ) space which is constrained by the cosmological data discussed above (14).
In Fig. 1, we plot the region that is constrained by the axinovae, for various assumptions.The gray regions are excluded by black hole superradiance constraints [43][44][45][46].The most conservative (weakest) constraint, shown in green, comes from assuming that the oscillation temperature is low and that the time to reach a critical star is given by (13).Over most of the green region the critical star mass is low and the growth (d log M/dt) is still in the constant regime.Given constraints on large scale structure we take the lowest possible oscillation temperature to ∼ 1 keV.The later an axion starts oscillating the larger the mass of the initial axion miniclusters, which leads to a longer axion star production time τ , suppressing the resulting appearance of axinova.
In the red region we again assume the lowest possible oscillation but now assume that the star growth is power law, M ∼ ρ s r 3 s (t/τ ) 2/3 , for all star masses, as discussed in Appendix B. At masses below ρ s r 3 s the power law predicts faster growth than the constant growth assumed in (13) and the green region.This makes the bound stronger.For f a > ∼ 10 15 GeV the axion star critical mass is larger than where exponential growth transitions to power law in the green and the two constraints coincide.
Finally, the blue region is the strongest constraint and is found by optimizing over the oscillation temper-FIG. 1.The exclusion region from axinovae for different assumptions for axion parameters, see text for more details.Existing limits from the black hole superradiance are shown in grey.The green region is the most conservative bound using a constant rate of dlogM/dt at M max * ≤ M * and a late Tosc, with a formation timescale given in Eq. ( 13) The red region uses a power-law growth (PL) with M ∝ t 2/3 over all the mass ranges and it also assumes the lowest oscillation temperature.The blue region presents the bound after optimizing over oscillation temperature.
ature.The maximum possible oscillation temperature arises when the axion starts oscillating with its zero temperature mass, T osc = m a M pl .These high temperatures will lead to the lightest axion miniclusters and the shortest star production times, but such miniclusters may not be massive enough to contain a critical star.At each point in the parameter space, we select the highest possible T osc that leads to a massive enough minicluster.Since, M h ∼ T −3 osc this selected temperature is still close to m a M pl .In Fig. 2, in Appendix D we show the constraint for M * ∼ t 1/5 , when the leading order the decay rate is independent of T osc .
In the excluded regions an O(1) fraction of all dark matter has passed through an axionova.This may lead to other observables in axion experiments or in cosmological observations.Given the high powers that appear in (18) if the constraints on decaying dark matter are improved in the future the region of parameter space excluded will not be greatly altered.

IV. CONCLUSIONS
We obtain new bounds on axion dark matter parameters m a , f a assuming the formation of dense axion minihalos, motivated by the post-inflationary scenario.Axion perturbations in the post-inflationary scenario will lead to the formation of dense substructures known as axion miniclusters or minihalos after matter-radiation equality, which can subsequently form coherent objects known as axion stars at the core of axion minihalos.Low mass dilute axion stars, supported by gradient pressure, can be cosmologically stable.However, they will accrete more axions from minihalos and continue to grow in mass until the axion self-coupling becomes important and the gradient pressure can no longer stop them from collapsing and emitting relativistic axions, in an axinova.The remnant of an axinova is a less massive star which will again grow, leading to recurrent axinova.
If the recurrent formation rate is large enough and axinovae are active, they can convert a significant fraction of dark matter into radiation which can be constrained by measurements of large scale structure formation.Our constraint only depends on the axion self-coupling and gravity.The self coupling can be mapped to axionphoton and axion-neutron couplings in specific models.Those constraints are obtained by requiring the population of dense axion stars formed in axion minihalos at high redshifts shall not dominate the mass of dark matter.If the axion is only a fraction of dark matter or only a few percent of axion dark matter is decaying, the conversion to dark radiation may be cosmologically significant in future observations but consistent with the current data.Alternatively, if the axinova has a branching fraction into standard model states there may be observables in the region or parameter space close to our bound.We leave a more detailed study of the cosmological evolution or possible visible signals to future work.smaller velocity.If a star formed within a small radius, the mass contained in this region is small.Therefore, lighter objects always start to form with a greater rate.For an NFW profile, the mass contained within r is The formation timescale given by self-interactions is τ self ∝ v 2 /ρ 2 .At small radius of an NFW halo, the density and velocity scale as ρ ∝ 1/r and v ∝ √ r.Therefore, M (t) ∝ t 2/3 .Similarly, if the gravity dominates the axion star formation, τ gr ∝ v 6 /ρ 2 ∝ r 5 at small radius and we obtain the mass growth power law M (t) ∝ t 2/5 .Since this scaling is active at short distance scales within the minihalo we consider a scenario where α = 3/2 at all axion star masses, see Fig. 1.

Appendix C: Press-Schechter with White Noise-like Power at Short Distances
We consider the density perturbations, δ ≡ δρ/ρ, to consist of two contributions, conventional ΛCDM adiabatic perturbations that are present at all scales and isocurvature perturbations which are only become important over a finite range of scales.We take the isocurvature contribution to be a power low with a cut-off at very small scales, corresponding to a wavenumber k 0 .For the case of the axion it is believed the short-scale behavior has a power spectrum that is approximately that of white noise, corresponding to n = 3 below.Modes from these two contributions have different growth behaviors after they enter the horizon, in particular the adiabatic perturbations have logarithmic growth until matter-radiation equality while the isocurvature modes do not.At late times, in the matter dominated era, they have similar growth.Taking into account these different growth behaviors the two-point function of the density perturbations is For a ΛCDM-like power spectrum A s ≈ 2 × 10 −9 , n s ≈ 0.97, and the pivot scale is k s = 5 × 10 −3 Mpc −1 .At late times D adi ≈ D iso ≈ a/a eq = (1 + z eq )/(1 + z) and the exact forms can be found in standard references e.g.[59,60].The constant I 1 ≈ 9.1 and L ≈ log(0.1aeq /a).
The Press-Schechter formalism assumes spherical collapse of over-densities and that the probability for these collapses follows a Gaussian distribution whose variance, smoothed at some scale R, is given by where W (kR) is the window function and can take various forms.Here we focus on the so-called sharp k-filter where W (z) = Θ(1 − z).For this choice of window function there is not a well defined mass, M , associated with the co-moving filter scale R, since the real space form of W does not have local support [61].However, we will follow the oft-used relation M = 6π 2 ρ 0 R 3 [62], where ρ 0 is the present day cosmological axion density.Note that for (C2) to be well defined we have to introduce an IR cut-off k IR and we define M 0 = 6π 2 ρ 0 k −3 0 .We are typically interested in halo masses and formation redshifts where the adiabatic perturbations are subdominant to the isocurvature perturbations, A s A 0 .In this regime, once structures can form i.e. z < z eq , the variance has the simple form In the Press-Schechter approach the halo mass function is related to the probability to find δ > δ c ≈ 1.686, with the fraction of matter in objects of mass M given by The exponential suppression means that the most massive objects, with mass M peak , to have formed are those for which σ(z, M peak ) = δ c .If the isocurvature perturbations were large enough, A 0 > nδ c , these objects would form at z eq .Instead, for more typical isocurvature perturbations of A 0 ≈ 0.1, the first halos to form are of mass M 0 and they form at and subsequently grow, with the peak mass of the halo mass function being Appendix D: Axion Relic Abundance from Misalignment We consider the relic abundance from the misalignment mechanism for an axion coupled to a dark confining gauge group "DarkQCD", which is taken to be SU (N C ) with N F vector-like quarks.The temperature dependence of the mass is understood in two limits.At low temperature the axion mass is independent of temperature and at high temperature the dilute instanton gas approximation is valid, leading to a power law dependence.In between there could be a first or second order transition or a smooth cross over depending on N F , N C [63,64].For simplicity we take the temperature dependence mass to have the form Here we take the critical temperature to be the same as the confinement scale of DarkQCD, T c = Λ = √ m 0 f a .The dilute instanton gas approximation gives b = (11N C + N F − 12)/6.Taking b large for temperatures in the vicinity of T c also approximates the form of a first order phase transition.After PQ symmetry breaking, and before the instantons generate a potential for the axion, the misalignment angle θ = a/f has a flat potential and is free to take on any initial value in each causal patch.The equation of motion for this angle is θ + 3H θ + m 2 (T )θ = 0 . (D2) Assuming the cosmology is governed by a fluid with equation of state p = ωρ (RD is ω = 1/3) then the scale factor a ∼ t 2 3(1+ω) and H = This equation can be solved exactly by noting that y = x α J n (βx γ ) with J n the n-th Bessel function, satisfies the equation Thus, the solution to (D3) takes the form . (D5) Requiring that the argument of the Bessel function changes by an O(1) amount before oscillation is deemed to have set in, and identifying various powers of t with H and m a (T ), the oscillation temperature is implicitly defined by Notice that for large b, an axion mass that rapidly changes from zero to m 0 as can arise in a first order phase transition, this is different from the usual m ∼ 3H/2 requirement since the rapid evolution of the axion mass provides its own "friction".From now on we consider the case of RD and thus m osc ∼ (2+b)H osc .We also consider the possibility that the dark sector and the are at different temperatures.Assuming there are no thermalizing interactions between them, and ignoring the complication of different thresholds in the two sectors we take the ratio of temperatures to be a constant, T D = ξT SM ≡ ξT .Thus, the oscillation temperature and mass are found by solving 8π 3 g * (T ) 90 If the oscillation begins while the mass is temperature dependent then where c(T ) = 90/8π 3 g * (T ) and, assuming the SM dominates the energy density of the Universe, 0.06 < ∼ c(T ) < ∼ 0.33.This solution is only consistent if T D > Λ which places the restriction b < ∼ ξ 2 M pl /f a [65].For b, ξ in violation of this bound the oscillation starts after the axion has attained its zero-temperature mass and Once the oscillation temperature is known, and using the fact that ratio of axion number density to entropy density is constant, the present day axion mass fraction can be determined: If the dark sector has roughly the same temperature as the standard model sector (ξ ∼ 1), the confinement scale corresponds to a Hubble of H ∼ m a f a /M pl , which is always smaller than m a because we require f a < M pl and the axion self-coupling is stronger than gravity.The axion mass will not be turned on until the dark confinement occurs.Therefore, T osc is greatly delayed, which enhances the relic abundance since it is less diluted.The blue dashed curve in Fig. 2 shows the axion parameters that give the dark matter relic abundance assuming a slightly colder dark sector (T DS = 0.5T SM ) and axion mass to be turned on as m a ∝ T −b .A large dark gauge group or a first-order phase transition in the dark sector will be needed for a large b.We also presented the T osc independent constraint in Fig. 2 which assumes axion star mass grows like M ∝ t 0.2 , corresponding to α = 5.For this value of α the decay rate ( 18) is independent of T osc in the region of parameter space dominated by self interactions.
While we have been focusing on a QCD-like axion model to study the relic abundance, there are other models that can enhance the self-coupling of axions while giving the correct relic abundance, such as a clockwork axion [66] (discussed in Appendix E), friendship axion [67], axions from dilute domain walls [68,69],and kinetic misalignment mechanism [70,71].In a clockwork axion scenario, a large field range is naturally produced for the axion field in the low-energy theory.The axion potential can have two confinement scales and two effective decay constants which can give the relic abundance that is needed while keeping the self-coupling strong.The friendship axion can resonantly convert the energy density in the axion sector with a larger decay constant to that with a lower decay constant if the mass ratio of two axions is close to 1. Therefore the relic density of axions with a low decay constant is greatly enhanced.Axion relic density can be greatly enhanced if the Peccei Quinn symmetry is followed by a period of inflation such that axion string networks are inflated away but will eventually reenter the horizon [68,69].In this scenario, the decay of diluted domain walls occurs very late, enhancing the relic density of axions.In kinetic misalignment [70,71] the axion field does not start at rest but instead has a nonzero initial velocity.The process of the axion settling into a minimum of the periodic potential, and generating an axion number density, is delayed since it can only occur after its initial kinetic energy has red-shifted away.The initial velocity, θi , for the field is proportional to the net PQ charge and its generation requires an explicit breaking of the PQ symmetry at some scale.This breaking should not be present at later times when the axion potential should be determined solely by instanton effects as can occur, for instance, if the breaking is from higher dimensional operators or arises from another scalar field acquiring a VEV.The kinetic energy of the field becomes comparable to the potential energy when θi (a i /a) 3 ≈ m a (T ), so large initial velocity and late generation both delay the onset of oscillations and increase the relic abundance.Kinetic misalignment tends to produce denser minihalos than conventional misalignment [72,73] due to a parametric resonance that enhances fragmentation [74].If the fragmentation is not complete the power spectrum of axion density perturbations has features at many scales and our power law ansatz will not be a good approximation.However, if the fragmentation 2. The exclusion plot of axion parameters from axinovae.The colored region represents the exclusion region assuming the axion star growth M ∝ t 0.2 , where the exclusion is independent of Tosc.The gray regions are the existing limits from the black hole superradiance.Dashed curves are the axion parameters giving the correct relic abundance assuming axion mass behaves as ma(T ) = ma(Λ/T ) b , where Λ = √ mafa is the dark confinement scale.A large b can naturally come from a first-order dark QCD phase transition.
completes before the kinetic motion is depleted the power spectrum is well approximated by white noise [72].In both cases the late-time halo mass function is peaked such that most of the mass is in mini-halos of mass M peak .While there have been many models that can enhance either the axion relic abundance or the self-coupling, diluting the relic abundance is also possible in scenarios such as nonstandard thermal histories that lead to entropy production [75].

Appendix E: Enhanced Axion Self-Coupling
The axion self-coupling is given by |λ| ∼ m 2 a /f 2 a ∼ Λ 4 /f 4 a , assuming a cosine instanton potential.To obtain the right relic abundance for axion dark matter, f a is usually large since the relic abundance of axions is proportional to f 2 a .However, axion self-couplings can be enhanced without affecting the standard misalignment mechanism or the formation of axion miniclusters.If the axion couples to two confining sectors, which can be naturally achieved with clockwork mechanism [66], the axion potential is Here Λ 1 , Λ 2 are the confinement scales of the two strongly coupled sectors and f 1 , f 2 are the corresponding decay constants.
We consider the situation where the vacuum misalignment mechanism is mostly set by V 1 (a) and so we require V 1 (a) V 2 (a) and V 1 (a) V 2 (a) which corresponds to the requirements Satisfying these constraints will guarantee that the misalignment mechanism and the axion mass term and the rolling of axion field are solely determined by the strong sector with a confinement scale of Λ 1 and breaking scale f 1 , which will be responsible for the relic abundance of the axion particles.However, this does not fully determine the axion self-couplings.If f 1 f 2 , the self-coupling can be dominated by the other strong sector, as long as the following condition is satisfied The conditions (E2) and (E3) can be consistent with each other provided f 1 f 2 .For instance, if Λ 2 /Λ 1 ≡ 1 then f 2 /f 1 ∼ ζ , with 1 < ζ < 2, will satisfy the conditions.Assuming the strong coupling sectors satisfy these requirements then m a ∼ Λ 4 1 /f 2 1 and |λ| ∼ Λ 4 2 /f 4 2 and the effective decay constant that labels the self-coupling strength is Therefore, the effective decay constant of an axion model that gives the self-coupling strength can be much smaller than the decay constant that is responsible for the relic abundance.They can be considered as two independent parameters.