Axion Star Explosions: A New Source for Axion Indirect Detection

If dark matter is composed of axions, then axion stars form in the cores of dark matter halos. These stars are unstable above a critical mass, decaying to radio photons that heat the intergalactic medium, offering a new channel for axion indirect detection. We recently provided the first accurate calculation of the axion decay rate due to axion star mergers. In this work we show how existing data concerning the CMB optical depth leads to strong constraints on the axion photon coupling in the mass range $10^{-14}\,{\rm eV}\lesssim m_a\lesssim 10^{-8}\,{\rm eV}$. Axion star decays lead to efficient reionization of the intergalactic medium during the dark ages. By comparing this non-standard reionization with Planck legacy measurements of the Thompson optical width, we show that couplings in the range $10^{-14}\,{\rm GeV}^{-1} \lesssim g_{a\gamma\gamma} \lesssim 10^{-10}\,{\rm GeV}^{-1}$ are excluded for our benchmark model of axion star abundance. Future measurements of the 21cm emission of neutral hydrogen at high redshift could improve this limit by an order of magnitude or more, providing complementary indirect constraints on axion dark matter in parameter space also targeted by direct detection haloscopes.


I. INTRODUCTION
Numerical simulations of axion-like dark matter cosmologies have shown that a solitonic core forms in the center of every dark matter halo, see [1,2] for the first simulations and [3][4][5][6][7][8] for further studies.In Ref. [9] we have explicitly calculated the number and energy density evolution of axion stars from hierarchical structure formation using semi-analytic models starting from an initial adiabatic perturbation spectrum.Ref. [2] demonstrated a power law relation between the mass of an axion star and its host halo: where z is the redshift, M h is the virial halo mass, α is a power law exponent, and [2]: where m a is the axion-like dark matter mass and M min is the smallest halo mass that could host a soliton of mass M S at a given redshift (given by the Jeans scale) [10,11].* miguel.escudero@cern.ch;-0000-0002-4487-8742 † charis.pooni@kcl.ac.uk; -0000-0002-6658-3478 ‡ malcolm.fairbairn@kcl.ac.uk; -0000-0002-0566-4127 § diego.blas@gmail.com;-0000-0003-2646-0112 ¶ xdu@carnegiescience.edu; -0000-0003-0728-2533 * * david.j.marsh@kcl.ac.uk; -0000-0002-4690-3016 Ref. [2] found a core-halo mass relation of the form M S ∝ M 1/3 h while Refs.[4,6,7] have found M S ∝ M 3/5 h .Refs. [8,12] suggest that there is an intrinsic diversity in the core-halo mass relation.Importantly, all numerical simulations of axion-like dark matter cosmologies show that the core-halo mass relation is bounded between α = 1/3 and α = 3/5.The Monte Carlo merger tree model we employed in Ref. [9] shows that an average slope α = 2/5 captures this diversity well for both the soliton mass function and, more importantly, the merger rate.In Figure 1 we show isocontours of constant M h for various values of M S and m a for the case α = 2/5.For the other cases we refer the reader to Appendix A. We note that while the simulations of Refs.[1][2][3][4][5][6][7][8] are performed for m a ∼ 10 −22 eV their results should apply to any axion mass for three reasons: 1) the axion stars appear to form during the gravitational timescale of the halo, 2) the primordial power spectrum expected from inflation is almost scale-invariant, and 3) the evolution of non-relativistic axion dark matter features a scaling symmetry that allows extrapolation to any axion mass, see e.g.[2].

II. AXION STAR EXPLOSIONS
Solitonic cores, also referred to as axion stars have ultrahigh phase-space occupation numbers that can trigger collective processes which cannot occur in vacuum.Parametric resonance can lead to exponentially fast decays of axion stars into photons [13][14][15][16][17][18][19] or into relativistic axions [20][21][22], meaning axion stars become unstable above a certain mass [13,19] (vector dark matter solitons have where g aγγ is the axion-photon coupling.Super-critical solitons explode into photons with E γ = m a /2.This happens by a collective process where photons produced by axion decays stimulate other axions to decay such that the axion star decay happens on a short timescale given by the light crossing time of the soliton: The energy released in the explosion of a supercritical soliton is E = M S − M decay S .There are, however, two factors that could prevent this explosion: 1) If the axion-like particle possesses self interactions, then axion stars can decay into axions leading a "Bosenova" above a critical mass M B−Nova set by the quartic self coupling λ [20][21][22][24][25][26].For a quartic coupling typical of a cosine instanton potential, λ = m 2 a /f 2 a , and axion-photon coupling g aγγ ≃ α EM /(2πf a ), one finds that M decay S ≃ 600M B−Nova .This means that for nominal values of g aγγ the axion stars will not actually decay into photons but rather into relativistic axions, with possible limits from such a decay explored in Ref. [27], see also [28].However, there are many models that feature significantly enhanced g aγγ interactions, see e.g.[29][30][31][32][33][34][35], or suppressed quartic interactions, see e.g.[36,37]. In hat follows, we will consider that the relevant effect is decay into photons and thus our results will apply to scenarios such as those previously mentioned.
2) In the parameter space of interest for our study, axion stars are always hosted by halos of mass M h < 10 5 M ⊙ which means that they lie below the baryonic Jeans scale and are too light to capture any significant amount of baryons from the intergalactic medium (IGM) [38,39].This means that on average the axion stars formed will simply be surrounded by a baryon environment with average IGM properties.Nevertheless, the IGM contains ionized particles, which could kinematically block the axion decay by plasma effects.More concretely, the photon gets an effective mass-squared that is proportional to the number of free electrons in the plasma [40]: Axion star decay is blocked until the plasma frequency drops to ω p (z decay ) < m a /2, resulting in an explosion of all the axion stars with M S > M decay S at specific redshift (see Appendix C): Once the plasma frequency has dropped such that parametric resonance decay is no longer blocked, then supercritical axion stars will explode as soon as they are formed.In Ref. [9] we calculated the cosmological evolution of the number/energy density of critical axion stars using the extended Press-Schechter formalism and Monte Carlo merger trees.We considered the case that an axion star will only explode either if it is plasma blocked until it is super critical, or if it is formed by a major merger defined by the merger of two axion stars of comparable mass, such that mass increase happens rapidly.
Soliton explosions lead to axion dark matter decay, and can inject energy into the IGM comparable to or greater than energy emission from astrophysical processes, such as core collapse supernovae [9], and thus offer a new method of indirect detection of axions.

III. HEATING THE INTERGALACTIC MEDIUM
Figure 2 outlines qualitatively the mechanism just described: axion stars start to form and grow when halos form and merge, which happens appreciably only at z ≲ 100.Once an axion star becomes massive enough and the plasma frequency is low enough, parametric resonance can take place and the star rapidly explodes into low energy photons of E γ = m a /2.These photons are absorbed by the IGM, heat it, and in turn via collisional ionizations lead to an increased ionization fraction of the Universe.
Axion star decay results in the release of a huge number of low energy photons with E γ = m a /2.We will be interested in m a < 10 −8 eV, whose associated low energy photons are efficiently absorbed via inverse Bremsstrahlung on ionized particles in the intergalactic medium, namely FIG. 2. Schematic of reionization caused by axion star explosions.Appreciable formation of dark matter halos hosting sufficiently massive axion stars occurs for z ≲ 100, see [9].When an axion star decays, it releases a huge number of low-energy photons, which are absorbed by inverse Bremsstrahlung, leading to heating of the IGM.If the IGM becomes hot enough, reionization occurs.Energetically, reionization requires a fraction of around f decay DM ≈ 10 −9 of the dark matter to decay.For the lowest mass axion-like particles, ma ≲ 5 × 10 −13 eV, axion star decay is kinematically blocked at early times, leading to a population of super critical stars which decay all at once in a burst once the plasma frequency falls low enough to allow the decay.This leads to patchy reionization.At higher axion masses, plasma blocking is not efficient at the relevant redshifts, and instead super-critical stars formed by major mergers decay immediately as they form, leading to a more uniform and continuous reionization history.
via γ e − p → e − p processes [41].The net result of this absorption is to heat the IGM.In turn, once the temperature of the IGM goes above T ∼ 1 eV, collisional ionization processes e − H → 2e − p ionize the Universe.As a result, the decay of axion stars may result into a period of early reionization, which is strongly constrained by Planck legacy data [42,43].This forms the basis of the constraints derived in this paper, see also [44][45][46][47][48][49][50][51] for constraints of this type for other dark matter candidates.
Importantly, depending upon the absorption length scale of the photons produced by axion stars the effect on the IGM can be different.Considering the γ e − p → e − p absorption length, see [41] and Appendix D, we find that for m a ≲ 5 × 10 −13 eV the photons are absorbed within very small volumes and this generates a shockwave very much like the one formed by supernova explosions as the amount of energy released is very similar [52].On the other hand, for m a ≳ 5 × 10 −13 eV the absorption length-scale is larger than the inter-separation of axion stars results in a homogeneous cosmological increase in the global IGM temperature.It is important to highlight that in principle a fraction as small as f DM ∼ 3 × 10 −9 ≃ 13.6 eV n b /ρ DM of dark matter converted into heat in the IGM would be enough to fully reionize the Universe, see Appendix B.

IV. METHODOLOGY
We model the IGM temperature by including the heating generated by axion star decays and accounting for cooling from Compton processes, collisional excitations and ionizations, and the expansion of the Universe.Collisional excitations and ionizations dominate in our scenario which we implement using the fitted rates in [53,54].We then track the free electron density x e using the effective 3-level atom approximation [55,56], including free protons and He + .We use the same rates and expressions as the recombination code RECFAST including their fudge factor [57,58] -modulo that we change T b → T γ in the terms in Eqs.(1) and (2) of [58] that account for photon ionization as the approximation T b ≃ T γ breaks down for z < 200.Finally, in the regime m a ≲ 5×10 −13 eV the photons from axion star explosions are absorbed on scales smaller then their typical hosting halo and generate shockwaves which ionize small patches of the Universe.We track such explosions until they are unable to ionize more IGM.We follow standard methods used for supernova remnants [52,[59][60][61], and find results very close to those expected from the Sedov-Taylor solutions [62,63].We refer the reader to Appendices E and F where we describe in detail all the relevant heating and ionization for the two regimes of interest, respectively.

V. CONSTRAINTS FROM THE CMB OPTICAL DEPTH
We derive conservative constraints on reionization driven by axion star explosions by using the 95% upper limit on the integrated Thomson optical depth to reionization from Planck CMB observations: τ reio < FIG. 3. Parameter space excluded by Planck measurements of the Thomson optical depth.We highlight also reach of future 21cm surveys and constraints from X/γ-Ray observations [64][65][66][67][68][69] and CAST [70], see [71].The evolution of the free electron fraction and the baryon temperature is shown for the points highlighted by colored symbols in Figure 4.
0.068 [42,43].In practice, and as customarily done for other types of energy injections [49][50][51], since observations show that the Universe should be fully ionized by redshift z ≃ 6 [73], and this implies a minimum contribution of τ reio = 0.0384, we then consider a region of parameter space excluded if the contribution to the optical width to reionization from axion star explosions at 50 > z > 6 is τ reio > 0.03.We note that this is a conservative approach as it assumes the Universe was not reionized until z = 6 by standard astrophysical sources.
The axion parameter space excluded as a result of this constraint is shown in the two regions in dark red in Figure 3, and the resulting evolution of the ionization fraction and IGM temperature for representative models is shown in Figure 4. Constraints in the region m a ≲ 5 × 10 −13 eV are a result of patchy reionization generated by the decay of all critical axion stars in the Universe once the axion decay is kinematically allowed, namely when ω p (z decay ) = m a /2.For m a ≳ 5 × 10 −13 eV the photon plasma mass is small compared to the axion mass, and the constraints arise instead from supercritical axion stars formed via mergers providing continuous heating and uniform reionization.The intermediate region of parameter space is constrained only at 1σ because in this regime the plasma frequency rises and blocks further axion star decay, see green curve of the upper panel of Figure 4.This region of parameter space is, however, expected to be tested with large scale CMB polarization data from LiteBIRD [74].
We note that the Planck constraints do not extend to arbitrarily large g aγγ -the larger the value of g aγγ , the smaller the mass of the halos that host critical soliton stars.In axion-like dark matter cosmologies structure formation is suppressed at M h < M min , see Eq. ( 2), and couplings g aγγ ≳ 10 −10 GeV correspond to axion stars For the scenarios we show there are no astrophysical sources of heating or reionization.We also show in black the ΛCDM evolution for xe from the Planck best fit cosmology and T b from the fiducial model of [72].In the lower panel we highlight scenarios that can be tested with future 21cm observations, conservatively (long dashed) and optimistically (dashed).
hosted in halos whose abundance is suppressed (which may in turn lead to further cosmological constraints, see e.g.[75,76]).In addition, we notice that the bounds change their shape for m a ≳ 5 × 10 −11 eV.For this region of parameter space the photon absorption efficiency in the IGM is smaller than one and this weakens the constraints.

VI. LYMAN-α AND SPECTRAL DISTORTIONS
We considered two other existing measurements that can limit this scenario for IGM heating.Since axion star explosions occur when Compton scattering is inefficient (z ≲ 10 4 ), they generate y-type distortions of the CMB energy spectrum.The change to the photon energy density from such an effect must satisfy δρ γ /ρ γ ≲ 5 × 10 −5 based on COBE/FIRAS measurements [47,77], which translates to f decay DM ≲ 2 × 10 −7 for decaying DM: significantly weaker than the Planck optical depth constraint.
At redshifts 2 ≲ z ≲ 7 the Lyman-α forest can measure the IGM temperature see e.g.[78][79][80].The IGM temperature at these redshifts satisfies T b ≲ 1.5 × 10 4 K which can also be used to constrain axion star decays.This is shown as the upper contour in Figure 3, with an accompanying T b (z) in Figure 4 (bottom panel).Lyman-α data could be competitive with the Planck optical depth constraint, although a dedicated study including other sources of heat and cooling mechanisms would be needed to set definitive constraints.

VII. 21 CM COSMOLOGY
Measurements of the hyperfine 21cm transition are expected to revolutionize our understanding of the cosmic dawn and the epoch of reionization, see [81,82] for reviews.There have already been several interesting bounds on the emission power of the 21cm line at various redshifts [83,84], and even a putative detection [85] (which is disputed see, e.g.[84]).The observational status of the cosmological 21cm line is still in its infancy but it is expected that ongoing and upcoming experiments such as HERA [86] and SKA [87] among others should be able to robustly detect the 21cm line and provide a view of the thermal state of the Universe at redshifts 4 ≲ z ≲ 30 [88], see [89] for a compilation of experiments targeting this line.
As highlighted in Figure 4, axion star explosions greatly enhance the IGM temperature at the redshifts where the 21cm line will be targeted z ≲ 35.In particular, the rise of T b at large redshifts z ≳ 20 is a feature that is not easy to achieve astrophysically, see e.g.[90].The sky-averaged 21 cm brightness temperature can be written as [91]: The sensitivity of SKA in this redshift range is expected to be ∆T 21 ∼ 10 mK, which means that SKA should be able to clearly differentiate between these two cases and thus future 21cm observations will test scenarios where axion star explosions heat the IGM.
In Figure 3 we highlight in light red the region of parameter space that could be tested by an experiment such as SKA by demanding that T 21 < 30 mK at z ≃ 20 (equivalent to T b ≳ 500 K by z ≃ 20).A tighter limit can be arrived at if one could differentiate any model from ΛCDM that causes T b > T γ at z ≳ 10.This is shown by the bottom light line in Figure 3.

VIII. CONCLUSIONS
Parametric-resonance instability of axion stars leads to an enhanced decay rate of axion dark matter into low energy radio photons.We have shown how the production of such photons heats the IGM, leads to reionization, and alters the optical depth of CMB photons.Planck legacy measurements of the optical depth form the basis of our new constraints on g aγγ highlighted in Figure 3, which is the strongest limit on axion dark matter in the relevant mass range by more than an order of magnitude.Future 21cm measurements of the IGM during the dark ages could improve this limit by more than an order of magnitude.The overlap of these indirect probes with the target regions of the DM-Radio haloscope program [93] invites further careful consideration of axion star explosions as a new tool in the hunt for dark matter.

IX. OUTLOOK
Aspects of this model that require further exploration include primarily a detailed study of soliton merger rates in simulations displaying a core-halo diversity in the redshift and particle mass ranges of interest.While our Monte Carlo merger trees in [9] suggest that M S ∝ M 2/5 h accounts well for diversity in the soliton merger rate only simulations in the redshifts of interests can fully support or modify this conclusion.Furthermore, the cosmological simulations of axion stars are typically carried out for m a ∼ 10 −22 eV.In these simulations axion stars appear to form upon halo collapse, and given the almost scale-invariant primordial power spectrum and the fact that the Schrodinger-Poisson equations solved in them features a scaling symmetry, their results can then be extrapolated to any axion mass.Nevertheless, it would be very interesting to perform explicit simulations of axion star formation and growth within the mass range of interest for our study, 10 −14 eV ≲ m a ≲ 10 −6 eV.
On the side of particle physics, the strongest assumption of our analysis is that either axion quartic couplings are suppressed or g aγγ couplings are enhanced compared with canonical expectations.While models of these types exist, construction of further explicit axion-like dark matter models above the QCD line would be required to address how natural such couplings are.Nevertheless, simulations of axion star explosions including both quartic and photon couplings may yet find that Bosenova triggered photon parametric resonance can occur for more standard coupling ratios.Furthermore, determining the true reach of 21cm measurements requires a full simulation of the 21cm anisotropies including anomalous lowenergy photon injection from axion stars.We suspect that anisotropies in this model will differ strongly from ΛCDM, and thus bounds may even be more powerful than our estimates.
In the main text we have shown all results for the benchmark core-halo mass relation in Eq. ( 1) with α = 2/5.Here we show the results for other two benchmark core-halo mass relations: α = 1/3 and α = 3/5.The former was the first one found in the literature [1,2].The case α = 3/5 was found in Refs.[4,6,7].In addition, Refs.[8,12] have shown evidence not for a strict corehalo mass relation but rather for diversity.Importantly, all numerical simulations show that axion stars lie somewhere in between α = 1/3 and α = 3/5.In addition, our Monte Carlo merger tree [9] shows that α = 2/5 captures well the effect of diversity in the merger rate of axion stars and in consequence represents our nominal choice for the core-halo mass relation.For completeness, in this appendix we show the resulting Planck constraints for these other two scenarios.In particular, Figure 5 shows the M S -M h relation for the two cases while Figure 6 shows the Planck CMB constraints for each of them.Planck constraints on the axion parameter space are strongly related to the amount of decaying dark matter that is injected as heat in the IGM.There are two distinct regions of parameter space.For m a ≲ 5 × 10 −13 eV there are three key effects: 1) axion stars cannot decay until the photon plasma mass is low enough, 2) the explosion thus occurs as a burst of energy when ω p (z) < m a /2, and 3) the photons generated from axion star explosions are absorbed on very small volumes and this generates shockwaves and bubbles of ionized IGM.This leads to patchy ionization.In this case, the constraints on the parameter space closely resemble regions of constant decaying dark matter density evaluated at the redshift of decay, z decay .This is explicitly shown in the left panel of Figure 7.There we can see that the region of parameter space excluded by Planck closely follows a region of parameter space where f decay DM ≃ 10 −6 .It is important to note that this is almost three orders of magnitude larger than what would energetically be needed if all the energy were to be injected homogeneously as heat in the IGM.However, since the observable τ is strongly sensitive to the volume of ionized IGM in the Universe this patchy reionization is much less efficient in fully ionizing the entire Universe.This is analogous to what happens in standard reionization where the most efficient sources of the global ionization of the Universe are UV and X-ray photons with mean free paths similar to the size of the observable Universe.
In the region of parameter space with m a ≳ 5 × 10 −13 eV the absorption of photons from axion star explosions takes place over distances that are larger than the typical inter-separation of axion stars.That means that the process can effectively be seen as homogeneous.In addition, in this case the emission is continuous because the photon plasma mass is small compared to m a /2.In this regime what matters is how much energy is released into heat and on which timescale does the IGM cool.The former we explicitly calculated in Ref. [9].The latter is simply given by Compton cooling and reads (see Eq. (E1c)): Thus, given energy arguments we expect regions of parameter space where df decay DM /dt/τ Compton ≳ 3 × 10 −9 can lead to reionization of the Universe.This is explicitly shown in the right panel of Figure 7 where we see that the Planck exclusion region is parallel and very close to this line.We notice that 21 cm observations will be sensitive to rates that are ∼ 2 orders of magnitude smaller than those currently tested by Planck measurements of τ reio .h .We can clearly see that for a fixed axion star mass the mass of a halo hosting it is much smaller for the case α = 3/5 than for the one with α = 1/3.This in turn means that significantly smaller gaγγ couplings can be probed.We highlight in grey regions of parameter space for which such axion stars cannot have possibly formed by z = 20 as their masses would lie below the effective Jeans mass generated by quantum pressure, see Eq. ( 2).

Appendix C: Cosmological Photon Plasma Mass
There are always some charged particles in the early Universe and this changes the global propagation properties of photons in cosmology.In particular, the photon gets a mass-squared that is proportional to the number of free electrons in the plasma [40]: where α EM is the electromagnetic fine structure constant, n e is the (free) electron number density and m e is the electron mass.Photons with energy E γ < ω p (z) cannot propagate.This in turn has important implications for the axion decay into photons because by pure kinematics the stars cannot explode into photons unless m a > 2ω p (z).In Figure 8 we show twice the plasma frequency as a function of redshift for the Planck best fit ΛCDM cosmology (as well as for considering a fully ionized universe).We see that axions with masses m a < 2×10 −14 eV cannot decay and thus our constraints are restricted to sufficiently massive axion-like particles.
The main cosmological implication of the plasma mass is that all the axion stars will not be able to explode until ω p (z) ≤ m a /2.For the range of redshifts z reio ≲ z ≲ 800 the relation between the redshift of decay and the axion DM /dt/τCompton.We can clearly appreciate how the constraints follow the same shape and follow closely the case 3×10 −9 as expected from energy arguments.The region ma ≳ 10 −10 eV deviates from this expectations because the probability of photon absorption in the IGM is smaller than 1.In addition, in dashed red we show regions of parameter space with τreio = 0.13, namely, with an optical depth that exceeds the Planck bound by a factor of 2. mass is approximately given by: This means that for axion masses in the range 2×10 −14 ≲ m a ≲ 10 −12 eV all axion stars that are massive enough will explode and convert all their energy into very low energy photons that can subsequently heat the intergalactic medium.For m a ≳ 10 −12 eV we will instead expect that axion stars to explode as soon as they form.In this regime, we expect a continuous emission of energy across redshift.
Effect of Over/Under-densities: We note that the plasma frequency displayed in Figure 8 is the one corresponding to the average electron density in the Universe.However, one could imagine that axion stars are at some point surrounded by either overdense or underdense baryonic environments.First, axion stars do form in the centers of dark matter halos and thus one could at first sight expect them to be surrounded by an overdense baryonic medium [94].However, the halos that host the axion stars that we consider in our study have masses M h < 10 5 M ⊙ .This in turn means that their potential wells are not deep enough to capture any significant amount of baryonic gas [38,39] and would thus we expect them to be surrounded by the average baryon density.It is also possible that some axion stars form in regions where the baryon density is below the cosmic average.This would in turn mean that one could in principle extend the constraints we discuss to slightly lower axion masses and this has been discussed in the context of light dark photon dark matter in [45,46,95,96].However, it is important to notice that the explosion mechanism of axion stars is independent on the value of the plasma frequency surrounding the star.The only effect of the plasma frequency would be to potentially delay the explosion itself.Furthermore, noting that the plasma frequency is only mildly sensitive to the baryon density, ω p ∝ √ n e , this means that we do not expect a significant effect from considering the small effect of over or under densities.
Appendix D: Absorption of photons from axion star explosions in the IGM Axion star explosions generate a huge number of very low energy photons with E γ = m a /2.For m a ≲ 10 −8 eV the most efficient absorption process of these photons in the IGM is free-free absorption also known as inverse Bremsstrahlung, γ e − p → e − p.The rate at which these photons are absorbed by the plasma is [41]: where σ T is the Thomson cross section, T b is the temperature of the electron-baryon fluid and where .This means that as the temperature of the baryons raises the free free absorption becomes less effective.Importantly, it also scales as x 2 e which means that as the Universe becomes ionized the absorption becomes highly efficient.
This rate should be compared with the Hubble expansion H(z) to see if these photons will be absorbed, which is well described around the redshifts of interest by Figure 9 explicitly shows the absorption lengthscale of photons produced from axion star decays for two characteristic baryon temperatures, T b = 9 K, and T b = 1 eV, which correspond to a cool IGM and one where the IGM is would be fully ionized.From this figure we notice three distinctive regions of parameter space: 1.For m a ≳ 10 −8 eV the photons have mean free paths larger than the size of the observable Universe and thus axion star explosions will not lead to relevant cosmological implications.
2. For 5 × 10 −13 eV ≲ m a ≲ 10 −8 eV the photons produced from axion star explosions will be absorbed by the IGM on length scales larger than the one it will take for this region to reach T b ∼ 1 eV, which is: FIG. 9. Typical length scale over which photons from axion star decays are absorbed in the IGM via inverse Bremsstrahlung on free electron-proton pairs with a plasma at two different temperatures, T b = 9 K, 1 eV.We show it for the reference value of z = 20 and with a pre-reionization free electron fraction xe = 2 × 10 −4 .We can appreciate three regimes: at ma ≳ 10 −8 eV the photons have mean free paths larger than the size of the Universe and therefore there are no cosmological signatures, for ma ≲ 5 × 10 −11 eV the lengthscale of absorption is always smaller than 1/H which mean that all of the photons will be absorbed even if T b ∼ 1 eV which is approximately the temperature of an ionized IGM.Finally, for ma ≲ 5 × 10 −13 eV the photons from axion stars are absorbed on very small lengthscales and thus generates an intense shock-wave which in turn generates a sort of patchy reionization.
That means that the absorption of photons in this region of parameter space can be seen as homogeneous, so it will take several nearby axion star explosions to heat up the IGM.This scale is labelled as typical bubble size in Figure 9.
3. For m a ≲ 5 × 10 −13 eV the photons are absorbed on very small lengthscales.Since a large amount of energy is released per axion star decay E ∼ M S ∼ 10 −4 M ⊙ the injection of these photons will lead to shockwaves in the IGM very similar to supernova remnants.The lower limit on m a of 2 × 10 −14 eV arises due to the photon plasma mass in the Universe.
It is important to note that these considerations have been made with x e ∼ 2 × 10 −4 as expected prereionization in ΛCDM.As soon as x e grows photons from axion star decays will be absorbed faster in the IGM, see Eq. (D3).
In what follows, in Appendix E we outline our modeling of energy injections for m a ≳ 5 × 10 −13 eV and in Appendix F we consider the case of m a ≲ 5 × 10 −13 eV.For m a ≳ 5×10 −13 eV axion stars explode into photons that are absorbed across lengthscales which are larger than the typical interseparation of axion stars.In this regime, the injection of energy can be seen as homogeneous.As discussed in the main text, we track the evolution of the baryon temperature including the new heating source from axion star explosions and also account for Compton and adiabatic cooling.We also consider for the net baryon cooling generated by collisional ionizations e H → 2e p, as well as collisional excitations e H → eH * → eHγ.We track the free electron fraction using the effective 3-level atom and follow the evolution of protons (p) and HeII ≡ He + .x i ≡ n i /n H so that x e ≡ n e /n H = x p + x HeII .The evolution equations that describe them are: which represent the rate equations for the free electron fraction contribution due to Hydrogen, Eq. (E1a) Helium, Eq. (E1b), recombination (α), photoionization (β) and collisional processes (coll).Both recombination collisional processes depend on the baryon temperature, T b , which is solved via the rate equation, Eq. (E1c), where the first term is simply adiabatic cooling, the second corresponds to Compton cooling, the third to the net gas cooling generated by possible collisions in the plasma, and the last term is a result of heating by axion star explosions.The energy deposited as heat in the IGM by axion star explosions is explicitly given by: where here the Θ function ensures that no emission is generated if the photon plasma mass blocks the decay, the last factor is an efficiency factor that takes into account that photons may not be absorbed if their mean free path is large, and finally the fraction of dark matter decaying into photons from axion star mergers is obtained from [9] and approximately reads as: The relevant energy transitions in the effective-three level model correspond to the ground state (1s), the second level with two quantum states (2s and 2p), and the third level (c), which denotes the continuum.Direct recombination and photoionization are prohibited, and the Lyman-α decay channel is heavily dampened by the optically thick plasma during the recombination epoch.Therefore, the only possible route to recombination (the ground state, 1s) is to go via the 2s− →1s γγ-photon decay channel.This occurs at a much slower rate compared to Lyman-α, Λ H,2s1s = 8.22458 s −1 which has a transitional energy, E H,2s1s = E H,2p1s = 10.2 eV, for Hydrogen, and Λ He,2s1s = 51.3s −1 with an energy transition of E He,2s1s = 20.62 eV, for Helium.This creates an overall effective decay rate for energy transitions between 2− →1, where γγ-decay channel dominates at earlier times.This effective decay rate is encapsulated by the Peebles Cfactor, which scales the rate equations in Eq. (E1a) and (E1b) due to these physical effects.The C-factor for Hydrogen, C H , and Helium, C He , are given by, C He = 1 + K HeI Λ He,2s1s n H (f He − x HeII )e −E2s,2p/Tγ (1 + K HeI (Λ He,2s1s + β HeI )n H (f He − x HeII )e −E2s,2p /T γ ) , (E5) respectively, where K H = λ 3 H,Lyα /8πH(z), λ 3 H,Lyα = 121.56nm, and for Helium, E 2s,2p = 0.60 eV and K HeI = λ 3 HeI, Lyα /8πH(z), where λ HeI, Lyα = 58.43nm [58].The recombination rates for Hydrogen, α H (T b ) and Helium, α HeI (T b ) in Eq. (E1a) and (E1b), are taken from, [97,98] and [99], respectively.
Assuming local thermal equilibrium and applying the fact the absorption (α) and emission (β) processes are in detailed balance, the photoionization rates are given by β H (T γ ) = (2πm e T γ ) 1.5 α H (T γ )e where f He = n He /n H , and where ⟨σv⟩ H and Helium, ⟨σv⟩ HeII , are the collisional ionization rates for Hydrogen and Helium taken from [54].The net cooling generated by collisional ionizations and excitations is given by where E H = 13.6 eV and E HeII = 24.6 eV.The first two terms arise due to cooling generated by ionization while the last term corresponds to cooling generated by collisional excitations.Namely, by processes of the type eH → eH * → eHγ, this is, collisions that are not able to ionize the neutral atom can nevertheless excite it and it will subsequently decay back to the ground state by emitting a photon.Since these photons are not absorbed by the gas the whole process cools it.We take these rates from [53] and which read: x e (f He − x He−II ) .
To cross check our code we considered also the decaying dark matter scenario studied in Ref. [47].There the authors consider a scenario with dark matter decaying into very low energy photons as we do.In their case the decay is modelled by a simple exponential law and in Figure 10 we show the excellent agreement for the evolution of the free electron fraction when we use their model of dark matter.In this reference, the authors did not include the cooling from collisional excitations and it is in that case where we find good agreement.In our study of axion star explosions we do include it since it has a relevant impact (as shown also in Figure 10).
Thomson CMB Optical Depth: Finally, after solving Equations (E1) we then calculate several integrated FIG. 10.Here we compare the free electron fraction evolution in a scenario of very light decaying dark matter that decays homogeneously and following a typical exponential decay law [47].For the comparison we have used the very same values as in Figure 15 of [47].We can appreciate the excellent agreement when collisional excitations are not considered, as is the case of [47].When collisional excitation cooling is considered, the effect on xe is somewhat reduced.
quantities.In particular, the Thomson optical depth to recombination which reads: In practice we choose z max = 50 and since we know the Universe was reionized by z = 6 we use x e = 1 for z ≤ 6 irrespective of the reionization evolution that axion star explosions lead to.Solutions for this reionization evolution are shown in Figure 4. Spectral Distortions: Additional heating of the IGM due to axion star explosions can generate a y-type distortion of the CMB at the redshifts of interest.The maximum amount of energy that could go into the CMB from axion star explosions is given by, where dQ/dz is the heating rate and ρ γ is the energy density of CMB photons.In Eq. (E11), we have chosen to integrate from z min = 0 to z max = 400, given z = 400 is a high enough redshift before axion stars have started to decay.With this expression, we have found that for the parameter space that is currently excluded by Planck constraints on τ the y distortion is much smaller than the one currently tested by COBE/FIRAS.In particular, for the red circle benchmark point in Figure 3 we find |y| = 3 × 10 −7 .This is two orders of magnitude smaller than the current sensitivity.For the diamond point in purple that can be tested by 21cm observations, we find |y| = 1 × 10 −8 , which could be within the sensitivity of future CMB observations.
Appendix F: Patchy Reionization: Growth of Bubbles and Initial Ionization of the Plasma For m a ≲ 5 × 10 −13 eV photons generated by axion star explosions in the early Universe are absorbed on very small lengthscales (see Figure 9).Axion star explosions release a huge amount of energy M S ∼ 10 −4 M ⊙ (see Figure 1) and will generate a shockwave very similar to those generated by supernova explosions that would be very violent and fully reionize the IGM around it.In this appendix we describe how we track the evolution of the bubbles generated by these explosions and calculate the volume average ionization fraction of the Universe.This scenario leads to patchy ionization.
The expansion of post explosionary shock-waves in a dense medium was first studied in the context of nuclear bombs [62,63].It was later realised that this could also be applied to astrophysical environments such as around supernova explosions [100] and the equations were subsequently developed to include the internal and external pressure changes, the heating of the interior of the bubble [101], energy lost to the shockwave due to the emission of radiation [102], the self gravity of the shell and the expansion of the Universe [52,59,103].
In order to model the behaviour of the expanding bubble with radius R and enclosed mass M , we use the following two coupled equations (see, e.g., [59][60][61]): where a dot represents derivative w.r.t time, Ω m is the total density parameter, Ω b the baryon density parameter, and H is the Hubble expansion rate.L tot represents the total luminosity and p is the bubble pressure resulting from this luminosity.The first term in equation (F1a) represents the driving pressure of the outflow, the second is drag due to accelerating the IGM from velocity HR to velocity Ṙ, the third is the self gravity of the expanding shell the fourth is the self gravity of the entire halo.The first term in equation (F1b) represents the increase in pressure caused by injection of energy while the second term is the drop in pressure caused by adiabatic expansion.In particular, the total luminosity is given by: where here L Explosion is the luminosity that is generated by the explosion of the axion star into photons, L Compton is the luminosity lost via Compton cooling against the CMB, and L Ionization takes into account the energy lost in ionising the swept IGM by the shockwave.Since the timescale for the explosion is very short compared to any other time scale, see Eq. ( 3), and the one for absorption as well, see Eq. (D3), in practice we assume that the blastwave expands freely until the mass contained in the bubble is comparable to the energy ejected by the axion star explosion.At that point we start the integration of Eqs.(F1a) and (F1b) neglecting the explosion luminosity but taking into account the fact that a large pressure has been generated by the explosion.This is similar to what is done for supernova explosions, see Section 8.6.1 of [104], and the starting radius, the velocity and the pressure read as follows: R 0 = 0.7 pc We thus assume that the pressure density is at that point essentially given by the energy output in the axion star explosion in the volume filled by the bubble and that the bubble is moving relativistically at a speed of v = 0.5 c.
From the initial conditions in Eq. (F3) we can then solve Eqs.(F1) with L Explosion = 0.The last thing needed is the other luminosities which can be written as [59]: where T γ is the CMB temperature at a given redshift, and where here n b is the background baryon number density, E H ≃ 13.6 eV is the energy it takes to ionize a hydrogen atom, and f m is the fraction of the baryonic mass kept in the interior of the bubble which by construction is f m ≪ 1.In our calculations, for concreteness we take f m = 0.1 but we have checked that smaller values yield very similar results.Finally, we stop the integration when the pressure of the interior of the bubble is p ≃ 2 T crit n b with T crit ≃ 15000 K which roughly corresponds to the temperature at which the IGM becomes is fully ionized.For smaller pressures the bubble will not be able to ionize the swept IGM and the radius at this point will tell us the region that has been ionized by the axion star explosion.Figure 11 shows the numerical solution to Eqs. (F1) for one energy injection example.For concreteness, we show the result for z decay = 20 and for two different axion star masses.There we can clearly see the steep decrease of the pressure as a result of the increase in the bubble radius.We note that we do not find any significant dependence of the results with respect to the redshift of injection.
By numerically solving for all the axion-star masses of interest we find that the comoving size that the bubbles

FIG. 1 .
FIG. 1.In blue dashed: isocontours of halo masses that host a given axion star mass at redshift z = 20 depending upon the value of the axion mass for the benchmark core halo mass relation MS ∝ M 2/5 h .In red we show isocontour lines of critical axion stars (M decay S ) as a function of several gaγγ values.

FIG. 4 .
FIG.4.Evolution of xe (upper) and T b (lower) for the six scenarios highlighted in Figure3with symbols.For the scenarios we show there are no astrophysical sources of heating or reionization.We also show in black the ΛCDM evolution for xe from the Planck best fit cosmology and T b from the fiducial model of[72].In the lower panel we highlight scenarios that can be tested with future 21cm observations, conservatively (long dashed) and optimistically (dashed).
where x HI is the fraction of neutral hydrogen in the Universe and T S is the spin temperature.At z ≲ 25 the spin temperature is T S ≃ T b[92] and we then see that for T b > T γ the 21cm signal could reach a maximum in emission of ∼ 35 mK at z ≃ 20 provided T b ≫ T γ .This is in strong contrast with what is expected in most cold IGM models where T b (z = 20) ∼ 7 K and which would lead to a strong absorption feature with T 21 ∼ −200 mK.
Appendix B: Decaying dark matter fraction and axion star merger evolution

FIG. 5 .
FIG. 5. Isocontours of halo masses (dashed blue) that host an axion star mass (y axis) as a function of axion mass (x axis) at a redshift z = 20, see Eq. (1).In red we show the value of a critical axion star above which the star can explode into photons given an axion-photon coupling in red, see Eq. (3).The left panel corresponds to the Schive relation, with MS ∝ M 1/3 h while the right panel corresponds to the case MS ∝ M 3/5

FIG. 6 .
FIG. 6. Equivalent to Figure 3 but for two other core-halo mass relations MS ∝ M 1/3 h (left) and MS ∝ M 3/5 h (right).In light red dashed lines we show the benchmark case with MS ∝ M 2/5 h for comparison.

FIG. 7 .
FIG. 7. Left: isocontours of f decay DM compared with the constraints obtained from Planck.Right: isocontours of df decayDM /dt/τCompton.We can clearly appreciate how the constraints follow the same shape and follow closely the case 3×10 −9 as expected from energy arguments.The region ma ≳ 10 −10 eV deviates from this expectations because the probability of photon absorption in the IGM is smaller than 1.In addition, in dashed red we show regions of parameter space with τreio = 0.13, namely, with an optical depth that exceeds the Planck bound by a factor of 2.

FIG. 8 .
FIG.8.Minimum mass of an axion that can decay into photons as a function of redshift, ma > 2 ωp(z).In red we show the evolution according to the Planck ΛCDM best fit cosmology and in purple we show the evolution in a fully ionized universe.

/2 T b m e − 7 / 2 . 3 π 2 × 10 −4 2 × 10
(D2) g BR is the Gaunt factor that for E γ ≪ T b can be approximated by g BR = √ log(2.25Tb /E γ ).Plugging in numerical values we find: Γ abs ≃ 2.6 × 10 −22 eV x e where here we have normalized the rate to the thermodynamic values of T b and x e as expected pre-reionization in ΛCDM at z = 20.We note that the rate has a strong redshift dependence and that it scales as m −2 a and T −3/2 b

10 −14 10 − 6 m
Photons are not absorbed in the Universe Photon absorption can be treated homogeneously Appendix E: Continuous Axion Star Mergers andExplosions: Homogeneous IGM Heating and Reionization