DESY 21-011
IFT-UAM/CSIC-20-144
FTUAM-20-21

Dark matter from an even lighter QCD axion: trapped misalignment

Luca Di Luzioa, Belen Gavela${}^{b,\,c}$, Pablo Quileza, Andreas Ringwalda

aDeutsches Elektronen-Synchrotron DESY,
Notkestraße 85, D-22607 Hamburg, Germany
bDepartamento de Fisica Teorica, Universidad Autonoma de Madrid,
Cantoblanco, 28049, Madrid, Spain
cInstituto de Fisica Teorica, IFT-UAM/CSIC,
Cantoblanco, 28049, Madrid, Spain

###### Abstract

We show that dark matter can be accounted for by an axion that solves the
strong CP problem, but is much lighter than usual due to a $Z_{\cal N}$
symmetry. The whole mass range from the canonical QCD axion down to the ultra-
light regime is allowed, with $3\leq{\cal N}\lesssim 65$. This includes the
first proposal of a “fuzzy dark matter” QCD axion with $m_{a}\sim 10^{-22}$
eV. A novel misalignment mechanism occurs –trapped misalignment– due to the
peculiar temperature dependence of the $Z_{{\cal N}}$ axion potential. The
dark matter relic density is enhanced because the axion field undergoes two
stages of oscillations: it is first trapped in the wrong minimum, which
effectively delays the onset of true oscillations. Trapped misalignment is
more general than the setup discussed here, and may hold whenever an extra
source of Peccei-Quinn breaking appears at high temperatures. Furthermore, it
will be shown that trapped misalignment can dynamically source the recently
proposed kinetic misalignment mechanism. All the parameter space is within
tantalizing reach of the experimental projects for the next decades. For
instance, even Phase I of CASPEr-Electric could discover this axion.

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>,
<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Contents

  1. 1 Introduction
  2. 2 Cosmological constraints on mirror worlds
    1. 2.1 On asymmetric SM/mirror temperatures
  3. 3 $Z_{\cal N}$ axion dark matter
    1. 3.1 Temperature dependence of the $Z_{\cal N}$ axion potential
    2. 3.2 $Z_{\cal N}$ axion misalignment
      1. 3.2.1 One stage of oscillations: simple ALP regime
      2. 3.2.2 Two stages of oscillation: trapped misalignment
        1. 3.2.2.1 Pure trapped misalignment: $\dot{\theta}_{\rm tr}\sim 0$
        2. 3.2.2.2 Trapped+kinetic misalignment: sizeable $\dot{\theta}_{\rm tr}$
      3. 3.2.3 Solving both the strong CP problem and the nature of DM
    3. 3.3 Trapped vs. non-trapped mechanisms
  4. 4 Implications for axion dark matter searches
    1. 4.1 Axion coupling to photons
    2. 4.2 Axion coupling to nucleons
    3. 4.3 Axion coupling to electrons
    4. 4.4 Axion coupling to nEDM
    5. 4.5 Ultra-light QCD axion dark matter
  5. 5 Conclusions
  6. A Ratios of $g_{*}$ and $g_{s}$
  7. B Anharmonicity function
  8. C Maximal relic density in the trapped+kinetic regime

##  1 Introduction

Were an axion-like particle (ALP) to be ever discovered, it would be
compelling to explore whether it has something to do with the strong CP
problem [1, 2, 3, 4], and/or whether it can be a viable dark matter (DM)
candidate [5, 6, 7]. The canonical QCD axion (aka “invisible axion”) satisfies

$m_{a}^{\rm
QCD}f_{a}=m_{\pi}\,f_{\pi}\,\frac{\sqrt{z}}{(1+z)}\,,\qquad\text{where}\quad
z=\frac{m_{u}}{m_{d}}\,$ (1.1)

and $m_{a}^{\rm QCD}$, $f_{a}$, $f_{\pi},m_{\pi},m_{u}$ and $m_{d}$ denote
respectively the axion mass, the axion scale, the pion decay constant, and the
pion, up and down quark masses [8, 9, 10, 11, 12, 13]. Eq. (1.1) holds
whenever QCD is the only confining group of the theory, irrespectively of the
ultraviolet (UV) model details. A departure from this $m_{a}$-$f_{a}$ relation
always requires to extend the gauge confining sector beyond the Standard Model
(SM) QCD group.

Axions that solve the strong CP problem but are heavier than the canonical QCD
axion have been explored since long and revived in the last years [14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. In contrast,
solutions to the strong CP problem with lighter axions were uncharted
territory until very recently.

The goal of this work is to determine the implications for DM of a freshly
proposed dynamical –and technically natural– scenario [31, 32], which solves
the strong CP problem with an axion much lighter than the canonical QCD one.
Its key point was to assume that Nature is endowed with a $Z_{\cal N}$
symmetry realized non-linearly by the axion field [31, 32]. ${\cal N}$ mirror
and degenerate worlds linked by the axion field would coexist with the same
coupling strengths as in the SM, with the exception of the effective
$\theta_{k}$-parameters,

$\mathcal{L}=\sum_{k=0}^{{\cal
N}-1}\left[\mathcal{L}_{\text{SM}_{k}}+\frac{\alpha_{s}}{8\pi}\left(\theta_{a}+\frac{2\pi
k}{{\cal N}}\right)G_{k}\widetilde{G}_{k}\right]\,+\dots\,.$ (1.2)

Here, $\mathcal{L}_{\text{SM}_{k}}$ denotes exact copies of the SM total
Lagrangian excluding the topological $G_{k}\tilde{G}_{k}$ term,111The SM is
identified from now on with the $k=0$ sector and this label on SM quantities
will be often dropped in the following.  $\theta_{a}\in[-\pi,\pi)$ is defined
in terms of the axion field $a$, $\theta_{a}\equiv a/f_{a}$, and the dots
stand for $Z_{\cal N}$-symmetric portal couplings among the SM copies.

The solution to the strong CP problem of this $Z_{\cal N}$ scenario required
${\cal N}$ to be odd. Overall, the $\sim 10$ orders of magnitude tuning
required by the SM strong CP problem is traded by a $1/{\cal N}$ adjustment,
where ${\cal N}$ could be as low as ${\cal N}=3$ (viable solutions were found
with $3\leq{\cal N}\lesssim 47$ for $f_{a}\lesssim 2.4\times 10^{15}$ GeV, and
any ${\cal N}$ for larger values of $f_{a}$ [32]). The resulting axion is
exponentially lighter than the canonical QCD one, because the non-perturbative
contributions to its potential from the ${\cal N}$ degenerate QCD groups
conspire by symmetry to suppress each other. This can be intuitively
understood from the large ${\cal N}$ limit of a non-linearly realized $Z_{\cal
N}$ shift symmetry, which is a continuous global $U(1)$ symmetry: the axion
acts as the Goldstone boson of the discrete symmetry [33]. Its mass thus
vanishes asymptotically for large ${\cal N}$ as befits a $U(1)$ Goldstone
boson. Indeed, it has been shown [32] that in the large ${\cal N}$ limit the
total axion potential is given in all generality by

$V_{\mathcal{N}}\left(\theta_{a}\right)\simeq-\frac{m_{a}^{2}f_{a}^{2}}{{\cal
N}^{2}}\,\cos({\cal N}\theta_{a})\,,$ (1.3)

where now the axion mass in vacuum obeys

$m_{a}^{2}\simeq\frac{m_{\pi}^{2}f_{\pi}^{2}}{f_{a}^{2}}\frac{1}{\sqrt{\pi}}\sqrt{\frac{1-z}{1+z}}\,\,\mathcal{N}^{3/2}\,z^{\mathcal{N}}\,,$
(1.4)

which is exponentially suppressed ($\propto z^{\cal N}$) in comparison to
$(m_{a}^{\rm QCD})^{2}$.222Note that the value of $m_{a}$ for ${\cal N}=1$
should be equal to $m_{a}^{\rm QCD}$ given in Eq. (1.1), while Eq. (1.4) only
holds in the large ${\cal N}$ limit.

The crucial properties of such a light axion are generic and do not depend on
the details of the putative UV completion. An important byproduct of this
construction is an enhancement of all axion interactions which is universal,
that is, model-independent and equal for all axion couplings, at fixed
$m_{a}$. The detailed exploration of the $Z_{\cal N}$ paradigm and of the
phenomenological constraints which do not require the axion to account for DM
can be found in Ref. [32].

We will determine here instead the parameter space that solves both the strong
CP problem and the nature of DM within the $Z_{\cal N}$ framework under
discussion. The question of DM is of strong experimental interest and very
timely, in particular given the plethora of projects targeting light DM
candidates, down to the region of “fuzzy” DM [34]. For instance, CASPEr-
Electric [35, 36, 37] probes directly the anomalous gluonic coupling of the
axion via an oscillating neutron electric dipole moment (nEDM): the strength
of that coupling cannot be modified for the canonical QCD axion irrespective
of the model details, unlike other axion-to-SM couplings that can be
selectively enhanced [38, 39, 40, 41, 42, 43, 44]. We show here that, in
contrast, a hypothetical early discovery at CASPEr-Electric could be
interpreted as a solution to the strong CP problem via a $Z_{\cal N}$ reduced-
mass axion, because of the same-size enhancement of all axion couplings in
such a scenario.

It will be shown that the cosmology of the $Z_{\cal N}$ axion exhibits several
novel aspects. The cosmological impact of hypothetical parallel “mirror”
worlds has been studied at length in the literature (for a review, see e.g.
[45]). In particular, the constraints on the number of effective relativistic
species $N_{\text{eff}}$ present in the early Universe imply that the mirror
copies of the SM must be less populated –cooler– than the ordinary SM world
[46]. This requires in turn that the SM has never been in thermal contact with
its mirror replica. Fortunately, mechanisms that source this world-asymmetric
initial temperatures while preserving the $Z_{\cal N}$ symmetry may arise
naturally in the cosmological evolution [47, 48, 49]. It has also been
suggested that DM could be simply constituted by mirror matter, for which
relevant constraints apply given the differences in temperatures. Note,
nevertheless, that in most cases only one replica of the SM was considered,
while large ${\cal N}$ could significantly modify those analyses. Furthermore,
an axionic nature of DM has not been previously considered in the mirror world
setup with a reduced-mass axion.333An axion heavier than the QCD one has been
contemplated as DM, within a $Z_{2}$ setup which realized linearly the
symmetry [50]. While gravity and axion-mediated interactions are naturally
small enough, the impact of other possible ($Z_{\cal N}$-symmetric)
interactions on the thermal communication between the SM and its mirror copies
will be discussed.

The evolution of the $Z_{\cal N}$ axion field and its contribution to the DM
relic abundance will be shown to depart drastically from both the standard
case and the previously considered mirror world scenarios. Due to the peculiar
temperature dependence of the $Z_{{\cal N}}$ axion potential, the production
of DM axions in terms of the misalignment mechanism is modified. The scenario
results in a novel type of misalignment, with a large value of the
misalignment angle. In particular, we will show that the relic density is
enhanced because the axion field undergoes two stages of oscillations, that
are separated by an unavoidable and drastic –non-adiabatic– modification of
the potential. The axion field is first _trapped_ in the wrong minimum (with
$\theta=\pi$), which effectively delays the onset of the true oscillations and
thus enhances the DM density.444While this manuscript was being written, Ref.
[51] appeared where a trapping mechanism was used in a different context, i.e.
to reduce the axion DM abundance. Conversely, an alternative enhancement
mechanism has also been proposed recently [52]. We will call this new
production mechanism _trapped misalingment_. Note that inflation played a
crucial role in previous dynamical setups advocated to drive the initial
misalignment to $\pi$ [53, 54, 55], while our mechanism results directly from
temperature effects.

Furthermore it will be shown that, in some regions of the parameter space,
trapped misalignment will automatically source the recently proposed kinetic
misalignment mechanism [56]. In the latter, a sizeable initial axion velocity
is the source of the axion relic abundance as opposed to the conventionally
assumed initial misalignment angle. The original kinetic misalignment proposal
[56] required an ad-hoc Peccei Quinn (PQ)-breaking non-renormalizable
effective operator suppressed by powers of the Planck mass, $M_{\rm Pl}$. In
contrast, we will show that the early stage of oscillations in the $Z_{\cal
N}$ axion framework naturally flows out into kinetic misalignment.

The interplay of the different mechanisms is studied in detail. Moreover,
although general phenomenological consequences of a large misalignment angle
were studied in Ref. [57], we will identify the main novel consequences that
follow from the scenario under discussion. On the phenomenological arena, we
discuss the implications of the $Z_{\cal N}$ reduced-mass axion for axion DM
searches, namely those experiments that rely on the hypothesis that an axion
or ALP sizeably contributes to the DM density. More specifically, the
experimental prospects to probe its coupling to photons, nucleons, electrons
and the nEDM operator are considered. The analysis sweeps through the whole
mass range down to ultra-light axions (with masses $m_{a}\ll 10^{-10}$ eV),
within the $Z_{\cal N}$ axion framework under discussion. The present and
projected sensitivity to the number of possible mirror world ${\cal N}$ will
be determined, and the constraints obtained will be compared and combined with
those stemming from experiments which are independent of whether axions
account or not for DM [32].

The structure of the paper can be easily inferred from the Table of Contents.

##  2 Cosmological constraints on mirror worlds

The existence of a hypothetical parallel mirror world, with microphysics
identical to that of the observable world and connected to the latter only
gravitationally, has a non-trivial impact on cosmology (for a review, see e.g.
[45]). In particular, extra relativistic species in the early Universe (mirror
photons, electrons and neutrinos) affect the number of effective neutrino
species $N_{\rm eff}$ that can be measured through: i) the abundances of light
elements –in particular Helium– at the time of Big Bang nucleosynthesis (BBN),
i.e. at $T\sim 1\,{\rm MeV}$; ii) through the damping tail in the cosmic
microwave background (CMB) power spectrum at recombination, $T\sim 0.26\,{\rm
eV}$. Present data yield [58, 59, 60]555For BBN, the constraint corresponds to
$BBN+Yp+D$ in Table V of Ref. [58], translating their 68.3% CL result for the
number of neutrino species $N_{\nu}$ to a 95% CL value $N_{\nu}=2.85\pm 0.28$,
and $N_{\mathrm{eff}}=1.015\,N_{\nu}$. For CMB, the combination of
TT+TE+EE+lowE+lensing+BAO data by Planck 2018 [59, 60] is considered here;
note that due to the disagreement between local and CMB measurements of the
Hubble expansion parameter $H_{0}$, the constraint on $\Delta N_{\rm eff}$
would weaken to $N_{\mathrm{eff}}=3.27\pm 0.30$ if $H_{0}$ measurements were
included.

$\text{BBN: }N_{\text{eff}}=2.89\pm 0.57\,\,,\qquad\text{CMB: }N_{\text{eff
}}=2.99_{-0.33}^{+0.34}\,.$ (2.1)

In order to satisfy these constraints in the presence of one mirror world, its
temperature $T^{\prime}$ needs to be smaller than that of the SM, $T_{0}$. In
the $Z_{\cal N}$ scenario under discussion the different $k>0$ copies of the
SM could have different temperatures $T_{k}<T_{0}$, and hence different energy
densities $\rho_{k}$ and entropies $s_{k}$,

$\displaystyle\rho_{k}=\frac{\pi^{2}}{30}{g}_{*}(T_{k})T_{k}^{4}\,,\qquad
s_{k}=\frac{2\pi^{2}}{45}{g}_{s}(T_{k})T_{k}^{3}\,,$ (2.2)

where $g_{*}$ and $g_{s}$ denote the effective degrees of freedom related to
energy and entropy density, respectively. In principle, the functional forms
$g_{*},g_{s}$ could have an extra $k$-dependence through the respective
$\theta_{k}$ parameters, but the overall impact is expected to be minor and
will be disregarded in the following. Moreover, we assume here that thermally
produced axions give a negligible contribution to $N_{\rm eff}$. This is
typically the case for $f_{a}\gtrsim 10^{7}$ GeV [61, 62, 63].

The different sectors evolve in time with separately conserved entropies,
therefore the ratio of entropy densities $\gamma^{3}_{k}\equiv(s_{k}/s)$ is
constant, while the ratio of temperatures is given by

$\displaystyle\frac{T_{k}}{T}=\gamma_{k}\cdot\bigg{[}\frac{g_{s}(T)}{g_{s}\left(T_{k}\right)}\bigg{]}^{1/3}\equiv\gamma_{k}\cdot
b_{k}(T,T_{k}),\qquad\text{for }\,k\equiv 1,\dots,{\cal N}-1\,,$ (2.3)

where from now on $T$ will denote the temperature of the SM copy, $T\equiv
T_{0}$. The Hubble expansion rate depends on the total energy density of the
Universe. In a radiation dominated era with ${\cal N}$ mirror worlds it reads

$H^{2}(t)=\frac{4\pi^{3}}{45}\sum_{k=0}^{{\cal
N}-1}{{g}_{*}(T_{k})}\frac{T_{k}^{4}}{M^{2}_{\rm
Pl}}\equiv\frac{4\pi^{3}}{45}\frac{T^{4}}{M^{2}_{\rm
Pl}}\,{\overline{g}_{*}(T)}\,,$ (2.4)

with the number of effective degrees of freedom $\overline{g}_{*}(T)$ given by

$\overline{g}_{*}(T)=g_{*}(T)\left(1+\sum_{k=1}^{{\cal
N}-1}c_{k}\gamma_{k}^{4}\right)\,,$ (2.5)

and
$c_{k}\left(T,T_{k}\right)\equiv\left[g_{*}\left(T_{k}\right)/g_{*}(T)\right]\cdot\left[g_{s}(T)/g_{s}\left(T_{k}\right)\right]^{4/3}$.

Let us discuss next the implications for $N_{\rm eff}$. In the SM,
$g_{*}(T=1\,{\rm MeV})=5.5+\frac{7}{4}N^{\rm SM}_{\rm eff}=10.83$ where
$N^{\rm SM}_{\rm eff}=3.046$ is the effective number of neutrino species.
Analogously, at recombination $g_{*}(T=0.26\,{\rm eV})=3.93$. It follows from
Eq. 2.5 that the contribution of the ensemble of ${\cal N}$ worlds to $\Delta
N_{\rm eff}=N_{\rm eff}-N^{\rm SM}_{\rm eff}$ at BBN and CMB is given by

$\displaystyle\text{BBN: }\quad\Delta N_{\rm eff}=\frac{4}{7}g_{*}(T=1\,{\rm
MeV})\sum_{k=1}^{{\cal
N}-1}c_{k}\gamma_{k}^{4}\,\quad\Longrightarrow\quad\sum_{k=1}^{{\cal
N}-1}c_{k}\gamma_{k}^{4}<0.067\,,$ (2.6) $\displaystyle\text{CMB: }\quad\Delta
N_{\rm eff}=\frac{4}{7}g_{*}(T=0.26\,{\rm eV})\sum_{k=1}^{{\cal
N}-1}c_{k}\gamma_{k}^{4}\,\quad\Longrightarrow\quad\sum_{k=1}^{{\cal
N}-1}c_{k}\gamma_{k}^{4}<0.13\,,$ (2.7)

where the bounds on the right-hand side stem from the constraints in Eq.
(2.1). These expressions illustrate an interesting prediction of the mirror
world(s) scenario: the deviation in the number of measured effective neutrino
species with respect $N^{\rm SM}_{\rm eff}$ is not constant but varies with
temperature. That is, it may vary with time as certain mirror species will
become non-relativistic at different times in the evolution, because of the
different temperatures of the mirror worlds.

####  $c_{k}\simeq 1$ approximation

The temperature dependence of the function ${g}_{s}$ is quite flat in most
regions of interest, ${g}_{s}(T_{k\neq 0})\simeq{g}_{s}(T)$. Within this
approximation $b_{k\neq 0}\simeq c_{k\neq 0}\simeq 1$, and the parameter
$\gamma_{k}$ directly gives the ratios of temperatures

$\gamma_{k}\simeq\frac{T_{k}}{T}\,,$ (2.8)

which remain constant throughout the Universe evolution, with $\gamma_{k}<1$
for all $k\neq 0$. This approximation is valid as far as the temperatures
involved correspond to the same plateau of the $N_{\rm eff}$ distribution (see
e.g. [64]). In other words, as long as the number of relativistic degrees of
freedom does not change between $T_{k}$ and $T$ and thus $c_{k}\simeq 1$
holds. This is the case at the CMB temperature, where Eq. (2.8) holds for all
$k$ because all exotic worlds are cooler than the SM one. For temperatures
around the BBN region, we have checked that the approximation is still good as
long as $\gamma_{k}$ is not too small, e.g. better than $25\%$ for
$\gamma_{k\neq 0}\gtrsim 0.1$, which does not change noticeably the analytical
expressions. We will work in this framework all through the rest of the paper,
except for the numerical results. More details can be found in Appendix A.

In general, the ${\cal N}-1$ copies of the SM may have different temperatures.
In this case Eqs. (2.6)-(2.7) allow to constrain the temperature of the
hottest of them, ${T_{\rm max}}$,

$\text{BBN:}\quad\frac{T_{\rm
max}}{T}<0.51\,,\qquad\text{CMB:}\quad\frac{T_{\rm max}}{T}<0.60\,.$ (2.9)

If all mirror images of the SM are instead assumed to have the same
temperature, $T_{k}\sim T^{\prime}$ for $k\neq 0$, the most restrictive bounds
follow,

$\text{BBN:}\quad\frac{T^{\prime}}{T}<\frac{0.51}{({\cal
N}-1)^{1/4}}\,,\qquad\text{CMB:}\quad\frac{T^{\prime}}{T}<\frac{0.60}{({\cal
N}-1)^{1/4}}\,.$ (2.10)

It is interesting that BBN data set the most constraining bound on $N_{\rm
eff}$, in spite of CMB measurements being more precise. This is due to the
temperature dependence of the mirror worlds contribution. Moreover, this non-
trivial dependence represents a smoking gun for the existence of the latter,
as it generically predicts incompatible measurements of $N_{\rm eff}$ from BBN
and CMB. Specifically, the scenario predicts in all generality the following
discrepancy:

$\displaystyle N^{\rm BBN}_{\rm eff}-N^{\rm CMB}_{\rm
eff}=3.92\sum_{k=0}^{{\cal
N}-1}c_{k}\gamma_{k}^{4}\,\simeq\,3.92\sum_{k=0}^{{\cal
N}-1}\gamma_{k}^{4}\,.$ (2.11)

If such a difference were to be experimentally established, it would allow to
predict the temperature of the mirror worlds e.g. in the two limiting cases in
Eqs. (2.9)–(2.10).

#### Constraints on portal couplings

The SM must avoid thermal contact with its mirror copies all through the
(post-inflation) history of the Universe so as to fullfil the condition
$T_{k\neq 0}\ll T$. This implies that the interactions between the SM and its
copies need to be very suppressed. Non-renormalizable interactions such as
gravity and axion-mediated ones are naturally small enough, while the Higgs
and hypercharge kinetic portal couplings can potentially spoil the condition
$T_{k\neq 0}\ll T$. For instance, in the $Z_{2}$ mirror case
$T^{\prime}/T\lesssim 0.5$ requires both portal couplings, defined as
$\mathcal{L}\supset\kappa|H|^{2}|H^{\prime}|^{2}+\epsilon
B^{\mu\nu}B^{\prime}_{\mu\nu}$, where $H$ and $H^{\prime}$ ($B^{\mu\nu}$ and
$B^{\prime}_{\mu\nu}$) denote respectively the SM Higgs doublet (hypercharge
field strength) and its mirror copy, and $\kappa$ and $\epsilon$ are
dimensionless couplings, to respect $\kappa\,,\epsilon\lesssim 10^{-8}$ [65].
Even smaller couplings are needed in the $Z_{\cal N}$ case with ${\cal N}>2$.
This can suggest a ‘naturalness’ issue for the Higgs and kinetic portal
couplings, as they cannot be forbidden in terms of internal symmetries.
Nevertheless, such small couplings may be technically natural because of an
enhanced _Poincare_ symmetry [66, 67]: in the limit where non-renormalizable
interactions are neglected, the $\kappa,\epsilon\to 0$ limit for the ensemble
of portals corresponds to an enhanced $\mathcal{P}^{\cal N}$ symmetry (namely,
an independent space-time Poincaré transformation $\mathcal{P}$ in each
sector). This protects those couplings from radiative corrections other than
those induced by the explicit $\mathcal{P}^{\cal N}$ breaking due to
gravitational and axion-mediated interactions, which are presumably small.

###  2.1 On asymmetric SM/mirror temperatures

The microphysics responsible for the evolution of the early SM Universe and of
its mirror copies is almost the same. Which mechanisms can then source
different temperatures for the SM and its replicae?

One difference in the microphysics of our setup is the axion coupling to the
$G_{\mu\nu}\tilde{G}^{\mu\nu}$ pseudo-scalar density in Eq. (1.2): the
effective value of the $\theta$ parameter differs for each sector $k$,
$\theta_{k}=2\pi k/{\cal N}$ (and thus relaxing to zero in the SM with
probability $1/{\cal N}$ – see Sect. 3). This implies that nuclear physics
would be drastically different for the SM and its mirror copies. Indeed, the
one-pion scalar exchange parametrized by the effective baryon chiral
Lagrangian,

$\mathcal{L}_{\chi\text{PT}}\supset c_{+}\sum_{k=0}^{{\cal
N}-1}\frac{m_{u}m_{d}\sin\theta_{k}}{[m_{u}^{2}+m_{d}^{2}+2m_{u}m_{d}\cos\theta_{k}]^{1/2}}\frac{\pi_{k}^{a}}{f_{\pi}}\overline{N}_{k}\tau^{a}N_{k}\,,$
(2.12)

where $\pi_{k}$ and $N_{k}$ denote respectively the pion and nucleon fields
and $c_{+}$ is an $\mathcal{O}(1)$ low-energy constant, induces long-range
forces (i.e. not spin-suppressed in the non-relativistic limit) among nucleons
in all worlds but the SM one [68]. This may lead to different thermalization
histories of the mirror copies, modulated by their ‘distance’ $k$ from the SM
sector ($k=0$). Qualitatively, nucleons are kept longer in thermal equilibrium
in the mirror replicas of the SM due to this extra interaction channel, which
could play a role in generating different temperatures. However, a
quantitative estimate is complicated by the fact that the interaction in Eq.
(2.12) becomes relevant only around the QCD phase transition.

Previously in the literature, mirror world scenarios leading to sufficiently
different temperatures relied on specific inflation implementations. Some of
them could also hold in our $Z_{\cal N}$ scenario. For instance, in the so-
called _asymmetric inflation_ it is argued that the initial condition after
inflation may be very different for the SM and the various mirror copies, even
if they have very similar microphysics dynamics. This mechanism was proposed
in the context of an exactly $Z_{2}$ symmetric mirror world [47, 48]. Two
inflaton fields would be present: $\phi$ for the SM and $\phi^{\prime}$ for
its $Z_{2}$ mirror copy, respectively slow-rolling towards the minimum of
their potentials (the generalization to the case of $Z_{\cal N}$ is
straightforward). Inflation ends when the potential steepens and the inflaton
field picks up kinetic energy so that it starts to oscillate towards the
minimum of its potential. In this process, quantum fluctuations can be such
that the two inflaton fields do not necessarily reach their minimum at the
same time, even in the same spatial region. In regions where $\phi^{\prime}$
reaches its minimum first, any mirror particles produced due to its
oscillations are diluted by the remaining inflation driven by the second field
$\phi$. Hence, by the time when $\phi$ reheats the SM sector, the density of
mirror matter will be very small, and in consequence $T^{\prime}\ll T$.
Quantitatively, whether such asymmetric inflation scenario can be realized or
not depends on the shape of the inflaton potential [48]. An alternative idea
to achieve different temperatures is that the inflaton itself transforms non-
linearly under the $Z_{\cal N}$ symmetry [49]: this automatically implies
different values for the inflaton fields in each mirror world, and thus
different thermal histories. In the present paper we will not commit ourselves
to a specific mechanism realizing different SM/mirror temperatures, but simply
assume as a working hypothesis $T_{k\neq 0}<T$.

Mirror baryonic matter has been advocated to explain partially or even totally
the nature of DM, in the context of the $Z_{2}$-symmetric mirror world (see
e.g. [45]). This may be achieved, in spite of the lower temperature of the
mirror world, by means of a proper baryogenesis mechanism in the mirror sector
[46]. In the following, we will not dwell with mirror baryonic DM (which is
assumed to give a negligible contribution due to the lower temperature of the
SM mirror copies) and focus instead on axionic DM.

##  3 $Z_{\cal N}$ axion dark matter

One of the key ingredients in order to compute the axion relic density is the
temperature at the onset of oscillations around the true minimum, since it
determines for how long axion oscillations are damped (and therefore to what
extent the relic density is diluted). In the traditional scenario, that
corresponds to the temperature at which the Hubble parameter is of the order
of the axion mass. We analyze here the intrinsic differences between the
$Z_{\cal N}$ scenario with a reduced-mass axion and the canonical QCD axion.

An important question is whether the $Z_{\cal N}$ axion can account for the
observed DM density, without requiring any further adjustment other than the
inherent $1/{\cal N}$ tuning needed to solve the strong CP problem, i.e. to
choose $\theta_{a}=0$ as the vacuum for the SM world among the ${\cal N}$
possible $\theta_{a}$ values, $\theta_{a}=\\{\pm 2\pi\ell/\mathcal{N}\\}\text{
with }\ell=0,1,\ldots,(\mathcal{N}-1)/2$. It will be demonstrated in this
section that this is indeed the case.

We will focus here on the study of the axionic zero mode, assuming spatial
homogeneity and the pre-inflationary PQ breaking scenario.

###  3.1 Temperature dependence of the $Z_{\cal N}$ axion potential

At temperatures well below the critical temperature of the QCD phase
transition, $T\ll T_{\rm QCD}\simeq 150$ MeV, the customary QCD axion
potential is well approximated by the zero temperature chiral potential.
However, the chiral expansion breaks down close to $T_{\rm QCD}$. Lattice
computations of the topological susceptibility $\chi(T)$ are instead available
for the cross-over regime, from which the QCD axion mass is estimated. As the
exact temperature dependence is not known, we will parametrize the QCD axion
potential at finite temperature $V(\theta_{a},T)$ by weighing the zero
temperature expression of the chiral axion potential $V(\theta_{a})$ by a
$h(T)$ factor

$V(\theta_{a},T)\equiv V(\theta_{a})\,\times
h(T)=-m_{\pi}^{2}f_{\pi}^{2}\sqrt{1-\frac{4m_{u}m_{d}}{\left(m_{u}+m_{d}\right)^{2}}\sin^{2}\left(\frac{\theta_{a}}{2}\right)}\,\times
h(T)\,,$ (3.1)

where

$h(T)\equiv\left\\{\begin{array}[]{ll}\,1&\text{ for }T<T_{\rm QCD}\\\
\left(\frac{T_{\rm QCD}}{T}\right)^{\alpha}&\text{ for }T>T_{\rm
QCD}\end{array}\right.\,.$ (3.2)

Constant terms in the potential have been left implicit in this expression, as
irrelevant for the axion mass and cosmological evolution. The coefficient
$\alpha$ parametrizes the power-law behaviour at high temperatures. Unless
otherwise stated, $\alpha=8$ will be assumed, as suggested by the dilute
instanton gas approximation (DIGA) and in agreement with some lattice
computations [69]. Note, however, that the potential in Eq. (3.1) differs from
the DIGA-motivated cosine potential commonly assumed for the QCD case: our
parametrization ensures that $V(\theta_{a},T)$ is a continuous function of the
temperature not only around the minimum but for any field value $\theta_{a}$.

The generalization of the ansatz above to the total $Z_{\cal N}$ axion
potential is simply given by the sum over the ${\cal N}$ different
potentials,666This potential depends on all the temperatures of the different
mirror worlds, but the simplified notation $V_{\cal N}(\theta_{a},T,T_{1}\dots
T_{{\cal N}-1})\equiv V_{\cal N}(\theta_{a},T)$ is used here for simplicity.
This notation also reflects the fact that the SM temperature $T$ –being the
largest– is the dominant driver of the Universe evolution.

$V_{\cal N}(\theta_{a},T)=-\sum_{k=0}^{{\cal
N}-1}m_{\pi}^{2}f_{\pi}^{2}\sqrt{1-\frac{4m_{u}m_{d}}{\left(m_{u}+m_{d}\right)^{2}}\sin^{2}\left(\frac{\theta_{a}}{2}+\frac{\pi
k}{\mathcal{N}}\right)}\times h(T_{k})\,,$ (3.3)

where the power-law temperature suppression is assumed to be $k$-independent
and as given in Eq. (3.2), that is, $\alpha_{k}=\alpha$ for all $k$.
Furthermore, $h(T_{k})=h(b_{k}\gamma_{k}T)\simeq h(\gamma_{k}T)$, see Eq. 2.3
and Appendix A. Irrespective of the specific $\gamma_{k\neq 0}$ values, three
temperature regimes can be distinguished which are qualitatively different.
They are determined by the minimum and maximum value of $\gamma_{k\neq 0}$,

$\gamma_{\text{min}}\,\leq\gamma_{k}\,\leq\gamma_{\text{max}}\,,$ (3.4)

corresponding respectively to the coolest ($\gamma_{\text{min}}T$) and the
hottest ($\gamma_{\text{max}}T$) mirror copies of the SM, see Fig 1.

Figure 1: Range of temperatures corresponding to qualitatively different
behaviours of the $Z_{\cal N}$ axion potential. The acronyms $HT,MT$ and $LT$
stand for high, medium and low temperatures, respectively.

#### Low temperatures $T<T_{\rm QCD}$

The potential for all ${\cal N}$ replicas is well approximated by the
–exponentially suppressed– zero-temperature potential in Eqs. (1.3) and (1.4),

$\displaystyle
V^{LT}_{\mathcal{N}}\left(\theta_{a},T\right)\,=\,V_{\mathcal{N}}\left(\theta_{a}\right)\,,$
(3.5)

with minima at $2\pi k/{\cal N}$ for ${\cal N}$ odd.

#### Medium temperatures $T_{\rm QCD}<T<T_{\rm QCD}/\gamma_{\text{max}}$

In this regime, only the SM contribution to the total potential is temperature
suppressed. It is weighed down by a factor $\left(T_{\rm
QCD}/T\right)^{\alpha}$, while all other contributions are still well
approximated by their zero-temperature potential,

$\displaystyle V^{MT}_{\mathcal{N}}\left(\theta_{a},T\right)$
$\displaystyle\simeq\left(\frac{T_{\rm
QCD}}{T}\right)^{\alpha}V(\theta_{a})+\sum_{k=1}^{{\cal N}-1}V(\theta_{a}+2\pi
k/{\cal N})$ $\displaystyle=\left[\left(\frac{T_{\rm
QCD}}{T}\right)^{\alpha}-1\right]V(\theta_{a})+\sum_{k=0}^{{\cal
N}-1}V(\theta_{a}+2\pi k/{\cal N})$
$\displaystyle\simeq\left[\left(\frac{T_{\rm
QCD}}{T}\right)^{\alpha}-1\right]V(\theta_{a})-m_{\pi}^{2}f_{\pi}^{2}\frac{1}{\sqrt{\pi}}\sqrt{\frac{1-z}{1+z}}\mathcal{N}^{-1/2}(-1)^{\mathcal{N}}z^{\mathcal{N}}\cos\left(\mathcal{N}\theta_{a}\right)$
$\displaystyle\xrightarrow{\,\,T\gg T_{\rm QCD}\,\,}-V(\theta_{a})\,,$ (3.6)

where in the last step the term $\propto z^{\cal N}$ has been neglected
compared to the size of $V(\theta_{a})$. It follows that the total axion
potential is well approximated in this regime by minus the QCD zero-
temperature axion potential, that is, the single $k=0$ world contribution with
the opposite sign. The potential thus depends only on $f_{a}$ and its minimum
is located at $\theta_{a}=\pm\,\pi$.

#### High temperatures $T>T_{\rm QCD}/\gamma_{\text{min}}$

The temperature effects are now important for all ${\cal N}$ replicas. In the
simplified case with all ${\cal N}-1$ mirror worlds having the same
temperature ($T_{k>0}\equiv T^{\prime}$, $\gamma_{k>0}\equiv\gamma^{\prime}$,
$b_{k>0}\equiv b^{\prime}$)

$\displaystyle V_{\mathcal{N}}^{HT}\left(\theta_{a},T\right)$
$\displaystyle\simeq\left(\frac{T_{\rm
QCD}}{T}\right)^{\alpha}\,V(\theta_{a})+\left(\frac{T_{\rm
QCD}}{T^{\prime}}\right)^{\alpha}\,\,\sum_{k=1}^{{\cal N}-1}V(\theta_{a}+2\pi
k/{\cal N})$ $\displaystyle\xrightarrow{\,\,{\cal N}\gg
1\,\,}-\left(\frac{T_{\rm
QCD}}{T\,\gamma^{\prime}\,b^{\prime}}\right)^{\alpha}V\left(\theta_{a}\right)\left[1-(\gamma^{\prime}\,b^{\prime})^{\alpha}\right]\xrightarrow{\,\,\gamma^{\prime}\ll
1\,\,}-\left(\frac{T_{\rm
QCD}}{T\,\gamma^{\prime}\,\,b^{\prime}}\right)^{\alpha}V\left(\theta_{a}\right)\,,$
(3.7)

where in the second step we neglected again the exponentially suppressed
$Z_{\cal N}$ contribution at $T=0$ and the definition $T^{\prime}\equiv
T\gamma^{\prime}b^{\prime}$ was used. As for the medium temperature regime,
the total axion potential depends essentially only on $f_{a}$ and its minimum
is located at $\theta_{a}=\pm\,\pi$. In the more general case in which the
mirror copies of the SM have different temperatures, only the coldest world
will give a sizeable contribution to the potential, and the minimum will no
longer be at $\pm\,\pi$ although typically not too far from it, see Fig. 2.

Figure 2: Temperature dependence of the $Z_{3}$ axion potential. Left panel:
all mirror worlds have the same temperature, $\gamma_{k}=\gamma=0.5$. Right
panel: mirror worlds have different temperatures,
$\gamma_{1}=0.5,\,\gamma_{2}=0.25$. The low and medium temperature potentials
do not depend on $\gamma_{k}$, as it can be appreciated by comparing the two
plots.

In summary, temperature effects lead to regimes in which the minimum of the
total axion potential is displaced to $\theta_{a}=\pm\,\pi$ or similar large
field values. This high-temperature shift of the minimum parallels that at
finite density, previously identified in Refs. [70, 71, 32] and recently
employed as well in Ref. [72]. This phenomenon is rich in consequences. In
particular, it follows that the axion mass at finite temperature,

$\displaystyle m_{a}(T)=\frac{1}{f_{a}}\sqrt{\frac{d^{2}V_{\cal
N}(\theta_{a},T)}{d\theta_{a}^{2}}\Bigg{|}_{\rm min}}\,\,,$ (3.8)

is given –in the three temperature regimes discussed above– by

$\displaystyle
m_{a}(T)\simeq\frac{m_{\pi}f_{\pi}}{f_{a}}\left\\{\begin{array}[]{lll}\displaystyle\frac{1}{\sqrt[4]{\pi}}\sqrt[4]{\frac{1-z}{1+z}}\,\,\mathcal{N}^{3/4}\,z^{\mathcal{N}/2}&\qquad\text{
for }\ \ \ T<T_{\rm QCD}\qquad\qquad\quad\quad\quad\quad\text{LT}\\\\[12.0pt]
\displaystyle\frac{\sqrt{z}}{1-z^{2}}&\qquad\text{ for }\quad T_{\rm
QCD}<T<T_{\rm QCD}/\gamma_{\text{max}}\qquad\\!\text{MT}\\\\[12.0pt]
\displaystyle\frac{\sqrt{z}}{1-z^{2}}\left(\frac{T_{\rm
QCD}}{\gamma_{\text{min}}\,T\,b_{\rm min}}\right)^{\alpha}&\qquad\text{ for
}\quad T_{\rm QCD}/\gamma_{\text{min}}<T\qquad\qquad\quad\ \
\text{HT}\,.\end{array}\right.$ (3.12)

The LT expression for $m_{a}(T)$ equals that for the zero-temperature mass
$m_{a}$ in Eq. (1.4), as expected. In contrast, for the intermediate
temperature case it corresponds to the tachyonic mass of the QCD axion at the
top of its potential,

$\displaystyle m_{a,\pi}^{\rm
QCD}\equiv\frac{1}{f_{a}}\sqrt{\frac{d^{2}V(\theta_{a})}{d\theta_{a}^{2}}}\Bigg{|}_{\theta_{a}=\pi}\,=\,\sqrt{\frac{1+z}{1-z}}\,\,m_{a}^{\rm
QCD}\,,$ (3.13)

which differs from $m_{a}^{\rm QCD}$ by a factor $\sim 1.7$ because the QCD
chiral axion potential is not a cosine, and its curvature around maxima and
minima are slightly different. Thus Eq. (3.12) can be rewritten as

$\displaystyle m_{a}(T)\simeq\left\\{\begin{array}[]{lll}\displaystyle
m_{a}&\qquad\text{ for }\ \ \ T<T_{\rm
QCD}\qquad\qquad\quad\quad\quad\quad\text{LT}\\\\[12.0pt] \displaystyle
m_{a,\pi}^{\rm QCD}&\qquad\text{ for }\quad T_{\rm QCD}<T<T_{\rm
QCD}/\gamma_{\text{max}}\qquad\\!\text{MT}\\\\[12.0pt] \displaystyle
m_{a,\pi}^{\rm QCD}\,\left(\frac{T_{\rm QCD}}{\gamma_{\text{min}}\,T\,b_{\rm
min}}\right)^{\alpha}&\qquad\text{ for }\quad T_{\rm
QCD}/\gamma_{\text{min}}<T\qquad\qquad\quad\ \ \text{HT}\,.\end{array}\right.$
(3.17)

In summary, at low temperatures the $Z_{\cal N}$ axion mass depends both on
$f_{a}$ and ${\cal N}$ and is exponentially suppressed by a factor $z^{{\cal
N}/2}$. In contrast, at medium and high temperatures it only depends on
$f_{a}$ in the large ${\cal N}$ limit, and is of the order of the canonical
QCD axion mass in vacuum and at high temperature respectively.

###  3.2 $Z_{\cal N}$ axion misalignment

We analyze here the consequences of the temperature-dependent $Z_{\cal N}$
axion potential for axion DM production. The axion relic density can be
estimated from the solution to the equation of motion (EOM) for a classical,
non-relativistic and homogeneous scalar field $\theta_{a}$ in an expanding
Universe:

$\displaystyle\ddot{\theta}_{a}+3H\dot{\theta}_{a}+\frac{1}{f_{a}^{2}}V^{\prime}(\theta_{a})=0\,,$
(3.18)

where $V^{\prime}(\theta_{a})={dV(\theta_{a})}/{d\theta_{a}}$ and the spatial
gradients have been neglected. For a field oscillating near its minimum,
$V^{\prime}(\theta_{a})\simeq m_{a}^{2}(T)f_{a}^{2}\,\theta_{a}$ is a good
approximation and the potential matches that of a damped harmonic oscillator.
The behaviour of the solution depends then on the interplay between two
variables, the axion mass $m_{a}(T)$ and the Hubble parameter $H(T)$, which
both vary with temperature.

We discuss next the evolution of the axion mass as a function of $T$, $f_{a}$
and ${\cal N}$. As the Universe cools down, the axion mass progressively grows
from its initial vanishingly small value at high temperatures, until it
suddenly drops down around $T_{\rm QCD}$ (feeling the $z^{\cal N}$ exponential
suppression of the potential) and reaches then a constant value. This is shown
in Fig. 3, where the $Z_{\cal N}$ trajectories are depicted in blue for
different values of ${\cal N}$. The figure illustrates the ${\cal
N}$-independence of $m_{a}(T)$ at high and medium temperatures (in the large
${\cal N}$ limit), while it becomes ${\cal N}$-dependent at low-temperature,
see also Eq. (3.12) and Fig. 2. At some point in that trajectory, the mass
overcomes the Hubble friction –depicted in  green– and the axion starts to
oscillate.

(a) $f_{a}=4\times 10^{18}\,\text{GeV}>f_{a}^{\rm crit}$: Simple ALP

(b) $f_{a}=1\times 10^{14}\,\text{GeV}<f_{a}^{\rm crit}$: Trapped mis.

Figure 3: Axion mass vs. Hubble temperature dependence. Left panel: “simple
ALP” misalignment regime (i.e. scalar with constant mass) corresponding to one
stage of oscillations. Right panel: trapped misalignment regime with two
stages of oscillations. Note the degeneracy in ${\cal N}$ of the trajectories
in the medium and high-temperature regimes. The value
$\gamma_{k>0}\equiv\gamma^{\prime}=1/10$ has been used.

There are two qualitatively distinct regimes, depending on whether this
crossing happens after or before the axion mass abruptly drops at $T_{\rm
QCD}$ (as illustrated in Fig. 3). Let us denote by $f_{a}^{\rm crit}$ the
value of $f_{a}$ for which the canonical QCD axion would start to oscillate
precisely at $T_{\rm QCD}$,

$f_{a}^{\rm crit}\equiv\frac{m_{\pi}f_{\pi}}{m_{\rm
crit}}\frac{\sqrt{z}}{1+z}\simeq 2\times
10^{17}\,\text{GeV},\quad\text{where}\quad m_{\rm crit}=1.67\,H(T_{\rm
QCD})\simeq 2.9\times 10^{-11}\,\text{eV}\,.$ (3.19)

Whenever $m_{a}>m_{\rm crit}$ ($m_{a}<m_{\rm crit}$) Hubble crossing takes
place just (more than) once in the cosmic trajectory. The factor $1.67$ for
$m_{\rm crit}$ stems from our criteria for the onset oscillations that better
fits the numerical solution to the EOM, as argued further below. Strictly
speaking, the value of $f_{a}^{\rm crit}$ for the $Z_{\cal N}$ axion is
slightly larger than that in Eq. (3.19), see Eq. (3.13); this small difference
will be taken into account in the numerical results and figures, but
disregarded in the discussion below for the different regimes of oscillation.
In Fig. 3, $f_{a}=f_{a}^{\rm crit}$ corresponds to the case in which the (
green) Hubble line would cross the canonical QCD axion trajectory (i.e. ${\cal
N}=1$,  orange) exactly at the flattening point of the latter, i.e. $\sim
T_{\rm QCD}$.

####  3.2.1 One stage of oscillations: simple ALP regime

For $f_{a}>f_{a}^{\rm\bf{crit}}$ the $T=0$ potential is already developed when
the axion starts to oscillate, see left panel of Fig. 3. Oscillations only
take place around the true minimum $\theta_{a}=0$, and the relic density
corresponds to that of a “simple ALP” with constant mass. This single regime
of oscillations with constant mass applies as well to the standard QCD axion
scenario (i.e. ${\cal N}=1$ in Fig. 3) for large enough $f_{a}$. The evolution
of the axion field is straightforward and goes as follows:

#####  $T>T_{1}$.

At high temperatures the axion field is frozen at a given and arbitrary
initial misalignment angle $\theta_{1}$ due to Hubble friction, $H\gg
m_{a}(T)$. As the Universe cools down the Hubble parameter diminishes, until
the friction term intercepts the $m_{a}(T)$ trajectory at a temperature
$T_{1}$ such that

$1.67H(T_{1})=m_{a}(T_{1})\equiv m_{1}\,.$ (3.20)

Since $T_{1}<T_{\rm QCD}$, the true vacuum potential has already developed
when the axion field starts to oscillate. In consequence, its mass at $T_{1}$
is already the zero-temperature mass:

$m_{1}=m_{a}\,,\qquad\text{with}\qquad\theta_{a}(T_{1})=\theta_{1}\,,\quad\dot{\theta}_{1}=0\,.$
(3.21)

#####  $T_{1}>T$.

The axion fields keeps oscillating around $\theta_{a}=0$ until nowadays and
its mass is almost constant, $\dot{m_{a}}\simeq 0$. It is adequate to apply
the WKB approximation provided
$m_{a}\gg\big{\\{}\dot{m}_{a}/{m}_{a},\dot{H}/H\big{\\}}$, i.e. provided the
oscillations are much faster than both the expansion rate of the Universe and
the rate at which the mass changes, which is the case here. The adiabatic
approximation predicts the existence of a conserved quantity [5, 6, 7], an
adiabatic invariant $N_{a}$ that can be interpreted as the comoving number of
non-relativistic axion quanta, defined for a generic axion mass $m_{a}$ as

$N_{a}\equiv\frac{\rho_{a}a^{3}}{m_{a}}={\rm
const.}\,,\qquad\text{with}\qquad{\rho_{a}}\equiv\frac{1}{2}{\dot{\theta}}_{a}^{2}{f_{a}^{2}}+\frac{1}{2}m_{a}^{2}{f_{a}^{2}}\theta_{a}^{2}\,,$
(3.22)

where $a$ denotes the scale factor and $\rho_{a}$ the axion energy density.
For the regime under discussion, it follows from $N_{a}$ conservation between
$T_{1}$ and today that the current energy density $\rho_{a,\,0}$ can be
expressed as

$\displaystyle\rho_{a,\,0}=m_{a}\frac{N_{a}}{a_{0}^{3}}\simeq\frac{1}{2}\,m_{a}^{2}\,(\theta_{1}\,f_{a})^{2}\,\left(\frac{a_{1}}{a_{0}}\right)^{3}\,,$
(3.23)

where $a_{0}$ and $a_{1}$ denote respectively the value of the scale factor
today and at the onset of oscillations, $a_{1}\equiv a(T_{1})$.777Whenever the
initial misalingment angle is large, the relic density in Eq. 3.23 needs to be
corrected by the so-called anharmonicity factor, see Appendix B. The latter
can be expressed in turn in terms of the axion mass using Friedmann equations
and the conservation of entropy, $a_{1}\propto 1/\sqrt{m_{a}M_{\rm Pl}}$,
allowing us to obtain the ratio of the axionic relic density to the current DM
abundance $\rho_{\rm DM}$ as a function of $m_{a}$, $\theta_{1}$ and
$f_{a}$,888The slight discrepancy with respect to the result in Ref. [73] is
due to the use of the updated value of $\rho_{\rm DM}\simeq
1.26\,{\mathrm{keV}}/{\mathrm{cm}^{3}}$ from Planck 2018 data [60] and a
different choice for the onset of the oscillations that better fits the
numerical result ($1.67\,H=m_{a}$ instead of $3H=m_{a}$, see e.g. Refs. [74,
75]).

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm DM}}\simeq
5.7\sqrt{\frac{m_{a}}{\mathrm{eV}}}\left(\frac{\theta_{1}\,f_{a}}{10^{12}\,\mathrm{GeV}}\right)^{2}\mathcal{F}(T_{1})\,,$
(3.24)

where
$\mathcal{F}(T_{1})\equiv\left(g_{*}(T_{1})/3.38\right)^{\frac{3}{4}}\left(g_{s}(T_{1})/3.93\right)^{-1}$
is an $\mathcal{O}(1)$ factor and $\rho_{\rm DM}\simeq
1.26\,{\mathrm{keV}}/{\mathrm{cm}^{3}}$ from Planck 2018 data [60]. For the
$Z_{\cal N}$ scenario under discussion, it follows from Eq. 1.4 that this
ratio can be rewritten in the large ${\cal N}$ limit as

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm DM}}\simeq 0.48\,{\cal
N}^{3/8}{z^{{\cal
N}/4}}\,\theta_{1}^{2}\,\left(\frac{f_{a}}{10^{13}\mathrm{GeV}}\right)^{3/2}\mathcal{F}(T_{1})\,.$
(3.25)

The initial axion misalignment $\theta_{1}$ can a priori take any value in the
$[-\pi,\pi)$ interval and the axion field will roll down to the closest
minimum999Note that the relevant axion field displacement that should enter in
Eq. (3.25) corresponds to the field distance to the closest minimum and it is
only when the closest minimum is $\theta_{a}=0$ that the field distance
coincides with $\theta_{1}$. among the ${\cal N}$ possibilities. Thus it is
with probability $1/{\cal N}$ that the final $\theta$-parameter will
correspond to the CP-conserving point $\theta=0$. In other words, in order for
the simple ALP regime to solve the strong CP problem and also account for DM
it is necessary and sufficient that the initial misalignment angle lies in the
interval $\theta_{1}\in[-\pi/{\cal N},\pi/{\cal N})$ within the simple ALP
regime, for any $f_{a}>2\times 10^{17}\,\text{GeV}$.

Eq. (3.23) also indicates that the simple ALP solutions which account for the
relic DM obey (see Fig. 7)

$\displaystyle m_{a}f_{a}^{4}\propto{\rm const.}\,,$ (3.26)

to be compared with the $Z_{\cal N}$ axion mass relation, that for a given
${\cal N}$ predicts $m_{a}f_{a}\propto{\rm const.}$ This behaviour is depicted
in the $\\{m_{a},1/f_{a}\\}$ plane in Fig. 6 by continuous superimposed lines
(for different values of $\theta_{1}$), together with some experimentally
excluded regions as a reference (while the complete set of present and
projected sensitivity regions is depicted in Fig. 11). It follows from Fig. 6
that the simple ALP $Z_{\cal N}$ solutions can solve both the strong CP
problem and account for DM mass down to the region of fuzzy DM ($m_{a}\sim
10^{-22}$ eV) and for any value of ${\cal N}\geq 3$.

####  3.2.2 Two stages of oscillation: trapped misalignment

For $f_{a}<f_{a}^{\rm crit}$, two different potential minima develop during
the Universe evolution from high to low temperatures, see Fig. 3 (right
panel). In a nutshell:

  * •

A first phase of oscillations takes place around $\theta_{a}=\pi$ at
temperatures above $T_{\rm QCD}$, because temperature effects are important
then and the axion mass is unsuppressed.

  * •

At $T_{\rm QCD}$ the true minimum $\theta_{a}=0$ develops (while $\pi$ becomes
a maximum). The axion mass suddenly becomes exponentially suppressed due to
the $Z_{\cal N}$ symmetry. Oscillations around $\theta_{a}=0$ will start
whenever the kinetic energy becomes smaller than the height of the potential
barrier.

These two stages are thus separated by a drastic (non-adiabatic) modification
of the potential. This novel axion DM production mechanism is named _trapped
misalignment_. Its cosmic evolution is illustrated in Figs. 4 and 5 for toy-
model trapped trajectories ( blue), versus that of a QCD-like axion ( orange)
with the same zero temperature mass. Let us expatiate next on the evolution
details.

#####  $T>T_{1}$.

In this period, the behaviour is alike to that for a simple ALP, see Eq.
(3.20). Before Hubble crossing at $T_{1}$, the axion field is _frozen_ and
remains constant at an arbitrary initial misalignment angle $\theta_{1}$.

#####  $T_{1}>T>T_{\rm QCD}$.

The finite temperature effective potential is relevant in this range, see Eqs.
(3.1) and (3.7), and the initial axion velocity is $\dot{\theta}_{1}=0$. Thus,
the axion “is trapped” in this _first stage of oscillations_ around
$\theta_{a}\sim\pi$ from $T_{1}$ until $T_{\rm QCD}$. Its mass is unsuppressed
due to the explicit $Z_{\cal N}$ breaking by the thermal background:
$m_{a}(T)$ raises progressively until it stabilizes at $\sim
m_{a}^{\rm{QCD}}$.

Trapping causes a delay of the start of oscillations around the true minimum
$\theta_{a}=0$, as compared to $T_{1}$ for the canonical QCD axion and the
simple ALP scenarios. The consequence is an enhancement of the DM density, as
the dilution time is thus shortened.

In Figs. 4 and 5 the evolution for a QCD-like axion (i.e. ${\cal N}=1$,
orange) is shown to oscillate around $\theta_{a}=0$ since well before $T_{\rm
QCD}$. In contrast, for ${\cal N}\geq 3$ ( blue) the oscillations around
$\theta=\pi$ take place from $T_{1}$ down to $T_{\rm QCD}$. The crucial impact
of the trapped stage is thus to delay the onset of oscillations around the
true minimum, and to set the initial value of $\theta_{a}$ and
$\dot{\theta_{a}}$ at $T_{\rm QCD}$,

$\theta_{\rm tr}\equiv\theta_{a}(T_{\rm QCD})\sim\pi,\,\quad\dot{\theta}_{\rm
tr}\equiv\dot{\theta}_{a}(T_{\rm QCD})\,.$ (3.27)

It is not possible to predict the exact value of $\dot{\theta}_{\rm tr}$.
However, an order of magnitude estimate stems from the energy density when the
temperature approaches $T_{\rm QCD}$, applying $N_{a}$ conservation to the
adiabatic interval $[T_{1},T_{\rm QCD})$,

$\displaystyle\rho_{a,\,{\rm QCD}}\equiv\rho_{a}(T_{\rm QCD})=m_{a,\pi}^{\rm
QCD}\frac{N_{a}}{a_{\rm QCD}^{3}}=\frac{1}{2}m_{1}m_{a,\pi}^{\rm
QCD}\left(\frac{a_{1}}{a_{\rm
QCD}}\right)^{3}(\theta_{1}-\pi)^{2}f_{a}^{2}\,,$ (3.28)

which shows that the precise value of $\rho_{a,\,{\rm QCD}}$ depends on the
arbitrary initial misalignment angle with respect to the high-temperature
minimum in $\pi$: $(\theta_{1}-\pi)$. This dependence is inherited by the mean
velocity at the end of this period, which follows from the equality of the
mean kinetic and potential energies for a harmonic oscillator,

$\displaystyle\sqrt{\langle\dot{\theta}^{2}_{\rm
tr}\rangle}=\frac{1}{\sqrt{2}}\sqrt{m_{1}m_{a,\pi}^{\rm
QCD}}\left(\frac{a_{1}}{a_{\rm QCD}}\right)^{3/2}|\theta_{1}-\pi|\,.$ (3.29)

#####  $T=T_{\rm QCD}$.

At this point, the low-temperature potential for the $Z_{\cal N}$ axion
scenario develops (see Eqs. (1.3), (1.4) and (3.5)), and thus $\theta_{a}=\pi$
becomes a maximum. The abrupt exponential suppression of the axion mass was
shown in Fig. 3 for different values of ${\cal N}$. It illustrates that the
two stages of oscillation are separated by a sudden, inherently non-adiabatic,
modification of the potential as $\dot{m}_{a}/m_{a}\gg m_{a}$. The energy
density of the axion just after the transition at $T_{\rm QCD}$ reads

$\displaystyle{\rho_{a,{\rm tr}}}=\frac{1}{2}{\dot{\theta}}_{\rm
tr}^{2}{f_{a}^{2}}+\frac{m_{a}^{2}{f_{a}^{2}}}{{\cal
N}^{2}}\big{[}1-\cos({\cal
N}\theta_{tr})\big{]}\,\simeq\frac{1}{2}{\dot{\theta}}_{\rm
tr}^{2}{f_{a}^{2}}+2\frac{m_{a}^{2}{f_{a}^{2}}}{{\cal N}^{2}}\,,$ (3.30)

that is, the total energy density at this point is larger or equal than the
maximum potential energy (the height of the barrier
$2{m_{a}^{2}{f_{a}^{2}}}/{{\cal N}^{2}}$). What happens next depends crucially
on the value of $\dot{\theta}_{\rm tr}$, which acts as initial condition for
the subsequent period. It determines at which temperature the _second stage of
oscillations_ begins:

  * •

For very small axion velocity, $\dot{\theta}_{\rm tr}\ll 2m_{a}/{\cal N}$,
oscillations around one minimum may start closely after $T_{\rm QCD}$. This
case will be denoted below “pure trapped misalignment”.

  * •

For large enough $\dot{\theta}_{\rm tr}$ so that the kinetic energy dominates
${\rho_{a,{\rm tr}}}$, $\dot{\theta}_{\rm tr}\gg 2m_{a}/{\cal N}$, the axion
field will roll over the top of the barriers for a long time before its
oscillations start at a temperature sensibly lower than $T_{\rm QCD}$. This
case will be denoted as “trapped+kinetic misalignment”.

#####  3.2.2.1 Pure trapped misalignment: $\dot{\theta}_{\rm tr}\sim 0$

#####  $T_{\rm QCD}>T$.

Were the axion velocity exactly zero at $T_{\rm QCD}$, $\dot{\theta}_{\rm
tr}=0$, the setup would not be viable: the axion would simply roll to the
closest minimum $\theta_{a}=\frac{\pi}{{\cal N}}({\cal N}-1)$. This is
precluded by the experimental limit on the nEDM that implies $\theta\lesssim
10^{-10}$.

Nevertheless, $\dot{\theta}_{\rm tr}$ may be very small but not zero. Even a
tiny kinetic energy may suffice to allow the axion to go over some potential
barriers, because the misalignment set by the trapping corresponds to a
maximum of the zero temperature potential. In Sect. 3.2.3 the minimum kinetic
energy necessary to roll over $\mathcal{O}({\cal N})$ potential
barriers101010Strictly speaking, about ${\cal N}/2$ barriers must be overflown
to approach the $\theta_{a}=0$ minimum, but we will obviate this finesse as it
does not change the order of magnitude estimates intended here. will be
estimated and shown to be attained in most of the parameter space. The axion
field will then end up in the CP-conserving minimum $\theta_{a}=0$ with
$1/{\cal N}$ probability.

An estimation of the time required for the axion to go from one maximum to the
next is $\delta t\sim 2\pi/m_{a}$. To overfly only $\sim{\cal N}$ barriers a
very short time is required because in this regime $m_{a}\ll H$, i.e. the
difference between $T_{\rm QCD}$ and the temperature at the onset of
oscillations around the true minimum is negligible. In conclusion, for
$\dot{\theta}_{\rm tr}\ll 2m_{a}/{\cal N}$ the kinetic energy is not large
enough to extend beyond $\sim T_{\rm QCD}$ the delay time accumulated during
the trapped period (see Eq. 3.30). The value of $\theta_{a}$ when it starts to
oscillate from the top of the barrier closest to $\theta_{a}=0$ is

$\theta_{2}\sim\frac{\pi}{{\cal N}}\,.$ (3.31)

The evolution is illustrated in Fig. 4, which depicts the numerical solution
for a toy example exhibiting pure trapped misalignment ( blue) and compares it
with that for a QCD-like axion ( orange) with the same zero temperature mass.
The figure reflects the enhancement of the relic DM density in the trapped
case. The latter can be easily computed analogously to that for the usual
misalignment mechanism. After the transition at $T_{\rm QCD}$, the adiabatic
approximation is valid again since $\dot{m}_{a}\simeq 0$, $m_{a}\gg H$. From
the $N_{a}$ conservation law –Eq.(3.22), and the initial conditions at the
onset of the second stage of oscillations $\big{\\{}\theta_{\rm
2}\simeq\pi/{\cal N},\,\dot{\theta}_{\rm 2}\simeq 0\big{\\}}$, it follows that
the current axionic relic density $\rho_{a,0}$ is given in the pure trapped
case by

$\displaystyle\rho_{a,0}\big{|}_{\rm
tr}=m_{a}\frac{N_{a}}{a_{0}^{3}}=\frac{1}{2}\,m_{a}^{2}\,\left(\frac{\pi}{{\cal
N}}f_{a}\right)^{2}\,\left(\frac{a_{\rm QCD}}{a_{0}}\right)^{3}\,,$ (3.32)

where $a_{\rm QCD}\equiv a(T_{\rm QCD})$. For a given final vacuum mass
$m_{a}$, it demonstrates an enhancement of the relic density in the scenario
with a $Z_{\cal N}$ trapped axion vs. that with a simple ALP, for two main
reasons: i) the relative scale dependence is larger by a factor $({a_{\rm
QCD}}/{a_{1}})^{3}$ and does not depend on the axion mass; ii) the initial
misalignment angle is fixed to $\theta_{\rm tr}\sim\pi/{\cal N}$, while for
the simple ALP it is arbitrary within the interval $[-\pi/{\cal N},\pi/{\cal
N})$. Note that the current relic density in the pure trapped scenario becomes
independent of the initial misalignment $\theta_{1}$, unlike the case of the
simple ALP and QCD-like pre-inflationary scenarios.

Figure 4: Axion evolution for the Pure Trapped misalingment vs. that for a
canonical QCD axion, illustrating the relic density enhancement for the same
final zero-temperature mass $m_{a}$.

For the $Z_{\cal N}$ reduced-mass axion under discussion, the product
$m_{a}f_{a}$ –and thus $\rho_{a,0}$ in Eq. (3.32)– is a function of just one
unknown: ${\cal N}$, see $m_{a}$ in Eqs. (1.4) and (3.12). Substituting this
dependence in the ratio

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm DM}}\bigg{|}_{\rm tr}\simeq
5.6\sqrt{\frac{m_{a}}{\mathrm{eV}}}\left(\frac{m_{a}}{m_{\rm
crit}}\right)^{3/2}\left(\frac{\pi}{{\cal
N}}\,\frac{f_{a}}{10^{12}\,\mathrm{GeV}}\right)^{2}\mathcal{F}(T_{\rm
QCD})\,,$ (3.33)

it follows that

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm DM}}\bigg{|}_{\rm tr}\sim 1.5\times
10^{7}\,\frac{z^{{\cal N}}}{\sqrt{{\cal N}}}\,.$ (3.34)

Thus in the pure trapped regime the value of ${\cal N}$ compatible with
simultaneously solving the strong CP problem and accounting for the ensemble
of the DM density is

${\cal N}\sim 21\,.$ (3.35)

This line of constant relic density is highlighted in Fig. 6 ( purple) for
relatively small $f_{a}$ values. Indeed, Eq. (3.32) predicts that the relic
density solution in the $\\{m_{a},f_{a}\\}$ plane has to be a line parallel to
those defining the canonical QCD and $Z_{\cal N}$ axion bands, because in all
these cases the parametric dependence is $m_{a}f_{a}={\rm const.}$ The
estimate in Eq. (3.35) disregards anharmonicity corrections and in order to
qualitatively account for the latter ones the axion DM region will be
represented as a band in the figures of Sect. 4.111111The size of this
correction will be different in general for the pure trapped case and the
simple ALP regime since in the former $H\ll m_{a}$, see Appendix B.

#####  3.2.2.2 Trapped+kinetic misalignment: sizeable $\dot{\theta}_{\rm tr}$

#####  $T_{\rm QCD}>T>T_{2}$.

If the axion velocity at $T_{\rm QCD}$ is large enough, $\dot{\theta}_{a}\gg
2m_{a}/{\cal N}$,121212Equivalently, if the axion kinetic energy at the end of
the trapped period $\dot{\theta}_{\rm tr}^{2}f_{a}^{2}/2$ is much larger than
the maximum potential energy $2m_{a}^{2}f_{a}^{2}/{\cal N}^{2}$. the axion
field has enough kinetic energy to roll many times over the barriers before it
starts to oscillate around some minimum, triggering the kinetic misalignment
mechanism. The onset of oscillations is thus delayed until much later than
$T_{\rm QCD}$.

In this regime, the WKB (adiabatic) approximation is not valid anymore and the
axion energy density is diluted as $\propto a^{-6}$ due to the expansion of
the Universe (as illustrated in Fig. 5). This can be shown by noting that
$\dot{\theta}_{a}f_{a}^{2}$ is the Noether charge density associated to the PQ
symmetry [56] and therefore it decays as $\dot{\theta}_{a}\propto a^{-3}$ as
long as the axion mass can be neglected.131313Alternatively, the dilution of
the axion velocity can be obtained from the axion EOM in the massless limit:
$\ddot{\theta}_{a}+3H\dot{\theta}_{a}=0\,\implies\dot{\theta}_{a}a^{3}={\rm
const}$. In consequence, a comoving PQ charge $q_{\rm kin}$ can be defined
which is a conserved quantity,

$\displaystyle q_{\rm kin}=\dot{\theta}_{a}a^{3}={\rm const}\,.$ (3.36)

This regime of rapidly decreasing energy density ends up at a temperature
$T_{2}$ at which the kinetic energy becomes of the order of the height of the
barrier, and thus the axion can no longer overcome it,

$\displaystyle\frac{1}{2}\dot{\theta}^{2}_{a}(T_{2})f_{a}^{2}=2\,\frac{m_{a}^{2}f_{a}^{2}}{{\cal
N}^{2}}\quad\implies\quad\dot{\theta}_{a,\,2}\equiv\dot{\theta}_{a}(T_{2})=2\,\frac{m_{a}}{{\cal
N}}\,.$ (3.37)

At this point the axion will start a second stage of oscillations, see Fig. 5.
Enforcing the conservation of $q_{\rm kin}$ from $T_{\rm QCD}$ until $T_{2}$,
and using Eq. (3.37), it follows that the scale factor at the onset of the
second stage of oscillations $a_{2}\equiv a(T_{2})$ can be expressed as

$\displaystyle q_{\rm kin}=\dot{\theta}_{\rm tr}\,a_{{\rm
QCD}}^{3}=2\,\frac{m_{a}}{{\cal N}}a_{2}^{3}\quad\implies\quad
a_{2}=\left(\frac{{\cal N}\,\dot{\theta}_{\rm
tr}}{2\,m_{a}}\right)^{1/3}a_{{\rm QCD}}\,.$ (3.38)

Figure 5: Axion evolution for the Trapped+Kinetic misalignment vs. that for a
canonical QCD axion, illustrating the strong enhancement of the relic density
enhancement for the same final zero-temperature mass $m_{a}$.

#####  $T_{2}>T$.

The axion finally starts the second stage of oscillations around one of the
${\cal N}$ zero-temperature minima at random. Thus, the $1/{\cal N}$
probability of the $Z_{\cal N}$ axion scenario [31, 32] to solve the strong CP
problem (i.e. the final $\theta$-parameter to corresponding the CP-conserving
point $\theta=0$) will be maintained when accounting for DM.

In order to determine the axion relic density today, the adiabatic
approximation is again appropriate after $T_{2}$. The comoving number of
axions $N_{a}^{\prime}$ is constant during this final period (and different
from that during the first oscillation period $N_{a}$). Using $q_{\rm kin}$ in
Eq. (3.38) and the mean axion velocity at $T_{\rm QCD}$ in Eq. 3.29,
$N_{a}^{\prime}$ is well approximated by [56]

$\displaystyle
N^{\prime}_{a}=\frac{\rho_{a,2}a_{2}^{3}}{m_{a}}=2\,\frac{m_{a}f_{a}^{2}a_{2}^{3}}{{\cal
N}^{2}}=\frac{\dot{\theta}_{\rm tr}}{{\cal N}}\,f_{a}^{2}\,a_{\rm QCD}^{3}\,.$
(3.39)

Finally the relic density today,

$\displaystyle\rho_{a,0}\big{|}_{\rm
tr+kin}\simeq\,\frac{m_{a}}{a_{0}^{3}}N_{a}^{\prime}\,,$ (3.40)

can be written as

$\displaystyle\rho_{a,0}\big{|}_{\rm
tr+kin}=C\,\frac{m_{a}\sqrt{m_{1}m_{a,\pi}^{\rm QCD}}}{\sqrt{2}\,{\cal
N}}\left(\frac{\sqrt{a_{1}a_{\rm
QCD}}}{a_{0}}\right)^{3}|\theta_{1}-\pi|\,f_{a}^{2}\,,$ (3.41)

where $C$ denotes a deviation from the analytical –exactly adiabatic and
harmonic– estimations above, to be computed numerically. For our setup, we
find numerically $C\simeq 2$, which confirms the result first obtained in Ref.
[56]; this correction will be included in all figures below for the
kinetic+trapped mechanism.

The comparison in Fig. 5 between the trapped+kinetic evolution ( blue) with
that for a QCD-like axion ( orange) shows how the delay produced by the
kinetic-dominated period results in an enhancement of the current axionic
relic density, which is even stronger than in the pure trapped scenario. The
period of kinetic misalignment has been shortened in the figure for
illustration purposes.

The comparison of the resulting axionic relic density with that for the pure
trapped scenario in Eq. 3.32 shows two competing effects. In fact, for the
kinetic+trapped scenario:

  * •

the dilution factor is smaller, $(\sqrt{a_{1}a_{\rm QCD}}/{a_{0}})^{3}$ vs.
$(a_{\rm QCD}/{a_{0}})^{3}$;

  * •

the relevant mass scale is larger, $\sqrt{m_{1}m_{a,\pi}^{\rm QCD}}$ vs.
$m_{a}$.

Hence, it will depend on the specific point of the parameter space whether the
kinetic+trapped or the pure trapped mechanism gives the dominant contribution
(see Section 3.3).

The ratio of the predicted axion relic density to the observed DM density can
be obtained following same the procedure applied for the simple ALP regime,
yielding

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm DM}}\bigg{|}_{\rm tr+kin}\simeq
7.9\sqrt{\frac{m_{a}}{\mathrm{eV}}}\,\frac{\sqrt{m_{a}\,m_{a,\pi}^{\rm
QCD}}}{\big{(}m_{\rm
crit}^{3}m_{1}\big{)}^{1/4}}\left[\frac{\,f_{a}}{10^{12}\,\mathrm{GeV}}\right]^{2}\,\frac{\left|\theta_{1}-\pi\right|}{{\cal
N}}\,\mathcal{F}_{\rm kin,1}(T_{1})\,,$ (3.42)

where $\mathcal{F}_{\rm kin,1}(T_{1})$ is a function of the degrees of freedom
at play, which can be found in Appendix A together with further details. The
relic density depends on the temperatures of all the different SM copies via
the value of $m_{1}$ (i.e. the $\gamma_{k}$ values, see Eqs. (2.3) and (2.8)).
In the simplified case in which all the copies of the SM have the same
temperature $T_{k\neq 0}\sim T^{\prime}$, Eq. (3.42) can be written –for large
${\cal N}$ and for the largest possible kinetic misalignment
contribution141414Here the maximum axion velocity at the end of the trapped
period is considered, i.e. twice the mean axion velocity in Eq. 3.29.– as

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm DM}}\bigg{|}_{\rm
tr+kin,\,max}\simeq
0.64\left(\frac{\,f_{a}}{10^{9}\,\mathrm{GeV}}\right)^{5/6}\,\frac{z^{{\cal
N}/2}}{{\cal N}^{1/3}}\,\left|\theta_{1}-\pi\right|\,\mathcal{F}_{\rm
kin,\,2}(T_{1})\,,$ (3.43)

where $\mathcal{F}_{\rm kin,\,2}(T_{1})$ is another function of the degrees of
freedom at play, see Appendix A.

It also follows from the results above that the trapped+kinetic DM solutions
behave in the $\\{m_{a},1/f_{a}\\}$ plane as

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm DM}}\bigg{|}_{\rm tr+kin}\propto
m_{a}f_{a}^{19/12}\,,$ (3.44)

in contrast to $m_{a}f_{a}={\rm const.}$ for the pure trapped case in Eq.
(3.32). The trapped+kinetic trajectories that account for DM are illustrated
by double lines in Fig. 6, for different values of the initial misalignment
angle $\theta_{1}$.

####  3.2.3 Solving both the strong CP problem and the nature of DM

We estimate here the minimum kinetic energy necessary in the trapped regime
for the axion to roll over at least $\mathcal{O}({\cal N})$ maxima after
$T_{\rm QCD}$. The result will be of general interest whenever
$f_{a}<f_{a}^{\rm crit}$, and of practical relevance mainly for the pure
trapped misalignment.

Since the trapping forces $\theta_{\rm tr}\simeq\pi$, the corresponding total
energy density is always larger than the maximum potential energy (the height
of the barrier $2{m_{a}^{2}{f_{a}^{2}}}/{{\cal N}^{2}}$). As the axion rolls
over $\sim{\cal N}$ different maxima, the energy density is diluted due to the
expansion of the Universe,

$\displaystyle\rho_{a}(a_{{\rm QCD}+{\cal N}})=\rho_{a,{\rm
tr}}\left(\frac{a_{\rm QCD}}{a_{{\rm QCD}+{\cal N}}}\right)^{p}\,,$ (3.45)

where $a_{{\rm QCD}+{\cal N}}$ is the scale factor when the axion has rolled
over just ${\cal N}$ maxima after $T_{\rm QCD}$, and $p$ parametrizes how
strong is the dilution (we find numerically that $p\sim 1.5$). The minimal
kinetic energy $K_{\rm min}$ necessary for the axion to be able to roll over
just ${\cal N}$ maxima before oscillating must saturate the condition

$\displaystyle\left(K_{\rm min}+2\frac{m_{a}^{2}{f_{a}^{2}}}{{\cal
N}^{2}}\right)\left(\frac{a_{\rm QCD}}{a_{{\rm QCD}+{\cal
N}}}\right)^{p}=2\,\frac{m_{a}^{2}{f_{a}^{2}}}{{\cal N}^{2}}\,.$ (3.46)

Given the Hubble parameter in a radiation dominated Universe
$H(t)=\frac{1}{2t}$ and the Friedmann equation in Eq. 2.4, it follows that the
temperature $T_{{\rm QCD}+{\cal N}}$ after overcoming ${\cal N}$ maxima
(corresponding to a time difference of $\delta t\sim 2\pi{\cal N}/m_{a}$, see
discussion before Eq. (3.31)) is

$\displaystyle T_{{\rm QCD}+{\cal N}}=\left[\frac{4\pi^{3}}{45}g_{*}(T_{\rm
QCD})\right]^{1/4}\sqrt{\frac{M_{\rm Pl}}{2(t_{\rm QCD}+{2\pi{\cal
N}}/{m_{a}})}}\simeq T_{\rm QCD}-8.2\,\frac{2\pi{\cal N}}{m_{a}}\frac{T_{\rm
QCD}^{3}}{M_{\rm Pl}}\,,$ (3.47)

where $t_{\rm QCD}$ is the time corresponding to $T_{\rm QCD}$, and the
approximation ${2\pi{\cal N}}/{m_{a}}\ll t_{\rm QCD}$151515Or equivalently
${2\pi{\cal N}\,T^{2}_{\rm QCD}}\ll{m_{a}\,M_{\rm Pl}}$. has been used, as in
this regime $H\ll m_{a}$. The corresponding dilution factor reads

$\displaystyle\left(\frac{a_{\rm QCD}}{a_{{\rm QCD}+{\cal
N}}}\right)^{p}\simeq\left(\frac{T_{{\rm QCD}+{\cal N}}}{T_{\rm
QCD}}\right)^{p}\simeq 1-8.2\,p\,\frac{2\pi{\cal N}}{m_{a}}\frac{T_{\rm
QCD}^{2}}{M_{\rm Pl}}\,.$ (3.48)

Substituting it in Eq. 3.46, the minimum kinetic energy required is
determined,

$\displaystyle K_{\rm min}\simeq 8.2\,p\,\frac{4\pi}{{\cal
N}}\,\frac{m_{a}T_{\rm QCD}^{2}}{M_{\rm Pl}}f_{a}^{2}\,.$ (3.49)

To evaluate how probable it is to simultaneously solve the strong CP problem
and explain DM, $K_{\rm min}$ needs to be compared with the mean kinetic
energy generically expected at $T_{\rm QCD}$. It follows from Eq. (3.29) that

$\displaystyle K_{\rm tr}\equiv\frac{1}{2}{\langle\dot{\theta}^{2}_{\rm
tr}\rangle}\,f_{a}^{2}=\frac{1}{4}{m_{1}m_{a,\pi}^{\rm
QCD}}\left(\frac{a_{1}}{a_{\rm
QCD}}\right)^{3}(\theta_{1}-\pi)^{2}f_{a}^{2}\,.$ (3.50)

In the case in which all mirror copies but the SM have the same temperature
$T^{\prime}$ and taking for example $\gamma^{\prime}=0.2$, the latter equation
can be rewritten as

$\displaystyle K_{\rm tr}\simeq\frac{1}{4}\left(\frac{T_{\rm
QCD}^{14}(m_{a,\pi}^{\rm QCD})^{5}}{M_{\rm
Pl}^{7}}\right)^{1/6}\,(\theta_{1}-\pi)^{2}\,f_{a}^{2}\times\kappa\,,$ (3.51)

where $\kappa$ is a numerical factor that ranges in the interval $(4-8)$. This
implies that the axion field will roll over at least ${\cal N}$ barrier tops
whenever

$\displaystyle m_{a}\lesssim\frac{{\cal N}}{4\pi\,8.2\,p}\left(\frac{T_{\rm
QCD}}{\sqrt{M_{\rm Pl}}}\right)^{1/3}(m_{a,\pi}^{\rm
QCD})^{5/6}|\theta_{1}-\pi|^{2}\,.$ (3.52)

We compare next the parametric dependence of all DM axion solutions –which
also solve the strong CP problem– discussed up to this point. They are
depicted in the $\\{m_{a},1/f_{a}\\}$ plane of Fig. 6 and for several values
of $\theta_{1}$:

Figure 6:  $Z_{\cal N}$ axion solutions that solve both the strong CP problem
and account for DM. The superimposed single continuous lines indicate either
simple ALP solutions or pure trapped ones, while double lines denote
trapped+kinetic misalignment regions. The tenuous oblique lines for various
values of ${\cal N}$ indicate $Z_{\cal N}$ axion solutions to the strong CP
problem. Solid colors denote regions experimentally excluded by axion-nEDM
[76] or astrophysical data. Further constraints and prospects can be seen in
Fig. 11.

  * •

In the excluded grey area the strong CP problem is not solved since the
condition in Eq. 3.52 is not satisfied. This precludes any DM solution with
${\cal N}\lesssim 9$, although only for $m_{a}\gtrsim m_{\rm crit}\simeq
4\times 10^{-11}\,\text{eV}$.

  * •

The stars indicate for each trajectory the point below which $f_{a}$ is large
enough ($f_{a}>f_{a}^{\rm{crit}}$) so that the oscillations around the true
minimum start after the zero temperature potential is fully developed. This is
the “simple ALP” case with DM solutions obeying $m_{a}f_{a}^{4}\propto{\rm
const.}$

  * •

Above the stars, double (single) continuous lines indicate trapped+kinetic
(pure trapped) trajectories. For large values of $f_{a}$ within the
$f_{a}<f_{a}^{\rm{crit}}$ range, the trapped+kinetic solution tends to
dominate the parameter space, while for small $f_{a}$ the pure trapped
mechanism sets the relic density.

  * •

The kinetic kick received by a given trapped+kinetic trajectory when its
$f_{a}$ value is smaller than $f_{a}^{\rm crit}$ (but very close to it) can be
clearly seen. Above that point, the slope of the trajectory obeys
$m_{a}f_{a}^{19/12}={\rm const.}$, as explained around Eq. (3.44).

  * •

The pure trapped regime, that populates the top of the figure for small enough
$f_{a}$, corresponds to the continuous ${\cal N}\sim 21$ line, whose slope
obeys $m_{a}f_{a}={\rm const.}$

Overall, Fig. 6 shows that good solutions to both the strong CP problem and
the nature of DM exist within the $Z_{\cal N}$ axion scenario whenever
$m_{a}\lesssim 10^{-4}$ eV, down to the fuzzy DM region ($m_{a}\sim 10^{-22}$
eV). In the range $m_{\rm crit}\lesssim m_{a}\lesssim 10^{-4}$ eV the number
of worlds must be ${\cal N}\sim 21$. For the lighter axion scenarios,
$m_{a}\lesssim m_{\rm crit}$, values of ${\cal N}$ up to $\sim 65$ are
possible.

Finally, note that this analysis has focused only on the axion zero mode as a
first step. Within the trapped mechanism the axion field self-interactions may
become important whenever the axion is not close to the minimum. Non-
linearities will induce the production of higher momentum axion quanta, i.e.
axion fragmentation [77]. The consequences of the fragmentation of the axion
field within the trapped misalignment mechanism is left for a future work.

###  3.3 Trapped vs. non-trapped mechanisms

The trapped mechanism identified in this paper is actually more general than
the axion setup with a non-linearly realized $Z_{\cal N}$ shift symmetry
discussed above. Trapped misalignment could arise in fact in a large variety
of QCD axion or ALP scenarios.

In QCD axion models, any additional PQ-breaking source at low energies needs
to be extremely suppressed in order to comply with the nEDM bounds. However,
the high-temperature axion potential is largely unconstrained. Indeed, any PQ-
breaking source which is active only at high energies may induce a
displacement of the potential minima at high temperatures. The axion field can
then get trapped in its displaced minimum for a certain period of time, until
the low-energy QCD potential develops. This can dramatically change the relic
density predicted today, due to the trapped misalignment phase. An enhancement
can be expected whenever the transition between the two phases of the
potential is non-adiabatic161616In regimes which are mainly adiabatic, a
reduction of the axion relic density may result instead [51]. and takes place
at a smaller temperature than the oscillation temperature in absence of
trapping.

$\displaystyle\vline\begin{array}[]{l|ll}\<EMAIL_ADDRESS>\centering\text{Misalignment mechanisms}&\qquad\rho_{a,0}&\quad
g_{aXX}\propto{1}/f_{a}\\\\[5.0pt] \hline\cr\\\\[-10.0pt] \text{Simple
ALP}&\quad\sqrt{m_{a}}f_{a}^{2}&\qquad{\centering
m_{a}^{1/4}\<EMAIL_ADDRESS>\text{Trapped}&\quad
m_{a}^{2}f_{a}^{2}&\qquad m_{a}\\\\[8.0pt] \text{Trapped+Kinetic}&\quad
m_{a}f_{a}^{19/12}&\qquad m_{a}^{12/19}\\\\[5.0pt]
\hline\cr\hline\cr\\\\[-10.0pt] \text{QCD-like }&\quad
m_{a}^{5/6}f_{a}^{2}&\qquad m_{a}^{5/12}\\\\[5.0pt]
\hline\cr\end{array}\vline$ (3.61) Table 1: Parametric dependence of axion
relic density and axion couplings $g_{aXX}$ (for a fixed axionic relic
density), for the different types of misalignment mechanisms discussed in the
text.

The parametric dependence of the axion relic density on the
$\\{m_{a},f_{a}\\}$ plane is a hallmark of the different misalignment
mechanisms discussed above. Table 1 summarizes the results presented earlier
on (see Eqs. (3.26), (3.32) and (3.44)). It also shows the generic $m_{a}$
dependence of the axion-SM interaction couplings $g_{aXX}\sim 1/f_{a}$ (with
$X$ denoting a generic SM field) for each DM production scenario discussed
above. For comparison, the corresponding results for the canonical QCD axion
scenario are shown as well (within the $Z_{\cal N}$ axion there is no regime
where the QCD-like relic density formula applies other than for ${\cal N}=1$).
In Table 1 (and in all figures), for the QCD-like and trapped+kinetic cases
the DIGA value $\alpha=8$ has been assumed (see Eq. 3.2).171717 For arbitrary
$\alpha$, $\rho_{a}\propto m_{a}f_{a}^{{(14+3\alpha)}/{(2(4+\alpha))}}$ in the
trapped+kinetic case, while $\rho_{a}\propto
m_{a}^{(\alpha+2)/(\alpha+4)}f_{a}^{2}$ in QCD-like misalignment, from which
the corresponding $m_{a}$ dependence of $g_{aXX}\propto{1}/f_{a}$ can be
obtained.

The differences in slope shown in Table 1 are illustrated qualitatively in
Fig. 7. The crossing point of all lines but the trapped+kinetic misalignment
trajectory corresponds to $m_{\rm crit}$. The figure shows intuitively how the
experimental parameter space is populated differently in the trapped
mechanisms (blue tones) compared to that of the QCD axion ( orange). For a
given $m_{a}$ to which an experiment is sensitive, the value of the signal
strength expected $g_{aXX}\sim 1/f_{a}$ is generically larger in the trapped
regimes, and thus the experimental reach is higher.

Figure 7: Sketch of the parametric dependence of $g_{aXX}(m_{a})$ for a fixed
axionic relic density, for the different types of misalignment discussed. For
the QCD-like and the pure trapped cases, the DIGA-inspired power-law
temperature suppression is assumed, i.e. $\alpha=8$. The dashed line reminds
that the pure trapped $Z_{\cal N}$ scenario can only take place for
$m_{a}>m_{\rm crit}$.

##  4 Implications for axion dark matter searches

In this section we discuss the implications of the above $Z_{\cal N}$ scenario
for axion DM searches, namely those experiments that rely on the hypothesis
that the axion sizeably contributes to the DM relic density. More
specifically, the $Z_{\cal N}$ axion coupling to photons, nucleons, electrons
and to the nEDM will be discussed. The whole mass range from the canonical QCD
axion mass down to the ultra-light QCD axion regime (with masses $m_{a}\ll
10^{-10}$ eV) is considered, including the fuzzy DM regime down to $m_{a}\sim
10^{-22}$ eV.

All the couplings considered but the axion-nEDM coupling are model dependent:
they can be enhanced or suppressed in specific UV completions of the axion
paradigm. We will explore for each of those the parameter space coupling vs.
$m_{a}$. On the right-hand side of the figures, though, the naive expectations
for $f_{a}$ when assuming $\mathcal{O}(1)$ couplings will be indicated as
well. This is done for pure illustrative purposes as the relation between
$f_{a}$ and axion couplings is model-dependent.

In contrast, the value of the axion-nEDM coupling only depends on $f_{a}$,
that is, it only assumes that the axion solves the strong CP problem. The same
model independence holds for previous analyses of highly-dense stellar objects
and gravitational wave prospects, which in addition did not need to assume an
axionic nature of DM. These efforts led to strong constraints on the
$\\{m_{a},f_{a}\\}$ parameter space [70, 71, 32], which can be thus directly
compared with the prospects for axion-nEDM searches.

The figures below depict with solid (translucent) colors the experimentally
excluded areas of parameter space (projected sensitivities). The blue tones
are reserved exclusively for experiments which do rely on the assumption that
axions account sizeably for DM, while the remaining colors indicate searches
which are independent of the nature of DM. In case the axion density
$\rho_{a}$ provides only a fraction of the total DM relic density $\rho_{\rm
DM}$, the sensitivity to couplings of axion DM experiments should be rescaled
by $(\rho_{a}/\rho_{\rm DM})^{1/2}$.

Crucially, we have shown in Sect. 3 that in some regions of the
$\\{m_{a},f_{a}\\}$ plane the $Z_{\cal N}$ axion can realize DM via the
misalignment mechanism and variants thereof, depending on the possible
cosmological histories of the axion field evolution in the early Universe. The
identified regions will be superimposed on the areas of parameter space
experimentally constrained/projected by axion DM experiments.

In order to compare the experimental panorama with the predictions of a
benchmark axion model, the figures will also show the expectations for the
coupling values within the $Z_{\cal N}$-KSVZ axion model developed in Ref.
[32]. The canonical KSVZ QCD axion solution (i.e. ${\cal N}=1$) is shown as a
thick yellow line, embedded into a faded band encompassing the model
dependency of the KSVZ axion [38, 40]. Oblique orange lines will signal
instead the center of the displaced yellow band that corresponds to solving
the strong CP problem with other values of ${\cal N}$, that is, for a $Z_{\cal
N}$ reduced-mass axion, see Eq. (1.4).

The entire DM relic density can be accounted for within the $Z_{\cal N}$ axion
paradigm in the regions encompassed by the purple band in the
figures.181818The band’s width accounts qualitatively for corrections to the
analytic solutions obtained in the previous section, see e.g. comment after
Eq. (3.35) and Fig. 6. For instance, in the trapped regime a factor of 2
uncertainty on $f_{a}$ has been applied.  These correspond to initial values
of $\theta_{a}$ (from the misalignment mechanism) which range from
$\theta_{1}=3/{\cal N}$ down to $\theta_{1}=0.003/{\cal N}$. The figures
illustrate that in the pure trapped regime the relic density is independent of
the initial misalignment angle. In contrast, for the simple ALP and the
trapped+kinetic mechanisms it does depend on the value of $\theta_{1}$ (in the
latter case, through its dependence on the axion velocity at $T_{\rm QCD}$).

###  4.1 Axion coupling to photons

The effective axion-photon-photon coupling $g_{a\gamma}$ is defined via the
Lagrangian

$\delta\mathcal{L}\equiv\frac{1}{4}g_{a\gamma}aF\tilde{F}\,,$ (4.1)

with [78]

$g_{a\gamma}=\frac{\alpha}{2\pi f_{a}}(E/N-1.92(4))\,,$ (4.2)

where $E/N$ takes model-dependent values ($E$ and $N$ denote the model-
dependent anomalous electromagnetic and strong contributions). This coupling
is being explored by a plethora of experiments, as illustrated in Fig. 8. For
reference, predictions of the benchmark $Z_{\cal N}$-KSVZ axion model [32] are
depicted for $E/N=0$. The conclusion is that, for a large range of values of
${\cal N}$, an axion-photon signal can be expected in large portions of the
parameter space of upcoming axion DM experiments, if the $Z_{\cal N}$ axion
accounts for the entire DM relic density.

Figure 8: Axion-photon coupling vs. axion mass. Axion limits readapted from
[79] include: laboratory axion experiments and helioscopes ( dark green) [80,
81, 82, 83], axion DM experiments ( blue) [84, 85, 86, 87, 88, 89, 90, 91, 92,
93] and astrophysical bounds ( green) [94, 95, 96, 97, 98]. Projected
sensitivities appear in translucent colors. The  orange oblique lines
represent the theoretical prediction for solutions to the strong CP problem
within the $Z_{\cal N}$ paradigm, for the benchmark axion-photon couplings
$E/N=0$ and for different (odd) numbers of SM copies ${\cal N}$. In the
purple band area the $Z_{\cal N}$ axion can account for the entire DM density,
in addition to solve the strong CP problem.

###  4.2 Axion coupling to nucleons

The axion coupling to nucleons $N$ ($N=p,n$),

$g_{aN}=C_{aN}m_{N}/f_{a}\,,$ (4.3)

is defined via the Lagrangian term

$\delta\mathcal{L}\equiv
C_{aN}\frac{\partial_{\mu}a}{2f_{a}}\overline{N}\gamma^{\mu}\gamma_{5}N\,,$
(4.4)

with [78]

$\displaystyle C_{ap}$
$\displaystyle=-0.47(3)+0.88(3)\,c^{0}_{u}-0.39(2)\,c^{0}_{d}-C_{a,{sea}}\,,$
(4.5) $\displaystyle C_{an}$
$\displaystyle=-0.02(3)+0.88(3)\,c^{0}_{d}-0.39(2)\,c^{0}_{u}-C_{a,{sea}}\,,$
(4.6)

where
$C_{a,{sea}}=0.038(5)\,c^{0}_{s}+0.012(5)\,c^{0}_{c}+0.009(2)\,c^{0}_{b}+0.0035(4)\,c^{0}_{t}$
denotes a “sea-quark” contribution in the nucleon, while in the reference
$Z_{\cal N}$-KSVZ axion model $c^{0}_{q}=0$ for all quark flavours $q$. The
parameter space of the axion-nucleon coupling vs. axion mass is displayed in
Fig. 9. It demonstrates that no signals can be expected in the axion DM
experiments foreseen up to now.

Figure 9: Axion-neutron coupling vs. axion mass. Axion limits readapted from
[79] include: axion DM experiments ( blue) [99, 100, 101, 102, 103] and
astrophysical bounds ( green) [104, 105]. Projected sensitivities appear in
translucent colors, delimited by dashed lines. The  orange oblique lines
represent the theoretical prediction for the $Z_{\cal N}$ axion-neutron
couplings, for different (odd) numbers of SM copies ${\cal N}$. The purple
band encompasses the region where the $Z_{\cal N}$ axion can account for the
entire DM density. The two theoretical predictions are represented for the
benchmark value $c^{0}_{q}=0$.

###  4.3 Axion coupling to electrons

The axion coupling to electrons,

$g_{ae}=C_{ae}m_{e}/f_{a}\,,$ (4.7)

is defined via the Lagrangian term

$\delta\mathcal{L}\equiv
C_{ae}\,\frac{\partial_{\mu}a}{2f_{a}}\,\overline{e}\gamma^{\mu}\gamma_{5}e\,,$
(4.8)

with [106, 61]

$C_{ae}=c_{e}^{0}+\frac{3\alpha^{2}}{4\pi^{2}}\left[\frac{E}{N}\log\left(\frac{f_{a}}{m_{e}}\right)-1.92(4)\log\left(\frac{{\rm
GeV}}{m_{e}}\right)\right]\,,$ (4.9)

where in the reference $Z_{\cal N}$-KSVZ axion model $c^{0}_{e}=0$ and
$E/N=0$. The parameter space of the axion-electron coupling vs. axion mass is
displayed in Fig. 10. It shows that axion-magnon conversion techniques [107,
108] could barely detect a $Z_{\cal N}$ reduced-mass axion which would account
for DM. It would be interesting to explore whether experiments which probe the
axion electron coupling with different techniques, such as QUAX [109], will be
able to test this scenario in the future.

Figure 10: Axion-electron coupling vs. axion mass. Axion limits readapted from
[79] include: solar axion experiments ( dark green) [110, 111], axion DM
experiments ( blue) [107, 108] and astrophysical bounds ( green) [112] and
hints (grey band) [113]. Projected sensitivities appear in translucent colors,
delimited by dashed lines. The  orange oblique lines represent the theoretical
prediction for the $Z_{\cal N}$ axion-electron couplings assuming
$c^{0}_{e}=0$ and $E/N=0$, for different (odd) numbers of SM copies ${\cal
N}$. The purple band encompasses the region where the $Z_{\cal N}$ axion can
account for the entire DM density.

###  4.4 Axion coupling to nEDM

The axion coupling to the nEDM operator provides a crucial test of the
$Z_{\cal N}$ axion setup: it only depends on $f_{a}$, that is, on the
assumption that the axion solves the strong CP problem, see Eq. (4.11). Unlike
the couplings discussed above, it does not depend on the details of the UV
completion of the axion model. The signal thus depends solely on $f_{a}$ and
$m_{a}$.

In fact, standard mechanisms to enhance QCD axion couplings via large Wilson
coefficients do not modify the size of the nEDM coupling, which is fixed
purely by the QCD $m_{a}$-$f_{a}$ relation. In contrast, the $Z_{\cal N}$
axion setup under discussion is responsible for a “universal” enhancement of
all axion couplings (including that to the nEDM operator), for any given
$m_{a}$ mass, following the rescaling with ${\cal N}$ of the modified
$m_{a}$-$f_{a}$ relation in Eq. (1.4).

Were the relic DM density to be entirely made up of axions, an oscillating
nEDM signal $d_{n}(t)$ could be at reach [76],

$d_{n}(t)=g_{d}\,\frac{\sqrt{2\rho_{\text{DM,
local}}}}{m_{a}}\cos(m_{a}\,t)\,,$ (4.10)

where $g_{d}=C_{an\gamma}/(m_{n}f_{a})$ is the coupling of the axion to the
nEDM operator, defined via the Lagrangian term

$\delta\mathcal{L}\equiv-\frac{i}{2}\frac{C_{an\gamma}}{m_{n}}\frac{a}{f_{a}}\overline{n}\sigma_{\mu\nu}\gamma_{5}nF^{\mu\nu}\,,$
(4.11)

where $C_{an\gamma}=0.011(5)\,e$ [114] and $\rho_{\rm DM,\,local}\approx 0.4\
\text{GeV}/\text{cm}^{3}$ is the local energy density of axion DM.191919Note
that the local relic density differs from the _mean_ relic density as measured
through the CMB and used in Eq. 3.33. For instance, the axion DM experiment
CASPEr-Electric [35, 99, 115] employs nuclear magnetic resonance techniques to
search for an oscillating nEDM. Its reach in the $\\{m_{a},f_{a}\\}$ parameter
space is displayed in Fig. 11. The reach of other techniques to probe an
oscillating nEDM based on storage rings [116] are also illustrated there,
together with present constraints which exclude ${\cal N}\gtrsim 65$.

Figure 11: Model-independent constraints in the $\\{m_{a},f_{a}\\}$ plane,
including axion DM experiments ( blue) [76, 116, 99] astrophysical ( green)
and finite-density constraints ( orange) [70, 71, 32]. The nEDM bounds have
been relaxed by a factor of 3 to take into account the stochastic nature of
the axion field, see Ref. [117]. Projected sensitivities appear in translucent
colors (those relative to finite-density effects can be seen in Fig. 8 of
[32]). The purple band encompasses the region where the $Z_{\cal N}$ axion can
account for the entire DM density.

Fig. 11 demonstrates the brilliant prospects of CASPEr-Electric Phase I and II
to discover a $Z_{\cal N}$ axion in a large range of light masses and for
$3\leq{\cal N}\lesssim 65$ if such axion accounts for DM. In particular, the
$Z_{\cal N}$ axion paradigm with reduced mass is –to our knowledge– the first
axion model that could explain a positive signal in CASPEr-Electric Phase I
(and large regions of Phase II), and at the same time solve the strong strong
CP problem. Canonical QCD axion models that exhibit large enough nEDM
couplings to be detectable at these experiments predict automatically too
heavy axions and therefore out of their reach. The figure illustrates in
addition that future proton storage ring facilities may access the region of
interest, albeit only in a small parameter region with large ${\cal N}$ and
$\theta_{1}$ values.

Furthermore, data on highly dense stellar objects share with the oscillating
nEDM signal its model-independence character: they directly provide
constraints on the $\\{m_{a},f_{a}\\}$ parameter space. Those data have the
added value of not relying on an axionic nature for DM. They have already set
strong constraints on the $\\{m_{a},f_{a}\\}$ parameter space [70, 71, 32].
For the $Z_{\cal N}$ axion scenario under discussion, they implied a strict
limit on the number of mirror worlds: $3\leq{\cal N}\lesssim 47$ for
$f_{a}\lesssim 2.4\times 10^{15}$ GeV. They also provide tantalizing discovery
prospects for any value of $f_{a}$ and down to ${\cal N}\sim 9$ –possibly even
lower– with future neutron star and gravitational wave data, down to the
ultra-light mass region, see Fig. 8 in Ref. [32]. It is thus of crucial
interest to compare them with the oscillating nEDM projects described above.
Fig. 11 tends to this task and shows that present finite-density data leave
wide open the complete parameter space that CASPEr-Electric I, II and III can
cover, with exciting synergy prospects between both type of approaches.

###  4.5 Ultra-light QCD axion dark matter

The very light pseudoscalar DM candidates usually referred to as ultra-light
axions, with mass in the range $m_{a}\in\left[10^{-22},10^{-10}\right]$ eV, do
not customarily address the strong CP problem. This is because, for the
canonical QCD axion, such light masses would imply trans-Planckian decay
constants, $f_{a}>M_{\rm Pl}$. In contrast, the $Z_{\cal N}$ axion scenario is
–to our knowledge– _the first axion model of fuzzy DM that can also solve the
strong CP problem_.

Ultra-light axions are typically searched via gravitational interactions, with
no reference to their highly suppressed couplings to the SM (for a recent
White Paper see [118]). The $Z_{\cal N}$ scenario offers instead an
interesting complementarity between strong CP related experiments and purely
gravitational probes of fuzzy dark matter.

Whenever the ultra-light DM axions obey cosinus-like potentials, they are also
subject to CMB constraints on the matter-power spectrum, due to the presence
of anharmonic terms which impact the growth of perturbations (see e.g. Refs.
[119, 120]). This implies a bound around the recombination era on any
additional energy density fluctuations beyond cold DM,
$\delta\rho/\rho\lesssim 10^{-3}$, corresponding to $a_{0}(z_{\rm
rec})/f_{a}\lesssim 10^{-3}$ [121], with $a_{0}(z_{\rm rec})$ denoting the
axion field amplitude at the recombination redshift $z_{\rm rec}$. Expanding
the axion potential around the minimum,
$V(a)\simeq\frac{1}{2}m_{a}^{2}a^{2}-\frac{1}{4!}\lambda_{a}a^{4}+\ldots\,$,
the bound translates into a condition on the quartic coupling [121]

$\left.\frac{\lambda_{a}a_{0}^{4}}{m_{a}^{2}a_{0}^{2}}\right|_{\rm
eq}\sim\frac{\lambda_{a}\,\text{eV}^{4}}{m_{a}^{4}}\lesssim 10^{-3}\,.$ (4.12)

Ref. [121] considered such constraint in the context of the $Z_{\cal N}$ axion
model of Ref. [31]. We here revisit the latter analysis with the formulae
derived in Ref. [32] for the $Z_{\cal N}$ scenario under discussion. For the
latter, the potential in the large ${\cal N}$ limit –Eqs. 1.3 and 1.4–
corresponds to a quartic coupling given by

$\lambda_{a}=m_{a}^{2}f_{a}^{2}/{\cal N}^{2}\,,$ (4.13)

and the bound in Eq. (4.12) translates into ${\cal N}\lesssim 85$. This
constraint does not impact the DM prospects discussed in this paper, as ${\cal
N}\gtrsim 65$ values are already excluded by either finite density constraint
or nEDM data, see Fig. 11.

In summary, the $Z_{\cal N}$ scenario under discussion allows us to account
for DM and solve the strong CP problem in a sizeable fraction of the ultra-
light axions parameter space: $m_{a}\in\left[10^{-22},10^{-10}\right]$ eV. For
instance, Fig. 8 shows that a photon-axion signal could be at reach for the
upper masses in that range, provided $\theta_{1}$ takes sizeable values and
the kinetic misalignment regime takes place. More importantly, in that entire
mass range a model-independent discovery signal is open for observation at the
oscillating nEDM experiments such as CASPEr-Electric, see Fig. 11. The latter
figure also points to the exciting prospects ahead, because the near-future
data from high-density stellar systems and gravitational wave facilities
should cover that entire mass region, for any value of $\theta_{1}$ and down
to ${\cal N}\sim 9$.

##  5 Conclusions

This work constitutes a proof-of-concept that an axion lighter –or even much
lighter– than the canonical QCD one may both solve the SM strong CP problem
and account for the entire DM density of the Universe. Large regions of the
$\\{m_{a},\,f_{a}\\}$ parameter space to the left of the of the canonical QCD
axion band can accomplish that goal.

While the implications of a $Z_{\cal N}$ shift symmetry to solve the strong CP
problem (with a $1/{\cal N}$ probability and ${\cal N}$ degenerate worlds)
have been previously analyzed [32], leading to a lighter-than-usual axion, the
question of DM and the cosmological evolution was left unexplored. We showed
here that the evolution of the axion field through the cosmological history
departs drastically from both the standard one and from previously considered
mirror world scenarios.

In particular, we identified a novel axion production mechanism which holds
whenever $f_{a}\lesssim 3.2\times 10^{17}\,\text{GeV}$: trapped misalignment,
which is a direct consequence of the temperature dependence of the axion
potential. Two distinct stages of oscillations take place. At large
temperatures the minimum of the finite-temperature potential shifts from its
vacuum value, i.e. $\theta=0$, to large values, e.g. $\theta=\pi$, where the
axion field gets trapped down to a temperature $T\sim T_{\rm QCD}$. The axion
mass is unsuppressed during this trapped period and thus of the order of the
canonical QCD axion mass. The underlying reason is that the SM thermal bath
explicitly breaks the $Z_{\cal N}$ symmetry, because its temperature must be
higher than that of the other mirror worlds. This trapped period has a major
cosmological impact: the subsequent onset of oscillations around the true
minimum at $\theta=0$ is delayed as compared with the standard QCD axion
scenario. The result is an important enhancement of the DM relic density. In
other words, lower $f_{a}$ values can now account for DM.

We have determined the minimum kinetic energy $K_{\text{min}}$ required at the
end of trapping for the axion to roll over $\sim{\cal N}/2$ maxima before it
starts to oscillate around the true minimum (so as to solve the strong CP
problem). We showed that the axion kinetic energy is of
$\cal{O}(K_{\text{min}})$ in sizeable regions of the parameter space, fuelled
by the (much larger than in vacuum) axion mass at the end of the trapped
period. In this pure trapped scenario, the final oscillations start at
temperatures smaller but close to $T\sim T_{\rm QCD}$.

In fact, the axion kinetic energy at the end of trapping is shown to be in
general much larger than $K_{\text{min}}$. Trapped misalignment then
automatically seeds kinetic misalignment [56] between $T\sim T_{\rm QCD}$ and
lower temperatures. The axion rolls for a long time over the low-temperature
potential barriers before final oscillations start at $T\ll T_{\rm QCD}$,
extending further the delay of oscillations around the true minimum ensured by
the trapped period. In consequence, the trapped+kinetic misalignment mechanism
enhances even more strongly the DM relic density.

Our novel trapped mechanism is more general than the $Z_{\cal N}$ framework
considered here. It could arise in a large variety of ALP or QCD axion
scenarios. For instance, it may apply to axion theories in which an explicit
source of PQ breaking is active only at high temperatures and the transition
to the true vacuum is non-adiabatic. Note also that in our scenario kinetic
misalignment does not rely on the presence of non-renormalizable PQ-breaking
operators required in the original formulation [56]. It is instead directly
seeded by trapped misalignment, which is itself a pure temperature effect.

For values of the $Z_{\cal N}$ axion scale $f_{a}\gtrsim 3.2\times
10^{17}\,\text{GeV}$, the trapped mechanism does not take place, since there
is only one stage of oscillations. The $T=0$ potential is already developed
when the Hubble friction is overcome, and the axion oscillates from the start
around the true minimum $\theta_{a}=0$. The relic density corresponds then to
that of a simple ALP regime with constant axion mass, alike to the standard
QCD axion scenario.

We have determined the current axion relic density stemming from the various
misalignment mechanisms, analyzing their different dependence on the
$\\{m_{a},\,f_{a}\,,{\cal N}\\}$ variables. The ultimate dependence on the
arbitrary initial misalignment angle has been determined as well for the
simple ALP and trapped+kinetic scenarios. For the pure trapped scenario, the
relic density turns out to be independent of the initial misalignment, which
results in a band centered around ${\cal N}\sim 21$ to account for the
ensemble of DM. Overall, DM solutions are found within the $Z_{\cal N}$
paradigm for any value of $3\leq{\cal N}\lesssim 65$.

The results above have been next confronted with the experimental arena of the
so-called axion DM searches. As a wonderful byproduct of the lower-than-usual
$f_{a}$ values allowed in the $Z_{\cal N}$ axion paradigm to solve the strong
CP problem, all axion-SM couplings are equally enhanced for a given $m_{a}$.
This increases the testability of the theory in current and future
experiments. In consequence, many axion DM experiments which up to now only
aimed to target the nature of DM, are simultaneously addressing the SM strong
CP problem, provided mirror worlds exist. We have studied the present and
projected experimental sensitivity to the axion coupling to photons, electrons
and nucleons, as a function of the axion mass and ${\cal N}$. It follows that
an axion-photon signal is at reach in large portions of the parameter space of
upcoming axion DM experiments, while no such prospects result for the coupling
to nucleons, and only marginally for the coupling to electrons.

A different and crucial test is provided by the $aG\tilde{G}$ coupling (that
fixes the value of $1/f_{a}$), which can be entirely traded by an axion-nEDM
coupling. The signal has two remarkable properties, for any given $m_{a}$: i)
in all generality, it does not depend on the details of the putative UV
completion of the axion model, unlike all other couplings considered; ii) its
strength is enhanced in the $Z_{\cal N}$ paradigm, which is impossible in any
model of the canonical QCD axion. It follows that the $Z_{\cal N}$ paradigm is
–to our knowledge– the only true axion theory that could explain a positive
signal in CASPEr-Electric phase I and in a large region of the parameter space
in phase II. The reason is that a traditional QCD axion with an nEDM coupling
in the range to be probed by that experiment would be automatically heavier,
and therefore outside its reach. Such a signal could instead account for DM
and solve the strong CP problem within the $Z_{\cal N}$ scenario. The same
applies to the Storage Ring projects that aim to detect oscillating EDMs.

Furthermore, our results demonstrate a beautiful synergy and complementarity
between the expected reach of axion DM experiments and axion experiments which
are independent of the nature of DM. For instance, oscillating nEDM
experiments on one side, and data expected from highly dense stellar objects
and gravitational waves on the other, have a wide overlap in their sensitivity
reach. Their combination will cover in the next decades the full range of
possible ${\cal N}$ and $m_{a}$ values, in the mass range from the standard
QCD axion mass down to $\sim 10^{-22}$ eV, that is, down to the fuzzy DM
range. To our knowledge, the $Z_{\cal N}$ axion discussed here is the first
model of fuzzy DM which also solves the strong CP problem.

## Acknowledgments

We thank Gonzalo Alonso-Álvarez, Gary Centers, Victor Enguita, Yann
Gouttenoire, Benjamin Grinstein, Lam Hui, David B. Kaplan, Ryosuke Sato,
Geraldine Servant, Philip Sørensen, Luca Visinelli and Neal Weiner for
illuminating discussions. The work of L.D.L. is supported by the Marie
Skłodowska-Curie Individual Fellowship grant AXIONRUSH (GA 840791). L.D.L.,
P.Q. and A.R. acknowledge support by the Deutsche Forschungsgemeinschaft under
Germany’s Excellence Strategy \- EXC 2121 Quantum Universe - 390833306. M.B.G.
acknowledges support from the “Spanish Agencia Estatal de Investigación” (AEI)
and the EU “Fondo Europeo de Desarrollo Regional” (FEDER) through the projects
FPA2016-78645-P and PID2019-108892RB-I00/AEI/10.13039/501100011033. M.B.G. and
P. Q. acknowledge support from the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreements 690575
(RISE InvisiblesPlus) and 674896 (ITN ELUSIVES), as well as from the Spanish
Research Agency (Agencia Estatal de Investigación) through the grant IFT
Centro de Excelencia Severo Ochoa SEV-2016-0597. This project has received
funding/support from the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie grant agreement No 860881-HIDDeN.

##  Appendix A Ratios of $g_{*}$ and $g_{s}$

Throughout the paper, different ratios of the number of degrees of freedom in
the SM and the mirror worlds appear. They need to be computed numerically and
correspond to $\mathcal{O}(1)$ corrections that in large regions of the
parameter space can be neglected. This allows the derivation of simple
analytical formulas for the different observables, such as the contributions
to $N_{\rm eff}$, or the axion relic density. In this Appendix we review the
definition of these quantities and compute them numerically as a function of
both the SM and mirrork temperatures for different values of $\gamma_{k}$. To
this aim, the tabulated effective number of degrees of freedom (for the SM
case) as a function of temperature provided in Ref. [122] will be used (see
also Fig. 12).

Figure 12: Effective number of degrees of freedom.

#### Ratio of temperatures vs. $\gamma_{k}$: $b_{k}(T,T_{k})$

The function $b_{k}(T,T_{k})$ defined in Eq. (2.3),

$\displaystyle
b_{k}(T,T_{k})\equiv\bigg{[}\frac{g_{s}(T)}{g_{s}\left(T_{k}\right)}\bigg{]}^{1/3}=\frac{T_{k}/T}{\gamma_{k}}\,,$
(A.1)

is shown to approach $1$ for $T_{k}\gg T_{\text{QCD}}$ in the right plot in
Fig. 13. Therefore the approximation $b_{k}(T,T_{k})\sim 1$ can be safely
adopted in the computation of the high-temperature axion potential in Eq. 3.7.

Figure 13: Numerical evaluation of $b_{k}(T,T_{k})$ as a function of the SM
temperature $T$ (left) and of the k-th mirror world temperature $T_{k}$
(right), for different values of $\gamma_{k}$.

#### Bounds on $\gamma_{k}$ through $N_{\rm eff}$: $c_{k}(T,T_{k})$

The ratio

$\displaystyle
c_{k}\left(T,T_{k}\right)\equiv\left[\frac{g_{*}\left(T_{k}\right)}{g_{*}(T)}\right]\left[\frac{g_{s}(T)}{g_{s}\left(T_{k}\right)}\right]^{4/3}\,,$
(A.2)

which appears in Eq. (2.5), has been depicted in Fig. 14. The left plot shows
that its value is $c_{k}(T=1\text{ MeV},T_{k})\lesssim 1.2$ whenever the SM
BBN is taking place, i.e. $T\sim 1\text{ MeV}$. In consequence, the
approximation $c_{k}(T,T_{k})\sim 1$ can be safely applied when using the
experimental limits in Eq. 2.1 to extract bounds on $\gamma_{k}$ from Eq. 2.6.

Figure 14: Numerical evaluation of $c_{k}(T,T_{k})$ as a function of the SM
temperature $T$ (left) and of the k-th mirror world temperature $T_{k}$
(right), for different values of $\gamma_{k}$.

#### Trapped+kinetic misalignment: $\mathcal{F}_{\rm
kin,\,2}(T_{1},T_{1}^{\prime})$

The prediction for the relic density in the kinetic+trapped misalignment
scenario in Eq. 3.42 involves the ratio

$\displaystyle\mathcal{F}_{\rm
kin,1}(T_{1})\equiv\Bigg{(}\frac{\sqrt{g_{*}(T_{1})g_{*}(T_{\rm
QCD})}}{3.38}\Bigg{)}^{\frac{3}{4}}\left(\frac{3.93}{\sqrt{g_{s}(T_{1})g_{s}(T_{\rm
QCD})}}\right)\,,$ (A.3)

which feeds into the final result in Eq. 3.43 (see also Appendix C) through
$\mathcal{F}_{\rm kin,\,2}$:

$\mathcal{F}_{\rm kin,\,2}(T_{1},T_{1}^{\prime})\equiv\mathcal{F}_{\rm
kin,1}(T_{1})\times b^{\prime}(T_{1},T^{\prime}_{1})^{1/3}\times
c^{\prime}(T_{\rm BBN},T^{\prime}_{\rm
BBN})^{-1/12}\times\left(3.38/g_{*}(T_{1})\right)^{1/12}=\\\
\left[\frac{\sqrt{g_{*}(T_{1})g_{*}(T_{\text{QCD}})}}{3.38}\right]^{\frac{3}{4}}\left[\frac{3.93}{\sqrt{g_{s}(T_{1})g_{s}(T_{\text{QCD}})}}\right]\left[\frac{g_{s}(T_{1})}{g_{s}\left(T^{\prime}_{1}\right)}\right]^{1/9}\left[\frac{g_{*}\left(T_{\rm
BBN}\right)}{g_{*}(T^{\prime}_{\rm BBN})}\right]^{1/12}\\\
\times\left[\frac{g_{s}(T^{\prime}_{\rm BBN})}{g_{s}\left(T_{\rm
BBN}\right)}\right]^{1/9}\left(\frac{3.38}{g_{*}(T_{1})}\right)^{1/12}\,,$
(A.4)

where $T^{\prime}_{\rm BBN}$ denotes the temperature of all the mirror copies
of the SM when the temperature of the latter is $T_{\rm BBN}\simeq 1\text{
MeV}$.

The numerical evaluation of $\mathcal{F}_{\rm kin,\,2}$ can be found in Fig.
15. It shows that its values lie in the range $0.45-0.85$, depending on the
temperature at the onset of the first stage of oscillations, $T_{1}$, and on
the value of $\gamma_{k}$.

Figure 15: Numerical evaluation of $\mathcal{F}_{\rm
kin,\,2}(T_{1},T_{1}^{\prime})$ as a function of the SM temperature $T$ (left)
and of the k-th mirror world temperature $T_{k}$ (right), for different values
of $\gamma_{k}$.

##  Appendix B Anharmonicity function

The relic density of a scalar field with a harmonic potential can be computed
analytically using the adiabatic approximation. However, pseudo Goldstone
bosons are in general described by periodic potentials that only resemble a
harmonic potential close to the minimum. If the initial misalignment angle is
large $\theta\gtrsim\pi/2$, the onset of oscillations gets delayed leading to
an enhancement of the relic density. This enhancement can be parametrized via
the so-called anharmonicity factor $f_{\rm anh}$. For an ALP-like regime, one
can use an empirical formula proposed in Ref. [57] which encodes this
deviation from the harmonic result:

$\displaystyle\rho_{a,\,0}=\frac{1}{2}m_{a}^{2}\left(\frac{a_{1}}{a_{\rm
0}}\right)^{3}\theta_{1}^{2}f_{a}^{2}f_{\rm anh}(\theta_{1})\,.$ (B.1)

The definition of $f_{\rm anh}(\theta_{1})$ thus depends on the criteria
chosen to define the onset of oscillations at temperature $T_{1}$. In our
case, $T_{1}$ correspond to the temperature in which the axion mass term
overcomes Hubble friction as given in Eq. (3.20).

Figure 16: Comparison of the empirical formula in Eq. (B.4) for the
anharmonicity factor with the numerical result in the simple ALP case.

For the case of a simple ALP regime with constant mass and a cosine potential,
and for large misalignment angles $\pi-\theta_{1}\lesssim 10^{-2}$, our
numerical result is in excellent agreement with the Ref. [57] proposal.
However, the latter does not apply for small misalignment angles, as it does
not have the expected limiting behaviour $f_{\rm
anh}\left(\theta_{1}\right)\xrightarrow{\,\,\theta_{1}\rightarrow 0\,\,}1$.
For this reason, we will employ the following anharmonicity prescription:

$\displaystyle f_{\rm
anh}(\theta_{1})=\left\\{\begin{array}[]{ll}\left[1-\log\left({1-\left({\theta_{1}}/{\pi}\right)^{3}}\right)\right]^{1.3}&\text{
for }\quad\pi-\theta_{1}>10^{-2}\,,\\\\[12.0pt]
0.06\left[\log\left(\frac{2.32}{\pi-|\theta_{1}|}\right)+4\,\log\left(\log\left(\frac{2.32}{\pi-|\theta_{1}|}\right)\right)\right]&\text{
for }\quad\pi-\theta_{1}<10^{-2}\,,\end{array}\right.$ (B.4)

where the first expression is inspired by Ref. [123], which obtained an
expression for the anharmonicity factor for the QCD axion (we have adapted it
so as to fit our numerical result for the constant mass ALP); and the second
expression is taken from Ref. [57].

##  Appendix C Maximal relic density in the trapped+kinetic regime

The predicted relic density in the trapped+kinetic scenario of Eq. 3.41 can be
rewritten as

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm DM}}\bigg{|}_{\rm tr+kin}\simeq
7.9\sqrt{\frac{m_{a}}{\mathrm{eV}}}\bigg{(}\frac{m_{a}^{2}\,m_{a,\pi}^{\rm
QCD}}{m_{\rm crit}^{3}}\bigg{)}^{1/4}\left(\frac{m_{a,\pi}^{\rm
QCD}}{m_{1}}\right)^{1/4}\left(\frac{\,f_{a}}{10^{12}\,\mathrm{GeV}}\right)^{2}\,\frac{\left|\theta_{1}\right|}{{\cal
N}}\,\mathcal{F}_{\rm kin,1}(T_{1})\,,$ (C.1)

where $\mathcal{F}_{\rm kin,1}(T_{1})$ was given in Eq. (A.3). The relic
density depends via $m_{1}$ on the full temperature-dependent potential in Eq.
3.3 and therefore on the temperatures of all the different worlds (i.e. $T$
and the values of $\gamma_{k}$, see Eq. 2.3). In the simplified case in which
all the mirror copies of the SM have the same temperature $T_{k\neq
0}=T^{\prime}$ one can use the potential in Eq. 3.7 to express the factor
$({m_{a,\pi}^{\rm QCD}}/{m_{1}})^{1/4}$ as a function of the parameters of the
finite-temperature potential,

$\displaystyle\left(\frac{m_{a,\pi}^{\rm
QCD}}{m_{1}}\right)^{1/4}\simeq\left(\frac{\gamma^{\prime}\,b^{\prime}(T_{1},T_{1}^{\prime})\,\sqrt{m_{a,\pi}^{\rm
QCD}M_{\rm pl}}}{T_{\rm
QCD}}\right)^{\frac{\alpha}{2\alpha+8}}\left(1.67\times
1.66\,\sqrt{g_{*}(T_{1})}\right)^{\frac{-\alpha}{4\alpha+16}}\,.$ (C.2)

For $\alpha=8$ it results in a parametric dependence on the mass scales of the
scenario of the form

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm
DM}}\propto\sqrt{\frac{m_{a}}{\mathrm{eV}}}\,\Bigg{[}\frac{m_{a}^{6}\,\big{(}m_{a,\pi}^{\rm
QCD}\big{)}^{5}M_{\rm Pl}^{2}}{m_{\rm crit}^{9}T_{\rm
QCD}^{4}}\Bigg{]}^{1/12}\left(\frac{\,f_{a}}{10^{12}\,\mathrm{GeV}}\right)^{2}\,\frac{\left|\theta_{1}\right|}{{\cal
N}}\,,$ (C.3)

from which the parametric dependences shown in Table 1 follow.

Finally, the maximal relic density that can be generated by the
trapped+kinetic misalignment mechanism can be obtained as a function of
$f_{a}$ and ${\cal N}$, by choosing the value of $\gamma^{\prime}$ that
saturates the $N_{\rm eff}$ bound in Eq. 2.6.202020Corresponding to
$\gamma^{\prime}={0.51}/{({\cal N}-1)^{1/4}}c^{\prime-1/4}(T_{\rm
BBN},T^{\prime}_{\rm BBN})$, with $c^{\prime}\equiv c_{k\neq 0}$.  In the case
of the $Z_{\cal N}$ axion mass in Eq. 1.4 we obtain

$\displaystyle\frac{\rho_{a,0}}{\rho_{\rm DM}}\bigg{|}_{\rm
tr+kin,\,max}\simeq
0.64\left(\frac{\,f_{a}}{10^{9}\,\mathrm{GeV}}\right)^{5/6}\,\frac{z^{{\cal
N}/2}}{{\cal N}^{1/3}}\,\left|\theta_{1}\right|\,\mathcal{F}_{\rm
kin,\,2}(T_{1})\,,$ (C.4)

where $\mathcal{F}_{\rm kin,\,2}(T_{1})$ can be found in Eq. (A.4).

## References

  * [1] R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [2] R. Peccei and H. R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D 16 (1977) 1791–1797.
  * [3] S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978) 223–226.
  * [4] F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * [5] J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B 120 (1983) 127–132.
  * [6] L. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B 120 (1983) 133–136.
  * [7] M. Dine and W. Fischler, “The Not So Harmless Axion,” Phys. Lett. B 120 (1983) 137–141.
  * [8] J. E. Kim, “Weak Interaction Singlet and Strong CP Invariance,” Phys. Rev. Lett. 43 (1979) 103.
  * [9] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, “Can Confinement Ensure Natural CP Invariance of Strong Interactions?,” Nucl. Phys. B166 (1980) 493–506.
  * [10] A. R. Zhitnitsky, “On Possible Suppression of the Axion Hadron Interactions. (In Russian),” Sov. J. Nucl. Phys. 31 (1980) 260. [Yad. Fiz.31,497(1980)].
  * [11] M. Dine, W. Fischler, and M. Srednicki, “A Simple Solution to the Strong CP Problem with a Harmless Axion,” Phys. Lett. B104 (1981) 199–202.
  * [12] J. E. Kim, “A Composite Invisible Axion,” Phys. Rev. D31 (1985) 1733.
  * [13] K. Choi and J. E. Kim, “Dynamical Axion,” Phys. Rev. D32 (1985) 1828.
  * [14] V. A. Rubakov, “Grand unification and heavy axion,” JETP Lett. 65 (1997) 621–624, arXiv:hep-ph/9703409 [hep-ph].
  * [15] Z. Berezhiani, L. Gianfagna, and M. Giannotti, “Strong CP problem and mirror world: The Weinberg-Wilczek axion revisited,” Phys. Lett. B500 (2001) 286–296, arXiv:hep-ph/0009290 [hep-ph].
  * [16] L. Gianfagna, M. Giannotti, and F. Nesti, “Mirror world, supersymmetric axion and gamma ray bursts,” JHEP 10 (2004) 044, arXiv:hep-ph/0409185.
  * [17] S. D. H. Hsu and F. Sannino, “New solutions to the strong CP problem,” Phys. Lett. B605 (2005) 369–375, arXiv:hep-ph/0408319 [hep-ph].
  * [18] A. Hook, “Anomalous solutions to the strong CP problem,” Phys. Rev. Lett. 114 no. 14, (2015) 141801, arXiv:1411.3325 [hep-ph].
  * [19] H. Fukuda, K. Harigaya, M. Ibe, and T. T. Yanagida, “Model of visible QCD axion,” Phys. Rev. D92 no. 1, (2015) 015021, arXiv:1504.06084 [hep-ph].
  * [20] C.-W. Chiang, H. Fukuda, M. Ibe, and T. T. Yanagida, “750 GeV diphoton resonance in a visible heavy QCD axion model,” Phys. Rev. D93 no. 9, (2016) 095016, arXiv:1602.07909 [hep-ph].
  * [21] S. Dimopoulos, A. Hook, J. Huang, and G. Marques-Tavares, “A Collider Observable QCD Axion,” JHEP 11 (2016) 052, arXiv:1606.03097 [hep-ph].
  * [22] T. Gherghetta, N. Nagata, and M. Shifman, “A Visible QCD Axion from an Enlarged Color Group,” Phys. Rev. D93 no. 11, (2016) 115010, arXiv:1604.01127 [hep-ph].
  * [23] A. Kobakhidze, “Heavy axion in asymptotically safe QCD,” arXiv:1607.06552 [hep-ph].
  * [24] P. Agrawal and K. Howe, “Factoring the Strong CP Problem,” JHEP 12 (2018) 029, arXiv:1710.04213 [hep-ph].
  * [25] P. Agrawal and K. Howe, “A Flavorful Factoring of the Strong CP Problem,” JHEP 12 (2018) 035, arXiv:1712.05803 [hep-ph].
  * [26] M. K. Gaillard, M. B. Gavela, R. Houtz, P. Quilez, and R. Del Rey, “Color unified dynamical axion,” Eur. Phys. J. C78 no. 11, (2018) 972, arXiv:1805.06465 [hep-ph].
  * [27] M. A. Buen-Abad and J. Fan, “Dynamical axion misalignment with small instantons,” JHEP 12 (2019) 161, arXiv:1911.05737 [hep-ph].
  * [28] A. Hook, S. Kumar, Z. Liu, and R. Sundrum, “High Quality QCD Axion and the LHC,” Phys. Rev. Lett. 124 no. 22, (2020) 221801, arXiv:1911.12364 [hep-ph].
  * [29] C. Csáki, M. Ruhdorfer, and Y. Shirman, “UV Sensitivity of the Axion Mass from Instantons in Partially Broken Gauge Groups,” JHEP 04 (2020) 031, arXiv:1912.02197 [hep-ph].
  * [30] T. Gherghetta and M. D. Nguyen, “A Composite Higgs with a Heavy Composite Axion,” arXiv:2007.10875 [hep-ph].
  * [31] A. Hook, “Solving the Hierarchy Problem Discretely,” Phys. Rev. Lett. 120 no. 26, (2018) 261802, arXiv:1802.10093 [hep-ph].
  * [32] L. Di Luzio, B. Gavela, P. Quilez, and A. Ringwald, “An even lighter QCD axion,” IFT-UAM/CSIC-20-143 FTUAM-20-21 DESY 21-010 To appear, 2020.
  * [33] S. Das and A. Hook, “Non-linearly realized discrete symmetries,” JHEP 10 (2020) 071, arXiv:2006.10767 [hep-ph].
  * [34] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, “Ultralight scalars as cosmological dark matter,” Phys. Rev. D 95 no. 4, (2017) 043541, arXiv:1610.08297 [astro-ph.CO].
  * [35] D. Budker, P. W. Graham, M. Ledbetter, S. Rajendran, and A. Sushkov, “Proposal for a Cosmic Axion Spin Precession Experiment (CASPEr),” Phys. Rev. X4 no. 2, (2014) 021030, arXiv:1306.6089 [hep-ph].
  * [36] D. F. Jackson Kimball et al., “Overview of the Cosmic Axion Spin Precession Experiment (CASPEr),” arXiv:1711.08999 [physics.ins-det].
  * [37] D. Jackson Kimball et al., “Overview of the Cosmic Axion Spin Precession Experiment (CASPEr),” Springer Proc. Phys. 245 (2020) 105–121, arXiv:1711.08999 [physics.ins-det].
  * [38] L. Di Luzio, F. Mescia, and E. Nardi, “Redefining the Axion Window,” Phys. Rev. Lett. 118 no. 3, (2017) 031801, arXiv:1610.07593 [hep-ph].
  * [39] M. Farina, D. Pappadopulo, F. Rompineve, and A. Tesi, “The photo-philic QCD axion,” JHEP 01 (2017) 095, arXiv:1611.09855 [hep-ph].
  * [40] L. Di Luzio, F. Mescia, and E. Nardi, “Window for preferred axion models,” Phys. Rev. D 96 no. 7, (2017) 075003, arXiv:1705.05370 [hep-ph].
  * [41] P. Agrawal, J. Fan, M. Reece, and L.-T. Wang, “Experimental Targets for Photon Couplings of the QCD Axion,” JHEP 02 (2018) 006, arXiv:1709.06085 [hep-ph].
  * [42] G. Marques-Tavares and M. Teo, “Light axions with large hadronic couplings,” JHEP 05 (2018) 180, arXiv:1803.07575 [hep-ph].
  * [43] L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, “The landscape of QCD axion models,” Phys. Rept. 870 (2020) 1–117, arXiv:2003.01100 [hep-ph].
  * [44] L. Darmé, L. Di Luzio, M. Giannotti, and E. Nardi, “Selective enhancement of the QCD axion couplings,” arXiv:2010.15846 [hep-ph].
  * [45] Z. Berezhiani, “Mirror world and its cosmological consequences,” Int. J. Mod. Phys. A 19 (2004) 3775–3806, arXiv:hep-ph/0312335.
  * [46] Z. Berezhiani, D. Comelli, and F. L. Villante, “The Early mirror universe: Inflation, baryogenesis, nucleosynthesis and dark matter,” Phys. Lett. B503 (2001) 362–375, arXiv:hep-ph/0008105 [hep-ph].
  * [47] E. W. Kolb, D. Seckel, and M. S. Turner, “The Shadow World,” Nature 314 (1985) 415–419.
  * [48] V. S. Berezinsky and A. Vilenkin, “Ultrahigh-energy neutrinos from hidden sector topological defects,” Phys. Rev. D62 (2000) 083512, arXiv:hep-ph/9908257 [hep-ph].
  * [49] G. Dvali, E. Koutsangelas, and F. Kuhnel, “Compact Dark Matter Objects via $N$ Dark Sectors,” Phys. Rev. D 101 (2020) 083533, arXiv:1911.13281 [astro-ph.CO].
  * [50] M. Giannotti, “Mirror world and axion: Relaxing cosmological bounds,” Int. J. Mod. Phys. A 20 (2005) 2454–2458, arXiv:astro-ph/0504636.
  * [51] S. Nakagawa, F. Takahashi, and M. Yamada, “Trapping Effect for QCD Axion Dark Matter,” arXiv:2012.13592 [hep-ph].
  * [52] R. Alonso and J. Scholtz, “Kick-alignment: matter asymmetry sourced dark matter,” arXiv:2012.14907 [hep-ph].
  * [53] R. T. Co, E. Gonzalez, and K. Harigaya, “Axion Misalignment Driven to the Hilltop,” JHEP 05 (2019) 163, arXiv:1812.11192 [hep-ph].
  * [54] F. Takahashi and W. Yin, “QCD axion on hilltop by a phase shift of $\pi$,” JHEP 10 (2019) 120, arXiv:1908.06071 [hep-ph].
  * [55] J. Huang, A. Madden, D. Racco, and M. Reig, “Maximal axion misalignment from a minimal model,” arXiv:2006.07379 [hep-ph].
  * [56] R. T. Co, L. J. Hall, and K. Harigaya, “Kinetic Misalignment Mechanism,” arXiv:1910.14152 [hep-ph].
  * [57] A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, “Large-misalignment mechanism for the formation of compact axion structures: Signatures from the QCD axion to fuzzy dark matter,” Phys. Rev. D 101 no. 8, (2020) 083014, arXiv:1909.11665 [astro-ph.CO].
  * [58] B. D. Fields, K. A. Olive, T.-H. Yeh, and C. Young, “Big-Bang Nucleosynthesis After Planck,” JCAP 03 (2020) 010, arXiv:1912.01132 [astro-ph.CO].
  * [59] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. I. Overview and the cosmological legacy of Planck,” Astron. Astrophys. 641 (2020) A1, arXiv:1807.06205 [astro-ph.CO].
  * [60] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” arXiv:1807.06209 [astro-ph.CO].
  * [61] S. Chang and K. Choi, “Hadronic axion window and the big bang nucleosynthesis,” Phys. Lett. B 316 (1993) 51–56, arXiv:hep-ph/9306216.
  * [62] S. Hannestad, A. Mirizzi, and G. Raffelt, “New cosmological mass limit on thermal relic axions,” JCAP 07 (2005) 002, arXiv:hep-ph/0504059.
  * [63] L. Di Luzio, G. Martinelli, and G. Piazza, “Axion hot dark matter bound, reliably,” arXiv:2101.10330 [hep-ph].
  * [64] D. Baumann, “Primordial Cosmology,” PoS TASI2017 (2018) 009, arXiv:1807.03098 [hep-th].
  * [65] Z. Berezhiani, “Through the looking-glass: Alice’s adventures in mirror world,” arXiv:hep-ph/0508233 [hep-ph].
  * [66] R. R. Volkas, A. J. Davies, and G. C. Joshi, “Naturalness of the invisible axion model,” Phys. Lett. B215 (1988) 133–138.
  * [67] R. Foot, A. Kobakhidze, K. L. McDonald, and R. R. Volkas, “Poincare protection for a natural electroweak scale,” Phys. Rev. D 89 no. 11, (2014) 115018, arXiv:1310.0223 [hep-ph].
  * [68] L. Ubaldi, “Effects of theta on the deuteron binding energy and the triple-alpha process,” Phys. Rev. D81 (2010) 025011, arXiv:0811.1599 [hep-ph].
  * [69] S. Borsanyi et al., “Calculation of the axion mass based on high-temperature lattice quantum chromodynamics,” Nature 539 no. 7627, (2016) 69–71, arXiv:1606.07494 [hep-lat].
  * [70] A. Hook and J. Huang, “Probing axions with neutron star inspirals and other stellar processes,” JHEP 06 (2018) 036, arXiv:1708.08464 [hep-ph].
  * [71] J. Huang, M. C. Johnson, L. Sagunski, M. Sakellariadou, and J. Zhang, “Prospects for axion searches with Advanced LIGO through binary mergers,” Phys. Rev. D 99 no. 6, (2019) 063013, arXiv:1807.02133 [hep-ph].
  * [72] D. Brzeminski, Z. Chacko, A. Dev, and A. Hook, “A Time-Varying Fine Structure Constant from Naturally Ultralight Dark Matter,” arXiv:2012.02787 [hep-ph].
  * [73] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, “WISPy Cold Dark Matter,” JCAP 06 (2012) 013, arXiv:1201.5902 [hep-ph].
  * [74] G. Alonso-Álvarez, T. Hugle, and J. Jaeckel, “Misalignment \& Co.: (Pseudo-)scalar and vector dark matter with curvature couplings,” JCAP 02 (2020) 014, arXiv:1905.09836 [hep-ph].
  * [75] D. J. E. Marsh, “Axion Cosmology,” Phys. Rept. 643 (2016) 1–79, arXiv:1510.07633 [astro-ph.CO].
  * [76] P. W. Graham and S. Rajendran, “New Observables for Direct Detection of Axion Dark Matter,” Phys. Rev. D88 (2013) 035023, arXiv:1306.6088 [hep-ph].
  * [77] N. Fonseca, E. Morgante, R. Sato, and G. Servant, “Axion fragmentation,” JHEP 04 (2020) 010, arXiv:1911.08472 [hep-ph].
  * [78] G. Grilli di Cortona, E. Hardy, J. Pardo Vega, and G. Villadoro, “The QCD axion, precisely,” JHEP 01 (2016) 034, arXiv:1511.02867 [hep-ph].
  * [79] C. O’Hare, “cajohare/axionlimits: Axionlimits,” July, 2020.  https://doi.org/10.5281/zenodo.3932430.
  * [80] R. Bähre et al., “Any light particle search II —Technical Design Report,” JINST 8 (2013) T09001, arXiv:1302.5647 [physics.ins-det].
  * [81] E. Armengaud et al., “Conceptual Design of the International Axion Observatory (IAXO),” JINST 9 (2014) T05002, arXiv:1401.3233 [physics.ins-det].
  * [82] CAST Collaboration, V. Anastassopoulos et al., “New CAST Limit on the Axion-Photon Interaction,” Nature Phys. 13 (2017) 584–590, arXiv:1705.02290 [hep-ex].
  * [83] BabyIAXO Collaboration, A. Abeln et al., “Conceptual Design of BabyIAXO, the intermediate stage towards the International Axion Observatory,” arXiv:2010.12076 [physics.ins-det].
  * [84] C. Hagmann, P. Sikivie, N. Sullivan, and D. Tanner, “Results from a search for cosmic axions,” Phys. Rev. D 42 (1990) 1297–1300.
  * [85] 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–72, arXiv:1706.00209 [physics.ins-det].
  * [86] D. Alesini, D. Babusci, D. Di Gioacchino, C. Gatti, G. Lamanna, and C. Ligi, “The KLASH Proposal,” arXiv:1707.06010 [physics.ins-det].
  * [87] HAYSTAC Collaboration, L. Zhong et al., “Results from phase 1 of the HAYSTAC microwave cavity axion experiment,” Phys. Rev. D 97 no. 9, (2018) 092001, arXiv:1803.03690 [hep-ex].
  * [88] D. J. Marsh, K.-C. Fong, E. W. Lentz, L. Smejkal, and M. N. Ali, “Proposal to Detect Dark Matter using Axionic Topological Antiferromagnets,” Phys. Rev. Lett. 123 no. 12, (2019) 121601, arXiv:1807.08810 [hep-ph].
  * [89] ADMX Collaboration, T. Braine et al., “Extended Search for the Invisible Axion with the Axion Dark Matter Experiment,” Phys. Rev. Lett. 124 no. 10, (2020) 101303, arXiv:1910.08638 [hep-ex].
  * [90] S. Lee, S. Ahn, J. Choi, B. Ko, and Y. Semertzidis, “Axion Dark Matter Search around 6.7 $\mu$eV,” Phys. Rev. Lett. 124 no. 10, (2020) 101802, arXiv:2001.05102 [hep-ex].
  * [91] S. Beurthey et al., “MADMAX Status Report,” arXiv:2003.10894 [physics.ins-det].
  * [92] D. Alesini et al., “Search for Invisible Axion Dark Matter of mass m${}_{a}=43~{}\mu$eV with the QUAX–$a\gamma$ Experiment,” arXiv:2012.09498 [hep-ex].
  * [93] M. Lawson, A. J. Millar, M. Pancaldi, E. Vitagliano, and F. Wilczek, “Tunable axion plasma haloscopes,” Phys. Rev. Lett. 123 no. 14, (2019) 141802, arXiv:1904.11872 [hep-ph].
  * [94] H.E.S.S. Collaboration, A. Abramowski et al., “Constraints on axionlike particles with H.E.S.S. from the irregularity of the PKS 2155-304 energy spectrum,” Phys. Rev. D 88 no. 10, (2013) 102003, arXiv:1311.3148 [astro-ph.HE].
  * [95] A. Ayala, I. Domínguez, M. Giannotti, A. Mirizzi, and O. Straniero, “Revisiting the bound on axion-photon coupling from Globular Clusters,” Phys. Rev. Lett. 113 no. 19, (2014) 191302, arXiv:1406.6053 [astro-ph.SR].
  * [96] A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi, and A. Ringwald, “Revisiting the SN1987A gamma-ray limit on ultralight axion-like particles,” JCAP 1502 no. 02, (2015) 006, arXiv:1410.3747 [astro-ph.HE].
  * [97] Fermi-LAT Collaboration, M. Ajello et al., “Search for Spectral Irregularities due to Photon–Axionlike-Particle Oscillations with the Fermi Large Area Telescope,” Phys. Rev. Lett. 116 no. 16, (2016) 161101, arXiv:1603.06978 [astro-ph.HE].
  * [98] M. D. Marsh, H. R. Russell, A. C. Fabian, B. P. McNamara, P. Nulsen, and C. S. Reynolds, “A New Bound on Axion-Like Particles,” JCAP 12 (2017) 036, arXiv:1703.07354 [hep-ph].
  * [99] D. F. Jackson Kimball et al., “Overview of the Cosmic Axion Spin Precession Experiment (CASPEr),” Springer Proc. Phys. 245 (2020) 105–121, arXiv:1711.08999 [physics.ins-det].
  * [100] T. Wu et al., “Search for Axionlike Dark Matter with a Liquid-State Nuclear Spin Comagnetometer,” Phys. Rev. Lett. 122 no. 19, (2019) 191302, arXiv:1901.10843 [hep-ex].
  * [101] A. Garcon et al., “Constraints on bosonic dark matter from ultralow-field nuclear magnetic resonance,” Sci. Adv. 5 no. 10, (2019) eaax4539, arXiv:1902.04644 [hep-ex].
  * [102] I. M. Bloch, Y. Hochberg, E. Kuflik, and T. Volansky, “Axion-like Relics: New Constraints from Old Comagnetometer Data,” JHEP 01 (2020) 167, arXiv:1907.03767 [hep-ph].