The Cosmic Axion Background

Existing searches for cosmic axions relics have relied heavily on the axion being non-relativistic and making up dark matter. However, light axions can be copiously produced in the early Universe and remain relativistic today, thereby constituting a Cosmic $\textit{axion}$ Background (C$a$B). As prototypical examples of axion sources, we consider thermal production, dark-matter decay, parametric resonance, and topological defect decay. Each of these has a characteristic frequency spectrum that can be searched for in axion direct detection experiments. We focus on the axion-photon coupling and study the sensitivity of current and future versions of ADMX, HAYSTAC, DMRadio, and ABRACADABRA to a C$a$B, finding that the data collected in search of dark matter can be repurposed to detect axion energy densities well below limits set by measurements of the energy budget of the Universe. In this way, direct detection of relativistic relics offers a powerful new opportunity to learn about the early Universe and, potentially, discover the axion.


I. INTRODUCTION
The existence of an axion with mass well below the electroweak scale could resolve the strong CP puzzle [1][2][3][4], and is entirely in line with UV expectations given the ubiquity of axions in string theory, where they arise from the deconstruction of extra-dimensional gauge fields [5][6][7]. The discovery of cosmologies where such a particle produced in the early Universe could constitute dark matter [8][9][10] has motivated a broad program to detect non-relativistic axions, and the development of instruments that will cover enormous swaths of unexplored parameter space in the coming decades. Yet the axion need not be dark matter, and the mere existence of an axion in the spectrum implies the possibility that a relic population of these states was produced in the early history of the Universe. Generically, such a population could be relativistic -a characteristic feature of the axion is an approximate shift symmetry, leading to a potential suppressed by powers of the axion decay constant, f a , and correspondingly the axion mass, m a , is expected to be small. Accordingly, the Universe may be awash in a sea of relativistic axions, a Cosmic axion Background (CaB).
In this work we will broadly discuss the production and detection of such a CaB. The possibility of a relativistic axion population is not new, and has been discussed in several contexts, including axion contributions to ∆N eff [11], axions with keV energies motivated by the prospect of moduli decaying into axions through Plancksuppressed higher dimensional operators [12][13][14][15], or constraints on the conversion of relativistic axions off primordial magnetic fields in the early Universe [16,17]. However, our focus here is to systematize the study of the * jdror1@ucsc.edu † hitoshi@berkeley.edu; hitoshi.murayama@ipmu.jp ‡ nrodd@berkeley.edu CaB and demonstrate the terrestrial detection prospects, thereby opening new paths to discovery. In addition to outlining a number of distinct scenarios where relativistic axions could be produced in the early, or even late, Universe, we will demonstrate that such a population can leave a detectable fingerprint in instruments designed to search for dark matter. Studies of an additional relativistic component to the Universe are particularly relevant in light of recent discrepancies in measurements of the Hubble constant between the early (H 0 = 67.4 ± 0.5 km/s/Mpc) and late (H 0 = 73.3 ± 0.8 km/s/Mpc) Universe [18]. An additional contribution to N eff of around 0.4 -which relativistic axions could providemay play a role in resolving the discrepancy, as they can relax the uncertainties in the early Universe measurement, giving a value of 66.3 ± 1.4, which would reduce, although far from resolve, the tension [19]. This provides an experimental target which we will denote by "H 0 Preferred" throughout. A simplified representation of the CaB landscape discussed in this work, is provided in Fig. 1. The black dashed curves show the differential axion energy density, Ω a (ω) (a precise definition is provided below), as a function of the energy, ω, for the CaB variants discussed in this work. The colored and shaded regions show the reach of two existing (solid curves) and future (dotted curves) instruments in this same space. We will explain this figure in more detail later in the introduction, but already we emphasize that dark-matter searches will probe interesting CaB parameters, particularly at lower frequencies.
In Fig. 1, and throughout this work, we will focus on the axion-photon coupling, In general, the coupling of the axion to the Standard Model (SM) is highly uncertain and there exist experiments targeting a number of different axion-SM couplings  (2), as a function of energy. The black dashed curves show four different realizations of the CaB, corresponding to thermal production (with Ta = T0, the CMB temperature), a Gaussian distribution representative of parametric-resonance production (with ρa = ργ,ω = 0.3 µeV, and σ/ω = 0.1), dark-matter decay (χ → aa), and cosmic-string production (fa = 10 15 GeV, T d = 10 12 GeV). For the dark-matter decay distribution the parameters are set to parameters already accessible to ADMX, as justified later in this work. In particular, we take mDM 5.4 µeV and τ 2 × 10 3 tU , with tU the age of the universe. While the thermal distribution will always peak roughly where shown and the cosmic-string production is dominant at lower frequencies, the parametric resonance and dark-matter decay signals can populate the full energy energy range shown. In all cases we set the axion photon coupling to the largest allowed value consistent with star-emission bounds over this energy range, g SE aγγ = 0.66 × 10 −10 GeV −1 . The colored regions denote the sensitivity in this same space that could be obtained by reanalyzing existing ADMX and HAYSTAC data, or with the future sensitivities of DMRadio and MADMAX. In determining the sensitivities, we have assumed that the CaB axion-photon coupling saturates star-emission bounds. We show the region of parameter space where the CaB could partially alleviate the Hubble tension, labelled H0 Preferred. Finally, the gray dotted line depicts the approximate boundary, to the left of which the CaB has sufficient number densities to be treated as a classical wave.
(for a review, see e.g. [20,21]). While we restrict our discussion to g aγγ , many aspects of the CaB extend to more general couplings.
At present, there are two primary classes of searches for axion backgrounds using the coupling in (1). The two strategies are broadly distinguished by where the axions are produced: a relativistic population produced in the cores of compact astrophysical objects or a nonrelativistic dark-matter population. For existing relativistic searches, the axions are produced in compact objects, such as stars like the Sun, which act as a bright source of axions with energies in the ∼keV range. Avoiding excess cooling of these objects from axion emission already puts a stringent bound on g aγγ [22] with comparable limits obtained by directly searching for the emitted axions in helioscopes [23] or absorption in direct detection experiments [24]. Together these searches, which we collectively refer to as "star-emission" bounds, are able to set strong bounds on axions with m a 1 keV, with the strongest limits across this full range given by g aγγ 0.66 × 10 −10 GeV −1 as determined by the CAST helioscope [25] and observations of Horizontal Branch stars [26,27]. For m a 10 −10 eV, these bounds can be strengthened by X-ray searches from conversion of axions emitted by SN-1987A [28] (assuming the supernova remnant is a proto-neutron star [29]), NGC 1275 [30], and super star clusters [31], reaching g aγγ 3.6 × 10 −12 GeV −1 . We collectively denote these existing star-emission bounds by g SE aγγ . These will be relevant as in the current work we will only consider axions with m a 1 keV, and therefore the same axions constituting the cosmic background could also be produced in these astrophysical objects, and must therefore satisfy g aγγ ≤ g SE aγγ . Dark-matter searches instead look for non-relativistic axions with a much larger local number density [23]. Tra-ditionally, axion dark matter has been searched for in microwave cavity haloscopes [32,33]. In the presence of a large magnetic field, axions in the 1 − 50 µeV mass range can resonantly excite the modes of an O(m) sized cavity (as m −1 a ∼ µeV −1 ∼ m). This detection principle underlies many of the strongest existing bounds on axion dark matter, as determined by the ADMX [34][35][36] and HAYSTAC [37] collaborations (see also Ref. [38]), which already require dark-matter axions in this mass range to have g aγγ orders of magnitude below g SE aγγ . Ideas are currently being developed to extend the accessible axion dark-matter mass window to both higher and lower values. For m a ≤ 1 µeV, resonant conversion can still be obtained when the axion power is read out through a high quality-factor lumped-element resonator [39][40][41][42]. A broadband readout of the signal in this mass range has already been used to set limits comparable to g SE aγγ by the ABRACADABRA [43,44] and SHAFT [45] instruments, and in the future DMRadio will aim to significantly improve on these pathfinding results [40,42,46,47]. At higher masses the MADMAX Collaboration will search for dark matter using a dielectric haloscope, which searches for the electromagnetic emission that an axion generates at dielectric boundaries in the presence of a magnetic field [48][49][50]. Other proposed instruments searching for dark matter through the axion-photon coupling include resonant frequency conversion in superconducting cavities [51,52], looking for a phase difference in locked lasers [53,54], exciting quasidegenerate modes in a superconducting cavity [55], detection of small energy deposits in crystals [56,57], and matching the axion mass to a plasma frequency [58,59], although this list is far from exhaustive. In summary, it is likely that in the coming decades the axion dark matter hypothesis will either be confirmed, or required to satisfy g aγγ g SE aγγ in the mass range 1 neV m a 1 meV. Let us now sketch how this progress in the search for axion dark matter can be repurposed to search for the CaB. The detectable power deposited by an axion population via (1) is naively ∝ g 2 aγγ ρ a up to the details of the experimental readout. Taking the experimental factors as constant between the two scenarios, we can obtain an estimate of the sensitivity for a dark-matter instrument to the CaB by matching the power between the two, i.e. (g 2 aγγ ρ a ) CaB = (g 2 aγγ ρ a ) DM . Assuming axions fully constitute dark matter, astrophysical observations determine that ρ a = ρ DM 0.4 GeV/cm 3 , see e.g. [60]. 1 The unknown parameter being searched for is then g aγγ , which beyond g aγγ ≤ g SE aγγ is a free parameter, although in specific scenarios like the QCD axion sharper predictions are possible. Regardless, for a given instrument we can project the reach in g aγγ to determine the associated reach in deposited axion power. For the CaB both g aγγ ≤ g SE aγγ and ρ a are free parameters. If the 1 We take ρDM = 0.4 GeV/cm 3 throughout.
CaB is a relic of the early Universe, measurements of ∆N eff further require the energy density to be less than that of the Cosmic Microwave Background (CMB) [19], ρ a ρ γ , although the density may be predicted in certain scenarios. This poses an immediate challenge: for equal g aγγ , the power deposited by the CaB will be at least a factor of ρ DM /ρ γ 10 9 smaller. The situation is even more dire. The detectability of power deposited by axion dark matter is enhanced by the exceptionally long coherence time of the signal, originating from the narrow energy distribution associated with non-relativistic dark matter in the Milky Way. Indeed, for dark matter we expect ∆ω/ω ∼ 10 −6 , whereas generically the CaB will have a broad distribution in energy, ∆ω/ω ∼ 1. As we will review, this typically enhances sensitivity to the dark matter signal by a further three orders of magnitude relative to the CaB. Accordingly, for equal g aγγ , a relativistic axion that is a relic of the early Universe is at most a trillionth as detectable as dark matter.
As we will show, the challenge is not insurmountable. Upcoming axion dark-matter instruments will have a sensitivity that such a CaB will be detectable. This is demonstrated in Fig. 1, where we recast the existing results of ADMX and HAYSTAC and the expected future reach of DMRadio and MADMAX onto the equivalent CaB parameter space, assuming g CaB aγγ = g SE aγγ . In detail, we define Ω a (ω) as the relic density per unit log (angular) frequency of the axion with ρ c = 3M 2 Pl H 2 0 the critical density. 2 Fixing the coupling, we can recast the stated dark-matter sensitivity to a sensitivity on ρ a and hence Ω a (ω). The figure demonstrates that DMRadio will be sensitive to scenarios where Ω a (ω) 5 × 10 −5 , roughly corresponding to ρ a < ρ γ -a target cavity instruments may also reach in the futureand further the H 0 preferred parameter space discussed earlier.
The sensitivity to such small energy densities suggests that the data collected by axion direct detection experiments can be repurposed to probe a range of cosmic sources of axions beyond non-relativistic dark matter. The axion distribution can be narrow or broad and have a peak frequency over a wide range of energies, depending on how and when they were produced, which motivates the discussion in this work on mechanisms for generating the CaB. In particular, we discuss a thermal axion background, emission from cosmic strings, and production from a parametric resonance in the early Universe, which is expected to produce a roughly Gaussian distribution, all of which are shown in Fig. 1. For such cosmic relics, the axions will free-stream over long distances and their spectrum ultimately depends on the cosmic history of the Universe. In this sense, axion experiments looking for a stochastic axion background are in close analogy with searches for a stochastic gravitational wave background (only axions may have a much larger coupling). 3 While ADMX is close, no existing instrument is currently sensitive to cosmological relics. This motivates scenarios where the CaB is produced in the late Universe, where ρ a can be larger than ρ γ , and in particular we discuss dark matter decaying to two relativistic axions, χ → aa, with m a m χ /2. The resulting spectrum of axions receives two contributions: one from the decay of dark matter within the Milky Way, which generates a sharp ∆ω/ω ∼ 10 −3 spectrum, and the broader spectrum resulting from dark-matter decays throughout the Universe. Both contributions can be seen in the spectrum shown in Fig. 1. As we will show, dark-matter instruments can be repurposed into axion telescopes to search for this dark-matter indirect-detection channel. Such searches can further exploit the fact that the Milky Way signal will undergo a daily modulation in microwave cavity instruments, as the relative direction of the signal, primarily from the Galactic Center, and the experimental magnetic field vary throughout the day. Indeed, we will show that ADMX is currently sensitive to unexplored parameter space -a reanalysis of their existing data may already reveal a signal of the CaB.
In the remainder of this work we will expand the above discussion as follows. To begin with, in Sec. II, we introduce different possible CaB sources focusing on thermal production, dark-matter decay, parametricresonance production, and emission from topological defects. Then, in Sec. III we study the viability of detecting a relativistic axion background with instruments designed to search for dark-matter axions through the axion-photon coupling. As already mentioned, we focus on the axion-photon coupling, and further will restrict our attention to the sensitivity with resonant cavity instruments such as ADMX and HAYSTAC, and lumpedcircuit readout approaches such as DMRadio. Our analysis will justify the sensitivities shown for these instruments in Fig. 1. We will not, however, return to carefully consider the sensitivity of instruments focused on higher mass axion dark matter m a ∼ ω > 100 µeV, such as MADMAX. Detecting a relic CaB requires sensitivity to g aγγ many orders of magnitude below g SE aγγ . This simply will not be achieved in any proposed high mass instrument. 4 In Sec. IV we then combine the results to determine projected limits on various CaB scenarios, and finally present our outlook in Sec. V.

II. CaB SOURCES
We now turn to a discussion of specific production mechanisms for the CaB. As mentioned already, axions can be produced in the early and late Universe, and we will consider examples of both. In each scenario, our goal will be to characterize the associated axion energy spectrum, which will be a central ingredient when we come to detection. For this purpose we will again use Ω a (ω) as defined in (2). We emphasize that the present discussion is not intended to be an exhaustive consideration of all scenarios from which a CaB could emerge, rather, we simply demonstrate that there are many possibilities. Nonetheless, the analysis will reveal a common theme that emerges across production mechanisms, in particular that the CaB will generically be a broad distribution, ∆ω/ω ∼ 1. When contrasted with the highly coherent signal predicted for dark matter, this expectation will represent a fundamental difference when approaching searches for relativistic axions.

A. Thermal Relic
We begin by studying the simplest example of a CaB source, thermal production during the early Universe. Early studies of thermal axion production can be found in [65][66][67][68][69] with a more detailed analysis performed in [70][71][72]. Fundamentally, if an axion was ever in thermal contact with the SM bath at high temperatures, then a residual thermal population is expected to exist to the present day, generating a CaB with the closest resemblance to the CMB. Indeed, a thermal axion relic will also be described by a blackbody spectrum, so that with a total energy density comparable to that of the CMB. The above distribution is defined by a single parameter, the present day axion temperature, T a . 5 Remaining agnostic as to the exact axion-SM interaction that brought the axion into thermal equilibrium initially, at some temperature, T d , the two will decouple. If we assume that since axion freeze-out there has not been any entropy dilutions beyond those in the SM, and further that there was not an early period of matter domination, then as entropy is approximately conserved, we can relate the present and decoupling axion temperatures as follows, Here T 0 2.7 K is the present day CMB temperature and g * ,S (T ) is the number of entropic degrees as a function of temperature. Accordingly, we can specify the thermal axion in terms of T a or T d . The spectrum for different decoupling temperatures is shown in Fig. 2. Generically, a thermal distribution is associated with a number density of axions of O(1 − 100 cm −3 ) and a peak energy of ∼ 10 −4 eV. Again, both are comparable to the CMB. As the figure shows, T d 1 MeV is excluded by ∆N eff measurements, which would include the case where T a = T 0 as shown in Fig. 1. Nevertheless, the range 1 MeV T d 1 GeV is somewhat favored as a solution to the present H 0 tension, a possibility that was studied in detail in Ref. [73].
Ultimately, T d is determined by the microphysics responsible for the axion coming into thermal contact. Axions coupled to photons have a maximum possible decoupling temperature since processes such as, γe → ae, will keep it in equilibrium. Equating this rate with Hubble leads to an estimate of the decoupling temperature, We conclude that axions saturating the star-emission bounds would have a decoupling temperature of around a TeV, while additional interactions can keep the axion thermally coupled at lower temperatures. This motivates a range of decoupling temperatures. Lastly, we emphasize that a thermal abundance of axions will always form as long as the temperature of the Universe was ever above the decoupling temperature, making this population a robust prediction for any theory without a low reheating temperature after inflation.

B. Dark-Matter Decay
Dark matter need not be absolutely stable, and axions offer one possible decay channel. Provided that the dark-matter mass is significantly larger than m a , then the axions produced through this process will be relativistic. If these same axions have a sizable photon coupling then they are in principle detectable in terrestrial experiments, opening up a new channel for the indirect detection program. An important aspect of this scenario, is that since the dark-matter abundance is considerably larger than the CMB (ρ DM /ρ γ 10 9 ), decaying dark matter can result in a CaB energy density that is significantly larger than allowed for a relic population, given bounds from ∆N eff . 6 Accordingly, in the short term decaying dark matter represents the most accessible CaB candidate. Currently, this scenario is only constrained indirectly as a consequence of the fact that a significant fraction of dark matter decaying into radiation would modify the expansion history of the Universe [74,75]. Qualitatively, these bounds require the decay rate to be less than the current Hubble rate Γ H 0 , so that the dark-matter lifetime is longer than the age of the Universe, τ = 1/Γ t U . As dark matter decaying to a relativistic species would modify the expansion history, this possibility has been suggested as a potential resolution to the Hubble tension [76] (though this was later refuted [77]). Given the existing tension in H 0 measurements between the early and late Universe, the current bounds depend noticeably on which data set is used. Recently, using only local measurements, the Dark Energy Survey constrained the dark-matter lifetime to τ 50 Gyr ∼ 3.6 t U [78] and we consider this as our nominal bound. Although this value is likely to be revised with developments on the Hubble tension, this will not qualitatively impact our discussion.
The first goal of this section will be to describe the CaB that results from decaying dark matter. We will then outline an explicit decaying dark-matter model, and lastly, we discuss the feasibility of detecting axions arising from the related mechanism of neutrino decays.

The Axion Spectrum from Decaying Dark Matter
Axions produced from dark-matter decay will have a spectrum that results from two distinct sources: the decay of galactic dark matter within the Milky Way and contribution from decays of the extragalactic dark matter throughout the Universe. While in both cases the fundamental process will be dark matter, which we denote χ, decaying to axions, the resulting spectra will be significantly different. Nonetheless, the contributions produce similar axion abundances.
Consider first the extragalactic contribution resulting from dark matter decaying to axions throughout the isotropic, homogeneous, and expanding Universe. The number density of axions observed today produced per unit time and per unit energy is given by the product of several factors. The first of these is the number density of dark-matter particles at a given time t, ρ DM (t)/m χ . We must also weight this by the rate at which dark matter decays at this time, which is Γe −Γt (we assume that the decay rate Γ is constant through cosmic history). Each decay is associated with a differential energy spectrum of the emitted axions, dN/dω , normalized such that its integral over all ω gives the number of emitted axions. As the emitted axions are assumed to be relativistic, the axion energy as observed today will be suppressed by a ratio of scale factors, ω = ω a, where a is the scale factor at time t and we take a 0 = 1. Finally, the present number density will be diluted as compared to the density emitted at t, as the Universe is now larger by a factor of 1/a 3 . Combining these factors and then integrating over all time from t = 0 to the present t = t 0 , we obtain the total extragalactic differential number density as Changing integration variables to the scale factor, we can write this as with Ω DM 0.27 the cosmological dark-matter density and t(a) the age of the Universe as a function of scale factor, so that t(1) = t U .
The above result is appropriate for a general darkmatter decay axion spectrum. Throughout this work, however, we will specialize to the simple example of a two-body decay, χ → aa, again assuming m a m χ /2, so that the produced axions are relativistic. In this case, the spectrum takes the following particularly simple form, Inserting this into the above gives, where Θ is the Heaviside step-function and we have exchanged the decay rate for the lifetime. The axion energy spectrum as observed today is shown in Fig. 3 for different dark matter lifetimes and masses. The sharp peak is associated with galactic decays, described shortly, but the broad continuum arises from the above expression. The redshifting of the axions produced throughout the Universe smooths the sharp two-body spectrum into a continuum. We now compute the axion abundance created within the Milky Way. For this, the number density of axions per unit energy around Earth is determined by the conventional indirect detection expression for decaying dark matter, and given by 7 Note we have assumed the observable decays within the Milky Way occur at t = t U . The integral at the end of this is commonly referred to as the D-factor in the indirect detection literature (for details see e.g. [79]). For a canonical Milky Way dark-matter profile, the Dfactor has a full sky integrated value of D MW 2.7 × 10 32 eV/cm 2 · sr. 8 The expression in (9) again holds for a general spectrum. If we specialize to the case of χ → aa, then the local axion energy density per unit log frequency is approximately given by, 7 For τ t U , we have e −t U /τ 1, and so this factor is commonly neglected in indirect detection analyses. 8 To obtain this value we assumed a canonical Navarro-Frenk-White profile [80,81], took the Earth-Galactic Center distance as 8.127 kpc [82], and a local dark-matter density of 0.4 GeV/cm 3 .
Integrating (8) and (10), we can determine the total energy density for the two contributions. This is maximized for τ ∼ t U , where we have ρ MW a 2ρ EG a 10 3 ρ γ , so that, as claimed, the energy densities from the two contributions are comparable, and a combined density larger than the CMB can be obtained for a range of lifetimes (ρ a ρ γ for τ 10 4 t U ).
The reason (10) is an approximation is because it assumes the observed axion spectrum is the same as that in the dark-matter rest frame. While this is often a reasonable approximation, axion experiments are often sensitive to extremely narrow energy distributions -recall, for dark matter, ∆ω/ω ∼ 10 −6 . This motivates a more detailed consideration of the axion energy distribution. There are two contributions that will resolve the distribution in (10) to have a finite width: the velocity dispersion of dark matter in the Milky Way and the finite velocity of the Earth through the dark-matter halo. Both velocities result in a net motion between the observer and source of axions, and therefore the axion energies will be Doppler shifted by a factor of v ∼ 10 −3 , which is the magnitude of both velocity components. In this work, we will simply replace δ(ω − m χ /2) in (10) with a Gaussian of width 10 −3 centered at half the dark-matter mass. The actual distribution is more complex, indeed it depends on the dark-matter distribution and varies across the sky given the motion of the Earth in the halo frame (for further details, see Ref. [83]). Nevertheless, the main aspect of the distribution relevant for forecasting sensitivities is the width, and the Gaussian approximation adequately accounts for this.

An Explicit Model: Decaying Scalar Dark Matter
Above we considered the axion abundance produced through dark-matter decay, with all model dependence in the axion spectrum, dN/dω, and the lifetime τ . We now study a decaying dark-matter model which predicts a detectable CaB, and generates the simple χ → aa spectrum used above.
Consider a theory with a complex scalar field, Φ, with potential, The theory has a spontaneously broken U(1) and we identify the Goldstone boson with the axion and the radial mode with dark matter, decomposing the field as Φ = (χ + f a )e ia/fa / √ 2.
In the broken phase, the relevant axion dark-matter couplings are, 9 from which we identify the dark-matter mass as m χ = √ 2λf a . Further, the axion dark-matter coupling allows us to compute the rate of dark-matter decay as, with corresponding axion spectrum as given in (7). In order for Γ χ→aa to be, at least, comparable to the age of the Universe, we then require f a to be well below the weak scale. This may seem hard to reconcile given the stringent bounds on the axion-photon coupling, g SE aγγ 1 TeV −1 , however, this can be natural if the axion obtains its photon coupling through axion-axion or photon-dark-photon mixing, and as we demonstrate in App. A this does not require any elaborate model building. Nevertheless, this does require forbidding any significant terms in the scalar potential mixing between the SM Higgs and Φ. In generating g aγγ , χ may also obtain a coupling directly to photons. Even though searches for χ → γγ are significantly more stringent than the axion searches discussed in this work, these constraints are not significant in the parameter space of interest, as we work in the limit of f a g −1 aγγ .

Cosmic Neutrino Background Decay
A CaB may also be produced as a byproduct of neutrino decays of the cosmic neutrino background (CνB). Flavor off-diagonal couplings of axions to neutrinos can be the result of a global lepton number broken by multiple scalar fields with the axion playing the role of a Majoron [84] or, more generally, the familon [85]. A generic axion can have a coupling to neutrinos given by (Q = Q † ), Here we assume the neutrinos are Majorana and work with two-component fermion notation. The neutrino decay rates were recently calculated in Ref. [86] for L i − L j gauge bosons in the high-energy limit and the results can be translated to axions with the relation, 1/f a ↔ g X /m X , Parametrically, Γ νi→νj a ∼ m 3 ν Q 2 /f 2 a and for the decay to be comparable to Hubble while avoiding the staremission bounds requires f a 1 TeV while keeping g aγγ 1 TeV −1 . As mentioned previously, this can occur naturally for axions that inherit a photon interaction through axion-axion or photon-dark photon mixing, see also App. A.
The strongest constraints on the neutrino lifetime are from observations of neutrino free-streaming in the CMB [87]. Current limits allow for neutrino lifetimes well below the age of the Universe; indeed, a recent reanalysis of the bounds in Ref. [88] found a conservative limit that is on the order of several days. Accordingly there is a significant possibility that neutrino decays populate the CaB. For decays while the neutrinos are still relativistic, the axions are produced with an energy comparable with the neutrino temperature, and hence this results in a spectrum similar to the thermal background considered in Sec. II A. Axions produced from late-time neutrino decays -after neutrinos have become non-relativistic -will have a peaked spectrum around the neutrino mass. In either event, the resulting spectrum will be subject to the same challenge as the thermal background, in that it is located at high frequencies where it is unlikely to be observable in the near future due to the lack of sensitive experiments at these energies. As such, we will not evaluate this case in detail, but note should the thermal CaB become accessible, then likely so too would this scenario.

C. Parametric Resonance
A CaB can also be produced through the process of parametric resonance in the early Universe [89][90][91][92]. In order for this process to occur, the axion must be coupled to a scalar field which is heavily displaced from its minimum after inflation. This will occur by quantum fluctuations for any scalar field, unless it is fixed to the origin by an effective mass larger than the Hubble scale at inflation. When such a scalar field begins to oscillate about its minimum it will produce axions with a bose-enhanced rate that will typically deplete its energy density within an e-fold. The characteristic axion energy as observed today will have redshifted dramatically and can be much lower than the energy of axions produced by perturbative decay of a scalar field. As the parametricresonance phenomena is a non-perturbative process that occurs out of equilibrium, computing the spectrum in detail requires evolving multiple scalar fields on the lattice. We will not attempt such a calculation here but instead perform qualitative estimates. Earlier work on relativistic axion production considered potential modifications to ∆N eff and parallel production of gravitational waves [93]. In this subsection we review the dynamics of axion production and explore the parameter space leading to a detectable axion background. We follow the notation and discussion of Refs. [94,95] which studied the prospect of ultralight bosonic dark matter produced through parametric resonance.
We now focus on an explicit realization of the parametric-resonance phenomena, which can be achieved using the same model introduced in Sec. II B. Recall, there we had the axion arise from a global symmetry breaking complex scalar, Φ, with a radial mode χ playing the role of dark matter. Our starting point will again be the potential given in (11) and as our initial condition we take χ to have a large field value, χ i f a . At early times we assume the second derivative of the potential with respect to the field is greater than Hubble squared, V (χ i ) H 2 , such that the scalar is stuck and the field redshifts as vacuum energy. When V (χ i ) ∼ H 2 the field begins to oscillate, resulting in exponential production of both the radial and axion modes. Since χ i f a , m 2 χ χ 2 is small relative to λχ 4 and can be neglected during the oscillations. This leads to a broad resonance that rapidly depletes the energy density stored in the original scalar field. Furthermore, since the axion energy at the time of production is set by the effective mass of χ, m eff χ (χ) ≡ λχ, it is independent of the temperature of the SM bath and often considerably smaller. Subsequent redshift until today can lead to relativistic axions over a wide range of energies, well below the temperature of the CMB and potentially within reach of low-frequency axion haloscopes.
We now estimate the abundance and energy spectrum of the axion and radial mode. Axions are emitted with energy ω a ∼ m eff χ (χ i ) during radiation domination at a temperature of oscillation, T osc ∼ √ ΓM Pl , where Γ is the oscillation timescale. For perturbative production, Γ cannot be arbitrarily large, in detail Γ m 3 χ /f 2 a . For parametric resonance, the particle production will occur within a few oscillations, so Γ ∼ λχ i . The characteristic axion energy, as measured today, is redshifted using T osc and given by, where s(T ) is the entropy of the SM bath. Accordingly, provided m eff χ (χ i ) M Pl , we haveω a T 0 , and the axion energy will be well below the cosmic photon temperature.
With an estimate for where the spectrum will peak, next we consider the a and χ comoving number densities after the oscillations have concluded. These can be parameterized as, where f denotes the fraction of energy transferred to axions,ω a andω χ are the mean energies of each particle at the time of production, and ρ χ,osc 1 4 λ 2 χ 4 i is the total energy density in the radial direction prior to oscillations. Since the vacuum mass of radial mode can be neglected in this limit both a and χ are produced with comparable energy densities, f 1/2 and comparable energies,ω a ω χ λχ i . As noted above, the energy is determined by the effective mass, which is driven by the quartic. These determine the comoving number densities to be so that the axion energy density today is fixed by the initial scalar field value, We conclude for χ i < M Pl , parametric resonance produces a maximum axion relic density a few orders of magnitude below that of the CMB. The characteristic frequency of the axions may span many orders of magnitude, and depends on the initial field value as well as the quartic. To explore the parameter space it is helpful to focus on a specific case where χ makes up dark matter. 10 Requiring χ to have the observed dark-matter abundance provides an additional constraint, To study the resulting model space, we take the free parameters to be the three parameters {m χ , χ i , λ}, one combination of which is restricted by the requirement of (20). Note, the value of the vacuum expectation value is a dependent parameter, f a = m χ / √ 2λ. This leads to a prediction for the axion energy today, The relative spread in the axion spectrum must be determined using lattice simulations though we expect it to be O(1). Simulations for a similar theory have found the spectrum to be roughly a Gaussian with a relative width of order unity [96]. Here we have only included the influence of redshift on the axion energy spectrum. It is known that there are additional number-changing processes that tend to move the axion spectrum toward a thermal distribution. These are slow and not expected to effectively thermalize axions on a cosmological timescale, but may shift the peak axion frequency by an order of magnitude [96].
FIG. 4. The parameter space of parametric resonance producing axions and χ, where everywhere in the plot χ constitutes all of dark matter. The constraints shown arise from darkmatter stability, warmness, efficient parametric resonance, and isocurvature. Further, we show two contours of axion energy density, and separately the mean CaB energy. The higher mean energy ofωa = 10 −9 eV falls entirely in the region excluded by DM stability, however this constraint is removed if we no longer assume χ constitutes the dark matter of the Universe.
The parameter space of PR production, assuming χ makes up dark matter, has several constraints summarized below: 11 1. DM Unstable: As discussed in section II B, χ may also decay (perturbativity) into axions with a rate given by (13). As discussed earlier in the context of dark-matter indirect detection (see Sec. II B), we use the nominal bound of τ 50 Gyr ∼ 3.6 t U [78].

2.
Warmness: This bound arises from the requirement that χ are cold enough to constitute dark matter today, roughly taken to be p χ (T )/m χ 10 −3 at recombination. In detail, we require 3. Inefficient PR: Parametric resonance is efficient at producing axions when the initial field value is much larger than its vacuum value. Otherwise the resonance is a narrow feature and is unable to convert all the energy density in χ into field excitations. For this condition we take the rough bound, 4. Isocurvature: During inflation we assume χ has an effective mass below the Hubble scale at inflation, λχ i H inf such that fluctuations during inflation displace the field away from the minimum. These isocurvature perturbations can be observed in the CMB, placing a bound χ i /H inf (π 2 βP R (k * )) −1/2 , where β ≤ 0.011 is the isocurvature fraction and P R (k * ) 2.1 × 10 −9 is the observed amplitude of the curvature power spectrum at the pivot scale [97]. Combining these results places a limit on the scalar quartic, λ.
The parameter space of relativistic axions produced through parametric resonance in light of these bounds is shown in Fig. 4, where we have fixed the abundance of χ to match that of dark matter. We see that for χ to constitute dark matter, we require m χ 1 keV and the condition of a detectable CaB further restricts χ to large initial field values and smaller masses.
These results demonstrate that a consistent CaB produced through parametric resonance can occur over an enormous range of frequencies. From (19), detectability will be maximized for χ i ∼ M Pl (up to consistency of the warmness criteria). The spectrum is then expected to be roughly an O(1) width Gaussian [96], with peak frequency determined from (21) asω a 5 × 10 −7 eV (m χ /1 eV). In the scenario where χ constitutes dark matter, this allows a mean energy as low asω a ∼ 10 −28 eV, when m χ ∼ 10 −22 eV is in the fuzzy dark-matter regime, and as high as 5 × 10 −12 eV when saturating the dark-matter stability criteria, shown in Fig. 4. Removing the requirement that χ be dark matter, the frequency range can then be extended even further, particular to higher frequencies which may be accessible by DMRadio, or even ADMX and HAYSTAC.

D. Topological Defect Decay
A CaB may also be produced through the decay of topological defects. In this section we study the abundance and energy spectrum of axions emitted from a network of cosmic strings formed during a thermal phase transition. Axions produced during the phase transition itself are in thermal contact with the SM bath and will contribute to the thermal background discussed in Sec. II A. We focus on the axions produced after chemical decoupling, when the cosmic-string network has entered the time where the network has a constant number of strings per Hubble volume (up to log-violations, which we also take into account), commonly referred to as the "scaling regime". In computing the spectrum we work in the limit m a → 0 where strings remain until late times and there are no domain walls. If there is a finite mass, and the domain-wall number is equal to 1, the network will quickly collapse when H ∼ m a . This will produce an additional burst of axion production and a sharp drop in the axion spectrum at a characteristic frequency. We do not consider these effects but they may produce additional distinctive signals.
The spectrum of axions emitted by cosmic strings is still an active area of discussion in the literature with the debate centered on whether the typical axion energy emitted by a string is of order the inverse length or inverse thickness of the string (in particular, see Refs. [98,99] and [100,101]). We estimate the abundance and spectrum following numerical simulations done for the QCD axion in Refs. [98,99], where the simulations suggest that the spectrum is dominated by low-energy axions.
To begin with, the energy density of cosmic strings can be parameterized using the average length of string within a Hubble length, ξ, and the energy of that string, given by the product of its tension, µ eff , and a Hubble length, 1/H. The total energy is then averaged over Hubble volume, 1/H 3 . Following Ref. [98] we write this as, This form is convenient since both ξ and µ eff only evolve logarithmically with time. Their evolution and can be parameterized as [98] Here m r ∼ f a is the string width, γ is roughly a constant in time which we will approximate as unity, α 0.24 ± 0.02 [99] is also a constant, and finally we take β 0 since we are interested in late times, where the log term is the dominant contribution.
The rate of axion energy emission during the scaling regime per unit volume, Γ, is given by the difference of the energy density of "free" strings (strings without inter-commutation and radiation) and the energy density stored in strings, We assume that both energy densities are equal at an initial time, t i . The free energy density at a later time t is then, This follows as ρ free s ∝ t −1 and we require ρ free s and ρ s to match at t = t i . Inserting (23) and (26) into (25) and working in the large log limit gives, FIG. 5. The energy density (left) and spectrum (right) of a CaB produced from the emission of cosmic strings. The relic density is shown as a function of the axion decay constant for different decoupling temperatures, T d . The spectrum is shown as a function of energy, however, the shape at high frequencies is sensitive to the decoupling temperature, which we take to be equal to fa (solid) or fa/10 3 (dashed). In both figures we also show bounds from ∆N eff and the region preferred to mildly alleviate the Hubble tension. On the left, we further show the region excluded by the maximum possible reheating scale of the Universe, derived by assuming instantaneous reheating from the maximum scale of inflation.
where ρ SM is the total energy density in the SM and the combination has only a logarithmic time dependence. The relic density in axions today is then given by, where a d is the scale factor at the time the network enters the scaling regime. This expression applies during both radiation and matter domination, however we note that the simulations to estimate ξ were only performed for radiation domination, and we focus on axions produced during this epoch. The relic density for different decoupling temperatures is shown on the left of Fig. 5. In addition to bounds on the energy density, cosmic strings have a constraint on the maximum value of the decay constant. The energy scale of inflation is given in terms of the tensor-to-scalar ratio as [102], Using the upper bound on r < 0.056 [97] and setting V = 3M 2 Pl H 2 I we can derive an upper bound on the Hubble scale of inflation, H I < 6 × 10 13 GeV. Relating H I to the reheat temperature of the universe, gives a maximum possible T RH . To have a thermal phase transition in the early universe requires T RH to be above the critical temperature for a phase transition, ∼ f a , putting an upper bound on the decay constant. Lastly, we note that if the cosmic string network remains until recombination there is an additional bound from the string energy density imprinted on the cosmic microwave background [103]. The spectrum of axions produced after this epoch corresponds to frequencies, ω 10 −31 eV, and are not observable with the experiments considered in this work. In presenting our axion spectrum and experimental projections, we assume the network collapses before this time.
We now move on to calculate how this energy is distributed. The emission spectrum of axions from strings has been a source of uncertainty in the literature with the debate centered around whether axion emission is dominated by coherent motion of the string producing axions with wavelength of order the string length ("IR-dominated") or by small loops and kinks along the string producing axions with wavelength of order the string width ("UV-dominated"). This has profound consequences for QCD axion dark matter as it predicts a relic abundance produced from topological defect decay with an uncertainty of a few orders of magnitude (see Refs. [98,99] when compared to Ref. [100]). Fundamentally, the enormous separation of scales between the string length and its width make this a challenging problem to resolve. In either case, the spectrum can likely be approximated by a power-law parameterized by a spec-tral index q with a high and low energy cut-off [98], normalizes F such that the integral over all x is unity. Here x is the appropriately normalized energy, while x 1 and x 2 are IR and UV cutoffs. In Fig. 6 we show xF (x) for different values of q, demonstrating that for q < 1 the spectrum is UV-dominated while for q > 1 the spectrum is IRdominated. While the debate is yet to be settled (in particular, see Ref. [101]), a recent analysis [99] suggests that the spectrum is IR dominated during the scaling regime with estimates for the IR and UV cutoffs of x 1 10 and x 2 m r /H respectively. Furthermore, the authors of Ref. [99] find a best fit for the spectral index over time of, and we assume this form in our results. Given a spectral function F (x), the axion spectrum as observed today is the appropriately weighted timeintegral of this expression, where ω is the axion energy as measured today. The spectrum for different f a is shown on the right of Fig. 5. At low energies the spectrum is roughly a constant (a consequence of the network being in the scaling regime) while at high energies the spectrum falls off as it relies on producing axions with energies much larger then the Hubble scale from cosmic-string oscillations. The frequency where the drop begins depends on the decoupling temperature of the axion with the SM bath, T d , with solid lines denoting T d = f a and dashed curves showing T d = f a /10 3 . In all cases, the abrupt change just below ω = 10 −22 eV is associated with a drop in the cosmic string energy density after the QCD phase transition, where the number of relativistic degrees of freedom in the SM drops considerably. As for all CaB candidates, in addition to the dependence on the energy density and spectrum, axion detection from cosmic strings is sensitive to the axion-SM coupling. For generic axions, the axion-photon coupling is g aγγ α/2πf a , which when combined with the densities above would be challenging to observe. Accordingly, when we discuss the experimental prospects, we will again consider g aγγ larger than this simplest expectation, which can be induced by mechanisms including the clockwork [104] (see Refs. [105] and [106] for recent summaries of such mechanisms in the context of the QCD axion and ultralight axion dark matter, respectively). Nonetheless, we note that if experiments could probe the scenario where g aγγ ∼ α/2πf a , it would be possible to FIG. 6. The function in (31), which determines the spectrum of axion emission from cosmic strings for different spectral indices, q. If q < 1 the emission spectrum is dominated by axions with a wavelength of order the string width while if q > 1 the spectrum is dominated by modes of order the string length. These two scenarios are referred to as UV and IR dominated, respectively. We assume an IR dominated spectrum in this work as suggested in Refs. [98,99].
probe all value of f a . This is because the detectable CaB power is ∝ g 2 aγγ ρ a , and for cosmic strings we have ρ a ∝ µ eff ∝ f 2 a , resulting in g 2 aγγ ρ a being independent of f a .

III. DETECTING THE CaB
Having motivated the possibility of a local CaB, we now turn to the question of how that population could be detected. We focus on detection at a few of the many instruments constituting the burgeoning program to detect ultralight dark matter. Our central conclusion will be that experiments designed with axion dark matter in mind are generally also sensitive to a relativistic population. Indeed, it is possible that ADMX has already collected a detectable signal that would have been missed by an analysis focused on the non-relativistic axion.
Qualitatively, detection of relativistic axions proceeds as for their non-relativistic counterparts. In both cases, the axion can be described as an oscillating classical wave, 12 which through a coupling to the SM induces a 12 The classical wave description holds in the limit of a large number of states per de Broglie volume, naλ 3 dB 1. For dark-matter axions, using the mean expected dark-matter density and speed, this is satisfied for ma 10 eV. For the CaB, we instead have naλ 3 dB ∼ (ρa/ργ )(ω/1 meV) −4 . In the present work we will consider detection exclusively in scenarios withω 1 meV and sufficient densities that classicality applies. Nevertheless, our description will not apply for arbitrarily large mean energies or detectable time-varying signal in, for example, electromagnetic waves or nuclear spins. A central difference is the signal bandwidth. For dark matter, the expectation is that the signal power will be deposited in an extremely narrow range of frequencies centered around its unknown mass m a . In general, the oscillation frequency is set by the axion energy. For a non-relativistic particle, the energy is ω m a (1 + v 2 /2), and given our expectation for the local dark matter is that the speeds vary over a range ∆v ∼ 10 −3 and take a mean valuē v ∼ 10 −3 , axion dark matter carries a large quality factor of Q DM a =ω/∆ω ∼ 10 6 , whereω is the average energy. For a CaB the expectation is that the local axion field has a wide distribution of energies, such that generically Q CaB a ∼ 1. An exception is dark matter decaying to axions within the Milky Way, where we expect Q CaB a ∼ 10 3 . Regardless, in either case Q CaB a Q DM a , and this will represent a challenge to detection. There are additional important differences between the relativistic and darkmatter cases -for instance, the relativistic signal can exhibit a unique daily-modulation signal even at a single detector -and we will explore these as well. As is the case for the bulk of this work, we restrict our attention to experiments focused on the axion coupling to electromagnetism, g aγγ , although much of our formalism can be lifted for other SM couplings.
We divide our discussion of the CaB detection into five parts. Firstly we outline several basic features of a relativistic axion population -its expected amplitude and distribution in both time and frequency -applicable to any detection strategy. We next use these results to sketch our expected sensitivity to the CaB by comparing the experimentally detectable power associated with the relativistic and dark-matter axion field. Having provided a general sensitivity estimate sufficient for understanding Fig. 1, we then focus specifically on the axion-photon coupling, detailing axion electromagnetism with a specific focus on the differences in the relativistic case. Finally, we apply these lessons to existing axion dark-matter detection strategies, and discuss representative examples of both broadband and resonant detection strategies approaches. In the following section we will use these results to set estimated limits on a number of different CaB scenarios discussed in Sec. II.

A. Properties of the Relativistic Axion
To start our discussion, we will outline general properties of the relativistic axion field relevant to its detection. In particular, we treat the CaB as the superposition of many non-interacting axion particles with energies ω drawn from probability distribution p(ω). 13 Within small densities. The approximate boundary between the two regimes is shown in Fig. 1. 13 Formally we can define p(ω) = (1/na)dna/dω, as it is the differential number density that controls the probability of observing this framework, we will derive the expected amplitude of the axion field a -more specifically of a 2 -in both the time and frequency domain, and further quantify the fluctuations around the central value. Any experimental detection will involve a coupling to the axion field, and therefore the measurements will inherit these average values and fluctuations. We will consider a general energy distribution, and show that our results contain the nonrelativistic limit as a special case. Indeed, our results are a direct generalization of the non-relativistic field, which will allow us to bootstrap known dark-matter results to the CaB. For the energies and densities considered in this work, the axion field will always contain an enormous number of particles per de Broglie volume. Consequently, the field can be described in terms of an emergent classical wave. In this respect, the CaB directly mirrors non-relativistic dark matter for m a 10 eV, where the associated statistics were derived in Ref. [107], and we will generalize a number of results from that reference. We imagine the classical axion wave as constructed from a large number, N a , of non-interacting waves, Each element of this sum is associated with a random variable ω i , an energy drawn from p(ω). Beyond their energy, however, there is no reason to imagine the various states are phase coherent, and this is ensured by the uniform random variable φ i ∈ [0, 2π). The amplitude is fixed by ensuring the field carries energy density ρ a , which would be equal to ρ DM for non-relativistic dark matter. In the discretized picture,ω = N −1 a i ω i , but more generally we takeω = dω ωp(ω).
In principle, there is an additional contribution to the phase neglected in (34): the spatial variation controlled by −k i · x, where k i is the particle momentum. If we imagine measuring the axion field at a single point, this contribution is irrelevant at the level of the phase, as we can always center our coordinates such that x = 0. 14 Yet where it can be relevant is through effects sensitive to the spatial gradients of the axion field. As we will discuss in Sec. III C, whilst these gradients are usually neglected for a non-relativistic field, they are parametrically important for a relativistic population. The amplitude of these effects is fixed by the massive dispersion relation, |k i | = ω 2 i − m 2 a , and therefore also controlled by the energy distribution, p(ω). The direction of k, however, is not, and will itself be drawn from a distribution on the celestial sphere. For instance, if the dark matter is made of axions, the direction of k may point an axion at a given energy. 14 If the axion field is measured at multiple spatially separated locations, however, the k · x contribution to the phase is physical, and can be used to perform interferometry on the wave [108]. towards a dark-matter stream incident on the Earth, or in the relativistic case the CaB would be biased towards the center of the Milky Way if it originates from darkmatter decay. As shown in [108], the angular distribution can be fully incorporated into the description of the nonrelativistic axion field, and the arguments there can be generalized to relativistic axions. 15 We will not pursue this direction in the current work, however. While the directional distribution will be relevant for the fine details of the relativistic signal -in particular as it relates to daily-modulation effects unique to the relativistic axion, discussed in Sec. III C -it unnecessarily complicates an estimate of the experimental reach, which is our focus.
To make progress in our description of the axion field, we re-organize the sum in (34) such that states with nearby energy are combined. Specifically, we partition the particles into sets, indexed by j, containing all those with ω ∈ [ω j , ω j + ∆ω], within which the states are distinguished only by the random phase. Combining the particles within a given energy cell then amounts to a random walk in the complex plane (see Ref. [107]), leaving The end stage of the walk is a new random phase φ j , with the distance traveled dictated by the Rayleigh random variable α j , drawn from p(α) = α e −α 2 /2 , and the density of states at that energy, controlled by p(ω j ). Through its 15 In the non-relativistic case, it is convenient to express the energy and momentum both in terms of the particle velocity, v. In [108], it was shown how the statistics of the axion field can be described in terms of p(v) (often written f (v)), which includes directional information. This approach can be generalized to the relativistic case by describing the field in terms of k and p(k), rather than the energy as we do in the text.
dependence on α j and φ j , a(t) is itself a random variable. Although a = 0, we expect a 2 > 0; indeed, a 2 is an exponentially random variable, with mean (approximating ∆ω as differential) For a non-relativistic axion, 1/ω −1 ω m a , and we recover the familiar dark-matter result, a 2 = ρ a /m 2 a . For a general energy distribution, however, there is no such simplification (to be clear ω/ω = 1). To exemplify this point, consider a p(ω) which is log flat over [ω 1 , ω 2 ]. If ω 1 ∼ ω 2 , we have ω/ω ∼ 1. However, if they are parametrically separated, ω 1 = ω 2 for 1, then ω/ω ∼ ( ln 2 ) −1 1. To complete our discussion of the relativistic axion in the time domain, in Fig. 7 we show three realizations of the axion field, determined directly from (34), for three different p(ω). In the left two figures we take m a ω, in order to depict examples of the CaB, which are then contrasted with the expected dark-matter axion on the right. On the left, we take p(ω) to be a positive-definite normal distribution, which can be considered an example of a CaB emerging from parametric-resonance production. We takeω = 10 neV, CMB energy density, ρ a = ρ γ , and further set the distribution to be wide, specifically σ =ω. The fact that a number of frequencies are contributing is visible in the realizations. In the middle, we take an even broader p(ω), corresponding to the CaB as predicted from cosmic-string production. In detail, the distribution is determined by (33) with f a = 10 15 GeV and T d = 10 12 GeV, such that ρ a ∼ ρ γ . Nevertheless, we only draw frequencies in a restricted range of ω ∈ [5 neV, 1 µeV], over which we haveω ∼ 10 neV. The presence of both high and low-frequency contributions in p(ω) can be seen in the realizations. Finally, on the right, we show the conventional dark-matter axion scenario, withω m a = 10 neV, and small variations around this as predicted by the standard halo model. The variations are not visible in the time domain, with the period of the realizations highly regular. The statistical nature of the amplitude discussed above, can be seen.
While we can understand the time dependence of the CaB, its properties are more transparent in the frequency domain. As such, consider the Fourier transform of (35). We imagine making measurements of the axion field at a frequency f = 1/∆t for a total integration time T , thereby collecting a set of N = T /∆t discrete measurements of the field, which we denote by {a n = a(n∆t)}. We then calculate the power spectral density (PSD), which quantifies the power in the field at a given frequency, as Technically ω is a discrete variable, given by 2πk/T , with k = 0, 1, . . . , N − 1 the relevant Fourier mode, although we will often assume a sufficiently long integration time that we can approximate ω as continuous. As in the time domain, the PSD is an exponentially distributed random variable, and therefore specified entirely by its mean, Once more, this result reduces to the correct dark-matter expression in the non-relativistic limit. The energy of the non-relativistic wave is specified by its speed, drawn from a distribution f (v). Changing variables to v ω = 2ω/m a − 2, in the non-relativistic limit, we have This agrees with the dark-matter case in [107], demonstrating that the general expression in (38) contains the non-relativistic limit as a special case. Nevertheless, the scaling in (38) is misleading. While it quantifies the power distribution in fluctuations of the axion, experiments can only measure the induced fluctuations in SM fields derivatively coupled to the axion. Accordingly, it is more appropriate to consider the power in g aSM ∂a, where g aSM is the axion-SM coupling. Taking ∂a ∼ ωa, we can determine the parametrics of the accessible power by approximating p(ω) as a uniform distribution over a range of widthω/Q a . Doing so, the power scales as

B. Rough Sensitivity
We will now use (40) to determine the parametric sensitivity of dark-matter experiments to the CaB, leaving a detailed calculation of the sensitivities to the following subsections. We begin with the following simple estimate: assume the CaB can be detected if the power it deposits atω matches the power produced by the darkmatter axion at the sensitivity threshold. For the moment, we assume that experiments extract power from a relativistic and non-relativistic axion wave identically, although we will later justify this assumption up to O(1) factors. In order to compute the power matching, we assume optimistically that the coupling saturates the existing bounds, g aSM = g SE aSM , so that for a fixed Q a andω we can constrain ρ a . For the dark-matter power, we need the dark-matter equivalent of (40), which is obtained by setting ρ a = ρ DM , Q a = Q DM a ∼ 10 6 ,ω ∼ m a , and fixing the coupling to an existing sensitivity threshold, denoted by g lim aSM . Equating the powers at the frequencyω = m a , we expect sensitivity to an axion background that constitutes the following fraction of the CMB energy density, This scaling is overly pessimistic. The CaB will deposit its power over a much wider range than dark matter, so there is more information than can be gleaned by comparing power at a single frequency. In principle, our sensitivity depends on how this additional information is obtained, either through a broadband or resonant readout strategy, and so we will consider the two cases separately. As we do so, however, we emphasize a fundamental challenge: the broad nature of the signal will make it harder to distinguish from backgrounds. There are handles, for instance as we will show for the axionphoton coupling, the signal power will continue to scale quadratically with the magnetic field, and further for the case of the CaB from dark-matter decay, there can be a unique daily-modulation signal. Beyond such remarks, we will not attempt to determine the optimal analysis for a relativistic signal here, although we note it will likely require a more accurate characterization of the background than in dark-matter searches. Indeed, both ADMX and HAYSTAC usually remove features much broader than that expected of dark matter, see for example [35,109], which raises the possibility that a signal of the CaB may already be hiding in existing data, albeit in the most optimistic scenarios.
For a broadband readout of the axion power, we integrate over a range of energies, and therefore at the level of the integrated signal power the distribution p(ω) would seem irrelevant. Yet even when the entire spectrum is resolved, the width of p(ω) still determines an important physical property of the axion field: the coherence time. The coherence time has a straightforward interpretation in the frequency domain. Recall that the measurement time T determines the frequency resolution of the associated discrete Fourier transform, according to ∆ω = 2π/T . For sufficiently small T , the entire signal will fit within a single bin, and the signal amplitude will be associated with one draw from the exponential distribution as outlined in Sec. III A. As T is increased, eventually the resolution will be sufficient to resolve the structure in p(ω). At this stage the signal will occupy multiple bins, each of which will have an independent exponential draw that then combine incoherently. The transition between these two cases defines the coherence time, which we can quantify by 2π/τ = σ ω , with σ ω the width of p(ω). Parametrically, we expect σ ω ∼ω/Q a , so that τ ∼ 2πQ a /ω, or numerically, Consequently, the coherence time of the CaB will generally be short on the timescale of experimental measurements, 16 implying that we will operate in the regime T > τ , where we expect the sensitivity to the signal power to be impeded by the increased background that enters when the signal is distributed over a broader range.
As the CaB has a coherence time that is smaller than for dark matter by a factor of Q DM a /Q CaB a 1, the background will be enhanced by this same scale, and subsequently there is a reduction in signal power sensitivity of Q DM a /Q CaB a . This leads to a refined estimate for the broadband sensitivity of Turning to a resonant detection strategy, the estimate in (41) will be modified by the experimental quality factor Q, associated with the cavity or readout circuit of the instrument, in three ways. Firstly, the power recorded by the resonator is controlled by min(Q, Q a ), so that for Q < Q DM a , we have overestimated the deposited darkmatter power. For the moment, we will assume that we are in this limit, for instance ADMX and HAYSTAC currently operate with Q ∼ 10 5 and Q ∼ 10 4 , respectively. We also expect that Q CaB a Q, so that for equal couplings and density, we expect the CaB power to be suppressed by a factor of Q/Q CaB a , rather than Q DM a /Q CaB a as assumed in (41). The experimental quality factor will enter a second time in defining the instrumental bandwidth of ω 0 /Q, where ω 0 is the resonant frequency. The bandwidth conventionally dictates the range over which the signal can be analyzed. For dark matter, the signal is narrower than the bandwidth by a factor of Q DM a /Q. This implies that when searching for dark matter, the background can be restricted to a smaller range, suppressing its contribution by Q DM a /Q. For the CaB there is no such suppression -the signal extends over the full bandwidth -so using the Dicke radiometer equation [110], our sensitivity will suffer due to the increased background by a further factor of Q DM a /Q, similar to the broadband consideration. The third consideration is in the CaB's 16 Throughout we will always assume that we are operating in the T > τ regime so that the energy distributions can be resolved.
favor. Its broad nature implies that the signal will deposit power over many bandwidths collected during a dark-matter search. At most the number of bins can be Q/Q CaB a , producing an enhancement in the sensitivity of Q/Q CaB a . 17 Taken together, these three factors modify (41) to As Q has dropped out, we are left with the same parametric scaling as for broadband detection.
Given an experimental limit on or sensitivity for dark matter, we can estimate our sensitivity using either (43) or (44). For several specific instruments, we have already shown the results in Fig. 1. We can also consider the expected reach more generally. Taking Q DM a = 10 6 and Q CaB a = 1, all that remains is to fix the couplings. Focussing on the axion-photon coupling, we have g aSM = g aγγ ∼ α/(2πf a ). For the CaB, we take the optimistic value of g SE aSM = g SE aγγ = 0.66 × 10 −10 GeV −1 , whereas for dark matter, we exploit the fact that a broad goal of the axion dark-matter program is to probe g aγγ values at the scale of the QCD axion, which satisfies m a f a m π f π . Assuming this goal is achieved across a wide range of masses, then the corresponding sensitivity to a relativistic population is given by so that for energies below 1 µeV, we could be sensitive to a CaB with an energy density below that of the CMB.
In what follows we will refine this estimate.

C. Relativistic Axion E&M
We now specialize our discussion to detecting the CaB through a coupling to electromagnetism. As is well known, the coupling in (1) leads to the following classical equations of motion [23] ∇ · E = ρ − g aγγ B · ∇a , with ρ and J the charge and current densities, respectively. For non-relativistic dark-matter axions the momentum is parametrically smaller than the energy (∇a ∼ 17 N additional measurements can be thought of as scaling the experimental measurement time T → N T . Assuming T > τ , our sensitivity to the power will scale as √ T → √ N T [111]. k ω ∼ ∂ t a) and this justifies neglecting the two terms involving ∇a. This leaves a single modification to the Ampére-Maxwell equation, which by analogy enters as an effective current J eff = g aγγ B ∂ t a. The axion field thereby converts magnetic field lines into oscillating currents, identifying large magnetic fields as a central ingredient in the detection of axion dark matter.
For the CaB spatial gradients cannot be neglected, yielding two additional sources in (46). The first of these is the generation of an additional effective current, J eff = −g aγγ E × ∇a. As our focus is on searching for the CaB with existing axion dark-matter detectors, which rely on large magnetic fields, this term will not be relevant. The second gradient term, which provides a contribution to Gauss' law cannot be immediately discarded. In detail, a relativistic axion field generates an effective charge density ρ eff = −g aγγ B · ∇a; again by analogy, the axion converts magnetic field lines into oscillating lines of charge. This effect is proportional to B·∇a ∼ B·ka, and is therefore dependent on the incident direction of the axion relative to the experimentally established magnetic field. Nevertheless, we will see that for all the experiments we consider the effective charge does not significantly contribute to the signal. Yet the incident direction of the axion remains detectable: the relativistic field can undergo appreciable spatial oscillations over the instrument, leading to an interference pattern that depends on the incoming angle of the axion wave. We will show explicitly how this effect arises for resonant cavity instruments. For a true cosmological relic, the signal, like the CMB, will be almost completely isotropic (up to ∼ 10 −3 variations associated with our peculiar velocity with respect to the Hubble flow), resulting in an effectively time independent signal. In the scenario where the CaB arises from dark-matter decay, the galactic component will be far from isotropic, instead pointing preferentially towards the Galactic Center. Given that decaying dark-matter will emerge as a case that can be probed already by existing datasets (as it can generate ρ a ρ γ ), this modulation will be an important fingerprint of a genuine CaB signal.
The final equation in (46) allows for backreaction of the electromagnetic fields on the axion itself. To determine when this is effective, consider a particularly simple experimental configuration with B · k = |E| = 0, but with a DC magnetic field of strength B 0 . From Ampére-Maxwell, the axion will induce an AC electric field oscillating parallel to the magnetic field, and with amplitude E ∼ g aγγ B 0 a, generating a backreaction of g aγγ E · B ∼ g 2 aγγ B 2 0 a. For a relativistic field we take ( + m 2 a )a a, and so the condition for backreaction to be irrelevant is parametrically g 2 aγγ B 2 0 /ω 2 1, with The size of the backreaction is negligible for the parameter space considered in this work, though might be of phe-nomenological interest for experiments looking for significantly lower frequency axions. More generally, effects subleading in g aγγ can be neglected. This motivates studying the fields in (46) in powers of g aγγ (see also [112]), Fields carrying a subscript 0 are the dominant fields generated by the experiment, for example a large static magnetic field in ADMX or HAYSTAC, whereas a subscript a denotes axion-induced effects, which are O(g aγγ ). To simplify the discussion, we will assume the large fields are DC (i.e. static), as is the case for many axion darkmatter proposals, although not all, see e.g. [51,52,55]. Under this assumption, the equations for the DC fields reduce to those of electro-and magneto-statics, so that all the physics of interest is contained in the equations for the axion-induced AC fields, 18 We can separate the equations as follows, In the following subsections we will solve the equations in (50) for different experimental configurations. Before doing so, we consider the equations parametrically. Firstly, in the non-relativistic limit, we can drop all axion gradients, and the equations reduce significantly, There are two relevant spatial scales in the problem. The first is the experimental size L, which through the boundary conditions dictates the scale over which the primary fields vary. The second is the Compton wavelength of the axion field, λ a ∼ 1/ω, or in the non-relativistic case λ a ∼ 1/m a . 19 Accordingly, in order to understand the relevance of various terms in (51) qualitatively we can substitute ∇ → 1/L and ∂ t → 1/λ a , and then determine the relevance of each term for specific experiments. A more careful discussion of these scalings is provided in Ref. [112].
To begin with, resonant cavity instruments are designed with a principle that λ a ∼ L, and therefore all terms are relevant. Experiments searching for lighter dark matter, such as ABRACADABRA or DMRadio, have λ a L, suppressing time derivatives with respect to spatial gradients. In particular, (51) then implies that E a ∼ (L/λ a )B a , so that the induced electric fields are parametrically suppressed with respect to the magnetic fields, a point that has been widely discussed [112][113][114][115][116]. In the high mass regime considered by, for instance MAD-MAX [48,117], where λ a L, we instead neglect the spatial gradients. Accordingly, B a ∼ (λ a /L)E a , so that the dominant effect is now the induced electric fields.
A similar analysis can be performed in the relativistic case. We first reinstate the terms containing gradients of the axion field, and then replace those gradients with their parametric scaling of 1/λ a . For simplicity, we consider DC field configurations that are purely magnetic. Then, (50) reduces to Compared to (51), we see for B a there is a single additional term that depends on the relative angle between B 0 and k. For a resonant cavity, once more all terms are in principle relevant. In the low-frequency limit (λ a L), we have We now see that E a ∼ B a , so that the induced electric field is no longer parametrically suppressed. Nevertheless, the origin of the two effects is different. The AC magnetic field is generated by J eff , whereas the AC electric field originates from ρ eff . An identical analysis in the high-frequency regime (λ a L), results in 19 Even in the non-relativistic case, the relevant spatial scale is the Compton wavelength, and not the distance over which the phase of axion field itself varies, which is set by the coherence length. The rationale is that the axion field will drive oscillations in the electromagnetic fields, which having a lightlike dispersion will vary over a spatial scale set by the timescale of their oscillations. and again, neither field is suppressed. Going forward, we will specialize to two specific scenarios from which we can largely infer how the CaB could appear in experiments designed to search for dark matter. In particular, we will consider broadband and resonant detection for λ a L, and resonant detection in the cavity regime λ ∼ L. We will not consider the regime where λ a L, relevant for experiments such as MADMAX. For a CaB, such experiments will not be able to reach axion energy densities relevant for cosmic sources, given the scaling in (45) (see also Fig. 1), though they may have promise in looking for dark-matter decay. Nonetheless, this is the parameter range relevant for the thermal CaB, and therefore it may be interesting to consider dedicated experiments searching for such a background. We will not pursue this direction here.

D. Low-Frequency Detection (λa L)
Armed with the expressions for the induced electric and magnetic fields, we now compute the CaB sensitivity of axion dark-matter instruments focusing on the frequencies well below a µeV (λ a 1 m). For the CaB, our estimated sensitivity in (45) suggests that such experiments are ideal for probing cosmic relics, for which measurements of ∆N eff bound ρ a < ρ γ . While existing instruments only have sensitivity for dark-matter axions with a coupling comparable to the star-emission bounds [43][44][45]118], these results are paving the way for future experiments that will probe the couplings predicted for the QCD axion, such as DMRadio [46,47]. In this mass range, both broadband and resonant search strategies have been proposed. As such, in this section we will consider both types of detection, and to be concrete envision a large scale realization of the DMRadio (or equivalently ABRACADABRA) instrument. 20 Our starting point is the geometry of DMRadio and ABRACADABRA: a large toroidal magnet with field strength B 0 . To understand the effects generated by the CaB on such an instrument, we can use the equations of axion electrodynamics in the λ a L limit, as stated in (53). From the second equation, we see that the axion field will convert the DC magnetic field into an oscillating toroidal current, which will then induce an AC magnetic field in the center of the torus. In the presence of an axion field, a pickup loop placed in the center of the torus would see a varying magnetic flux in a region where conventionally there should be none. This detection principle is identical to the conventional strategy for detecting axion dark matter with such an instrument; that the effect would be the same is clear from the fact the equation for B a in (53) does not involve any gradients of the axion field. Even though the CaB modes have a significantly smaller de Broglie wavelength than dark matter, there is no issue of this leading to an incoherent effect across the instrument, as that would only occur for λ a L, outside the range considered by these instruments.
There are, however, differences for the CaB detection from the conventional dark-matter axion search. Firstly, the range of frequencies over which the AC magnetic field will be excited are significantly larger than for dark matter, as we have emphasized many times already. A second difference is that unlike in the non-relativistic case, there is now an unsuppressed electric field generated from the effective charge. Recalling that the relativistic axion converts magnetic field lines into oscillating charge lines, the instrument would behave like a torus of oscillating charge, inducing an axial AC electric field near the center of the torus. Supposing the pickup loop used to search for B a is perfectly perpendicular to this field, the above detection scheme is unaffected. Nevertheless, the oscillating electric field is present, and its detection could provide a confirmation of any magnetic field excess.
Turning to the actual detector response, integrating the effective current over the torus, the flux induced in the pickup loop will be [41] where B 0 is the magnetic field at the inner radius of the torus and V B is the magnetic field volume. In ABRA-CADABRA this flux is read out by inductively coupling the pickup loop to a SQUID, which will observe a flux Φ a proportional to Φ pickup , with a proportionality constant β. For a 100 m 3 instrument, β 0.5% [41]. As discussed in Sec. III A, this signal can be searched for by collecting a time series data set of the SQUID flux, taking the discrete Fourier transform of the measurements, and finally forming the PSD, S Φ (ω). Going through these steps and recalling that the PSD of the CaB is exponentially distributed, with mean given in (38), the result is that the power will be an exponentially distributed quantity, with mean where we have introduced the time derivative of (38), S ∂ta = ω 2 S a . Here λ B (ω) is the contribution to the flux from background sources, which in general will be frequency dependent. For a broadband strategy, λ B is limited by noise within the SQUID, numerically given by λ B 1.6 × 10 5 eV −1 . The steps leading to this result parallel closely those in the non-relativistic case considered in [107], and we refer there for additional details.
In particular, recalling that in the non-relativistic limit p(ω) = f (v ω )/(m a v ω ), (56) contains the dark-matter result as a special case.
Importantly, the exponential nature of the PSD implies that we can exploit the full likelihood framework of [107].
We can then analytically determine the expected sensitivity to our signal through the use of the Asimov data set [121], where instead of considering the distribution of our sensitivity on a set of simulated data, we instead replace the data with its asymptotic expectation. Doing so, for a given g aγγ , our sensitivity to the CaB density is given by Here T is the data collection time, and TS is the test statistic associated with the sensitivity threshold (for 95% expected limits TS 2.71, whereas for an n-σ discovery, TS = n 2 up to the look elsewhere effect). To arrive at this result we assumed that T τ , such that p(ω) is well resolved, and also that we are in the background dominated regime. Further, we have left the terminals on the frequency integration unspecified, but these can at largest be the range over which the experiment has sensitivity.
The result in (57) is consistent with the simple estimate claimed in Sec. III B. To see this, we model the broad CaB as a log-uniform distribution p(ω) = Q CaB a /ω defined over a rangeω/Q a . Then, treating λ B as independent of frequency, our sensitivity can be written (dropping ρ γ ) For dark matter, we instead approximate the distribution as a uniform p(ω) = Q DM a /m a over a narrow range, which results in an identical expression but with Q CaB a → Q DM a andω → m a . Taking m a =ω, the ratio of the two expressions yields (44).
The sensitivity in (57) represents the quantitative result for broadband sensitivity, however to provide additional intuition -especially for why the width of p(ω) plays a fundamental role -we can rederive the result qualitatively from a signal-to-noise ratio [41,111]. The signal strength is simply the average signal flux, which is given by (using (36)) whereas the flux noise is given by √ λ B . Given the signal and background flux magnitudes, if we perform a measurement for a time T , the naive expectation is our sensitivity will grow as √ T , and indeed for a time it will. Nevertheless, as described already, the CaB has a finite coherence time τ ∼ 2πQ a /ω, and for T > τ the signal will no longer combine coherently, leading to the signal-to-noise only scaling with the parametrically reduced (T τ ) 1/4 [111]. Assuming T > τ , then our sensitivity is given by |Φ a |(T τ ) 1/4 / √ λ B = 1. Together, these scalings provide which is parametrically identical to (58), and demonstrates that the appearance of the width of p(ω) is a consequence of measuring the CaB over times longer than the field is coherent. The alternative qualitatively different readout strategy proposed in this frequency range is resonant detection. For the toroidal geometry described above, this can be achieved by reading out the pickup loop through a resonant circuit, which for our purposes can be characterized by four parameters: the quality factor Q, resonant frequency ω 0 , total circuit inductance L T , and thermal noise temperature T 0 . As in the broadband case, we can generalize the known dark-matter result (see e.g. Ref. [107]) to the CaB as follows, which, up to experimental factors, is identical to (57). This expression determines the expected CaB energy density sensitivity for a given set of experimental parameters. Yet we can recast this result in the spirit of Sec. III B, and forecast CaB sensitivity in terms of the expected dark-matter reach. Indeed (57) holds equally well for dark matter, taking ρ a = ρ DM , and evaluating where in the final step we assumed f (v) follows the canonical standard halo model. For a resonant instrument, the frequency range is commonly restricted to a single bandwidth of size ∆ω = ω 0 /Q around ω 0 . For Q < Q DM a this detail is irrelevant in the dark-matter computation, as for m a = ω 0 , then the dark-matter distribution is contained entirely within the bandwidth. For the CaB however, we will generically have Q CaB a Q, and therefore expect p(ω) to be constant over the range of integration, leading to Having computed the result for both dark matter and relativistic axions, we can take the ratio to determine the CaB sensitivity as a function of the experimentally achieved dark-matter coupling, g lim aγγ , to be .
If we approximate the energy distribution as log flat, so p(m a ) = Q a /m a , takeω = m a , and recall that the CaB sensitivity can in principle be enhanced across at most N ∼ Q/Q CaB a bandwidths, then this result reproduces the parametric scaling given in (44).

E. Resonant Cavity Detection (λa ∼ L)
We next consider detection with a physical resonator, as pursued by both ADMX and HAYSTAC, where a microwave cavity is constructed in order to resonantly enhance power produced at a frequency tuned to the cavity dimension, ω 0 ∼ 1/L. The design principle for these instruments is to optimize the search for the narrow spectral feature dark matter predicts when its mass is such that 1/m a ∼ L ∼ m. In particular, in the presence of a large static magnetic field B 0 , an axion background will source oscillating electromagnetic fields, thereby generating potentially detectable power in the cavity. As we will show in this section, this statement is true for both dark matter and the CaB. Importantly, the existing reach of ADMX is already sufficient to probe open parameter space, although a reanalysis of the data would be required, and the same will soon be true of HAYSTAC.
Parametrically our final sensitivity will be identical to the resonant circuit expression (64). There will, however, be important differences. According to the expressions in (52), the CaB will source AC electric and magnetic fields in the presence of a large B 0 , and if we work in a regime where the DC field is spatially uniform, we have whilst B a can be determined from Faraday's law. This result is identical to the case of dark matter: we have dropped the contribution from the effective charge as it cannot excite resonant cavity modes. 21 Unlike for a non-relativistic axion, when λ a ∼ L the axion wave can undergo O(1) spatial oscillations across the detector, which will induce a new daily modulation effect we will explore. Taking a simplified picture of the CaB where a(t, x) ∝ cos(ωt − k · x), we will find an explicit dependence on the direction of k. In particular, when we compute the cavity form factor -which captures the overlap between the axion source and the relevant mode of the cavity being measured -a dependence on the direction of k and hence the CaB will appear. For a cosmic relic, the CaB will be incident approximately isotropically on the detector, and the effect will take on a sky-averaged, and effectively time-independent, value. 22 However, for local sources of relativistic axions, the distribution can be anisotropic. This will certainly be the case for darkmatter decays in the Milky Way, where the flux predominantly originates from the Galactic Center. In the rest frame of the Earth, where the orientation of the cavity is time invariant, the direction of the Galactic Center will vary throughout the day, leading to a daily variation in the cavity form factor. Accordingly, in general the CaB power, which is proportional to the form factor, will undergo O(1) variations throughout the day. These daily modulations provide a novel handle that can be used to distinguish a CaB with a local origin from potential backgrounds, although in the present work we will not quantitatively calculate this effect. Our focus is instead on determining the CaB sensitivity of these instruments, which will require a determination of the power a relativistic axion field deposits in a cylindrical cavity as used by both ADMX and HAYSTAC. We do so by repeating the analogous non-relativistic computation [32,33] (see [122] for a recent review) while accounting for the three important modifications in the relativistic case. Two of these we have already discussed: the additional gradient term in (65) and the fact that the CaB carries power over a much broader frequency range. The third relativistic novelty arises from the fact the axion field need not be spatially coherent across the instrument when |k| = ω ∼ 1/L, which can suppress the integrated power deposited over the cavity.
We begin with (65). We will solve this equation within a cylindrical cavity, where a DC magnetic field B 0 = B 0ẑ has been established along the cavity axis, and we will assume the field is spatially uniform. The cavity will have a set of normal modes for the electric fields it can support, e mn (x), indexed by three integers ( , m, n), and satisfying with ω mn the resonant frequency of the mode. On dimensional grounds we expect ω mn ∼ 1/L, however the exact value will be determined by the spatial variation of the modes. As these basis modes must satisfy the electric boundary conditions, O(1) changes to the resonant frequency can be achieved by modifying these boundary conditions, as ADMX achieves by mechanically varying the location of tuning rods within the instrument. We can normalize the basis modes as follows, The integral is performed over the cavity volume V , which reveals that the modes carry dimension e ∼ V −1/2 . As the orthonormal modes are a complete basis for the electric field in the cavity, we can write E a = α mn e mn and solve for the coefficients. Using this expansion in (65), and transforming to the frequency domain, we obtain ,m,n (ω 2 − ω 2 mn )α mn e mn = − ω 2 B 0 g aγγ a(ω)cos(k · x)ẑ , where we have assumed a simple form for the axions spatial dependence, a(ω, x) = a(ω) cos(k · x). Using the orthonormality condition, we can then isolate the electric field coefficients as, The first line of this result identifies ω mn as the resonant frequencies, with the apparent divergence a remnant of our idealized treatment of the problem. In particular we have neglected the fact that the electric field will penetrate the walls of the cavity and dissipate energy due to the finite resistance of the material, which is the origin of the cavity quality factor, Q. This dissipation can be accounted for as in the standard dark-matter calculation. More novel is the second line of (69), and we define which we will use to define a generalized cavity form factor shortly, and note dimensionally κ mn ∼ V 1/2 . The energy density in the AC cavity fields for this mode is given by U mn = (|E a | 2 + |B a | 2 )/2 |E a | 2 , as for frequencies near the resonant frequency we will have |B a | |E a | from Faraday's law. Collecting our expressions above, the energy density has the form Here we have defined a transfer function, T (ω), which accounts for the resonant response of the circuit when dissipation is included, which is sharply peaked around its maximum, T (ω mn ) = 4Q 2 /ω 2 mn . The energy density is further expressed in terms of a cavity form factor C mn = |κ mn | 2 /V , where cavity volume has been introduced to ensure this is an intensive quantity. In detail, we define a relativistic cavity form factor, which can be contrasted with the conventional result used for dark matter Let us discuss several details of these form factors, which parameterize the overlap between the induced E a and the static B 0 . In both cases, numerically we have C < 1. Further, we expect C to be O(1) for the lowest lying mode -higher modes correspond to basis functions of shorter wavelength, which generically suppress the integral. In the non-relativistic case, contributions to e perpendicular toẑ do not contribute, which identifies the cavity transverse magnetic (TM) modes as relevant.
We have yet to specify e, however for a cylindrical cavity they are well known (and our results are qualitatively similar for other geometries). Taking the cylinder to have height L, and radius R ∼ L, we have where J n are Bessel functions of the first kind, and the resonant frequency is given by ω 00 = j 0 /R, with j 0 the th zero of J 0 . Explicitly evaluating the cavity form factor, we obtain Numerically, C 100 0.69, and C 00 ∼ C 100 / 2 , so that the response of higher order modes is rapidly suppressed. Accordingly, it is common to focus on the lowest lying mode, defining ω 0 = ω 100 and C = C 100 .
The relativistic cavity factor is complicated by a dependence on the incident angle of the axion throughn (where k = ωn) and a frequency dependence through k. With a view to searching for the CaB through a repurposed dark-matter search, we will only consider the response induced for the lowest TM mode, although we note from (73) that in this case transverse electric modes could also contribute. Combining (73) and (75) for = 0, where we employ a shorthand s α = sin α and c α = cos α, with α = arccos(n·ẑ). The additional relativistic novelty is K(ω, α), which accounts for the incoherence of the axion field over the experimental volume, and is shown in Fig. 8. As shown there, for ω ω 0 we have K(ω, α) → 1, corresponding to the limit where the axion is spatially coherent over the instrument, whereas for ω ω 0 instead K(ω, α) → 0 as a result of destructive interference across the cavity, and for L ∼ R the result is only weakly dependent on α. 23 This factor effectively removes the contribution of frequency with ω > ω 0 , and we will approximate it by a step function, K(ω, α) Θ(ω 0 − ω).
Having determined the modified form factor, we can return to determining the measurable signal power, which  (77). We fix the cavity radius to be equal to the height, R = L, and demonstrate the approximate insensitivity of the result to the relative angle of the incident axions and the detector magnetic field, α = arccos(n ·ẑ). For R = L, a large dependence on α can arise, generating a potential daily modulation in the CaB signal. Note that this factor enters in the relativistic case in addition to the transfer function in (72), which is sharply peaked at ω0.
is related to the cavity energy density by P = ω 0 U/Q. Using (71), we have The form of |a(ω)| 2 has already been extensively discussed, it will be an exponentially distributed variable with mean given in (38). Accordingly, the power will also be exponentially distributed, however the total power on average will be given by where in the final step we assumed that Q CaB a Q, so that p(ω) only varies slowly over the range where the transfer function has appreciable support.
In (79) we have an expression for the signal power the relativistic axion will deposit in the cavity when analyzing a single resonant frequency. Combined with an expected background contribution and a set of experimental parameters, this result is sufficient to forecast the CaB sensitivity. Here, however, we will instead use a matched power approach to obtain the projected reach. In particular, existing ADMX limits are a combination of individual experimental runs, so rather than combining these run-by-run, we will simply recast the combined dark-matter g aγγ limits. To do so, we require the darkmatter analogue of (79), which is obtained by integrating over the full frequency range, and taking a mean |a(ω)| 2 as given in (39). Re-evaluating the integral, we find We cannot directly match CaB and dark-matter powers above, as comparing the signal strength will only determine the sensitivity assuming the background is equal in the two cases. However, as already discussed in Sec. III B, it is not, and the larger background for the CaB will reduce its sensitivity by Q DM a /Q. Once this is accounted for, we can combine (79) and (80) (evaluated at ω 0 = m a ), to obtain our estimated sensitivity in this regime of which, as promised, is identical to (64) up to a numerical prefactor. Several aspects of the above discussion are overly idealized. 24 In particular, we considered a configuration where the scale at which the relativistic power suppression occurs -as encoded in K(ω, α) -matches the resonant frequency scale, both set by the characteristic size of the cavity. In practice, however, microwave cavity instruments introduce tuning rods in order to vary the resonant frequency, which will shift ω 0 by a factor of a few above the characteristic value ∼ 1/L. Yet as the relativistic suppression occurs for ω 1/L, when the axion field is no longer spatially coherent across the cavity, K(ω, α) as it appears in Fig. 8 will remain qualitatively unchanged. Combining the two effects, naively this would suggest a significant reduction in sensitivity as the modes that are resonantly enhanced also experience the incoherent suppression. Nonetheless, in an actual instrument, the geometry is such that L = R, and once that is accounted for K(ω, α) now varies considerably with angle, and for α = 90 • the suppression is postponed till higher frequencies, so that the power can still be absorbed, and with the appearance of a significant daily-modulation effect that can be exploited.

IV. PROJECTED LIMITS
Having outlined various forms the CaB can take, and having determined the experimental sensitivity to it, in this section we combine these results to sketch projected sensitivities. As we will show, detecting the CaB will not necessarily require dedicated instruments, instead the rapid progress in the search for axion dark matter will simultaneously open enormous swaths of relativistic axion parameter space. The present discussion will 24 We thank Jonathan Ouellet for this observation. not be exhaustive. Instead we will show the estimated reach in three cases to demonstrate various aspects of the detection schemes we have proposed, and the interplay with specific CaB candidates. Firstly, we will discuss the reach of the existing resonant cavity instruments HAYSTAC and ADMX for a simple Gaussian p(ω) as predicted in the parametric resonance scenario, showing that an order of magnitude improvement in the g aγγ sensitivity of ADMX would translate to sensitivity to ρ a < ρ γ , although this would still be short of the prediction of parametric-resonance production discussed in Sec. II C. Secondly, we demonstrate that a large scale broadband ABRACADABRA style instrument could probe relativistic axions originating from cosmic strings in the parameter space where they could help alleviate the Hubble tension. Finally, we will consider the case of most immediate interest: indirect detection of dark matter decaying to axions. As this scenario allows ρ a > ρ γ , we will see that ADMX is already sensitive to unexplored parameter space, and the situation will improve dramatically with future instruments.
Recall that the CaB signal is determined by three quantities: ρ a , g aγγ , and p(ω) (equivalently, the form of Ω a (ω) fixes ρ a and p(ω)). As already discussed, here we will take the approach of fixing g aγγ to g SE aγγ , the largest value consistent with star-emission constraints. When considering detection with frequencies ω 10 −9 eV, we take g SE aγγ = 0.66 × 10 −10 GeV −1 for consistency with CAST [25] and Horizontal Branch [26,27] constraints. In the future, this limit may be tightened by IAXO [123] (see also Refs. [124,125] for future searches using the CMB). Should these experiments detect an axion signal, the same axion could be produced in the early Universe, and would strongly motivate further searches for the CaB. At lower frequencies, the bounds strengthen further, and we will adopt g SE aγγ = 3.6 × 10 −12 GeV −1 as determined from super star clusters [31].

A. Gaussian
To begin with, we consider searching for a Gaussian energy distribution in existing resonant cavity instruments. Such a distribution can be motivated by the parametricresonant production mechanism reviewed in Sec. II C, however here we can also envision it as providing an opportunity to explore our formalism. In detail, we consider a positive definite Gaussian with mean µ =ω, and variable width σ = κω, where we will explore several values of κ. In detail, we take Fixing g aγγ , we can then determine a limit on ρ a for a given κ.
To do so, we will recast existing bounds on axion dark matter collected by the ADMX and HAYSTAC instruments. ADMX is already probing the couplings predicted Sensitivity for a Gaussian CaB, as predicted by parametric-resonance scenario, that can be obtained by searching for a relativistic axion in data already collected by ADMX and HAYSTAC in the search for axion dark matter. In particular we take p(ω) as in (82), and consider three values of the width κ. If the CaB is a cosmological relic, then it must have ρa < ργ, which ADMX would reach with an order of magnitude improvement in its gaγγ reach.
for the QCD axion for m a ∼ 2 − 3 µeV [34][35][36]. Given (45), we would therefore expect the instrument to be on the verge of ρ a ∼ ρ γ sensitivities for the CaB. HAYSTAC, on the other hand, is already within a factor of a few from the QCD prediction for m a ∼ 23 − 24 µeV [37]. Here we take these existing limits and recast them using (81).
Our forecast sensitivity is provided in Fig. 9. In determining the plotted sensitivities, we combined the single bandwidth sensitivities in (81) across multiple bins, accounting for the spread of p(ω). In doing so we assumed the frequency range scanned by the instruments was divided into bins of width ω 0 /Q, taking Q = 10 5 for ADMX and 10 4 for HAYSTAC. In detail, we compute a limit ρ a,i in each bin, indexed by i, and then determine the combined limit as ρ −2 a = ρ −2 a,i . Note that in the event of the limit being identical across N bins, this returns the expected ρ a = ρ a,i / √ N . As the figure demonstrates, at present neither instrument is sensitive to the ρ a < ρ γ required for a cosmological relic. Nonetheless, ADMX is within two orders of magnitude of the relevant parameter space, which from (81), would be achieved with an order of magnitude improvement in their sensitivity to the dark matter g aγγ . A factor of 30 improvement would allow ADMX to access the parameter space required to reduce the H 0 tension. The situation is more challenging for HAYSTAC, which would require at least three orders of magnitude improvement in their coupling sensitivity to reach ρ a ∼ ρ γ , although widening the range of masses considered by a factor α would also enhance there sensitivity by √ α. Note that in Fig. 9, the peak sensitivity depends on the width of the distribution, κ, and is not always located within the range of masses directly probed by ADMX and HAYSTAC for κ > 1. To understand this, recall from (40) that the signal power is determined by ρ a /ω = n a , the number density. For a fixed ρ a , we can increase n a by decreasingω, and still obtain a constraint as long as p(ω) has support.

B. Cosmic Strings
Next we consider sensitivity to a cosmic-string origin of the CaB as discussed in Sec. II D. In this work we will use the specific results provided in Refs. [98,99], as already discussed, however we note that further improvements in the predicted string spectrum will impact the sensitivities we present. Regardless, as demonstrated in Fig. 5, the cosmic-string spectrum is expected to be especially broad. As such, we will use it as an example to forecast sensitivity with a futuristic broadband instrument operating in the low-frequency regime, defined by λ a L. The spectrum and energy density of cosmic-string axions is determined determined by the symmetry breaking scale, f a , and subsequent temperature at which the string network enters the scaling regime, T d ≤ f a . In terms of their impact on the predicted CaB spectrum, T d < f a provides an effective cutoff on the spectrum at higher frequencies, whereas the energy density in the spectrum is controlled by ∼ f 2 a . Accordingly, for a given T d , we can construct our sensitivity to f a by determining where a detectable axion power is produced. Recall our broadband sensitivity to ρ a was given in (57) for an ABRA-CADABRA type instrument. Assuming a frequency independent background, we can rearrange that result to obtain dω ω dn a dω For a fixed g SE aγγ , the left hand side determines the signal strength, so that the result determines our sensitivity. The differential energy density can be determined from (33) as For a given T d and f a , this provides all the ingredients for the signal prediction, and what remains is to set the experimental parameters. For this purpose we use the parameters adopted in the most optimistic scenario provided in the original ABRACADABRA proposal [41], which involved a 100 m 3 volume instrument with a 5 T magnetic field operating for a year. In order to determine the expected 95% sensitivity we further take TS = 2.71. Lastly, we need to specify the frequency range of the search, which will enter in the terminals for the signal integral in (83). Here we take 2×10 −13 eV < ω < 10 −7 eV,  where the lower limit is set at 50 Hz, where the 1/f noise is expected to begin dominating, and the upper limit is determined by the physical size of the instrument. Going to such low frequencies requires us to take the enhanced star emission bounds of g SE aγγ = 3.6 × 10 −12 GeV −1 . The end result of the above discussion is the forecast sensitivity shown in Fig. 10. In the same figure we also depict the parameter space where the cosmic-string spectrum over this frequency range would be in tension with ∆N eff bounds, and once more where the model could reduce the Hubble tension. This entire parameter space would be covered by the large scale broadband instrument. The sensitivity is dominated by the contribution at the low-frequency end of the instrument, and the flattening of the sensitivity at T d ∼ 10 9 GeV arises when the decoupling induced cutoff in the spectrum occurs within the experimental frequency range.

C. Dark-Matter Indirect Detection
In the examples thus far, the relevant CaB parameter space will only be reached in future instruments. This in a consequence of bounds from observational cosmology limiting ρ a ρ γ . However, if the CaB is produced in the late Universe, ρ a > ρ γ is allowed, and in principle may already be detectable. Axions produced from darkmatter decay, discussed in Sec. II B, is one such scenario. Searching for these axions would open up a new channel in the broader search for dark-matter indirect detection.
Recall that relativistic axions produced from darkmatter decay will receive a contribution from both decays in the local Milky Way halo and extragalactic dark matter. The latter will arrive approximately isotropically at the Earth, whereas the former will preferentially originate from the Galactic Center given its higher darkmatter density. For detection at resonant cavity instruments, we reiterate that the local decays will be associated with an observable daily modulation, as the signal is proportional to sin 4 α, with α the relative angle of the incident axions and the cavity magnetic field. Effectively this search uses the resonant cavities as dark-matter telescopes, although with peak sensitivity obtained when the instrument is perpendicular to the source. Whilst this effect can be used as a fingerprint of a genuine signal, for the sensitivity estimates we present here we will simply take the sky averaged value.
Doing so, our sensitivity for resonant cavity instruments was provided in (81), which we can re-express as where the signal has been written in terms of the differential number density evaluated at ω = m a . This number density is a combination of the local and extragalactic densities given in (9) and (5), respectively. As discussed, the local distribution will be broadened by the Doppler shifts arising from both the finite velocity distribution of the dark matter, and also the Earth's motion relative to the halo, although here we will simply model the distribution as a Gaussian with relative width 10 −3 . With these distributions, for a given m DM , the flux is dictated by the dark-matter lifetime τ . While a search for indirect detection with axions has not yet been performed, there are constraints on dark matter decaying to a relativistic species from cosmology, which roughly limit the lifetime to be longer than the age of the Universe [75,77,78]. Repurposed axion darkmatter searches will be able to do considerably better. To begin with, by repeating the approach of combining experimental bandwidths as we did for the Gaussian case above, we determine that ADMX is already sensitive to open parameter space, as we show on the left of Fig. 11. The strongest sensitivity is obtained when 2m DM falls within the ADMX search window of 2 − 3 µeV, as this corresponds to detecting the local population of axions, which have energies peaking at m DM /2. For higher darkmatter masses, the decay can still be detected thanks to the redshifted extragalactic spectrum. We emphasize once more that this figure shows the sensitivity ADMX can obtain to this scenario using existing data. Even for local decays, the signal is significantly broader than that expected of dark matter, and thus in existing searches the effect would have been discarded as background.
The future prospects for this search are considerable, as we demonstrate on the right of Fig. 11. Performing a naive forecast of ADMX and HAYSTAC by assuming they improve their reach for g aγγ by one and two orders of magnitude, respectively, they will both be able to probe open parameter space. At lower masses, were a future instrument such as DMRadio able to reach the projected sensitivity of the QCD axion g aγγ prediction from 0.5 neV FIG. 11. Sensitivity to relativistic axions resulting from dark-matter decay using existing (left) and future (right) resonant instruments. We emphasize that the forecast sensitivities are likely to be obtained on significantly different timescales for the various instruments. The H0 exclusion region is a result of the decay of dark matter to a relativistic species leading to an observable modification of cosmology, and we show the constraints τ 3.6 tU obtained by the DES collaboration [78].
to 0.8 µeV [46,47], their reach would be considerable. In particular, repurposing the low-frequency resonant result in (64), and assuming DMRadio operates with Q = 10 6 (although their actual readout strategy will be more complex), we see that the instrument will be sensitive to an enormous range of parameters. We note that the future projections shown here may have different timelines and do not necessarily represent a fair comparison between instruments.

V. DISCUSSION
In this work, we considered the possibility of detecting a Cosmic axion Background, an ultrarelativistic background of relic axions. Being naturally light, axions produced in the early Universe are typically relativistic and may exist over a large range of possible energies. Its detection would have imprints of the history of the Universe in its energy spectrum, in close analogy with the program searching for a stochastic gravitational wave background. In particular, we discuss sources sensitive to the reheat temperature of the Universe (thermal production), the origin of dark matter (dark-matter decay), inflation (parametric resonance), and early Universe phase transitions (cosmic-string emission). Most of these production mechanisms predict axions with energies in the range detectable by current and future axion dark-matter experiments. The exception is the well-motivated thermal CaB, where new ideas will be required to probe the predicted energies and densities. For relativistic axion production before recombination, probes of the expansion of the Universe limit the axion energy density to be below that of the CMB. Moreover, recent measurements have found a persistent discrepancy between the early and late measurements. An additional source of radia-tion is the simplest solution to partially alleviate the tension, motivating a value for the CaB energy density and we discuss the production of axions in light of this target. If an axion is discovered from star emission searches such as hinted by the recent Xenon-1T excess [126] or by future experiments such as IAXO [123], there would additionally be a clear target coupling for the axion, greatly increasing the urgency of conducting CaB searches.
For detection, we have focused on axions coupled to photons. Relativistic axion relics contribute terms to Maxwell's equations which are negligible for axion dark matter. The influence of the new terms depends on the experimental design and we study both resonant cavity experiments (e.g., ADMX, HAYSTAC) and lumpedcircuit readout instruments (e.g., ABRACADABRA and DMRadio). For resonant cavities, the power deposited in the cavity becomes sensitive to the axion propagation direction and will exhibit a daily modulation for nonisotropic sources (such as dark-matter decay). For oscillating magnetic field searches, we show that the power deposited is independent of the incoming axion direction, however oscillating electric fields develop within the experiment that might be useful in confirming a positive signal. In either case, we develop simple relations to estimate the prospective sensitivity of axion experiments to a CaB. While present searches for dark-matter axions can have sensitivity to a CaB, this typically requires a dedicated search due to the presence of a broad energy spectra. In particular, with a search for a relatively broad axion energy distribution, current data from the ADMX experiment may be able to discover dark matter decaying into axions, thereby turning the instrument into an indirect-detection axion telescope. In the future, axion detection experiments will be able to discover ambient axion densities well below that of the CMB and probe a wide range of motivated sources.
Throughout this work, we have largely neglected the influence of a prospective axion mass, m a . If, for a given axion source, m a becomes comparable to the axion energy as the Universe expands, then axions would cluster in galaxies. This would increase the local density, analogous to the CνB. If axions become highly non-relativistic they make up a form of dark matter as currently searched for by axion haloscopes. Alternatively, if the axion is still mildly relativistic, it can have a relatively wide energy distribution and still be overlooked by present analyses. The detection prospects of axions with intermediate masses is an interesting question we leave for future work. Experiments could also look for CaB with other interaction terms, such as the axion-nucleon coupling. Since the projections of experiments such as CASPEr [111] have sensitives well below star-emission bounds, they will also be able to probe cosmological sources of relativistic axions. In addition, one could consider other relativistic bosons such as dark photons. We leave the study of their prospective sources and detection capabilities for future work.
The axion-photon coupling has a natural value related to its decay constant. If the interactions arises from integrating out charged fermions, f , of charge Q f , then the coupling generated is, In the absence of parametrically small charges, the expectation is that g aγγ ∼ α/2πf a . However, it is conceptually simple to envision scenarios where the coupling is well below this scale either through axion-axion [127] or photondark photon [128] kinetic mixing (see also Refs. [105,106] for related discussions). To see this we first consider a kinetic mixing between an axion, a, and a second axion, b. We take the second axion to have a decay constant f b that couples to electromagnetism, whereas a does not couple to any charged fields, so that it has no direct photon coupling. In detail, we take the following Lagrangian, To diagonalize the axion kinetic terms to leading order in ε, we send a → a, b → b − εa. Once axion masses are introduced, this transformation leaves the rest of the Lagrangian approximately unchanged only if m a m b . With this shift the Lagrangian for a can now be written as, Importantly, the coupling between a and electromagnetism is determined by the free parameters ε and f b -it can be small even if f a is well below the weak scale. While axion-axion mixing breaks the relationship between g aγγ and f a , it also requires an additional light axion that may influence the phenomenology. Alternatively, one can consider a kinetic mixing between the photon and a dark photon. Consider the Lagrangian, where α is the dark gauge coupling and the photon-dark photon mixing is set by . If the dark photon, A , has a mass below the photon plasma mass, then the photondark photon mixing can be eliminated with the transformation A → A , A → A − A . This leaves the dark photon approximately massless, but generates an axionphoton coupling, The coefficient can be arbitrarily small even for f a below the electroweak scale, breaking its conventional relationship to g aγγ . Interestingly, there may also be an opportunity to discover additional light states in this scenario.
If f a 1 TeV then in the limit that the dark-photon mass is negligible, there are millicharged particles in the spectrum. These particles can have a mass of, at most, ∼ 4πf a . For f a 100 eV solar cooling bounds on millicharged particles would require 10 −13 [129] and may place a bound on the parameter space.