Probing Relativistic Axions from Transient Astrophysical Sources

Emission of relativistic axions from transient sources, such as axion star explosions, can lead to observable signatures. We show that axion bursts from collapsing axion stars can be detectable over the wide range of axion masses $10^{-15} \, \textrm{eV} \lesssim m \lesssim 10^{-7} \, \textrm{eV}$ in laboratory experiments, such as ABRACADABRA, DMRadio and SHAFT. The detection of axion bursts could provide new insights into the fundamental axion potential, which is challenging to probe otherwise. An ensemble of bursts in the distant past would give rise to a diffuse axion background distinct from the usual cold axion DM. Coincidence with other signatures would provide a new window into multi-messenger astronomy.


I. INTRODUCTION
With the discovery of the Higgs boson, the prominent role that scalar particles play in nature has become apparent. The quantum chromodynamics (QCD) axion [1-7] is a leading solution to the strong CP problem and also constitutes an attractive dark matter (DM) candidate [8][9][10]. Axions, as well as axion-like particles (ALPs) in general, are expected to be ubiquitous in string theories (e.g., [11,12]). Significant efforts have been devoted and are ongoing to explore the axion/ALP parameter space (see [13,14] for a recent overview).
Axions associated with the Galactic DM halo contributing to cold DM are understood to be nonrelativistic and are the prime target of experimental efforts [13,14]. Relativistic axions can also be readily produced through a variety of mechanisms in the early universe, potentially contributing to a cosmic axion background (e.g., [15] and references therein). Experimental searches for continuouslyproduced relativistic axions from astrophysical environments like stars have been also widely explored [16].
Here we explore the possibility of detecting relativistic axions from transient astrophysical sources using terrestrial detectors. We discuss both the simple case of a transient axion burst signal in the limiting case of minimal wave spreading, as well as the well-motivated case of axion star bosenovae based on a concrete production mechanism and the detailed numerical simulations of Ref. [27].
We show how transient sources can lead to signatures that can be used to explore the fundamental axion potential, which is otherwise difficult to probe via conventional searches for cold axion DM. Further, we identify that the historic accumulation of axions from transient sources would give rise to a diffuse axion background. Our results are general, apply to many models and can be extended to other new physics searches from transient sources.
When the axion star mass M exceeds the critical value M c , the leading attractive self-interaction destabilizes the star, triggering collapse. The collapse is self-similar and well-described by a nonrelativistic evolution, until the final moments, when self-interactions induce relativistic decay processes in the star that deplete its mass [26] as the star expands again. The star collapses and expands O(few) more times before settling into a gravitationallybound and non-relativistic configuration [27]. This clump may relax again into a new dilute axion star.
The differential energy spectrum dE/dk of emitted axions from the collapse process has been studied via numerical simulations in [27] (see Fig. 3 in A). The leading relativistic axion momentum (k) peak was found to be at k/m ≈ 2.4, with less prominent peaks at higher values of k. The collapsing axion stars are found to lose ∼ 30 − 60% of their initial mass into outgoing axion radiation, which we will show to be potentially detectable. The duration of the axion burst is approximately δt burst ≈ 400/m ≈ 30 ns/m 5 1 .
Further details about axion stars and the collapse simulations of [27] are reviewed in A.

III. AXION BURST PROPERTIES
Consider an astrophysical burst of total energy output E, a distance R from the detector, which emits axions of fixed momentum k 0 = q m, where q 1 for relativistic axion emission. If the corresponding frequency ω 0 of the outgoing axions is in the sensitive range of a given axion DM experiment, then the burst may be detectable.
As the burst travels toward the detector, wave spreading will dilute the total energy density ρ as ∝ (δt) −1 , where δt is the apparent burst time as seen at the detector. In the absence of wave spreading, this energy dilution is dictated solely by the burst duration at the source, δt = δt burst , or in terms of the equivalent length scale, δx = δx burst . The burst energy density will also be diluted as ∝ R −2 due to the propagation of spherical waves away from the source. Altogether, in the limiting case when 1 Though not explicit in [27], based on dimensional analysis we argue that the burst duration δt burst must be approximately independent of f . The scaling with m is given, and aside from f , the only other relevant dimensionful parameter present is M P ; however, gravity is practically decoupled during the final stages of collapse, as long as f M P . Therefore, δt burst will be practically independent of M P and must therefore also be nearly independent of f .

R
δx, the energy density at the detector takes the form For a QCD axion bosenova of the kind simulated in Ref. [27], we can characterize the burst by the total energy emitted around the main relativistic momentum peak for a single explosion, which is In the limit of minimal wave spreading, this implies ρ * ∼ 10 7 ρ DM f 2 12 (100 AU/R) 2 , where ρ DM ≈ 0.4 GeV/cm 3 is the local DM density. See A for further details.
The temporal properties of the burst are also affected by wave spreading effects. In the absence of wave spreading, δt = δt burst , whereas when wave spreading dominates, δt ≈ δk/m × R/(q 2 q 2 + 1). Additionally, the apparent coherence time τ * of the axion burst (as seen at the detector) is τ * = 2π/δω ∼ 2π q 2 + 1/(qm) in the absence of wave spreading, but increases to τ * ≈ 2πR/(q 3 mδt burst ) in the limiting case of strong wave spreading. The latter follows from the consideration that, in the case of strong wave spreading, at a given time the detector is only immersed in relativistic axion waves with a small dispersion of δk ≈ mq 2 q 2 + 1δt burst /R. For further details concerning wave spreading, see B. Wave spreading can affect the resulting sensitivity, as we discuss below.
In general, if objects of mass M constitute a f DM fraction of the DM, and explode on an average timescale of τ , then the number of exploding objects within a distance R of a detector on Earth will be N star (R) = (f DM ρ DM /M )(4πR 3 /3), which for an axion star with mass M = M c gives N star (R) ≈ f DM (R/100 AU) 3 (m 5 /f 12 ). An experiment running for t int = 1 yr could thus detect at least one axion star explosion on average at a distance R if In what follows, we set τ = 10 Gyr (comparable to the present age of our Galaxy), and we ensure that our choice of f DM satisfies gravitational lensing bounds in the mass range 10 −11 M M 10 M [40]. See C for further discussion of axion star explosions, including the distribution of masses and frequency of explosions.

IV. AXION BURST SIGNALS
The transient signal from a relativistic axion burst differs from the usual DM signal associated with a cold oscillating axion field in several important aspects. DM axions are expected to oscillate coherently on the timescale τ DM ≈ 2π(mv 2 DM ) −1 ∼ 2π × 10 6 /m, where the typical virial velocity within the Milky Way DM halo is v DM ∼ 10 −3 , implying an oscillator quality factor of Q DM ∼ 10 6 . On the other hand, a relativistic axion burst is composed of incoherent waves and will only have a quality factor of Q * = O(1) when the effects of wave spreading are negligible. As a result, the signal-to-noise ratio (SNR) in our model will generally be suppressed by some power of the ratio Q * /Q DM 1. However, the burst energy density reaching a detector, ρ * , can greatly exceed the local DM density. In contrast, for a relativistic cosmic axion background from the early universe, like that considered in Ref. [15], the energy density is typically well below that of the cold DM component.
The sensitivity to temporally coherent oscillating signals improves with the integration time t int in an experiment as ∝ t 1/2 int , up to a maximum number of oscillations dictated by the coherence time.
In the temporally incoherent regime, the SNR still grows with the integration time, albeit more slowly, as ∝ t 1/4 int [41][42][43]. We focus on signals that are linear in both the axion field and the coupling constant g. Since the quality factor for the axion burst signal is expected to be Q * = O(1), there is no expected benefit from resonant-type experiments, such as those described in [44], and so we focus solely on broadband-type experiments.
Traditional searches for axion DM usually exploit the axion's possible coupling to the electromagnetic field: where F is the electromagnetic field tensor andF is its dual. The coupling in Eq. (5) is derivative in nature and gives rise to an axion-induced effective , where E and B denote applied electric and magnetic fields, respectively. The time derivative and spatial gradient associated with a spinless field φ have the typical sizes |∂ t φ| ∼ εφ 0 and |∇φ| ∼ kφ 0 , respectively, where ε = γm is the typical particle energy (γ = 1/ √ 1 − v 2 is the Lorentz factor, with v being the particle speed), k = γmv is the typical particle momentum, and φ 0 is the typical field amplitude. The energy density associated with the field is given by ρ ∼ ε 2 φ 2 0 . Broadband laboratory searches for axion DM via the electromagnetic coupling in Eq. (5) mainly search for an axion-induced oscillating magnetic flux in the presence of an applied static magnetic field [45][46][47][48][49][50][51]. In this case, the sensitivity to a relativistic axion burst (with the axion-photon coupling of Eq. (5) denoted g * ), relative to the standard cold DM search paradigm (denoted g DM ), is given by: where we have made use of the inequality τ * < δt, with τ * being the coherence time of the axion burst as seen at the detector. In deriving Eq. (6), we have assumed that the apparatus is capable of sampling data points at a rate of at least O(m), which allows one to optimally search for transient signals in the collected data.
Axions can also couple to fermions via derivativetype couplings. In this case, Eq. (6) is modified by the presence of an extra factor of v DM /v * ≈ v DM , which implies an enhanced sensitivity to a relativistic axion burst by the factor of v * /v DM ∼ 10 3 compared to Eq. (6). However, we are not aware of sufficiently sensitive broadband techniques for the relevant mass range in this case. We discuss the axionfermion coupling and some other types of couplings in more detail in D and E.
As a consequence of the relations above, if δt < t int , then the SNR remains the same regardless of the degree of wave spreading, since the signal can be captured in its entirety. On the other hand, if δt > t int , then the experimental sensitivity will degrade in accordance with Eq. (6). We also remark that the comparison in Eq. (6) involves a fixed value of ω 0 , rather than m. Since cold DM axions have ω 0 ≈ m, whereas axions in bursts are relativistic with ω 0 m, experiments searching for relativistic axion bursts can therefore be sensitive to lower axion masses compared to cold DM searches at the same signal frequency.
In Fig. 1, we estimate the regions of parameter space where the sensitivity ratio g * /g DM < 1 using Eq. (6) for the photon coupling, and including a v * /v DM ∼ 10 3 enhancement factor for the fermion coupling (see D). We focus on the optimal case of a minimal-uncertainty burst [δk ≈ (δx burst ) −1 ] of very short duration (δt burst = 2π/m), and we assume t int = 1 yr. Since this is a ratio of sensitivities, the result is independent of the details of a particular broadband-type experiment. It is intriguing that, for both types of couplings, the sensitivity ratio favors the detection of relativistic axion bursts over a wide swath of feasible parameter space (e.g., for E M and R pc). For the specific case of the QCD axion star collapse [27], we use E as given by Eq. (3) and depict the result by the light blue shaded region in Fig. 1.
[Note that in practice, for the axion star bursts we consider here, δt 1 yr for the main relativistic peak at q ≈ 2.4, as long as R O(10) pc.] When the ratio g * /g DM < 1 at some value of ω 0 , this implies the possibility of the detection of an axion star bosenova with enhanced sensitivity over analogous axion cold DM searches, provided that at least one axion bosenova occurs during the experimental measurements. In fact, we will see that QCD axion stars on the KSVZ line can potentially be probed by searching for relativistic axion bursts.
In Fig. 2(a), we illustrate the axion-photon coupling g φγ sensitivity reach of a future ABRACADABRA-type 2 instrument to our axion star bosenova signal (using Eq. (3) for the emitted energy as before), along with existing constraints, for the source distances R = 10 3 AU, 1 pc, 10 pc. 2 ABRACADABRA is in the process of merging with the DMRadio collaboration [50,51]; for simplicity, we use the long-term projection in [47] to estimate the sensitivity to our signal of interest. For the current ABRA-CADABRA limits, see [48,50]; see also the related limits from SHAFT [49].
The translation of the parameter g φγ to 1/f in Fig. 2(a) is model dependent; for concreteness, we have used |g φγ | = 1.92 α/(2πf ), which is consistent with the KSVZ axion model when the ratio of the electromagnetic and colour anomaly coefficients is given by E/N = 0. The sensitivity of conventional cold axion DM searches falls off as 1/f (or faster), because the axion-to-standard-model couplings are proportional to 1/f . On the other hand, the energy emitted in axion star bosenovae generally grows in size as ∝ f , since E ∝ M c ∝ f ; therefore, the product of the two factors is independent of f for axion star burst signals in ABRACADABRA-type detectors, and hence the sensitivity regions in Fig. 2(a) are vertical lines in the f − m plane. Axion star bosenovae may thus be a preferred method of discovery for large values of the decay constant f . In Fig. 2(b), we consider the sensitivity of a future ABRACADABRA-type instrument to explosions of arbitrary E, taken equal to the DM object mass M . Fixing the expected number of axion star explosions in 1 year, N , implies a relationship between R and M = E at fixed τ and f DM ; we vary E and show the most optimistic reach of an ABRACADABRA-type detector. In the estimation, we assume a minimumuncertainty and short-duration burst, as in Fig. 1. We observe that generic axion star bursts can be both detectable (SNR > 1) and occur frequently (N > 1) over a wide parameter space. We include for reference a region for N = 10 −4 , which could be relevant for a scenario in which τ is much longer than our assumption of 10 Gyr.

V. IMPLICATIONS
In models of cold oscillating axion DM with the cosine potential (1), the leading-order quadratic term of the potential gives rise to a cosinusoidal timevarying function of amplitude φ LO = φ 0 and angular frequency set by m. Higher-order terms in the cosine potential (1) modify the purely cosinusoidalin-time function into a Jacobi elliptical function due to the effects of anharmonicity. The next-to-leadingorder quartic term in Eq. (1) induces a correction to the leading-order cosinusoidal time-varying solution of the order of δφ NLO /φ LO ∼ φ 2 0 /f 2 ∼ ρ/(mf ) 2 ∼ 10 −37 /(m 5 f 12 ) 2 , assuming that the oscillating axion field saturates the local cold DM density ρ DM . Unlike traditional cold DM searches, the signal from axion star bursts can only arise in the presence of axion self-interactions of at least next-to-leading order (φ 4 ) in the axion potential, such as in Eq. (1). Furthermore, in the final stages of axion star collapse, the underlying processes are relativistic, and so one would expect next-to-next-to-leading-order (φ 6 ) and  [47]) to axion star bosenovae using Eq. (6) (blue regions), for the different explosion distances R = 10 3 AU, 1 pc, 10 pc. The QCD axion band is illustrated in cyan. The red lines mark N ≥ 1 for different choices of the three relevant parameters {R, τ, fDM} as labelled, chosen to satisfy gravitational lensing bounds [40], and the black dashed contours mark fixed values of the axion star mass M . The green regions are excluded by astrophysical observations [52] or by the non-observation of black hole superradiance (BHSR) in rapidly-rotating black holes [30]. For concreteness, we have use the model-dependent expression |g φγ | = 1.92 α/(2πf ) to translate coupling strength in the experimental limits to the vertical axis scale 1/f . The gray regions denote f > 10 17 GeV (light gray) from the consideration of axion star instability [53], and f > MP = 2.4 × 10 18 GeV (dark gray) from theoretical considerations [54] (though see also [55,56]). Finally, the yellow regions indicate where parametric resonance conversion of axions to photons would be relevant, either during axion star collapse (dark yellow) or for stable axion stars at M = Mc (light yellow) [57]. higher-order terms in the axion potential to have a non-negligible effect on the structure of the axion emission spectrum in Fig. 3. Therefore, the detection of such bursts could provide insight into the fundamental axion potential, which is challenging to probe via conventional cold DM searches.
In analogy with the diffuse neutrino background (e.g., [58][59][60]), relativistic axions from historic transients will accumulate into a diffuse axion background. The transient event rate would depend on the considered source formation and axion emission model. Highly redshifted axions could become nonrelativistic, while the boosted diffuse axion background component is expected to have a distinct phase space distribution.
Depending on the specific axion coupling, a variety of multi-messenger signatures accompanying relativistic axion bursts could be potentially expected from transient sources. This includes radiophoton emission via the axion-photon coupling (e.g., [61]), see the yellow regions in Fig. 2(a), as well as gravitational waves in the case of binary mergers or asymmetric axion star explosions. Such coincidence signatures would provide a complementary handle for exploring relativistic axions from transients. Searching for correlated signatures using a network of spatially-separated detectors would allow for the localization of the source. Some aspects of the detection prospects of relativistic axions and the multi-messenger aspect have been considered recently in the context of neutron star and black hole mergers [62], as well as neutron star -axion star collisions [63]. Further, it is expected that a typical axion star explodes several times before finally settling in a dilute, gravitationally-bound configuration [27]. Additionally, relativistic axions may convert into photons in Earth's ionosphere (see, e.g., [64]). We leave the detailed exploration of such signatures for future work.
Finally, we have pointed out that resonant-type experiments are less advantageous for axion burst signals than are broadband searches. However, this can be partly ameliorated by the use of simultaneous resonant searches in a narrow range of frequencies, which allows for the capture of a greater fraction of the burst energy and for probing a larger region of the emission spectrum. This constitutes a middle-ground approach between purely resonant and purely broadband search strategies. More generally, the broadness of the axion burst signal provides an important way to distinguish it from a cold DM signal (which would be very narrow in frequency space). A full treatment of the signal shape for axion bosenova signals, including the effect of multiple subsequent explosions as found in [27], is left for future work.

Appendix A: Basic Properties of Axion Stars
Axion stars are characterized by a classical field wavefunction in the non-relativistic limit. For a Klein-Gordon field φ, we define [65] φ(r, t) = 1 √ 2m e −i m t ψ(r, t) + h.c. , (A1) and the equation of motion for ψ(r, t) is a nonlinear Schrödinger equation at leading order in the self-interaction potential. The gravitational potential V g (|ψ| 2 ) is defined by and the total mass is defined by M ≡ m d 3 r |ψ| 2 . At low densities, the self-interaction term is negligible and gravity can balance the gradient energy, giving rise to a stable configuration known as a "dilute axion star" [36][37][38][39]. Since axion star wavefunctions are non-compact, we define the effective radius as R = R 99 , inside of which a mass 0.99M is contained. Solutions to the equations of motion in the dilute limit satisfy with the maximum mass (minimum radius) allowed by structural stability 3 of the star given by When axion stars exceed their critical mass M c , they collapse and explode in a bosenova of relativistic particles, a process simulated in Ref. [27]. We reproduce the emitted axion spectrum in Fig. 3. In [27], it was shown that after N * = O(few) explosions, the axion star settles into a diffuse, gravitationally-bound configuration; in the end, the total energy loss of the star is well-fit by the linear To get a rough estimate of the signal from this bosenova, suppose that a collapse and bosenova of the kind simulated in Ref. [27] occurs a distance R from a terrestrial detector. The integral over the main relativistic peak in Fig. 3, centered around k/m ≈ 2.4, gives which we use in the calculation of Eq. (3) in the Main Text, with N * = 1. Higher-momentum peaks are relevant to the signal as well, though we leave a full analysis of the signal shape for future work. Note that the largest proportion of emitted axions are found in the non-relativistic region of the spectrum, particularly at k/m 0.5. It would be interesting to consider such axions as a transient, non-relativistic, DM-like signal, although this is beyond the scope of the present work, which focuses on relativistic axion bursts. Note however that very non-relativistic axions emitted from a burst would, given enough time, become indistinguishable from the usual cold DM. The duration of the axion star burst for f = 2 × 10 −4 M P = 2.4 × 10 15 GeV in the simulation of [27] was approximately δt burst ≈ 400/m ≈ 30 ns/m 5 . Since the emitted axions are relativistic, the axion burst will have a corresponding intrinsic spread in space of δx burst ≈ δt burst ≈ 10 m/m 5 along the direction of propagation. In this approximation, substituting Eq. (3) into Eq. (2), we obtain: which does not depend on the size of the terrestrial detector. This will be modified by wave spreading effects, see B.

Appendix B: Wave Spreading Effects
Our signal depends on the duration of the burst, δt, as seen at the detector, which in turn depends on the spread in momentum of the emitted axion particles at the source. A true δ-function power spectrum in frequency space is impossible, as the power spectrum cannot be infinitesimally narrow without violating the uncertainty principle; the spread in momentum δk will be at least of order (δx burst ) −1 . Even a small spread in the emitted axion momenta can be important, as it leads to a further dilution of the axion energy density at the position of Earth due to wave spreading (in the particle picture, the fastest axions reach Earth sooner than slower ones, even if they are emitted at the same time). This is particularly important if an explosion occurs sufficiently far away. Here we derive the relevant time/distance scales needed for the analysis in the Main Text.
Suppose the axion momentum has a spread δk around some central value k 0 . The relativistic energy-momentum relation is ω = √ k 2 + m 2 = m q 2 + 1, with q = k/m, which implies v 2 = q 2 /(q 2 + 1). Taking the differentials of both sides, Eq. (B1) is relevant because the wave spreading effect is proportional to δv/v. In particular, for an instantaneous burst with small momentum spread δk/k 1 around the central peak at k 0 , we find This implies using δt = δx q 2 + 1/q; this expression is used in the results of the Main Text in Eq. (6). The smallest momentum spread at the source that is allowed by the uncertainty principle is δk ≈ (δx burst ) −1 , which implies So the intrinsic burst duration is negligible as long as R/δx burst q(q 2 + 1) m δx burst , which holds unless q is very large or R is very small.
Let us now consider the special case of axion star bosenovae. As we have just demonstrated, wave spreading in flight can often dominate the intrinsic signal duration. When the burst duration is set by the axion mass scale m, we can write δx burst = ξ/m, where ξ is a dimensionless constant in natural units, in which case Eq. (B4) becomes: For the specific case of a bosenova explosion of a QCD axion star, for which ξ ≈ 400, with δk/m ∼ 1 and q ≈ 2.4 for the main relativistic peak in Fig. 3, the wave spreading effect dominates over the intrinsic spread for any burst outside of a small radius R 0.03 pc when m 10 −15 eV.

Appendix C: Axion Star Collapse Frequency
A thorough treatment of the distribution of axion stars and collapse rate is beyond the scope of this work. Here, we have parametrised the relevant effects via three constants: (a) the fraction f DM of DM axions contained in axion stars at the present day; (b) the peak value f c = M/M c of the (unknown) mass distribution of stable axion stars in our Galaxy; and (c) the typical timescale τ for an axion star to grow in mass to the critical point and collapse. In lieu of a detailed study of axion cosmology, we comment on each of these points below: 1. Cosmological simulations of the axion field have recently been performed by several independent groups [25,[66][67][68]. For the QCD axion, it has been shown that if the axion global symmetry (e.g., the Peccei-Quinn symmetry [1]) is broken after inflation, overdensities known as axion miniclusters naturally form from field fluctuations shortly before matterradiation equality. Quantitative results, including the spectrum of overdensities, depend on very sensitive details of global string and domain wall networks in the early universe, a topic of heated debate of late [69][70][71][72][73]. Still, a qualitative picture is emerging in which these miniclusters form in the early universe and subsequently undergo mergers and possibly tidal destruction inside of galactic DM halos, and those that survive seed the formation of axion stars [24,25,74,75].
In direct simulations of axion minicluster halos, it has been shown that a large fraction of order ∼ 0.75 of the axion DM density remains in miniclusters down to redshift z ≈ 100 [76]. At lower redshift, uncertainties about merger histories and tidal stripping make it very difficult to determine the presentday distribution of these miniclusters (for recent efforts, see [77,78]). There are also so-called hybrid DM scenarios in which axions and primordial black holes both constitute a sizeable fraction of DM; in that case, axions can form axion stars in the gravitational potential of black holes directly from the galactic background [79]. On the basis of the above, it is not implausible to expect axion stars to constitute a few percent or more of the total DM abundance in galaxies.
To estimate the maximum possible reach, we would like to consider the case when f DM = 1. On the other hand, gravitational lensing experiments rule out axion stars with large values of f DM in some range of axion star masses. In the mass range 10 −11 M M 10 M , the constraint is as strong as f DM < 5 × 10 −3 , but can also be significantly weaker over that mass range [40]. In our rate estimates, we check to ensure that the abundance assumed is consistent with the gravitational lensing bounds at each value of M that we consider.
2. The initial mass distribution of axion stars depends on the mechanism of their formation. For example, axion stars forming inside of miniclusters have a mass at formation determined by the density and virial velocity of their host minicluster [25].
Furthermore, once axion stars form (either from relaxation of miniclusters or some other mechanism), they can merge [25,[80][81][82] or accrete further mass from the axion background [83]. Therefore, the mass distribution of axion stars will be ever-changing.
However, axion star masses do not grow without bound; as we have described, when the mass reaches M c , the axion star collapses and emits a large fraction E/M c ∼ 0.3 − 0.6 of its mass as relativistic particles [27]. Therefore, it is plausible that axion star masses will be clustered in some narrow range within an order of magnitude of M c . In this work, we make the simplifying assumption that all axion stars have f c ≡ M/M c = 1, in order to estimate the number of nearby axion stars and their mass.
3. In order for axion stars to be seen exploding in the sky near enough to Earth, it is necessary to ask how long it would take a typical axion star to grow enough mass (through accretion and mergers) to trigger collapse. The situation is actually more complex than mere mass growth: after an initial bosenova, simulations suggests the non-relativistic remnant of the original axion star remains gravitationally bound, though it retains only ∼ 40 − 70% of its original mass and the axions in the resulting diffuse configuration are not in their ground state [27]. We speculate, therefore, that the final axion configuration may relax again into a sub-critical axion star; if this is true, then this new star will begin to accrete mass again, growing towards M = M c , until it collapses and explodes again. These processes might occur repeatedly in galactic haloes, leading to a higher frequency of axion star bosenovae.
All of these complex dynamical processes deserve a dedicated treatment, but for the purposes of the present paper, we parametrise the overall result by a single timescale τ , the axion star lifetime, defined as the average total time for an average star to accrete mass, explode, reform into a gravitationally-bound relic, and relax again into a stable star. If τ is much greater than the lifetime of the Universe, then axion stars would be effectively stable and bosenovae will be exceedingly rare. On the other hand, if τ is too small, then DM haloes could become unstable; every explosion would convert an O(1) fraction of the axion star's mass-energy content (which may be associated with the cold DM) to relativistic axions that escape the galaxy 4 . Therefore, in the present paper, we consider the interesting intermediate value τ = 10 Gyr, which is comparable to the present age of our Galaxy.

FIG. 4.
Regions of parameter space for the axionnucleon coupling g φn where g * /gDM < 1 using Eq. (D2) (purple regions), for the different explosion distances R = 1 AU, 10 3 AU, 1 pc, 1 kpc. The red and black lines, as well as the cyan and horizontal gray shaded regions, are the same as in Fig. 2(a) in the Main Text. The green regions denote current astrophysical constraints on g φn [30,52].
It is crucial to emphasize that τ is highly modeldependent, and is worthy of a more complete treatment than we have offered here. For axion stars forming within axion miniclusters, for example, it has been suggested recently [25,83] that mass growth of an axion star by accretion saturates at some mass M sat < M c , constraining this possible mechanism for collapse at M M c (see also Section V of [84] for the case of ultralight axions). However, this too is somewhat model-dependent; for QCD axions with f 10 11 GeV (or ALPs with comparable self-interaction strength), the stars saturated in denser configurations, effectively driving M → M c prior to saturation and triggering collapse as we suggest. We leave further discussion of these issues for future work.

Appendix D: Derivative Coupling of Axion to Fermions
In the Main Text, we discussed effects of the axion-photon coupling that involve the time derivative of the axion field. Here, we discuss effects of the axion-fermion coupling that involve the spatial gradient of the axion field. Let us consider the derivative-type coupling of an axion field to fermions: where f is a fermion field andf = f † γ 0 is its Dirac adjoint. The spatial components in Eq. (D1) give rise to an interaction of the form H ∝ g φf ∇φ · Σ, where Σ is the Dirac spin matrix vector of the fermion. Broadband searches for axion DM via the fermion couplings in Eq. (D1) rely on the precession of polarised fermion spins about ∇φ [85][86][87][88][89][90]. In this case, the sensitivity to a relativistic axion burst, relative to the standard cold DM search paradigm, is given by: (D2) where we have used the inequality τ * < δt, analogously to Eq. (6) in the Main Text. In this case, the sensitivity to a relativistic axion burst is enhanced by the factor of v * /v DM ∼ 10 3 compared to Eq. (6).
In Fig. 4, we show the regions of parameter space where g * /g DM < 1 for the axion-nucleon coupling g φn , if a bosenova occurred within a distance R of the detector during an experimental search duration of t int = 1 yr. At present, the proposals for a broadband search for the axion-nucleon couplings in the relevant mass range are limited 5 . However, the sensitivity of such a search would be enhanced for the relativistic axion signal, relative to the cold DM case, by the factor of v * /v DM ∼ 10 3 , see Eq. (D2). We see again that relativistic axion bursts can outperform cold axion DM searches in the large-f region, even for quite distant bursts.

Appendix E: Non-derivative Axion Couplings
In the Main Text and in D, we have considered the primary couplings searched for in axion experiments, g φγ and g φn in Eqs. (5) and (D1), respectively, to estimate the sensitivity ratios in Eqs. (6) and (D2). For completeness, we also consider nonderivative couplings below.
The axion's coupling to the gluon field tensor G and its dualG, where b is the colour index, gives rise to time-varying electric dipole moments of nucleons [92] and atoms [86,93], which are non-derivative in nature. Additionally, scalar-type non-derivative couplings such as may arise in models of axions with CP violation in the quark sector [94] or in models of scalar-field DM that give rise to apparent variations of the fun-damental constants [95,96]. In the case of nonderivative couplings, the sensitivity to a relativistic axion burst, relative to the standard cold DM search paradigm, is given by: