The QCD axion sum rule

Axions and axion-like particles (ALPs) are the subject of intense exploration at present, with both novel experimental proposals and theoretical work constantly being put forward. Being pseudo-Goldstone bosons (pGBs), their search constitutes an epic and generic quest to uncover symmetries hidden in Nature and awaiting discovery.

The paradigmatic example of this type of particle is the true axion that solves the strong CP problem of the Standard Model of Particle Physics (SM), that is, the fact that the value of one parameter of the strong force (QCD) - the $\bar{\theta}$ parameter- has to be adjusted by over ten orders of magnitude to comply with experimental bounds. In the successful wake of explaining small parameters through symmetries, the axion would be the pseudo-Goldstone boson of a dynamical solution via a global axial $U(1)_{A}$ symmetry. This “Peccei-Quinn” (PQ) symmetry must be classically exact although hidden (aka spontaneously broken), and explicitly broken only by QCD quantum effects which give the axion $a$ a tiny mass $m_{a}$Peccei:1977hh ; Peccei:1977ur ; Weinberg:1977ma ; Wilczek:1977pj . Indeed, the anomalous $a\,G_{\mu\nu}\tilde{G}^{\mu\nu}$ coupling to the gluon field strength $G_{\mu\nu}$ is directly responsible for its mass $m_{a}$, which is then necessarily linked to the axion scale $f_{a}$ by the relation¹¹1Which takes into account the mixing with the $\eta^{\prime}$ and subsequently the $\pi^{0}$.

where $\chi_{\rm QCD},m_{\pi},f_{\pi},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) defines the “canonical QCD axion”, also often called the “invisible axion” or “standard axion”. It corresponds to a straight line in the $\{m_{a},1/f_{a}\}$ parameter space: the QCD axion band. Its importance stems from its presumed universality, i.e. it is presumed to hold for whatever ultraviolet (UV) axion model, as far as QCD is the only source of Peccei-Quinn breaking. What lies beyond that band is considered to be ALP territory. ALPs are pGBs similar to axions but a priori they are not associated to a solution of the strong CP problem. Indeed the interest in axions extends far beyond the QCD axion; ALPs arise in a variety of theories (e.g.  Dienes:1999gw ; Gelmini:1980re ; Davidson:1981zd ; Wilczek:1982rv ; Svrcek:2006yi ; Arvanitaki:2009fg ; Cicoli:2013ana ; Alexander:2023wgk ) and are also attractive dark-matter candidates Preskill:1982cy ; Abbott:1982af ; Dine:1982ah . In practice, the mass and scale of generic ALPs do not need to abide by Eq. (1), and they are treated as independent free parameters.

If an ALP signal is ever detected, an immediate and compelling question will be whether it can nevertheless be interpreted as a true axion that solves the strong CP problem. Indeed, much theoretical effort is being dedicated to modify the relation in Eq. (1) so as to obtain for instance a heavier-than-QCD or lighter-than-QCD true axion. In other words, to obtain solutions in which the QCD band is displaced significantly towards the right rubakov:1997vp ; Berezhiani:2000gh ; Gianfagna:2004je ; Hsu:2004mf ; Hook:2014cda ; Fukuda:2015ana ; Chiang:2016eav ; Dimopoulos:2016lvn ; Gherghetta:2016fhp ; Kobakhidze:2016rwh ; Agrawal:2017ksf ; Agrawal:2017evu ; Gaillard:2018xgk ; Hook:2019qoh ; Gherghetta:2020ofz ; Csaki:2019vte ; Kivel:2022emq or the left  Hook:2018jle ; DiLuzio:2021pxd ; DiLuzio:2021gos of its standard position. This achievement requires in general the enlargement the gauge group of the forces in Nature beyond the SM ones, providing new PQ violating contributions.

In this paper we explore in depth a much simpler possibility. The physical axion mass eigenstate does not need to coincide with the scalar field that couples to $G_{\mu\nu}\tilde{G}^{\mu\nu}$ (apart from the mixing with the $\eta^{\prime}$ and $\pi^{0}$ via QCD). In other words, the axion-gluon interaction basis and the mass basis do not need to be simultaneously diagonal. Here, generic mixing of the axion with extra scalars is allowed –alike to the mismatch of electroweak interaction versus mass bases for fermions. The axion field (or combination of fields) is thus allowed to mix with other $N-1$ arbitrary SM singlet scalars. As long as the mixing potential is classically PQ invariant, the solution to the strong CP problem holds. The condition that has to be fulfilled by a general mixing potential in order to be PQ invariant will be identified. It will be shown in all generality that the axion properties are shared among the $N$ eigenstates, with each of them obeying an equation which departs from Eq. (1) because the extra scalar potential provides extra contributions to their masses. That is, one single axion field (or combination of fields) will reveal itself in Nature as multiple axion signals. Specifically, $N$ non-standard but true QCD axions will result, each of them with a different mass $m_{i}$ and different effective axion scale $f_{i}$ as determined from its coupling to gluons. We will denote by $g_{i}$ the factor by which each of them deviates from the standard relation,

The system is tightly constrained, though, since the PQ symmetry dictates how the mass induced by the topological susceptibility of QCD $\chi_{\rm QCD}$ is shared among the $N$ axion eigenstates. Denoting by $\beta_{i}\equiv 1/g_{i}$ the “axionness” of each eigenstate $a_{i}$, that is, the fraction of its mass which is due to the QCD topological susceptibility, it will be demonstrated that

This sum rule, together with other mathematically exact conditions and relations will be proven below in all generality.

The constraint in Eq. (3) is very rich in consequences, as it links the properties of a given axion eigenstate to the fate of the other $N-1$ axions. Among other aspects of key importance for the experimental search, we will address and solve the question of what is the maximal possible displacement from the standard QCD band for the axion eigenstate closest to that band. The area in $\{m_{a},1/f_{a}\}$ parameter space where these multiple QCD axion signals can be detected will be determined. While for low values of $N$ the displacement for a given axion may be drowned on the error bands of the projections, the prediction of up to $N$ different axion signals is a direct smoking gun, for instance for experiments which measure the axion coupling to the neutron electric dipole moment (nEDM) operator with very high precision in frequency, such as CASPER-electric Budker:2013hfa ; JacksonKimball:2017elr , as the axion-to-nEDM coupling follows directly from the axion gluonic coupling.

The strength of other axion couplings to Standard Model (SM) fields is model-dependent. An important one is the axion-photon coupling, explored in a plethora of experiments. The modification of the latter within the $N$ axion framework under discussion will be addressed as well and shown to also follow a displacement pattern, albeit subject to model-dependences.

An underlying condition for the results in this paper to have strong experimental impact is that the contribution to the masses induced by the generic mixing potential should not be vastly different from the QCD-induced mass. Whenever this is not the case for some axion eigenstates, these will be of no impact (i.e. either too heavy or very light but in both cases decoupled from $G_{\mu\nu}\tilde{G}^{\mu\nu}$). That is, they will decouple from the sum rule in Eq. (3) and thus will not significantly contribute to the solution to the strong CP problem. In this perspective, the usual single canonical QCD axion is just one particular case of the vast parameter space of solutions, in which all $N-1$ extra eigenstates have decoupled.

A mixing of the canonical QCD axion with other singlet scalars in Nature has previously appeared in past publications Kim:2004rp ; Choi:2014rja ; Kaplan:2015fuy ; Giudice:2016yja ; DiLuzio:2017ogq ; Chen:2021hfq ; Agrawal:2022lsp ; Fraser:2019ojt ; Darme:2020gyx but, either by choice or by construction, the limit in which all but one axion decouples was taken and the features discovered here were not discussed.

The problem will be formulated within the model-independent framework of EFTs and for a generic scalar potential, and the main mathematical tool to be used is the eigenvector-eigenvalue theorem JacobiDeBQ ; Denton:2019pka . Some UV complete toy examples will be also shown for illustration, though, in App. A.

Before proceeding to formulate the problem for $N$ singlet scalars and the most general PQ-invariant potential, we illustrate some results in this section within a $N=2$ EFT toy model. The reader interested in the general case and solutions can go directly to Sec. III. In what follows, the notation $m_{i}$, $f_{i}$ will be reserved to denote the mass and the scale of a given mass eigenstate $a_{i}$, and $F$ will denote a certain combination of all $f_{i}$, see below.

Let us consider a combination of two fields with the following effective interactions:²²2The pattern of results to be obtained next will hold as well for a general mixing potential, i.e. $V(\hat{a}_{1},\hat{a}_{2})$ instead of the second term in this equation, see Sect. III.

Note that in this Lagrangian there is only one combination of fields that couples to $G\widetilde{G}$, ${{a}_{G\widetilde{G}}}/{F}\equiv({{\hat{a}}_{1}}+{{\hat{a}}_{2}})/{\hat{f}}$, which is not necessarily a mass eigenstate, as it mixes with the orthogonal scalar combination, $a_{\perp}$, via a PQ invariant potential, i.e.

It is easy to see that, in this particular example, it is the shift of $\hat{a}_{1}$ that implements the PQ symmetry and ensures CP conservation in the QCD vacuum. Expanding the axions around their minima we obtain the following potential,

which corresponds to the mass matrix,

where $r\equiv\hat{m}_{2}^{2}\hat{f}^{2}/\chi_{\rm QCD}$. Given the simplicity of the system it can be exactly diagonalized. The axion masses read

and the corresponding physical eigenstates are given by

Generically, two axion mass eigenstates $a_{1}$ and $a_{2}$ result with different masses $m_{1}$ and $m_{2}$ and both coupled to $G\widetilde{G}$ with scales $f_{1}$ and $f_{2}$,

For each eigenstate, we can also compute exactly the factor $g_{i}$ in Eq. 2,

which describes how much the mass-decay constant relation deviates with respect to the single axion case, i.e. it measures the distance to the standard QCD band. For a single, canonical, QCD axion $g=1$.

Now we are ready to study two widely known limits in which one of the axion eigenstates ends up behaving like a single QCD axion, i.e. either $g_{1}=1$ or $g_{2}=1$. In the limit $r\to\infty$, one of the eigenstates, $a_{2}\simeq\hat{a}_{2}$, becomes infinitely massive, and the other one, $a_{1}\simeq\hat{a}_{1}$, reduces to the standard single QCD axion. In the opposite limit $r\to 0$, the determinant of $\mathbf{M}^{2}$ vanishes: one massive and standard QCD axion results $a_{2}\simeq(\hat{a}_{1}+\hat{a}_{2})/\sqrt{2}$ , while there is also a massless scalar but decoupled from the QCD sector $a_{1}\simeq(\hat{a}_{1}-\hat{a}_{2})/\sqrt{2}$. For intermediate $r$, however, both axion eigenstates have similar coupling to gluons and their mass-scale relation deviates from that for the standard QCD scenario (Eq. (1)).

This behaviour can be understood as given by the axionness of each eigenstate $\beta_{i}\equiv 1/g_{i}$, and noting that a sum rule is strictly obeyed:

This constraint illustrates how the QCD topological susceptibility is shared among the two axion eigenstates. It also implies that, whenever one of the two eigenstates decouples, the other converges towards the standard QCD axion. An intuitive interpretation of the $\beta_{i}$-factors can be found in App. B, where it is shown that for a physical axion to decouple from the sum rule it should either have a negligible PQ charge or a negligible projection onto the field (or combination of fields) coupling to $G\widetilde{G}$, ${a}_{G\widetilde{G}}$.

Examples of multi-axion solutions to the strong CP problem are illustrated in Fig. 2. If a experiment detects a signal at the location indicated by the star, and it corresponds to a pure QCD axion solution, it is predicted that

• •

  No signal should be found inside all the dashed area –delimitated by a grey line, including in particular the standard -yellow- QCD line.

• •

  A second axion signal should be found in the undashed area. For the $N=2$ and the toy model above, its location is marked in the figure by an orange diamond.

The first point stems from the sum rule in Eq. (13), and in fact it will hold whatever the total number of axion eigenstates, see next section. There, we will generalize that sum rule to an arbitrary number of axion eigenstates for generic PQ-invariant potentials.

In addition, Fig. 1 indicates that, for the eigenstate closest the the standard QCD band, the point of largest departure corresponds to $g_{1}=g_{2}=2$ (case not illustrated in Fig. 2); this fact will turn out to extend to arbitrary $N=2$ mixing potentials ³³3In the particular example above, the $r$ value which allows maximal departure of the eigenstate closest to the standard QCD band is $r=2$. Using the same simple potential, but for different $\hat{f}_{1}$ and $\hat{f}_{2}$, that maximal deviation is given by $r=(\hat{f}_{1}^{2}+\hat{f}_{2}^{2})/\hat{f}_{1}^{2}$..

The results above for the $\{m_{i},1/f_{i}\}$ QCD band apply directly to experiments which are independent of the axion UV complete models, such as CASPEr-electric JacksonKimball:2017elr or future experiments based on the piezoaxionic effect Arvanitaki:2021wjk . It is easy to prove, though, that the couplings to photons and nucleons get analogously displaced under some simplifying assumptions, see Section III.2 for this analysis in the generic $N$ axion scenario.

It will be shown next that whenever the axion field (or field combination) that couples to ${G\widetilde{G}}$ mixes with other $N-1$ scalar singlets via a PQ-invariant potential, it will lead in all generality to $N$ axion eigenstates, i.e. $N$ distinct massive pseudoscalars, each of them coupling to ${G\widetilde{G}}$ in a very specific way.

Let us first formulate the general case of a system of $N$ scalar fields $\{{\hat{a}}_{k}\}_{k=1}^{N}$ which couple to the QCD topological term through one arbitrary linear combination, and are in addition subject to an extra scalar potential $V_{B}({\hat{a}}_{k})$ which mixes them,

where, without loss of generality, the decay constants ${\hat{f}}_{k}$ have been redefined to absorb any coefficient in the coupling to gluons. Without the extra potential, the system would have exhibited $N$ independent global U(1) symmetries at the classical level, with one combination of the $N$ GBs becoming the standard massive QCD axion. The mixing potential reduces the symmetry. To be precise, note that neither the number of fields that couple to $G\widetilde{G}$ nor the dimension of the mixing potential $V_{B}$ needs to coincide with the number of scalars in Nature $N$ which mix, but the results below are independent of it, as shown in the basis of Eq. (19).

The condition to solve the strong CP problem in the presence of the potential is simply that the latter must exhibit a $U(1)_{PQ}$ symmetry, i.e. solely broken by QCD. Imposing this true PQ symmetry ensures that the minimum of the complete potential corresponds to the CP conserving point, $\bar{\theta}_{\text{eff}}=0$, and it will have important consequences on the properties of the different axions. From now on we will assume that all axion fields are already redefined and expanded around their minimum and thus we will drop the $\bar{\theta}$-term.

It is instructive to consider the complete scalar mass matrix which stems from Eq. (14). At energies below the QCD confinement scale, the gluonic coupling induces a non-zero potential⁴⁴4The leading order QCD axion potential DiVecchia:1980yfw ; Leutwyler:1992yt ; diCortona:2015ldu reads $V({\hat{a}}_{k})=-m_{\pi}^{2}f_{\pi}^{2}\sqrt{1-\beta\sin^{2}\left(\frac{1}{2}\sum\frac{{\hat{a}}_{k}}{{\hat{f}}_{k}}-\bar{\theta}\right)}$, where ${\beta\equiv\frac{4m_{u}m_{d}}{(m_{u}+m_{d})^{2}}}$. which corresponds to the following mass term,

where the QCD topological susceptibility $\chi_{\text{QCD}}$ was defined on the right-hand side (RHS) of Eq. (1). The complete mass matrix $\hat{\mathbf{M}}$ can then be defined as

where

has rank $1$ and encodes the QCD contribution to the masses of the axion eigenstates, while ${\big{(}\hat{{\mathbf{M}}}_{B}^{2}\big{)}_{kl}\equiv\partial^{2}V_{B}/\partial a_{k}\partial a_{l}}$ stems from the extra potential.

If the system has a PQ symmetry, the determinant of the complete mass matrix $\hat{{\mathbf{M}}}$ needs to vanish in the limit where the QCD topological susceptibility goes to zero, which in turn implies that the determinant of $\hat{{\mathbf{M}}}_{B}$ vanishes⁵⁵5In order for the result to apply both ways, one needs to add on the RHS the condition that the eigenvector of $\hat{\mathbf{M}}_{B}^{2}$ with vanishing eigenvalue, $a_{\rm PQ}$ (see App. B), has a non-zero projection to the linear combination that couples to gluons, i.e. $\langle a_{\rm PQ}|{a}_{G\widetilde{G}}\rangle\neq 0$. ,

A rotated basis. The problem can be equivalently formulated rotating to a field basis in which the combination of fields coupled to $G\widetilde{G}$ is redefined as a single field ${\hat{a}}_{G\widetilde{G}}$,

and as a consequence ${\hat{a}}_{G\widetilde{G}}$ mixes with the remaining $N-1$ fields via the (rotated) external potential $V^{\textrm{R}}_{B}$,

The relation of the scale $F$ to the $\{\hat{f}_{k}\}$ set in Eq. (18) follows from field normalization. In this perspective, the field coupling to ${G\widetilde{G}}$ is allowed to mix with other singlet scalars in Nature through a PQ-invariant potential. This is the generalization of the rotation performed earlier for $N=2$ between Eqs. (4) and Eqs. (5). The rotated basis will be preferred further below, as it renders particularly straightforward mathematical demonstrations.

Even though the maxion conditions require the presence of mass scales in the extra potential which are of the same order as the QCD-induced mass, and restrict the shape of the matrices, there are still infinitely many maxion matrices. Indeed, the class of matrices that verify these QCD Maxion conditions generate a $m$-parameter family of matrices with $m=N(N+1)/2$. Fixing all $\hat{f}_{i}=\hat{f}$, the number of families reduces to $m=N(N-1)/2+1$. Consequently, in the latter case all $N=2\,(3)$ maxion potentials are functions of $1\,(3)$ free parameter(s), plus the overall scale.

In a scenario with 2 axion fields produced with the same decay constant ($\hat{f}_{1}=\hat{f}_{2}=\hat{f}$), the only family of maxions is characterized by the following mass matrix:

where $r\in[0,1)\cup(1,2]$. For $r=2$, the toy example presented in Sec. II is recovered, while $r=1/5$ reproduces the example in App. A.

The eigenvalues obtained from this matrix lead to ${\Delta\equiv|m_{1}^{2}-m_{2}^{2}|/(m_{1}^{2}+m_{2}^{2})}$ ranging from $\sim 0.7-1$. In the limiting case in which $\Delta=1$, corresponding to $r=1$, the system is not PQ-invariant, as the massless eigenstate has no mixing with $a_{G\widetilde{G}}$; see Footnote 5. Nonetheless, values of $\Delta=1+\epsilon$ (with $\epsilon\ll 1$) are possible, corresponding to an almost massless maxion and another one with $m_{2}^{2}\approx 2\chi_{\rm QCD}/F^{2}$. This shows that highly hierarchical masses are possible. Assuming a different hierarchy for $\hat{f}_{1}$ and $\hat{f}_{2}$, values of $\Delta\sim 0$ can be attained in addition.

In scenarios with 3 axions, the maxion families are instead characterized by 3 free parameters, assuming equal decay constants for all fields. Identifying $\hat{a}_{1}$ with the PQ combination, an illustrative example is given by the mass matrix

with $r\in\left[\frac{1}{2}(1-\sqrt{13}),\frac{1}{2}(1+\sqrt{13})\right]$. The eigenvalues sourced by this mass matrix are shown in Fig. 4. If it was furthermore assumed that the extra source of interactions $V_{B}$ is diagonal in the scalar fields, it would follow that:

In this particular example, the entries of $\mathbf{M}^{2}$ in the preferred basis turn out to be simply the roots of Laguerre polynomials, as explained next.

Laguerre maxions. If the matrix $\hat{\mathbf{M}}_{B}^{2}$ in Eq. 16 is diagonal and all decay constants are equal $\hat{f}_{k}=\hat{f}$, one can show that the maxion conditions are automatically fulfilled if the characteristic polynomial of the mass matrix and its minor –in the preferred rotated basis– correspond to the Laguerre polynomials $\mbox{\large$p$}_{{\,{\mathbf{M}}^{2}}}(\lambda)=L_{N}(\lambda)$ and $\mbox{\large$p$}_{{\,{\mathbf{M}}_{1}^{2}}}(\lambda)=L_{N-1}(\lambda)$. We call the solutions produced in this scenario Laguerre maxions. Indeed, one can check that for the zeroes $\lambda_{i}$ of the $N$-th order Laguerre polynomial $L_{N}(\lambda)$, the following identity holds⁹⁹9This can be shown using the expression of $g_{i}$ in Eq. 90 and the property $\lambda\,dL_{N}(\lambda)/d\lambda=NL_{N}(\lambda)-NL_{N-1}(\lambda)$. ,

ensuring that all physical axions depart from the standard QCD axion line by a factor of $\sqrt{N}$. The Laguerre Maxion scenario is particularly predictive since the spread in the values of the QCD axion masses will be fixed and determined by the zeroes of the Laguerre polynomials (see example the cases of $N=3,30$ in Fig. 4).

The main experimental impact of our results is that multiple axion signals –instead of a single one– are expected from the generic QCD axion solution to the strong CP problem. The customary solution with a single axion QCD corresponds instead to a certain limit of the general parameter space. These multiple solutions can only be displaced towards the right of the customary QCD axion band.

In case of a positive ALP signal to the right of that band (or even within the single QCD band given its error bars), experiments are encouraged to look for similar signals within their frequency range. In addition, signals widely separated and thus accessed by different experiments are possible. This establishes a beautiful synergy between experiments sensitive to vastly different masses and scales, as the complete reconstruction of the solution to the strong CP problem may require this complementary search.

An illustration of this experimental impact on the $\{m_{a},1/f_{a}\}$ parameter space of axion-gluon couplings (which are directly tackled by axion-nEDM experiments) is provided in Fig. 2. Specific $N=2$ and $N=3$ QCD axion solutions are depicted, corresponding to mixing matrices in Eqs. 7 and 66 and a generic 3-axion scenario. A star indicates the putative location of the first signal detected in the axion/ALP parameter space. For it to be linked to a pure axion QCD solution of the strong CP problem, it follows that:

• •

  No other axion signal can be found anywhere in the greyed area (which includes the single QCD axion band). This holds independently of the total number of axion eigenstates $N$.

• •

  For $N=2$ only a second signal should appear exactly in the grey line, so as to saturate the sum rule; an example is indicated by an orange diamond.

• •

  For $N=3$, two other signals await discovery in the undashed area, see e.g. the clear blue diamonds.

• •

  An $N=3$ maxion solution is also possible in the example chosen, signaled by the dark blue points, as the distance of the star signal to the single QCD axion band is $\sqrt{3}$.

We have shown that the same pattern holds for arbitrary $N$, for which we have determined in all generality the condition for an arbitrary scalar potential to be PQ invariant. Furthermore, a precise sum rule for the QCD axion has been found, together with other exact results.

For the axion-photon couplings, we have illustrated the analogous multiple signals and displacements in Fig. 3, assuming universal $E/{\mathcal{N}}$ couplings and focusing on maxion signals: both a $N=3$ and a $N=30$ maxion multiplet are shown, for the mixing provided by the Laguerre matrices described in Section IV. The displacements of axion couplings to elementary fermions will be subject to model-dependences analogous to those of axion-photon couplings, to be further explored elsewhere future .

In practice, for signals close to the standard QCD band the possible displacement of a given signal will be obscured by the incertitude on the width of the band itself. For experiments measuring directly the coupling to gluons, there is a large uncertainty on the theoretical prediction of $n$EDM$(f_{a})$ JacksonKimball:2017elr . For axion-photon searches, the experimental measurements are very precise with uncertainties stemming instead from the UV model dependence, i.e. the different matter content of axion models and different $E_{k}/{\mathcal{N}}_{k}$ factors. In summary, the smoking gun may turn out to be the multiplicity of signals, as experiments have reached impressive precision on frequencies.

Furthermore, most of the experiments with best reach rely on the assumption that axions account for the dark matter (DM) of the universe. Haloscopes are sensitive to $\sqrt{\rho_{{\rm DM},local}}\times g_{aXX}$, where $\rho_{{\rm DM},local}$ denotes the local DM density and $g_{aXX}$ a generic axion coupling to the visible world. Uncertainties stem from $\rho_{{\rm DM},local}$ (which could be largely modified by $e.g.$ miniclusters, focusing effects, etc.) and, for the multiple QCD axion under discussion, on how DM is distributed among them. If democratically distributed among the $N$ axion eigenstates, the flux of each species would diminish by a factor of $N$ and consequently the bounds would be weakened. We have illustrated this weakening of the projected bounds for CASPEr-electric in Fig. 2, for the $N=3$ scenario. Nevertheless, the thermal evolution of the Universe with a multiple QCD axion may lead to very different spectra and evolution scenarios (e.g. Kitajima:2014xla ; Ho:2018qur ; Cyncynates:2021xzw ; Adams:2022pbo ), which deserve a dedicated future study.

The quest to unravel the fundamental symmetries of the visible and dark sectors of the Universe is inspiring. We have proven that the PQ solution to the strong CP problem within pure QCD leads in all generality to multiple axion signals with well delimitated properties. This stems from lifting the requirement that the basis of axion-gluon interactions and the axion mass basis are simultaneously diagonal. In this perspective, the usual single QCD axion is just one particular case of the landscape of solutions, corresponding to zero mixing. The results open novel experimental territory and blur the distinction between the search for the QCD axion and that for ALPs in a huge region of the parameter space.

Acknowledgments.—We are indebted to Daniel Alvarez-Gavela for useful discussions and input. We also thank Aneesh Manohar and Benjamin Grinstein for interesting discussions. This project has received funding /support from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska -Curie grant agreement No 860881-HIDDeN, and under the Marie Sklodowska-Curie Staff Exchange grant agreement No 101086085-ASYMMETRY. The work of M. R. is supported by the Marie Sklodowska -Curie grant agreement No 860881-HIDDeN. The work of P.Q. is supported in part by the U.S. Department of Energy Grant No. DE-SC0009919. B. G. acknowledge as well partial financial support from the Spanish Research Agency (Agencia Estatal de Investigación) through the grant IFT Centro de Excelencia Severo Ochoa No CEX2020-001007-S and through the grant PID2019-108892RB-I00 funded by MCIN/AEI/ 10.13039/501100011033. B.G. thanks very much the Particle Physics group of the University of California San Diego, where part of this work was carried out. Likewise, P. Q. thanks the Galileo-Galilei Institute for theoretical Physics in Florence (GGI) for their warm hospitality, where part of this work was carried out. We are also thankful to Anna Lewis for helping to improve the writing of the manuscript.