###### Abstract

We study axion dark matter production from a misalignment mechanism in
scenarios featuring a general nonstandard cosmology. Before the onset of Big
Bang nucleosynthesis, the energy density of the universe is dominated by a
particle field $\phi$ described by a general equation of state $\omega$. The
ensuing enhancement of the Hubble expansion rate decreases the temperature at
which axions start to oscillate, opening this way the possibility for axions
heavier than in the standard window. This is the case for kination, or in
general for scenarios with $\omega>1/3$. However, if $\omega<1/3$, as in the
case of an early matter domination, the decay of $\phi$ injects additional
entropy relative to the case of the standard model, diluting this way the
preexisting axion abundance, and rendering lighter axions viable. For a
misalignment angle $0.5<\theta_{i}<\pi/\sqrt{3}$, the usual axion window
becomes expanded to $4\times 10^{-9}$ eV $\lesssim m_{a}\lesssim 2\times
10^{-5}$ eV for the case of an early matter domination, or to $2\times
10^{-6}$ eV $\lesssim m_{a}\lesssim 10^{-2}$ eV for the case of kination.
Interestingly, the coupling axion-photon in such a wider range can be probed
with next generation experiments such as ABRACADABRA, KLASH, ADMX, MADMAX, and
ORGAN. Axion dark matter searches may therefore provide a unique tool to probe
the history of the universe before Big Bang nucleosynthesis.

New opportunities for axion dark matter searches
in nonstandard cosmological models

Paola Ariasa, Nicolás Bernalb, Dimitrios Karamitrosc,
Carlos Maldonadoa, Leszek Roszkowskic,d and Moira Venegasa
a Departamento de Física, Universidad de Santiago de Chile, Casilla 307,
Santiago, Chile
b Centro de Investigaciones, Universidad Antonio Nariño
Carrera 3 Este # 47A-15, Bogotá, Colombia
c National Centre for Nuclear Research, ul. Pasteura 7, 02-093 Warsaw, Poland
d AstroCeNT, Nicolaus Copernicus Astronomical Center Polish Academy of
Sciences,
ul. Rektorska 4, 00-614 Warsaw, Poland

PI/UAN-2021-693FT

###### Contents

  1. 1 Introduction
  2. 2 Axion DM in the standard cosmological scenario: misalignment essentials
  3. 3 Nonstandard cosmologies
  4. 4 Axion oscillation in NSC
    1. 4.1 Early matter domination and $\beta<4$
      1. 4.1.1 Region 1: $R_{\text{osc}}\ll R_{\text{eq}}$
      2. 4.1.2 Region 2: $R_{\text{eq}}\ll R_{\text{osc}}\ll R_{\text{c}}$
      3. 4.1.3 Region 3: $R_{\text{c}}\ll R_{\text{osc}}\ll R_{\text{end}}$
    2. 4.2 Kination and $\beta>4$
  5. 5 Axion coupling to two photons
  6. 6 Conclusions
  7. A Solutions for the $\phi-$radiation system
  8. B Axion oscillations during $\phi$ domination
    1. B.1 Oscillations in Region 1
    2. B.2 Oscillations in Region 2
    3. B.3 Oscillations in Region 3
  9. C Adiabatic invariant and the anharmonic factor

##  1 Introduction

The axion is a by-product of an elegant solution of the strong CP problem of
the QCD sector [1, 2, 3] of the standard model of particle physics (SM).
Interestingly, the same mechanism that solves the strong CP problem, leads to
an efficient mechanism of non-thermal production of a cold population in the
early universe, making the axion an excellent cold dark matter (CDM) candidate
[4, 5, 6]. Axions appear as pseudo-Nambu-Goldstone bosons when the so-called
Peccei-Quinn symmetry (PQS) is spontaneously broken at some energy $f_{a}\gg
v_{\rm EW}\simeq 247$ GeV. As the universe cools down, the axion potential
energy changes during the QCD phase transition epoch, due to instanton effects
that break explicitly the PQS, acquiring a small mass. During this process,
the value of the axion field gets realigned (the process is known as the
“misalignment mechanism”), changing from an arbitrary initial value $a_{i}$,
to the true vacuum value of the field, $\langle a\rangle=0$. This process is
of extreme relevance, since on the one hand, it solves the so-called strong CP
problem of the QCD sector and fills the universe with cold axion particles
produced during the oscillation of the field around the minimum.

Because of the realignment of the axion potential at the QCD scale, the mass
before the QCD phase transition is thermally suppressed and after the
transition it is set to a constant value. A distinctive feature of the axion
is that its mass has a definite relationship with the PQ scale by
$f_{a}\,\tilde{m}_{a}(T)=\sqrt{\chi(T)}$, where $\chi$ is the topological
susceptibility in QCD, which has been evaluated in the chiral limit [7, 8],
next-to-next-to leading order in chiral perturbation theory [9] and directly
via QCD lattice simulations [10], all coinciding with the central limit

$m_{a}\simeq 5.69\,\mbox{meV}\left(\frac{10^{9}\,\mbox{GeV}}{f_{a}}\right)\,,$
(1)

where $m_{a}$ corresponds to the axion mass at zero temperature.

Only two parameters are needed to fix the relic density of axions: the scale
$f_{a}$ (or equivalently its mass) and the initial value of the field at the
moment when it starts to oscillate, $a(T_{\text{osc}})\simeq a_{i}$. In a pre-
inflationary scenario, the PQS breaking occurs before or during inflation and
it is not restored afterwards. In this scenario, the axion field takes
different initial values in different disconnected patches that are eventually
stretched during inflation, homogenising its initial value in our universe to
$a_{i}=f_{a}\,\theta_{i}$, where the initial angle is randomly selected in the
range $\theta_{i}\in\left[-\pi,\pi\right]$. For
$\theta_{i}\sim\mathcal{O}(1)$, $m_{a}$ should be in the range of $\mu$eV in
order for the axion density not to exceed the observed CDM density
$\Omega_{\rm CDM}h^{2}\simeq 0.12$ [11, 12, 10, 13, 14]. Smaller masses
(equivalently, higher $f_{a}$ scales) can be reached if small-tuned values of
$\theta_{i}$ are considered. This is the so-called anthropic axion window [15,
12]. Nonetheless, in this pre-inflationary scenario, high values of $f_{a}$
are constrained by isocurvature perturbations [16, 15, 14], requiring
$f_{a}\lesssim 10^{16}$ GeV. On the other hand, in a post-inflationary
scenario, the PQS breaks after inflation, so $\theta_{i}$ takes different
values in different patches of the present universe, so its value gets
averaged to $\bar{\theta}_{i}=\pi/\sqrt{3}$, which then sets the DM abundance
to be correctly satisfied for $m_{a}\simeq 30~{}\mu$eV. However, in this
scenario topological defects such as cosmic strings and domain walls emerge,
contributing to the CDM abundance in a still disputed quantity due to the
uncertainties in their energy loss process [17, 18, 19, 20, 9]. Thus, in the
post-inflationary PQS breaking the accepted axion mass is in the range
$30~{}\mu$eV $\lesssim m_{a}\lesssim 5$ meV.

The QCD axion has been realised in many scenarios, where the most popular ones
correspond to the DFSZ model [21, 22] and the KSVZ model [23, 24, 25]. The
axion can feature couplings to nucleons, electrons and more importantly to two
photons. Stellar evolution arguments constrain the previously mentioned
interactions, where the most stringent one comes from the coupling of axions
to nucleons, leading to an axion flux from a supernova explosion that
otherwise would quench the observed neutrino pulse from SN1987a, unless the
axion decay constant is constrained to be [26, 27, 28]

$f_{a}\gtrsim 4\times 10^{8}~{}\mbox{GeV},$ (2)

corresponding to a mass of $m_{a}\lesssim 1.6\times 10^{-2}$ eV.

The uncertainty in the astrophysical models and assumptions on axion coupling
to fermions makes laboratory based searches crucial to test axions in the
$\mu$eV mass ballpark. So far, only the so-called haloscope experiments [29],
ADMX [30, 31] and HAYSTAC [32], using a highly tuned microwave cavity that
converts DM axions into photons in the presence of a static magnetic field,
have achieved enough sensitivity to touch the axion band in the sub-eV regime.
Several operating and proposed haloscopes aim to explore masses in between the
$\mu$eV to the meV range, such as CULTASK [33, 34] KLASH [35], ORGAN [36, 37],
RADES [38, 39] and QUAX [40]. A related concept makes use of dielectric
plates, such as the MADMAX [41] experiment. The LC circuit based detection
[42] is expected to reach axion masses in the $10^{-7}$ to $10^{-9}$ eV range.
The ABRACADABRA experiment [43, 44] and ADMX SLIC [45] have already released
promising results. In a higher mass range, the helioscope experiment CAST has
explored the $m_{a}\sim$ eV ballpark and it is expected to improve with the
next generation experiment IAXO [46]. For details on axion searches,
prospects, and their schedules, we refer, e.g., to Refs. [47, 48, 49].

So far, no signal of axion DM has been found, but given the many forthcoming
experiments that will access to probe the axion CDM prediction, it is expected
that in the following years either a discovery is made, or the axion is not
found in that parameter space. In the second pessimistic scenario, there are
several ways to go around the preferred DM parameter space, such that the
axion is either even more weakly or strongly coupled than thought. The most
straightforward approach is to consider that in the case of the pre-
inflationary scenario, there is a very small initial misalignment angle
$\theta_{i}$, to compensate for a higher $f_{a}$ scale shifting it to be as
high as $\sim 10^{16}$ GeV, as allowed by bounds on isocurvature
perturbations. Another way to open up the axion DM window to smaller masses or
higher PQ scales invokes the coupling of the axion to some other field, for
instance, to a hidden photon [50, 51], or if the axion potential was much
larger in the very early universe [52, 53, 54, 55]. A different possibility to
expand the range masses has been to consider a nonstandard cosmological (NSC)
history in the early universe, such that entropy is injected into the thermal
bath, diluting the axion energy density.

In this work we pursue the last alternative, that is, consider a NSC history
of the early universe prior to Big-Bang Nucleosynthesis (BBN) by means of a
new field that eventually dominates the expansion of the universe. There are
several works that have previously studied a NSC scenario in the context of
axion physics. Besides, pioneer papers that considered an early matter
dominance (EMD) period [56, 57, 58, 59], in Ref. [60, 61] the authors studied
the cases of low temperature reheating (LTRH) and kination, including the
anharmonicities in the axion potential. In Refs. [62, 63] thermally produced
axions in LTRH and kination cosmologies were analysed, whereas in Ref. [64] it
was considered the misalignment production of axion-like particles in early
matter domination and kination cosmologies. In Ref. [65] a full scan of
cosmological histories was performed, together with including the contribution
from the decay of an axionic string network. In Refs. [66, 67] the impact of
NSC on the formation of axion miniclusters was analysed.

Our approach is to consider a very detailed analysis of the DM production
during a NSC scenario. In particular, we study analytically the axion relic
density and we find the mass that can account for the whole DM abundance
observed today as a function of the NSC parameters. We also solve numerically
to check our results, getting a good agreement. We find the features and
characteristics of the NSCs where the misalignment mechanism can lead to an
extended range in the axion mass and we map their impact on the coupling to
two-photon plane. We point the experimental efforts and prospects that can
test those cosmological scenarios.

The manuscript is organised as follows: in the next section 2 we will review
the standard misalignment mechanism for axions in the framework of a standard
cosmological scenario. We will derive the analytical expressions of the relic
density and compare them with the numerical results. Then in section 3 we
introduce the characteristics and parameters that define the models of NSC
considered. In section 4 we derive the axion relic density in terms of the NSC
parameters and then we contrast with the numerical results. We divide our
analysis, starting with cosmologies with an equation of state $\omega<1/3$,
where we choose as benchmark values $\omega=0$ and $\omega=-1/3$. Then, we
study the case $\omega>1/3$, using as benchmarks $\omega=1$ and $\omega=5/3$.
In section 5 we map the results into the exclusion plot for the coupling of
axions to two photons. We show the opportunities to extend the DM parameter
space for different cosmologies, and we mention on the new generation of
experiments that will be able to test some of that parameter space. Finally,
in section 6 we summarise and conclude.

##  2 Axion DM in the standard cosmological scenario: misalignment essentials

Let us start by reviewing the essentials of the misalignment production and
write down the parameter space where the axion can explain the whole DM
observed today. The temporal evolution of the SM temperature $T$ can be
obtained from the evolution of the SM entropy density
$s(T)=\frac{2\pi^{2}}{45}\,g_{\star S}\,T^{3}$ given by

$\frac{ds}{dt}+3\,H\,s=0\,,$ (3)

where $g_{\star S}(T)$ accounts for the number of relativistic degrees of
freedom present in the SM entropy [68], and $H$ corresponds to the Hubble
expansion rate which, in a radiation dominated universe, is given by

$H(T)=\sqrt{\frac{\rho_{R}(T)}{3\,M_{P}^{2}}}=\frac{\pi}{3}\sqrt{\frac{g_{\star}(T)}{10}}\,\frac{T^{2}}{M_{P}}\,,$
(4)

where $\rho_{R}(T)$ is the SM radiation energy density, $g_{\star}(T)$
accounts for the relativistic degrees of freedom contributing to $\rho_{R}$
[68], and $M_{P}\simeq 2.4\times 10^{18}$ GeV is the reduced Planck mass.
Equation (3) allows to extract the temperature evolution in terms of the scale
factor $R$ as

$\frac{dT}{dR}=-\left(1+\frac{T}{3\,g_{\star S}}\frac{dg_{\star
S}}{dT}\right)^{-1}\left(\frac{T}{R}\right).$ (5)

Furthermore, the scale $f_{a}$ at which the PQS is spontaneously broken
determines the axion mass through the topological susceptibility of QCD,
$\chi(T)$, as

$\tilde{m}_{a}^{2}(T)=\frac{\chi(T)}{f_{a}^{2}}\,,$ (6)

where $\chi(T)$ has been estimated from lattice QCD simulations and found a
zero-temperature value of $\chi_{0}\equiv\chi(0)\simeq 0.0245$ fm-4, in the
symmetric isospin case [10]. For the numerical calculations the results of
Ref. [10] will be used, however, for analytical estimations we will take
instead an approximate expression, which has to be cut off by hand once the
mass reaches the zero-temperature value [15]

$\tilde{m}_{a}(T)\simeq m_{a}\times\begin{dcases}1&\text{ for }T\leq T_{\rm
QCD}\,,\\\ \left(\frac{T}{T_{\rm QCD}}\right)^{-4}&\text{ for }T\geq T_{\rm
QCD}\,,\end{dcases}$ (7)

with $m_{a}$ given in Eq. (1), and $T_{\rm QCD}\simeq 150$ MeV the approximate
temperature at the QCD phase transition.

To track the evolution of the axion field during the early universe, let us
write down the axion Lagrangian density as

$\mathcal{L}=\frac{1}{2}\,\partial_{\mu}a\,\partial^{\mu}a-\tilde{m}_{a}^{2}(t)\,f_{a}^{2}\left(1-\cos\frac{a}{f_{a}}\right).$
(8)

For an homogeneous axion field, we only track down the evolution of the zero
mode, which is found to be

$\ddot{\theta}+3\,H(t)\,\dot{\theta}+\tilde{m}_{a}^{2}(t)\,\sin\theta=0\,,$
(9)

where $\theta(t)\equiv a(t)/f_{a}$. At high energies (i.e., $T\gg T_{\rm
QCD}$), the last term can be ignored and the field is stuck on a constant
value $\theta_{i}$, which is in principle randomly selected. At temperature
$T=T_{\text{osc}}$, defined as the temperature where the equality

$3\,H(T_{\text{osc}})=\tilde{m}_{a}(T_{\text{osc}})$ (10)

holds, the axion field starts to oscillate around its true minimum,
$\theta=0$. Within the WKB approximation, i.e. $\theta\ll 1$ and slowly
varying $H$ and $\tilde{m}_{a}$ (adiabatic evolution of the axion),111For
larger values of the angle, one has to include the “anharmonic” contribution
to the potential. In our numerical results, the anharmonic factor derived in
Appendix C are included.

$a(t)=a_{i}\left[\frac{\tilde{m}_{a}(t_{\text{osc}})}{\tilde{m}_{a}(t)}\left(\frac{R_{\text{osc}}}{R(t)}\right)^{3}\right]^{1/2}\cos\left(\int\tilde{m}_{a}(t)\,dt\right).$
(11)

We should note that this we also have assumed that that the value of the field
at $t=t_{\text{osc}}$ is close to its initial value, i.e. $a(t=t_{\rm
osc})\simeq a_{i}$. It is convenient to rewrite Eq. (9) in terms of the scale
factor as

$\theta^{\prime\prime}+\left(\frac{4}{R}+\frac{H^{\prime}(R)}{H(R)}\right)\theta^{\prime}+\left(\frac{\tilde{m}_{a}(R)}{H(R)\,R}\right)^{2}\sin\theta=0\,,$
(12)

where primes denote derivatives with respect to $R$.

The axion energy density $\rho_{a}$ is given by

$\rho_{a}(t)=\frac{\dot{a}^{2}}{2}+\tilde{m}_{a}^{2}(t)\,f_{a}^{2}\left(1-\cos\frac{a}{f_{a}}\right)\simeq\frac{\dot{a}^{2}}{2}+\frac{\tilde{m}_{a}^{2}(t)\,a^{2}}{2}\,,$
(13)

where in the last step $a/f_{a}\ll 1$ was assumed. It is interesting to note
that even if $\rho_{a}$ is not conserved since its mass (and the scale factor)
vary with time, its comoving number $N_{a}\equiv n_{a}\,R^{3}$ is conserved.
The energy density for non-relativistic axions is
$\rho_{a}(t)=\tilde{m}_{a}(t)\,n_{a}(t)$, and at present (i.e., at $t=t_{0}$)
becomes

$\rho_{a}(t_{0})=\rho_{a}(t_{\text{osc}})\,\frac{m_{a}}{\tilde{m}_{a}(t_{\text{osc}})}\left(\frac{R(t_{\text{osc}})}{R_{0}}\right)^{3}\simeq\frac{\theta_{i}^{2}}{2}\,f_{a}^{2}\,m_{a}\,\tilde{m}_{a}(T_{\text{osc}})\,\frac{s(T_{0})}{s(T_{\text{osc}})}\,,$
(14)

with $R_{0}\equiv R(t_{0})$ the scale factor at present, and using entropy
conservation and
$\rho_{a}(T_{\text{osc}})\simeq\frac{1}{2}\,\tilde{m}_{a}^{2}(T_{\text{osc}})\,f_{a}^{2}\,\theta_{i}^{2}$.
The oscillation temperature $T_{\text{osc}}$ is therefore

$\displaystyle
T_{\text{osc}}\simeq\begin{dcases}\left(\frac{1}{\pi}\sqrt{\frac{10}{g_{\star}(T_{\text{osc}})}}\,m_{a}\,M_{P}\right)^{1/2}&\text{for
}\,T_{\text{osc}}\leq T_{\rm QCD}\,,\\\
\left(\frac{1}{\pi}\sqrt{\frac{10}{g_{\star}(T_{\text{osc}})}}\,m_{a}\,M_{P}\,T_{\rm
QCD}^{4}\right)^{1/6}&\text{for }\,T_{\text{osc}}\geq T_{\rm
QCD}\,.\end{dcases}$ (15)

The dependence of $T_{\text{osc}}$ with $m_{a}$ is depicted in the left panel
of Fig. 1, where the change in the slope at $T_{\text{osc}}=T_{\rm QCD}$ is
shown with a horizontal line.

Figure 1: Left panel: Oscillation temperature as function of the axion mass.
The bend around $m_{a}\simeq 4.8\times 10^{-11}$ eV corresponds to the QCD
transition. Right panel: Initial misalignment angle as a function of the axion
mass. In both panels, the blue solid line is for the numerical solution and
dotted blue line is for the analytical solution. The mismatch for
$\theta_{i}>\pi/\sqrt{3}$ is due to the fact that we do not consider
anharmonic effect in the analytical solution. In both panels the dark grey
areas correspond to the scale $f_{a}$ bigger than the Planck scale.

The corresponding axion mass is

$\tilde{m}_{a}(T_{\rm QCD})=\pi\,\sqrt{\frac{g_{\star}(T_{\rm
QCD})}{10}}\,\frac{T_{\rm QCD}^{2}}{M_{P}}\simeq 4.8\times
10^{-11}~{}\text{eV}.$ (16)

Thus, for masses above $\tilde{m}_{a}(T_{\rm QCD})$, the temperature effects
of QCD are important, and below they can be safely ignored. Now we can insert
$T_{\text{osc}}$ back into Eq. (14) in order to have the axion relic density
today

$\Omega_{a}\equiv\frac{\rho_{a}(T_{0})}{\rho_{c}}\simeq\begin{dcases}0.006\left(\frac{\theta_{i}}{1}\right)^{2}\left(\frac{m_{a}}{5.6~{}\mu\mbox{eV}}\right)^{-3/2}&\text{
for }m_{a}\lesssim\tilde{m}_{a}(T_{\rm QCD})\,,\\\
0.17\left(\frac{\theta_{i}}{1}\right)^{2}\left(\frac{m_{a}}{5.6~{}\mu\mbox{eV}}\right)^{-7/6}&\text{
for }m_{a}\gtrsim\tilde{m}_{a}(T_{\rm QCD})\,,\end{dcases}$ (17)

where we have used the fact that $\rho_{c}\simeq 1.1\times 10^{-5}~{}h^{2}$
GeV/cm3 is the critical energy density, $h\simeq 0.674$, and $s(T_{0})\simeq
2.9\times 10^{3}$ cm-3 [14]. In order to match the observed DM relic abundance
$\Omega_{a}\simeq 0.26$ [14], the initial misalignment angle $\theta_{i}$
could be tuned. The right panel of Fig. 1 shows the required angle
$\theta_{i}$ in order to generate the whole observed DM relic abundance, for
different axion masses. It is worth pointing-out that the analytical
approximation deviates significantly from the numerical result for
$\theta_{i}\gtrsim 1$, since the former uses the exact potential as well as
the corresponding anharmonic factor derived in Appendix C.

The initial angle is deeply connected to the cosmological history of the
axion: the realignment of the axion potential happens in different causally
disconnected patches in the universe with (in principle) different initial
misalignment angles. If the PQS is broken before or during inflation (and not
restored afterwards) one of these patches gets inflated away, homogenising the
initial angle over the observable universe and the angle is randomly selected
from $\left[-\pi,\pi\right]$. On the other hand, if the PQS breaking occurs
after the inflationary epoch, different initial angles in the different
disconnected patches are arbitrarily selected between $\left[-\pi,\pi\right]$,
thus, the angle is averaged to
$\langle\theta_{i}^{2}\rangle=\pi^{2}/3\simeq(1.8)^{2}$ [15]. Therefore, the
parameter space where the axion can explain the DM paradigm can be different,
depending on whether $f_{a}$ is higher or lower than $T_{I}\equiv H_{I}/2\pi$,
where $H_{I}$ corresponds to the inflationary scale. For the pre-inflationary
scenario, where $f_{a}\geq T_{I}$, there is a strong constraint on high
$f_{a}$ scales from isocurvature perturbations [15]. If the axion is present
during inflation, it develops quantum fluctuations which are nearly scale
invariant and uncorrelated with the adiabatic perturbations. Using Planck
constraints on the primordial isocurvature fraction [14], it is found that
[69]

$H_{I}\lesssim 1.3\times
10^{9}~{}\mbox{GeV}\left(\frac{f_{a}}{10^{16}\,\mbox{GeV}}\right)^{0.42},$
(18)

assuming axions make the whole DM abundance. On the other hand, an axion with
a PQ scale $f_{a}\simeq 10^{-16}$ GeV can produce the whole abundance for a
misalignment angle of $\theta_{i}\simeq 3\times 10^{-3}$ (cf. Fig. 1), which
is considered tuned or anthropic. That is, a pre-inflationary axion hints to a
low-scale inflation in the standard cosmological scenario.

In the post-inflationary scenario, $f_{a}<T_{I}$, the bound on the PQ scale
comes from axion DM not overclosing the universe, which requires
$f_{a}\lesssim 10^{11}~{}$GeV [18]. On the other hand, in this scenario
topological defects such as cosmic strings and domain walls emerge, radiating
cold axions. Their contribution it is still in dispute [17, 18, 19, 20, 9,
70], nonetheless it is usually assumed they enlarge the DM parameter space
down to $f_{a}\gtrsim 10^{9}$ GeV, translated into a mass range of
approximately $56~{}\mu\mbox{eV}\lesssim m_{a}\lesssim 5$ meV.

For the rest of our discussion – since we are only focused on studying the
impact of NSC on the misalignment mechanism – we will consider as the standard
cosmology axion mass range, the one obtained from the misalignment mechanism,
Eq. (17), by considering $\theta_{i}$ in a natural range of
$\left[0.5,\pi/\sqrt{3}\right]$, that leads to masses $1.6\times 10^{-6}$ eV
$\lesssim m_{a}\lesssim 1.4\times 10^{-5}$ eV.

##  3 Nonstandard cosmologies

It is typically assumed that the energy density of the universe was dominated
by radiation between the end of the reheating era after inflation and the
onset of matter domination at redshift $z\simeq 3400$, i.e., the standard
cosmological scenario. However, there is no evidence that the universe was
radiation dominated before the BBN. In that context, NSCs have been widely
studied, especially in the context of a fluid with an equation-of-state
parameter $\omega$, with $-1\leq\omega\leq 1$, that eventually dominates the
energy density of the universe.222Has been argued that cosmologies with
$\omega>1$ could feature superluminal propagation, although in Ref. [71] they
show that is not the case. Moreover, as discussed in Refs. [72, 73, 74],
scenarios with $\omega\to-1$ are disfavoured. A more concrete bound on
$\omega$ (and the other parameters related to $\phi$) would depend on the
underlying inflationary model, which is beyond our scope. Extensively studied
NSCs are the EMD, firstly explored in the context of supersymmetry and string
theory (see e.g. Refs. [75, 76, 77]), with $\omega=0$, and the kinetic energy
domination, known as kination dominance, with $\omega=1$ [78, 79, 80]. For a
detailed review on NSC, see Ref. [81]. In this work, we will consider the
existence of a particle field $\phi$ with a general equation of state,
$\omega=p_{\phi}/\rho_{\phi}$, that during a certain period has an impact on
the expansion of the universe. Prior to BBN, the field could decay into SM
radiation, with a total decay rate $\Gamma_{\phi}$. The Boltzmann equations
that govern the evolution of the energy density $\rho_{\phi}$ of $\phi$ and
the SM entropy density $s$ are333It is interesting to note that, in the case
of $\omega>1/3$, $\rho_{\phi}$ gets diluted faster than radiation, and
therefore $\phi$ could be stable, with $\Gamma_{\phi}=0$.

$\displaystyle\frac{d\rho_{\phi}}{dt}+3(1+\omega)\,H\,\rho_{\phi}$
$\displaystyle=-\Gamma_{\phi}\,\rho_{\phi}\,,$ (19)
$\displaystyle\frac{ds}{dt}+3\,H\,s$
$\displaystyle=+\frac{\Gamma_{\phi}}{T}\,\rho_{\phi}\,.$ (20)

The latter equation can be recast to find the relationship between the plasma
temperature $T$ and scale factor $R$ in the case of an entropy injection due
to the decay of $\phi$, as

$\frac{dT}{dR}=\left(1+\frac{T}{3\,g_{\star S}}\frac{dg_{\star
S}}{dT}\right)^{-1}\left[-\frac{T}{R}+\frac{\Gamma_{\phi}\,\rho_{\phi}}{3\,H\,s\,R}\right].$
(21)

The Hubble expansion rate is given by

$H=\sqrt{\frac{\rho_{\phi}+\rho_{R}+\rho_{a}}{3M_{P}^{2}}}\,,$ (22)

featuring the contribution from $\phi$. As the axion contribution is always
subdominant, the evolution of the $\phi$-radiation system is decoupled from
the DM evolution and can be solved separately. It is customary to define the
end of the $\phi$ domination with the temperature $T_{\text{end}}$. For
$\omega>1/3$, as $\phi$ does not have to decay, $T_{\text{end}}$ corresponds
to the equality $\rho_{R}(T_{\text{end}})=\rho_{\phi}(T_{\text{end}})$.
However, for $\omega<1/3$, it corresponds to the temperature at which the
field has mostly decayed away, i.e. when $H(T_{\text{end}})=\Gamma_{\phi}$,
and therefore [82, 59]

$T_{\text{end}}^{4}\equiv\frac{90}{\pi^{2}\,g_{\star}(T_{\text{end}})}\,M_{P}^{2}\,\Gamma_{\phi}^{2}\,.$
(23)

Equations (19) and (20) can be analytically solved. Here we present simplified
approximations that will be used in our analytical analysis later. However,
the complete details are outlined in the Appendix A. Before the decay of
$\phi$, the energy densities evolve as

$\rho_{R}(R)=\rho_{\text{eq}}\left(\frac{R_{\text{eq}}}{R}\right)^{4},\,\,\,\,\,\,\rho_{\phi}(R)=\rho_{\text{eq}}\left(\frac{R_{\text{eq}}}{R}\right)^{\beta},$
(24)

where we have defined $\beta\equiv 3(\omega+1)$ and $\rho_{\text{eq}}$ is the
energy density of $\phi$ and radiation at $T=T_{\text{eq}}$, with a scale
factor $R=R_{\text{eq}}$. Eventually, the decay of $\phi$ begins to affect the
evolution of the plasma temperature, at some scale factor $R_{c}$. Thus, we
solve analytically Eqs. (19) and (20) at first order in
$\Gamma_{\phi}/H_{\text{eq}}$, where $H_{\text{eq}}\equiv H(R_{\text{eq}})$,
following closely Ref. [74]. We find that the energy of $\phi$ and radiation
is well described by

$\displaystyle\rho_{\phi}(R)$
$\displaystyle\simeq\rho_{\text{eq}}\,\left[\left(\frac{R_{\text{eq}}}{R}\right)^{\beta}-\frac{2}{\beta}\frac{\Gamma_{\phi}}{H_{\text{eq}}}\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/2}\right],$
(25) $\displaystyle\rho_{R}(R)$
$\displaystyle\simeq\rho_{\text{eq}}\,\left[\left(\frac{R_{\text{eq}}}{R}\right)^{4}+\frac{2}{8-\beta}\frac{\Gamma_{\phi}}{H_{\text{eq}}}\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/2}\right].$
(26)

For $R\ll R_{c}$ the first terms of the rhs of Eqs. (25) and (26) are the
dominant ones. Eventually, as the decay of $\phi$ starts to be important, the
second terms of the rhs can not be neglected. We find, by equating the first
and second term of the rhs of Eq. (26) and using the expressions of Eq. (78),
that

$R_{\text{c}}\simeq
R_{\text{eq}}\left(\frac{(8-\beta)}{2}\left(\frac{T_{\text{eq}}}{T_{\text{end}}}\right)^{2}\right)^{\frac{2}{8-\beta}}.$
(27)

And by requiring the second and third term of the rhs of Eq. (25) are
comparable, we find an expression for when the decays start to be important,
meaning $R=R_{\text{end}}$, to be

$R_{\text{end}}\simeq
R_{\text{eq}}\left(\frac{\beta}{2}\left(\frac{T_{\text{eq}}}{T_{\text{end}}}\right)^{2}\right)^{2/\beta}.$
(28)

where we can also find the temperature at $R=R_{\text{c}}$ to be

$T_{\text{c}}\simeq
T_{\text{eq}}\left(\frac{2}{8-\beta}\frac{T_{\text{end}}^{2}}{T_{\text{eq}}^{2}}\right)^{\frac{2}{8-\beta}}.$
(29)

From this analysis it can be also extracted that deep during the $\phi$
domination, the relation between temperature and scale factor is

$T(R)\simeq
T_{\text{eq}}\left[\frac{2}{8-\beta}\frac{T_{\text{end}}^{2}}{T_{\text{eq}}^{2}}\right]^{1/4}\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/8}.$
(30)

Figure 2 shows an example of the evolution of the SM and $\phi$ energy
densities, for $T_{\text{eq}}=1$ GeV, $T_{\text{end}}=4$ MeV and $\beta=3$
(EMD).

Figure 2: Evolution of the energy densities for radiation (black) and $\phi$
field (dashed red) as a function of the scale factor $R$, for $\beta=3$,
$T_{\text{eq}}=1$ GeV, and $T_{\text{end}}=4$ MeV. The vertical lines depict
$R=R_{\text{eq}}$ (solid blue), $R=R_{c}$ (dashed blue), and
$R=R_{\text{end}}$ (dash-dotted line), respectively.

For the case of an EMD, or in general for $\beta<4$, $R_{\text{eq}}$
corresponds to the scale factor at which $\rho_{\phi}$ starts to dominate over
$\rho_{R}$, $R_{c}$ to the scale factor where effectively $\rho_{\phi}$ starts
to dominate the evolution of $\rho_{R}$, and $R_{\text{end}}$ is a proxy of
the scale factor where $\phi$ decays completely.

Figure 3: The evolution of the Hubble expansion rate $H$ normalised to the
case dominated by radiation $H_{R}$, as a function of the temperature for
different NSC. The left panel corresponds to $T_{\text{eq}}=10^{8}$ GeV and
$T_{\text{end}}=4$ MeV, for $\beta=1$, 2 and 3, whereas the right panel to
$T_{\text{end}}=4$ MeV for $\beta=6$ and 8.

Taking into account the contributions from radiation and $\phi$, Fig. 3 shows
the evolution of the Hubble expansion rate $H$ normalised to case dominated by
radiation $H_{R}$, as a function of the temperature for different NSCs. The
left panel corresponds to $T_{\text{eq}}=10^{8}$ GeV and $T_{\text{end}}=4$
MeV, for $\beta=1$, 2 and 3, whereas the right panel to $T_{\text{end}}=4$ MeV
and $\beta=6$ and 8. The different regions can be understood analytically. In
the case with $\beta<4$ (left panel)

$H(T)\simeq\begin{dcases}H_{R}(T)&\text{ for }T\gg T_{\text{eq}}\,,\\\
H_{R}(T_{\text{eq}})\left(\frac{T}{T_{\text{eq}}}\right)^{\frac{\beta}{2}}&\text{
for }T_{\text{eq}}\gg T\gg T_{\text{c}}\,,\\\
H_{R}(T_{\text{end}})\left(\frac{T}{T_{\text{end}}}\right)^{4}&\text{ for
}T_{\text{c}}\gg T\gg T_{\text{end}}\,,\\\ H_{R}(T)&\text{ for
}T_{\text{end}}\gg T\,,\end{dcases}$ (31)

whereas in the opposite case with $\beta>4$ (right panel)

$H(T)\simeq\begin{dcases}H_{R}(T_{\text{end}})\left(\frac{T}{T_{\text{end}}}\right)^{\frac{\beta}{2}}&\text{
for }T\gg T_{\text{end}}\,,\\\ H_{R}(T)&\text{ for }T_{\text{end}}\gg
T\,.\end{dcases}$ (32)

##  4 Axion oscillation in NSC

In order to compute the relic density at present, we follow the same lines as
in the standard case, taking into account the entropy injection due to the
decay of $\phi$. Thus, the axion energy density today gets the form

$\rho_{a}(T_{0})=\rho_{a}(T_{\text{osc}})\frac{m_{a}}{\tilde{m}_{a}(T_{\text{osc}})}\left(\frac{R_{\text{osc}}}{R_{0}}\right)^{3}=\rho_{a}(T_{\text{osc}})\frac{m_{a}}{\tilde{m}_{a}(T_{\text{osc}})}\frac{s(T_{0})}{s(T_{\text{osc}})}\times\frac{S_{\rm
osc}}{S_{\rm end}}.$ (33)

which is similar to Eq. (14), multiplied by a factor that takes into account
the dilution of the axion energy density due to entropy injection between
$T_{\text{osc}}$ and $T_{\text{end}}$, given by $S_{\rm osc}/S_{\rm
end}$444Let us precise that in this equation, an instantaneous decay of $\phi$
has been assumed at $T=T_{\text{end}}$, which of course is not true.
Numerically we are taking this into account, defining a temperature/scale
factor such that the field has effectively stopped to affect the SM radiation.
For sake of clarity we will use $T_{\text{end}}$ in our analytical
expressions.. The dilution factor can be expressed in terms of the parameters
of the NSC, depending on which stage of the $\phi$ dominance the axion starts
to oscillate. Moreover, this factor appears only for $\beta<4$, as for
$\beta>4$ $\phi$ does not deposit any energy in the plasma and then is equal
to one.

Let us stress that the axion energy density – in absence of dilution – is
inversely proportional to $T_{\text{osc}}$. In a NSC scenario, the oscillation
temperature will always be smaller or equal to the oscillation temperature in
a radiation dominated expansion, because $H>H_{R}$ (see Figs. 3 and 4). Then,
the effect of lowering the oscillation temperature in a NSC leads to an
increased axion energy density (as will be found when analysing cosmologies
with no entropy injection, i.e. $\beta>4$). On the other hand, the dilution
factor in Eq. (33) causes the axion relic abundance to decrease by diluting
its energy density.

In the rest of this section we will obtain analytical expressions for the
axion relic density in terms of the parameters of the NSC and we will find the
extreme range of masses that could be achieved for a given cosmology. For
simplicity of the analytical expressions, we drop the dependence on the
degrees of freedom (i.e., they will be taken to one), as we are only
interested in the order of magnitude. For the numerical results they are fully
considered, the same as the anharmonicities in the potential, as outlined in
Appendix C.

###  4.1 Early matter domination and $\beta<4$

Cosmologies with $\beta<4$ correspond to equations of state with $\omega<1/3$.
It turns useful for our detailed analysis to divide the NSC into three stages,
which will lead to different outcomes for the axion relic density, depending
in which one the oscillation of the axion field takes place:

  * $i)$

Region 1, $R\ll R_{\text{eq}}$: radiation dominated universe.

  * $ii)$

Region 2, $R_{\text{eq}}\ll R\ll R_{\text{c}}$: dominance of $\phi$, prior its
decay becomes significant.

  * $iii)$

Region 3, $R_{\text{c}}\ll R\ll R_{\text{end}}$: dominance of $\phi$, which
significantly decays into radiation, modifying the evolution of the SM
temperature. After this stage, the universe returns to a radiation dominance.

Figure 4: The oscillation temperature as a function of the axion mass for a
non-standard (solid black) and standard (dotted black) cosmology. The
horizontal blue lines denote the boundaries of the regions studied in section
4.1. The left panel corresponds to $\beta=3$, $T_{\text{eq}}\simeq 10^{3}$
GeV, $T_{\text{end}}=4$ MeV, whereas the right panel to kination $\beta=6$
with $T_{\text{end}}=4$ MeV.

In Fig. 4 we show the oscillation temperature $T_{\text{osc}}$ vs the axion
mass $m_{a}$ in a NSC scenario with $\beta=1$ and we have highlighted the
boundaries of the three NSC regions, and can be clearly recognised that for
each of them $T_{\text{osc}}$ features a different dependence on $m_{a}$ that
we will find analytically.

Our aim is to find the smallest axion mass – equivalently, the highest
possible $f_{a}$ scale – that can accommodates the whole cold DM observed
today for a given NSC. To do that, we start writing the entropy dilution
factor for the axion in each region. Secondly, we find the requirements on the
axion mass in order for the oscillation of the axion field to take place
during any of the stages of the NSC delineated above. The oscillation
temperature can be divided into two regimes, depending on whether it takes
place before or after the QCD phase transition. Therefore, we divide our
analysis of each region for the case of constant axion mass (oscillations
below the QCD transition) or thermal axion mass (oscillations above the QCD
transition). For both regimes, we find the smallest possible mass that gives
the whole relic density, as a function of the NSC parameters, by requiring the
highest possible dilution taking some NSC benchmark values. Our main results
for $\beta<4$ are summarised in Fig. 5, where we show the parameter space that
generates the whole DM abundance for $\beta=3$ (red band, top) and $\beta=2$
(blue band, bottom), for $T_{\text{end}}=4$ MeV. We have considered as a width
a range of initial angles $\theta_{i}\in\left[0.5,\pi/\sqrt{3}\right]$. The
regions to the left of the bands produce a DM overabundance and the ones to
the right a sub-production. We compare our analytical estimates with these
numerical results, finding an excellent agreement. The region between the
vertical dashed lines depicts the case of a standard cosmology.

Figure 5: Parameter space corresponding to the whole observed DM abundance,
for an early matter dominance case (EMD) with $\beta=3$ (top) and $\beta=2$
(bottom), for $T_{\text{end}}=4$ MeV and $0.5\leq\theta_{i}\leq\pi/\sqrt{3}$.
Between the vertical dashed lines the results for the standard cosmology are
recovered. The dotted lines correspond to $T_{\text{osc}}=T_{\text{eq}}$,
$T_{\text{osc}}=T_{c}$ and $T_{\text{osc}}=T_{\rm QCD}$. The grey bands show
the regions where $T_{\text{eq}}<T_{\text{BBN}}$ or $f_{a}>M_{P}$.

####  4.1.1 Region 1: $\boldsymbol{R_{\text{osc}}\ll R_{\text{eq}}}$

In this region, the oscillations of the axion field occur in a period of
radiation domination, before the influence of $\phi$ and only differs from the
standard cosmological scenario in that there is entropy transfer to the SM,
not to axions, causing a dilution of the DM density. We can get an analytical
estimation of the dilution factor, valid when the fluid gets to dominate the
expansion of the universe,555Meaning valid only when the chosen parameters
$\beta$, $T_{\text{eq}}$ and $T_{\text{end}}$ lead to a solid fluid
domination. as

$\gamma_{R_{1}}=\frac{S_{\rm osc}}{S_{\rm
end}}\simeq\left[\left(\frac{4}{\beta^{2}}\right)\left(\frac{T_{\text{end}}}{T_{\text{eq}}}\right)^{4-\beta}\right]^{3/\beta},$
(34)

for a more detailed derivation, see Appendix B.1, Eq. (B.1). Therefore, in
this regime,

$\Omega_{a}=\Omega_{a}^{\rm std}\,\gamma_{R_{1}}\,,$ (35)

with $\Omega_{a}^{\rm std}$ corresponding to the axion density in the standard
cosmological scenario. By inspecting $\gamma_{R_{1}}$ it can be seen that the
dilution is more important for cosmologies with low $\beta$. Moreover, a
higher dilution is obtained for the small $T_{\text{end}}$ and large
$T_{\text{eq}}$ values.

The condition for the oscillation of the axion field to occur during this
period is given by $T_{\text{eq}}\ll T_{\text{osc}}$. Let us analyse it for
both regimes.

#### Constant axion mass

The constant axion mass regime corresponds to $T_{\text{eq}}\ll
T_{\text{osc}}\lesssim T_{\rm QCD}$, which can also be written in terms of the
axion mass, by using the expression for the oscillation temperature in Eq.
(15), as

$4\times
10^{-12}~{}{\rm{eV}}\left(\frac{T_{\text{eq}}}{0.1~{}\mbox{GeV}}\right)^{2}\ll
m_{a}\lesssim 10^{-11}~{}\mbox{eV},$ (36)

independently on $\beta$. The axion relic abundance in this region is obtained
from Eq. (35). Imposing that the axion constitutes all DM of the universe, the
axion mass in terms of the NSC parameters obeys

$m_{a}\simeq\begin{dcases}5\times
10^{-10}\,\mbox{eV}\,\,\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{2}\left(\frac{0.1\,\mbox{GeV}}{T_{\text{eq}}}\right)^{2}\,\theta_{i}^{4/3}&\text{for
}\beta=2\,,\\\ 2\times
10^{-8}\,\mbox{eV}\,\,\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{2/3}\left(\frac{0.1\,\mbox{GeV}}{T_{\text{eq}}}\right)^{2/3}\,\theta_{i}^{4/3}&\text{for
}\beta=3\,,\end{dcases}$ (37)

for $m_{a}\lesssim 10^{-11}$ eV. In order to maximise the impact of the
dilution, we have chosen $T_{\text{end}}\sim T_{\text{BBN}}$.

For the two equations of state considered above, we have to satisfy that the
axion mass is below $\sim 10^{-11}$eV, but that would need of a
$T_{\text{eq}}$ which is (much) higher than $T_{\rm QCD}$ or the initial
misalignment angle should be taken unnaturally small. Therefore, we find that
if the oscillation of the axion field happens during this regime, it can not
accommodate the whole DM density and could only be a subcomponent.

#### Thermal axion mass

On the other hand, for $T_{\text{osc}}\gtrsim T_{\rm QCD}$, the axion mass
receives important thermal corrections. Then, since $T_{\text{osc}}\gg
T_{\text{eq}}$, we find

$\displaystyle
10^{-12}\,\mbox{eV}\left(\frac{T_{\text{eq}}}{0.1\,\mbox{GeV}}\right)^{6}\ll
m_{a},\,\,\,\,\,\,\,\,\,\,\,\,\mbox{for}\,\,\,m_{a}\gtrsim 10^{-11}\rm{eV}.$
(38)

Note that $T_{\text{eq}}$ appears to the sixth power, so raising its value
also increases the value of the axion mass that could oscillate in this
period. Then, it seems difficult to open up the DM axion window much if the
oscillation happens during this regime, since in order to get the highest
entropy dilution, a high $T_{\text{eq}}$ is needed, but also this rises the
axion mass range, that already in the SC subproduces the DM density. Next, we
find the mass that can accommodate the whole DM relic density

$m_{a}\simeq\begin{dcases}2\times
10^{-10}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{18/7}\left(\frac{0.2\,\mbox{GeV}}{T_{\text{eq}}}\right)^{18/7}\theta_{i}^{12/7}&\text{for
}\beta=2\,,\\\ 3\times
10^{-6}\,{\rm{eV}}\,\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{6/7}\left(\frac{0.1\,\mbox{GeV}}{T_{\text{eq}}}\right)^{6/7}\,\theta_{i}^{12/7}&\text{for
}\beta=3\,,\end{dcases}$ (39)

for $m_{a}\gtrsim 10^{-11}$ eV. We have considered $T_{\text{end}}\sim
T_{\text{BBN}}$, to get the maximum dilution possible, thus we also require
the highest $T_{\text{eq}}$, but recalling from Eq. (38), that increasing that
temperature, makes the axion mass to rise, a compromise has to be found. Thus,
by incorporating that requirement, we find for $\beta=3$ that around
$T_{\text{eq}}\lesssim 0.6$ GeV for $\theta_{i}\sim 1$, an axion of mass
$m_{a}\sim 6\times 10^{-7}$ eV could explain the whole DM density. For
$\beta=2$, it is found that the smallest axion mass that can satisfy the DM
abundance is around $m_{a}\sim 10^{-10}$ eV, with a $T_{\text{eq}}\sim 0.2$
GeV. For those estimations we have assumed $\theta_{i}\simeq\mathcal{O}(1)$.
These findings are supported by Fig. 5, where a plot of the temperature
$T_{\text{eq}}$ vs $m_{a}$ is shown, such that the red band corresponds to the
mass range that fits the whole DM relic density observed today, with more
details in the caption. Region 2 for the case of thermal axion masses
contributes to the lower part of the plot. From that figure we confirm that in
order for an axion to oscillate during the period defined by Region 2, the NSC
has to be quite short, with $T_{\text{BBN}}\lesssim T_{\text{eq}}\lesssim
0.1~{}$GeV and $T_{\text{end}}\sim T_{\text{BBN}}$. Moreover, this region
opens up to the range of mass just below the standard cosmological scenario,
even overlapping with that cosmology, allowing to reach PQ scales as high as
nearly $f_{a}\lesssim 10^{14}$ GeV for $\beta=3$ and $f_{a}\sim 10^{16}~{}$GeV
for $\beta=2$, with initial misalignment angles $\theta_{i}\sim 1$.

####  4.1.2 Region 2: $\boldsymbol{R_{\text{eq}}\ll R_{\text{osc}}\ll
R_{\text{c}}}$

In this regime, the axion oscillation happens during the domination of $\phi$,
but their decays still do not affect the evolution of the temperature in the
SM sector, so $T\propto R^{-1}$. The Hubble parameter – assuming the $\phi$
field completely dominates the expansion – in this case is given by

$H(R)\simeq\sqrt{\frac{\rho_{\phi}(R)}{3\,M_{P}^{2}}}\,,$ (40)

where the expression for $\rho_{\phi}$ to be considered is as in Eq. (24). All
analytical results from this region are collected in Appendix B.2. The
oscillation temperature is now dependent on the NSC parameters, so their
analytical expressions can be found in Eq. (94). Also the axion mass at the
QCD transition it is different than in the standard cosmological scenario. For
our analysis, it will be useful to write it down for some benchmark values,
see Eq. (95),

$m_{R_{2}}\simeq\begin{dcases}10^{-10}\,\mbox{eV}\,\left(\frac{T_{\text{eq}}}{2\,\mbox{GeV}}\right)&\beta=2,\\\
2.4\times
10^{-9}\,\mbox{eV}\left(\frac{T_{\text{eq}}}{10^{4}\,\mbox{GeV}}\right)^{1/2}&\beta=3.\end{dcases}$
(41)

#### Constant axion mass

In the case of a constant axion mass, the requirement for the oscillation to
occur in Region 2 can be written as $T_{\text{c}}\ll T_{\text{osc}}\lesssim
T_{\rm QCD}$, because $T_{\text{c}}<T_{\rm QCD}<T_{\text{eq}}$ in the whole
range. Using the expression for the oscillation temperature Eq. (94), we
convert the oscillation temperature range into a mass range for the axion,
given in Eq. (97). Let us write it down for some benchmark values

$\displaystyle\beta=2;\qquad 2.7\times
10^{-11}\,\mbox{eV}\,\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{2/3}\left(\frac{T_{\text{eq}}}{2~{}\mbox{GeV}}\right)^{4/3}\ll
m_{a}\ll m_{R_{2}},$ (42) $\displaystyle\beta=3;\qquad 9\times
10^{-10}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{6/5}\left(\frac{T_{\text{eq}}}{10^{4}~{}\mbox{GeV}}\right)^{4/5}\ll
m_{a}\ll m_{R_{2}}.$ (43)

Thus, for both $\beta=2$ and $\beta=3$, the mass range can lie below the one
for standard cosmology, and it gets lower for $\beta=2$, but, they have a
rather small range that satisfies $m_{a}\lesssim m_{R_{2}}$, by comparing with
Eq. (41). Note that in the case $\beta=2$, we have to choose a fairly small
$T_{\text{eq}}$ to meet the correct mass range.

The axion relic density is computed evaluating the energy density with the
corresponding oscillation temperature in this region, and then multiplying by
the dilution factor, which is the same as Region 1, since the fluid has not
yet started to decay. The detailed expressions for the axion relic density are
computed in the appendix, Eq. (98), from those expressions, we can find the
axion mass that accommodates the whole DM density in a NSC if it happens that
the axion oscillates in Region 2 with constant mass. As done before, we equate
the relic density in this region to $\Omega_{\rm CDM}$ and we again consider
$T_{\text{end}}\sim T_{\text{BBN}}$ to assure the maximum dilution and the
initial misalignment angle is considered $\mathcal{O}(1)$. We find

$m_{a}\simeq\begin{dcases}4.7\times
10^{-11}~{}\text{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right){\theta_{i}}^{2/3}&\beta=2\,,\\\
2.7\times
10^{-9}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{1/2}\,{\theta_{i}}&\beta=3\,,\end{dcases}$
(44)

for $m_{a}\lesssim m_{R_{2}}$. Let us first stress that the relic density is
independent on $T_{\text{eq}}$. Secondly, we require they are within the
limits of Eq. (43) which leads to the condition
$10^{4}\,\mbox{GeV}<T_{\text{eq}}<3.6\times 10^{4}$ GeV for $\beta=3$ and
$0.7\,\mbox{GeV}<T_{\text{eq}}<4$ GeV for $\beta=2$. Therefore, we have found
that the smallest axion mass that that could fit the DM relic abundance in the
case of $\beta=3$ is

$m_{a}\sim\mbox{few}\times
10^{-9}\,\mbox{eV},\,\,\,\,\,\,\,\mbox{with}\,\,\,\,\,\,10^{4}\,\mbox{GeV}<T_{\text{eq}}<3.6\times
10^{4}\,\mbox{GeV},$ (45)

which can be seen from Fig. 5 that it corresponds to the region in the red
band that it is roughly between the green dotted line ($T_{\text{osc}}=T_{\rm
QCD}$) and the blue line ($T_{\text{osc}}=T_{\text{c}}$). Analogously, for
$\beta=2$ the smallest axion mass is found roughly

$m_{a}\sim\mbox{few}\times 10^{-11}\mbox{eV},\quad\text{with}\quad
0.7\,\mbox{GeV}<T_{\text{eq}}<4\,\mbox{GeV},$ (46)

which also corresponds to the region in the lower panel of Fig. 5 that lies
between the dotted green line ($T_{\text{osc}}=T_{\rm QCD}$) and the dotted
blue line ($T_{\text{osc}}=T_{\text{c}}$). As we shown in our analysis, this
region has a short range, meaning that only for cosmologies with narrow range
of temperatures $T_{\text{eq}}$ can give the right dilution to the relic
density, such that the whole DM relic density can be produced.

#### Thermal axion mass

In case the oscillation happens above the QCD phase transition temperature,
the range to be considered is $T_{\rm QCD}\ll T_{\text{osc}}\ll
T_{\text{eq}}$, using expression Eq. (94) from the Appendix and converting the
above requirement into an axion range, we get

$\displaystyle m_{R_{2}}\ll m_{a}\ll 4\times
10^{-8}\mbox{eV}\left(\frac{T_{\text{eq}}}{0.6\,\mbox{GeV}}\right)^{6},$ (47)

valid for any equation of state. From the relic density found in Eq. (98) we
obtain the axion mass that can fit the whole DM relic density to be of the
order of

$m_{a}\simeq\begin{dcases}10^{-10}\,\mbox{eV}\left(\frac{0.1\,\rm{GeV}}{T_{\text{eq}}}\right)^{8/7}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{15/7}\theta_{i}^{10/7}&\text{for
}\beta=2\,,\\\ 2\times
10^{-9}\,\mbox{eV}\left(\frac{10^{4}~{}\rm{GeV}}{T_{\text{eq}}}\right)^{2/7}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{11/14}\theta_{i}^{11/7}&\text{for
}\beta=3\,,\end{dcases}$ (48)

for $m_{a}\gtrsim m_{R_{2}}$. Requiring they belong to the mass range in Eq.
(47), we get that the equilibrium temperature has to satisfy
$0.6\,\mbox{GeV}<T_{\text{eq}}<10^{4}$ GeV for $\beta=3$ and
$0.2\,\mbox{GeV}<T_{\text{eq}}<0.7\,\mbox{GeV}$ for $\beta=2$. Comparing with
our numerical result from Fig. 5 we find a good agreement. In the case of EMD
($\beta=3$) there is a much broader range of cosmologies that can lead to the
right amount of DM production. But, because the relic has a steeper dependence
on the mass than in the standard cosmological scenario ($m_{a}^{-14/11}$ vs.
$m_{a}^{-7/6}$), the range of masses covered is not as broad as in the
previous region. Instead, for $\beta=2$ our analysis supported by the lower
panel of Fig. 5 shows that the parameter space of NSC such that the axion DM
is produced during this regime it is extremely narrow (between black and green
dotted lines) and the axion mass is set to be between $10^{-10}-10^{-9}$ eV,
which actually opens up the PQ scale more than for a EMD era, to near
$f_{a}\sim 10^{17}$ GeV.

####  4.1.3 Region 3: $\boldsymbol{R_{\text{c}}\ll R_{\text{osc}}\ll
R_{\text{end}}}$

For this regime, oscillations of the axion field happen during the domination
and decay of $\phi$. The Hubble parameter is taken to be
$H\simeq\sqrt{\frac{\rho_{\phi}}{3M_{P}}}$ and temperature and scale factor
are related by $T\propto R^{-\beta/8}$, see Eq. (30). Therefore, this time
results more convenient to impose the condition $R_{\text{c}}\ll
R_{\text{osc}}\ll R_{\text{end}}$, than using the temperatures. Full
analytical expressions used for this region can be found in Appendix B.3. Let
us also remind that during this period, due to the entropy injection, the
expansion of the universe follows the behaviour $H\propto T^{4}$. The dilution
to the axion energy is different than in the previous regions, because $\phi$
is already decaying. Using $R_{\text{end}}$ given in Eq. (28) we find

$\gamma_{R_{3}}=\frac{S_{\rm osc}}{S_{\rm
end}}=\left(\frac{T_{\text{osc}}}{T_{\text{end}}}\frac{R_{\text{osc}}}{R_{\text{end}}}\right)^{3}\simeq\left(\frac{T_{\text{end}}}{T_{\text{osc}}}\right)^{24/\beta-3}.$
(49)

For smaller $\beta$, as expected, the dilution is more important. On the other
hand, since $T_{\text{end}}/T_{\text{osc}}$ is bigger than the ratio
$T_{\text{end}}/T_{\text{eq}}$ from Eq. (34), we obtain the relation
$\gamma_{3}>\gamma_{1}$, since in Region 3 the field $\phi$ is already
decaying, the entropy injection into the thermal bath is smaller than in the
previously considered regions. The highest dilution of the axion energy
density will be attained for the smallest $T_{\text{end}}$, so we shall again
choose $T_{\text{end}}=T_{\text{BBN}}$, for our analysis.

For the oscillation temperatures, the expressions have a very mild dependence
on $\beta$, that does not change the order of magnitude of the expressions
(see Eq. (103)), so we write them as valid for any equation of state

$\displaystyle
T_{\text{osc}}^{R_{3}}\simeq\begin{dcases}0.04~{}\mbox{GeV}\left(\frac{m_{a}}{10^{-4}\,\mu\mbox{eV}}\right)^{1/4}\left(\frac{T_{\text{end}}}{10~{}\mbox{MeV}}\right)^{1/2}&m_{a}\lesssim
m_{R_{3}}\,,\\\
1~{}\mbox{GeV}\left(\frac{m_{a}}{1\,\mu\mbox{eV}}\right)^{1/8}\left(\frac{T_{\text{end}}}{10~{}\mbox{GeV}}\right)^{1/4}&m_{a}\gtrsim
m_{R_{3}}\,,\end{dcases}$ (50)

where $m_{R_{3}}$ is the mass at the QCD phase transition, given by Eq. (104)

$m_{R_{3}}\simeq 3\times
10^{-8}~{}\rm{eV}\left(\frac{T_{\text{BBN}}}{T_{\text{end}}}\right)^{2}.$ (51)

Let us make a few comments about these expressions. Firstly, there is no
dependence of $T_{\text{eq}}$. Secondly, since we have considered
$T_{\text{end}}\sim T_{\text{BBN}}$, the mass in Eq. (51) is the highest at
which the transition could happen. Recalling that in the SC the transition
mass is $\sim 10^{-11}$ eV, the QCD phase transition is shifted to a slightly
higher mass in this scenario.

#### Constant axion mass

The mass range for the oscillation to happen during Region 3 can be obtained
from $R_{\text{c}}\ll R_{\text{osc}}\ll R_{\text{end}}$.666In this case we do
not include the scale factor at the QCD phase transition temperature, $R_{\rm
QCD}$, because it is found is always smaller than $R_{\text{osc}}$ for both,
constant and thermal mass regimes. We relate the scale factor to the
temperature using Eq. (30) and using the oscillation temperature found above,
we get

$\displaystyle\beta$ $\displaystyle=$ $\displaystyle 2;\quad 7\times
10^{-15}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{2}\ll
m_{a}\ll 5\times
10^{-11}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{2/3}\left(\frac{T_{\text{eq}}}{4\,\mbox{GeV}}\right)^{4/3},$
(52) $\displaystyle\beta$ $\displaystyle=$ $\displaystyle 3;\quad 4\times
10^{-15}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{2}\ll
m_{a}\ll 3\times
10^{-9}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{6/5}\left(\frac{T_{\text{eq}}}{10^{5}\,\mbox{GeV}}\right)^{4/5},$
(53)

for $m_{a}\lesssim m_{R_{3}}$. The mass range is in principle quite broad and
it moves to higher masses by increasing $T_{\text{end}}$ and only the upper
bound depends on $T_{\text{eq}}$. All the masses from the range are restricted
to be smaller than the mass at the QCD temperature, Eq. (51), which limits
$T_{\text{eq}}<10^{6}$ GeV for $\beta=3$ and $T_{\text{eq}}<12~{}$GeV for
$\beta=2$.

The analytical expression that accounts for the relic density during this
period is given in Eq. (110). From there, we find the axion masses that
accounts for the whole DM density for our two benchmark cosmologies, $\beta=2$
and $\beta=3$ are

$\displaystyle m_{a}\simeq\begin{dcases}5\times
10^{-11}\,\mbox{eV}\,\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)\,\theta_{i}^{2/3},&\qquad\beta=2,\\\
2.7\times
10^{-9}\,\mbox{eV}\,\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{1/2}\,\theta_{i},&\qquad\beta=3.\end{dcases}$
(54)

And requiring they are consistent with the mass ranges found in Eqs. (52) and
(52), respectively, it is found that $T_{\text{eq}}>4$ GeV for $\beta=2$ and
$T_{\text{eq}}>8\times 10^{4}~{}$GeV for $\beta=3$. These estimates are
supported by Fig. 5, where region 3 for constant axion mass starts above the
blue line $T_{\text{osc}}=T_{\text{c}}$. There, as found in our analysis, the
relic density is independent of $T_{\text{eq}}$ and the smallest axion mass
that can fulfil the whole DM density is, for $\beta=3$

$m_{a}\sim 10^{-8}\,\mbox{eV}\qquad\mbox{for}\quad T_{\text{eq}}\gtrsim
10^{4}\,\mbox{GeV},$ (55)

and for $\beta=2$

$m_{a}\sim 10^{-10}\,\mbox{eV}\qquad\mbox{for}\quad T_{\text{eq}}\gtrsim
4\,\mbox{GeV}.$ (56)

Numerically in Fig. 5, we have found that Region 3 starts to happen for
constant mass axions above the blue line $T_{\text{osc}}=T_{\text{c}}$, in
good agreement with our analysis and the numbers estimated above. Thus, for
approximately $T_{\text{eq}}\gtrsim 10^{3}$ GeV, the relic density is
independent on $T_{\text{eq}}$, where the first part corresponds to the
oscillation of the axion prior to the decay of $\phi$ (Region 2) and for
higher $T_{\text{eq}}$ when the decays are already happening, in a mass range
of few $\times 10^{-8}$ eV in mass and below, for a NSC with
$T_{\text{end}}\sim T_{\text{BBN}}$ and $\theta_{i}\sim 1$.

#### Thermal mass

If the oscillation of the axion field happens before the QCD phase transition,
the range to be considered is $R_{\text{c}}\ll R_{\text{osc}}\ll
R_{\text{end}}$, and it is mapped into the mass range by

$\displaystyle\beta$ $\displaystyle=$ $\displaystyle 2;\qquad m_{R_{3}}\ll
m_{a}\ll
10^{-7}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{10/3}\left(\frac{T_{\text{eq}}}{10^{3}\,\mbox{GeV}}\right)^{8/3},$
(57) $\displaystyle\beta$ $\displaystyle=$ $\displaystyle 3;\qquad
m_{R_{3}}\ll m_{a}\ll 5\times
10^{-7}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{22/5}\left(\frac{T_{\text{eq}}}{10^{7}\,\mbox{GeV}}\right)^{8/5}.$
(58)

Now using the equation for the relic density in this region, given by Eq.
(110) we obtain the axion mass that accounts for the whole DM relic density to
be approximately

$\displaystyle m_{a}\simeq\begin{dcases}2\times
10^{-11}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{5/2}\theta_{i}&\qquad\beta=2\,,\\\
10^{-8}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{4/3}\theta_{i}^{4/3}&\qquad\beta=3\,.\end{dcases}$
(59)

For $m_{a}\gtrsim m_{R_{3}}$.These expressions again do not depend on
$T_{\text{eq}}$ so if $T_{\text{end}}\sim T_{\text{BBN}}$, they correspond to
the lowest possible mass that can account for the whole DM relic density.
Comparing with the intersection mass given in Eq. (51), they are slightly
higher, meaning the relic density can not be completely fulfilled when axions
oscillate before the QCD phase transition for the range of masses found in
Eqs. (57) and (58). We conclude that axions in those mass ranges can only
account for the whole DM density if higher values of $T_{\text{end}}$ are
chosen, $T_{\text{end}}>1.4~{}T_{\text{BBN}}$ for $\beta=3$ and
$T_{\text{end}}>5~{}T_{\text{BBN}}$ for $\beta=2$.

###  4.2 Kination and $\beta>4$

In this case, the equation-of-state parameter $\omega>1/3$ (and hence
$\beta>4$), implying that the energy density of $\phi$ dilutes faster than the
one of radiation. Therefore, even if $\phi$ does not decay, it will naturally
become subdominant, giving rise to the onset of a radiation-dominated era, at
a temperature $T=T_{\text{end}}$ defined by the equality
$\rho_{\phi}(T_{\text{end}})=\rho_{R}(T_{\text{end}})$. It is interesting to
note that as $\phi$ does not decay, the SM entropy is always conserved. A
typical example of $\omega>1/3$ corresponds to kination [83, 84], where
$\omega=1$ (i.e. $\beta=6$), which implies that $H(T)\propto T^{3}$, or to the
case where $\omega=5/3$ (i.e. $\beta=8$) [85, 71] for which $H(T)\propto
T^{4}$.

Additionally, the oscillation temperature takes the analytical expression

$\frac{T_{\text{osc}}}{T_{\text{end}}}\simeq\begin{dcases}\left[\frac{g_{\star
S}(T_{\text{end}})}{g_{\star
S}(T_{\text{osc}})}\right]^{1/3}\left[\frac{m_{a}}{3\,H_{R}(T_{\text{end}})}\right]^{2/\beta}&\text{
for }T_{\text{osc}}\leq T_{\rm QCD}\,,\\\ \left[\frac{g_{\star
S}(T_{\text{end}})}{g_{\star
S}(T_{\text{osc}})}\right]^{\frac{\beta}{24+3\beta}}\left[\frac{m_{a}}{3\,H_{R}(T_{\text{end}})}\left(\frac{T_{\text{end}}}{T_{\rm
QCD}}\right)^{-4}\right]^{\frac{2}{8+\beta}}&\text{ for }T_{\text{osc}}\geq
T_{\rm QCD}\,,\end{dcases}$ (60)

for $T_{\text{osc}}\gg T_{\text{end}}$. In the opposite case where
$T_{\text{osc}}\ll T_{\text{end}}$, the oscillation takes place during the
standard cosmology, and therefore the NSC has no effect on the axion dynamics.

Note that the oscillation temperature for this scenario is in correspondence
with the one found for $\beta<4$ in region 2, with the replacement
$T_{\text{eq}}\rightarrow T_{\text{end}}$. This comes as no surprise, since
cosmologies with $\beta>4$ are characterised by the dominance of $\phi$,
without decaying, up to the temperature $T_{\text{end}}$. The same thing
happens for region 2, but the equality temperature $T_{\text{eq}}$ has to be
identified with $T_{\text{end}}$. Then, the mass at the QCD transition is
given by the same expression as the one found for region 2 in Eq. (95) with
the replacement $T_{\text{eq}}\rightarrow T_{\text{end}}$, which is

$m_{\rm{QCD}}\simeq\begin{dcases}4\times
10^{-10}\,\mbox{eV}\left(\frac{T_{\text{BBN}}}{T_{\text{end}}}\right)&\mbox{for
}\beta=6\,,\\\
10^{-8}\,\mbox{eV}\,\left(\frac{T_{\text{BBN}}}{T_{\text{end}}}\right)^{2}&\mbox{for
}\beta=8\,.\\\ \end{dcases}$ (61)

To find the axion mass range for which the oscillation is possible during this
period, we only have to impose $T_{\text{end}}\ll T_{\text{osc}}$, using the
appropriate expression for $T_{\text{osc}}$, depending if the scale is higher
or lower than the QCD temperature, $T_{\rm QCD}$.

#### Constant mass

In the case of constant axion mass, i.e. $T_{\text{osc}}\lesssim T_{\rm QCD}$,
we find the allowed mass range for oscillation in this regime to be

$\displaystyle\qquad 7\times
10^{-15}\,\mbox{eV}\left(\frac{T_{\text{end}}}{T_{\text{BBN}}}\right)^{2}\ll
m_{a}\lesssim m_{\rm QCD},$ (62)

valid for any $\beta$. Secondly, we again can borrow the expression found in
the analysis of region 2 for the relic density Eq. (98), with the dilution
factor of regions 1 and 2, Eq. (34) equal to one. We find the expression

$\displaystyle\Omega_{a}=\Omega_{\rm
std,-3/2}\left(\frac{T_{\text{end}}^{2}}{m_{a}M_{p}}\right)^{\frac{12-3\beta}{2\beta}},\qquad
m_{a}\lesssim m_{\rm QCD},$ (63)

where $\Omega_{\rm std,-3/2}$ is the relic density in the standard
cosmological scenario for the constant mass regime, defined in Eq. (91). From
there, it is straightforward to find the axion mass that can account for the
whole DM relic density, imposing $\Omega_{a}=\Omega_{\rm CDM}$ for our
benchmark values $\beta=6$ and $\beta=8$

$\displaystyle m_{a}\simeq\begin{dcases}2.4\times
10^{-3}\,\mbox{eV}\left(\frac{T_{\text{BBN}}}{T_{\text{end}}}\right)\theta_{i}^{2},&\qquad\beta=6\,,\\\
17\,\mbox{eV}\left(\frac{T_{\text{BBN}}}{T_{\text{end}}}\right)\theta_{i}^{8/3},&\qquad\beta=8\,,\end{dcases}$
(64)

valid for $m_{a}\lesssim m_{\rm QCD}$. We can immediately note that they are
significantly higher than the requirement on the mass for both cosmologies.
Increasing the value of $T_{\text{end}}$ does not help, because this also
increases the mass at the QCD phase transition Eq. (61). Thus, it is concluded
that an axion of constant mass can not account for the full DM relic density
for these cosmologies.

#### Thermal mass

For the range of temperatures $T_{\text{osc}}\gtrsim T_{\rm QCD}$, has to be
required on the one hand that $T_{\text{osc}}\gg T_{\text{end}}$, which leads
to a requirement on $T_{\text{end}}$

$T_{\text{end}}\ll
1.5\,\mbox{GeV}\,\left(\frac{m_{a}}{10^{-5}\,\mbox{eV}}\right)^{-1/6}.$ (65)

On the other hand, the axion mass has to be above the mass at the QCD
transition, thus $m_{a}\gtrsim m_{\rm QCD}$. To find an expression for the
relic density, we can again use the result found for region 2, Eq. (98), with
the dilution factor reduced to one, which gives

$\Omega_{a}=\Omega_{\rm std,-7/6}\left[\frac{\left(T_{\rm
QCD}^{4}\,m_{a}\,M_{P}\right)^{1/6}}{T_{\text{end}}}\right]^{7\frac{\beta-4}{\beta+8}}.$
(66)

With $\Omega_{\rm std,-7/6}$ the axion relic density in the standard cosmology
for the case of thermal axion mass, defined in Eq. (92). Then, the mass
required to produce the whole DM relic density is given by

$\displaystyle m_{a}\simeq\begin{dcases}2.4\times
10^{-3}\,\mbox{eV}\left(\frac{T_{\text{BBN}}}{T_{\text{end}}}\right)\,\theta_{i}^{2}&\qquad\beta=6,\\\
0.9\,\mbox{eV}\left(\frac{T_{\text{BBN}}}{T_{\text{end}}}\right)^{2}\,\theta_{i}^{16/7}&\qquad\beta=8.\end{dcases}$
(67)

Let us remark that these masses – especially in the case of $\beta=8$ – are
higher than the standard axion DM window. In the absence of dilution, a lower
oscillation temperature compared to the one in the standard scenario can be
reached, increasing the axion relic density, and therefore allowing to explore
higher axion masses. On the other hand, for higher values of $T_{\text{end}}$
it is possible to reach smaller masses, but recalling the requirement imposed
by Eq. (65), we find that the maximum $T_{\text{end}}$ temperature compatible
with the axion as the DM solution is

$T_{\text{BBN}}\lesssim T_{\text{end}}\lesssim 1.4~{}\mbox{GeV},$ (68)

valid for both $\beta=6$ and $\beta=8$ and for $\theta=1$.

Figure 6: Parameter space corresponding to the whole observed DM abundance,
for kination $\beta=6$ (red) and $\beta=8$ (blue), for
$0.5\leq\theta_{i}\leq\pi/\sqrt{3}$. Between the vertical dashed lines the
results for the standard cosmology are recovered. The dotted line corresponds
to $T_{\text{osc}}=T_{\text{end}}$, whereas in the grey band
$T_{\text{end}}<T_{\text{BBN}}$.

Figure 6 shows the parameter space generating the whole observed DM abundance,
for kination $\beta=6$ (red) and $\beta=8$ (blue). The widths of the bands
correspond to $0.5\leq\theta_{i}\leq\pi/\sqrt{3}$. Above the dotted black line
$T_{\text{osc}}<T_{\text{end}}$, which implies that the axion is produced in a
standard cosmology, and therefore is independent on $\omega$ and
$T_{\text{end}}$. However, below the dotted line (i.e., for
$T_{\text{osc}}>T_{\text{end}}$) bigger values for the axion mass could be
explored, especially for large values of $\beta$ and small $T_{\text{end}}$.
The white region in the left of the bands generates a DM overabundance which
overcloses the universe, whereas the region in the right can only account for
a fraction of the total DM relic density. In agreement with our analytical
estimations, from Fig. 6 it can be seen that both cosmologies are spawn
approximately over the region for $T_{\text{end}}$ found in Eq. (68). Note
that the whole parameter space can be explained from the analysis made for the
thermal mass scenario, $T_{\text{osc}}\gtrsim T_{\rm QCD}$. The apparent
change in slope of both bands around $T_{\text{end}}\sim 0.1$ GeV is due to
the change in the relativistic degrees of freedom (ignored in our analytical
treatment). A cosmology with $\beta=8$ can accommodate the DM relic density
over a higher range of masses, because the dependence on the mass is softer
than for kination and in the standard scenario, as can be found from the
expression found for the relic density, Eq. (66). In particular, for kination
($\beta=6$) the classical window is shifted to $1.6\times 10^{-6}$ eV
$\lesssim m_{a}\lesssim 10^{-2}$ eV, whereas for $\beta=8$ one can reach
$1.6\times 10^{-6}$ eV $\lesssim m_{a}\lesssim 6$ eV. Such higher values can
be explored in the near future by experiments like MADMAX [86, 41, 87, 88] and
(baby)IAXO [89, 90, 46].

##  5 Axion coupling to two photons

The coupling of axions to two photons is one of the most exploited to look for
signatures in observations and experimental searches. The interaction takes
the form

$\mathcal{L}_{a\gamma}=-\frac{1}{4}g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}=g_{a\gamma}\,a\,\vec{E}\cdot\vec{B}\,,$
(69)

where the coupling constant $g_{a\gamma}$ is model dependent and is related to
the PQ scale $f_{a}$ as

$g_{a\gamma\gamma}=\frac{\alpha}{2\pi
f_{a}}\left(\frac{E}{N}-\frac{2}{3}\frac{4+z}{1+z}\right)\simeq
10^{-13}\,\mbox{GeV}^{-1}\left(\frac{10^{10}\mbox{GeV}}{f_{a}}\right),$ (70)

where $z\equiv m_{u}/m_{d}$ and $E$ and $N$ are the electromagnetic and colour
anomalies associated with the axion anomaly. For KSVZ models $E/N=0$ [91, 92],
whereas for DFSZ models $E/N=8/3$ [22, 93]. Our findings can be mapped into
the $(g_{a\gamma},m_{a})$ to see in a comprehensive way the impact of a NSC
period over the axion relic density and the consequences for the axion search.
Figure 7 shows the parameter space of QCD axions (yellow band) and the
exclusion bounds from several astrophysical observations, cosmology based
constraints and laboratory searches (brownish coloured area). The axion DM
relic density generated in a standard cosmology is bounded between the dashed
lines. That is, the relic abundance obtained from the misalignment mechanism
for initial angles in the range $\theta_{i}\in\left[0.5,\pi/\sqrt{3}\right]$,
such that both, pre-inflationary and post-inflationary scenarios are included
without invoking a fine-tuned solution. The light blue area corresponds to the
experimental concepts with their corresponding projections of sensitivity, all
with particular emphasis to cover the QCD axion band.

Figure 7: The axion parameter space for the axion-photon coupling in the SM
and in NSC scenarios considered here, with benchmark scenarios: $\beta=2$
($\omega=-1/3$), $\beta=3$ (EMD, $\omega=0$), $\beta=6$ (kination, $\omega=1$)
and $\beta=8$ ($\omega=5/3$). The green regions correspond to current
cosmological, astrophysical and laboratory constraints. The light purple area
are projections for experimental searches, as described in the text. For the
standard cosmology window of axion mass we have assumed contributions to the
relic abundance from the misalignment mechanism only, with the initial angle
on the range $\theta_{i}\in\left[0.5,\pi/\sqrt{3}\right]$, in order to include
both, the pre- and post-inflationary scenarios. Figure adapted from [94].

We have delineated with brackets the range of masses – equivalently, the range
in the scale $f_{a}$ – that a NSC with the given equations of state could
access. Our benchmark values are $\beta=2$, equivalent to $\omega=-1/3$,
$\beta=3$, corresponding to an EMD,  i.e. $\omega=0$, $\beta=6$, known as a
kination regime, with $\omega=1$ and finally $\beta=8$, which corresponds to
an equation of state with $\omega=5/3$. We have fixed
$T_{\text{end}}=T_{\text{BBN}}=4$ MeV and we scan in the equality temperature
$T_{\text{eq}}$ according to the results obtained in subsection 4.1 for
$\beta<4$ and the scan is in $T_{\text{end}}$ for $\beta>4$, as shown in
subsection 4.2, all summarised in Figs. 5 and 6. There, it was found that for
a EMD cosmology, the same mass range as the standard axion CDM window can be
reached – for cosmologies with very low $T_{\text{eq}}$ – and then down to
$m_{a}\gtrsim 4\times 10^{-9}$ eV, for cosmologies with $T_{\text{eq}}\gtrsim
10^{3}~{}$ GeV. The latter is translated into a $g_{a\gamma}\sim 10^{-19}$
GeV-1. Meanwhile, for cosmologies with $\beta=2$ it is possible to get into
even smaller masses, thanks to the higher dilution of the axion relic density.
According to our findings, the smallest mass attainable for this cosmology is
of $m_{a}\gtrsim 10^{-10}~{}$eV, for $T_{\text{eq}}\gtrsim 0.1$ GeV, which
corresponds to a coupling of the order of $g_{a\gamma}\sim 10^{-20}$ GeV-1.
Both scenarios will be in the reach of the ABRACADABRA, KLASH and the next
generation of ADMX.

On the other hand, cosmologies with $\beta>4$ open up the axion DM window to
the right of the standard cosmological scenario, because there is no dilution
of the energy density, but a smaller $T_{\text{osc}}$, delaying the onset of
DM production, increasing the relic density. These cosmologies will be on the
reach of experiments such as ADMX, CULTASK, MADMAX, and ORGAN. For $\beta=6$
we have found the highest axion mass that can accommodate the whole DM relic
density is $m_{a}\lesssim 10^{-2}$ eV, with a temperature $T_{\text{end}}$
very close to the BBN value. The corresponding scale for the coupling to
photons is $g_{a\gamma}\sim 10^{-13}$ GeV-1. For $\beta=8$ and the same scale
on $T_{\text{end}}$, it is even possible to reach $m_{a}\lesssim 6$ eV, which
gives $g_{a\gamma}\sim 10^{-9}$ GeV-1. Cosmologies with $\beta=8$ will be
already probed with the future helioscope experiment IAXO. We also note that
cosmologies with $\beta=8$ and $T_{\text{end}}\lesssim 10^{-2}$ GeV conflict
with the horizontal branch bound and other astrophysical observations [95, 96,
97]. Let us note that the SN 1987A bound on the QCD axion rules out the
parameter space 15 meV $\lesssim m_{a}\lesssim 1$ keV [98, 99]. However, in
Ref. [100] has been pointed out that pion-related processes have been
previously underestimated and a better quantitative understanding of their
contribution is needed. Thus, it still seems worth to explore this region of
parameter space independently with experimentally well-known and controlled
sources. For an exhaustive review on the experimental programme of direct
detection searches and expected timescales, we refer to Refs. [48, 47, 49].
Finally, let us remark that for all the equations of state studied here, it is
possible to partially cover the standard axion DM mass range, although this
happens for a very short-lived NSC.

##  6 Conclusions

Throughout this paper we have made a comprehensive and exhaustive study of the
implications of a nonstandard cosmological (NSC) period – somewhere prior to
the BBN era – on the onset of the axion field oscillation and the
corresponding dark matter (DM) relic density generated in the so-called
process of misalignment. We have considered a generically new field $\phi$
whose energy content gets to dominate the expansion of the universe, until it
decays into SM radiation for $\omega<1/3$ or dilutes away in the case of
$\omega>1/3$, with $\omega$ the equation-of-state parameter. As examples, we
have considered benchmark cosmological scenarios that have been motivated
previously in the literature: $\omega=-1/3$, $\omega=0$ (early matter
domination), $\omega=1$ (kination) and $\omega=5/3$. Although we have also
found analytical expressions for a generic equation of state.

In the case of $\omega<1/3$, we have divided the period of NSC in three
regions: i) a radiation driven period before the dominance of $\phi$, ii)
domination of $\phi$ prior to its decay and iii) the era of decay and
corresponding entropy injection into the SM plasma. For all these regions, we
have found the corresponding oscillation temperature, DM relic density and the
axion mass as a function of the NSC parameters that can account for the whole
DM density observed today, distinguishing the regime of temperatures above the
QCD phase transition (where the mass gets important thermal effects) and below
(where it can be considered constant). Our findings have shown that in
principle, a wide range of cosmologies, with their characteristic parameters:
$\omega$, $T_{\text{eq}}$ (temperature that sets the beginning of the
dominance of $\phi$ over radiation) and $T_{\text{end}}$ (temperature where
$\phi$ has mostly decayed away) can lead to the right production of axion DM
in a wide range of masses below the standard axion window. For the case of
$\omega=-1/3$, we have found that the smallest possible axion mass (without a
tuning of the initial misalignment angle) is around $m_{a}\sim 10^{-10}~{}$eV
and can be achieved for cosmologies with $T_{\text{eq}}\gtrsim 0.1$ GeV,
$T_{\text{end}}\sim 4$ MeV. In that case, the oscillation of the axion field
happens either before or during the epoch of decay of $\phi$ and below the QCD
phase transition, so the axion mass can be considered constant. Something
similar is found for a NSC with matter domination $\omega=0$, where the
smallest axion mass is around $m_{a}\sim 10^{-8}$ eV, for cosmologies with
$T_{\text{eq}}\gtrsim 10^{2}$ GeV, $T_{\text{end}}\sim 4$ MeV, as shown in
Fig. 5.

For cosmologies with $\omega>1/3$, we have considered a stable $\phi$, with an
energy density that fades away due to the expansion of the universe. In that
case there is no dilution, and we find that the axion relic density is higher
than in the standard cosmological scenario, due to a lower oscillation
temperature – consequence of a bigger Hubble parameter – than in the standard
radiation-dominated scenario. Thus, higher axion masses than in a radiation-
dominated universe can lead to the explanation of the DM paradigm. For the
equations of state $\omega=1$ and $\omega=5/3$, we have found that cosmologies
where the extinction of $\phi$ takes place near the BBN epoch can lead to
higher axion masses, $m_{a}\sim 10^{-2}$ eV and $6$ eV, respectively. In this
case, the axion oscillation happens above the QCD phase transition and
therefore the mass receives thermal corrections, as shown in Fig. 6.

Concerning the impact of NSC on axion searches, we have projected our results
into the exclusion plot of the axion coupling to two photons, the most
exploited interaction to search for these particles. We have highlighted that
a handful of NSCs can spawn the axion DM window in the range
$10^{-20}\,\mbox{GeV}\lesssim g_{a\gamma}$, all the way up to $g_{a\gamma}\sim
10^{-9}$ GeV, see Fig. 7, the upper part being in slight conflict with
astrophysical searches. Remarkably, most of the parameter space is expected to
be probed by experimental prospects in the near future. These findings are
expected to have a positive impact on the motivation to search with
laboratory-based experiments for axion DM outside the so-called classical
window of standard cosmology. Then, a potential discovery of the QCD axion
will shed some light into the cosmological history of the early universe,
prior to the BBN era.

Let us finally stress that we have only considered the production of DM
through the misalignment mechanism. In the case where the PQS is broken after
the inflationary epoch, the potential contribution from topological defects
such as strings and domain walls could move the mass ranges found here.
Nonetheless, given the still ongoing uncertainty on the fraction contributed
by these mechanisms, we have only focused on the misalignment production.

## Acknowledgments

We are thankful to Luca Visinelli for useful comments and feedback. PA and MV
acknowledge support from FONDECYT project 1161150 and Proyecto POSTDOC-DICYT,
042131AR-AYUDANTE, VRIDEI. PA acknowledges support from DICYT 042131AR
project. NB received funding from the Spanish FEDER/MCIU-AEI under grant
FPA2017-84543-P, and the Patrimonio Autónomo - Fondo Nacional de
Financiamiento para la Ciencia, la Tecnología y la Innovación Francisco José
de Caldas (MinCiencias - Colombia) grant 80740-465-2020. DK and LR are
supported in part by the National Science Centre, Poland, research grant No.
2015/18/A/ST2/00748. LR is also supported by the project AstroCeNT: Particle
Astrophysics Science and Technology Centre, carried out within the
International Research Agendas programme of the Foundation for Polish Science
financed by the European Union under the European Regional Development Fund.
This project has received funding /support from the European Union’s Horizon
2020 research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 860881-HIDDeN. PA and MV gratefully thank AstroCeNT for their
hospitality, where part of this work was developed. PA is deeply thankful to
Laura Pejcha for her support and company throughout this project.

##  Appendix A Solutions for the $\phi-$radiation system

Let us start by solving the evolution Eqs. (19) and (20) during the dominance
of $\phi$, taking into account the decays. The Hubble parameter can be written
as $H\simeq\sqrt{\rho_{\phi}/(3M_{P}^{2})}$. By changing variables from time
to scale factor, and defining $H_{\text{eq}}\equiv H(T_{\text{eq}})$,
$u(r)=\frac{\rho_{\phi}}{\rho_{\rm eq}}r^{\beta}$, and $r=R/R_{\text{eq}}$,
the evolution for the $\phi$ field now reads

$\frac{du}{u^{1/2}}=-\frac{\Gamma_{\phi}}{H_{\text{eq}}}\,r^{\beta/2-1}\,dr.$
(71)

Whose solution to first order in $\Gamma_{\phi}/H_{\text{eq}}$ is

$\rho_{\phi}(R)=\rho_{\text{eq}}\left[\frac{\Gamma_{\phi}}{\beta\,H_{\text{eq}}}-\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/2}\left(1+\frac{\Gamma_{\phi}}{H_{\text{eq}}\,\beta}\right)\right]^{2}.$
(72)

Replacing into Eq. (20), using $s=(\rho+p)/T$, we find

$\rho_{R}(R)=\rho_{\text{eq}}\left[\left(\frac{R_{\text{eq}}}{R}\right)^{4}+\frac{2\Gamma_{\phi}}{(8-\beta)H_{\text{eq}}}\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/2}-\frac{2\Gamma_{\phi}\,}{(8-\beta)H_{\text{eq}}}\left(\frac{R_{\text{eq}}}{R}\right)^{4}\right].$
(73)

For $\beta<4$, the universe is radiation dominated before $R_{\text{eq}}$. It
is possible to get an educated guess for $R_{\text{c}}$ and $R_{\text{end}}$
by considering that when the decays start to influence the SM radiation, both
terms in $\rho_{R}$ should be of the same order, meaning

$\left(\frac{R_{\text{eq}}}{R_{\text{c}}}\right)^{4}\simeq\frac{2\Gamma_{\phi}}{(8-\beta)H_{\text{eq}}}\left(\frac{R_{\text{eq}}}{R_{\text{c}}}\right)^{\beta/2},$
(74)

from here, we get

$R_{\text{c}}=R_{\text{eq}}\left(\frac{(8-\beta)}{2}\left(\frac{T_{\text{eq}}}{T_{\text{end}}}\right)^{2}\right)^{\frac{2}{8-\beta}}.$
(75)

Similarly, the decays of $\phi$ will stop, when the two terms inside the
brackets of Eq. (72) become comparable, thus

$\left(\frac{R_{\text{eq}}}{R_{\text{end}}}\right)^{\beta}\simeq\frac{2\Gamma_{\phi}}{\beta\,H_{\text{eq}}}\left(\frac{R_{\text{eq}}}{R_{\text{end}}}\right)^{\beta/2},$
(76)

from where we find

$R_{\text{end}}=R_{\text{eq}}\left(\frac{\beta}{2}\left(\frac{T_{\text{eq}}}{T_{\text{end}}}\right)^{2}\right)^{2/\beta}.$
(77)

Furthermore, the corresponding temperatures can be obtained from

$\displaystyle\rho_{\text{eq}}\left(\frac{R_{\text{eq}}}{R_{\text{c}}}\right)^{4}\simeq\rho_{R}(T_{\text{c}}),\quad\mbox{and}\quad{{\frac{\Gamma_{\phi}}{H_{\text{eq}}}=\frac{T_{\text{end}}^{2}}{T_{\text{eq}}^{2}},}}$
(78)

which result in

$T_{\text{c}}\simeq
T_{\text{eq}}\left(\frac{2\Gamma_{\phi}\,}{H_{\text{eq}}(8-\beta)}\right)^{\frac{2}{8-\beta}}=T_{\text{eq}}\left(\frac{4}{(8-\beta)^{2}}\frac{T_{\text{end}}^{4}}{T_{\text{eq}}^{4}}\right)^{\frac{1}{8-\beta}}.$
(79)

From this analysis it can be also extracted that deep during the $\phi$
domination, the relation between temperature and scale factor is

$T\simeq
T_{\text{eq}}\left[\frac{2}{8-\beta}\frac{T_{\text{end}}^{2}}{T_{\text{eq}}^{2}}\right]^{1/4}\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/8}.$
(80)

We emphasise the equation above tells us that for smaller $\beta$ the
temperature decreases more slowly with the universe expansion and therefore it
implies that the universe is bigger when reaching the temperature
$T_{\text{end}}$. Let us also remark that the range of Eqs. (25)-(30) are only
valid if the fluid dominates strongly over radiation. There can be certain
cosmologies where this is not the case - and still have an impact on the axion
relic density - and thus, the above referred expressions have to be handled
carefully.

The case $\beta=0$ should be treated separately. In this case, we start with
Eq. (19) and defining $u=\sqrt{\rho_{\phi}/\rho_{\text{eq}}}$ and
$r=\ln(R/R_{\text{eq}})$. Thus, the equation to solve is

$du=-\frac{\Gamma_{\phi}}{2H_{\text{eq}}}dr,$ (81)

whose solution is given by

$\rho_{\phi}(R)\simeq\rho_{\text{eq}}\left(1+\frac{\Gamma_{\phi}}{H_{\text{eq}}}\ln\left(\frac{R_{\text{eq}}}{R}\right)\right).$
(82)

By replacing into the radiation equation, we get the evolution is

$\rho_{R}(R)\simeq\rho_{R_{\text{eq}}}\left(\left(\frac{R_{\text{eq}}}{R}\right)^{4}+\frac{\Gamma_{\phi}}{4H_{\text{eq}}}\right).$
(83)

We can find the moment when the radiation energy density starts to be
influenced by the $\phi$ field

$R_{\text{c}}=R_{\text{eq}}\left(\frac{4H_{\text{eq}}}{\Gamma_{\phi}}\right)^{1/4}.$
(84)

Additionally the moment when the decays of the $\phi$ field are significant is
represented by

$R_{\text{end}}=R_{\text{eq}}\times\exp\left(-\frac{H_{\text{eq}}}{\Gamma_{\phi}}\right).$
(85)

##  Appendix B Axion oscillations during $\phi$ domination

###  B.1 Oscillations in Region 1

To analyse this epoch, first we find the oscillation temperature. Since for
this period it is assumed the universe is radiation dominated, the Hubble
parameter is $H=\sqrt{\rho_{R}/(3M_{P}^{2})}$, so the oscillation temperature
is the usual found for a QCD axion, depicted in Fig. 1. There are two regimes
for the oscillation temperature, depending on whether the effects of
temperature on the axion mass, are important or not, Eq. (7), and they are the
same as the standard cosmological scenario, namely

$T_{\text{osc}}\simeq\begin{dcases}\left(m_{a}\,M_{P}\right)^{1/2},&T_{\text{osc}}\lesssim
T_{\rm QCD}.\\\ \left(m_{a}M_{P}T_{\rm
QCD}^{4}\right)^{1/6},&T_{\text{osc}}\gtrsim T_{\rm QCD}\end{dcases}$ (86)

These two temperatures intersect at the QCD epoch, which can be expressed as a
particular mass of the axion, which is the same as the standard cosmological
scenario, found in Eq. (16), but we rewrite here, for completeness, this time
forgetting about the degrees of freedom, as we are only interested in the
order of magnitude, so

$\tilde{m}_{a}(T_{\rm QCD})=m_{R_{1}}\simeq\frac{T_{\rm QCD}^{2}}{M_{P}}.$
(87)

The condition for $R_{\text{osc}}$ to occur in this regime can be re-written
to be $T_{\rm{osc}}\gg T_{\text{eq}}$ and by using the expressions for
$T_{\text{osc}}$ we can write it as a condition for the axion mass

$m_{a}\gg\begin{dcases}\frac{T_{\text{eq}}^{2}\,}{M_{P}},&m_{a}\lesssim
m_{R_{1}},\\\ \frac{T_{\text{eq}}^{6}}{M_{P}\,T_{\rm QCD}^{4}},&m_{a}\gtrsim
m_{R_{1}},\end{dcases}$ (88)

where the first equation correspond to the limit of the axion mass in the case
that the temperature has significant effects on the mass of the axion (before
QCD transition) and the second is the limit after the QCD transition.

The relic abundance, if the oscillation happens in this region, is given by
$\Omega_{a}(T_{0})=\Omega_{a}^{\rm std}(T_{0})\gamma_{R_{1}}$, such that
$\Omega_{a}^{\rm std}$ is the axion relic energy density today in the standard
cosmology, and $\gamma_{R_{1}}=S_{\rm osc}/S_{\rm end}$ the entropy dilution
factor. To find the latter, we take

$\displaystyle\gamma_{R_{1}}$ $\displaystyle=$ $\displaystyle\frac{S_{\rm
osc}}{S_{\rm end}}=\frac{S_{\rm c}}{S_{\rm
end}}=\left(\frac{T_{\text{c}}\,R_{\text{c}}}{T_{\text{end}}\,R_{\text{end}}}\right)^{3},$
$\displaystyle=$
$\displaystyle\left(\frac{\beta^{2}}{4}\left(\frac{T_{\text{eq}}}{T_{\text{end}}}\right)^{4-\beta}\right)^{-\frac{3}{\beta}}.$

Now, replacing the dilution factor in the axion relic density we find

$\displaystyle\Omega_{R_{1}}=\begin{dcases}\Omega_{\rm
std,-3/2}\left(\frac{\beta^{2}}{4}\right)^{-3/\beta}\left(\frac{T_{\text{end}}}{T_{\text{eq}}}\right)^{12/\beta-3},&\mbox{for}\,\,\,m_{a}\lesssim
m_{R_{1}}\\\ \Omega_{\rm
std,-7/6}\left(\frac{\beta^{2}}{4}\right)^{-3/\beta}\left(\frac{T_{\text{end}}}{T_{\text{eq}}}\right)^{12/\beta-3},&\mbox{for}\,\,\,m_{a}\gtrsim
m_{R_{1}},\end{dcases}$ (90)

where we have defined the axion relic density for a standard cosmology, when
temperature effects on the mass are unimportant and important, respectively,
as

$\displaystyle\Omega_{\rm std,-3/2}\equiv 5\times
10^{-11}\left(\frac{m_{a}}{1\,\mbox{eV}}\right)^{-3/2}\,\theta_{i}^{2},$ (91)
$\displaystyle\Omega_{\rm std,-7/6}\equiv 2.4\times
10^{-7}\left(\frac{m_{a}}{1\,\mbox{eV}}\right)^{-7/6}\,\theta_{i}^{2}.$ (92)

###  B.2 Oscillations in Region 2

We start by writing the Hubble parameter as

$H\simeq\sqrt{\frac{\rho_{\phi}}{3M_{P}^{2}}}\simeq
H_{\text{eq}}\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/2}=H_{\text{eq}}\left(\frac{T_{\text{osc}}}{T_{\text{eq}}}\right)^{\beta/2},$
(93)

where $H_{\text{eq}}\equiv H(R_{\text{eq}})$. We find the oscillation
temperature for both cases, where the temperatures effects on the mass are and
are not important, as

$\displaystyle
T_{\text{osc}}^{R_{2}}=\begin{dcases}T_{\text{eq}}\left(\frac{M_{P}m_{a}}{T_{\text{eq}}^{2}}\right)^{2/\beta}&m_{a}\lesssim
m_{R_{2}},\\\ \left(T_{\text{eq}}^{\frac{\beta-4}{2}}M_{P}\,m_{a}\,T_{\rm
QCD}^{4}\right)^{\frac{2}{\beta+8}}&m_{a}\gtrsim m_{R_{2}}.\end{dcases}$ (94)

These two temperatures intersect at the axion mass at the QCD temperature,
$\tilde{m}_{a}(T_{\rm QCD})$, in this case given by

$m_{R_{2}}=\left(\frac{T_{\rm
QCD}}{T_{\text{eq}}}\right)^{\beta/2}\,\frac{T_{\text{eq}}^{2}}{M_{P}}.$ (95)

By requiring $T_{\text{c}}\ll T_{\text{osc}}\ll T_{\text{eq}}$ we find the
mass range of the axion to have the oscillation happening in this era. For the
case of constant axion mass, the requirement is $T_{\text{c}}\ll
T_{\text{osc}}\lesssim T_{\rm QCD}$ so we find

$\displaystyle m_{R_{2}}\ll m_{a}\ll\frac{T_{\text{eq}}^{6}}{M_{P}\,T_{\rm
QCD}^{4}}.$ (96)

While for the case the oscillation happens before the QCD phase transition, it
is found

$\frac{T_{\text{eq}}^{2}}{M_{P}}\left(\frac{2\,T_{\text{end}}^{2}}{T_{\text{eq}}^{2}(8-\beta)}\right)^{\frac{\beta}{8-\beta}}\ll
m_{a}\ll m_{R_{2}}.$ (97)

The axion abundance can be found by taking the expression for the energy
density in the NSC, Eq. (33), where $S_{\rm osc}/S_{\rm end}$ is the entropy
dilution, which is the same found in Eq. (34). We get

$\Omega_{R_{2}}\simeq\begin{dcases}\Omega_{\rm
std,-3/2}\left(\frac{T_{\text{end}}^{2}}{m_{a}M_{P}}\right)^{\frac{3}{2\beta}(4-\beta)}&\text{for
}m_{a}\lesssim m_{R_{2}}\,,\\\ \Omega_{\rm
std,-7/6}\left(\frac{\beta^{2}}{4}\right)^{-3/\beta}\left(\frac{T_{\text{end}}}{T_{\text{eq}}}\right)^{\frac{3}{\beta}(4-\beta)}\left[\frac{T_{\text{eq}}^{7}}{\left(T_{\rm
QCD}^{4}\,m_{a}M_{P}\right)^{7/6}}\right]^{\frac{4-\beta}{\beta+8}}&\text{for
}m_{a}\gtrsim m_{R_{2}}\,.\end{dcases}$ (98)

Note that in the case of a constant axion mass, the dependence on
$T_{\text{eq}}$ cancels.

###  B.3 Oscillations in Region 3

We will assume $R_{\text{osc}}$ is still far from $R_{\text{end}}$, such that
the Hubble parameter, for our purposes, can still be considered as

$H\sim\sqrt{\frac{\rho_{\phi}}{3\,M_{P}^{2}}}=H_{\text{eq}}\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/2}.$
(99)

Nonetheless, the decays of $\phi$ are affecting radiation, such that its
energy density it is better described as

$\rho_{R}\sim\rho_{\text{eq}}\,\frac{2\Gamma_{\phi}}{(8-\beta)H_{\text{eq}}}\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/2}.$
(100)

Now we can equate the above expression with the energy density as a function
of the temperature, the well known $\rho_{R}=(\pi^{2}/30)g_{\star}(T)T^{4}$,
such that we can obtain the relation between scale factor and temperature in
this region to be

$\left(\frac{R_{\text{eq}}}{R}\right)^{\beta/2}=\left(\frac{8-\beta}{2}\right)\frac{T^{4}}{T_{\text{end}}^{2}T_{\text{eq}}^{2}}.$
(101)

Replacing back into the Hubble parameter, we find777Dropping the degrees of
freedom.

$H(T)=H_{\text{eq}}\left(\frac{8-\beta}{2}\right)\frac{T^{4}}{T_{\text{end}}^{2}T_{\text{eq}}^{2}}.$
(102)

As in the previous regions, we can find the oscillation temperature

$\displaystyle
T_{\text{osc}}^{R_{3}}=\begin{dcases}\left(\frac{2m_{a}M_{P}\,T_{\text{end}}^{2}}{8-\beta}\right)^{1/4}&\mbox{for}\,\,\,m_{a}\lesssim
m_{R_{3}}\\\ \left(\frac{2m_{a}M_{P}\,T_{\rm
QCD}^{4}\,T_{\text{end}}^{2}}{8-\beta}\right)^{1/8}&\mbox{for}\,\,\,m_{a}\gtrsim
m_{R_{3}}.\end{dcases}$ (103)

So for this region the dependence of $T_{\text{osc}}$ on $\beta$ is quite
mild, and on $T_{\text{eq}}$ nonexistent. The intersection between these
temperatures gives the moment of the QCD phase transition, which can be used
to obtain the mass of the axion at this point. We get

$m_{R_{3}}=\frac{(8-\beta)}{2}\frac{T_{\rm
QCD}^{4}}{M_{P}\,T_{\text{end}}^{2}}\simeq 2\times
10^{-13}\,\mbox{eV}\left(\frac{\mbox{GeV}}{T_{\text{end}}}\right)^{2}.$ (104)

We now can write down the mass ranges expected for the oscillation of the
axion field to happen in region 3. Here it is better to work with the
requirement $R_{\text{c}}\ll R_{\rm osc}\ll R_{\text{end}}$, but it can be put
into temperatures by using Eq. (101). Firstly, let us assume the oscillation
happens completely before the QCD transition, where temperature effects on the
axion mass are important. In that case

$m_{R_{3}}\ll m_{a}\ll\frac{T_{\text{eq}}^{6}}{M_{P}\,T_{\rm
QCD}^{4}}\left(\frac{2\,T_{\text{end}}^{2}}{(8-\beta)T_{\text{eq}}^{2}}\right)^{\frac{8+\beta}{8-\beta}},$
(105)

On the other hand, for the constant axion mass regime, the requirement maps
into a mass range of as

$\frac{2}{\beta}\frac{T_{\text{end}}^{2}}{M_{P}}\ll
m_{a}\ll\frac{T_{\text{eq}}^{2}}{M_{P}}\left(\frac{2\,T_{\text{end}}^{2}}{(8-\beta)T_{\text{eq}}^{2}}\right)^{\frac{\beta}{8-\beta}}.$
(106)

The dilution factor, $S_{\rm osc}/S_{\rm end}$, is different in the previous
regions, because the axion is produced in an era where the SM entropy is
increasing. We write

$\gamma_{R_{3}}=\frac{S_{\rm osc}}{S_{\rm
end}}=\left(\frac{T_{\text{osc}}R_{\text{osc}}}{T_{\text{end}}R_{\text{end}}}\right)^{3}.$
(107)

From Eqs. (101) and (28) we can easily find the ratio
$R_{\text{osc}}/R_{\text{end}}$ to be

$\frac{R_{\text{osc}}}{R_{\text{end}}}=\left[\frac{4}{\beta(8-\beta)}\left(\frac{T_{\text{end}}}{T_{\text{osc}}}\right)^{4}\right]^{2/\beta},$
(108)

and then the dilution factor is

$\gamma_{R_{3}}=\left(\frac{4}{\beta(8-\beta)}\right)^{6/\beta}\left(\frac{T_{\text{end}}}{T_{\text{osc}}}\right)^{24/\beta-3}.$
(109)

Now the axion abundance is easily found

$\displaystyle\Omega_{R_{3}}=\begin{dcases}\Omega_{\rm
std,-3/2}\left(\frac{2}{\beta}\right)^{6/\beta}\left(\frac{T_{\text{end}}^{2}}{m_{a}M_{P}}\right)^{\frac{3}{2\beta}(4-\beta)}&\mbox{for}\,\,m_{a}\lesssim
m_{R_{3}}\,,\\\ \Omega_{\rm
std,-7/6}\,(8-\beta)^{\frac{\beta+6}{2\beta}}\left(\frac{T_{\text{end}}^{6}}{T_{\rm
QCD}^{4}m_{a}M_{P}}\right)^{3/\beta-2/3}&\mbox{for}\,\,m_{a}\gtrsim
m_{R_{3}}\,.\end{dcases}$ (110)

Both expressions are independent on $T_{\text{eq}}$, showing that for this
region the main parameters are $T_{\text{end}}$ and $\beta$.

##  Appendix C Adiabatic invariant and the anharmonic factor

In oscillatory systems with varying period, the energy is not conserved, and
it is usually useful to define an “adiabatic invariant”,

$J\equiv\oint p\ d\theta\;,$ (111)

which is an approximate constant of motion.

To calculate the adiabatic invariant of the axion, we note that the
Hamiltonian that results in the EOM of Eq. (9) is

$\mathcal{H}=\dfrac{1}{2}\dfrac{p^{2}}{f_{a}^{2}\ R^{3}}+V(\theta)\ R^{3}\;,$
(112)

with

$\displaystyle p=f_{a}^{2}\ R^{3}\ \dot{\theta}$ (113) $\displaystyle
V(\theta)=\tilde{m}_{a}^{2}(t)f_{a}^{2}(1-\cos\theta)\;.$ (114)

Notice that the explicit time dependence comes from $R$ and $\tilde{m}_{a}$.
That is, the Hamiltonian varies slowly if
$\dot{\tilde{m}}_{a}/\tilde{m}_{a}\ll\tilde{m}_{a}$ and $H\ll\tilde{m}_{a}$,
which are the adiabatic conditions. When these conditions are met, the
adiabatic invariant of this system becomes

$J=2\sqrt{2}\ f_{a}^{2}\ \tilde{m}_{a}(t)\,R^{3}\ \int_{-\theta_{\rm
max}}^{\theta_{\rm max}}\sqrt{\cos\theta-\cos\theta_{\rm max}}\ d\theta\;,$
(115)

where we have defined $\theta_{\rm max}$ the maximum $\theta$ during its
oscillation, which corresponds to $p=0$. Moreover, we have assumed negligible
change of $\tilde{m}_{a}$ and $R$ during one period, due to the adiabatic
conditions. For the sake of consistency with the literature [101, 11], we
define $I\equiv J/(\pi\ f_{a}^{2})$; a rescaled adiabatic invariant

$I\equiv R^{3}\ \tilde{m}_{a}(t)\ \theta_{\rm max}^{2}\,f(\theta_{\rm
max})\;,$ (116)

where

$f(\theta_{\rm max})=\dfrac{2\sqrt{2}}{\pi\theta_{\rm
max}^{2}}\int_{-\theta_{\rm max}}^{\theta_{\rm
max}}d\theta\sqrt{\cos\theta-\cos\theta_{\rm max}}\;,$ (117)

is called the anharmonic factor, with $0.5\lesssim f(\theta_{\rm max})\leq 1$,
as shown in Fig. 8.

Figure 8: The anharmonic factor for $0\leq\theta_{\rm max}<\pi$.

The adiabatic invariant allows us to calculate the maximum value of the angle
$\theta$ at late times from its corresponding value at some point after the
adiabatic conditions where met (at some $T\ll T_{\text{osc}}$). Numerically,
this means that one can stop the integration once the axion starts to evolve
adiabatically – at some temperature $T_{\rm ad}$ with $\theta=\theta_{\rm
ad}$. Then, since $J$ is conserved between $T_{\rm ad}$ and $T_{0}$, the axion
angle at $T=T_{0}$ becomes

$\theta_{0}^{2}=\dfrac{s_{0}}{s_{\rm ad}}\dfrac{\tilde{m}_{a,{\rm
ad}}}{m_{a}}\theta_{\rm ad}^{2}\ f(\theta_{\rm ad})\times\dfrac{S_{\rm
ad}}{S_{0}}\;,$ (118)

where the subscript “${\rm ad}$” indicates that the corresponding quantity is
evaluated at $T=T_{\rm ad}$. Therefore, the energy density of the axion
becomes

$\rho_{a}(T_{0})=\dfrac{1}{2}\ \dfrac{s_{0}}{s_{\rm ad}}\ f_{a}\ m_{a}\
\tilde{m}_{a,{\rm ad}}\ \theta_{\rm ad}^{2}\ f(\theta_{\rm
ad})\times\dfrac{S_{\rm ad}}{S_{0}}\;.$ (119)

Although this form of the axion energy density is similar to Eq. (33) – with
$T_{\text{osc}}\to T_{\rm ad}$, multiplied with the anharmonic factor in Eq.
(117). It is valid for any value of $\theta_{i}$, as it takes into account the
exact form of the potential of the axion field.

## References

  * [1] R. Peccei and H.R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * [2] S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * [3] F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279.
  * [4] J. Preskill, M.B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 127.
  * [5] M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137.
  * [6] L. Abbott and P. Sikivie, _A Cosmological Bound on the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 133.
  * [7] R. Crewther, _Chirality Selection Rules and the U(1) Problem_ , _Phys. Lett. B_ 70 (1977) 349.
  * [8] P. Di Vecchia and G. Veneziano, _Chiral Dynamics in the Large N Limit_ , _Nucl. Phys. B_ 171 (1980) 253.
  * [9] M. Gorghetto and G. Villadoro, _Topological Susceptibility and QCD Axion Mass: QED and NNLO corrections_ , _JHEP_ 03 (2019) 033 [1812.01008].
  * [10] S. Borsanyi et al., _Calculation of the axion mass based on high-temperature lattice quantum chromodynamics_ , _Nature_ 539 (2016) 69 [1606.07494].
  * [11] K.J. Bae, J.-H. Huh and J.E. Kim, _Update of axion CDM energy_ , _JCAP_ 09 (2008) 005 [0806.0497].
  * [12] O. Wantz and E. Shellard, _Axion Cosmology Revisited_ , _Phys. Rev. D_ 82 (2010) 123508 [0910.1066].
  * [13] G. Ballesteros, J. Redondo, A. Ringwald and C. Tamarit, _Standard Model-axion-seesaw-Higgs portal inflation. Five problems of particle physics and cosmology solved in one stroke_ , _JCAP_ 08 (2017) 001 [1610.01639].
  * [14] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [15] M.P. Hertzberg, M. Tegmark and F. Wilczek, _Axion Cosmology and the Energy Scale of Inflation_ , _Phys. Rev. D_ 78 (2008) 083507 [0807.1726].
  * [16] M. Beltran, J. Garcia-Bellido and J. Lesgourgues, _Isocurvature bounds on axions revisited_ , _Phys. Rev. D_ 75 (2007) 103507 [hep-ph/0606107].
  * [17] C. Hagmann, S. Chang and P. Sikivie, _Axion radiation from strings_ , _Phys. Rev. D_ 63 (2001) 125018 [hep-ph/0012361].
  * [18] O. Wantz and E. Shellard, _The Topological susceptibility from grand canonical simulations in the interacting instanton liquid model: Chiral phase transition and axion mass_ , _Nucl. Phys. B_ 829 (2010) 110 [0908.0324].
  * [19] T. Hiramatsu, M. Kawasaki, T. Sekiguchi, M. Yamaguchi and J. Yokoyama, _Improved estimation of radiated axions from cosmological axionic strings_ , _Phys. Rev. D_ 83 (2011) 123531 [1012.5502].
  * [20] M. Kawasaki, K. Saikawa and T. Sekiguchi, _Axion dark matter from topological defects_ , _Phys. Rev. D_ 91 (2015) 065014 [1412.0789].
  * [21] M. Dine, W. Fischler and M. Srednicki, _A simple solution to the strong cp problem with a harmless axion_ , _Physics letters B_ 104 (1981) 199.
  * [22] A.R. Zhitnitsky, _On Possible Suppression of the Axion Hadron Interactions. (In Russian)_ , _Sov. J. Nucl. Phys._ 31 (1980) 260.
  * [23] J.E. Kim, _Weak-interaction singlet and strong cp invariance_ , _Physical Review Letters_ 43 (1979) 103.
  * [24] M.A. Shifman, A. Vainshtein and V.I. Zakharov, _Can confinement ensure natural cp invariance of strong interactions?_ , _Nuclear Physics B_ 166 (1980) 493.
  * [25] L. Di Luzio, F. Mescia and E. Nardi, _Redefining the Axion Window_ , _Phys. Rev. Lett._ 118 (2017) 031801 [1610.07593].
  * [26] J.R. Ellis and K.A. Olive, _Constraints on Light Particles From Supernova Sn1987a_ , _Phys. Lett. B_ 193 (1987) 525.
  * [27] G. Raffelt and D. Seckel, _Bounds on Exotic Particle Interactions from SN 1987a_ , _Phys. Rev. Lett._ 60 (1988) 1793.
  * [28] M.S. Turner, _Axions from SN 1987a_ , _Phys. Rev. Lett._ 60 (1988) 1797.
  * [29] P. Sikivie, _Experimental Tests of the Invisible Axion_ , _Phys. Rev. Lett._ 51 (1983) 1415.
  * [30] ADMX collaboration, _A SQUID-based microwave cavity search for dark-matter axions_ , _Phys. Rev. Lett._ 104 (2010) 041301 [0910.5914].
  * [31] ADMX collaboration, _Extended Search for the Invisible Axion with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 124 (2020) 101303 [1910.08638].
  * [32] HAYSTAC collaboration, _Results from phase 1 of the HAYSTAC microwave cavity axion experiment_ , _Phys. Rev. D_ 97 (2018) 092001 [1803.03690].
  * [33] CAPP collaboration, _First Results from an Axion Haloscope at CAPP around 10.7 $\mu$eV_, _Phys. Rev. Lett._ 126 (2021) 191802 [2012.10764].
  * [34] S. Lee, S. Ahn, J. Choi, B.R. Ko and Y.K. Semertzidis, _Axion Dark Matter Search around 6.7 $\mu$eV_, _Phys. Rev. Lett._ 124 (2020) 101802 [2001.05102].
  * [35] D. Alesini, D. Babusci, D. Di Gioacchino, C. Gatti, G. Lamanna and C. Ligi, _The KLASH Proposal_ , 1707.06010.
  * [36] B.T. McAllister, G. Flower, E.N. Ivanov, M. Goryachev, J. Bourhill and M.E. Tobar, _The ORGAN Experiment: An axion haloscope above 15 GHz_ , _Phys. Dark Univ._ 18 (2017) 67 [1706.00209].
  * [37] B.T. McAllister and M.E. Tobar, _The ORGAN Experiment_ , _Springer Proc. Phys._ 245 (2020) 37.
  * [38] A.A. Melcón et al., _Axion Searches with Microwave Filters: the RADES project_ , _JCAP_ 05 (2018) 040 [1803.01243].
  * [39] CAST collaboration, _First results of the CAST-RADES haloscope search for axions at 34.67 $\mu$eV_, 2104.13798.
  * [40] D. Alesini et al., _Search for invisible axion dark matter of mass m ${}_{a}=43~{}\mu$eV with the QUAX–$a\gamma$ experiment_, _Phys. Rev. D_ 103 (2021) 102004 [2012.09498].
  * [41] MADMAX collaboration, _A new experimental approach to probe QCD axion dark matter in the mass range above 40 $\mu$eV_, _Eur. Phys. J. C_ 79 (2019) 186 [1901.07401].
  * [42] P. Sikivie, N. Sullivan and D.B. Tanner, _Proposal for Axion Dark Matter Detection Using an LC Circuit_ , _Phys. Rev. Lett._ 112 (2014) 131301 [1310.8545].
  * [43] Y. Kahn, B.R. Safdi and J. Thaler, _Broadband and Resonant Approaches to Axion Dark Matter Detection_ , _Phys. Rev. Lett._ 117 (2016) 141801 [1602.01086].
  * [44] J.L. Ouellet et al., _First Results from ABRACADABRA-10 cm: A Search for Sub- $\mu$eV Axion Dark Matter_, _Phys. Rev. Lett._ 122 (2019) 121802 [1810.12257].
  * [45] N. Crisosto, P. Sikivie, N.S. Sullivan, D.B. Tanner, J. Yang and G. Rybka, _ADMX SLIC: Results from a Superconducting $LC$ Circuit Investigating Cold Axions_, _Phys. Rev. Lett._ 124 (2020) 241101 [1911.05772].
  * [46] IAXO collaboration, _Physics potential of the International Axion Observatory (IAXO)_ , _JCAP_ 06 (2019) 047 [1904.09155].
  * [47] P. Sikivie, _Invisible Axion Search Methods_ , _Rev. Mod. Phys._ 93 (2021) 015004 [2003.02206].
  * [48] I.G. Irastorza and J. Redondo, _New experimental approaches in the search for axion-like particles_ , _Prog. Part. Nucl. Phys._ 102 (2018) 89 [1801.08127].
  * [49] J. Billard et al., _Direct Detection of Dark Matter – APPEC Committee Report_ , 2104.07634.
  * [50] P. Agrawal, G. Marques-Tavares and W. Xue, _Opening up the QCD axion window_ , _JHEP_ 03 (2018) 049 [1708.05008].
  * [51] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi and F. Takahashi, _Relic Abundance of Dark Photon Dark Matter_ , _Phys. Lett. B_ 801 (2020) 135136 [1810.07188].
  * [52] K. Choi, H.B. Kim and J.E. Kim, _Axion cosmology with a stronger qcd in the early universe_ , _Nuclear Physics B_ 490 (1997) 349–364.
  * [53] T. Banks and M. Dine, _The cosmology of string theoretic axions_ , _Nuclear Physics B_ 505 (1997) 445–460.
  * [54] T. Banks, M. Dine and M. Graesser, _Supersymmetry, axions, and cosmology_ , _Physical Review D_ 68 (2003) .
  * [55] L. Heurtier, F. Huang and T.M.P. Tait, _Resurrecting Low-Mass Axion Dark Matter Via a Dynamical QCD Scale_ , 2104.13390.
  * [56] P.J. Steinhardt and M.S. Turner, _Saving the Invisible Axion_ , _Phys. Lett. B_ 129 (1983) 51.
  * [57] G. Lazarides, R.K. Schaefer, D. Seckel and Q. Shafi, _Dilution of Cosmological Axions by Entropy Production_ , _Nucl. Phys. B_ 346 (1990) 193.
  * [58] M. Kawasaki, T. Moroi and T. Yanagida, _Can decaying particles raise the upper bound on the Peccei-Quinn scale?_ , _Phys. Lett. B_ 383 (1996) 313 [hep-ph/9510461].
  * [59] G.F. Giudice, E.W. Kolb and A. Riotto, _Largest temperature of the radiation era and its cosmological implications_ , _Phys. Rev. D_ 64 (2001) 023508 [hep-ph/0005123].
  * [60] L. Visinelli and P. Gondolo, _Axion cold dark matter in non-standard cosmologies_ , _Phys. Rev. D_ 81 (2010) 063508 [0912.0015].
  * [61] M. Venegas, _Relic Density of Axion Dark Matter in Standard and Non-Standard Cosmological Scenarios_ , 2106.07796.
  * [62] D. Grin, T.L. Smith and M. Kamionkowski, _Axion constraints in non-standard thermal histories_ , _Phys. Rev. D_ 77 (2008) 085020 [0711.1352].
  * [63] P. Carenza, M. Lattanzi, A. Mirizzi and F. Forastieri, _Thermal axions with multi-eV masses are possible in low-reheating scenarios_ , 2104.03982.
  * [64] N. Blinov, M.J. Dolan and P. Draper, _Imprints of the Early Universe on Axion Dark Matter Substructure_ , _Phys. Rev. D_ 101 (2020) 035002 [1911.07853].
  * [65] N. Ramberg and L. Visinelli, _Probing the Early Universe with Axion Physics and Gravitational Waves_ , _Phys. Rev. D_ 99 (2019) 123513 [1904.05707].
  * [66] A.E. Nelson and H. Xiao, _Axion Cosmology with Early Matter Domination_ , _Phys. Rev. D_ 98 (2018) 063516 [1807.07176].
  * [67] L. Visinelli and J. Redondo, _Axion Miniclusters in Modified Cosmological Histories_ , _Phys. Rev. D_ 101 (2020) 023008 [1808.01879].
  * [68] M. Drees, F. Hajkarim and E.R. Schmitz, _The Effects of QCD Equation of State on the Relic Density of WIMP Dark Matter_ , _JCAP_ 06 (2015) 025 [1503.03513].
  * [69] Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * [70] D. Harari and P. Sikivie, _The Dynamics of Global Strings_ , in _7 th Moriond Workshop: Searches for New and Exotic Phenomena_, pp. 513–517, 1987.
  * [71] F. D’Eramo, N. Fernandez and S. Profumo, _When the Universe Expands Too Fast: Relentless Dark Matter_ , _JCAP_ 05 (2017) 012 [1703.04793].
  * [72] R. Allahverdi, K. Dutta and A. Maharana, _Constraining Non-thermal Dark Matter by CMB_ , _JCAP_ 10 (2018) 038 [1808.02659].
  * [73] R. Allahverdi and J.K. Osiński, _Freeze-in Production of Dark Matter Prior to Early Matter Domination_ , _Phys. Rev. D_ 101 (2020) 063503 [1909.01457].
  * [74] P. Arias, D. Karamitros and L. Roszkowski, _Frozen-in fermionic singlet dark matter in non-standard cosmology with a decaying fluid_ , 2012.07202.
  * [75] A. Vilenkin and L.H. Ford, _Gravitational Effects upon Cosmological Phase Transitions_ , _Phys. Rev. D_ 26 (1982) 1231.
  * [76] G.D. Coughlan, W. Fischler, E.W. Kolb, S. Raby and G.G. Ross, _Cosmological Problems for the Polonyi Potential_ , _Phys. Lett. B_ 131 (1983) 59.
  * [77] T. Moroi and L. Randall, _Wino cold dark matter from anomaly mediated SUSY breaking_ , _Nucl. Phys. B_ 570 (2000) 455 [hep-ph/9906527].
  * [78] J.D. Barrow, _Massive Particles as a Probe of the Early Universe_ , _Nucl. Phys. B_ 208 (1982) 501.
  * [79] L.H. Ford, _Gravitational Particle Creation and Inflation_ , _Phys. Rev. D_ 35 (1987) 2955.
  * [80] B. Spokoiny, _Deflationary universe scenario_ , _Phys. Lett. B_ 315 (1993) 40 [gr-qc/9306008].
  * [81] R. Allahverdi et al., _The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe_ , _Open J. Astrophys._ 4 (2020)  [2006.16182].
  * [82] D.J.H. Chung, E.W. Kolb and A. Riotto, _Production of massive particles during reheating_ , _Phys. Rev. D_ 60 (1999) 063504 [hep-ph/9809453].
  * [83] R.R. Caldwell, R. Dave and P.J. Steinhardt, _Cosmological imprint of an energy component with general equation of state_ , _Phys. Rev. Lett._ 80 (1998) 1582 [astro-ph/9708069].
  * [84] V. Sahni and A.A. Starobinsky, _The Case for a positive cosmological Lambda term_ , _Int. J. Mod. Phys. D_ 9 (2000) 373 [astro-ph/9904398].
  * [85] J. Khoury, B.A. Ovrut, P.J. Steinhardt and N. Turok, _The Ekpyrotic universe: Colliding branes and the origin of the hot big bang_ , _Phys. Rev. D_ 64 (2001) 123522 [hep-th/0103239].
  * [86] MADMAX Working Group collaboration, _Dielectric Haloscopes: A New Way to Detect Axion Dark Matter_ , _Phys. Rev. Lett._ 118 (2017) 091801 [1611.05865].
  * [87] J. Egge, S. Knirck, B. Majorovits, C. Moore and O. Reimann, _A first proof of principle booster setup for the MADMAX dielectric haloscope_ , _Eur. Phys. J. C_ 80 (2020) 392 [2001.04363].
  * [88] S. Beurthey et al., _MADMAX Status Report_ , 2003.10894.
  * [89] E. Armengaud et al., _Conceptual Design of the International Axion Observatory (IAXO)_ , _JINST_ 9 (2014) T05002 [1401.3233].
  * [90] T. Dafni, C.A.J. O’Hare, B. Lakić, J. Galán, F.J. Iguaz, I.G. Irastorza et al., _Weighing the solar axion_ , _Phys. Rev. D_ 99 (2019) 035037 [1811.09290].
  * [91] J.E. Kim, _Weak Interaction Singlet and Strong CP Invariance_ , _Phys. Rev. Lett._ 43 (1979) 103.
  * [92] M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, _Can Confinement Ensure Natural CP Invariance of Strong Interactions?_ , _Nucl. Phys. B_ 166 (1980) 493.
  * [93] M. Dine, W. Fischler and M. Srednicki, _A Simple Solution to the Strong CP Problem with a Harmless Axion_ , _Phys. Lett. B_ 104 (1981) 199.
  * [94] C. O’Hare, _cajohare/axionlimits: Axionlimits_ , July, 2020.  10.5281/zenodo.3932430.
  * [95] M. Regis, M. Taoso, D. Vaz, J. Brinchmann, S.L. Zoutendijk, N.F. Bouché et al., _Searching for light in the darkness: Bounds on ALP dark matter with the optical MUSE-faint survey_ , _Phys. Lett. B_ 814 (2021) 136075 [2009.01310].
  * [96] D. Grin, G. Covone, J.-P. Kneib, M. Kamionkowski, A. Blain and E. Jullo, _A Telescope Search for Decaying Relic Axions_ , _Phys. Rev. D_ 75 (2007) 105018 [astro-ph/0611502].
  * [97] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, _WISPy Cold Dark Matter_ , _JCAP_ 06 (2012) 013 [1201.5902].
  * [98] G.G. Raffelt, _Astrophysical methods to constrain axions and other novel particle phenomena_ , _Phys. Rept._ 198 (1990) 1.
  * [99] P. Carenza, T. Fischer, M. Giannotti, G. Guo, G. Martínez-Pinedo and A. Mirizzi, _Improved axion emissivity from a supernova via nucleon-nucleon bremsstrahlung_ , _JCAP_ 10 (2019) 016 [1906.11844].
  * [100] P. Carenza, B. Fore, M. Giannotti, A. Mirizzi and S. Reddy, _Enhanced Supernova Axion Emission and its Implications_ , _Phys. Rev. Lett._ 126 (2021) 071102 [2010.02943].
  * [101] D.H. Lyth, _Axions and inflation: Sitting in the vacuum_ , _Phys. Rev. D_ 45 (1992) 3394.