Dark matter from an even lighter QCD axion: trapped misalignment

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_\mathcal{N}$ symmetry. The whole mass range from the canonical QCD axion down to the ultra-light regime is allowed, with $3\le\mathcal{N}\lesssim65$. This includes the first proposal of a"fuzzy dark matter"QCD axion with $m_a\sim10^{-22}$ eV. A novel misalignment mechanism occurs --{\it trapped misalignment}-- due to the peculiar temperature dependence of the $Z_{\mathcal{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.

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 and m QCD a , f a , f π , m π , 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.
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 N symmetry realized non-linearly by the axion field [31,32]. 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 θ k -parameters, Here, L SM k denotes exact copies of the SM total Lagrangian excluding the topological G kGk term, 1 θ a ∈ [−π, π) is defined in terms of the axion field a, θ a ≡ a/f a , and the dots stand for Z N -symmetric portal couplings among the SM copies. The solution to the strong CP problem of this Z N scenario required N to be odd. Overall, the ∼ 10 orders of magnitude tuning required by the SM strong CP problem is traded by a 1/N adjustment, where N could be as low as N = 3 (viable solutions were found with 3 ≤ N 47 for f a 2.4 × 10 15 GeV, and any N for larger values of f a [32]). The resulting axion is exponentially lighter than the canonical QCD one, because the nonperturbative contributions to its potential from the N degenerate QCD groups conspire by symmetry to suppress each other. This can be intuitively understood from the large N limit of a non-linearly realized Z 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 N as befits a U (1) Goldstone boson. Indeed, it has been shown [32] that in the large N limit the total axion potential is given in all generality by where now the axion mass in vacuum obeys (1.4) which is exponentially suppressed (∝ z N ) in comparison to (m QCD a ) 2 . 2 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 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 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 N reduced-mass axion, because of the same-size enhancement of all axion couplings in such a scenario. Note that an exceptionally light QCD axion has formerly been considered [45], albeit it required an ad-hoc fine-tunning of potential parameters.
It will be shown that the cosmology of the Z 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. [46]). In particular, the constraints on the number of effective relativistic species N eff present in the early Universe imply that the mirror copies of the SM must be less populated -cooler-than the ordinary SM world [47]. 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 N symmetry may arise naturally in the cosmological evolution [48][49][50]. 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 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. 3 While gravity and axion-mediated interactions are naturally small enough, the impact of other possible (Z N -symmetric) interactions on the thermal communication between the SM and its mirror copies will be discussed.
The evolution of the Z 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 2 Note that the value of m a for N = 1 should be equal to m QCD a given in Eq. (1.1), while Eq. (1.4) only holds in the large N limit. 3 An axion heavier than the QCD one has been contemplated as DM, within a Z 2 setup which realized linearly the symmetry [51]. Z 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 θ = π), which effectively delays the onset of the true oscillations and thus enhances the DM density. 4 We will call this new production mechanism trapped misalingment. While an early stage of oscillations has been previously proposed in the literature in order to adiabatically suppress the axionic energy density and/or to solve the isocurvature problem [54][55][56][57], our trapped misalingment differs from previous works in that the modification of the potential (from the high temperature large axion mass to the exponentially suppressed zero temperature mass) is strongly nonadiabatic and, when active, always leads to an enhancement of the relic density. Note that inflation played a crucial role in previous dynamical setups advocated to drive the initial misalignment to π [58][59][60], 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 [61]. 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 [61] required an ad-hoc Peccei Quinn (PQ)breaking non-renormalizable effective operator suppressed by powers of the Planck mass, M Pl . In contrast, we will show that the early stage of oscillations in the Z 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. [62], 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 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 10 −10 eV), within the Z N axion framework under discussion. The present and projected sensitivity to the number of possible mirror world 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.

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. [46] In order to satisfy these constraints in the presence of one mirror world, its temperature T needs to be smaller than that of the SM, T 0 . In the Z 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 ρ k and entropies s k , 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 kdependence through the respective θ 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 eff . This is typically the case for f a 10 7÷8 GeV [66][67][68] and if model-dependent axion couplings to SM fermions are neglected [69]. The different sectors evolve in time with separately conserved entropies, therefore the ratio of entropy densities γ 3 k ≡ (s k /s) is constant, while the ratio of temperatures is given by where from now on T will denote the temperature of the SM copy, T ≡ T 0 . The Hubble expansion rate depends on the total energy density of the Universe. In a radiation dominated era with N mirror worlds it reads with the number of effective degrees of freedom g * (T ) given by 5 For BBN, the constraint corresponds to BBN + Y p + D in Table V of Ref. [63], translating their 68.3% CL result for the number of neutrino species N ν to a 95% CL value N ν = 2.85 ± 0.28, and N eff = 1.015 N ν . For CMB, the combination of TT+TE+EE+lowE+lensing+BAO data by Planck 2018 [64,65] is considered here; note that due to the disagreement between local and CMB measurements of the Hubble expansion parameter H 0 , the constraint on ∆N eff would weaken to N eff = 3.27 ± 0.30 if H 0 measurements were included.
Let us discuss next the implications for N eff . In the SM, g * (T = 1 MeV) = 5.5 + 7 4 N SM eff = 10.83 where N SM eff = 3.046 is the effective number of neutrino species. Analogously, at recombination g * (T = 0.26 eV) = 3.93. It follows from Eq. (2.5) that the contribution of the ensemble of N worlds to ∆N eff = N eff − N SM eff at BBN and CMB is given by 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 SM 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.
The temperature dependence of the function g s is quite flat in most regions of interest, 1, and the parameter γ k directly gives the ratios of temperatures which remain constant throughout the Universe evolution, with γ k < 1 for all k = 0. This approximation is valid as far as the temperatures involved correspond to the same plateau of the N eff distribution (see e.g. [70]). 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 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 γ k is not too small, e.g. better than 25% for γ k =0 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 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 max , If all mirror images of the SM are instead assumed to have the same temperature, T k ∼ T for k = 0, the most restrictive bounds follow, It is interesting that BBN data set the most constraining bound on N 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 eff from BBN and CMB. Specifically, the scenario predicts in all generality the following discrepancy: 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 =0 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 =0 T . For instance, in the Z 2 mirror case T /T 0.5 requires both portal couplings, defined as L ⊃ κ|H| 2 |H | 2 + B µν B µν , where H and H (B µν and B µν ) denote respectively the SM Higgs doublet (hypercharge field strength) and its mirror copy, and κ and are dimensionless couplings, to respect κ , 10 −8 [71]. Even smaller couplings are needed in the Z N case with 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 [72,73]: in the limit where non-renormalizable interactions are neglected, the κ, → 0 limit for the ensemble of portals corresponds to an enhanced P N symmetry (namely, an independent space-time Poincaré transformation P in each sector). This protects those couplings from radiative corrections other than those induced by the explicit P N breaking due to gravitational and axion-mediated interactions. The former are presumably small being Planck suppressed, while axion-mediated corrections to κ scale like m 2 H /f 2 a [74] and hence they can be safely neglected for the standard high f a values considered. Moreover, axion-mediated interactions among the different sectors (leading to interaction rates scaling as ∼ T 5 /f 4 a ) are also small enough during the evolution of the Universe, such that they do not spoil the evolution of the independent thermal baths, as long as the PQ breaking pre-inflationary scenario is considered.

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 µνG µν pseudo-scalar density in Eq. (1.2): the effective value of the θ parameter differs for each sector k, θ k = 2πk/N (and thus relaxing to zero in the SM with probability 1/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, where π k and N k denote respectively the pion and nucleon fields and c + is an O(1) lowenergy constant, induces long-range forces (i.e. not spin-suppressed in the non-relativistic limit) among nucleons in all worlds but the SM one [75,76]. 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 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 [48,49]. Two inflaton fields would be present: φ for the SM and φ for its Z 2 mirror copy, respectively slow-rolling towards the minimum of their potentials (the generalization to the case of Z 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 φ reaches its minimum first, any mirror particles produced due to its oscillations are diluted by the remaining inflation driven by the second field φ. Hence, by the time when φ reheats the SM sector, the density of mirror matter will be very small, and in consequence T T . Quantitatively, whether such asymmetric inflation scenario can be realized or not depends on the shape of the inflaton potential [49]. An alternative idea to achieve different temperatures is that the inflaton itself transforms non-linearly under the Z N symmetry [50]: 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 =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. [46]). 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 [47]. Nevertheless, the abundance of mirror baryons strongly depends on the specific baryogenesis mechanism at play and thus in the following, we will not dwell with mirror baryonic DM (which is assumed to give a negligible contribution 6 due to the lower temperature of the SM mirror copies) and focus instead on axionic DM.

Z 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 N scenario with a reduced-mass axion and the canonical QCD axion.
An important question is whether the Z N axion can account for the observed DM density, without requiring any further adjustment other than the inherent 1/N tuning needed to solve the strong CP problem, i.e. to choose θ a = 0 as the vacuum for the SM world among the N possible θ a values, θ a = {±2π /N } with = 0, 1, . . . , (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.

Temperature dependence of the Z N axion potential
At temperatures well below the critical temperature of the QCD phase transition, T T QCD 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 QCD . Lattice computations of the topological susceptibility χ(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 (θ a , T ) by weighing the zero temperature expression of the chiral axion potential V (θ a ) by a h(T ) factor Constant terms in the potential have been left implicit in this expression, as irrelevant for the axion mass and cosmological evolution. The coefficient α parametrizes the powerlaw behaviour at high temperatures. Unless otherwise stated, α = 8 will be assumed, as suggested by the dilute instanton gas approximation (DIGA) and in agreement with some lattice computations [77]. 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 (θ a , T ) is a continuous function of the temperature not only around the minimum but for any field value θ a . The generalization of the ansatz above to the total Z N axion potential is simply given by the sum over the N different potentials, 7 where the power-law temperature suppression is assumed to be k-independent and as given in Eq. (3.2), that is, 3) and Appendix A. Irrespective of the specific γ k =0 values, three temperature regimes can be distinguished which are qualitatively different. They are determined by the minimum and maximum value of γ k =0 ,

Low temperatures T < T QCD
The potential for all N replicas is well approximated by the -exponentially suppressedzero-temperature potential in Eqs. (1.3) and (1.4), with minima at 2πk/N for N odd.

Medium temperatures T QCD < T < T QCD /γ max
In this regime, only the SM contribution to the total potential is temperature suppressed. It is weighed down by a factor T QCD /T α , while all other contributions are still well approximated by their zero-temperature potential, where in the last step the term ∝ z N has been neglected compared to the size of V (θ 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 θ a = ± π.

High temperatures T > T QCD /γ min
The temperature effects are now important for all N replicas. In the simplified case with all N − 1 mirror worlds having the same temperature ( where in the second step we neglected again the exponentially suppressed Z N contribution at T = 0 and the definition T ≡ T γ b was used. As for the medium temperature regime, the total axion potential depends essentially only on f a and its minimum is located at θ a = ± π. 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 ± π although typically not too far from it, see  In summary, temperature effects lead to regimes in which the minimum of the total axion potential is displaced to θ a = ± π or similar large field values. This hightemperature shift of the minimum parallels that at finite density, previously identified in Refs. [32,78,79] and recently employed as well in Ref. [80]. This phenomenon is rich in consequences. In particular, it follows that the axion mass at finite temperature, is given -in the three temperature regimes discussed above-by 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, which differs from m QCD a by a factor ∼ 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.9) can be rewritten as (3.11) In summary, at low temperatures the Z N axion mass depends both on f a and N and is exponentially suppressed by a factor z N /2 . In contrast, at medium and high temperatures it only depends on f a in the large N limit, and is of the order of the canonical QCD axion mass in vacuum and at high temperature respectively.

Z N axion misalignment
We analyze here the consequences of the temperature-dependent Z 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 θ a in an expanding Universe:θ where V (θ a ) = dV (θ a )/dθ a and the spatial gradients have been neglected. For a field oscillating near its minimum, V (θ a ) m 2 a (T )f 2 a θ 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 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 QCD (feeling the z N exponential suppression of the potential) and reaches then a constant value. This is shown in Fig. 3, where the Z N trajectories are depicted in blue for different values of N . The figure illustrates the N -independence of m a (T ) at high and medium temperatures (in the large N limit), while it becomes N -dependent at low-temperature, see also Eq. (3.9) and Fig. 2. At some point in that trajectory, the mass overcomes the Hubble friction -depicted in green-and the axion starts to oscillate.
There are two qualitatively distinct regimes, depending on whether this crossing happens after or before the axion mass abruptly drops at T QCD (as illustrated in Fig. 3). Let us denote by f crit a the value of f a for which the canonical QCD axion would start to oscillate precisely at T QCD , Whenever m a > m crit (m a < m crit ) Hubble crossing takes place just (more than) once in the cosmic trajectory. The factor 1.67 for m 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 crit a for the Z N axion is slightly larger than that in Eq. (3.13), see Eq. (3.10); 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 crit a corresponds to the case in which the (green) Hubble line would cross the canonical QCD axion trajectory (i.e. N = 1, orange) exactly at the flattening point of the latter, i.e. ∼ T QCD .

One stage of oscillations: simple ALP regime
For f a > f crit a 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 θ 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. 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 θ 1 due to Hubble friction, H 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 (3.14) Since T 1 < T 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: T 1 > T . The axion fields keeps oscillating around θ a = 0 until nowadays and its mass is almost constant,ṁ a 0. It is adequate to apply the WKB approximation provided m a ṁ a /m a ,Ḣ/H , 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 where a denotes the scale factor and ρ 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 ρ a, 0 can be expressed as where a 0 and a 1 denote respectively the value of the scale factor today and at the onset of oscillations, a 1 ≡ a(T 1 ). 8 The latter can be expressed in turn in terms of the axion mass using Friedmann equations and the conservation of entropy, a 1 ∝ 1/ √ m a M Pl , allowing us to obtain the ratio of the axionic relic density to the current DM abundance ρ DM as a function of m a , θ 1 and f a , 9 The initial axion misalignment θ 1 can a priori take any value in the [−π, π) interval and the axion field will roll down to the closest minimum 10 among the N possibilities. Thus it is with probability 1/N that the final θ-parameter will correspond to the CPconserving point θ = 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 θ 1 ∈ [−π/N , π/N ) within the simple ALP regime, for any f a > 2 × 10 17 GeV. Eq. (3.17) also indicates that the simple ALP solutions which account for the relic DM obey (see Fig. 7) 20) to be compared with the Z N axion mass relation, that for a given N predicts m a f a ∝ const. This behaviour is depicted in the {m a , 1/f a } plane in Fig. 6 by continuous superimposed lines (for different values of θ 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 N solutions can solve both the strong CP problem and account for DM mass down to the region of fuzzy DM (m a ∼ 10 −22 eV) and for any value of N ≥ 3.

Two stages of oscillation: trapped misalignment
For f a < f crit a , two different potential minima develop during the Universe evolution from high to low temperatures, see Fig. 3 (right panel). In a nutshell: 8 Whenever the initial misalingment angle is large, the relic density in Eq. (3.17) needs to be corrected by the so-called anharmonicity factor, see Appendix B. 9 The slight discrepancy with respect to the result in Ref. [81] is due to the use of the updated value of ρ DM 1.26 keV/cm 3 from Planck 2018 data [65] 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. [82,83]). 10 Note that the relevant axion field displacement that should enter in Eq. (3.19) corresponds to the field distance to the closest minimum and it is only when the closest minimum is θ a = 0 that the field distance coincides with θ 1 .
• A first phase of oscillations takes place around θ a = π at temperatures above T QCD , because temperature effects are important then and the axion mass is unsuppressed.
• At T QCD the true minimum θ a = 0 develops (while π becomes a maximum). The axion mass suddenly becomes exponentially suppressed due to the Z N symmetry.
Oscillations around θ 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.14). Before Hubble crossing at T 1 , the axion field is frozen and remains constant at an arbitrary initial misalignment angle θ 1 .
T 1 > T > T QCD . The finite temperature effective potential is relevant in this range, see Eqs. (3.6) and (3.7), and the initial axion velocity isθ 1 = 0. Thus, the axion "is trapped" in this first stage of oscillations around θ a ∼ π from T 1 until T QCD . Its mass is unsuppressed due to the explicit Z N breaking by the thermal background: m a (T ) raises progressively until it stabilizes at ∼ m QCD a . Trapping causes a delay of the start of oscillations around the true minimum θ 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. N = 1, orange) is shown to oscillate around θ a = 0 since well before T QCD . In contrast, for N ≥ 3 (blue) the oscillations around θ = π take place from T 1 down to T 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 θ a andθ a at T QCD , θ tr ≡ θ a (T QCD ) ∼ π,θ tr ≡θ a (T QCD ) . (3.21) It is not possible to predict the exact value ofθ tr . However, an order of magnitude estimate stems from the energy density when the temperature approaches T QCD , applying N a conservation to the adiabatic interval [T 1 , T QCD ), which shows that the precise value of ρ a, QCD depends on the arbitrary initial misalignment angle with respect to the high-temperature minimum in π: (θ 1 − π). 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, , and thus θ a = π becomes a maximum. The abrupt exponential suppression of the axion mass was shown in Fig. 3 for different values of N . It illustrates that the two stages of oscillation are separated by a sudden, inherently non-adiabatic, modification of the potential asṁ a /m a m a . The energy density of the axion just after the transition at T QCD reads ρ a,tr = 1 2θ that is, the total energy density at this point is larger or equal than the maximum potential energy (the height of the barrier 2m 2 a f 2 a /N 2 ). What happens next depends crucially on the value ofθ 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,θ tr 2m a /N , oscillations around one minimum may start closely after T QCD . This case will be denoted below "pure trapped misalignment".
• For large enoughθ tr so that the kinetic energy dominates ρ a,tr ,θ tr 2m a /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 QCD . This case will be denoted as "trapped+kinetic misalignment".
3.2.2.1 Pure trapped misalignment:θ tr ∼ 0 T QCD > T . Were the axion velocity exactly zero at T QCD ,θ tr = 0, the setup would not be viable: the axion would simply roll to the closest minimum θ a = π N (N − 1). This is precluded by the experimental limit on the nEDM that implies θ 10 −10 .
Nevertheless,θ 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 O(N ) potential barriers 11 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 θ a = 0 with 1/N probability.
An estimation of the time required for the axion to go from one maximum to the next is δt ∼ 2π/m a . To overfly only ∼ N barriers a very short time is required because in this regime m a H, i.e. the difference between T QCD and the temperature at the onset of oscillations around the true minimum is negligible. In conclusion, forθ tr 2m a /N the kinetic energy is not large enough to extend beyond ∼ T QCD the delay time accumulated during the trapped period (see Eq. (3.24)). The value of θ a when it starts to oscillate from the top of the barrier closest to θ a = 0 is 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 QCD , the adiabatic approximation is valid again sinceṁ a 0, m a H. From the N a conservation law -Eq. (3.16), and the initial conditions at the onset of the second stage of oscillations θ 2 π/N ,θ 2 0 , it follows that the current axionic relic density ρ a,0 is given in the pure trapped case by where a QCD ≡ a(T QCD ). For a given final vacuum mass m a , it demonstrates an enhancement of the relic density in the scenario with a Z N trapped axion vs. that with a simple ALP, for two main reasons: i) the relative scale dependence is larger by a factor (a QCD /a 1 ) 3 and does not depend on the axion mass; ii) the initial misalignment angle is fixed to θ tr ∼ π/N , while for the simple ALP it is arbitrary within the interval [−π/N , π/N ). Note that the current relic density in the pure trapped scenario becomes independent of the initial misalignment θ 1 , unlike the case of the simple ALP and QCD-like pre-inflationary scenarios.
For the Z N reduced-mass axion under discussion, the product m a f a -and thus ρ a,0 in Eq.   T QCD > T > T 2 . If the axion velocity at T QCD is large enough,θ a 2m a /N , 13 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 QCD .
In this regime, the WKB (adiabatic) approximation is not valid anymore and the axion energy density is diluted as ∝ a −6 due to the expansion of the Universe (as illustrated in 13 Equivalently, if the axion kinetic energy at the end of the trapped periodθ 2 tr f 2 a /2 is much larger than the maximum potential energy 2m 2 a f 2 a /N 2 . Fig. 5). This can be shown by noting thatθ a f 2 a is the Noether charge density associated to the PQ symmetry [61] and therefore it decays asθ a ∝ a −3 as long as the axion mass can be neglected. 14 In consequence, a comoving PQ charge q kin can be defined which is a conserved quantity, q kin =θ a a 3 = const . (3.30) 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, At this point the axion will start a second stage of oscillations, see Fig. 5. Enforcing the conservation of q kin from T QCD until T 2 , and using Eq. (3.31), it follows that the scale factor at the onset of the second stage of oscillations a 2 ≡ a(T 2 ) can be expressed as (3.32) T 2 > T . The axion finally starts the second stage of oscillations around one of the N zero-temperature minima at random. Thus, the 1/N probability of the Z N axion scenario [31,32] to solve the strong CP problem (i.e. the final θ-parameter to corresponding the CP-conserving point θ = 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 is constant during this final period (and different from that during the first oscillation period N a ). Using q kin in Eq. (3.32) and the mean axion velocity at T QCD in Eq. (3.23), N a is well approximated by [61] Finally the relic density today, can be written as 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 2,    which confirms the result first obtained in Ref. [61]; 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.26) shows two competing effects. In fact, for the kinetic+trapped scenario: • the dilution factor is smaller, ( √ a 1 a QCD /a 0 ) 3 vs. (a QCD /a 0 ) 3 ; • the relevant mass scale is larger, m 1 m QCD a,π vs. m a . Hence, it will depend on the specific point of the parameter space whether the ki-netic+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 where F 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 ρ a,0 ρ DM tr+kin ∝ m a f 19/12 a , (3.38) in contrast to m a f a = const. for the pure trapped case in Eq. (3.26). The trapped+kinetic trajectories that account for DM are illustrated by double lines in Fig. 6, for different values of the initial misalignment angle θ 1 .

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 O(N ) maxima after T QCD . The result will be of general interest whenever f a < f crit a , and of practical relevance mainly for the pure trapped misalignment. Since the trapping forces θ tr π, the corresponding total energy density is always larger than the maximum potential energy (the height of the barrier 2m 2 a f 2 a /N 2 ). As the axion rolls over ∼ N different maxima, the energy density is diluted due to the expansion of the Universe, where a QCD+N is the scale factor when the axion has rolled over just N maxima after T QCD , and p parametrizes how strong is the dilution (we find numerically that p ∼ 1.5). The minimal kinetic energy K min necessary for the axion to be able to roll over just N maxima before oscillating must saturate the condition Given the Hubble parameter in a radiation dominated Universe H(t) = 1 2t and the Friedmann equation in Eq. (2.4), it follows that the temperature T QCD+N after overcoming N maxima (corresponding to a time difference of δt ∼ 2πN /m a , see discussion before Eq. (3.25)) is where t QCD is the time corresponding to T QCD , and the approximation 2πN /m a t QCD 16 has been used, as in this regime H m a . The corresponding dilution factor reads Substituting it in Eq. (3.40), the minimum kinetic energy required is determined, To evaluate how probable it is to simultaneously solve the strong CP problem and explain DM, K min needs to be compared with the mean kinetic energy generically expected at T QCD . It follows from Eq. (3.23) that In the case in which all mirror copies but the SM have the same temperature T and taking for example γ = 0.2, the latter equation can be rewritten as where κ is a numerical factor that ranges in the interval (4 − 8). This implies that the axion field will roll over at least N barrier tops whenever  [84] 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.46) is not satisfied. This precludes any DM solution with N 9, although only for m a m crit 4 × 10 −11 eV.
• The stars indicate for each trajectory the point below which f a is large enough (f a > f crit a ) 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 4 a ∝ 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 crit a 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 crit a (but very close to it) can be clearly seen. Above that point, the slope of the trajectory obeys m a f 19/12 a = const., as explained around Eq. (3.38).  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.
Overall, Fig. 6 shows that good solutions to both the strong CP problem and the nature of DM exist within the Z N axion scenario whenever m a 10 −4 eV, down to the fuzzy DM region (m a ∼ 10 −22 eV). In the range m crit m a 10 −4 eV the number of worlds must be N ∼ 21. For the lighter axion scenarios, m a m crit , values of N up to ∼ 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 [85]. The consequences of the fragmentation of the axion field within the trapped misalignment mechanism is left for a future work.

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 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 hightemperature 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-adiabatic 17 and takes place at a smaller temperature than the oscillation temperature in absence of trapping. 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.20), (3.26) and (3.38)). It also shows the generic m a dependence of the axion-SM interaction couplings g aXX ∼ 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 N axion there is no regime where the QCD-like relic density formula applies other than for N = 1). In Table 1 (and in all figures), for the QCD-like and trapped+kinetic cases the DIGA value α = 8 has been assumed (see Eq. (3.2)). 18 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 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 ∼ 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. α = 8. The dashed line reminds that the pure trapped Z N scenario can only take place for m a > m crit .

Implications for axion dark matter searches
In this section we discuss the implications of the above Z 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 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 10 −10 eV) is considered, including the fuzzy DM regime down to m a ∼ 10 −22 eV. 18 For arbitrary α, ρ a ∝ m a f (14+3α)/(2(4+α)) a in the trapped+kinetic case, while ρ a ∝ m (α+2)/(α+4) a f 2 a in QCD-like misalignment, from which the corresponding m a dependence of g aXX ∝ 1/f a can be obtained.
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 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 [32,78,79], 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 ρ a provides only a fraction of the total DM relic density ρ DM , the sensitivity to couplings of axion DM experiments should be rescaled by (ρ a /ρ DM ) 1/2 .
Crucially, we have shown in Sect. 3 that in some regions of the {m a , f a } plane the Z 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 N -KSVZ axion model developed in Ref. [32]. The canonical KSVZ QCD axion solution (i.e. 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 N , that is, for a Z N reduced-mass axion, see Eq. (1.4).
The entire DM relic density can be accounted for within the Z N axion paradigm in the regions encompassed by the purple band in the figures. 19 These correspond to initial values of θ a (from the misalignment mechanism) which range from θ 1 = 3/N down to θ 1 = 0.003/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 θ 1 (in the latter case, through its dependence on the axion velocity at T QCD ).

Axion coupling to photons
The effective axion-photon-photon coupling g aγ is defined via the Lagrangian δL ≡ 1 4 g aγ aFF , (4.1) with [86] g aγ = α 2π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 N -KSVZ axion model [32] are depicted for E/N = 0. The conclusion is that, for a large range of values of N , an axion-photon signal can be expected in large portions of the parameter space of upcoming axion DM experiments, if the Z N axion accounts for the entire DM relic density.   [92][93][94][95][96][97][98][99][100][101] and astrophysical bounds (green) [102][103][104][105][106]. Projected sensitivities appear in translucent colors. The orange oblique lines represent the theoretical prediction for solutions to the strong CP problem within the Z N paradigm, for the benchmark axion-photon couplings E/N = 0 and for different (odd) numbers of SM copies N .
In the purple band area the Z N axion can account for the entire DM density, in addition to solve the strong CP problem.

Axion coupling to nucleons
The axion coupling to nucleons N (N = p, n), is defined via the Lagrangian term with [86] (see also [107] for similar results) 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 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 [87] include: axion DM experiments (blue) [108][109][110][111][112] and astrophysical bounds (green) [113,114]. Projected sensitivities appear in translucent colors, delimited by dashed lines. The orange oblique lines represent the theoretical prediction for the Z N axion-neutron couplings, for different (odd) numbers of SM copies N . The purple band encompasses the region where the Z N axion can account for the entire DM density. The two theoretical predictions are represented for the benchmark value c 0 q = 0.

Axion coupling to electrons
The axion coupling to electrons, is defined via the Lagrangian term δL ≡ C ae ∂ µ a 2f a eγ µ γ 5 e , (4.8) with [66,115] where in the reference Z 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 axionmagnon conversion techniques [116,117] could barely detect a Z 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 [118], will be able to test this scenario in the future.   [87] include: solar axion experiments (dark green) [119,120], axion DM experiments (blue) [116,117] and astrophysical bounds (green) [121] and hints (grey band) [122]. Projected sensitivities appear in translucent colors, delimited by dashed lines. The orange oblique lines represent the theoretical prediction for the Z N axion-electron couplings assuming c 0 e = 0 and E/N = 0, for different (odd) numbers of SM copies N . The purple band encompasses the region where the Z N axion can account for the entire DM density.

Axion coupling to nEDM
The axion coupling to the nEDM operator provides a crucial test of the Z 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 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 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 [84], where C anγ = 0.011(5) e [123] and ρ DM, local ≈ 0.4 GeV/cm 3 is the local energy density of axion DM. 20 For instance, the axion DM experiment CASPEr-Electric [35,108,124] 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 [125] are also illustrated there, together with present constraints which exclude N 65. Fig. 11 demonstrates the brilliant prospects of CASPEr-Electric Phase I and II to discover a Z N axion in a large range of light masses and for 3 ≤ N 65 if such axion accounts for DM. In particular, the Z 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 N and θ 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 [32,78,79]. For the Z N axion scenario under discussion, they implied a strict limit on the number of mirror worlds: 3 ≤ N 47 for f a 2.4 × 10 15 GeV. They also provide tantalizing discovery prospects for any value of f a and down to N ∼ 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] (BBN constraints [78] are not shown, as they do not overlap with the purple region where DM is accounted for). It is thus of crucial interest to compare them with the oscillating nEDM projects described above. Fig. 11 tends to this task and  [84,108,125] astrophysical (green) [126][127][128][129][130] and finite-density constraints (orange) [32,78,79]. The nEDM bounds have been relaxed by a factor of 3 to take into account the stochastic nature of the axion field, see Ref. [131]. 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 N axion can account for the entire DM density.
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.

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 ∈ 10 −22 , 10 −10 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 Pl . In contrast, the Z N axion scenario is -to our knowledgethe 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 [132]). The Z 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. [133,134]). This implies a bound around the recombination era on any additional energy density fluctuations beyond cold DM, δρ/ρ 10 −3 , corresponding to a 0 (z rec )/f a 10 −3 [135], with a 0 (z rec ) denoting the axion field amplitude at the recombination redshift z rec . Expanding the axion potential around the minimum, V (a) 1 2 m 2 a a 2 − 1 4! λ a a 4 + . . . , the bound translates into a condition on the quartic coupling [135] λ Ref. [135] considered such constraint in the context of the Z N axion model of Ref. [31]. We here revisit the latter analysis with the formulae derived in Ref. [32] for the Z N scenario under discussion. For the latter, the potential in the large N limit -Eqs. (1.3) and (1.4)-corresponds to a quartic coupling given by 13) and the bound in Eq. (4.12) translates into N 85. This constraint does not impact the DM prospects discussed in this paper, as N 65 values are already excluded by either finite density constraint or nEDM data, see Fig. 11.
Finally, we mention that the pre-inflationary PQ breaking scenario considered here is subject to potential CMB constraints on iso-curvature fluctuations. Since the dynamics of the axion field evolution is highly non-linear it is not obvious to track the evolution of the iso-curvature fluctuations from the end of inflation till recombination. We note, however, that these bounds are relaxed for low-scale inflation H I 10 11 GeV (see e.g. Ref. [45], for an analysis based on a fine-tuned model).
In summary, the Z 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 ∈ 10 −22 , 10 −10 eV. For instance, Fig. 8 shows that a photon-axion signal could be at reach for the upper masses in that range, provided θ 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 θ 1 and down to N ∼ 9.

Conclusions
This work constitutes a proof-of-concept that an axion lighter -or even much lighterthan 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 canonical QCD axion band can accomplish that goal.
While the implications of a Z N shift symmetry to solve the strong CP problem (with a 1/N probability and 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 3.2 × 10 17 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. θ = 0, to large values, e.g. θ = π, where the axion field gets trapped down to a temperature T ∼ T 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 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 θ = 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 min required at the end of trapping for the axion to roll over ∼ 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 O(K 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 ∼ T QCD .
In fact, the axion kinetic energy at the end of trapping is shown to be in general much larger than K min . Trapped misalignment then automatically seeds kinetic misalignment [61] between T ∼ T QCD and lower temperatures. The axion rolls for a long time over the low-temperature potential barriers before final oscillations start at T T 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 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 [61]. It is instead directly seeded by trapped misalignment, which is itself a pure temperature effect.
For values of the Z N axion scale f a 3.2 × 10 17 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 θ 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 , 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 N ∼ 21 to account for the ensemble of DM. Overall, DM solutions are found within the Z N paradigm for any value of 3 ≤ N 65.
The results above have been next confronted with the experimental arena of the socalled axion DM searches. As a wonderful byproduct of the lower-than-usual f a values allowed in the Z 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 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 aGG 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 N paradigm, which is impossible in any model of the canonical QCD axion. It follows that the Z N paradigm is -to our knowledgethe 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 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 N and m a values, in the mass range from the standard QCD axion mass down to ∼ 10 −22 eV, that is, down to the fuzzy DM range. To our knowledge, the Z N axion discussed here is the first model of fuzzy DM which also solves the strong CP problem.
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 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 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 mirror k temperatures for different values of γ k . To this aim, the tabulated effective number of degrees of freedom (for the SM case) as a function of temperature provided in Ref. [136] will be used (see also Fig. 12). is shown to approach 1 for T k T QCD in the right plot in Fig. 13. Therefore the approximation b k (T, T k ) ∼ 1 can be safely adopted in the computation of the high-temperature axion potential in Eq. (3.7).
Bounds on γ k through N eff : c k (T, T k ) The ratio which appears in Eq. (2.5), has been depicted in Fig. 14. The left plot shows that its value is c k (T = 1 MeV, T k ) 1.2 whenever the SM BBN is taking place, i.e. T ∼ 1 MeV. In consequence, the approximation c k (T, T k ) ∼ 1 can be safely applied when using the experimental limits in Eq. (2.1) to extract bounds on γ k from Eq. (2.6).   Trapped+kinetic misalignment: F kin, 2 (T 1 , T 1 ) The prediction for the relic density in the kinetic+trapped misalignment scenario in Eq. (3.36) involves the ratio  where T BBN denotes the temperature of all the mirror copies of the SM when the temperature of the latter is T BBN 1 MeV.
The numerical evaluation of F 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 γ k .

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 θ π/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 anh . For an ALP-like regime, one can use an empirical formula proposed in Ref. [62] which encodes this deviation from the harmonic result: The definition of f anh (θ 1 ) thus depends on the criteria chosen to define the onset of oscillations at temperature T 1 . In our case, T 1 corresponds to the temperature in which the axion mass term overcomes Hubble friction as given in Eq. (3.14). For the case of a simple ALP with constant mass and a cosine potential, and for large misalignment angles π − θ 1 10 −2 , our numerical result is in excellent agreement with the Ref. [62] proposal. However, the latter does not apply for small misalignment angles, as it does not have the expected limiting behaviour f anh (θ 1 ) θ 1 →0 − −−− → 1. For this reason, we will employ the following anharmonicity prescription: 3 1.3 for π − θ 1 > 10 −2 , 0.06 log 2.32 π−|θ 1 | + 4 log log 2.32 π−|θ 1 | for π − θ 1 < 10 −2 , where the first expression is inspired by Ref. [137], 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. [62].  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 N , by choosing the value of γ that saturates the N eff bound in Eq. (2.6). 21 In the case of the Z N axion mass in Eq. where F kin, 2 (T 1 ) can be found in Eq. (A.4).