An even lighter QCD axion

We explore whether the axion which solves the strong CP problem can naturally be much lighter than the canonical QCD axion. The $Z_\mathcal{N}$ symmetry proposed by Hook, with $\mathcal{N}$ mirror and degenerate worlds coexisting in Nature and linked by the axion field, is considered in terms of generic effective axion couplings. We show that the total potential is safely approximated by a single cosine in the large $\mathcal{N}$ limit, and we determine the analytical formula for the exponentially suppressed axion mass. The resulting universal enhancement of all axion interactions relative to those of the canonical QCD axion has a strong impact on the prospects of axion-like particle experiments such as ALPS II, IAXO and many others. The finite density axion potential is also analyzed and we show that the $Z_\mathcal{N}$ asymmetric background of high-density stellar environments sets already significant model-independent constraints: $3\le\mathcal{N}\lesssim47$ for an axion scale $f_a\lesssim 2.4\times10^{15}$ GeV, with tantalizing discovery prospects for any value of $f_a$ and down to $\mathcal{N}\sim9$ with future neutron star and gravitational wave data, down to the ultra-light mass region. In addition, two specific ultraviolet $Z_\mathcal{N}$ completions are developed: a composite axion one and a KSVZ-like model with improved Peccei-Quinn quality.


Introduction
The axion experimental program is in a blooming phase, with several new experiments and detection concepts promising the exploration of regions of parameter space thought to be unreachable until a decade ago. Many of those experiments are simply prototypes, awaiting the jump to become 'big-experiments', or, in the case of more consolidated techniques, they are still far from saturating their full physics potential. Nonetheless, they sometimes reach sensitivities which go well-beyond astrophysical limits, albeit often still far from the customary QCD axion window.
On the other hand, since axion couplings are inherently ultraviolet (UV) dependent, such early stage experiments already provide valuable probes of the QCD axion parameter space. Imagine for definiteness that ALPS II would detect a signal in 2021, would it be possible to interpret that as an axion that solves the strong CP problem? Since the strong CP problem is one of the strongest motivations for new physics, if an axion-like particle (ALP) will be ever discovered, there or elsewhere, it would be compelling to explore whether it had something to do with the strong CP problem. This work explores whether wide regions in the ALP parameter space, well outside the traditional QCD axion band, may correspond to solutions of the strong CP problem. This is a question of profound theoretical and experimental relevance.
In axion solutions to the strong CP problem 1 both the axion mass and the couplings to ordinary matter scale as 1/f a , where f a is the axion decay constant, denoting the scale of new physics. The precise relation between mass and decay constant depends on the characteristics of the strong interacting sector of the theory. When QCD is the only confining group to which the axion a couples, in which case we denote the axion mass as m QCD a , they are necessarily linked by the relation [3,4] m QCD where χ QCD , m π , f π , m u and m d denote respectively the QCD topological susceptibility, the pion mass, its decay constant, and the up and down quark masses. Equation (1.1) is completely model-independent as far as QCD is the only source of the axion mass, and it defines the "canonical QCD axion", also often called "invisible axion". For this axion the aG µνG µν coupling to the gluon strength G µν is directly responsible for the axion mass, since the only source of explicit breaking of the global axial Peccei-Quinn (PQ) symmetry U(1) PQ is its QCD anomaly. The strength of other axion couplings to Standard Model (SM) fields is instead model-dependent: it varies with the matter content of the UV complete axion model.
In recent years there have been many attempts to enlarge the canonical QCD axion window, by considering UV completions of the axion effective Lagrangian which departed from the minimal DFSZ [5,6] and KSVZ [7,8] constructions. Most approaches actually focussed on the possibility of modifying the Wilson coefficient of specific axion-SM effective operators [9][10][11][12][13][14][15]. That is, the size of the coupling coefficients, at fixed f a , is modified. This has for example allowed to populate new regions of the parameter space by moving vertically the axion band in the axion mass versus coupling plane, see Fig. 1 left. The results are then "channel specific": different couplings c are modified differently. The parameter space of solutions can be alternatively changed by varying the axion mass at fixed f a . This corresponds to horizontal displacements of the canonical axion band in the parameter space, see right panel in Fig. 1. It always requires that the magnitude of the relation between the axion mass m a and 1/f a departs from that in Eq. (1.1): the confining sector of the SM must be enlarged beyond QCD. New instanton sources give then additional contributions to the right-hand side of Eq. (1.1). The practical consequence is a universal modification of the parameter space of all axion couplings at a given m a , at variance with the vertical displacement scenarios. This feature could a priori allow for the two mechanisms in Fig. 1 to be distinguished. 2 Examples of horizontal enlargement of the parameter space towards the right of the canonical QCD axion band are heavy axion models that solve the strong CP problem at low scales (e.g. f a ∼ TeV) [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. The present work explores instead left horizontal shifts: true axions that solve the strong CP problem with m a m QCD a . This avenue is more challenging, since it requires a new source of PQ breaking aligned with QCD, whose contribution to the axion mass needs to almost cancel that from QCD without relying on fine-tunings.
A possible mechanism to achieve this lighter-than-usual true axion in a technically natural way was recently put forth by Hook [35], in terms of a discrete Z N symmetry. N mirror and degenerate worlds would coexist in Nature, linked by an axion field 3 which implements non-linearly the Z N symmetry. One of those worlds is our SM one. All the confining sectors contribute now to the right-hand side of Eq. (1.1), conspiring by symmetry to suppress the axion mass without spoiling the solution to the strong CP problem.
The direct consequence is, for fixed f a , a N -dependent reduced axion mass in spite of all confining scales being equal to Λ QCD . In other words, for a given value of m a it follows a universal enhancement of all axion interactions relative to those of the canonical QCD axion. In this paper, we expand on the mathematical properties of the implementation of the Z N symmetry and determine the analytic form of the exponential suppression of the axion mass and its potential in the large N limit. The phenomenological analysis of the number of possible mirror worlds N will be next carried out with present and projected data.
The study will also explore the Z N axion potential at finite density, to confront present constraints and prospects from very dense stellar objects and gravitational waves. It has been recently pointed out in [37,38] that a generic reduced-mass axion leads to strong effects on those systems, raising the effective mass in the dense media. In the scenario considered here, a stellar background made only of SM matter is by nature Z N -asymmetric: we will show analytically how such an asymmetric background breaks the cancellations which guaranteed an exponentially suppressed axion mass for the Z N symmetric vacuum potential. Limits on the number of possible worlds will be obtained in turn.
The theoretical framework to be used throughout the work described above is that of effective axion couplings. Nevertheless, two concrete UV completions of the Z N scenario under consideration will be developed as well: a modelà la KSVZ [7,39], and a composite modelà la Choi-Kim [40,41]. The status of the Peccei-Quinn (PQ) quality problem will be also addressed.
An important remark is that we will consider in this paper experiments that can test the solution to the strong CP problem without further assumptions. Indeed, it is most relevant to get a clear panorama on the strong CP problem by itself, given its fundamental character. In particular, we will not discuss axion or ALP experiments that do rely on the assumption that the DM of the Universe may be constituted by axions. The cosmological evolution of the axion field in the Z N scenario under discussion and its contribution to the DM relic abundance departs drastically from the standard case, and it is discussed in a companion paper [42].
The structure of the present paper can be easily inferred from the Table of Contents. 2 Down-tuning the axion mass In Ref. [35] it was shown how to naturally down-tune the axion mass from its natural QCD value in Eq. (1.1), exploiting the analyticity structure of the QCD axion potential in the presence of a Z N symmetry. For pedagogical purposes, before turning to the generic Z N case we analyze the (unsuccessful) case of a Z 2 symmetry: the SM plus one degenerate mirror world linked by an axion which realizes the symmetry non-linearly.

The Z 2 case
Consider the SM plus a complete copy SM , related via a Z 2 symmetry which exchanges each SM field with its mirror counterpart, while the axion field is shifted by π: Figure 2: Z 2 axion potential. The mirror contribution to the axion potential V SM (a/f a ) (in green) partially cancels that of the SM, V SM (a/f a ) (in blue), leading to a total shallower potential V 2 (a/f a ) (in orange). The total potential has a maximum in a/f a = 0 and thus this Z 2 axion does not solve the SM strong CP problem.
The Lagrangian, including the anomalous effective couplings of the axion to SM fields, then reads where θ parametrizes the anomalous QCD coupling, α s is the QCD fine-structure constant, the Lorentz indices of the field strength G µν have been obviated, and the dots stand for possible Z 2 -symmetric portals between the two mirror worlds (see Sect. 2.2.1). Without loss of generality, we can perform a uniform shift in a such that the θ term in Eq. (2.3) is set to zero. Therefore, the effective θ-parameter of the SM corresponds to θ eff ≡ a /f a , where a denotes the vacuum expectation value (vev) of the axion field. In the case of an exact Z 2 symmetry, all couplings and masses of the mirror world and the SM would coincide with the exception of the effective θ-parameter. It is this difference (namely the π shift in the effective θ-parameters of the SM and its mirror) the one responsible for displaced contributions to the total axion potential, with destructive interference effects. Were the QCD axion potential to be a simple cosine, the total potential would vanish because the two contributions (from QCD and mirror QCD) would have exactly the same size but opposite sign, i.e. ∝ cos(a/f a ) and ∝ cos(a/f a + π) = −cos(a/f a ) respectively. However, for the true chiral axion potential [43][44][45] the exact cancellation disappears and a residual potential -and thus a non-zero axion mass-remains, which at leading chiral order reads (keeping only two flavours)

(2.4)
This Z 2 -symmetric world would not solve the strong CP problem, though, because a/f a = 0 is a maximum of the axion potential, as illustrated in Fig. 2. Indeed, as already pointed out in Ref. [35], a/f a = 0 is a minimum of the potential only for odd values of N , while it is a maximum for N . Thefore, the simplest viable axion model that solves the strong CP problem with a reduced axion mass incorporates a Z 3 symmetry.

Z N axion
We consider now N copies of the SM that are interchanged under a Z N symmetry which is non-linearly realized by the axion field: with k = 0, . . . , N − 1. One of those worlds will be our SM one. The most general Lagrangian implementing this symmetry describes N mirror worlds whose couplings take exactly the same values as in the SM, with the exception of the effective θ-parameter: for each copy the effective θ value is shifted by 2π/N with respect to that in the neighbour k sector, where L SM k denotes exact copies of the SM total Lagrangian excluding the strong anomalous coupling, and the dots stand for Z N -symmetric portal couplings that may connect those different sectors (to be discussed in Sect. 2.2.1). In this equation θ a ≡ a/f a is the angular axion field defined in the interval [−π, π), and a universal (equal for all k sectors) bare theta parameter has been set to zero via an overall shift of the axion field. The SM is identified from now on with the k = 0 sector: to ease the notation, the label k = 0 on SM quantities will be often dropped below. Each QCD k sector contributes to the θ a potential, which in the 2-flavour leading order chiral expansion reads and denotes the chiral condensate [44], while χ 0 ≈ (75 MeV) 4 is the zero temperature QCD topological susceptibility [45,46]. Alternatively, the total Z N axion potential can be written as For any N , θ a = 0 is an extrema of the axion potential. Indeed, using the property sin(2π(N − k)/N ) = − sin (2πk/N ) it is straightforward to see that The same holds for any θ a = 2πn/N with n ∈ Z, because of the periodicity of the potential. For N odd the potential V (θ a ) has N minima located at which includes the origin θ a = 0, while for N even the origin becomes a maximum. This result -valid for any N -can be shown for instance using the exact Fourier series expansion of the potential in Eqs. (2.8)-(2.11) (see final part of Appendix C). It follows that N odd is required in order to solve the SM strong CP problem (albeit with a 1/N tuning in the cosmological evolution [35,42]). The k = 0 worlds have instead non-zero effective θ-parameters: θ k ≡ 2πk/N for θ a = 0, see Eq. (2.7). A typical shape of the axion potential for N = 3 is illustrated in Fig. 3.  The different effective θ k values translate into slightly different masses for the pion mass in each mirror world, m π (θ k ). At quadratic order in m π a reduction factor of up to ∼ √ 3 results [45], (2.14) Interestingly, nuclear physics would be drastically different in the different mirror copies. In particular, a new scalar pion (π k ) to nucleon (N k ) coupling is generated in all worlds but the SM one (see e.g. Refs. [47,48]): where c + is an O(1) low-energy constant of the baryon chiral Lagrangian. Its impact on the cosmological histories of the mirror worlds is discussed in Ref. [42] for the Z N scenario under discussion.
Overall, for our world to be that with vanishing effective θ, the ∼ 10 orders of magnitude tuning required by the SM strong CP problem has been traded by a 1/N adjustment, while N could a priori be as low as N = 3. 4

Renormalizable portals to the SM
Renormalizable portals between the SM and its mirror copies (left implicit in Eq. (2.7)) are allowed by the Z N symmetry. In the following, we classify for completeness the portal operators connecting the different k sectors.

Higgs portals
The most general Z N symmetric scalar potential for the Higgs doublets H k of the different mirror worlds includes terms of the form where v denotes the Higgs vev and κ i are dimensionless parameters. Note that the Z Nsymmetric mixings between different worlds may include next-neighbour, next-to-next neighbour etc. interactions. All κ i≥1 terms provide renormalizable portals between the mirror Higgs copies (H k =0 ) and the SM Higgs (H k=0 ).

Kinetic mixing
Terms mixing the U(1) k Y hypercharge field strengths of mirror worlds are a priori also allowed by the Z N symmetry, where F µν k denote here the k-hypercharge field strenghts and i are free dimensionless parameters.
The above renormalizable portals are subject to strong cosmological constraints, as discussed in Ref. [42]. This can suggest a 'naturalness' issue for the Higgs and the 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 Poincaré symmetry [52,53]: in the limit where non-renormalizable interactions are neglected, the κ i =0 and i =0 → 0 limit corresponds to an enhanced P N symmetry (namely an independent space-time Poincaré transformation P in each sector). Those couplings are then protected from receiving radiative corrections other than those induced by the explicit P N breaking due to gravitational and axion-mediated interactions, which are presumably small. In addition, other terms in the scalar potential which depend on the details of the UV completion of the Z N axion scenario may be present and strongly constrained; an example is given below in Sect. 3.1.

Axion potential in the large N limit
It is non-trivial to sum the series which defines the axion potential, Eq. (2.11). However, the presence of the Z N symmetry allows for the application of powerful mathematical tools related to its Fourier decomposition and holomorphicity properties, that lead to simplified expressions in the large N limit.

Holomorphicity bounds and convergence of Riemannian sums
As first noticed in Ref. [35], the fact that the potential in Eq. (2.11) corresponds to a Riemann sum allows one to express it as an integral plus subleading terms, where the definition of each single-world potential, V θ a + 2πk N , can be read off Eq. (2.8). Most importantly, the integral does not depend on the field θ a and the amplitude of the axion potential is thus solely contained in the subleading terms. The latter are nothing but the error E committed in approximating the Riemann sum by an integral, Powerful theorems exist that describe the fast convergence of this approximation. It can be shown, applying complex analysis, that if some conditions are satisfied the convergence of the rectangular rule is exponential (see e.g. Section 3 in Ref. [54]). More precisely, if V (θ a ) is a 2π-periodic function and it can be extended to a holomorphic function V (w) in a rectangle from 0 to 2π and from −ib to +ib, then the error of the rectangular rule is constrained as where M is an upper limit on V (w) in the rectangular region defined above. As a consequence, the axion mass will be exponentially suppressed for large N . More in detail, let us apply the theorem to the second derivative of the potential, which can be extended in the complex plane to a holomorphic function until the expression under the square root vanishes. Indeed, this function has branch points in 5 Naively, it is tempting to apply the theorem assuming b = log z in Eq. (2.20). This is not possible though, since V (w) is not bounded in the rectangular region, due to a divergence in the branch point. As we show in Appendix A, it is possible to optimize the bound obtained above on the axion mass (V (θ a )/f 2 a ) by allowing a departure from log z, b = log z + ∆b, which leads to where the factor 3/2 stems from the order of the divergence of Eq. (2.21) in the branch point w cut . Implementing this result in Eq. (2.20), it follows that In Fig. 4 we compare this analytical bound with the numerical result: our analytical bound captures the correct dependence on N of the Z N axion mass, although it misses the overall constant factor. The overall factor will be analytically determined in the following. Nevertheless, the discussion above has the two-fold interest of determining the correct exponential suppression and of being very general, as it only relies on the holomorphicity structure of the potential, and not on the specific form it takes. As a consequence, the exponential suppression of the axion mass is not spoiled when considering the subleading chiral corrections to Eq. (2.11).

Fourier expansion: axion mass from hypergeometric functions
It is possible to gain further physical insight on the origin of the cancellations in the potential by constructing its Fourier series expansion. As shown in Appendix B, the Fourier series of any scalar potential respecting the Z N shift symmetry only receives contributions from modes that are multiples of N . Moreover, if the potential can be written as a sum of shifted contributions, as it is the case for the Z N axion under discussion -see Eq. (2.18)-then the Fourier series of the total potential V N (θ a ) can be easily obtained in terms of the Fourier series of a single V (θ a ) term, leading to whereV (n) denotes the coefficient of the Fourier series for the single-world potential V (θ a ),V It is convenient to express this integral in terms of the Gauss hypergeometric function (see Appendix C and Ref. [55] for conventions and relevant properties), As shown in Appendix C, in the large N limit this expression further simplifies tô leading to the following expression for the total potential where in the second line we have kept only the first mode in the expansion, as the higher modes are exponentially suppressed with respect to it. The total potential is thus safely approximated by a single cosine. It trivially follows from Eq. (2.30) that θ a = 0 is a minimum of the total potential for N odd, and a maximum for N even. Here and all through this work purely constant terms in the potential are obviated, as they have no impact on the axion mass. Eq. (2.30) can be rewritten as where the Z N axion mass m a in the large N limit is finally given by a compact and analytical formula, The overall coefficient is thus determined, in addition to exhibiting the z N exponential suppression of the potential and the specific N dependence previously argued in Eq. (2.25) from holomorphicity arguments. In summary, in the large N limit the axion mass is reduced with respect to that of the QCD axion by a factor where m QCD a denotes the mass of the canonical QCD axion as given in Eq. (1.1). 6 This ratio is illustrated in Fig. 4, which compares the numerical behaviour with: a) the analytical dependence previously proposed in Ref. [35]; b) that from the holomorphicity bound in Eq. (2.24); c) the full analytical result in Eq. (2.33). Our analytical results improve on previous ones by Hook on a number of aspects: i) the explicit determination of the exponential behavior controlled by z N ∼ 2 −N ; ii) the improved N dependence from the factor N 3/2 ; iii) the z-dependence of the axion mass in 1−z 1+z ; iv) the determination of the prefactor 1/ √ π. In practice, the large N results in Eqs. (2.30)-(2.33) turn out to be an excellent approximation already for N = 3.

UV completions and alternative scenarios
Up to this point, the analysis has been largely independent from the precise UV completion of the Z N axion scenario. For the sake of illustration, in this section we provide two UV completions of the axion effective Lagrangian in Eq. (2.7). We also briefly discuss an alternative implementation of the Z N symmetry in which the resulting axion is heavier than usual (rather than lighter).

KSVZ Z N axion
Consider N copies of vector-like Dirac fermions Q k (k = 0, . . . , N − 1) transforming in the fundamental representation of QCD k , together with a gauge singlet complex scalar S. The action of the Z N symmetry on these fields is postulated to be while the SM Lagrangian and its copies obey Eq. (2.5) under Z N . The most general Lagrangian containing the new degrees of freedom then reads where P R ≡ (1 + γ 5 )/2. It exhibits an accidental U(1) PQ symmetry

4)
S → e iα S , (3.5) that is spontaneously broken by the vev of S, v S , via a proper 'mexican-hat' potential V(S, H k ), whose structure is discussed below. Decomposing the S field in a polar basis, in terms of canonically normalized radial (ρ) and axion modes, the latter can be rotated away from the Yukawa term in Eq. (3.3) via an axion-dependent axial transformation (3.7) The heavy quarks, with real mass 7 m Q k = y S v S √ 2 , can next be integrated out in order to obtain the low-energy axion effective field theory. Because the transformation in Eq. (3.7) is QCD k anomalous, with anomaly factor 2N k = 1, the resulting axion effective Lagrangian is given by which yields precisely Eq. (2.7), after the identification v S = f a . Furthermore, the presence of the singlet scalar S introduces new scalar portals between the SM and its mirror worlds, in addition to the generic ones in Eq. (2.16). The scalar potential in the latter equation should thus be enlarged by (3.10) Note that, because the Higgs vev v is the same in all k sectors due to the unbroken Z N symmetry, the required hierarchy of scales is obtained with a single fine-tuning between v and f a , as for elementary canonical QCD axions. It is also possible to choose the representations of the Q k fields to transform nontrivially under the electroweak k gauge groups, so that they could e.g. mix with SM k quarks in a Z N invariant way and decay efficiently in the early Universe, thus avoiding possible issues with colored/charged stable relics in the SM sector [9,11]. Depending on the Q k quantum numbers, this would change in turn the value of the electromagneticto-QCD anomaly ratio of the PQ current, usually denoted as E/N , which enters the axion-photon coupling.

Peccei-Quinn quality
The threat posed on traditional QCD axion models by quantum non-perturbative gravitational corrections [56][57][58][59][60][61][62][63][64][65] may also affect the models discussed here, as f a is not very far from the Planck scale. These contributions are usually parametrized via effective operators, suppressed by powers of the Planck mass, that could explicitly violate the PQ symmetry and thus spoil the solution to the strong CP problem [56][57][58][59]. 8 7 Note that we crucially removed also the k-dependent phases from the Yukawas, in order to properly integrate out the heavy Q k fields. 8 UV sources of PQ breaking can be avoided in some invisible axion constructions within a variety of extra assumptions or frameworks [66][67][68][69][70][71][72][73][74][75][76][77][78], or be arguably negligible under certain conditions [79]. In the context of the KSVZ Z N axion model above, the exponentially small axion mass could seem to worsen this threat, increasing the sensitivity to explicit PQ-breaking effective operators. Interestingly, promoting the in built Z N symmetry to a gauge symmetry leads to an accidental U(1) PQ invariance, that for large N is efficiently protected from those extra sources of explicit breaking. Indeed, the lowest-dimensional PQ-violating operator in the scalar potential compatible with the Z N symmetry is S N , leading to an explicitly PQ-breaking contribution to the potential of the form where M Pl = 1.22 × 10 19 GeV is the Planck mass and c is a dimensionless coefficient with phase δ ≡ Arg c. Considering now V N (θ a ) + V PQ−break. , expanding for small θ a the axion potential V N (θ a ) ≈ V N (0)+ 1 2 m 2 a f 2 a θ 2 a , and solving the tadpole equation, the induced effective θ parameter in the SM sector reads where m 2 a from Eq. (2.32) has been used, and in the last step we neglected the second term in the denominator in the first line of Eq. (3.12): this is always justified in the θ a 10 −10 regime. In summary, unlike the customary ad-hoc Z N protection mechanism for the standard KSVZ axion, in the Z N axion scenario under discussion the discrete symmetry is already present by construction. Note that the scaling with N is slightly different as compared to the standard KSVZ axion, due to the enhancement factor 1/z N . But eventually the (f a /M Pl ) N suppression dominates and provides an efficient protection mechanism, even though the axion mass is exponentially suppressed. For the sake of an estimate, Fig. 5 shows the regions in the {N , f a } plane that saturate the nEDM bound for |c| = 1 and sin δ = 1.

Composite Z N axion
It is also possible to construct a UV completion of the Z N scenario which corresponds to a dynamical (composite) axionà la Kim-Choi [40,41], without extending its exotic fermionic content. In the original version of that model, the SM fields are not charged under the PQ symmetry while two exotic massless quarks, ψ and χ, transform under an extra confining "axi-color" group SU( N ) a and one of them, ψ, is also a triplet of QCD. Upon confinement of the axi-color group at a large scale Λ a ∼ f a Λ QCD , pseudo-Goldstone bosons composed of the exotic quarks emerge. All but one of them are coloured under QCD and become safely heavy. The light remaining one is the composite axion, whose mass obeys the usual formula for QCD axions Eq. (1.1).
We implement the Kim-Choi idea in the framework of our Z N framework without increasing the number of massless exotic fermions representations. The fermion ψ is simply extended to be now a triplet under all QCD k mirror sectors, see Table 1. The axion field will thus be unique and will couple to all anomalous terms. Upon SU( N ) a confinement at the large scale of order f a , the QCD k couplings α k s can be neglected, and therefore a large global flavor symmetry arises in the exotic fermionic sector: SU(3 N + 1) L × SU(3 N + 1) R × U(1) V . 9 This symmetry is spontaneously broken down to SU(3 N + 1) L+R × U(1) V by the exotic fermion condensates. Among the resulting Goldstone bosons, the QCD k singlet corresponds to the composite axion. Its associated PQ current reads (with f PQ ≡ N f a ) which corresponds to the only element of the Cartan sub-algebra of SU(3 N + 1) that has a vanishing anomaly coefficient with respect to SU( N ) a , but a non-vanishing one with respect to all the QCD k gauge groups.
Without further elements the model would be viable, but all mirror worlds would have the same θ-parameter: a heavier than usual axion would result. A simple Z N implementation which leads instead to relatively shifted potentials, and thus to a reduced axion mass, is to have a relative phase between the argument of the determinant of the quark mass matrix of the mirror worlds, where Y u and Y d denote the Yukawa matrices for the up and down quark sectors, respectively. One of the many possible Z N charge assignments for the quarks that yield Eq. (3.15) is that in which only the right-handed up quarks would transform as 10 corresponding to a Yukawa quark Lagrangian of the form The resulting low-energy axion effective field theory is then the desired one as in Eq. (2.7). In this Z N composite axion model only the exotic fermions are charged under the PQ symmetry, while the Z N charges are carried solely by SM quarks. This means that the Z N and PQ symmetries are not directly linked. As a consequence, gauging the inbuilt Z N symmetry would not soften the PQ quality problem, contrary to the case of the KSVZ Z N axion model discussed earlier above. Our Z N composite axion model is then subject to the usual PQ quality threat. Standard softening solutions often applied to composite axion models could be explored, e.g. those based on a chiral gauging of the global symmetry of the coset space or on introducing a moose structure [66,69,73,80].

Ultra-light QCD axions
The term ultra-light axions usually refers to the mass range m a ∈ 10 −33 , 10 −10 eV (with the extrema of the interval corresponding respectively to an axion Compton wavelength of the size of the Hubble horizon and to the Schwarzschild radius of a stellar mass black hole). As a theoretical motivation for ultra-light axions, the so-called string Axiverse [81] is often invoked, according to which a plenitude of ultra-light axions populating mass regions down to the Hubble scale 10 −33 eV is a generic prediction of String Theory, although without a direct reference to the solution of the strong CP problem. 11 On the other hand, according to the usual QCD mass vs. f a relation, Eq. (1.1), axion masses below the peV correspond to axion decay constants larger than the Planck mass, and hence they are never entertained within canonical QCD axion models. The Z N axion framework discussed in the present work allows in contrast to populate the sub-peV axion mass region while keeping sub-Planckian axion decay constants, with the advantage of providing as well a direct solution to the strong CP problem. As shown in Sec. 4.2, the tantalizing prospects for testing the Z N scenario, through observational data on very dense stellar objects and gravitational waves, can sweep through the discovery region of the ultra-light axion range.

A heavier-than-QCD axion
A remark is in order regarding the Z N charge of the axion in the different sectors. If the implementation of the Z N symmetry would be such that the N world replicas are degenerate but the axion field is exactly the same in all of them, that is, if Eqs. (2.5)-(2.6) were replaced by a −→ a , (3.19) the potentials of the different mirror worlds would not be relatively shifted but exactly superpose. The axion would then be a factor √ N heavier than the usual QCD axion in Eq. (1.1). This scenario was proposed in Ref. [36] for a Z 2 symmetry with just one mirror world degenerate with the SM, but its generalization to N copies is trivial. Such a heavier-than-QCD axion solution is viable, and it would transform the ALP arena to the right of the canonical QCD axion band into solutions to the SM strong CP problem. The axion Z N charge assignment explored throughout this work, Eq. (2.6), results instead in lighter-than-QCD axions, that is, solutions located to the left of the QCD axion band. Note that this option induces a comparatively much larger impact: a natural exponential suppression of the axion mass ∝ z N as the byproduct of the cancellations between the mirror potentials, Eq. (2.32), instead of the mild √ N enhancement just discussed. All in all, to explore the right-hand side region of the QCD axion band for solutions to the strong CP problem, other heavy axion scenarios proposed in the literature seem more efficient and appealing (e.g. those with mirror worlds much heavier than the SM, or scenarios with novel confining scales much larger than Λ QCD , as mentioned in Sect. 1).

Experimental probes of down-tuned axions
The Z N axion with reduced mass can provide a solution to the SM strong CP problem, independently of whether it accounts or not for the DM content of the Universe. It is hence interesting to get a perspective on the experimental panorama that does not rely on the supplementary assumption that the axion may be the DM particle: all experimental bounds and prospects below will be independent of that hypothesis. On the other hand, Ref. [42] will focus on experimental probes that do rely on it.

Axion coupling to photons
From an experimental point of view, a highly relevant axion coupling is that to photons, defined via the Lagrangian term δL = 1 4 g aγ aFF as [45,103] where E and N denote model-dependent anomalous electromagnetic and strong contributions, respectively. Fig. 6 shows the parameter space of the reference Z N axion model (with E/N = 0) in the coupling vs. mass plane. Predictions for the axion photon coupling are obtained by rescaling the Z N axion mass in Eq. constrains [83][84][85][86][87][88][89][90] and astrophysical bounds [91][92][93][94][95][96][97][98][99][100][101] are shown in blue and green, respectively. Projected sensitivities appear in translucent colors delimited by dashed lines. The orange oblique lines represent the theoretical prediction for the Z N axion photon couplings assuming E/N = 0 for different (odd) number of worlds N . These lines are solid for the regions of the parameter space in which the KSVZ UV completion of the Z N axion is free from PQ quality problem and dashed otherwise. The secondary vertical axis shows the corresponding axion decay constant f a if E/N = 0 is assumed. Supplementary constraints in case the axion is assumed to account for DM can be found in Ref. [42]. Axion limits adapted from Ref. [102].
the left of the parameter space: each of those oblique lines can be considered to be the center of a displaced yellow band. It is particularly enticing that experiments set a priori to only hunt for ALPs may in fact be targeting solutions to the strong CP problem.

Finite density constraints on f a
This subsection summarizes the model-independent constraints on f a for the Z N scenario under discussion. The result of the analysis is illustrated in Fig. 8. Interestingly, apart from the usual constraints stemming from the SN1987A [104] and black hole superradiance measurements [105][106][107][108] (depicted in purple), novel bounds apply to the exceptionally light Z N axion due to finite density effects. Indeed, it has been recently pointed out in Refs. [37,38] that finite density media may have a strong impact on the physics of very light axions or ALPs. In those media, the minimum of the total potential may be shifted to π. This has a number of phenomenological consequences that span from the modification of the nuclear processes in stellar objects due to θ ∼ O(1), to modifications in the orbital decay of binary systems (and subsequently in the emitted gravitational waves). For the scenario considered here, the important point is that a background made only of ordinary matter breaks the Z N symmetry. This hampers the symmetry-induced cancellations in the potential which led to a reduced-mass axion in vacuum: the effective axion mass will be larger within a dense medium.
We will first elaborate on the Z N axion potential in a nuclear medium. Following Refs. [37,109], one can compute the finite density effects on the axion potential by considering the quark condensates in a medium made of non-relativistic neutrons and protons [110,111]. Applying the Hellmann-Feynman theorem, the quark condensate at a finite density n N of a given nucleon N can be expressed as where qq 0 = 1 2 uu + dd ≡ −Σ 0 is the quark condensate in vacuum -see Eq. (2.10)and σ N is defined by where m q ≡ 1 2 (m u +m d ) and M N is the mass of the nucleon N . Because the Z N potential  Figure 7: Example of the in-medium potential dependence as a function of the nuclear density for N = 5. For large densities (light green) the total potential develops a minimum in θ a ∼ π.
is proportional to the quark condensate, see Eq. (2.8), we can simply obtain the potential within a SM nuclear medium V f.d. N (θ a , n N ) by weighting the SM (i.e. k = 0) contribution in the vacuum potential by the factor in Eq. (4.2), that is, In the last step of these expressions the large N limit has been taken, which allowed us to neglect the term corresponding to the exponentially reduced axion potential in vacuum (see Eq. 2.30). This shows that, if the nucleon density is large enough, the Z N asymmetric background spoils the cancellations among the mirror world contributions to the potential, in such a way that the total potential in matter is proportional to minus the SM one in vacuum V (θ a ). Therefore, the minimum of the potential is located at θ a = π. More precisely, which requires for the minimum to sit at θ a = π. This is illustrated in Fig. 7.
A large value of the θ parameter inside dense stellar objects is rich in physical consequences, which translates into strong constraints for the Z N scenario. As it was pointed out in Ref. [37], θ ∼ O(1) inside the solar core is excluded due to the increased protonneutron mass difference (which would prohibit the neutrino line corresponding to the Be 7 -Li 7 mass difference observed by Borexino [112]). Similarly, for θ ∼ π in nearby neutron stars (NS), Co 56 would be lighter than Fe 56 [47,48] and therefore Fe 56 could have been depleted due to its β-decay to Co 56 . The presence of iron in the surface of neutron stars and its implications in terms of the allowed θ values could be explored through dedicated X-ray measurements [113]. The corresponding current and projected constraints that were derived in Ref. [37] (within the simplifying assumption z = 1) are translated here to the Z N scenario and further generalized for any z.
A conservative criterion consistent with θ = π inside the medium is to impose that the axion mass at θ a = 0 becomes tachyonic, i.e. −m 2 (4.7) Requiring this quantity to be positive, it directly follows a limit on the number of worlds allowed by the stellar bounds above: where the most recent estimation of σ N has been used. 12 This bound does not apply for the whole range of f a , though, because the argument only makes physical sense as long as the reduced Compton wavelength of the axion is smaller than the stellar object, ∼ 1/f a is the effective axion mass in the medium, (4.9) 12 We employ here σ N 59 MeV which is in agreement with recent determinations based on Roy-Steiner equations σ N = 59.1(3.5) MeV [114] and ChPT estimates σ N = 59 (7) MeV [115].
For the case of the sun, r core ∼ 139.000 km implies f a 2.4 × 10 15 GeV for the observational constraints to apply. Finally, the area in parameter space excluded is depicted in dark blue in Fig. 8. Analogously, the future sensitivity prospects from neutron star data are depicted in shaded blue. 13 Additional constraints which apply if the axion is assumed to account as well for DM are discussed in Ref. [42].
Even stronger bounds may be established by relaxing the requirement stemming from Eq. (4.7). Indeed, as it can be seen in Fig. 8, long before the mass in θ a = 0 becomes tachyonic, the absolute minimum of the potential corresponds to θ ∼ O(1). Therefore one could constrain larger masses or smaller N values in the Z N scenario than those obtained above. This would require, however, a dedicated analysis to ensure that the axion field would fall into the absolute minimum, so as to overcome the potential barrier; this development lies beyond the scope of the present work.
The fact that the position of the minimum of the axion potential depends on the nuclear density of the medium not only modifies the effective θ-parameter inside stellar 13 Our results are analogous to those in Eq. (1.7) of Ref. [37], with their generic parameter identified as = m 2 a /m 2 a (N = 1) π −1/2 1 − z 2 (1 + z) N 3/2 z N −1 . Note that the location of the QCD axion line, as well as our projected exclusion regions for neutron stars and gravitational waves, are shifted towards the left by a factor of five with respect to those in Refs. [37,38]. objects but may also source a long-range force between them. This has been studied in Refs. [37,38]. This new long-range force sourced by the axion can be constrained by the measurement of double pulsar or neutron star (NS) -pulsar binaries [116][117][118]. Moreover, the existence of these axionic long-range forces would modify the gravitational wave signal emitted by NS-NS mergers or black hole-NS mergers which will be probed in the future by LIGO/VIRGO and Advanced LIGO [38,119]. The projected constraints from Ref. [38] are show in green in Fig. 8. It is striking that the whole ultra-light DM region, included the so-called "fuzzy dark matter" region (m a ∼ 10 −22 eV) [120], will be within observational reach in the next decades, for a wide range of N values.

Conclusions
An axion which solves the strong CP problem may be much lighter than the canonical QCD axion, down to the range of ultra-light axions, provided Nature has a Z N symmetry implemented via N degenerate world copies, one of which is our SM. The axion field realizes the symmetry non-linearly, which leads to exponential cancellations among the contributions from each mirror copy to the total axion potential. For large N , we have shown that the total axion potential is given by a single cosine and we determined analytically the -exponentially suppressed-dependence of the axion mass on the number of mirror worlds, using the properties of hypergeometric functions and the Fourier expansion. In practice, the formula in Eq. (2.32) gives an excellent approximation even down to N = 3. We have also improved the holomorphicity bounds previously obtained.
We compared next the predictions of the theory with present and future data from experiments which do not rely on the additional assumption that an axion abundance may explain the DM content of the Universe. It is particularly enticing that experiments set a priori to hunt only for ALPs may in fact be targeting solutions to the strong CP problem. For instance, ALPS II is shown to be able to probe the Z N scenario here discussed down to N ∼ 25 for a large enough axion-photon coupling, while IAXO and BabyIAXO may test the whole N landscape for values of that coupling even smaller, see Fig. 6. In turn, Fermi SN data can only reach N 43 but are sensitive to smaller values of the coupling.
Highly dense stellar bodies allow one to set even stronger bounds in wide regions of the parameter space. These exciting limits have an added value: they avoid model-dependent assumptions about the axion couplings to SM particles, because they rely exclusively on the anomalous axion-gluon interaction needed to solve to the strong CP problem. A dense medium of ordinary matter is a background that breaks the Z N symmetry. This hampers the symmetry-induced cancellations in the total axion potential: the axion becomes heavier inside dense media and the minimum of the potential is located at θ a = π. From present solar data we obtain the bound N 47 provided f a 2.4 × 10 15 GeV, while larger N values are allowed for smaller f a . Furthermore, we showed that projected neutron star and pulsar data should allow to test the scenario down to N ∼ 9 -and possibly even below-for the whole range of f a , see Fig. 8. Furthermore, gravitational wave data from NS-NS and BH-NS mergers by LIGO/VIRGO and Advanced LIGO will allow to probe all values of N for the remaining f a range, up to the Planck scale and including the ultra-light axion mass range.
These analytical and phenomenological results have been derived within the modelindependent framework of effective couplings. Nevertheless, for the sake of illustration, we have developed two examples of UV completed models. One is a Z N KSVZ model, which is shown to enjoy an improved PQ quality behaviour: its Z N and PQ symmetries are linked and in consequence gauging Z N alleviates much the PQ quality problem. The other UV completion considered in this paper is a Z N version of the composite axionà la Kim-Choi. While this model is viable, its PQ quality is not improved with respect to the usual situation, because its Z N and PQ symmetries are independent.
This work is intended to be a proof-of-concept that a much-lighter-than usual axion is a viable solution to the strong CP problem, with spectacular prospects of being tested in near future data. It also pinpoints that experiments searching for generic ALPs have in fact discovery potential to solve the strong CP problem.
The down-tuned axion considered here could also explain the DM content of the Universe in certain regions of the parameter space. The impact of such a light axion on the cosmological history is significant and it will be discussed in a separate paper [42].

A Holomorphicity properties of Z N axion potential
In order to determine the parameter b in Eq. (2.20), which controls the exponential suppression of the axion mass, it is necessary to study the region in the complex plane where the extension of the potential V (w) is holomorphic. As the plots in Fig. 9 illustrate, both the potential and its second derivative have branch cuts starting in w cut = π ± i log z. However, the second derivative V (w) diverges at the branch point and thus b cannot be extended all the way to log z. In order to optimize the bound on the axion mass we allow b to depart from log z, b = | log z + ∆b |. Taking into account that V (w) for small ∆b can be approximated by the procedure amounts to minimize the function B(∆b) that determines the bound |E N (V )| ≤ B(∆b) (see Eq. (2.20), The requirement dB(∆b) d(∆b) = 0 shows that the bound is optimized for where the factor 3/2 comes form the order of the divergence in Eq. (A.1). Figure 9: Representation of the complex functions V (w) (left) and V (w) (right). Colors represent the phase of the corresponding complex function and the brightness represents the modulus. The singularities can be clearly spotted: branch cuts starting from w cut = π ± i log z in both functions and divergences in those same points for V (w).

B Fourier series of the Z N axion potential
We show here that the coefficients of the Fourier series of any Z N symmetric potential, such as the Z N axion potential in Eq. (2.8), vanish unless the corresponding Fourier mode is a multiple of N . Moreover it will be shown that, when the potential is expressed as the non-vanishing coefficients of the Fourier series can be expressed in terms of the Fourier transformation of a single term in the sum Eq. (B.1). Let us denote byV N (n) the coefficients of the Fourier series of the total potential, and byV 2πk/N (n) the coefficients of the Fourier series of each of the terms in the sum in Eq. (B.1), We will stick to the notation that omits the subindex for the first world (k = 0),V 0 (n) ≡V (n), where the factor of two comes from the negative modes and the constant term (i.e. θ a -independent) has been obviated.

C Analytical axion mass dependence from hypergeometric functions
We show here that the Fourier series coefficients of the single world axion potential in Eq. can be written for large Fourier modes, n 1, as a simple analytical formula that exponentially decays with n. Moreover, by applying the result in Appendix B, it will be shown that the potential for the Z N axion approaches a single cosine and a simple formula for the Z N axion mass follows.

(C.3)
For convenience, the hypergeometric function can be also expressed as (see Eq. (9.131) from Ref. [55]) so that 2 F 1 −1/2, n − 1/2 n + 1 z 2 = 1 − z 2 1/2 2 F 1 −1/2, 3/2 n + 1 The relation in Eq. (C.3) is exact. However, only the modes n which are multiples of N contribute to the potential (see Appendix B), and therefore it is pertinent to focus on the large n limit. While the limit of the Gauss hypergeometric function when one or more of its parameters become large is difficult to compute in general, some asymptotic expansions of the hypergeometric function are known in the literature. In particular, following Ref. [121,122], Putting all this together, it follows that, in the large n limit, the coefficient of the Fourier series for a single world is given byV (n) = (−1) n m 2 π f 2 π 2 √ π 1 − z 1 + z n −3/2 z n , (C. 8) which in turn implies in this limit that the total Z N potential in Eq. (B.10) can be written as This expression allows us to understand several properties of the total potential. Firstly, it can be shown now that the total potential approaches a single cosine in the large N , since all the other modes are then exponentially suppressed with respect to the first one, Secondly, we can also obtain an analytical expression for the axion mass that confirms the dependence obtained from the holomorphicity arguments in Section 2.3.1, and completes the expresion with the correct prefactor, Finally, it is now trivial to show that the potential in the large N limit has N minima (maxima) located at θ a = {±2π /N } for = 0, 1, . . . , N −1 2 , for odd (even) N . The results above assumed the large N limit. However, note that the conclusion about the location of the extrema is true for any N . This can be easily seeing after obtaining the exact Fourier expansion of the Z N axion potential in Eq. (2.8), which reads, (C. 13) For even N , it trivially follows that θ a = 0 is a maximum, as all factors in this expression are positive except for the factor (1 − 2 ) < 0. For odd N instead the (−1) t N factor fluctuates in sign, but the first term (t = 1) is positive and dominates the expansion (e.g. it is exponentially larger in magnitude than the t = 2 term which is negative). The periodicity of the potential extends these conclusions to the location of all extrema.