Dark Matter Targets for Axion-like Particle Searches

Many existing and proposed experiments targeting QCD axion dark matter (DM) can also search for a broad class of axion-like particles (ALPs). We analyze the experimental sensitivities to electromagnetically-coupled ALP DM in different cosmological scenarios with the relic abundance set by the misalignment mechanism. We obtain benchmark DM targets for the standard thermal cosmology, a pre-nucleosynthesis period of early matter domination, and a period of kination. These targets are theoretically simple and assume $\mathcal{O}(1)$ misalignment angles, avoiding fine-tuning of the initial conditions. We find that some experiments will have sensitivity to these ALP DM targets before they are sensitive to the QCD axion, and others can potentially reach interesting targets below the QCD band. The ALP DM abundance also depends on the origin of the ALP mass. Temperature-dependent masses that are generated by strong dynamics (as for the QCD axion) correspond to DM candidates with smaller decay constants, resulting in even better detection prospects.


I. INTRODUCTION
The particle nature of dark matter (DM) is unknown. One particularly well-motivated DM candidate is the QCD axion, which also provides a solution to the strong CP problem [1][2][3][4]. The QCD axion is a pseudoscalar boson with an approximate shift symmetry, and its mass and couplings are mostly controlled by a single parameter, the axion decay constant f a . A relic cosmological abundance can be obtained from the misalignment mechanism [5][6][7], and as a result, axion DM has phenomenological properties very different from thermally-produced weakly interacting massive particles.
More general light pseudoscalars are known as axionlike particles, or ALPs. ALPs arise as pseudo-Nambu Goldstone bosons (pNGBs) associated with the breaking of global U (1) symmetries, or as zero modes of higher dimensional gauge fields that are generic in string theory [8][9][10][11]. Unlike the QCD axion, ALPs do not have to interact via the strong force and therefore they are not associated with the strong CP problem. As a result, they exhibit a wider range of couplings and masses and offer a compelling class of DM candidates. For a review of ALP and axion model-building and cosmology see, e.g., [12,13].
Recent years have seen a resurgence of interest in searching for axions and ALPs, with a number of active experiments and new proposals under consideration (for reviews see e.g. [14,15]). These include resonant cavity experiments at various frequencies, such as ADMX [16,17], ORGAN [18] and HAYSTAC [19,20], and also new ideas including dielectric haloscopes (MAD-MAX [21,22] and photonic materials [23]), resonant LCcircuits [24,25], detection-induced magnetic flux oscillations (ABRACADABRA [26,27]) and NMR-based techniques (ARIADNE [28,29] and CASPEr [30][31][32]). Col-lectively these experiments cover many orders of magnitude of possible ALP mass, and are sensitive to ALP couplings to photons or nucleons depending on the experiment. In this work, we focus on cosmological relic populations of electromagnetically-interacting ALPs. The impact of resonant cavity searches on these ALPs has previously been considered in [8].
Recently, the Physics Beyond Colliders Working Group has forecast the sensitivity of future experiments to axions and ALPs [33], building on the review [15]. For the QCD axion, a number of groups have developed models to expand the parameter space, classifying the possibilities for UV-complete theories [34][35][36], model-building photophilic [37] and photophobic [38] axions, and extending the standard misalignment mechanism [39,40]. These analyses highlight the breadth of viable QCD axion models extending beyond the canonical KSVZ and DFSZ scenarios, and motivate continued experimental exploration of the axion mass m a and the axion-photon coupling g aγγ parameter space (the "ALP plane").
Our aim is to map cosmological models onto the ALP plane, identifying regions where the correct relic abundance is obtained from simple assumptions about the expansion history, the ALP model, and the initial conditions. These regions of parameter space are therefore compelling targets for experiments searching for electromagnetically-coupled ALP DM. Since ALPs do not necessarily couple to the strong interactions, and their relic density depends on the expansion rate at early times, these targets can differ significantly from the QCD axion with a standard radiation-dominated cosmological history.
In Section II we consider an ALP with relic density set by the misalignment mechanism. The final abundance strongly depends on the expansion history of the universe before Big Bang Nucleosynthesis (BBN). We study ALPs that begin to oscillate during radiation domination (as in the standard cosmology), during an epoch of early matter domination (EMD), or during a kination phase. ALPs in these alternative cosmologies have also been considered recently in e.g. [41][42][43][44][45], which we build upon in our work. In Section III we study the impact of the origin of the ALP mass on the relic abundance. We determine the parameter space favored by ALP DM with a fixed mass during and after the onset of oscillations, a mass derived from higher-dimensional Planck-scale suppressed operators, and ALPs with a mass that changes with temperature.
At low masses, we find that experiments will be sensitive to ALPs with O(1) initial misalignment angles well before they are able to probe the QCD axion. At higher masses, the cosmological models motivate continuing ALP searches to couplings below the QCD region. In some cases, existing proposals will have the required sensitivity, while in other scenarios -particularly EMD -new search strategies may be required. We present the theoretical targets, existing constraints and experimental projections in Section IV, with the main results collected in Figs. 3, 4 and 5. Our findings are summarized in Section V. Appendices A and B contain details of the relic abundance calculations for different cosmologies and ALP mass temperature-dependence.

II. ALP DARK MATTER
We take the ALP Lagrangian to be whereF µν is the dual electromagnetic field-strength tensor. The photon coupling g aγγ = rα/(2πf a ) is related to the ALP decay constant f a . r is a model-dependent constant generically expected to be O(1); we set r = 1 in this work. The free parameters are then the ALP mass m a and photon coupling g aγγ (or equivalently f a ). At early times, the ALP field is frozen. The relic abundance today depends on the distribution of initial values a 0 ≡ f a θ 0 of the field before it begins to evolve. Here θ 0 is the initial misalignment angle. One possibility is that the angle θ 0 is uniform across all initially causallydisconnected regions that make up the observable universe today; this is the case if the ALP exists prior to inflation. Typically we expect θ 0 ∼ 1, in which case saturating the observed dark matter density identifies favored regions of the ALP parameter.
An alternative initial condition is a stochastic distribution of θ 0 across all causally-disconnected regions. This occurs in pNGB models where the global symmetry is broken after inflation. This scenario can be modeled by considering an ALP with an effective average misalignment angle of θ 0 ∼ π/ √ 3 [46]. In this case, topological defects and other large inhomogeneities formed during global symmetry-breaking also contribute to the present day dark matter density; however, the magnitude of these contributions is still a subject of debate. According to different studies the inclusion of large fluctuations may increase [47,48] or decrease [49] the relic density relative to that of the misalignment estimate. In what follows, we assume that the relic density is reasonably well-approximated by the misalignment calculation. Therefore, both the uniform and stochastic initial conditions can be studied if we vary θ 0 over a sufficient range. To avoid fine-tuning and to capture both possibilities, we take θ 0 ∈ [0.1, 2] below. 1 Accounting for topological defects should not significantly affect our conclusions, provided their contribution is at or below the same order of magnitude as the misalignment contribution.
The equation of motion for the ALP zero mode in the early universe is 2ä where H is the Hubble parameter ρ tot is the total energy density of the universe, and M p is the reduced Planck mass. At early times the ALP field is fixed. Oscillations begin when the Hubble parameter becomes comparable to m a , for some O(1) value of q. In Appendix A, we give a detailed discussion of q and list the values that give the best fits of the analytic formulae to the results of numerical integration. For temperature-independent ALP masses, we find q ≡ q 0 = 1.6 provides good precision across the various cosmological scenarios we consider. At a given time, the ALP energy density is where again we keep only the quadratic part of the ALP potential as an approximation. The ALP number density at time t can be defined as where we have allowed for the possibility of a timevarying ALP mass.
1 Monodromy scenarios allow a much larger initial misalignment of the ALP field -Refs. [50,51] consider displacements of up to 10 3 fa. This leads to larger possible values of gaγγ for a given value of ma. 2 The m 2 a a term should be replaced by V (a) for field values larger than O(fa). In what follows, we use the approximation above, noting that for larger initial misalignment angles, going beyond this approximation can have O(1) effects on the predicted relic density [52].
Soon after oscillations begin, the ALP energy density redshifts as matter. Let us denote the corresponding temperature as T osc . We also define a reference temperature T ad below which the evolution of the universe is adiabatic and the ratio of the ALP number density to entropy density n/s is conserved. At the onset of oscillations, the ALP number density is Its value at T ad is given by redshifting n a (T osc ) by (R osc /R ad ) 3 (R is the FRW scale factor). The presentday ALP density then depends on the cosmological evolution between T osc and T ad . The ALP relic density today can be written Here ρ c ≈ 10 −5 h 2 GeV/cm 3 , h ≈ 0.68, and T 0 ≈ 2.7 K [53]. In these expressions we have assumed that the ALP mass is temperature independent. We will return to the temperature-dependent case in Sec. III. Different cosmological scenarios correspond to different T ad and n(T ad ). Below, we consider three well-motivated possibilities in which the ALP begins to oscillate during a "standard" period of radiation domination, during a period of early matter domination followed by reheating, or during kination. The schematic evolution of the ALP energy density for these cosmologies is shown in Fig. 1.

A. Standard cosmology
In the conventional case, the ALP starts to oscillate during radiation domination (RD). The total energy density is given by where g * (T ) is the effective number of relativistic degrees of freedom. Away from mass thresholds, T ∼ R −1 and so ρ tot ∝ R −4 , where R is the scale factor. In this scenario, our approximate criterion for the onset of oscillations is where q 0 = 1.6, as discussed in Appendix A. Below T osc , the evolution is assumed to be adiabatic, so we can set T ad = T osc . Using Eqs. (8) and (10), we find an approximate expression for the ALP relic density today, assuming a temperature-independent mass during and  . Schematic evolution of the ALP energy density relative to the total energy as a function of the scale factor R for different cosmologies. The scale factor is normalized to unity at the start of ALP oscillations. The lettuce, mustard and tomato lines correspond to a universe with early matter (EMD), radiation, or kination domination before primordial nucleosynthesis, respectively. The transition from EMD or kination to radiation domination is denoted by the vertical dashed line. The initial ALP density is fixed by requiring that ALP-radiation equality occurs at the same value of R/Rosc for all three cases, such that these models have equal DM densities at late times. Since the initial value of the ALP energy density depends on faθ0, cosmologies with an early period of early matter (kination) domination, require larger (smaller) values of faθ0 to saturate the observed dark matter relic density than in standard radiation domination, for fixed ma.
after oscillations: Eq. (11) typically reproduces the results from numerical solutions of the ALP equation of motion (see Appendix B) to within about 10-20%. The scaling with input parameters is straightforward to understand: at the onset of oscillations, the ALP constitutes a fraction ρ a /ρ tot ∼ f 2 a /M 2 Pl of the total energy density, which immediately starts growing since ρ a ∝ R −3 redshifts more slowly than radiation. Correspondingly, larger f a leads to larger relic abundances. Similarly, increasing m a corresponds to earlier onset of oscillations and therefore a longer period over which ρ a /ρ tot grows, so the relic density also grows with m a .

B. Early Matter Domination
A period of early matter domination (EMD) modifies the conventional calculation of the axion relic density [54] (for recent work, see Refs. [42][43][44][45]). EMD can be modeled by a heavy long-lived particle or an oscillating scalar field φ that dominates the energy density, such that ρ tot ≈ ρ φ ∝ R −3 . This scalar field can be a saxion or another scalar modulus with small couplings that lead to long lifetimes. The entropy injected by the decay of the scalar field dilutes the energy density of the ALP below the reheating scale, allowing for larger initial ALP energy densities and reducing the tuning required in the misalignment angle for large f a . This cosmology therefore favors a different region of the ALP parameter space compared to the RD case described above.
In the EMD scenario, the ALP is again initially displaced from the origin and begins to oscillate when m a ∼ q 0 H = 1.6H if the mass is independent of temperature. Assuming a reheating temperature around 10 MeV (near the lower limit allowed by BBN [55,56]), oscillation occurs during EMD for m a 10 −13 eV, and during RD for smaller masses. The initial energy fraction in the ALP at oscillation is again of order f 2 a θ 2 0 /(M 2 pl ). The key difference in the EMD scenario is that H 2 ∝ R −3 after the onset of ALP oscillations and prior to reheating, and so the ALP energy fraction ρ a /(3M 2 Pl H 2 ) remains constant during this epoch. Accordingly, the ALP comes to dominate the energy density later than in the radiationdominated case, allowing for larger f a consistent with the observed dark matter relic density -see Fig. 1.
Assuming adiabatic expansion below T RH , the presentday ALP density is given by Eq. (8) with T ad = T RH . Since ρ a (t RH ) ρ a (t osc )H 2 (t RH )/H 2 (t osc ) and ρ a (t osc ) ∼ θ 2 0 f 2 a H 2 (t osc ), ρ a (T RH ) entering Eq. (8) is independent of H(t osc ) and m a . The present-day ALP density is found to be for temperature-independent ALP mass (see Appendix A for more details). In this case Ω a is determined by f a and T RH , with larger f a and T RH corresponding to larger Ω a . Note that this expression is only valid if T RH is larger than T osc ; otherwise, the standard radiation-dominated scenario is obtained. An expression for T osc is given by Eq. (40) below, where a temperature-independent ALP mass corresponds to b → 0.
We have also solved for the corresponding relic abundance numerically, by considering a three fluid model describing the modulus, the radiation energy density, and the ALP as described in Appendix. B. These numerical solutions agree with Eq. (12) to within about 20-25% across the parameter space of interest and for the values of T RH we have checked.

C. Kination
The final cosmological scenario we consider is known as kination [57,58]. ALP physics with an early period of kination has previously been studied in [41,43]. As for the EMD case, the energy density at early times is dominated by a long-lived scalar field φ, but rolling in a steep potential such that its kinetic energy dominates ρ tot . For a polynomial potential V (φ) ∝ φ N , the energy density after some early time t 0 evolves as In the limit N → ∞, ρ φ dilutes as ∼ R −6 . We will consider this large N limit in what follows. Assuming that radiation comes to dominate when it reaches a temperature T kin , using Eq. (8) with T ad = T kin we find that The relic density depends linearly on m a and inversely on T kin , which must be larger than about 5 MeV. The expression above only applies if T osc > T kin ; otherwise, one reproduces the standard RD scenario given by Eq. (11). An expression for T osc in the kination case is given in Eq. (43), where b → 0 corresponds to a temperatureindependent ALP mass. Note that the fractional density ρ a /ρ tot ∝ R 3 grows rapidly during kination, allowing the ALP to saturate the DM relic abundance for smaller values of f a compared to RD and EMD scenarios considered above. Eq. (14) reproduces the numerically-obtained relic density (c.f. Appendix B) to within O(10)%.
As an illustration of the key differences between the three scenarios discussed so far, we sketch the evolution of the various relevant energy densities in Fig. 1 for the three cosmologies. Here various parameters are fixed for illustrative purposes. The qualitative picture is clear: the faster the dilution of the dominant energy component in the pre-BBN era, the larger the final ALP abundance for fixed f a θ 0 . In the kination cosmology, for example, the ALP energy fraction rises more rapidly than in radiation domination. In contrast, in the EMD case this energy fraction remains constant until reheating. Sincein order to saturate the observed DM relic density -ALPradiation equality must occur around T 1 eV, the EMD and kination cases require a larger and smaller initial energy fraction, respectively, than in radiation domination, corresponding to larger and smaller preferred values of f a θ 0 for a given mass. From the experimental standpoint, this means that the kination scenario will provide a compelling and more easy-to-reach target than in the standard ALP cosmology, while a period of early matter domination will make the ALP more difficult to access with terrestrial experiments. However, all relic density-preferred bands can lie above the QCD band (i.e. at stronger coupling) for sufficiently small ALP masses. We will detail this picture further in Sec. IV.

III. ORIGIN OF THE ALP MASS
In the previous section, we assumed that the ALP mass is independent of temperature at the onset of oscillations. This is the simplest class of models, and in general it seems reasonable to remain agnostic about the origin of the ALP mass. However, motivated by the QCD axion, we consider two further variations.
Famously, the QCD axion appears to conflict with the straightforward application of effective field theory principles and the expectation that quantum gravity violates global symmetries [59][60][61][62][63]. Adding Planck-suppressed PQ-violating higher-dimension operators to the action, one finds that the axion solution to the strong CP problem is inoperative unless the Wilson coefficients are strongly suppressed up to operator dimension d ∼ 12. Solutions to this problem are known; it might be the case that all quantum gravity-induced PQ-violation is exponentially small [9,10]. In the ALP case, it is also of interest to compare the masses and couplings for which a viable dark matter candidate is obtained with the typical mass generated by Planck-suppressed operators.
Secondly, the QCD axion relic abundance is nontrivially affected by the strong temperature dependence of the topological susceptibility of QCD. Similarly, it is imaginable that the ALP mass is controlled by infrared physics (e.g., a new strongly coupled gauge theory) that introduces temperature dependence. As in QCD this dependence can have important implications for the preferred regions of the mass-coupling parameter space.

A. ALP mass from UV physics
We consider the typical contribution to the ALP mass from a dimension-d operator, parametrizing Φ = f a exp(ia/f a ). For simplicity, we suppose that the Wilson coefficient c is real and that the full potential is minimized at a = 0. The contribution to the ALP mass from Eq. (15) is We relate the scale f a to the ALP-photon coupling g aγγ by assuming [36] where r is an anomaly coefficient that we expect to be O(1). 3 Combining Eqs. (15) and (17) we obtain In the results presented in Sec. IV, we set c = r = 1. The resulting mass-coupling relation for d = 8, 10, 12 is shown along with the preferred DM regions and the experimental limits and projections in Figs. 3 and 4. To summarize, we will find that Planck-suppressed operators below dimension 8 must be absent across all of the parameter space we consider. In the high-f a region, even more suppression is required. For example, almost all of the viable ALP parameter space in the standard RD scenario requires that the Planck-suppressed contributions to the ALP potential start at dimension 12. The viability of this possibility depends on the specific UV model. We will discuss the implications of these results further below.

B. T -dependent ALP masses: general considerations
We now turn to the complementary case where the ALP mass is set by T -dependent infrared (IR) physics. First, we outline generic properties, constraints, and requirements on these scenarios. We then define a simple family of T -dependent masses and compute the relic density in the different cosmological scenarios, providing simple analytic expressions that reproduce the results of a more complete numerical treatment to within a few tens of percent.
First, we note that the temperature controlling the ALP mass does not need to equal the temperature of the SM bath. This is generically the case if the ALP mass is generated by couplings to a hidden sector (HS) that is not in kinetic equilibrium with the SM. For a given SM temperature T we parametrize the temperature of the hidden sector, T HS , as In what follows, all temperatures will correspond to temperatures of the SM photon bath, unless otherwise stated, and factors of ξ will be used to convert to hidden sector temperatures. We assume that m a is primarily sensitive to the temperature above a scale Λ, corresponding to a Standard Model (SM) bath temperature T Λ . The ALP zero mode is initially frozen at f a θ 0 and starts to oscillate when T = T osc . In order for T -dependence to have an effect on Ω a , we require ξ osc T osc > Λ where ξ osc ≡ ξ(T osc ).
The scale Λ cannot be arbitrarily low. In order for m a to vary significantly with temperature, there must exist a population of relativistic degrees of freedom in the HS. The presence of additional relativistic degrees of freedom modifies the expansion rate of the Universe, and which can alter the predictions of light element abundances and the CMB power spectrum. These constraints can be avoided if T Λ T BBN , where T BBN ∼ 5 MeV is the temperature of the SM bath around the onset of BBN. Otherwise, we must ensure that the effects from radiation in the HS at temperatures above T Λ are consistent with the measurements of the primordial abundances and CMB. Modifications of the expansion rate are typically parametrized by the effective number of neutrino species, N eff . For the parameter space of interest, T Λ is always above the temperature of recombination, so the BBN limit is most relevant. These constraints, detailed in e.g. Ref. [64], can be satisfied at ∼ 2σ confidence level provided ∆N eff 0.5. In terms of the effective number of relativistic degrees of freedom g * HS (T BBN ) in the HS at T BBN , we have where we use the shorthand ξ BBN ≡ ξ(T BBN ). From this we see that BBN N eff constraints can be avoided if ξ BBN 0.1 for g * HS (T BBN ) 10 3 . While the considerations above are quite general, there may be additional model-dependent constraints in concrete realizations. For example, one must also ensure that the relic abundance of any heavy states in the HS makes up a small component of the matter density today. These can decay or annihilate into HS radiation; however, one must then verify that ξ remains small. Heavy hidden sector states can also decay or annihilate to the SM, but this may require connector particles between the HS and SM that may again increase ξ(T BBN ). Furthermore, their decays to the SM must not significantly disrupt BBN or the CMB. In an effort to be as model-agnostic as possible we will not consider these issues further, although we emphasize that they will likely be important in concrete ALP scenarios with T -dependent masses. For related discussions in specific strongly coupled hidden sector models, see, e.g., Refs. [45,[65][66][67][68][69][70].
Summarizing, in order for T -dependence to affect the ALP relic density and be in agreement with N eff constraints, we require either or If ξ osc T osc < Λ, the onset of ALP oscillations proceeds as in the T -independent case discussed earlier. The relationship between T osc , m a , and f a depends on the particular cosmological scenario, as we discuss below. EMD and kination cosmologies will have additional requirements in order for T -dependence to be relevant.
Nontrivial temperature dependence enhances Ω a relative to the T -independent prediction for a given m a , f a , and θ 0 . The resulting abundance can be computed numerically (as discussed in Appendix B), but simple analytic estimates can again be used to reproduce the full results to within O(10%) in most cases. The size of the enhancement for the different cosmological scenarios can be estimated as follows (see also Appendix A for more details). Let us define the enhancement factor where Ω a is the ALP relic density assuming a Tdependent ALP mass and Ω T −ind a is the corresponding T -independent result as computed in Sec. II. Both quantities are evaluated for the same set of m a (T = 0), f a θ 0 . We model the different cosmological scenarios by assuming that the Hubble parameter for temperatures above some scale T * evolves as where w is the equation of state parameter, ρ = wp.
Early matter domination, radiation domination, and kination correspond to w = 0, 1/3, and 1, respectively. Below T * , the evolution is assumed to be adiabatic and follows that of a standard radiation-dominated cosmology. We approximate the transition to radiation domination at T * (if it occurs) as instantaneous. As in Eq. (4), the ALP begins to oscillate when Here the subscript T indicates that the value of q for the T -dependent case can differ from q 0 = 1.6. Given these assumptions and provided the ALP begins oscillating while its mass is changing with temperature, a straightforward calculation discussed further in Appendix A shows that the enhancement factor is given approximately by The subscript in γ T is indicates that this expression applies if the mass at the onset of oscillations, m osc , differs from m a , the low-temperature ALP mass. If m osc m a , the relic density can be significantly enhanced in the RD and EMD cosmologies. The scaling with m a /m osc is a product of two counter-acting effects: the delay in the start of oscillations and the growth of ALP mass with time. This is made explicit in Eq. A15 below. In the kination case, 2/(w + 1) − 1 = 0 and these effects nearly cancel, so the relic density is only enhanced if q T = q 0 . Given our assumptions about the origins of m a (T ), discussed below, this enhancement is milder than in RD and EMD, and is at most an O(1) effect.
To proceed further, we focus on a class of models with ALP mass T -dependence similar to that of the QCD axion. We will assume that, for T HS = ξ(T )T > Λ, the ALP mass is given by where m a is the zero-temperature mass, taken to be of the form For ξ(T )T < Λ, m a (T ) = m a . In Eq. (27), b is a positive exponent. In QCD-like theories, b is related to the β-function of the gauge group and can be obtained analytically from the dilute instanton gas approximation where N c and N f are the number of colors and light flavors, respectively [71]. For QCD, b = 4, and the semiclassical approximation is in reasonable agreement with lattice results at high temperatures [72][73][74]. In these simulations the scaling predicted by DIGA appears to hold down to Λ QCD , where m a saturates to near its zerotemperature value and remains approximately constant at lower temperatures [72,73]. In this sense our model of the temperature dependence mimics QCD and generalizes it to arbitrary Λ, b, and ξ(T ). In our plots we will take b = 4 as an illustrative example, corresponding to the QCD-like case. As such, we assume g * HS (T ) = 52 2 1 + tanh 10 1 − Λ ξ(T )T (29) where the factor of 52 corresponds to the number of relativistic degrees of freedom for SU (3) with three light flavors and the temperature T is understood to be that of the SM radiation bath; the tanh function smoothly decouples these degrees of freedom at the transition temperature. The total number of relativistic degrees of freedom at a temperature T is then g * (T ) = g * SM (T ) + ξ(T ) 4 g * HS (T ). Note that for the T -independent predictions we take g * HS = 0, as the mass m a can be set by physics in the ultraviolet and does not necessarily require new light degrees of freedom present near the onset of oscillations.
With the form of T -dependence specified, one can show (c.f. Appendix A) that there is a maximum allowed enhancement factor, γ T ≤ γ max . Defining T Λ such that T Λ = Λ/ξ(T Λ ) (the SM temperature at which the ALP mass saturates to its low-temperature value), if the ALP has not started oscillating by T Λ and m a > q 0 H(T Λ ), oscillations will begin suddenly at T Λ and so q max = m a /H(T Λ ). Thus, In this case the oscillation temperature is simply T osc = T Λ = Λ/ξ osc . Note that the opposite limit in which the ALP is still frozen at T Λ and m a < q 0 H(T Λ ) corresponds to ξ osc T osc < Λ, which reproduces the Tindependent case. One can show that this occurs when min(γ T , γ max ) < 1. Summarizing these considerations, the enhancement factor can be written compactly as where γ T and γ max are defined in Eqs. (26) and (30), respectively. Explicit expressions for these quantities in concrete cosmological scenarios are given below. Denoting the predicted oscillation temperature for a given exponent b in Eq. (27) as T osc (b), the true oscillation temperature is given by Further details can be found in Appendix A. Finally, in the EMD and kination cosmologies, if T osc predicted by Eq. (32) is smaller than T RH or T kin , the results for RD should be used.

C. Radiation domination with T -dependence
Let us first apply these results to determine the effects of T -dependence on the standard calculation of the ALP relic abundance, where the ALP is assumed to oscillate during radiation domination. A related discussion can be found in Ref. [8], which we generalize to allow for a HS at a different temperature than the SM. In general, ξ(T ) changes across mass thresholds as particles in both the HS and SM annihilate. This heating typically changes ξ by at most O(1) factors unless the change in number of degrees of freedom is very large. Since the precise form of ξ(T ) is model-dependent, in the remainder of this study we assume for simplicity that ξ(T ) ≈ ξ osc for temperatures of interest and treat ξ osc as a free parameter. In concrete models ξ osc can be set by, e.g., the branching ratio of the inflaton into the HS relative to the SM [75], or be equal to one for a HS in kinetic equilibrium with the SM.
In Appendix A, we find that the predicted oscillation temperature T osc (b) and enhancement in this case can be estimated by Meanwhile, the maximum enhancement factor γ max is With these expressions, Eqs. (11), (31) and (32) can then be used to estimate Ω a across the ALP plane accounting for T -dependence in the ALP mass. In the parameter space we consider, Eq. (31) typically yields γ = γ T for RD.
The preferred ALP dark matter region in the Tdependent case for ξ osc = 1 is shown on the left in Fig. 2 for b = 4 and g * HS (T ) given by Eq. (29) with ξ(T ) = ξ osc . The region shaded gold features an ALP with h 2 Ω a = 0.12 for natural values of the initial misalignment angle, θ 0 ∈ [0.1, 2], obtained by solving the ALP equation of motion (EOM) numerically (c.f. Appendix B). The gold dotted contours correspond to the analytic estimates given above. The corresponding Tindependent preferred ALP DM region, obtained by numerically solving the ALP EOM to late times, lies between the dashed gold contours. ξ osc = 1 illustrates the maximum allowed enhancement of the relic abundance in a RD cosmology. Since the HS is at the same temperature as the SM, there are strong bounds from N eff for ξ osc = 1, and the shaded gray region on the left in Fig. 2 is excluded by requiring ∆N eff < 0.5 (corresponding to T Λ > T BBN ∼ 5 MeV). It is likely that in concrete models the lower bound on T Λ will need to be somewhat higher than 5 MeV to avoid BBN constraints, and so the results shown should be understood to correspond to the most optimistic case.
On the right in Fig. 2 we show corresponding results assuming a decoupled hidden sector with ξ osc = 0.1. The enhancement of h 2 Ω a is smaller, however the cooler HS in principle allows for T Λ < T BBN , since ξ BBN < 0.1 (assuming the number of relativistic degrees of freedom in the hidden sector does not change significantly between oscillation and the onset of BBN). Again, one must be mindful of additional model-dependent constraints on small-Λ scenarios, as well as those with decoupled hidden sectors with dark radiation or significant late-time abundances of stable relics.
Larger ξ osc and b can in principle increase Ω a further, however this often comes at the cost of additional entropy injection after oscillation in simple models. For example, it could be that ξ osc > 1, however the large corresponding amount of HS entropy needs to be transferred to the SM before BBN, erasing the resulting enhancement for ξ osc much larger than 1. A hidden sector predicting b > 4 could also increase h 2 Ω a somewhat, however as b increases one also expects g * HS (T osc ) to increase in a QCD-like theory, and so the resulting enhancement again gets washed out by the requisite HS entropy dump for large b before BBN. These effects are encapsulated in the g * dependence of γ T in Eq. (34).

D. Early Matter Domination with T -dependence
We proceed similarly for the case of early matter domination, deriving a set of analytic expressions that can be used to estimate the relic abundance. We again allow the HS to be at a different temperature than the SM bath and parametrize ALP mass temperature dependence as in Eq. (27). The evolution of the energy densities in φ (the field responsible for EMD), SM and HS radiation can be modeled bẏ where the dot indicates a derivative with respect to time, Γ φ = Γ φ→SM +Γ φ→HS and Γ φ→SM , Γ φ→HS are the partial widths of φ into to SM and HS radiation, respectively. These equations can be solved during φ domination (i.e. while Γ φ /H 1 and ρ φ R 3 ≈ const) and yield where we assumed that the initial energy densities are negligible compared to those produced by φ decays. Here i corresponds to either HS or SM subscripts, and H MD and R MD are the Hubble parameter and FRW scale factor at the onset of EMD, respectively. We can use Eq. (39) to compute the temperatures of the HS and SM radiation baths during the epoch of matter domination. Defining the reheating temperature as that for which ρ φ = ρ SM + ρ HS and assuming Eq. (39) holds down to that temperature allows one to relate Γ SM and Γ HS to T RH and ξ osc (again neglecting relative heating effects). We then find that the oscillation temperature and abundance enhancement factor can be approximated by (see Appendix A for more details). γ max in EMD is given by The preferred regions of the ALP parameter space in the EMD scenario are illustrated in Fig. 2 for ξ osc = 1 and ξ osc = 0.1 with b = 4. We show results for T RH = 10 and 500 MeV. The blue and purple shaded regions feature an ALP with Ω a h 2 0.12 for θ 0 ∈ [0.1, 2] as obtained from the numerical solution. The corresponding dotted contours show the analytic predictions of Eqs. (40) - (42) and are a good fit to the numerical results. For ξ osc = 1, the gray shaded region is excluded by the measured value of N eff at BBN. This constraint is alleviated for ξ osc = 0.1, however the enhancement factor γ is reduced as a result. Other model-dependent constraints are likely to apply in the region where T Λ < T BBN as discussed in Sec. III B.
The behavior illustrated in Fig. 2 is straightforward to understand. First, note that the oscillation temperature is reduced as m a and f a are increased. If the oscillations begin below T Λ , temperature-dependent effects are unimportant and the preferred regions are the same as discussed in Sec. II B (i.e. γ = 1). This occurs for larger ALP masses, and the preferred value of g aγγ is independent of m a . For small enough masses, T osc > T Λ so that T -dependence enhances the relic density for a fixed m a , f a relative to the T -independent case. Here, the ALP DM regions pick up dependence on f a , increasing the preferred values of g aγγ in Fig. 2. At even lower values of m a , oscillation occurs after reheating (T osc < T RH ), and the predictions reduce to those of the T -dependent RD scenario of Sec. II A (the gold shaded region). Fig. 2 shows that allowing for a T -dependent ALP mass interpolates between the T -independent EMD and T -dependent RD scenarios. Smaller values for b tilt the interpolating region towards the left, while larger values steepen it. Increasing T RH causes the EMD band to match onto RD predictions at larger m a . In all cases, the preferred ALP DM regions are bounded by the Tdependent RD and T -independent EMD contours for a given θ 0 .

E. Kination with T -dependence
Finally, we comment on the kination cosmology with a T -dependent ALP mass near the onset of oscillations. From Eq. (26), the only enhancement comes from the slightly different values of q defining the oscillation time. In other words, the gain in energy from the growth of the mass is almost completely canceled by the loss in energy from starting to oscillate later. As explained in Appendix A, we find  . Theoretical targets (colored bands) and current experimental constraints (filled regions) on the ALP-photon coupling gaγγ as a function of the ALP mass ma. The shaded bands show regions where the ALP saturates the observed DM relic abundance for the standard (yellow), early matter-dominated (green) and kination (red) cosmologies for initial misalignment angles of θ0 ∈ [0. 1,2]. For the latter cosmologies, we take the reheating/kination temperature to be 10 MeV. We also show the QCD axion band, which does not have a relic density requirement imposed, in blue. The gray dotted diagonal lines correspond to ALPs which get their mass from dimension 8, 10 and 12 Planck-suppressed operators. Further discussion can be found in Section IV.

Meanwhile,
In most realistic models, one expects b ∼ O(1), and so typically γ = γ T . The enhancement is milder than in RD and EMD as it only depends on the exponent b. Nevertheless, allowing for m osc = m a changes the ALP oscillation temperature. Since T osc > T kin in order for the period of kination to modify the ALP evolution, T -dependence will change the regions of the ALP plane where kination is relevant for fixed T kin .
In our results below we will take T kin = 10 MeV. Since the preferred regions on the ALP parameter space assuming kination with and without T -dependence are similar for T kin = 10 MeV, we do not show predictions for kination in Fig. 2. However, we provide the corresponding T -dependent predictions in Fig. 5. Again we find that the analytic estimates above provide a good fit to the numerics (c.f. Appendix B) across the parameter space considered.

IV. PROJECTIONS AND RESULTS
We now investigate the potential for current and future ALP direct detection experiment and astrophysical observations to explore these natural ALP dark matter targets. The present status is summarized in Fig. 3, and future prospects are shown in in Figs. 4 (for temperatureindependent masses) and 5 (for temperature-dependent masses). These figures also show the preferred regions in the three cosmological histories considered in Secs. II and III: the standard cosmology, early matter domination (EMD) with T RH = 10 MeV, and kination with T kin = 10 MeV. T RH and T kin are the temperatures at which the universe transitions to standard radiation-  For the latter cosmologies, we take the reheating temperature to be 10 MeV. We also show the standard QCD axion target in blue and existing experimental constraints in solid gray. The gray dotted diagonal lines correspond to ALPs which obtain mass from dimension 8, 10 or 12 Planck scale suppressed operators. Further discussion can be found in Section IV.
are generated by Planck-suppressed operators of various dimensions as discussed in Sec II. For the temperaturedependent results in Fig. 5 we have assumed b = 4 in Eq. (27) and taken g * HS as in Eq. (29) with ξ = ξ osc . The left hand-plot in Fig. 5 shows the case where the hidden sector is in thermal equilibrium the SM, and the right-hand plot the case where the HS is decoupled from the SM with a lower temperature, 0.1 × T SM . Axion-like particles can be constrained by a variety of astrophysical measurements. These limits include the results from CAST [76]; cooling of Horizontal Branch ("HB" in Fig. 3) stars, massive stars [77,78], and SN1987A [79][80][81]; non-observation of a γ-ray excess from SN1987A [82]; the extragalactic background light [83,84]; searches for spectral irregularities in γ rays with HESS [85] and Fermi-LAT [86], and in X-rays with Chandra [87]. 4 These limits are shown in Fig. 3 as pastelcolored shaded regions.
The ALP parameter space is also constrained by a number of resonant cavity experiments. We show the regions excluded by ADMX [16,17] and ADMX Sidecar [89], Phase 1 of HAYSTAC [20], the ORGAN Pathfinder [18] and the older UF [90] and RBF [91,92] experiments in dark blue in Fig. 3. These experiments target the classical QCD axion DM window for m a between 10 −6 and 10 −4 eV. We see in Fig. 3 that the resonant cavity experiments are already probing significant regions of the kination-favored parameter space and are just beginning to extend into the QCD axion window.
We turn now to near-term prospects for direct detection in the ALP parameter space. The past few years have seen a renaissance in ideas for searching very light DM, including coherent bosonic candidates like ALPs. We show in Fig. 4 a summary of the impact these new experiments will have on the ALP parameter space for temperature-independent ALP masses, and they are sensitive to the properties of the ALP self-interactions. Accordingly we omit them from our plots, noting that they impact the region of parameter space ma 10 −11 eV. T -dependent ALP mass, ξosc = 0.1 Figure 5. As in Fig. 4, but for a T -dependent ALP mass. In the left panel, the hidden sector responsible for generating the ALP potential is assumed to be in thermal equilibrium with the SM, while in the right panel we assume the hidden sector is decoupled with temperature given by 0.1 × TSM. In both cases b = 4 was assumed in Eq. (27). Note that there may be additional important constraints on the hidden sector, as discussed further in the text. On the left, the EMD band assumes TRH = 10 MeV, while on the right TRH = 500 MeV. In both cases T kin = 10 MeV. The parameter space on the left is constrained by ∆N eff at BBN, since there are necessarily new HS states in equilibrium with the SM bath at TΛ. The ALP target regions assume that g * HS = 52 above TΛ, corresponding to the expected value for a SU (3) hidden sector with 3 light flavors.
in Fig. 5 a similar summary for temperature-dependent ALP masses. In many cases, allowing for T -dependence the experimental prospects are even more promising, although constraints on new relativistic degrees of freedom generating the ALP potential can exclude some of the parameter space. We emphasize that these considerations are model-dependent and that specific scenarios could feature even more stringent constraints on the hidden sector than those considered. Each experiment is capable of ruling out the region above the corresponding solid line.
Some future experiments are extensions of resonant microwave cavities technique, as in upgrades to ADMX [93], CAPP [94], KLASH [95,96], and, at higher frequencies, ORGAN [18]. These experiments provide a broader sensitivity in the QCD axion region m a ∼ 10 −6 − 10 −5 eV extending to lower values of g aγγ . More recently, new ideas based on dielectric stacks have appeared which are sensitive to higher mass ALPs, as in MADMAX [21] and photonic materials [23] ("Dielectric Stack" in Figs. 4 and 5). We also show the sensitivity of the proposal for a large-scale helioscope, IAXO [97], which will extend the reach of CAST, and the projections for the ALPS-II light-shining-through-walls experiment [98], which is currently under construction at DESY and will have sensitivity above g aγγ ∼ 10 −11 GeV −1 . It is also possible that future measurements of radio emission lines from the magnetospheres of neutron stars could lead to constraints in the µeV mass-range [99]. We show the limits that could be obtained with 100 hours of observation of the magnetar SGR J1745-2900 under the assumptions of an NFW dark matter density profile ("NSM" in Figs. 4 and 5), and also a spike profile ("NSM Spike") which would lead to stronger bounds. At large ALP masses the intensity line-mapping experiment SPHEREx [100] will be able to probe to the bottom of the kination region.
At very low masses, the ABRACADABRA suite of experiments (the region we show is the union of the broadband and resonant searches) and DM-Radio promise to cover a large amount of parameter space down to very small values of g aγγ 5 .
Other recent proposals at low mass make use of birefringence in the presence of an ALP background and include the interferometer concept [110], ADBC [111], and an experiment based on optical ring cavities [112].
While these will be able to explore new parts of parameter space, we find that they will not be sensitive to the kinds of ALP dark matter we study in this paper. We find that DM-Radio will be able to probe ALP dark matter up to m a ∼ 10 −6 eV assuming a standard cosmology or a period of kination in the early Universe. ABRACADABRA will be able to discover (or rule out) ALP dark matter in all of the cosmological scenarios we have considered with masses below ∼ 4 × 10 −7 eV. If the ALP mass is generated by a new strongly coupled gauge sector, the signal at ABRACADABRA for a given mass is likely to be even larger.

V. SUMMARY AND CONCLUSIONS
We have investigated the implications of current and future direct detection experiments for ALP dark mat-ter with mass 10 −12 ≤ m a ≤ 1 eV in a variety of wellmotivated cosmological scenarios. We have presented simple analytic expressions for the corresponding relic density from misalignment in the standard cosmological scenario with radiation domination (RD), as well as allowing for a period of early matter domination (EMD) and kination, in Eqs. (11), (12), and (14), respectively. These results apply to ALPs for which the mass is independent of the temperature between the onset of oscillations and today. A T -dependent ALP mass of the form in Eq. (27) enhances the relic abundance relative to these predictions so that Ω a = γΩ T −ind a , with the enhancement factor γ given by Eq. (31) and Eqs. While ALP dark matter is currently relatively unconstrained, future experiments have the ability to probe much of the well-motivated ALP parameter space. ALPs that obtain their masses from Planck-suppressed operators will be thoroughly tested by future experiments (provided they can saturate the observed dark matter relic abundance). The amount of suppression required for a viable ALP dark matter candidate depends on the cosmological scenario under consideration. It is possible for an ALP associated with d = 12 operators to be consistent with the DM relic density in a standard cosmological scenario for masses above 10 −2 eV. In other cases, such as a period of kination down to temperatures of a few MeV, an ALP with mass set by d = 8 Planck-suppressed operators can provide a viable dark matter candidate, a significantly less stringent requirement than that for the QCD axion (which requires Planck-suppressed PQbreaking operators to arise at d = 12 or higher).
ALP dark matter can be easier to detect than the QCD axion. For low masses (below the standard QCD axion window for a fixed f a ) experiments such as ABRA-CADABRA and DM-Radio will have sensitivity to ALP dark matter before they are sensitive to the QCD axion. In particular, the ABRACADABRA experiments can constrain the existence of ALPs in the various cosmologies we have considered for m a 4 × 10 −7 eV down to 10 −12 eV. Below 10 −11 eV black-hole superradiance complements the ABRACADABRA and DM-Radio sensitivity, although the precise details are modeldependent. If there was a period of kination in the early Universe down to temperatures near the BBN scale, or if the ALP mass at the onset of oscillations is smaller than its present day value, ABRACADABRA and DM-Radio can be more sensitive to ALP dark matter than to the QCD axion across their entire mass sensitivity ranges. At higher m a , the experimental prospects are more positive in the kination and T -dependent RD cases as well. We find that resonant cavity experiments (such as ADMX, CAPP and ORGAN), as well as MADMAX, can also probe ALPs in these more optimistic scenarios before they reach the QCD axion window.
For even larger masses, m a 10 −4 eV, ALP dark matter becomes more difficult to detect than the QCD axion in all of the scenarios we have considered. However, some proposed experiments using terahertz frequency resonators, dielectric stacks, or line-intensity mapping targeting the QCD axion in this mass range can also probe ALP DM that begins oscillating during kination (for low T kin ) and come close to the standard ALP prediction with O(1) initial misalignment angles and Tdependent masses for m a up to an eV. For ALPs in the standard RD and EMD cosmologies with masses set in the UV, this high-m a region will be difficult to access with existing experimental proposals. However, other probes of this parameter space beyond direct detection experiments may exist in some cases. For example, a period of early matter domination can also lead to the formation of ALP miniclusters, which can have interesting astrophysical consequences [43,44,113]. Future inquiry along these lines, and new ideas to access this region experimentally, are worth continued investigation.
ALPs can provide a compelling and viable dark matter candidate, behaving much like the QCD axion in the early Universe, but in many cases allowing for larger couplings to photons. ALP dark matter, therefore, can be easier to detect than the QCD axion, especially at low masses. More generally we emphasize the importance of vigorously pursuing the axion direct detection program, targeting a wide range of masses and exploring the ALP parameter space beyond the canonical QCD axion window.

Temperature-Independent ALP Mass
First we consider cases where m a is fixed to its zero-temperature value at times before the onset of ALP oscillations. The relic density is then given by Eq. (8), reproduced here for convenience: T ad is a reference temperature below which the evolution of the universe is adiabatic and T 0 is the temperature today. T osc is the SM temperature at the onset of ALP oscillations, and we haveθ(t osc ) 0 and n a (T osc ) = 1/2m a f 2 a θ 2 0 . Below T osc the comoving ALP number density is assumed to be conserved, so that Combining Eqs. (A1) -(A2) yields a general expression for the relic density, To obtain Ω a we therefore need to determine T osc /T ad and/or R osc /R ad in the various cases.
We parametrize the relationship at oscillation between the Hubble parameter and ALP mass as Here q 0 is a positive number that should be chosen to accurately reproduce the numerical predictions, described in Appendix B. In our final estimates we will take q 0 = 1.6; the reasoning behind this choice is explained below. At the time when oscillations begin, it is assumed that the universe is dominated by a fluid with equation of state p = wρ (for which H 2 ∝ R −3(w+1) ), and then later transitions instantaneously to radiation domination at some T = T * : In all cases we consider (RD, EMD, and kination), we can take Henceforth we will refer only to T * . In RD, we can also set T * = T osc . In EMD, T * = T RH , and in kination T * = T kin . Combining Eqs. (A4) -(A5), we obtain R osc /R * in terms of T * and m a . Plugging in to Eq. (A3) and setting q 0 = 1.6 yields the final result for the relic density for each case, given in the main text as Eqs. (11), (12), and (14). Note that in the kination case we could have instead taken T ad = T osc and used the conservation of comoving entropy to relate T osc to T kin . This yields the same result.
We now discuss q 0 in more detail, which will be particularly relevant when the ALP mass is T -dependent. First we inspect the form of solutions to the ALP EOM. Consider the EOM during a period with equation of state p = wρ, such that R ∼ t 2 3(w+1) . With θ(t) ≡ a(t)/f a we havë which has solutions Here r = 1 2 − 1 1+w , c 1,2 are integration constants, and J r (x), Y r (x) are Bessel functions of the first and second kind. Given that the Bessel functions only exhibit oscillatory behavior when their arguments are O(1) or larger, we see that ALP oscillations begin when with A an O(1) number. Since H 2/(3t(1 + w), we define the onset of oscillations as For EMD (w = 0), RD (w = 1/3), and kination (w = 1), we obtain m a {3/2, 2, 3} × AH osc , respectively. Comparing our analytic and numerical solutions, we find that choosing A such that m a 1.6H osc reproduces the numerical results to within a few tens of percent across the parameter space considered in the various cosmologies. (This appears consistent with the discussion of Ref. [12], which found m a ∼ 2H osc is a better choice than m a ∼ 3H osc in the RD scenario.) This value of q 0 can be adjusted to yield slightly better agreement in each cosmology, but for simplicity we take a common value. Therefore, introducing the parametrization (A10) was not really necessary in this case; however, a similar parametrization is useful when considering temperature-dependent masses, so we keep it for comparison. Summarizing, we take q 0 = 1.6 for all cosmologies, or in the parametrization (A10), 16 15 , EMD (A11)

Temperature-Dependent ALP Mass
We can proceed similarly when m a (T ) varies near the onset of oscillations. As discussed in the main text, the temperature controlling the ALP mass does not need to equal the temperature of the SM bath. We parametrize the temperature of the hidden sector T HS as where T is the temperature of the SM photon bath. We take an instanton-motivated class of models in which m a (T HS ) = (Λ/T HS ) b m a for T > Λ. Subsequently, all temperatures will correspond to SM temperatures, unless otherwise stated, and factors of ξ will be used to convert to hidden sector temperatures. Denoting the ALP mass at the onset of oscillations as m osc and the mass today as m a , the relic density can be expressed as: Here we have again assumed that at the time when oscillations begin, the universe is dominated by a fluid with equation of state p = wρ (for which H 2 ∝ R −3(w+1) ), and then later transitions instantaneously to radiation domination at some T = T * so that Eq. (A5) applies. We parametrize the oscillation time in the case of temperature-dependent ALP masses via generalizing Eq. (A4). q T should again be a positive number chosen so that the analytic formulae provide a good approximation to the numerical results. Note that q T and q 0 can be different. (A15) Thus, to determine Ω a in the T -dependent case, we can use the T -independent results of Eqs. (11), (12), and (14), multiplying by γ (assuming that T -dependence is relevant at the time of oscillations) and using the value of T osc predicted for T -dependent masses. To do so, we must determine m osc and appropriate choices of q T in each case. How should q T be chosen? Consider the ALP EOM for times when m a (T ) = m a (Λ/ξ(T )T ) b . During a period of adiabatic evolution dominated by a fluid with equation of state ρ = wp, the temperature evolves as away from mass thresholds. As a result, the T -dependent mass term in the ALP EOM can be written in terms of t so that the EOM becomesθ where T i is a constant (neglecting the T -dependence of g * and ξ). The solutions of this equation are sightly more complicated than in the T -independent case, but can still be written in terms of Bessel functions with arguments (A18) The solution starts to oscillate when the above quantity ∼ O(1) ≡ A. Using H = 2/(3(1 + w)t) and m a (t) = m a t Since we have assumed adiabatic evolution, this result applies to our radiation-dominated and kination scenarios.
For the EMD cosmology, entropy is injected into the bath from decays. However, we can still derive an approximate time-temperature relation from Eq. (39) (again neglecting T -dependence in g * and ξ). The ALP EOM in this case is approximatelyθ Proceeding as before, we find that the onset of oscillations occurs when When comparing to our numerical solutions, we find that neglecting the b-dependence in q T for b = 4 typically results in ∼ 70 − 100% discrepancies from the numerics. The disagreement becomes worse for larger values of b, in which case choosing the correct value of q T becomes particularly important.
Finally, we obtain the oscillation temperature from m osc = q T H osc , assuming that m a (T ) = m a (Λ/ξT ) b at oscillation. We refer to this oscillation temperature as T osc (b) to distinguish it from the T -independent prediction. To obtain T osc (b), we must specify H(T ). In RD this is simple. In kination, the temperature of the radiation bath is given by conservation of the comoving entropy density once T kin is specified, and the scale factor evolution from T kin to T osc is determined by Eq. (A5). For EMD, one can use the late-time solution for the radiation energy density, Eq. (39), along with Eq. (A5), to write the Hubble parameter in terms of the temperature. Solving for T osc (b) in this way yields Eqs. (33), (40), and (43). Inserting these temperatures into m a (T ) and Eq. (A15) yields Eqs. (34), (41) and (44) for γ T . This is not quite the end of the story for the temperature-dependent scenarios. Our results for q T , T osc (b), and γ T followed from inspecting solutions to the ALP EOM with m a (T ) = m a (Λ/ξT ) b . Therefore, they only apply if q T H Λ ≤ m a , where H Λ is the Hubble parameter at T Λ = Λ/ξ(T Λ ), when the mass saturates to its T -independent value. If instead q T H Λ > m a , oscillations begin after the mass has already saturated to m a , and the ALP EOM solution is given by Eq. (A8). Then there are two possibilities we must consider. First, if q 0 H Λ > m a , we recover the T -independent case. Second, if q 0 H Λ < m a , at T Λ the ALP will start to rapidly oscillate across this threshold, and H osc is given by m a (T osc ) = m a , H osc H Λ , q = q max ≡ m a H Λ .
1/(3t 0 ). The resulting r(t) yields H(t) through Eq. (B5), while T (t) is defined from g * (T )T 4 = g * (T kin )T 4 kin r −4 (t). (B7) In modeling a period of early matter domination, the radiation and matter densities are coupled and given as a function of time by numerically solving the systeṁ with H = 1/( √ 3M Pl ) √ ρ φ + ρ R and Γ φ chosen so that ρ R = ρ φ at a temperature T RH . The temperature is defined from the radiation energy density ρ R = π 2 /30g * (T )T 4 in the usual way. The final time is chosen to correspond to temperatures below T RH so that Eq. (B2) can be used. It is worth noting that the numerical treatment of the radiation bath outlined above is technically not entirely correct: near the QCD phase transition (or when the relativistic HS DOFs annihilate) the radiation bath equation of state can deviate somewhat from ρ = 1/3p. Conservation of comoving entropy density can instead be used in this case to track the evolution of ρ R with t. However, we find that the corresponding effects are rather small and so we neglect them in our analysis.