Chiral models of composite axions and accidental Peccei-Quinn symmetry

We introduce a class of composite axion models that provide a natural solution to the strong CP problem, and possibly account for the observed dark matter abundance. The QCD axion arises as a composite Nambu-Goldstone boson (NGB) from the dynamics of a chiral gauge theory with a strongly-interacting and confining SU(N) factor and a weakly-interacting U(1), with no fundamental scalar fields. The Peccei-Quinn (PQ) symmetry is accidental and all the mass scales are generated dynamically. We analyze specific models where the PQ symmetry is broken only by operators of dimension 12 or higher. We also classify several other models where the PQ symmetry can be potentially protected up to the dimension 15 or 18 level. Our framework can be easily extended to a scenario where the Standard Model (SM) is unified into a simple gauge group, and we discuss the case of non-supersymmetric SU(5) unification. The GUT models predict the existence of additional pseudo NGBs, parametrically lighter than the GUT and PQ scales, which could have an impact on the cosmological evolution and leave observable signatures. We also clarify the selection rules under which higher-dimensional PQ-violating operators can generate a potential for the axion in the IR, and provide a discussion of the discrete symmetries in composite axion models associated to the number of domain walls. These results can be of general interest for composite axion models based on a QCD-like confining gauge group.


Introduction
The strong CP problem can be elegantly solved by the QCD axion [1][2][3][4], which assumes the existence of a spontaneously broken global symmetry -known as the Peccei-Quinn (PQ) symmetry -anomalous under SU(3) c . Crucially, QCD non-perturbative effects generate a potential for the axion such that the field at the minimum of its potential cancels the θ dependence exactly. In most of the axion models studied in the literature, the PQ symmetry is imposed by hand, leaving its UV origin unspecified.
In this work we introduce a class of models where the PQ symmetry arises as an accidental global invariance of a strongly-coupled chiral gauge theory, with no fundamental scalar fields in the PQ sector. The PQ symmetry is robustly protected by gauge invariance and the GUT dynamics does not induce PQ-breaking operators. The idea of an accidental PQ symmetry was first proposed in Ref. [5], while that of a composite axion was put forward by Kim [6] and then further studied in Refs. [7,8]. The setup we consider can be seen as an extension of Kim's model where however mass terms are forbidden by a chiral U(1) gauge symmetry, so that the PQ invariance is accidental. The UV dark sector is similar to that of the Composite Dark Matter models introduced in [9][10][11]; the crucial difference lies in the use of accidental symmetries to protect the PQ symmetry rather than species number. Moreover, the structural absence of mass terms automatically guarantees the vanishing of the CP-violating theta angle for the new strong dynamics.
The axion solution is notoriously fragile: additional UV contributions to the axion potential need to be extremely suppressed in order not to shift its minimum and spoil the cancellation of θ. Simple arguments suggest that ultimately all global symmetries are explicitly broken by quantum gravity at the non-perturbative level. Planck suppressed operators with order one coefficients generally induce too large a PQ breaking, giving rise to the so-called axion quality problem [12][13][14][15]. However, it is still unclear whether nonperturbative quantum gravity effects will be exponentially suppressed by a factor e −cM P l /fa , where c is an O(1) coefficient, as suggested by some semiclassical computations. In the latter case, models with an accidental Peccei-Quinn symmetry can provide a fully satisfactory solution to the strong CP problem for axion decay constants f a ≤ 1.2 · 10 16 GeV -see Ref. [16] for a review. Some of the models analyzed in this work, in particular the minimal models of Section 4, naturally have an accidental PQ symmetry but do not solve the quality problem. They have a simple structure and could be relevant if the gravitational UV completion is weakly coupled and the dominant source of PQ breaking besides QCD comes from exponentially suppressed gravitational effects. In Sections 4 and 5 we present models with a slightly more elaborate structure that can also provide a natural structural solution to the axion quality problem. We analyze one particular such construction in Section 5 in which the PQ symmetry is accidental up to dimension 12 operators.
Attempts to construct accidental composite axions solving the quality problem have been pioneered by Randall [17], who put forward a chiral model based on a new product gauge group SU (n) × SU (m), where SU (n) is assumed to get strong in the IR and confine, while SU (m) remains weakly coupled and in the Higgs phase. The parameter space is significantly constrained by the requirement of having a confining SU (n) dynamics and perturbative SM gauge couplings up to the Planck scale [18]. Additional models based on non-abelian product gauge groups, which generalize Randall's construction and improve its quality, have been put forward in Refs. [20,21]. Supersymmetric composite axion models, also based on product gauge groups, have been instead constructed in Refs. [22,23]. Other attempts to build composite axions and solve the quality problem by means of chiral gauge theories include the models of Ref. [24] and Ref. [25]. The first makes use of a chiral and confining simple gauge group, while the second exploits an axial U(1) gauge group and is more similar to our approach. Holographic constructions dual to composite axion models have been analysed in Refs. [26][27][28][29][30][31]. Finally, several alternative approaches to the quality problem where the axion is a fundamental field have been proposed in the literature; for an incomplete list of works, see for example Refs. [32][33][34][35][36][37][38][39][40][41][42].
In a string theoretic context, QCD axions can also arise through dimensional reduction of higher-form fields. Such axions are protected by an exact shift symmetry at the perturbative level, but the PQ symmetry is expected to be broken by non-perturbative effects such as non-QCD instantons. The size of these corrections depends on the value of the action of the relevant instanton configuration, and there are models in which they are small enough to have a viable QCD axion. This requirement, together with that of an axion decay constant parametrically smaller than the Planck scale, necessitates non generic UV features; nonetheless these models can provide fully consistent UV solutions to the axion quality problem [43,44] (see also the recent work [45]). Our approach differs in that the solution to the quality problem is obtained as an IR property of the low-energy theory, the PQ symmetry being robustly protected by gauge invariance irrespectively of the details of the dynamics taking place at the Planck scale, provided that the theory admits a quantum gravity UV completion.
On a different note, the possible connection between the axion and an SU(5) Grand Unification dynamics has been first proposed in [46,47] and recently re-examined in Refs. [48][49][50][51][52]. All these models differ from the SU(5) model here described in two aspects: the Peccei-Quinn symmetry is not accidental and the dynamics responsible for its breaking is the GUT dynamics itself. We shall see instead that our setup makes it easy to realize the PQ symmetry as an accidental invariance of the GUT theory. Along the same lines we follow, Refs. [20,21] considered SU(5) unified models with accidental composite axions. Additional attempts to build an accidental Peccei-Quinn invariance in theories of Grand Unification can be found in Refs. [5,53]. For an extensive review of QCD axion models we refer the reader to Ref. [54].
This article is organized as follows. The broad class of theories that will be considered in our analysis is introduced in Section 2. In Section 3, we study PQ violation from higher-dimensional operators in a general setting, and formulate a criterion to identify those operators which can affect the axion potential in the IR. Models with n f = 2 and n f = 3 dark flavors are classified and analyzed in Sections 4 and 5 respectively, and their phenomenology is discussed in Section 6. Appendix A provides a general analysis of discrete symmetries and the number of domain walls in composite axion models based on a vector-like confining gauge group. Appendix B presents a simple QCD model with n f = 3 in detail, while Appendix C contains a comprehensive list of charge assignments for n f = 3 yielding a high-quality axion. A comparison with Randall's model is also included in Appendix D.

Model building
We consider theories where the only UV degrees of freedom in the Peccei-Quinn sector are fermions and gauge fields 1 . The gauge group is assumed to be of the form G = G D × G SM , where G SM is either the SM gauge group or one of its grand unified counterparts, and G D contains a confining subgroup that generates the Peccei-Quinn scale dynamically, giving rise to the composite axion. The PQ sector fulfills the following conditions: i) it entails good theoretical control over the IR dynamics and its global symmetry breaking pattern; ii) it leaves the SM gauge symmetries unbroken; iii) the representations of the new fermions under G are chiral. In order to have control over the IR dynamics in the PQ sector, we assume that G D includes a confining dark color SU(N DC ) subgroup under which the new fermions transform as fundamental plus anti-fundamental vectorlike representations. We take N DC ≥ 3, so that the fundamental representation of SU(N DC ) is complex. To ensure that the SM gauge symmetries are not broken spontaneously by SU(N DC ) confinement we require that the new fermions transform also as vectorlike representations of G SM . The chirality of the model then rests on the existence of an additional subgroup G W ⊂ G D whose dynamics is perturbative at the SU(N DC ) confinement scale. The case with G W = SU(m) was first discussed in Ref. [17] and then further analysed and generalised in Refs. [18,20].
In this work we focus on the case in which the weak group is abelian, G W = U(1) D .
A general class of models fulfilling the conditions discussed above can be defined in terms of two sets of left-handed Weyl fermions transforming as Therefore, n f is the number of different U(1) D charges assigned to the ψ (or equivalently to the χ) fields, and the SM representations r i are in general reducible. The U(1) D charges are chosen so that the representations are overall chiral. This requirement and the cancellation of U(1) D gauge anomalies imply that the multiplicity of fermions n f has to be at least 2: n f ≥ 2.
1 The SM sector can include scalars, such as the Higgs field or the ordinary heavy scalars of a grand unified UV completion.

-4 -
The most general renormalizable Lagrangian built from this field content comprises only kinetic terms, including a possible kinetic mixing between the hypercharge gauge boson and the U(1) D gauge field: Here G a µν , F D µν and B µν are the SU(N DC ), U(1) D and SM hypercharge field strengths respectively. No gauge invariant mass term or Yukawa coupling can be written due to the chiral structure of the model, and one can therefore set the SU(N DC ) theta term to zero through an axial field redefinition. The label κ runs over all irreducible representations of SU(N DC ) × G SM × U(1) D , whose number we will denote by n irr . 2 The Lagrangian of Eq. (2.2) has a [U(1) V ] n irr × [U(1) A ] n irr accidental global symmetry at the classical level. 3 The n irr vectorial U(1) V factors are linearly realized and each of them implies an accidentally conserved charge. One of the n irr axial U(1) A factors is anomalous under SU(N DC ), and one is gauged by U(1) D , whose vector boson acquires a mass through the Higgs phenomenon; the others are spontaneously broken by the SU(N DC ) condensates ψ κ χ κ . The theory thus contains (n irr − 2) Nambu-Goldstone bosons (NGBs), one of which is identified with the axion. A necessary condition in order to have a composite axion is thus n irr ≥ 3. The NGBs can acquire mass through higher-dimensional operators or from the QCD dynamics if they have anomalous couplings to the gluon field.
If the U(1) D ×G SM weak gauging is switched off, the global invariance of the Lagrangian Below the SU(N DC ) confinement scale, the global symmetry group is spontaneously broken to the diagonal subgroup, and the low-energy theory is described by a set of (N 2 f −1) NGBs. By turning on U(1) D × G SM , one of these NGBs becomes the longitudinal polarization of the massive U(1) D gauge boson, and most of the other ones become pseudo and acquire a potential through perturbative loop effects. The only light degrees of freedom are the (n irr − 2) NGBs coming from the accidental symmetry discussed above.
This same model building scheme was thoroughly analyzed for n f = 2 in Ref. [11], where the construction of Refs. [9,10] was generalized to new fermions charged under the SM gauge group. In the minimal setup of these works the choice of representations was such that n irr = 2 and no light NGBs or axions were present in the infrared. The Dark Matter (DM) candidate was identified with the lightest pseudo NGB charged under an accidental species number, and the dynamical confinement scale was assumed to be of order 1−100 TeV, much smaller than a typical Peccei-Quinn scale. We shall discuss further 2 One has nirr = n f if all the ri are irreducible representations. 3 The accidental symmetry is bigger than that if any of the irreducible representations inside a given ri has a non-trivial multiplicity. the precise relation between this work and Ref. [11] in Section 4, where we classify and investigate models with n f = 2.
In the models studied in this work, the axion decay constant f a is determined in terms of the dark color confinement scale Λ DC . The observation of the Supernova 1987A, together with the sharp prediction of the low-energy axion-nucleon couplings that follows in our models (see Eq. (6.1) and the discussion in Section 6), provides the most stringent lower bound on f a , independently of its role as a cosmological relic: f a 4 · 10 8 GeV [55] (see however Ref. [56]). It is well known that the axion can naturally reproduce the observed Dark Matter abundance through the misalignment mechanism for values of its decay constant f a ∼ 10 12 GeV (for a recent numerical analysis see Refs. [57,58]). Motivated by these considerations, we will consider f a ∼ 10 9 GeV and f a ∼ 10 12 GeV as benchmark values when analyzing our models.

IR confinement and UV perturbativity
Our construction relies on the assumption that the SU(N DC ) dynamics confines in the infrared. Since the new fermions have no mass term, SU(N DC ) must be asymptotically free; at one loop this gives the bound N DC ≥ 2N f /11. A stronger constraint N f < N conf is obtained by considering the lower boundary of the conformal window. The value of N conf is a non-perturbative property that cannot be derived with perturbative computations.
Lattice simulations suggest that N conf ≈ 12 for an SU(3) gauge theory with fermions in the fundamental representation, see Ref. [59]. In selecting our models we will use a naive linear extrapolation of this result to generic numbers of dark colors (see for instance Refs. [60][61][62][63] for some evidence in support of this assumption) and require where y i,α is the hypercharge of the representation r (α) i and b (1) SM = −41/6. In the following, we will focus on models satisfying the conditions of Eqs. (2.3), (2.5) and (2.6).

Peccei-Quinn symmetry and axion quality
As highlighted in Section 1, the existence of higher-dimensional PQ-violating operators poses a significant challenge to the axion solution of the strong CP problem. A simple estimate gives an idea of how serious this problem can be. If the PQ breaking effect comes from a single insertion of a local operator generated at the Planck scale, one expects a shift of the value of θ at the minimum of the potential of order where ∆P Q is the dimension of the operator, c is a dimensionless effective coefficient and φ CP is a CP -violating phase. The current experimental bound on the electric dipole moment (EDM) of the neutron [65] requires ∆θ 10 −10 , which translates into a lower bound on the dimension of the PQ breaking operator Assuming |c|, φ CP ∼ 1 and using M Pl ∼ 10 19 GeV, Λ QCD ∼ 100 MeV, this gives ∆P Q 9 for f a ∼ 10 9 GeV. Values of the PQ scale of order f a ∼ 10 12 GeV, which allow the axion to reproduce the observed dark matter abundance naturally through the misalignment mechanism, lead to a stronger bound ∆P Q 13. In theories of composite axion, for 5 More precisely, one should integrate out the top quark at its mass scale mt and use b (3) SM = 7 + 2/3 to run from mZ to mt. This leads to a small correction to Eq. (2.5) that we neglect for simplicity. a purely fermionic PQ-violating operator generated through diagrams with loops, one expects the effective coefficient to scale as c ∼ g 2 U V (4π/g U V ) ∆P Q /3−2 , where g U V is a UV coupling. This corresponds to values of c larger than 1, hence to a stronger constraint on ∆P Q , unless > ∆P Q /6 − 1 and g U V is sufficiently small. 6 Such strong bounds on the dimension ∆P Q make it clear that obtaining a natural resolution of the axion quality problem is not an easy task.
The problem is even more serious if PQ-violating operators are generated at energies lower than the Planck scale by perturbative dynamics. In this case the minimum values of ∆P Q required to preserve the axion solution are larger. Theories of Grand Unification where the PQ symmetry does not commute with the GUT group are an example of this situation.
For GUT scales of order 10 17 GeV one in general needs ∆P Q 10 for f a ∼ 10 9 GeV, and ∆P Q 16 for f a ∼ 10 12 GeV. Compatibility of the axion solution with Grand Unification thus imposes non trivial constraints on the structure of the theory.
In our analysis of candidate models with accidental PQ symmetry below the Planck scale, we will take the minimum dimensionalities ∆P Q = 9 and ∆P Q = 12 as targets for a solution of the axion quality problem for f a = 10 9 GeV and f a = 10 12 GeV respectively.
We choose these values considering the high sensitivity of Eq. Throughout this work we will compute the dimensionality of operators by assuming that the theory is weakly coupled at the UV scale and that the dark color force gets strong near the confinement scale Λ DC . In this case, the scaling dimension of operators is equal to the classical dimension up to small corrections. On the other hand, if the dark color dynamics has a walking (as opposed to running) behavior, i.e. it remains strongly coupled all the way up to the UV cutoff scale or slightly below, then the RG evolution above Λ DC is characterized by scaling dimensions that might be significantly larger than the classical ones, thus ameliorating or even erasing the axion quality problem. A dynamics of this kind can for example arise if N DC is such that dark color confines but is close to the conformal window. In this case the RG evolution between the UV scale, where the PQ breaking effects are generated, and the infrared confined phase below Λ DC is characterized by a strongly-coupled flow lingering close to a pair of fixed points at complex coupling [66,67].
While the scaling dimension of PQ-violating operators cannot be computed analytically at strong coupling, and the predictivity of our theoretical construction is lost in this limit, 7 6 The more accurate estimate of ∆θ in composite axion theories performed by Ref. [21] suggests that the actual constraint might be even stronger due to the presence of very large multiplicity factors. These are smaller in our models compared to those of Ref. [21] though still quite large. For ∆P Q = 12, including such factors leads to an upper bound on fa stronger by roughly one order of magnitude. 7 One could however use supersymmetry to gain analytic control over the strong dynamics, see for example a scenario of this kind could be relevant for the resolution of the quality problem and is worth being investigated. Interestingly, realizations of this idea in the context of warped extra dimensions could provide calculability and lead to a concrete candidate model, see for example Ref. [28].

General remarks on PQ violation: a spurion analysis
As already noticed by Ref. [18], not all PQ-violating operators can generate a potential for the axion at low energy. Since identifying the dangerous operators is obviously crucial to single out viable theories, in this section we perform a detailed analysis based on selection rules and spurion arguments. Our results apply to any SU(N DC ) × G W × G SM gauge theory with Weyl fermions ψ and χ transforming as vectorlike representations of SU(N DC ) × G SM .
A generic PQ-violating operator OP Q induces a potential for the axion if the matrix element ψ a,ϕ |OP Q |0 is non vanishing. Here |ψ a,ϕ is an arbitrary state containing axions and possibly other SM-neutral NGBs, denoted by ϕ, with zero momentum. 8 Therefore, OP Q is an interpolating operator for the state |ψ a,ϕ and must have its same quantum numbers.
In order to determine the form of OP Q implied by this requirement, it is useful to first neglect the effects of the weak gauging of G W × G SM . In this limit, the global symmetry

under which the fermions transform as
Any G F -violating operator made of the ψ and χ fields can thus be characterized by a spurionic tensorial structure Ref. [40] for a model of this kind. 8 While we are ultimately interested in a potential for just the axion field, we conservatively consider the possibility of generating a mixed potential with other SM-neutral NGBs, which can in turn lead to an axion potential at energies smaller than the mass of the ϕ's.
right indices, as well as an equal number of lower left and upper right indices. In other words, the spurionic structure of an interpolating operator in the low-energy theory must be of the form T {i}{j} {j}{i} ; clearly, the interpolating operator in the UV theory must have the same spurionic structure. Furthermore, any such operator must have zero U(1) V charge, i.e. be a 'mesonic' operator. This implies that OP Q is necessarily a polynomial of the SU(N DC )-singlet bilinears (ψ r χr) and (ψ r χr) † (whose flavor indices are left understood), where ψ r and χ r denote fields transforming as the irreducible SM representation r. Such fermionic bilinears are exactly those acquiring a vacuum expectation value at the SU(N DC ) confinement scale. The same conclusion was obtained by Dobrescu in Ref. [18].
Let us now consider the effects of the weak gauging of G W × G SM . The latter breaks Operators satisfying this requirement are polynomials of the bilinears (ψ r χr), (ψ r χr) † , (ψ † r ψ r ) and (χ † r χr). Hence, operators built of (ψ † r ψ r ) and (χ † r χr) factors are also allowed as a consequence of the weak gauging. 9 The axion potential in that case will be suppressed by powers of the weak gauge couplings.
It is instructive to analyze some examples in the more familiar case of QCD, where a similar selection rule characterizes the operators that can generate a potential for the pions. Consider for example the local operator O = (ud c )(su c ) in three-flavor QCD where electromagnetism plays the role of the weak gauging. 10 In the limit of vanishing up and down quark masses, π 0 is an exact NGB at the renormalizable level, and it is interesting to ask whether it gets a potential and a mass from insertions of O. The latter explicitly breaks the U (1) 3A axial symmetry, but has non zero isospin (i.e. charge under the vectorial U (1) 3V ) and thus cannot generate any potential for π 0 . More in detail, one can see that in the low-energy theory O gives rise to non-derivative interactions among π 0 and the other (pseudo) NGBs; these interactions violate isospin by +1/2 units. Since the rest of the dynamics conserves isospin, there are no diagrams at any loop order in chiral perturbation 9 Notice that (ψ † r ψr) and (χ † r χr) transform as vectors under the Lorentz group. If in a given local operator they appear with gauge indices contracted to form a singlet of G SM , then a new operator can be constructed with lower dimensionality where each bilinear is replaced by a covariant derivative. 10 Here q c denotes the complex conjugate of the right-handed quark q. Therefore, q c is a left-handed field transforming as an anti-fundamental under color.
theory that generate a potential for π 0 . Another instructive example is that of the operator O = (uu c )(u †σµ u) 2 in two-flavor QCD with vanishing quark masses plus electromagnetism.
is generated at leading order by diagrams with two loops of charged pions. In the limit in which electromagnetism is switched off, the charged pion becomes massless and these loops vanish identically (one can for example use dimensional regularization to preserve chiral invariance). Therefore, the operator O generates a potential only in presence of the electromagnetic (weak) gauging; this is expected since it contains factors of (u † u).
So far we have implicitly assumed that the axion potential is generated by single This is smaller than the dimension of the operator built out of their product by 4(N − 1).
It is then possible for multiple insertions to give the dominant contribution to the axion potential. Indeed, while the product of local operators will have to satisfy the selection rule discussed above and thus be a polynomial of (ψ r χr), (ψ r χr) † , (ψ † r ψ r ) and (χ † r χr), the individual operators need only to violate flavor (otherwise they play no role) and at least one of them needs to break the PQ symmetry. Notice that G SM invariance of each individual operator can be achieved by adding SM or GUT fields. In sec. 5.2 we will discuss one example of a theory where the existence of dimension-6 operators prevents the resolution of the quality problem despite their single insertion cannot generate an axion potential.
From these considerations it should be clear that a complete classification of the higherdimensional operators capable of spoiling the axion quality is not an easy task, even for a specific model. One would like to determine, for each given model, the smallest dimension ∆P Q at which an axion potential is generated. This can be defined as the minimum over (1) the dimensions of all gauge-invariant, PQ-violating local operators that are polynomials of (ψ r χr), (ψ r χr) † , (ψ † r ψ r ), (χ † r χr); and (2) the effective dimensions of all gauge-invariant, PQ-violating polynomials of (ψ r χr), (ψ r χr) † , (ψ † r ψ r ), (χ † r χr) that can be formed as the product of gauge-invariant, flavorviolating local operators.  Non-local operators of type (2) arise from multiple insertions of local operators. We will thus consider a model able to resolve the axion quality problem for some given physically- allowed value of f a if the smallest PQ-violating dimension ∆P Q defined as above satisfies the bound of Eq. (3.2). In particular, we will require ∆P Q ≥ 9 for f a ∼ 10 9 GeV and ∆P Q ≥ 12 for f a ∼ 10 12 GeV.
Notice that multiple insertions of operators with dimension bigger than (d + 4)/2 are always more suppressed than a single insertion of a dimension d operator. This means that non-local operators of type (2) are not relevant to determine ∆P Q if they have dimension larger than (∆P Q +4)/2. The calculation of ∆P Q is therefore simplified if multiple insertions can be neglected based on the minimum dimensionality of flavor-violating operators. The number of n f = 2 theories of this kind is for example reported in Tables 7 and 8. 4 Models with n f = 2 Models with n f = 2 and r i irreducible have been studied and classified in Ref. [11]. With reference to the general structure of Table 1, the most general charge assignments that are vector-like under SU(N DC ) × G SM and overall chiral in the case of irreducible r i are as with r 1 =r 2 ≡ r and q ∈ [0, 1) (Type III). The only phenomenologically viable models consistent with values of the confinement scale Λ DC ∼ 1−100 TeV were identified in Ref. [11] to be those with SU Theories with reducible r i were found to be strongly constrained by the requirement of perturbativity. In practice, the only non-trivial candidate was a Type-II theory with r equal to the direct sum of one electroweak doublet plus one color triplet; by combining these two irreps into a fundamental representation of a unified SU(5) group one could obtain a GUT model. However, such a theory turns out to be not phenomenologically viable in the minimal scenario considered in Ref. [11]. Indeed, it contains two light pseudo NGBs, one of which is an axion-like particle. The corresponding PQ symmetry does not commute with SU (5)  When raising the dynamical scale to values above ∼ 10 8 GeV, however, new choices of reducible representations r i become allowed by perturbativity, and it is possible to construct realistic axion models. This requires having n irr ≥ 3, as discussed in Section 2, and one way to obtain this is to consider n f = 2 with r i reducible. The simplest possibility is to start from the minimal model of Ref. [11] with QCD triplets and enlarge each representation by the addition of one SM singlet. We follow this route and present this minimal QCD model first and then its unified version. These models turn out to have ∆P Q = 6 and therefore do not solve the quality problem, but it is nonetheless very instructive to analyze them. We then generalize this construction in Section 4.3, where we classify all possible models with n f = 2 and generic reducible representations. Our analysis shows that these additional models have a chance of solving the quality problem.

A Minimal QCD model
We consider a Type-II model with the charge assignment detailed in Table 2 and r = For this choice, the SM triplets can be assigned a non-zero hypercharge y, which modifies the low-energy axion-photon coupling as we shall see. Type I and Type III models can be also constructed for y = 0 and lead to the same low-energy phenomenology. 11 They differ from Type-II models only for the spectrum of heavy pseudo NGBs and other resonances with mass of order Λ DC , which appear to be inaccessible with current and near-future experiments. For any rational value of the parameter q in the interval (−1, 1), the dark fermions have chiral representations and no gauge invariant mass term or Yukawa coupling can be written. 12 At the classical level the model has a U(1) 4 L × U(1) 4 R global symmetry, corresponding to the group of phase transformations of each of the fields ψ i and χ c i (χ c is the conjugate of χ and a right-handed Weyl fermion). At confinement the pattern of symmetry breaking giving four spontaneously broken U(1) axial factors. One of them is anomalous under SU(N DC ) and one is gauged by U(1) D . We are thus left with two broken generators, associated with two Nambu-Goldstone bosons. These are neutral under the unbroken vectorial U(1) V symmetries, and will be denoted by s and a in the following.
, s and a are interpolated by the axial currents The current J µ a has an anomaly in the background of SM gauge fields: where G µν and B µν are the QCD gluon and hypercharge field strengths respectively, the subscripts 3 and Y refer to the QCD and hypercharge gauge couplings, and dual field strengths are defined asF µν = µνρσ F ρσ /2. From now on we suppress indices, writing FF in place of F a µνF a,µν . The anomaly coefficients are where t i and Y are the QCD and hypercharge generators respectively. We can thus identify U(1) a with the Peccei-Quinn symmetry and the corresponding NGB a will play the role of the QCD axion. The other spontaneously broken accidental symmetry, U(1) s , is free from gauge anomalies and its NGB s is exact, hence massless, at the level of the renormalizable Lagrangian. The appearance of such massless state in the IR can be also deduced from an 't Hooft anomaly matching argument. 13 At the level of discrete symmetries, one can define the action of charge conjugation C and parity P on the new fields as follows 14 C : where G a µ denotes the SU(N DC ) gauge fields, T a are the corresponding generators, and η µµ is equal to +1 for µ = 0 and −1 for µ = 1, 2, 3. Both C and P defined in this way are explicitly broken by the chiral U(1) D gauging, but CP is preserved thanks to the absence of a dark color theta angle. It is also possible to define an alternative charge conjugation symmetry, C , which is exact in the new sector and broken only by the EW interactions of SM fields. We define C : whereas the SU(N DC ) gauge fields do not transform. Under C the NGBs s and a have charge −1 and +1 respectively, while they are both odd under the action of CP .
After dark confinement the U(1) D gauge boson acquires a mass The pseudo NGBs, whose corresponding global symmetries are explicitly broken by the weak gauging, acquire a mass from loops of SM and/or U(1) D gauge bosons of order At energies below the PQ scale one can integrate out all these heavy states and write an effective theory describing the relevant degrees of freedom, which are the NGBs a and s, and the SM fields. Since the generators Q a and Q s commute, the effective theory is simply that of two shift-symmetric scalar fields, plus terms capturing the anomalies and SM interactions. At quadratic order in the derivatives, and defining the axion decay constant f a as customary in axion physics, the effective lagrangian at the PQ matching scale is: From the latter equation we estimate where N DC = 3 is at the edge of the conformal window [59]. Furthermore, Eq. (2.6) gives the bound y 2 0.1 for N DC in the range of Eq. (4.11).
As we discuss in detail in Appendix A, in composite axion models based on a confining vector-like group the number of distinct vacua, also known as the domain wall number N DW , is equal to the anomaly coefficient computed in the basis where all the PQ charges are coprime integers. In our notation, The UV contribution to the coupling of the axion to photons can be characterized, as usual, by the E/N ratio and its value is controlled in our model by the arbitrary hypercharge y; we can also generate ∂ µ aj µ , as well as ∂ µ sj µ , at the matching scale if present in the UV theory.
Since the Peccei-Quinn symmetry -identified with U(1) a -is accidental, it is important to determine at which level it can be broken by higher-dimensional operators. The lowest dimensional PQ-breaking operators consistent with Lorentz and gauge invariance have scaling dimension 6 and are of the form ψ 1 ψ 2 χ 1 χ 2 and ψ 3 ψ 4 χ 3 χ 4 . These operators are gauge singlets for arbitrary values of q and carry a non-vanishing Peccei-Quinn charge.
Moreover, they have the appropriate flavour structure to generate a potential for the axion, according to the criterion of Section 3.1.

A Minimal SU (5) GUT model
The existence of dimension-6 PQ-breaking operators should prompt one to worry that additional UV dynamics, in particular Grand Unification, may spoil the axion solution.
We shall show that it is actually easy to embed the QCD model in a GUT theory in such a way that U(1) a is an accidental symmetry of the unified dynamics below the Planck scale.
As long as gravitational effects are exponentially suppressed, the theory thus provides a solution to the strong CP problem. In the following we will assume that the GUT group is spontaneously broken at energies higher than the dark confinement scales. This is mostly motivated by the fact that f a ≤ 1.2 × 10 16 GeV is required (according to an estimate based on the results of Ref. [16]) to protect the axion from too large non-perturbative gravitational effects, and GUT scales smaller than ∼ 10 16 GeV are generically at odds with proton decay.
Extending the theory is also motivated by its chiral nature: if M GUT > Λ DC and the GUT dynamics does not break SU(N DC ) × U(1) D , as we are assuming, we expect dark fermions charged under the SM to come in complete GUT representations. For the case of (non-supersymmetric) Georgi-Glashow (GG) unification [19], this can be realized by the straightforward embedding of the color triplets into fundamentals of SU(5) GUT , as shown in Table 2 terms with GUT scalars are allowed by gauge and Lorentz invariance if these fields do not carry U(1) D charges, as we will assume in the following. Since all fermions transform in complete GUT multiplets, the one-loop differential running of the gauge couplings is the same as in the SM, leaving the accuracy to which unification is realized in the minimal GG model unaltered. This is no longer true at two-loop level, since the splitting among the radiatively generated masses of the pseudo NGBs and threshold corrections can change the differential running. A non-minimal GUT sector with additional scalars or fermions can improve the unification without interfering with the present set-up.
While most of the qualitative features of the GUT theory are carried over from the simpler QCD model, the additional fields introduce some differences at low energy. We assume that the GUT group is broken at scales higher than Λ DC , hence the symmetry breaking pattern of Eq. (4.1) is enlarged to whiles has an even smaller decay width that can be induced only by a mixing withs generated by higher-order C -breaking loops. The metastability ofã ands can have an impact on cosmology, as discussed in Section 6.
Apart fromã ands, the low-energy limit of the GUT theory is also characterized by the presence of the two light NGBs a and s. The generalization of Eq. (4.3) is straightforward, 5, 1, 1) , (4.17) and the same properties for a and s are derived as for the QCD model. At energies of order Λ DC the theory can be matched to an effective Lagrangian for the four SM-singlet NGBs. This is similar to that of Eq. (4.9), with f a and N DW given again by For g GUT 0.1 and Λ DC < M GUT , the estimated value of c j is always smaller than 10 −4 .

Generic Models with n f = 2
We now generalize the classification to the case of arbitrary G SM representations irreducible. To this aim, it is convenient to introduce the following short-hand notation:  Combining the first two equations we obtain The solution p 1 = −q 1 implies p 2 = −q 2 and leads to a vector-like charge assignment, which we exclude. We are thus interested in choices of representations such that T 1 d 2 = T 2 d 1 .
One class of solutions is obtained for d 1 = d 2 , T 1 = T 2 . This is the only possibility if the choice of the r i is restricted to irreducible representations of SU(N ) with rank 2 or lower, as done in Ref. [11] to maintain perturbativity up to the Planck scale. In that case, the possible charge assignments are the same as those of Ref. [11], which lead to Type I, and for given arbitrary q 1 and p 1 , we find It is always possible to redefine the U(1) D coupling constant in such a way that one of the charges is normalized to one, e.g. q 1 = 1.
The possible choices of G SM representations r 1 , r 2 are narrowed by the confinement and perturbativity constraints of Section 2.1. As in the case of the minimal QCD and GUT models, we require that there exist a non-vanishing range of N DC in which the dark color group confines and the SM couplings remain perturbative for values of f a as low as 10 9 GeV. This criterion selects a sufficiently large and representative set of models, as we will see. Furthermore, when imposing Eq. (2.5) we assume a relation Λ DC = 100 f a between the dark confinement scale and the axion decay constant. This makes it possible to simplify the analysis and it is a good approximation in all the models considered. A thorough and more accurate analysis would require deriving for each model the exact relation between f DC and f a (by determining the embedding of U(1) PQ in the global symmetry group) and varying f a in the whole allowed range. Some models with larger ratios f DC /f a that have been excluded by our selection could be reconsidered in this way. We leave such analysis of n f = 2 models for future work and will perform it only for the n f = 3 models of Section 5. Besides passing the requirements of confinement and perturbativity, viable models must have at least two fragments with different Dynkin indices in order to have an axion, i.e. a NGB with a non-zero anomalous coupling to QCD gauge fields. The models satisfying these conditions, and the dimension of the corresponding representations, are listed in Table 3.

5
(10 ⊕ 5 1) 5 15 12 We still need to assess to what level the PQ symmetry is protected in these theories.
From Eq. (4.21) it follows that an operator of the form is always a gauge singlet, since d 1 (p 1 + q 1 ) + d 2 (p 2 + q 2 ) = 0 and d 1 , d 2 are positive integers.
If gcd(d 1 , d 2 ) > 1, where gcd(d 1 , d 2 ) denotes the greatest common divisor of d 1 and d 2 , it is possible to reduce the dimension of the operator by sending d 1 → d 1 /gcd(d 1 , d 2 ) and . This class of operators has the right form to violate the PQ global symmetries associated to the models of Table 3 and generate a potential for the axion according to the criterion of Section 3.1. The PQ-violating dimension ∆P Q , as defined in

Sec. 3.1, has an upper bound ∆P
. It follows that models with d 1 = d 2 , among which are those analyzed in the previous sections, generically have dangerous dimension-6 operators and cannot solve the axion quality problem. Models with ∆ maẍ P Q ≥ 9 have instead a chance to solve the quality problem, for suitable values of f a , depending on the actual value of ∆P Q . Determining the latter requires specifying the U(1) D charges and the embedding of U(1) PQ into the global symmetry group. We leave such analysis for future work [68]. 5 Models with n f = 3 In Section 2 we showed that models featuring a composite axion in the low-energy spectrum are those with at least three irreducible representations. Here we analyze models with exactly three irreducible representations, hence with n f = n irr = 3. We will show that in this case, for appropriate choices of the U(1) D charges, it is possible to have a sufficiently protected accidental Peccei-Quinn symmetry and thus obtain candidate solutions to the axion quality problem.
As for n f = 2, we make no a priori assumptions on the SM representations of the dark fermions. We classify models of the form detailed in Table 4 Table 5, together with the allowed range of dark colors and the minimum value of f a . Starting from any of the models of this table, except the last two, one can obtain another viable model by conjugating the representation r 3 (or equivalently r 2 ); such new theory has a different spectrum of heavy resonances but the same low-energy axion phenomenology compared to its parent theory.
Computing the minimum PQ-violating dimension ∆P Q defined in Sec. 3.1 for all the models of Tab. 5 is beyond the scope of this work. However, similarly to what we have seen for n f = 2 models, it is possible to set an upper bound on ∆P Q by exploiting the conditions of anomaly cancellation, Eq (5.1). From the first two equations, in particular, it follows that operators of the form are neutral under U(1) D provided that 17 These operators are invariant also under SU(N DC ) × G SM , manifestly violate U(1) PQ , and are of the right form to generate an axion potential through a single insertion. Their dimensionality is equal to 3 i (κ i +κ i ). Therefore, for a given model of Tab. 5, there exists an upper bound on the dimension of the PQ-violating operators given by ∆ maẍ P Q = 3κ min , where κ min = min{ i (κ i +κ i )} for positive integer κ i ,κ i satisfying Eq. (5.3). The values of ∆ maẍ P Q are reported in Tab. 5 for the selected models. All the models of Tab. 5 can address the quality problem if their value of ∆P Q equals the upper value ∆ maẍ P Q computed as above. It is in fact plausible that this equality can be achieved with a suitable choice of the U(1) D charges. In the following we shall focus on GUT models with (r 1 , r 2 , r 3 ) = (1,5, 10), where f a can be as low as 4 · 10 8 GeV. We will present their relevant properties and show that choices of U(1) D charges exist for which ∆P Q = ∆ maẍ P Q = 12. QCD models with (r 1 , r 2 , r 3 ) = (1, 3, 6) will be briefly discussed in Appendix B.

Analysis of (1,5, 10) models
Let us consider models of Table 4  Yukawa operators, in particular, are especially dangerous and one naively expects that they lead to a too large correction to the axion potential unless their coefficient is extremely small. To simplify our analysis, in the following we will consider only sets of U(1) D charges that do not allow any operator with two dark fermions. In practice, this is equivalent to require q i = −p j for any i, j, and p i = p j , q i = q j for any i = j.
From Section 2 we know that the number of NGBs that are singlets of G SM is equal to (n irr − 2). For n irr = 3, this means that the only light degree of freedom is the axion itself. The Peccei-Quinn generator can be easily shown to be Q a = diag(5, 1 5 , −1 10 ) . (5.4) The corresponding current J µ a has an anomaly in the background of SU(5) gauge fields where G µν is the SU (5) Then, the anomalous couplings to the SM gauge fields can be written as where i runs over the three simple factors of the SM gauge group SU(3) c × SU(2) EW × U(1) Y , and the anomaly matrix has the form Since the matrix A Ii has rank 2, one might be tempted to conclude that there is a particular combination of the ϕ I which does not couple to any of the anomalies, and is thus stable.
However, the mass generated radiatively by SU(5) GUT interactions for the two pseudo NGBs contained in the 24 is different from that of the pseudo NGB coming from the 75, and furthermore the mass matrix is not diagonal in the chosen basis. We explicitly verified this by computing the masses at 1 loop following Ref. [69]. Therefore, the rotation required to decouple one of the pseudo NGBs from the anomalies would induce a mixing between the ϕ I 's, 18 and we can conclude that all three states can decay through the anomalies. 18 We thank Luca Vecchi for bringing this point to our attention.
From a qualitative viewpoint, the situation is very similar to that of theã in the previous section.
Finally, U(1) D is spontaneously broken at dark confinement and the dark photon obtains a mass To derive this formula we have used the anomaly conditions (5.1).

Explicit solutions and axion quality for (1,5, 10) models
To fully specify a model, one must assign U(1) D charges that satisfy the conditions of anomaly cancellation in Eq. (5.1). The choice of charges determines the structure of higherdimensional operators and therefore the degree to which the quality problem can be solved in a given model. Finding all the integer solutions to the system (5.1) does not appear to be an easy problem and it will be tackled in a separate work [68]. However, some explicit solutions are easy to find. A particularly interesting family of solutions is given by where the parameter q is an integer number that must be different than 1 to have chiral representations. That of Eq. With the correct normalization for α D , the two subgroups can be identified with SU(5) GUT and U(1) D , the ψ i fields can be organized into a single spinorial representation of SO(10), and the same is true for the χ i . Since the 16 is anomaly free, anomaly cancellation emerges automatically and is satisfied independently by each of these two sets of fermions.
The model can thus be viewed as the low-energy limit of a more minimal theory with an extended unification group SO(10) L × SO(10) R , and a structure that is highly reminiscent of the theories discussed in Ref. [21]. The precise dark fermion content is shown in Tab. 6, from which it is apparent that the doubling of the SO(10)'s is necessary to ensure a chiral structure. 20 As in [21], all SM generations can be taken to transform as Therefore, ∆P Q ≤ 12 in this case, as shown in Tab. 5. We have thus performed a search for models with ∆P Q = 12 (dubbed high-quality models in the following), which can resolve the quality problem for values of the axion decay constant as high as 10 12 GeV and explain the whole DM abundance in terms of axions. The result of this search are reported below.
A complete classification depends on the precise representations of the GUT scalars required by the unification dynamics. This is because these fields can help to render some PQ-violating operators gauge invariant, at the price of a slight increase in their dimension.
In general, only the following multiple insertions have d eff < 12: • Two 4-fermion operators, with at most three GUT scalars in their product.
• Three 4-fermion operators, with at most one GUT scalar in their product.
In principle these operators might contain SM fermions in flavor-singlet combinations. In practice, SU(N DC ) invariance and the absence of 2-fermion operators exclude this possibility for 4-fermion operators. One can instead construct gauge invariant 6-fermion operators with four dark fermions and two SM fermions. 21 However, for representations r 1 = 1, r 2 =5, r 3 = 10, it is possible to show that the product of four dark fermions can be always turned into a singlet of SU(5) by multiplying it with one or two GUT scalars. 22 In this case, the 6-fermion operator can be replaced by an operator with four dark fermions and up to two scalars that has lower dimension and same flavor structure. One can thus focus on operators made of dark fermions only.
Although the classification of dangerous operators depends on the scalar content of the model, it is possible to impose a stronger condition and perform the same analysis without imposing SU(5) gauge invariance on all operators that feature GUT scalars. This selects a smaller class of models, dubbed robust in the following, that have high quality independently of the choice of GUT scalars, provided these do not carry U(1) D charge. 23 Besides selecting robust models, we also performed a search for high-quality models with a minimal scalar sector, i.e. GUT scalars transforming as 5 + 24 of SU(5). Table 7 illustrates the result of our numerical study of the parameter space. However, it would need to be paired with another 6-fermion operator to be of the correct form, resulting in an effective dimension d eff = 14. 22 Given any such operator, one can always contract as many of its upper and lower indices as possible with delta tensors and remain with only one type of indices. Their corresponding representation is a product of antisymmetric tensors (each of rank 1 or 2) which always contains a fully antisymmetric tensor. In the case of SU(5) any antisymmetric tensor can be contracted with at most two fundamentals to give a singlet.
The presence of a GUT scalar transforming as a fundamental of SU(5) is instrumental to embed the Higgs boson. 23 In practice, robust models are selected by imposing SU(5) gauge invariance only on the 6-fermion and 4-fermion operators involved in a double insertion, and on two of the 4-fermion operators involved in a triple insertion. 24 To perform the classification, we have made use of the Mathematica package LieART [70,71]. with the SM at low energy. On the other hand, their presence can have an impact on the cosmological evolution and can thus be tested by celestial observations. The cosmology of axion models is complex, for a review see for instance Refs. [72,73]. A scenario where the scale of inflation H I is larger than the dark confinement (and PQ) scale, or where the PQ symmetry is restored after reheating, appears to be strongly disfavoured in our models. Indeed, in this case the degrees of freedom in the PQ sector would thermalize with the SM bath after reheating, thanks to the QCD interactions of the dark fermions.
After confinement, the heavy and accidentally-stable states of the PQ sector (such as dark baryons) would acquire a relic abundance much larger than the observed DM density, in conflict with observations (for a detailed analysis see [11]). Moreover, as we show in Appendix A, the domain wall number is N DW > 1 in all of the examples, allowing the formation of stable domain wall and string networks. If the PQ phase transition occurs after inflation, they will quickly dominate the energy density of the universe, in conflict with standard cosmology [74,75].
We focus, therefore, on the scenario where the PQ symmetry is broken during inflation and the strongly-interacting sector is in its confined phase. In this case, all the accidentallystable states with mass ∼ Λ DC are heavier than the scales of inflation and reheating, and are not populated at those temperatures. We thus consider only the light dark states, i.e.
the QCD axion a and the other SM-neutral NGBs. The latter include, in general, NGBs that are exact up to non-perturbative gravitational effects and thus almost massless, and heavier pseudo NGBs whose associated symmetry is explicitly broken by GUT corrections.
In the following we will analyze models with n f = 2, where these states are respectively the singlet s and, in GUT models, theã ands; GUT models with n f = 3 and irreducible representations have instead only pseudo NGBs, whose cosmology and phenomenology are similar to those ofã.
Let us discuss first the massive pseudo NGBs. The decay rate of Eq. in the range mã ,s < T RH < f a . In the first case, the statesã ands are not popu-lated at reheating, the cosmological history is the standard one and the DM abundance can be reproduced by the axion with the usual misalignment mechanism. In the second case,ã ands can be populated at reheating. If T RH is larger than the temperature GUT ) > mã at which processes gg ↔ gã go out of equilibrium,ã can reach thermal equilibrium with the SM bath. After an early freeze-out, its energy density will be large and can lead to an epoch of early matter domination, before decaying to photons prior to BBN. In this scenario the abundance of axions given by the misalignment mechanism can be depleted, increasing the value of f a required to match the observed DM abundance for fixed initial misalignment. The interactions ofs with the SM bath are more suppressed than those ofã, and the region of parameter space where it reaches thermal equilibrium, if any, will be small. When the reheating temperature is higher than their mass but lower than their freeze out temperature,ã ands have an abundance that depends on the reheating dynamics and cannot be robustly estimated. Considering that its decay can happen after BBN, a too large density ofs can be problematic and in conflict with observations. It would be particularly interesting to explore the possible role ofã ands in cosmic inflation and their imprints as primordial non-gaussianities in cosmological correlators, especially because their mass is parametrically separated from Λ DC and can be close to the scale of inflation. We leave a more detailed analysis of the role ofã ands in cosmology to a future work. In the following, we focus on the first scenario where T RH is smaller than the mass ofã ands, and give a brief overview of the low-energy phenomenology.
Let us start considering the QCD axion a. Its phenomenological properties have been thoroughly explored and are well predicted in terms of the UV parameters by combining chiral perturbation theory and Lattice QCD results [57,76]. Using the notation and the results of Ref. [76], we quote here the predictions in our models for the low-energy parameters relevant for experimental searches. At energies below the QCD confinement scale, and including low-energy QCD contributions, the axion mass, the axion-photon coupling and the axion-nucleon couplings are: m a = 5.70 (7) 10 12 GeV f a µeV, while the coupling with electrons is generated only through radiative electromagnetic corrections and is suppressed [77,78], c e 3 · 10 −4 . The vanishing of the leading-order UV contributions to the axion-fermion couplings in our model gives a sharp prediction for these quantities, common to all the so-called "hadronic" axion models. Corrections to these predictions from GUT thresholds can be generated, as mentioned in Section 4. values of the mass in the range 10 −20 eV < m s < 10 −10 eV, the existence of this particle may be tested through black hole observations [82,83]. An estimate along the lines of Ref. [16] seems to suggest f a 5 · 10 16 GeV for this to be the case, with m s exponentially smaller for lower f a . Since f a < 1.2 · 10 16 GeV is required in our models to have a highquality axion, the mass and the abundance of s should be extremely small, according to these estimates. In light of the exponential sensitivity to order one numerical factors, and keeping an open mind on the nature of the UV effects, it could be possible for s to have a mass in the right range to play the role of a fuzzy dark matter component with abundance generated through the misalignment mechanism [84].
To conclude, both the QCD axion a and the singlet s can have a thermal component, which behaves as dark radiation and can be parametrised by a correction in the effective number of neutrinos ∆N eff ν . An irreducible component for the QCD axion comes from the thermal production of axions from the SM bath [85][86][87]. Its value depends on the reheating temperature and on f a , but it is always bounded by ∆N eff ν ≤ 0.026. Additional contributions can be present, both for a and s, depending on the dynamics of reheating or possibly from the decay of a population of primordial black holes [88]. In our treatment these can be seen as initial conditions and do not give rise to sharp predictions, but can be a helpful probe of our models.

Conclusions
The QCD axion solution to the strong CP problem hinges on the existence of a spontaneously broken Peccei-Quinn symmetry and is realized in terms of a very simple EFT at low energies. Despite this apparent simplicity, the PQ symmetry is usually imposed by hand, and explicit models are extremely sensitive to PQ-violating effects in the UV. Global symmetries are ultimately broken by quantum gravity and can be at best approximate. This is especially true for the PQ symmetry, which must be also broken by QCD instantons. Approximate global symmetries emerge at low energy in the effective field theory if they are broken only by irrelevant operators, and in that case they are called accidental.
Proton stability and the smallness of neutrino masses can be elegantly explained as consequences of accidental symmetries of the SM. The stability of dark matter could also be the consequence of an accidental symmetry. In all these cases, it is enough to protect the relevant global symmetry up to operators of dimension six to be compatible with current experimental data. Interestingly, breaking effects at the dimension-6 level imply proton decay rates (for a cutoff of order of the GUT scale) and DM decay rates (for a cutoff of order of the Planck scale and DM masses in the O(100 TeV) range) that are not far from current experimental sensitivities. This means that the paradigm of accidental symmetries could be tested experimentally in the near future. It is hence natural to consider whether the Peccei-Quinn symmetry may also be accidental, that is, emerge at low energy in the effective field theory. While this is an attractive solution, the fact that the axion potential is generated by QCD dynamics at an energy vastly smaller than the cutoff scale, and the existence of stringent bounds from neutron EDM experiments put the PQ symmetry in a different class compared to the accidental symmetries mentioned above. In order for the QCD axion solution to work, PQ-breaking operators must have at least dimension 9, while the most attractive scenario where the axion can account for the whole DM abundance requires dimension 12 or larger. This makes the model building of accidental axions much harder than for other SM extensions.
It is clearly possible that the PQ symmetry is not accidental, at least in the sense that it does not emerge accidentally due to a separation of scales in the effective field theory.
One can imagine a situation where PQ-breaking effects arise only at the non-perturbative level in the full UV theory, such that any other contribution to the axion potential is much smaller than that from QCD instantons. Models explicitly realizing this scenario have been proposed in the context of string theory. In these constructions the 4D effective action has a non-linearly realised PQ symmetry which acts as a shift symmetry on a fundamental scalar field. The strong CP problem is solved as long as the UV sources of PQ breaking (such as string theory instantons) are small enough, which can be achieved in some of the proposed models. In theories of this kind the axion quality entirely relies on the properties of the UV dynamics and cannot be explained within the effective field theory.
In this work, we have introduced a class of composite axion models where the PQ symmetry is an accidental symmetry of the effective field theory.
Our constructions are compatible with an SU (5)  It is useful to compare the level of PQ symmetry protection achieved in the class of models studied in this work with that of previous works on composite axions. In the QCD models of Ref. [17], the PQ-violating dimension depends on the number m of colors of the weakly-gauged non-abelian group: ∆P Q = 3m. For this reason, a careful analysis of perturbativity of the QCD coupling is crucial to assess the viability of these models as solutions of the quality problem [18]. In the case of simple GUT extensions of the models of Ref. [17], we find that multiple insertions with effective dimension 10 can generate an axion potential, so that ∆P Q = min{3m, 10}, see Appendix D. The models of Ref. [24] can have an accidental PQ symmetry up to dimension 9, while ∆P Q = 12 in the models of Refs. [20,21]. The constructions of Ref. [25] can lead to an even higher level of protection, but do not seem complete, in that gauge anomaly cancellation is postulated through the inclusion of an unspecified set of massless fermions whose charges and low-energy dynamics have not been investigated yet.
The PQ-violating dimension in our models depends on the choice of U(1) D charges.
By performing a detailed scan of n f = 3 models with GUT representations (r 1 , r 2 , r 3 ) = (1,5, 10) or QCD representations (r 1 , r 2 , r 3 ) = (1, 3, 6), we have found several charge assignments that lead to ∆P Q = 12. These models are perturbative for values of the axion decay constant as low as 4 · 10 8 GeV and thus give a robust solution to the axion quality problem. For f a ∼ 10 12 GeV they can also explain the DM abundance entirely in terms of axions through the misalignment mechanism. Our classification revealed several other models with n f = 2 and n f = 3 that can also address the quality problem and where the PQ-violating dimension can be higher than 12. For example, we identified GUT models with n f = 3 and QCD models with n f = 2 or 3 where ∆P Q can be as large as 18 and 15 respectively, see Tabs. 5 and 3. Determining the actual level of protection of these constructions requires a more elaborate analysis that we leave for future work [68].
Some of the models that we have discussed are close to the edge of the conformal window (for instance those of Section 4.2 with the choice N DC = 3), and could feature a walking dynamics. In this case the PQ-violating could have large anomalous dimensions: it would be interesting to understand if such a scenario can improve the level of protection of the PQ symmetry in this class of models or not.
The low-energy dynamics of our models shares all the universal properties of hadronic axion theories. In particular, the GUT models imply an E/N ratio equal to 8/3, and will be thus tested by on-going and planned experiments aiming at observing QCD axions.
In order to probe the Peccei-Quinn dynamics and distinguish between different proposals, however, further experimental observables are needed. If f a ∼ 10 12 GeV and the breaking of PQ symmetry is induced by d = 12 local operators generated at M Pl , then the models of Section 5 can be tested through the next generation neutron EDM experiments. An experimental indication that the PQ symmetry is accidental could therefore come in the near future, similarly to the case of baryon number and DM symmetry. Another way to probe the PQ dynamics and distinguish our GUT models from other ones could come from the parametrically lighter states that acquire mass through non-renormalizable interactions.
These states may have an impact on cosmological observables and give the tantalizing opportunity of testing the GUT and PQ dynamics through celestial observations. Besides solving the strong CP problem, the QCD axion could also be the main component of DM, either through the ordinary misalignment mechanism or through a possible realization of the recently proposed kinetic misalignment mechanism [89,90] in the composite scenario.
Furthermore, models with n f = 2 predict the existence of an additional massless scalar s.
This could be a component of DM, if UV effects generate a mass in the appropriate range and it is populated through a primordial mechanism, and can also give a contribution to dark radiation that, once added to the axion contribution, could be tested by future CMB and large-scale structure observatories.
The global PQ symmetry is broken both explicitly (by the QCD anomaly) and spontaneously (by the SU(N DC ) condensate). As first pointed out by Sikivie in Ref. [75], these symmetry breaking effects often leave unbroken some discrete symmetry groups, which can have important conceptual and phenomenological consequences. In this Appendix we will provide a general analysis of the discrete symmetries related to U(1) PQ in composite axion models based on a QCD-like confining gauge group.
Let us denote by G V the exact global flavour symmetry group preserved by the vacuum.
In our models this is a vectorial group of the form G V = U(1) V,1 × · · · × U(1) V,n irr , where each factor U(1) V,i acts with charge Q V,i (ψ j ) = δ ij and Q V,i (χ j ) = −δ ij . As in the main text, we fix the normalization of the U(1) PQ charges of the dark fermions ψ i , χ i so that they are coprime integers; with this choice the continuous parameter associated to the symmetry transformation is defined in the range [0, 2π). In general, the U(1) PQ PQ group can have a non-vanishing intersection with the unbroken flavor group, In other words, the SU(N DC ) condensates induce the spontaneous breaking The local dark color gauge group SU(N DC ), which leaves invariant physical states of the Hilbert space, is associated to a global group that instead can act non-trivially on the Hilbert space (see for example [91]). Such global group is linearly realized in the confining phase of vector-like theories, and physical states are classified according to its irreducible representations. One may therefore ask whether this global symmetry plays any role in our analysis of the discrete symmetries related to U(1) PQ . Since U(1) PQ commutes with SU(N DC ), the transformations of such a global SU(N DC ) that are contained in U(1) PQ have to be part of the (global) center group Z N DC . The latter acts as a flavour symmetry, and since it is preserved by the vacuum, it must be a subgroup of the vectorial group: Therefore, there is no need to consider separately transformations of the center Z N DC .
Invariance under the symmetry Z V restricts the physical domain of the axion. Denoting and the domain of PQ shift transformations on the axion is correspondingly restricted.
On the other hand, the anomaly breaks explicitly the U(1) PQ symmetry down to a discrete group Z A . Since the vectorial group is anomaly free, it follows that Z A contains Z V as a subgroup, and in particular its order N A is an integer multiple of N V . The quotient group has order N PQ = N A /N V and is an exact global symmetry group. It is non-linearly realised on the vacuum and acts on the axion as a discrete shift symmetry: It is straightforward to show, using the effective chiral lagrangian approach reviewed in [76], that (A.4) is a symmetry of the non-perturbative axion potential. It follows that the axion potential has N PQ distinct minima and the so-called domain wall number N DW equals N PQ .
Let us now determine N A and N V in our models. As in the main text, we use a normalization where Tr(t a t b ) = δ ab /2 for fields in the fundamental representation and use a Dirac notation. Under an anomalous PQ transformation with parameter α, the QCD theta angle shifts by 25  with κ i , κ i arbitrary integers. Summing the equations it follows that α Q i PQ = π(κ i + κ i ), and since the charges Q i PQ are coprime integers, the only non trivial solution in the range α ∈ (0, 2π) is α = π. We find, therefore, that Z V = Z 2 and N V = 2.
We arrive at our final result which is always an integer multiple of N DC and, in particular, greater than 1.

B (1, 3, 6) QCD models
In this Appendix we give a brief overview of n f = 3 QCD models with (r 1 , r 2 , r 3 ) = (1, 3,6) or (1, 3, 6), and classify their charge assignments that can solve the quality problem. 25 The SU(5) GUT case is obtained by replacing the QCD generators with the corresponding GUT generators, obtaining the corresponding A5.
The constraints of Sec. 2.1 imply 3 ≤ N DC ≤ 9 for these models. The PQ generator takes the form High quality models with ∆P Q = 12 can be characterized as done in Section 5.2. The only difference with respect to the GUT case is the absence of SU(3) c charged scalars, since we only include the SM Higgs field. Therefore, there is no longer a distinction between model dependent and robust solutions to the quality problem. Following the same procedure as in Section 5.2, we find the results reported in Table 8. The three charge assignments with n max = 10 that solve the quality problem with (r 1 , r 2 , r 3 ) = (1,3, 6) and Here and in the following, the U(1) D charges are reported in the format (p 1 , p 2 , p 3 ; q 1 , q 2 , q 3 ) to keep the notation more compact. With (r 1 , r 2 , r 3  The first two are also a solution for (r 1 , r 2 , r 3 ) = (1, 3, 6) and N DC = 4. C High quality (1,5, 10) models Following the analysis of Section 5.2, we present the charge assignments for n f = 3 GUT models with (r 1 , r 2 , r 3 ) = (1,5, 10) or (1, 5, 10) and ∆P Q = 12. We report all the highquality solutions with n max = 10 classified in Table 7.

D Randall's model
The model by Randall [17] is based on a product gauge group G = SU(N DC ) × G W × G SM , with a non-abelian weak factor G W = SU(m), and is defined by the chiral field content of Table 9. The original work of Ref. [17] Table 9: Composite axion model of Ref. [17]. In the original work only QCD interactions were considered, with G SM = SU(3) c and r = 3; the minimal GUT extension is obtained by taking G SM = SU(5) and r = 5. The fermions ψ 2 , χ 1 , χ 2 appear in multiple copies, labelled by the indices i = 1, . . . , dim(r), j = 1, . . . , m, k = 1, . . . , m · dim(r).
As already noticed in Ref. [18], the leading PQ-violating operators that generate a po- GeV and for f a = 10 12 GeV. The more sophisticated analysis of perturbativity performed in Ref. [18] suggests that the axion solution is actually spoiled for values of the Peccei-Quinn scale f DC 10 11 GeV.
One can wonder whether multiple insertions can induce PQ-violating effects that are more dangerous than those from single insertions. We find that multiple insertions are less dangerous than single ones in the case of QCD models. On the other hand, in the GUT case the leading PQ-violating operators have dimension 7 and include a GUT scalar Φ in the fundamental representations of SU(5); they are of the form O 1 = ψ 1 ψ 2 (χ 1 ) 2 Φ and O 2 = ψ 1 ψ 2 (χ 2 ) 2 Φ † . According to the criteria discussed in Section 3.1, these operators are harmless at the level of a single insertion but give rise to a potential for the axion when combined in a double insertion. On dimensional ground this effect is equivalent to the insertion of a single operator of effective dimension 10, and thus provides the leading PQ-breaking effect when m ≥ 4. This suggests that GUT models can solve the quality problem only for low values of the axion decay constant, f a ∼ 10 9 GeV.