Dark Matter Targets for Axion-like Particle Searches

TL;DR

This paper identifies experimental targets for axion-like particle dark matter detection across various cosmological scenarios, highlighting some experiments' potential to detect ALPs before the QCD axion and emphasizing the importance of mass origin.

Contribution

It provides benchmark ALP dark matter targets for different cosmologies, considering temperature-dependent masses and misalignment angles, guiding future experimental searches.

Findings

Some experiments can detect ALP DM targets before the QCD axion.

Temperature-dependent ALP masses improve detection prospects.

Certain cosmological scenarios yield lower decay constant targets.

Abstract

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.

Equations

$${\mathcal{L}}_{\mathrm{ALP}}=\frac{1}{2}{\partial }_{\mu }a{\partial }^{\mu }a-\frac{1}{2}{m}_a^2{a}^2-\frac{1}{4}{g}_{a\gamma \gamma }a{F}_{\mu \nu }{\tilde{F}}^{\mu \nu },$$

$$\ddot{a}+3H\dot{a}+m_a^{2}a=0,$$

$$H^2=\frac{{\rho }_{\mathrm{tot}}}{3M_p^2},$$

$$m_a=qH,$$

$${\rho }_a=\frac{1}{2}\left({\dot{a}}^2+m_a^2{a}^2\right)$$

$$n_a(t)={\rho }_a(t)/m_a(t)$$

$$n_a(T_{\mathrm{osc}})=\frac{1}{2}{m}_a^2{f}_a^2{\theta }_0^2.$$

$${\Omega }_a=\frac{m_a{n}_a(T_{\mathrm{ad}})}{{\rho }_c}{\left(\frac{T_0}{T_{\mathrm{ad}}}\right)}^3\frac{g_{\ast S}(T_0)}{g_{\ast S}(T_{\mathrm{ad}})}.$$

$${\rho }_{\mathrm{tot}}=\frac{{\pi }^2}{30}{g}_{\ast }(T)T^4,$$

$$m_a=q_{0}H(T_{\mathrm{osc}})=q_0\sqrt{\frac{{\pi }^2}{90}}\frac{g_{\ast }^{1/2}(T_{\mathrm{osc}})T_{\mathrm{osc}}^2}{M_p}$$

$${\Omega }_a{h}^2$$

$$\times {\left(\frac{90}{g_{\ast }(T_{\mathrm{osc}})}\right)}^{1/4}.$$

$${\Omega }_a{h}^2\simeq 0.12\times {\left(\frac{f_a{\theta }_0}{9\times 10^{14}\mathrm{GeV}}\right)}^2\times \left(\frac{T_{\mathrm{RH}}}{10\mathrm{MeV}}\right)$$

$${\rho }_{\varphi }(t)\simeq {\rho }_{\varphi }(t_0)\times {\left(\frac{R(t_0)}{R(t)}\right)}^{\frac{6N}{N+2}}.$$

$$\begin{array}{rl}{\Omega }_a{h}^2\simeq 0.12& \times {\left(\frac{f_a{\theta }_0}{2.4\times 10^{11}\mathrm{GeV}}\right)}^2\times \left(\frac{m_a}{1\mu \mathrm{eV}}\right)\\ & \times {\left(\frac{11}{g_{\ast }(T_{\mathrm{kin}})}\right)}^{\frac{1}{2}}\times \left(\frac{10\mathrm{MeV}}{T_{\mathrm{kin}}}\right).\end{array}$$

$$\mathcal{L}\supset \frac{c{\Phi }^d}{M_p^{d-4}}+h.c.,$$

$$m_a^2=2c{d}^2{f}_a^2{\left(\frac{f_a}{M_p}\right)}^{d-4}.$$

$$g_{a\gamma \gamma }=r\frac{\alpha }{2\pi f_a},$$

$$g_{a\gamma \gamma }=\frac{r\alpha }{2\pi M_p}{\left(\frac{m_a^2}{2c{d}^2{M}_p^2}\right)}^{-1/(d-2)}.$$

$$T_{\mathrm{HS}}\equiv \xi (T)T$$

$$\Delta N_{\mathrm{eff}}=\frac{4}{7}{\left(\frac{11}{4}\right)}^{4/3}{g}_{\ast \mathrm{HS}}(T_{\mathrm{BBN}}){\xi }_{\mathrm{BBN}}^4$$

$${\xi }_{\mathrm{osc}}{T}_{\mathrm{osc}}>\Lambda ,T_{\Lambda }>T_{\mathrm{BBN}}$$

$${\xi }_{\mathrm{osc}}{T}_{\mathrm{osc}}>\Lambda ,{\xi }_{\mathrm{BBN}}\lesssim \mathrm{0.1.}$$

$$\gamma \equiv \frac{{\Omega }_a}{{\Omega }_a^{T-\mathrm{ind}}}$$

$$H^2\propto R^{-3(w+1)}$$

$$m_{\mathrm{osc}}=q_T{H}_{\mathrm{osc}}.$$

$${\gamma }_T\simeq {\left(\frac{q_T}{q_0}\right)}^{\frac{2}{w+1}}{\left(\frac{m_a}{m_{\mathrm{osc}}}\right)}^{\frac{2}{w+1}-1}.$$

$$m_a(T)=m_a{\left(\frac{\Lambda }{\xi (T)T}\right)}^b$$

$$m_a=\frac{{\Lambda }^2}{f_a}.$$

$$g_{\ast \mathrm{HS}}(T)=\frac{52}{2}\left(1+\tanh \left[10\left(1-\frac{\Lambda }{\xi (T)T}\right)\right]\right)$$

$${\gamma }_{\max }={\left(\frac{m_a}{H(T_{\Lambda })q_0}\right)}^{\frac{2}{w+1}}.$$

Full text

††preprint: FERMILAB-PUB-19-197-A-T

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 Peccei and Quinn (1977a, b); Wilczek (1978); Weinberg (1978). 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 Abbott and Sikivie (1983); Dine and Fischler (1983); Preskill et al. (1983), 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 axion-like 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 Arias et al. (2012); Svrcek and Witten (2006); Arvanitaki et al. (2010); Cicoli et al. (2012). 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., Marsh (2016); Hook (2018).

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. Graham et al. (2015); Irastorza and Redondo (2018)). These include resonant cavity experiments at various frequencies, such as ADMX Asztalos et al. (2010); Du et al. (2018), ORGAN McAllister et al. (2017) and HAYSTAC Brubaker et al. (2017); Zhong et al. (2018), and also new ideas including dielectric haloscopes (MADMAX Caldwell et al. (2017); Brun et al. (2019) and photonic materials Baryakhtar et al. (2018)), resonant LC-circuits Chaudhuri et al. (2015); Silva-Feaver et al. (2017), detection-induced magnetic flux oscillations (ABRACADABRA Kahn et al. (2016); Ouellet et al. (2018)) and NMR-based techniques (ARIADNE Arvanitaki and Geraci (2014); Geraci et al. (2018) and CASPEr Graham and Rajendran (2013); Budker et al. (2014); Garcon et al. (2017)). Collectively 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 Arias et al. (2012).

Recently, the Physics Beyond Colliders Working Group has forecast the sensitivity of future experiments to axions and ALPs Beacham et al. (2019), building on the review Irastorza and Redondo (2018). For the QCD axion, a number of groups have developed models to expand the parameter space, classifying the possibilities for UV-complete theories Di Luzio et al. (2017a, b); Agrawal et al. (2018a), model-building photophilic Farina et al. (2017) and photophobic Craig et al. (2018) axions, and extending the standard misalignment mechanism Agrawal et al. (2018b); Co et al. (2018). 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\gamma\gamma}$ 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. Visinelli and Gondolo (2010); Ramberg and Visinelli (2019); Visinelli and Redondo (2018); Nelson and Xiao (2018); Draper et al. (2018), 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 $\mathcal{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

[TABLE]

where $\tilde{F}^{\mu\nu}$ is the dual electromagnetic field-strength tensor. The photon coupling $g_{a\gamma\gamma}=r\alpha/(2\pi f_{a})$ is related to the ALP decay constant $f_{a}$. $r$ is a model-dependent constant generically expected to be ${\cal O}(1)$; we set $r=1$ in this work. The free parameters are then the ALP mass $m_{a}$ and photon coupling $g_{a\gamma\gamma}$ (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}\equiv f_{a}\theta_{0}$ of the field before it begins to evolve. Here $\theta_{0}$ is the initial misalignment angle. One possibility is that the angle $\theta_{0}$ is uniform across all initially causally-disconnected regions that make up the observable universe today; this is the case if the ALP exists prior to inflation. Typically we expect $\theta_{0}\sim 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 $\theta_{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 $\langle\theta_{0}\rangle\sim\pi/\sqrt{3}$ Kolb and Turner (1990). 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 Hiramatsu et al. (2012); Gorghetto et al. (2018) or decrease Klaer and Moore (2017) 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 $\theta_{0}$ over a sufficient range. To avoid fine-tuning and to capture both possibilities, we take $\theta_{0}\in[0.1,2]$ below.111Monodromy scenarios allow a much larger initial misalignment of the ALP field – Refs. Jaeckel et al. (2017); Berges et al. (2019) consider displacements of up to $10^{3}f_{a}$. This leads to larger possible values of $g_{a\gamma\gamma}$ for a given value of $m_{a}$. 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 is222The $m_{a}^{2}a$ term should be replaced by $V^{\prime}(a)$ for field values larger than $\mathcal{O}(f_{a})$. In what follows, we use the approximation above, noting that for larger initial misalignment angles, going beyond this approximation can have $\mathcal{O}(1)$ effects on the predicted relic density Turner (1986).

[TABLE]

where $H$ is the Hubble parameter

[TABLE]

$\rho_{\mathrm{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}$,

[TABLE]

for some ${\cal 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\equiv q_{0}=1.6$ provides good precision across the various cosmological scenarios we consider. At a given time, the ALP energy density is

[TABLE]

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

[TABLE]

where we have allowed for the possibility of a time-varying ALP mass.

Soon after oscillations begin, the ALP energy density redshifts as matter. Let us denote the corresponding temperature as $T_{\rm osc}$. We also define a reference temperature $T_{\mathrm{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

[TABLE]

Its value at $T_{\mathrm{ad}}$ is given by redshifting $n_{a}(T_{\rm osc})$ by $(R_{\rm osc}/R_{\rm ad})^{3}$ ($R$ is the FRW scale factor). The present-day ALP density then depends on the cosmological evolution between $T_{\mathrm{osc}}$ and $T_{\mathrm{ad}}$. The ALP relic density today can be written

[TABLE]

Here $\rho_{c}\approx 10^{-5}\,h^{2}\mathrm{GeV/cm}^{3}$, $h\approx 0.68$, and $T_{0}\approx 2.7\,K$ Tanabashi et al. (2018). 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_{\mathrm{ad}}$ and $n(T_{\mathrm{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.

II.1 Standard cosmology

In the conventional case, the ALP starts to oscillate during radiation domination (RD). The total energy density is given by

[TABLE]

where $g_{*}(T)$ is the effective number of relativistic degrees of freedom. Away from mass thresholds, $T\sim R^{-1}$ and so $\rho_{\mathrm{tot}}\propto R^{-4}$, where $R$ is the scale factor. In this scenario, our approximate criterion for the onset of oscillations is

[TABLE]

where $q_{0}=1.6$, as discussed in Appendix A. Below $T_{\rm osc}$, the evolution is assumed to be adiabatic, so we can set $T_{\mathrm{ad}}=T_{\rm osc}$. Using Eqs. (8) and (10), we find an approximate expression for the ALP relic density today, assuming a temperature-independent mass during and after oscillations:

[TABLE]

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 $\rho_{a}/\rho_{\mathrm{tot}}\sim f_{a}^{2}/M_{\rm Pl}^{2}$ of the total energy density, which immediately starts growing since $\rho_{a}\propto 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 $\rho_{a}/\rho_{\mathrm{tot}}$ grows, so the relic density also grows with $m_{a}$.

II.2 Early Matter Domination

A period of early matter domination (EMD) modifies the conventional calculation of the axion relic density Banks and Dine (1997) (for recent work, see Refs. Ramberg and Visinelli (2019); Visinelli and Redondo (2018); Nelson and Xiao (2018); Draper et al. (2018)). EMD can be modeled by a heavy long-lived particle or an oscillating scalar field $\phi$ that dominates the energy density, such that $\rho_{\mathrm{tot}}\approx\rho_{\phi}\propto 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}\sim 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 Kawasaki et al. (2000); Hannestad (2004)), oscillation occurs during EMD for $m_{a}\gtrsim 10^{-13}$ eV, and during RD for smaller masses. The initial energy fraction in the ALP at oscillation is again of order $f_{a}^{2}\theta_{0}^{2}/(M_{\rm pl}^{2})$. The key difference in the EMD scenario is that $H^{2}\propto R^{-3}$ after the onset of ALP oscillations and prior to reheating, and so the ALP energy fraction $\rho_{a}/(3M_{\rm Pl}^{2}H^{2})$ remains constant during this epoch. Accordingly, the ALP comes to dominate the energy density later than in the radiation-dominated case, allowing for larger $f_{a}$ consistent with the observed dark matter relic density – see Fig. 1.

Assuming adiabatic expansion below $T_{\rm RH}$, the present-day ALP density is given by Eq. (8) with $T_{\rm ad}=T_{\rm RH}$. Since $\rho_{a}(t_{\rm RH})\simeq\rho_{a}(t_{\rm osc})H^{2}(t_{\rm RH})/H^{2}(t_{\rm osc})$ and $\rho_{a}(t_{\rm osc})\sim\theta_{0}^{2}f_{a}^{2}H^{2}(t_{\rm osc})$, $\rho_{a}(T_{\rm RH})$ entering Eq. (8) is independent of $H(t_{\rm osc})$ and $m_{a}$. The present-day ALP density is found to be

[TABLE]

for temperature-independent ALP mass (see Appendix A for more details). In this case $\Omega_{a}$ is determined by $f_{a}$ and $T_{\rm RH}$, with larger $f_{a}$ and $T_{\rm RH}$ corresponding to larger $\Omega_{a}$. Note that this expression is only valid if $T_{\rm RH}$ is larger than $T_{\rm osc}$; otherwise, the standard radiation-dominated scenario is obtained. An expression for $T_{\rm osc}$ is given by Eq. (43) below, where a temperature-independent ALP mass corresponds to $b\rightarrow 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_{\rm RH}$ we have checked.

II.3 Kination

The final cosmological scenario we consider is known as kination Joyce (1997); Ferreira and Joyce (1998). ALP physics with an early period of kination has previously been studied in Visinelli and Gondolo (2010); Visinelli and Redondo (2018). As for the EMD case, the energy density at early times is dominated by a long-lived scalar field $\phi$, but rolling in a steep potential such that its kinetic energy dominates $\rho_{\mathrm{tot}}$. For a polynomial potential $V(\phi)\propto\phi^{N}$, the energy density after some early time $t_{0}$ evolves as

[TABLE]

In the limit $N\rightarrow\infty$, $\rho_{\phi}$ dilutes as $\sim R^{-6}$. We will consider this large $N$ limit in what follows. Assuming that radiation comes to dominate when it reaches a temperature $T_{\rm kin}$, using Eq. (8) with $T_{\rm ad}=T_{\rm kin}$ we find that

[TABLE]

The relic density depends linearly on $m_{a}$ and inversely on $T_{\rm kin}$, which must be larger than about 5 MeV. The expression above only applies if $T_{\rm osc}>T_{\rm kin}$; otherwise, one reproduces the standard RD scenario given by Eq. (11). An expression for $T_{\rm osc}$ in the kination case is given in Eq. (46), where $b\rightarrow 0$ corresponds to a temperature-independent ALP mass. Note that the fractional density $\rho_{a}/\rho_{\mathrm{tot}}\propto 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 $\mathcal{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}\theta_{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. Since – in order to saturate the observed DM relic density – ALP-radiation equality must occur around $T\simeq 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}\theta_{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 Barr and Seckel (1992); Holman et al. (1992a, b); Kamionkowski and March-Russell (1992a, b). 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\sim 12$. Solutions to this problem are known; it might be the case that all quantum gravity-induced PQ-violation is exponentially small Svrcek and Witten (2006); Arvanitaki et al. (2010). 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 non-trivially 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.

III.1 ALP mass from UV physics

We consider the typical contribution to the ALP mass from a dimension-$d$ operator,

[TABLE]

parametrizing $\Phi=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

[TABLE]

We relate the scale $f_{a}$ to the ALP-photon coupling $g_{a\gamma\gamma}$ by assuming Agrawal et al. (2018a)

[TABLE]

where $r$ is an anomaly coefficient that we expect to be ${\cal O}(1)$.333Ref. Agrawal et al. (2018a) constructs models with $r\gg 1$, leading to a large enhancement of $g_{a\gamma\gamma}$. Combining Eqs. (15) and (17) we obtain

[TABLE]

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.

III.2 $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_{\rm HS}$, as

[TABLE]

In what follows, all temperatures will correspond to temperatures of the SM photon bath, unless otherwise stated, and factors of $\xi$ will be used to convert to hidden sector temperatures.

We assume that $m_{a}$ is primarily sensitive to the temperature above a scale $\Lambda$, corresponding to a Standard Model (SM) bath temperature $T_{\Lambda}$. The ALP zero mode is initially frozen at $f_{a}\theta_{0}$ and starts to oscillate when $T=T_{\mathrm{osc}}$. In order for $T$-dependence to have an effect on $\Omega_{a}$, we require $\xi_{\rm osc}\,T_{\rm osc}>\Lambda$ where $\xi_{\rm osc}\equiv\xi(T_{\rm osc})$.

The scale $\Lambda$ 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_{\Lambda}\gtrsim T_{\rm BBN}$, where $T_{\rm BBN}\sim 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_{\Lambda}$ 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_{\rm eff}$. For the parameter space of interest, $T_{\Lambda}$ is always above the temperature of recombination, so the BBN limit is most relevant. These constraints, detailed in e.g. Ref. Aghanim et al. (2018), can be satisfied at $\sim 2\sigma$ confidence level provided $\Delta N_{\rm eff}\lesssim 0.5$. In terms of the effective number of relativistic degrees of freedom $g_{*{\rm HS}}(T_{\rm BBN})$ in the HS at $T_{\rm BBN}$, we have

[TABLE]

where we use the shorthand $\xi_{\rm BBN}\equiv\xi(T_{\rm BBN})$. From this we see that BBN $N_{\rm eff}$ constraints can be avoided if $\xi_{\rm BBN}\lesssim 0.1$ for $g_{*{\rm HS}}(T_{\rm BBN})\lesssim 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 $\xi$ 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 $\xi(T_{\rm 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. Feng and Shadmi (2011); Cline et al. (2014); Boddy et al. (2014); Hochberg et al. (2014); Forestell et al. (2018); Berlin et al. (2018); Draper et al. (2018).

Summarizing, in order for $T$-dependence to affect the ALP relic density and be in agreement with $N_{\rm eff}$ constraints, we require either

[TABLE]

or

[TABLE]

If $\xi_{\rm osc}T_{\rm osc}<\Lambda$, the onset of ALP oscillations proceeds as in the $T$-independent case discussed earlier. The relationship between $T_{\rm 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 $\Omega_{a}$ relative to the $T$-independent prediction for a given $m_{a}$, $f_{a}$, and $\theta_{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 $\mathcal{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

[TABLE]

where $\Omega_{a}$ is the ALP relic density assuming a $T$-dependent ALP mass and $\Omega_{a}^{T-{\rm ind}}$ 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}\theta_{0}$. We model the different cosmological scenarios by assuming that the Hubble parameter for temperatures above some scale $T_{*}$ evolves as

[TABLE]

where $w$ is the equation of state parameter, $\rho=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

[TABLE]

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

[TABLE]

The subscript in $\gamma_{T}$ is indicates that this expression applies if the mass at the onset of oscillations, $m_{\rm osc}$, differs from $m_{a}$, the low-temperature ALP mass. If $m_{\rm osc}\ll m_{a}$, the relic density can be significantly enhanced in the RD and EMD cosmologies. The scaling with $m_{a}/m_{\rm 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. 66 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}\neq 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 $\mathcal{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_{\rm HS}=\xi(T)T>\Lambda$, the ALP mass is given by

[TABLE]

where $m_{a}$ is the zero-temperature mass, taken to be of the form

[TABLE]

For $\xi(T)T<\Lambda$, $m_{a}(T)=m_{a}$. In Eq. (27), $b$ is a positive exponent. In QCD-like theories, $b$ is related to the $\beta$-function of the gauge group and can be obtained analytically from the dilute instanton gas approximation (DIGA). DIGA predicts $b=\left(11N_{c}-2N_{f}\right)/6+N_{f}/2-2$, where $N_{c}$ and $N_{f}$ are the number of colors and light flavors, respectively Gross et al. (1981). For QCD, $b=4$, and the semiclassical approximation is in reasonable agreement with lattice results at high temperatures Borsanyi et al. (2016a, b); Dine et al. (2017). In these simulations the scaling predicted by DIGA appears to hold down to $\Lambda_{QCD}$, where $m_{a}$ saturates to near its zero-temperature value and remains approximately constant at lower temperatures Borsanyi et al. (2016a, b). In this sense our model of the temperature dependence mimics QCD and generalizes it to arbitrary $\Lambda$, $b$, and $\xi(T)$. In our plots we will take $b=4$ as an illustrative example, corresponding to the QCD-like case. As such, we assume

[TABLE]

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_{*{\rm SM}}(T)+\xi(T)^{4}g_{*{\rm HS}}(T)$. Note that for the $T$-independent predictions we take $g_{*{\rm 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, $\gamma_{T}\leq\gamma_{\rm max}$. Defining $T_{\Lambda}$ such that $T_{\Lambda}=\Lambda/\xi(T_{\Lambda})$ (the SM temperature at which the ALP mass saturates to its low-temperature value), if the ALP has not started oscillating by $T_{\Lambda}$ and $m_{a}>q_{0}H(T_{\Lambda})$, oscillations will begin suddenly at $T_{\Lambda}$ and so $q_{\rm max}=m_{a}/H(T_{\Lambda})$. Thus,

[TABLE]

In this case the oscillation temperature is simply $T_{\rm osc}=T_{\Lambda}=\Lambda/\xi_{\rm osc}$. Note that the opposite limit in which the ALP is still frozen at $T_{\Lambda}$ and $m_{a}<q_{0}H(T_{\Lambda})$ corresponds to $\xi_{\rm osc}T_{\rm osc}<\Lambda$, which reproduces the $T$-independent case. One can show that this occurs when $\operatorname{min}(\gamma_{T},\gamma_{\rm max})<1$.

Summarizing these considerations, the enhancement factor can be written compactly as

[TABLE]

where $\gamma_{T}$ and $\gamma_{\rm 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_{\rm osc}(b)$, the true oscillation temperature is given by

[TABLE]

Further details can be found in Appendix A. Finally, in the EMD and kination cosmologies, if $T_{\rm osc}$ predicted by Eq. (35) is smaller than $T_{\rm RH}$ or $T_{\rm kin}$, the results for RD should be used.

III.3 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. Arias et al. (2012), which we generalize to allow for a HS at a different temperature than the SM. In general, $\xi(T)$ changes across mass thresholds as particles in both the HS and SM annihilate. This heating typically changes $\xi$ by at most $\mathcal{O}(1)$ factors unless the change in number of degrees of freedom is very large. Since the precise form of $\xi(T)$ is model-dependent, in the remainder of this study we assume for simplicity that $\xi(T)\approx\xi_{\rm osc}$ for temperatures of interest and treat $\xi_{\rm osc}$ as a free parameter. In concrete models $\xi_{\rm osc}$ can be set by, e.g., the branching ratio of the inflaton into the HS relative to the SM Adshead et al. (2016), 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_{\rm osc}(b)$ and enhancement in this case can be estimated by

[TABLE]

[TABLE]

Meanwhile, the maximum enhancement factor $\gamma_{\rm max}$ is approximately

[TABLE]

With these expressions, Eqs. (11), (31) and (35) can then be used to estimate $\Omega_{a}$ across the ALP plane accounting for $T$-dependence in the ALP mass. In the parameter space we consider, Eq. (31) typically yields $\gamma=\gamma_{T}$ for RD.

The preferred ALP dark matter region in the $T$-dependent case for $\xi_{\rm osc}=1$ is shown on the left in Fig. 2 for $b=4$ and $g_{*\mathrm{HS}}(T)$ given by Eq. (29) with $\xi(T)=\xi_{\rm osc}$. The region shaded gold features an ALP with $h^{2}\Omega_{a}=0.12$ for natural values of the initial misalignment angle, $\theta_{0}\in[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 $T$-independent preferred ALP DM region, obtained by numerically solving the ALP EOM to late times, lies between the dashed gold contours. $\xi_{\rm 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_{\rm eff}$ for $\xi_{\rm osc}=1$, and the shaded gray region on the left in Fig. 2 is excluded by requiring $\Delta N_{\rm eff}<0.5$ (corresponding to $T_{\Lambda}>T_{\rm BBN}\sim 5$ MeV). It is likely that in concrete models the lower bound on $T_{\Lambda}$ 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 $\xi_{\rm osc}=0.1$. The enhancement of $h^{2}\Omega_{a}$ is smaller, however the cooler HS in principle allows for $T_{\Lambda}<T_{\rm BBN}$, since $\xi_{\rm 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-$\Lambda$ scenarios, as well as those with decoupled hidden sectors with dark radiation or significant late-time abundances of stable relics.

Larger $\xi_{\rm osc}$ and $b$ can in principle increase $\Omega_{a}$ further, however this often comes at the cost of additional entropy injection after oscillation in simple models. For example, it could be that $\xi_{\rm osc}>1$, however the large corresponding amount of HS entropy needs to be transferred to the SM before BBN, erasing the resulting enhancement for $\xi_{\rm osc}$ much larger than 1. A hidden sector predicting $b>4$ could also increase $h^{2}\Omega_{a}$ somewhat, however as $b$ increases one also expects $g_{*{\rm HS}}(T_{\rm 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 $\gamma_{T}$ in Eq. (37).

III.4 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 $\phi$ (the field responsible for EMD), SM and HS radiation can be modeled by

[TABLE]

where the dot indicates a derivative with respect to time, $\Gamma_{\phi}=\Gamma_{\phi\to{\rm SM}}+\Gamma_{\phi\to{\rm HS}}$ and $\Gamma_{\phi\to{\rm SM}}$, $\Gamma_{\phi\to{\rm HS}}$ are the partial widths of $\phi$ into to SM and HS radiation, respectively. These equations can be solved during $\phi$ domination (i.e. while $\Gamma_{\phi}/H\ll 1$ and $\rho_{\phi}R^{3}\approx\mathrm{const}$) and yield

[TABLE]

where we assumed that the initial energy densities are negligible compared to those produced by $\phi$ decays. Here $i$ corresponds to either HS or SM subscripts, and $H_{\rm MD}$ and $R_{\rm MD}$ are the Hubble parameter and FRW scale factor at the onset of EMD, respectively. We can use Eq. (42) 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 $\rho_{\phi}=\rho_{\rm SM}+\rho_{\rm HS}$ and assuming Eq. (42) holds down to that temperature allows one to relate $\Gamma_{\rm SM}$ and $\Gamma_{\rm HS}$ to $T_{\rm RH}$ and $\xi_{\rm osc}$ (again neglecting relative heating effects). We then find that the oscillation temperature and abundance enhancement factor can be approximated by

[TABLE]

[TABLE]

(see Appendix A for more details). $\gamma_{\rm max}$ in EMD is given by

[TABLE]

These expressions can be inserted into Eqs. (31) – (35) and used along with Eq. (12) to estimate $\Omega_{a}$ in the $T$-dependent case. Again if $T_{\rm osc}<T_{\rm RH}$, the RD expressions should be used. The resulting predictions agree well with full numerical solutions of the three-fluid system of equations, discussed in Appendix B.

The preferred regions of the ALP parameter space in the EMD scenario are illustrated in Fig. 2 for $\xi_{\rm osc}=1$ and $\xi_{\rm osc}=0.1$ with $b=4$. We show results for $T_{\rm RH}=10$ and 500 MeV. The blue and purple shaded regions feature an ALP with $\Omega_{a}h^{2}\simeq 0.12$ for $\theta_{0}\in[0.1,2]$ as obtained from the numerical solution. The corresponding dotted contours show the analytic predictions of Eqs. (43) – (45) and are a good fit to the numerical results. For $\xi_{\rm osc}=1$, the gray shaded region is excluded by the measured value of $N_{\rm eff}$ at BBN. This constraint is alleviated for $\xi_{\rm osc}=0.1$, however the enhancement factor $\gamma$ is reduced as a result. Other model-dependent constraints are likely to apply in the region where $T_{\Lambda}<T_{\rm BBN}$ as discussed in Sec. III.2.

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_{\Lambda}$, temperature-dependent effects are unimportant and the preferred regions are the same as discussed in Sec. II.2 (i.e. $\gamma=1$). This occurs for larger ALP masses, and the preferred value of $g_{a\gamma\gamma}$ is independent of $m_{a}$. For small enough masses, $T_{\rm osc}>T_{\Lambda}$ 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\gamma\gamma}$ in Fig. 2. At even lower values of $m_{a}$, oscillation occurs after reheating ($T_{\rm osc}<T_{\rm RH}$), and the predictions reduce to those of the $T$-dependent RD scenario of Sec. II.1 (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_{\rm 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 $T$-dependent RD and $T$-independent EMD contours for a given $\theta_{0}$.

III.5 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

[TABLE]

[TABLE]

Meanwhile,

[TABLE]

In most realistic models, one expects $b\sim\mathcal{O}(1)$, and so typically $\gamma=\gamma_{T}$. The enhancement is milder than in RD and EMD as it only depends on the exponent $b$. Nevertheless, allowing for $m_{\rm osc}\neq m_{a}$ changes the ALP oscillation temperature. Since $T_{\rm osc}>T_{\rm 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_{\rm kin}$.

In our results below we will take $T_{\rm kin}=10$ MeV. Since the preferred regions on the ALP parameter space assuming kination with and without $T$-dependence are similar for $T_{\rm 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 temperature-independent 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_{\mathrm{RH}}=10\,\mathrm{MeV}$, and kination with $T_{\mathrm{kin}}=10\,\mathrm{MeV}$. $T_{\mathrm{RH}}$ and $T_{\mathrm{kin}}$ are the temperatures at which the universe transitions to standard radiation-dominated evolution; temperatures of $5-10\,\mathrm{MeV}$ correspond to the lowest values compatible with BBN. In each case the bands are obtained by varying the initial misalignment angle $\theta_{0}$ between $0.1$ (bottom edge of each band) and $2$ (upper edge of each band). The ALP regions can be extended to smaller values of $g_{a\gamma\gamma}$ at the cost of fine-tuning $\theta_{0}<0.1$. The gray dotted lines in Figs. 3, 4 and 5 show the ALP mass-coupling relation if the masses are generated by Planck-suppressed operators of various dimensions as discussed in Sec II. For the temperature-dependent results in Fig. 5 we have assumed $b=4$ in Eq. (27) and taken $g_{*{\rm HS}}$ as in Eq. (29) with $\xi=\xi_{\rm 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\times T_{\rm{SM}}$.

Axion-like particles can be constrained by a variety of astrophysical measurements. These limits include the results from CAST Anastassopoulos et al. (2017); cooling of Horizontal Branch (“HB” in Fig. 3) stars, massive stars Cadamuro and Redondo (2012); Friedland et al. (2013), and SN1987A Raffelt (1996); Dolan et al. (2017); Lee (2018); non-observation of a $\gamma$-ray excess from SN1987A Payez et al. (2015); the extragalactic background light Masso and Toldra (1997); Overduin and Wesson (2004); searches for spectral irregularities in $\gamma$ rays with HESS Abramowski et al. (2013) and Fermi-LAT Ajello et al. (2016), and in X-rays with Chandra Chen and Conlon (2018).444ALPs can also be constrained by observations of near-extremal black holes and the resulting constraints on superradiance Arvanitaki et al. (2015). However, these constraints are strongly model-dependent in that 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 $m_{a}\lesssim 10^{-11}$ eV. These limits are shown in Fig. 3 as pastel-colored shaded regions.

The ALP parameter space is also constrained by a number of resonant cavity experiments. We show the regions excluded by ADMX Asztalos et al. (2010); Du et al. (2018) and ADMX Sidecar Boutan et al. (2018), Phase 1 of HAYSTAC Zhong et al. (2018), the ORGAN Pathfinder McAllister et al. (2017) and the older UF Hagmann et al. (1990) and RBF De Panfilis et al. (1987); Wuensch et al. (1989) 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 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 Shokair et al. (2014), CAPP Petrakou (2017), KLASH Alesini et al. (2017); Gatti et al. (2018), and, at higher frequencies, ORGAN McAllister et al. (2017). These experiments provide a broader sensitivity in the QCD axion region $m_{a}\sim 10^{-6}-10^{-5}$ eV extending to lower values of $g_{a\gamma\gamma}$. More recently, new ideas based on dielectric stacks have appeared which are sensitive to higher mass ALPs, as in MADMAX Caldwell et al. (2017) and photonic materials Baryakhtar et al. (2018) (“Dielectric Stack” in Figs. 4 and 5). We also show the sensitivity of the proposal for a large-scale helioscope, IAXO Irastorza et al. (2013), which will extend the reach of CAST, and the projections for the ALPS-II light-shining-through-walls experiment Bähre et al. (2013), which is currently under construction at DESY and will have sensitivity above $g_{a\gamma\gamma}\sim 10^{-11}\leavevmode\nobreak\ \rm{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 $\mu$eV mass-range Hook et al. (2018). 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 Creque-Sarbinowski and Kamionkowski (2018) 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\gamma\gamma}$555We also note that there is the BEAST proposal McAllister et al. (2018) which could be relevant at low masses. However, since the BEAST projections are a topic of current discussion in the literature Ouellet and Bogorad (2018); Beutter et al. (2018); Kim et al. (2018) we do not show them on our plot. We also do not show other limits from other as-yet-unpublished proposals, such as Goryachev et al. (2018); Marsh et al. (2018); Bogorad et al. (2019); Janish et al. (2019); Edwards et al. (2019).. Other recent proposals at low mass make use of birefringence in the presence of an ALP background and include the interferometer concept DeRocco and Hook (2018), ADBC Liu et al. (2018), and an experiment based on optical ring cavities Obata et al. (2018).

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}\sim 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 $\sim 4\times 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 matter with mass $10^{-12}\leq m_{a}\leq 1$ eV in a variety of well-motivated 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 $\Omega_{a}=\gamma\Omega_{a}^{T-{\rm ind}}$, with the enhancement factor $\gamma$ given by Eq. (31) and Eqs. (37) – (38) for RD, Eqs. (44) – (45) for EMD, and Eqs. (47) – (48) for kination, respectively.

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 PQ-breaking 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 ABRACADABRA 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}\lesssim 4\times 10^{-7}$ eV down to $10^{-12}$ eV. Below $10^{-11}$ eV black-hole super-radiance complements the ABRACADABRA and DM-Radio sensitivity, although the precise details are model-dependent. 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}\gtrsim 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_{\rm kin}$) and come close to the standard ALP prediction with $\mathcal{O}(1)$ initial misalignment angles and $T$-dependent 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 Blinov et al. ; Nelson and Xiao (2018); Visinelli and Redondo (2018). 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.

Acknowledgements We thank Aaron Chou, Yonatan Kahn, Jeff Filippini, and Manuel Meyer for useful discussions. MJD is supported by the Australian Research Council. The work of PD and JK was supported by NSF grant PHY-1719642. This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.

Appendix A Analytic Estimates for the Relic Density

In this Appendix we derive analytic estimates of the ALP DM number density.

A.1 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:

[TABLE]

$T_{ad}$ is a reference temperature below which the evolution of the universe is adiabatic and $T_{0}$ is the temperature today. $T_{\rm osc}$ is the SM temperature at the onset of ALP oscillations, and we have $\dot{\theta}(t_{\rm osc})\simeq 0$ and $n_{a}(T_{\rm osc})=1/2m_{a}f_{a}^{2}\theta_{0}^{2}$. Below $T_{\rm osc}$ the comoving ALP number density is assumed to be conserved, so that

[TABLE]

Combining Eqs. (49) – (50) yields a general expression for the relic density,

[TABLE]

To obtain $\Omega_{a}$ we therefore need to determine $T_{\rm osc}/T_{\rm ad}$ and/or $R_{\rm osc}/R_{\rm ad}$ in the various cases.

We parametrize the relationship at oscillation between the Hubble parameter and ALP mass as

[TABLE]

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\rho$ (for which $H^{2}\propto R^{-3(w+1)}$), and then later transitions instantaneously to radiation domination at some $T=T_{*}$:

[TABLE]

In all cases we consider (RD, EMD, and kination), we can take

[TABLE]

Henceforth we will refer only to $T_{*}$. In RD, we can also set $T_{*}=T_{\rm osc}$. In EMD, $T_{*}=T_{\rm RH}$, and in kination $T_{*}=T_{\rm kin}$. Combining Eqs. (52) – (53), we obtain $R_{\rm osc}/R_{\rm*}$ in terms of $T_{*}$ and $m_{a}$. Plugging in to Eq. (51) 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_{\rm ad}=T_{\rm osc}$ and used the conservation of comoving entropy to relate $T_{\rm osc}$ to $T_{\rm 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\rho$, such that $R\sim t^{\frac{2}{3(w+1)}}$. With $\theta(t)\equiv a(t)/f_{a}$ we have

[TABLE]

which has solutions

[TABLE]

Here $r=\frac{1}{2}-\frac{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 $\mathcal{O}(1)$ or larger, we see that ALP oscillations begin when

[TABLE]

with $A$ an $\mathcal{O}(1)$ number. Since $H\simeq 2/(3t(1+w)$, we define the onset of oscillations as

[TABLE]

For EMD ($w=0$), RD ($w=1/3$), and kination ($w=1$), we obtain $m_{a}\simeq\left\{3/2,2,3\right\}\times AH_{\rm osc}$, respectively. Comparing our analytic and numerical solutions, we find that choosing $A$ such that $m_{a}\simeq 1.6H_{\rm 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. Marsh (2016), which found $m_{a}\sim 2H_{\rm osc}$ is a better choice than $m_{a}\sim 3H_{\rm 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 (58) 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 (58),

[TABLE]

A.2 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_{\rm HS}$ as

[TABLE]

where $T$ is the temperature of the SM photon bath. We take an instanton-motivated class of models in which $m_{a}(T_{HS})=(\Lambda/T_{HS})^{b}m_{a}$ for $T>\Lambda$. Subsequently, all temperatures will correspond to SM temperatures, unless otherwise stated, and factors of $\xi$ will be used to convert to hidden sector temperatures.

Denoting the ALP mass at the onset of oscillations as $m_{\rm osc}$ and the mass today as $m_{a}$, the relic density can be expressed as:

[TABLE]

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\rho$ (for which $H^{2}\propto R^{-3(w+1)}$), and then later transitions instantaneously to radiation domination at some $T=T_{*}$ so that Eq. (53) applies. We parametrize the oscillation time in the case of temperature-dependent ALP masses via

[TABLE]

generalizing Eq. (52). $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.

For a fixed $m_{a}$, $f_{a}$ and $\theta_{0}$, we can compare the relic densities in the $T$-dependent and $T$-independent cases. From Eqs. (51), (53), and (64) we have

[TABLE]

Thus, to determine $\Omega_{a}$ in the $T$-dependent case, we can use the $T$-independent results of Eqs. (11), (12), and (14), multiplying by $\gamma$ (assuming that $T$-dependence is relevant at the time of oscillations) and using the value of $T_{\rm osc}$ predicted for $T$-dependent masses. To do so, we must determine $m_{\rm 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}(\Lambda/\xi(T)T)^{b}$. During a period of adiabatic evolution dominated by a fluid with equation of state $\rho=wp$, the temperature evolves as

[TABLE]

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

[TABLE]

where $T_{i}$ is a constant (neglecting the $T$-dependence of $g_{*}$ and $\xi$). 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

[TABLE]

The solution starts to oscillate when the above quantity $\sim\mathcal{O}(1)\equiv A$. Using $H=2/(3(1+w)t)$ and $m_{a}(t)=m_{a}t^{\frac{2b}{3(1+w)}}(\Lambda/T_{i})^{b}$, oscillations begin when

[TABLE]

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. (42) (again neglecting $T$-dependence in $g_{*}$ and $\xi$). The ALP EOM in this case is approximately

[TABLE]

Proceeding as before, we find that the onset of oscillations occurs when

[TABLE]

Eqs. (70) and (72) reduce to the results from Eq. (58) when $b=0$, as they should. Taking the same values for $A$ as in Eq. (62), we obtain

[TABLE]

When comparing to our numerical solutions, we find that neglecting the $b$-dependence in $q_{T}$ for $b=4$ typically results in $\sim 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_{\rm osc}=q_{T}H_{\rm osc}$, assuming that $m_{a}(T)=m_{a}(\Lambda/\xi T)^{b}$ at oscillation. We refer to this oscillation temperature as $T_{\rm osc}(b)$ to distinguish it from the $T$-independent prediction. To obtain $T_{\rm 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_{\rm kin}$ is specified, and the scale factor evolution from $T_{\rm kin}$ to $T_{\rm osc}$ is determined by Eq. (53). For EMD, one can use the late-time solution for the radiation energy density, Eq. (42), along with Eq. (53), to write the Hubble parameter in terms of the temperature. Solving for $T_{\rm osc}(b)$ in this way yields Eqs. (36), (43), and (46). Inserting these temperatures into $m_{a}(T)$ and Eq. (66) yields Eqs. (37), (44) and (47) for $\gamma_{T}$.

This is not quite the end of the story for the temperature-dependent scenarios. Our results for $q_{T}$, $T_{\rm osc}(b)$, and $\gamma_{T}$ followed from inspecting solutions to the ALP EOM with $m_{a}(T)=m_{a}(\Lambda/\xi T)^{b}$. Therefore, they only apply if $q_{T}H_{\Lambda}\leq m_{a}$, where $H_{\Lambda}$ is the Hubble parameter at $T_{\Lambda}=\Lambda/\xi(T_{\Lambda})$, when the mass saturates to its $T$-independent value. If instead $q_{T}H_{\Lambda}>m_{a}$, oscillations begin after the mass has already saturated to $m_{a}$, and the ALP EOM solution is given by Eq. (56). Then there are two possibilities we must consider. First, if $q_{0}H_{\Lambda}>m_{a}$, we recover the $T$-independent case. Second, if $q_{0}H_{\Lambda}<m_{a}$, at $T_{\Lambda}$ the ALP will start to rapidly oscillate across this threshold, and $H_{\rm osc}$ is given by

[TABLE]

In this sudden-oscillations case, $T_{\rm osc}=T_{\Lambda}=\Lambda/\xi_{\rm osc}$ and $\gamma$ is given by

[TABLE]

which yields Eqs. (38), (45), and (48).

Summarizing, we have derived the enhancement factor relevant for $T$-dependent ALP masses of the form $m_{a}(T)=m_{a}(\Lambda/\xi T)^{b}$ for three qualitatively different regimes. If $q_{T}H_{\Lambda}\leq m_{a}$, $\gamma=\gamma_{T}$. If instead, $q_{T}H_{\Lambda}>m_{a}>q_{0}H_{\Lambda}$, oscillations begin suddenly at $T_{\Lambda}$, and $\gamma$ should instead be taken as $\gamma_{\rm max}$. The condition $q_{T}H_{\Lambda}>m_{a}>q_{0}H_{\Lambda}$ corresponds to $\gamma_{T}>\gamma_{\rm max}$. Finally, if $q_{T}H_{\Lambda},\,q_{0}H_{\Lambda}>m_{a}$, one reproduces the $T$-independent case, and $\gamma=1$. This case occurs when both $\gamma_{T},\gamma_{\rm max}<1$. Therefore, the correct value of $\gamma$ is given compactly by

[TABLE]

The temperature $T_{\rm osc}$ should be specified as in Eq. (35). If the resulting $T_{\rm osc}<T_{\rm RH}$ or $T_{\rm kin}$, the results for radiation domination should be used. Using the above expressions typically reproduces the numerical results to within a few tens of percent across the parameter space.

Appendix B Numerics

To fit coefficients in some of the analytical expressions, and to determine their accuracy, we solve the ALP EOM

[TABLE]

numerically for the EMD, RD, and kination scenarios. We consider times between $t_{0}$ and $t_{f}$ with $\theta(t_{0})=\theta_{0}$. The initial time $t_{0}$ is chosen sufficiently above $t_{\rm osc}$ (defined implicitly by Eq. (35)) to capture the start of ALP evolution, while $t_{f}$ is taken large enough so that we can subsequently approximate the ALP number density as redshifting like nonrelativistic matter. To obtain the relic density today, we use

[TABLE]

where $n_{a}(t_{f})$ is the ALP number density at the time $t_{f}$:

[TABLE]

The various cosmologies predict different relations for $H(t)$ and $m_{a}(t)$ entering Eq. (80). The ALP contribution to $H$ is negligible at times near oscillation in the parameter space of interest. We therefore neglect it and model the matter, radiation, and kination components as perfect fluids with equation of state parameters $w=0,\frac{1}{3}$, and $1$, respectively. Furthermore, we treat the SM and HS radiation (if any) as a single fluid with temperatures related by $\xi(T)$, which is assumed to be approximately constant and set by $\xi_{\rm osc}$, taken as an input parameter. In the standard RD cosmology, we can simply use

[TABLE]

to determine $T(t)$. In the kination scenario, we have

[TABLE]

where $r\equiv R/R_{\rm kin}$ is the ratio of the FRW scale factor to its value at $T_{\rm kin}$, defined as the temperature for which the radiation and kination field energy densities are equal. The parameter $r$ is related to $t$ through

[TABLE]

where $t_{0}$ should be taken before ALP oscillation and during the kination phase so that $r_{0}$ can be defined via $H(r_{0})=1/(3t_{0})$. The resulting $r(t)$ yields $H(t)$ through Eq. (84), while $T(t)$ is defined from

[TABLE]

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 system

[TABLE]

with $H=1/(\sqrt{3}M_{\rm Pl})\sqrt{\rho_{\phi}+\rho_{\rm R}}$ and $\Gamma_{\phi}$ chosen so that $\rho_{R}=\rho_{\phi}$ at a temperature $T_{\rm RH}$. The temperature is defined from the radiation energy density $\rho_{R}=\pi^{2}/30g_{*}(T)T^{4}$ in the usual way. The final time is chosen to correspond to temperatures below $T_{\rm RH}$ so that Eq. (81) 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 $\rho=1/3p$. Conservation of comoving entropy density can instead be used in this case to track the evolution of $\rho_{\rm R}$ with $t$. However, we find that the corresponding effects are rather small and so we neglect them in our analysis.