aainstitutetext: Università di Roma La Sapienza, Piazzale Aldo Moro 5, 00185
Roma, Italybbinstitutetext: Istituto Nazionale di Fisica Nucleare (INFN) -
Sezione di Pisa
Polo Fibonacci Largo B. Pontecorvo, 3, 56127 Pisa, Italyccinstitutetext:
Department of Physics, Center for Theoretical Physics, Columbia University,
538 West 120th Street, New York, NY 10027, U.S.A.ddinstitutetext: Rudolf
Peierls Centre for Theoretical Physics, Beecroft Building, Clarendon
Laboratory,
Parks Road, University of Oxford, OX1 3PU, UK

# Chiral models of composite axions and accidental Peccei-Quinn symmetry

Roberto Contino c Alessandro Podo d Filippo Revello

###### Abstract

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 $\rm 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 $\rm 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.

##  1 Introduction

The strong CP problem can be elegantly solved by the QCD axion Peccei:1977hh ;
Peccei:1977ur ; Wilczek:1977pj ; Weinberg:1977ma , which assumes the existence
of a spontaneously broken global symmetry - known as the Peccei-Quinn (PQ)
symmetry - anomalous under $\,{\rm 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 $\theta$ 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.
Georgi:1981pu , while that of a composite axion was put forward by Kim
Kim:1984pt  and then further studied in Refs. Kaplan:1985dv ; Choi:1985cb .
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 Harigaya:2016rwr ; Co:2016akw ;
Contino:2020god ; 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 $\theta$. 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 Kamionkowski:1992mf ; Barr:1992qq ; Holman:1992us
; Ghigna:1992iv . However, it is still unclear whether non-perturbative
quantum gravity effects will be exponentially suppressed by a factor
$e^{-cM_{Pl}/f_{a}}$, 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}\leq 1.2\cdot 10^{16}\,\rm GeV$ –
see Ref. Hebecker:2018ofv  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 Randall:1992ut , who put forward a chiral model
based on a new product gauge group $SU(n)\times 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 Dobrescu:1996jp .
Additional models based on non-abelian product gauge groups, which generalize
Randall’s construction and improve its quality, have been put forward in Refs.
Redi:2016esr ; Vecchi:2021shj . Supersymmetric composite axion models, also
based on product gauge groups, have been instead constructed in Refs.
Lillard:2017cwx ; Lillard:2018fdt . Other attempts to build composite axions
and solve the quality problem by means of chiral gauge theories include the
models of Ref. Gavela:2018paw  and Ref. Fukuda:2017ylt . The first makes use
of a chiral and confining simple gauge group, while the second exploits an
axial ${\rm U}(1)$ gauge group and is more similar to our approach.
Holographic constructions dual to composite axion models have been analysed in
Refs. Choi:2003wr ; Flacke:2006ad ; Cox:2019rro ; Bonnefoy:2020llz ;
Yamada:2021uze ; Lee:2021slp . 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.
Chun:1992bn ; Cheng:2001ys ; DiLuzio:2017tjx ; Duerr:2017amf ;
Bonnefoy:2018ibr ; Lee:2018yak ; Hook:2019qoh ; Ardu:2020qmo ; Nakai:2021nyf ;
Darme:2021cxx ; Bhattiprolu:2021rrj .

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 Svrcek:2006yi ; Conlon:2006tq  (see
also the recent work Demirtas:2021gsq ). 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 $\,{\rm
SU}(5)$ Grand Unification dynamics has been first proposed in Dine:1981rt ;
Wise:1981ry  and recently re-examined in Refs. Boucenna:2017fna ;
DiLuzio:2018gqe ; FileviezPerez:2019fku ; FileviezPerez:2019ssf ;
Quevillon:2020aij . All these models differ from the $\,{\rm 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.
Redi:2016esr ; Vecchi:2021shj  considered $\,{\rm 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.
Georgi:1981pu ; DiLuzio:2020qio . For an extensive review of QCD axion models
we refer the reader to Ref. DiLuzio:2020wdo .

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.

##  2 Model building

We consider theories where the only UV degrees of freedom in the Peccei-Quinn
sector are fermions and gauge fields111The SM sector can include scalars, such
as the Higgs field or the ordinary heavy scalars of a grand unified UV
completion.. The gauge group is assumed to be of the form $\mathcal{G}={\rm
G}_{\rm D}\times{\rm G}_{\rm SM}$, where ${\rm G}_{\rm SM}$ is either the SM
gauge group or one of its grand unified counterparts, and ${\rm G}_{\rm 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
$\mathcal{G}$ are chiral. In order to have control over the IR dynamics in the
PQ sector, we assume that ${\rm G}_{\rm D}$ includes a confining dark color
${\rm SU}(N_{\rm DC})$ subgroup under which the new fermions transform as
fundamental plus anti-fundamental vectorlike representations. We take $N_{\rm
DC}\geq 3$, so that the fundamental representation of ${\rm SU}(N_{\rm DC})$
is complex. To ensure that the SM gauge symmetries are not broken
spontaneously by ${\rm SU}(N_{\rm DC})$ confinement we require that the new
fermions transform also as vectorlike representations of ${\rm G}_{\rm SM}$.
The chirality of the model then rests on the existence of an additional
subgroup ${\rm G}_{\rm W}\subset{\rm G}_{\rm D}$ whose dynamics is
perturbative at the ${\rm SU}(N_{\rm DC})$ confinement scale. The case with
${\rm G}_{\rm W}={\rm SU}(m)$ was first discussed in Ref. Randall:1992ut  and
then further analysed and generalised in Refs. Dobrescu:1996jp ; Redi:2016esr
. In this work we focus on the case in which the weak group is abelian, ${\rm
G}_{\rm W}={\rm U}(1)_{\rm 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

$\psi_{i}\sim(\square,p_{i},r_{i}),\qquad\chi_{i}\sim(\bar{\square},q_{i},\bar{r}_{i}),\quad\quad
i=1,...,n_{f}$ (1)

under ${\rm SU}(N_{\rm DC})\times{\rm U}(1)_{\rm D}\times{\rm G}_{\rm SM}$.
The ${\rm U}(1)_{\rm D}$ charges $p_{i}$ and $q_{i}$ are defined such that
$p_{i}\not=p_{j}$, $q_{i}\not=q_{j}$ for $i\not=j$ and $p_{i}\not=-q_{i}$.
Therefore, $n_{f}$ is the number of different ${\rm U}(1)_{\rm D}$ charges
assigned to the $\psi$ (or equivalently to the $\chi$) fields, and the SM
representations $r_{i}$ are in general reducible. The ${\rm U}(1)_{\rm D}$
charges are chosen so that the representations are overall chiral. This
requirement and the cancellation of ${\rm U}(1)_{\rm D}$ gauge anomalies imply
that the multiplicity of fermions $n_{f}$ has to be at least 2: $n_{f}\geq 2$.

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 ${\rm U}(1)_{\rm D}$ gauge field:

$\mathcal{L}_{\rm BSM}=-\dfrac{1}{4}{\cal G}^{a}_{\mu\nu}{\cal
G}^{a,\mu\nu}-\dfrac{1}{4}F^{D}_{\mu\nu}F^{D\,\mu\nu}+\dfrac{\varepsilon}{2}F^{D}_{\mu\nu}B^{\mu\nu}+\sum_{\kappa}\psi^{\dagger}_{\kappa}iD_{\mu}\bar{\sigma}^{\mu}\psi_{\kappa}+\sum_{\kappa}\chi^{\dagger}_{\kappa}iD_{\mu}\bar{\sigma}^{\mu}\chi_{\kappa}\,.$
(2)

Here ${\cal G}^{a}_{\mu\nu}$, $F^{D}_{\mu\nu}$ and $B_{\mu\nu}$ are the ${\rm
SU}(N_{\rm DC})$, ${\rm U}(1)_{\rm 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 ${\rm
SU}(N_{\rm DC})$ theta term to zero through an axial field redefinition. The
label $\kappa$ runs over all irreducible representations of ${\rm SU}(N_{\rm
DC})\times{\rm G}_{\rm SM}\times{\rm U}(1)_{\rm D}$, whose number we will
denote by $n_{irr}$.222One has $n_{irr}=n_{f}$ if all the $r_{i}$ are
irreducible representations. The Lagrangian of Eq. (2) has a $[{\rm
U}(1)_{V}]^{n_{irr}}\times[{\rm U}(1)_{A}]^{n_{irr}}$ accidental global
symmetry at the classical level.333The accidental symmetry is bigger than that
if any of the irreducible representations inside a given $r_{i}$ has a non-
trivial multiplicity. The $n_{irr}$ vectorial ${\rm U}(1)_{V}$ factors are
linearly realized and each of them implies an accidentally conserved charge.
One of the $n_{irr}$ axial ${\rm U}(1)_{A}$ factors is anomalous under ${\rm
SU}(N_{\rm DC})$, and one is gauged by ${\rm U}(1)_{\rm D}$, whose vector
boson acquires a mass through the Higgs phenomenon; the others are
spontaneously broken by the ${\rm SU}(N_{\rm DC})$ condensates
$\langle\psi_{\kappa}\chi_{\kappa}\rangle$. 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}\geq 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 ${\rm U}(1)_{\rm D}\times{\rm G}_{\rm SM}$ weak gauging is switched
off, the global invariance of the Lagrangian of Eq. (2) gets enhanced to
$\,{\rm SU}(N_{f})_{L}\times\,{\rm SU}(N_{f})_{R}\times U(1)_{V}$, where
$N_{f}=\sum_{i=1}^{n_{f}}{\rm dim}(r_{i})$ is the multiplicity of Dirac
fermions transforming in the fundamental representation of ${\rm SU}(N_{\rm
DC})$. Below the ${\rm SU}(N_{\rm 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_{f}^{2}-1)$ NGBs. By turning on ${\rm
U}(1)_{\rm D}\times{\rm G}_{\rm SM}$, one of these NGBs becomes the
longitudinal polarization of the massive ${\rm U}(1)_{\rm 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.
Contino:2020god , where the construction of Refs. Harigaya:2016rwr ;
Co:2016akw  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 the precise relation between this work and
Ref. Contino:2020god  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 $\Lambda_{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. (52)
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}\gtrsim
4\cdot 10^{8}\;\rm GeV$ Raffelt:2006cw  (see however Ref. Bar:2019ifz ). 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}\sim 10^{12}\,$GeV (for a recent numerical analysis see Refs.
Bonati:2015vqz ; Borsanyi:2016ksw ). Motivated by these considerations, we
will consider $f_{a}\sim 10^{9}\,$GeV and $f_{a}\sim 10^{12}\,$GeV as
benchmark values when analyzing our models.

###  2.1 IR confinement and UV perturbativity

Our construction relies on the assumption that the ${\rm SU}(N_{\rm DC})$
dynamics confines in the infrared. Since the new fermions have no mass term,
${\rm SU}(N_{\rm DC})$ must be asymptotically free; at one loop this gives the
bound $N_{\rm DC}\geq 2N_{f}/11$. A stronger constraint $N_{f}<N_{\rm conf}$
is obtained by considering the lower boundary of the conformal window. The
value of $N_{\rm conf}$ is a non-perturbative property that cannot be derived
with perturbative computations. Lattice simulations suggest that $N_{\rm
conf}\approx 12$ for an $\,{\rm SU}(3)$ gauge theory with fermions in the
fundamental representation, see Ref. DeGrand:2015zxa . In selecting our models
we will use a naive linear extrapolation of this result to generic numbers of
dark colors (see for instance Refs. Appelquist:1996dq ; Appelquist:1999hr ;
Poppitz:2009uq ; Poppitz:2009tw  for some evidence in support of this
assumption) and require

$N_{f}\leq 4N_{\rm DC}\,.$ (3)

A second condition we impose on our theories is that the Standard Model gauge
couplings should remain perturbative until the Planck scale. In practice, we
require that there are no Landau poles below the Planck mass. Consider for
example the QCD coupling; it is useful to define an effective number of new
Dirac fermions in the fundamental of $\,{\rm SU}(3)_{\rm c}$ as

$\Delta N_{f}^{(\rm QCD)}=2N_{\rm
DC}\sum_{i=1}^{n_{f}}\sum_{\alpha}T(r_{i}^{(\alpha)})\,,$ (4)

where $r_{i}^{(\alpha)}$ are the irreducible representations of $\,{\rm
SU}(3)_{\rm c}$ contained in $r_{i}$, and $T(r_{i}^{(\alpha)})$ are their
corresponding Dynkin indices.444We normalize the Dynkin index so that it is
equal to $1/2$ for the fundamental representation. Running the QCD coupling up
to the Planck scale, the perturbativity requirement implies the following
condition:

$\Delta N_{f}^{(\rm QCD)}\leq\dfrac{3}{2}\left(\log\dfrac{M_{\rm
Pl}}{\Lambda_{\rm
DC}}\right)^{-1}\left(\dfrac{2\pi}{\alpha_{s}(m_{Z})}+b^{(3)}_{\rm
SM}\log\dfrac{M_{\rm Pl}}{m_{Z}}\right),$ (5)

where $\alpha_{s}(m_{Z})=0.118$ is the QCD coupling constant at the Z pole
Zyla:2020zbs , $m_{Z}=91\,\rm GeV$ is the $Z$ mass and $b^{(3)}_{\rm SM}=7$ is
the SM one-loop beta function coefficient.555More precisely, one should
integrate out the top quark at its mass scale $m_{t}$ and use
$b^{(3)\,\prime}_{\rm SM}=7+2/3$ to run from $m_{Z}$ to $m_{t}$. This leads to
a small correction to Eq. (5) that we neglect for simplicity. Numerically,
this gives $\Delta N_{f}^{(\rm QCD)}\lesssim 26.7$ for $\Lambda_{\rm
DC}=10^{11}\,\rm GeV$ and $\Delta N_{f}^{(\rm QCD)}\lesssim 35.6$ for
$\Lambda_{\rm DC}=10^{13}\,\rm GeV$. More in general, since $\Delta
N_{f}^{(\rm QCD)}\propto N_{\rm DC}$ and $\Lambda_{\rm DC}\propto\sqrt{N_{\rm
DC}}f_{a}$ (as it will be shown in the following), the perturbativity
condition of Eq. (5) sets a lower bound on $f_{a}$ for fixed $N_{\rm DC}$, or
equivalently an upper bound on $N_{\rm DC}$ for a given $f_{a}$.
Perturbativity of the hypercharge coupling up to the Planck scale similarly
requires a bound

$2N_{\rm
DC}\sum_{i}\sum_{\alpha}y_{i,\alpha}^{2}\text{dim}(r_{i}^{(\alpha)})\leq\frac{3}{4}\left(\log\dfrac{M_{\rm
Pl}}{\Lambda_{\rm
DC}}\right)^{-1}\left(\dfrac{2\pi}{\alpha_{s}(m_{Z})}+b^{(1)}_{\rm
SM}\log\dfrac{M_{\rm Pl}}{m_{Z}}\right),$ (6)

where $y_{i,\alpha}$ is the hypercharge of the representation
$r_{i}^{(\alpha)}$ and $b^{(1)}_{\rm SM}=-41/6$.

In the following, we will focus on models satisfying the conditions of Eqs.
(3), (5) and (6).

##  3 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 $\theta$ at
the minimum of the potential of order

$\Delta\theta\approx|c|\,\phi_{CP}\left(\frac{M_{\rm Pl}}{\Lambda_{\rm
QCD}}\right)^{4}\left(\frac{f_{a}}{M_{\rm Pl}}\right)^{\Delta_{\cancel{PQ}}},$
(7)

where $\Delta_{\cancel{PQ}}$ is the dimension of the operator, $c$ is a
dimensionless effective coefficient and $\phi_{CP}$ is a $CP$-violating phase.
The current experimental bound on the electric dipole moment (EDM) of the
neutron nEDM:2020crw  requires $\Delta\theta\lesssim 10^{-10}$, which
translates into a lower bound on the dimension of the PQ breaking operator

$\Delta_{\cancel{PQ}}>\dfrac{10+\log_{10}\\!\left(|c|\,\phi_{CP}\dfrac{M_{\rm
Pl}^{4}}{\Lambda_{\rm QCD}^{4}}\right)}{\log_{10}\\!\left(\dfrac{M_{\rm
Pl}}{f_{a}}\right)}\,.$ (8)

Assuming $|c|,\phi_{CP}\sim 1$ and using $M_{\rm Pl}\sim 10^{19}\,$GeV,
$\Lambda_{\rm QCD}\sim 100\,$MeV, this gives $\Delta_{\cancel{PQ}}\gtrsim 9$
for $f_{a}\sim 10^{9}\,$GeV. Values of the PQ scale of order $f_{a}\sim
10^{12}\,{\rm GeV}$, which allow the axion to reproduce the observed dark
matter abundance naturally through the misalignment mechanism, lead to a
stronger bound $\Delta_{\cancel{PQ}}\gtrsim 13$. In theories of composite
axion, for a purely fermionic PQ-violating operator generated through diagrams
with $\ell$ loops, one expects the effective coefficient to scale as $c\sim
g_{UV}^{2}(4\pi/g_{UV})^{\Delta_{\cancel{PQ}}/3-2\ell}$, where $g_{UV}$ is a
UV coupling. This corresponds to values of $c$ larger than 1, hence to a
stronger constraint on $\Delta_{\cancel{PQ}}$, unless
$\ell>\Delta_{\cancel{PQ}}/6-1$ and $g_{UV}$ is sufficiently small.666The more
accurate estimate of $\Delta\theta$ in composite axion theories performed by
Ref. Vecchi:2021shj  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. Vecchi:2021shj  though still
quite large. For $\Delta_{\cancel{PQ}}=12$, including such factors leads to an
upper bound on $f_{a}$ stronger by roughly one order of magnitude. Such strong
bounds on the dimension $\Delta_{\cancel{PQ}}$ 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 $\Delta_{\cancel{PQ}}$ 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
$\Delta_{\cancel{PQ}}\gtrsim 10$ for $f_{a}\sim 10^{9}\,$GeV, and
$\Delta_{\cancel{PQ}}\gtrsim 16$ for $f_{a}\sim 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
$\Delta_{\cancel{PQ}}=9$ and $\Delta_{\cancel{PQ}}=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. (8) on the value of the PQ scale and the possibility of a
(loop) suppression of the effective coefficient $c$ and of the $CP$-violating
phase $\phi_{CP}$.

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 $\Lambda_{\rm 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 $\Lambda_{\rm 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_{\rm 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
$\Lambda_{\rm DC}$ is characterized by a strongly-coupled flow lingering close
to a pair of fixed points at complex coupling Kaplan:2009kr ; Gorbenko:2018ncu
. 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,777One could however use supersymmetry to
gain analytic control over the strong dynamics, see for example Ref.
Nakai:2021nyf  for a model of this kind. 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. Cox:2019rro .

###  3.1 General remarks on PQ violation: a spurion analysis

As already noticed by Ref. Dobrescu:1996jp , 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 ${\rm SU}(N_{\rm DC})\times{\rm
G}_{\rm W}\times{\rm G}_{\rm SM}$ gauge theory with Weyl fermions $\psi$ and
$\chi$ transforming as vectorlike representations of ${\rm SU}(N_{\rm
DC})\times{\rm G}_{\rm SM}$.

A generic PQ-violating operator ${\cal O}_{\cancel{PQ}}$ induces a potential
for the axion if the matrix element $\langle\psi_{a,\varphi}|{\cal
O}_{\cancel{PQ}}|0\rangle$ is non vanishing. Here $|\psi_{a,\varphi}\rangle$
is an arbitrary state containing axions and possibly other SM-neutral NGBs,
denoted by $\varphi$, with zero momentum.888While 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 $\varphi$’s. Therefore, ${\cal O}_{\cancel{PQ}}$ is an interpolating
operator for the state $|\psi_{a,\varphi}\rangle$ and must have its same
quantum numbers.

In order to determine the form of ${\cal O}_{\cancel{PQ}}$ implied by this
requirement, it is useful to first neglect the effects of the weak gauging of
${\rm G}_{\rm W}\times{\rm G}_{\rm SM}$. In this limit, the global symmetry is
$G_{F}=\,{\rm SU}(N_{f})_{L}\times\,{\rm SU}(N_{f})_{R}\times\,{\rm
U}(1)_{V}$, under which the fermions transform as

$\psi^{i_{L}}\rightarrow
e^{i\alpha}\,L^{i_{L}}_{j_{L}}\,\psi^{j_{L}},\quad\quad\quad\chi_{i_{R}}\rightarrow
e^{-i\alpha}\,(R^{*})_{i_{R}}^{j_{R}}\,\chi_{j_{R}}\,.$ (9)

Any $G_{F}$-violating operator made of the $\psi$ and $\chi$ fields can thus
be characterized by a spurionic tensorial structure

$T^{\\{i_{L}\\}\\{j_{R}\\}}_{\,\\{k_{L}\\}\\{l_{R}\\}}\,,$ (10)

where $\\{n\\}$ denotes a set of $n$ indices. At energies below the ${\rm
SU}(N_{\rm DC})$ confinement scale, $G_{F}$ is spontaneously broken to $\,{\rm
SU}(N_{f})_{V}\times\,{\rm U}(1)_{V}$, and the corresponding NGBs are
conveniently described by the field $\Sigma=\exp(i\pi^{\hat{a}}T^{\hat{a}})$.
The latter transforms as $\Sigma\to L\Sigma R^{\dagger}$ and thus carries one
upper index of $\,{\rm SU}(N_{f})_{L}$ plus one lower index of $\,{\rm
SU}(N_{f})_{R}$. A state of zero-momentum NGBs can be excited by the chiral
field $\Sigma$ and, more in general, by effective operators in the low-energy
theory that are polynomials of $\Sigma$ and $\Sigma^{\dagger}$. Since we are
interested in states with axions and light NGBs, we need to consider
exponentials of their fields. Because $\Sigma$ is unitary, each one of its
indices will contribute to determine the spurionic structure of the
interpolating effective operator (no $G_{F}$ singlet contractions are
possible). The latter must therefore have an equal number of upper left and
lower 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 $\,{\rm U}(1)_{V}$ charge, i.e. be a ‘mesonic’
operator. This implies that ${\cal O}_{\cancel{PQ}}$ is necessarily a
polynomial of the ${\rm SU}(N_{\rm DC})$-singlet bilinears
$(\psi_{r}\chi_{\bar{r}})$ and $(\psi_{r}\chi_{\bar{r}})^{\dagger}$ (whose
flavor indices are left understood), where $\psi_{r}$ and $\chi_{r}$ denote
fields transforming as the irreducible SM representation $r$. Such fermionic
bilinears are exactly those acquiring a vacuum expectation value at the ${\rm
SU}(N_{\rm DC})$ confinement scale. The same conclusion was obtained by
Dobrescu in Ref. Dobrescu:1996jp .

Let us now consider the effects of the weak gauging of ${\rm G}_{\rm
W}\times{\rm G}_{\rm SM}$. The latter breaks explicitly the flavor group
$G_{F}$ down to $[{\rm U}(1)_{V}]^{n_{irr}}\times[{\rm U}(1)_{A}]^{n_{irr}}$,
as we have already discussed, where the axial $[{\rm U}(1)_{A}]^{n_{irr}}$ is
non-linearly realized. The state $|\psi_{a,\varphi}\rangle$ is neutral under
all the vectorial ${\rm U}(1)$’s and under the SM group, while it transforms
non-trivially under the PQ symmetry. Since the same quantum numbers must also
characterize ${\cal O}_{\cancel{PQ}}$, we conclude that

A gauge-invariant PQ-violating operator ${\cal O}_{\cancel{PQ}}$ can generate
a potential for the axion at low energy only if it has vanishing vectorial
charges.

Operators satisfying this requirement are polynomials of the bilinears
$(\psi_{r}\chi_{\bar{r}})$, $(\psi_{r}\chi_{\bar{r}})^{\dagger}$,
$(\psi^{\dagger}_{r}\psi_{r})$ and $(\chi^{\dagger}_{\bar{r}}\chi_{\bar{r}})$.
Hence, operators built of $(\psi^{\dagger}_{r}\psi_{r})$ and
$(\chi^{\dagger}_{\bar{r}}\chi_{\bar{r}})$ factors are also allowed as a
consequence of the weak gauging.999Notice that $(\psi^{\dagger}_{r}\psi_{r})$
and $(\chi^{\dagger}_{\bar{r}}\chi_{\bar{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_{\rm SM}$, then a new operator can be
constructed with lower dimensionality where each bilinear is replaced by a
covariant derivative. 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.101010Here $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. In the limit of vanishing up and down
quark masses, $\pi^{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 $\pi^{0}$. More in detail, one can see that in the
low-energy theory $O$ gives rise to non-derivative interactions among
$\pi^{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 theory that generate a
potential for $\pi^{0}$. Another instructive example is that of the operator
$O=(uu^{c})(u^{\dagger}\bar{\sigma}^{\mu}u)^{2}$ in two-flavor QCD with
vanishing quark masses plus electromagnetism. By conveniently parametrizing
$\Sigma=\exp[i(\pi^{+}\sigma^{+}+\pi^{-}\sigma^{-})]\exp(i\pi^{0}\sigma^{3})$,
it is simple to see that at low energy
$O\sim\cos\pi^{0}\cos\sqrt{\pi^{+}\pi^{-}}\big{(}\partial_{\mu}\pi^{0}+i\pi^{+}\\!\overleftrightarrow{\partial}_{\\!\\!\mu}\pi^{-}/2+\dots\big{)}^{2}$.
A potential for $\pi^{0}$ 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^{\dagger}u)$.

So far we have implicitly assumed that the axion potential is generated by
single insertions of a local PQ-violating operator. However, multiple
insertions of local operators are also relevant and their importance was not
emphasized in the previous literature. On dimensional grounds, the insertion
of $N$ local operators of dimensions $d_{i}$ gives a contribution equivalent
to the insertion of a single operator of effective dimension

$d_{\text{eff}}=\sum_{i=1}^{N}d_{i}-4(N-1)\,.$ (11)

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 $(\psi_{r}\chi_{\bar{r}})$,
$(\psi_{r}\chi_{\bar{r}})^{\dagger}$, $(\psi^{\dagger}_{r}\psi_{r})$ and
$(\chi^{\dagger}_{\bar{r}}\chi_{\bar{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 ${\rm G}_{\rm 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 higher-dimensional 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 $\Delta_{\cancel{PQ}}$ 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 $(\psi_{r}\chi_{\bar{r}})$,
$(\psi_{r}\chi_{\bar{r}})^{\dagger}$, $(\psi^{\dagger}_{r}\psi_{r})$,
$(\chi^{\dagger}_{\bar{r}}\chi_{\bar{r}})$; and

  * (2)

the effective dimensions of all gauge-invariant, PQ-violating polynomials of
$(\psi_{r}\chi_{\bar{r}})$, $(\psi_{r}\chi_{\bar{r}})^{\dagger}$,
$(\psi^{\dagger}_{r}\psi_{r})$, $(\chi^{\dagger}_{\bar{r}}\chi_{\bar{r}})$
that can be formed as the product of gauge-invariant, flavor-violating 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 $\Delta_{\cancel{PQ}}$ defined as above satisfies the
bound of Eq. (8). In particular, we will require $\Delta_{\cancel{PQ}}\geq 9$
for $f_{a}\sim 10^{9}\,$GeV and $\Delta_{\cancel{PQ}}\geq 12$ for $f_{a}\sim
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 $\Delta_{\cancel{PQ}}$ if they have dimension larger than
$(\Delta_{\cancel{PQ}}+4)/2$. The calculation of $\Delta_{\cancel{PQ}}$ 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. Contino:2020god . With reference to the general structure of Table 1,
the most general charge assignments that are vector-like under ${\rm
SU}(N_{\rm DC})\times{\rm G}_{\rm SM}$ and overall chiral in the case of
irreducible $r_{i}$ are as follows: $p_{1}=-p_{2}=1$, $q_{2}=-q_{1}=q$, with
$r_{1}=r_{2}\equiv r$ and $q\in[0,1)$ (Type I); $p_{1}=-p_{2}=1$,
$q_{2}=-q_{1}=q$, with $r_{1}=\bar{r}_{2}\equiv r$ and $q\in[-1,1)$ (Type II);
$p_{1}=-q_{2}=1$, $p_{2}=-q_{1}=q$, with $r_{1}=\bar{r}_{2}\equiv r$ and
$q\in[0,1)$ (Type III). The only phenomenologically viable models consistent
with values of the confinement scale $\Lambda_{\rm DC}\sim 1-100\,\rm TeV$
were identified in Ref. Contino:2020god  to be those with $\,{\rm
SU}(3)_{c}\times\,{\rm SU}(2)_{EW}\times\,{\rm U}(1)_{Y}$ representations
$r=(\mathbf{1},\mathbf{2})_{0}$ and $r=(\mathbf{3},\mathbf{1})_{0}$.

| ${\rm SU}(N_{\rm DC})$ | ${\rm U}(1)_{\rm D}$ | ${\rm G}_{\rm SM}$
---|---|---|---
$\psi_{1}$ | $\square$ | $p_{1}$ | $r_{1}$
$\psi_{2}$ | $\square$ | $p_{2}$ | $r_{2}$
$\chi_{1}$ | $\bar{\square}$ | $q_{1}$ | $\bar{r}_{1}$
$\chi_{2}$ | $\bar{\square}$ | $q_{2}$ | $\bar{r}_{2}$
Table 1: General structure of accidental axion models with $n_{f}=2$. The
possible charge assignments $\\{p_{1},p_{2},q_{1},q_{2}\\}$ are restricted by
the cancellation of gauge anomalies as discussed in the main text. The group
${\rm G}_{\rm SM}$ is either the SM gauge group or one of its GUT extensions,
and the possible representations $r_{i}$ are classified in Section 4.3.

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 $\,{\rm 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. Contino:2020god . Indeed, it contains
two light pseudo NGBs, one of which is an axion-like particle. The
corresponding PQ symmetry does not commute with $\,{\rm SU}(5)$ and is
explicitly broken by dimension-6 four-fermion operators generated at the GUT
scale. Such effects induce an extra contribution to the potential of the
axion-like NGB that dominates over the QCD one for values of the strong
dynamical scale above $\sim 10^{8}\,$GeV, thus spoiling the resolution of the
strong CP problem. Lower values of the strong dynamical scale are instead
excluded by astrophysical bounds on axion-like particles.

When raising the dynamical scale to values above $\sim 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}\geq 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. Contino:2020god  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
$\Delta_{\cancel{PQ}}=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.

###  4.1 A Minimal QCD model

We consider a Type-II model with the charge assignment detailed in Table 2 and
$r=(\mathbf{3},\mathbf{1})_{+y}$ under ${\rm G}_{\rm SM}=\,{\rm
SU}(3)_{c}\times\,{\rm SU}(2)_{EW}\times\,{\rm U}(1)_{Y}$.

| ${\rm SU}(N_{\rm DC})$ | ${\rm U}(1)_{\rm D}$ | $\rm{\rm G}_{\rm SM}$
---|---|---|---
$\psi_{1}$ | $\square$ | $+1$ | $1$
$\psi_{2}$ | $\square$ | $-1$ | $1$
$\psi_{3}$ | $\square$ | $+1$ | $r$
$\psi_{4}$ | $\square$ | $-1$ | $\bar{r}$
$\chi_{1}$ | $\bar{\square}$ | $-q$ | $1$
$\chi_{2}$ | $\bar{\square}$ | $+q$ | $1$
$\chi_{3}$ | $\bar{\square}$ | $-q$ | $\bar{r}$
$\chi_{4}$ | $\bar{\square}$ | $+q$ | $r$
Table 2: Minimal accidental axion models with $n_{f}=2$. All the fields
$\psi_{i}$, $\chi_{i}$ are left-handed Weyl fermions, and the parameter $q$ is
a rational number in the interval $(-1,1)$. The QCD and GUT models discussed
in the text correspond to ${\rm G}_{\rm SM}=\,{\rm SU}(3)_{c}\times\,{\rm
SU}(2)_{EW}\times\,{\rm U}(1)_{Y}$ with $r=(\mathbf{3},\mathbf{1})_{+y}$, and
to ${\rm G}_{\rm SM}=\,{\rm SU}(5)_{\rm GUT}$ with $r=\mathbf{5}$,
respectively.

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.111111For non-vanishing hypercharge, Type I and Type
III models have mixed gauge anomalies, see Ref. Contino:2020god . They differ
from Type-II models only for the spectrum of heavy pseudo NGBs and other
resonances with mass of order $\Lambda_{\rm 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.121212For $q=-1$ the fields singlet under ${\rm G}_{\rm SM}$ form
vectorlike representations. Therefore, differently from the case with
$n_{f}=n_{irr}=2$ studied in Ref. Contino:2020god , this value of $q$ must be
discarded.

At the classical level the model has a ${\rm U}(1)_{L}^{4}\times{\rm
U}(1)_{R}^{4}$ global symmetry, corresponding to the group of phase
transformations of each of the fields $\psi_{i}$ and $\chi_{i}^{c}$
($\chi^{c}$ is the conjugate of $\chi$ and a right-handed Weyl fermion). At
confinement the pattern of symmetry breaking is

${\rm U}(1)_{L}^{4}\times{\rm U}(1)_{R}^{4}\longrightarrow{\rm
U}(1)_{V}^{4}\,,$ (12)

giving four spontaneously broken ${\rm U}(1)$ axial factors. One of them is
anomalous under ${\rm SU}(N_{\rm DC})$ and one is gauged by ${\rm U}(1)_{\rm
D}$. We are thus left with two broken generators, associated with two Nambu-
Goldstone bosons. These are neutral under the unbroken vectorial ${\rm
U}(1)_{V}$ symmetries, and will be denoted by $s$ and $a$ in the following.

Defining $\Psi_{L}=(\psi_{1}\,\psi_{2}\,\psi_{3}\,\psi_{4})$ and
$\Psi_{R}=(\chi^{c}_{1}\;\chi^{c}_{2}\;\chi^{c}_{3}\;\chi^{c}_{4})$, $s$ and
$a$ are interpolated by the axial currents

$J_{s}^{\mu}=\overline{\Psi}\gamma^{\mu}\gamma^{5}Q_{s}\Psi\,,\qquad
J_{a}^{\mu}=\overline{\Psi}\gamma^{\mu}\gamma^{5}Q_{a}\Psi\,,$ (13)

with

$Q_{s}={\rm diag}(3,-3,-\mathbb{1},\mathbb{1})\,,\qquad Q_{a}={\rm
diag}(-3,-3,\mathbb{1},\mathbb{1})\,.$ (14)

The current $J^{\mu}_{a}$ has an anomaly in the background of SM gauge fields:

$\langle\partial_{\mu}J_{a}^{\mu}\rangle=-\dfrac{\alpha_{3}}{4\pi}A_{3}\,G\tilde{G}-\dfrac{\alpha_{Y}}{4\pi}A_{Y}\,B\tilde{B},$
(15)

where $G_{\mu\nu}$ and $B_{\mu\nu}$ 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 as
$\tilde{F}^{\mu\nu}=\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}/2$. From now on
we suppress indices, writing $F\tilde{F}$ in place of
$F^{a}_{\mu\nu}\tilde{F}^{a,\mu\nu}$. The anomaly coefficients are

$\begin{split}A_{3}\,\delta_{bc}&=N_{\rm DC}\,{\rm
Tr}\left[\\{t_{b},t_{c}\\}Q_{a}\right]=2N_{\rm DC}\,\delta_{bc}\\\\[5.69046pt]
A_{Y}&=N_{\rm DC}\,{\rm Tr}\left[\\{Y,Y\\}Q_{a}\right]=12y^{2}N_{\rm
DC}\,,\end{split}$ (16)

where $t_{i}$ and $Y$ are the QCD and hypercharge generators respectively. We
can thus identify ${\rm 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, $\rm 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.131313For example,
one can think of weakly gauging the group $\rm U(1)_{V}$ generated by the
vectorlike charge $Q_{V}={\rm diag}(1,0,0,\mathbb{1})$, which is free from
gauge anomalies. The mixed $\rm[U(1)_{V}]^{2}U(1)_{s}$ anomaly is reproduced
at low energy by $s$.

At the level of discrete symmetries, one can define the action of charge
conjugation $C$ and parity $P$ on the new fields as follows141414This
definition corresponds to standard transformation rules for the Dirac field
$\Psi$: $\Psi(x)\to-i\gamma^{2}\Psi^{*}(x)$ under $C$ and
$\Psi(\vec{x},t)\to\gamma^{0}\Psi^{*}(-\vec{x},t)$ under $P$.

$\displaystyle C:\begin{cases}\psi_{i}(x)\leftrightarrow\chi_{i}(x)\\\ {\cal
G}^{a}_{\mu}(x)T_{a}\rightarrow-{\cal
G}^{a}_{\mu}(x)T_{a}^{*}\,,\end{cases}\qquad
P:\begin{cases}\psi_{i}(\vec{x},t)\leftrightarrow\chi^{c}_{i}(-\vec{x},t)\\\
{\cal G}^{a}_{\mu}(\vec{x},t)\rightarrow\eta_{\mu\mu}{\cal
G}^{a}_{\mu}(-\vec{x},t)\,,\end{cases}$ (17)

where ${\cal G}^{a}_{\mu}$ denotes the ${\rm SU}(N_{\rm DC})$ gauge fields,
$T_{a}$ are the corresponding generators, and $\eta_{\mu\mu}$ is equal to $+1$
for $\mu=0$ and $-1$ for $\mu=1,2,3$. Both $C$ and $P$ defined in this way are
explicitly broken by the chiral ${\rm U}(1)_{\rm 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^{\prime}$,
which is exact in the new sector and broken only by the EW interactions of SM
fields. We define

$C^{\prime}:\begin{cases}\psi_{1}\leftrightarrow\psi_{2}\,,\quad\psi_{3}\leftrightarrow\psi_{4}\,,\quad\chi_{1}\leftrightarrow\chi_{2}\,,\quad\chi_{3}\leftrightarrow\chi_{4}\\\
A_{\mu}^{D}\rightarrow-A_{\mu}^{D}\,,\quad
G^{a}_{\mu}\lambda_{a}\rightarrow-G^{a}_{\mu}\lambda_{a}^{*}\,,\quad
B_{\mu}\rightarrow-B_{\mu}\,,\end{cases}$ (18)

whereas the ${\rm SU}(N_{\rm DC})$ gauge fields do not transform. Under
$C^{\prime}$ 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 ${\rm U}(1)_{\rm D}$ gauge boson acquires a mass

$m_{\gamma_{D}}=2(1-q)e_{D}f_{\rm DC}\,.$ (19)

The pseudo NGBs, whose corresponding global symmetries are explicitly broken
by the weak gauging, acquire a mass from loops of SM and/or ${\rm U}(1)_{\rm
D}$ gauge bosons of order $\delta m^{2}\sim(\alpha_{i}/4\pi)\Lambda_{DC}^{2}$,
where the couplings are evaluated at $\Lambda_{\rm DC}$. Other baryonic and
mesonic resonances of the strong dynamics have mass of order $\Lambda_{\rm
DC}$. If $\Lambda_{\rm DC}$ is identified with the Peccei-Quinn scale, all
these states are extremely heavy; for values of the ${\rm U}(1)_{\rm D}$
coupling $\alpha_{D}\gtrsim 10^{-2}$ and $q\sim O(1)$, their masses are always
larger than $10^{8}\,\rm GeV$. All resonances can be classified according to
their transformation rules under the unbroken vectorial ${\rm U}(1)$
symmetries and some of them are accidentally stable. They are however out of
the reach of current experiments and can be neglected if their cosmological
abundance is zero.

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:

$\mathcal{L}_{\rm
eff}=\dfrac{1}{2}(\partial_{\mu}a)^{2}+\dfrac{1}{2}(\partial_{\mu}s)^{2}+\dfrac{\alpha_{3}}{8\pi}\dfrac{a}{f_{a}}\,G\tilde{G}+c_{B}\dfrac{\alpha_{Y}}{8\pi}\,\dfrac{a}{f_{a}}\,B\tilde{B},$
(20)

with

$f_{a}=\frac{f_{\rm DC}}{2A_{3}}=\dfrac{f_{\rm DC}}{4N_{\rm DC}}\,,\qquad
c_{B}=6y^{2}\,.$ (21)

From the latter equation we estimate $\Lambda_{\rm DC}\sim 4\pi f_{\rm
DC}/\sqrt{N_{\rm DC}}=16\pi\sqrt{N_{\rm DC}}f_{a}$ in this model.151515Notice
that the value of the PQ decay constant $f_{\rm DC}$ depends on how it is
defined. A redefinition of the PQ charges, for example, also changes $f_{\rm
DC}$ if the latter is defined in terms of the matrix element of the PQ
current. The dark color scale $\Lambda_{\rm DC}$ and the axion decay constant
$f_{a}$ set, respectively, the mass of the dark hadrons and the axion coupling
to photons. As such, they do not depend on a redefinition of the PQ charges.
This means that they are related to $f_{\rm DC}$ in a way that depends on the
convention chosen. The naive estimate $\Lambda_{\rm DC}\sim 4\pi f_{\rm
DC}/\sqrt{N_{\rm DC}}$ does not take this dependence into account and is thus
necessarily a very rough approximation. Knowing that the number of dark
flavors is $N_{f}=8$ and the effective number of additional QCD flavors is
$\Delta N_{f}^{(\rm QCD)}=2N_{\rm DC}$, we can derive the constraints implied
by the requirements of dark confinement and perturbativity of SM couplings. We
demand that the theory remains perturbative for values of $f_{a}$ as low as
$10^{9}\,$GeV. Equations (3) and (5) then imply the bound on the number of
dark colors

$3\leq N_{\rm DC}\leq 13\,,$ (22)

where $N_{\rm DC}=3$ is at the edge of the conformal window DeGrand:2015zxa .
Furthermore, Eq. (6) gives the bound $y^{2}\lesssim 0.1$ for $N_{\rm DC}$ in
the range of Eq. (22).

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_{\rm{DW}}$, is equal to the anomaly coefficient computed
in the basis where all the PQ charges are coprime integers. In our notation,

$N_{\rm{DW}}=|A_{3}|=2\,N_{\rm{DC}}.$ (23)

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 find $E/N=c_{B}=6y^{2}$. Higher-
derivative interactions, such as $B^{\mu\nu}\partial_{\mu}s\partial_{\nu}a$,
$(\partial_{\mu}s)^{2}B_{\rho\sigma}B^{\rho\sigma}$,
$(\partial_{\mu}a)^{2}(\partial_{\nu}s)^{2}$ etc., are suppressed by higher
powers of $f_{\rm DC}$ and can be safely neglected at the energies we are
considering. By matching the renormalizable UV dynamics to the effective
theory at $\Lambda_{\rm DC}$, no derivative operator of the form
$\partial_{\mu}aj^{\mu}$, where $j^{\mu}$ is an axial current made of SM
fermions, is generated at leading order in the SM couplings. This is because
SM fermions communicate with the new fermions only through SM gauge fields.
The same feature characterizes many models of composite axions. Of course,
since $\partial_{\mu}j^{\mu}$ mixes with $G\tilde{G}$ and $B\tilde{B}$ under
the renormalization group, operators $\partial_{\mu}aj^{\mu}$ are generated
through loop diagrams with one insertion of the anomalous couplings in Eq.
(20). Furthermore, non-renormalizable operators of the form

$\bar{\Psi}\gamma_{\mu}\gamma^{5}\,\Psi\bar{\Psi}_{SM}\gamma^{\mu}\gamma^{5}\Psi_{SM}$
(24)

can also generate $\partial_{\mu}aj^{\mu}$, as well as
$\partial_{\mu}sj^{\mu}$, at the matching scale if present in the UV theory.

Since the Peccei-Quinn symmetry – identified with ${\rm 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 $\psi_{1}\psi_{2}\chi_{1}\chi_{2}$ and
$\psi_{3}\psi_{4}\chi_{3}\chi_{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.

###  4.2 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 ${\rm 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}\leq 1.2\times 10^{16}\,$GeV is
required (according to an estimate based on the results of Ref.
Hebecker:2018ofv ) to protect the axion from too large non-perturbative
gravitational effects, and GUT scales smaller than $\sim 10^{16}\,$GeV are
generically at odds with proton decay.

Extending the theory is also motivated by its chiral nature: if $M_{\rm
GUT}>\Lambda_{\rm DC}$ and the GUT dynamics does not break $\,{\rm SU}(N_{\rm
DC})\times{\rm U}(1)_{\rm 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 Georgi:1974sy , this can
be realized by the straightforward embedding of the color triplets into
fundamentals of ${\rm SU}(5)_{\rm GUT}$, as shown in Table 2. With this
choice, the number of ${\rm SU}(N_{\rm DC})$ flavors is $N_{f}=12$, while the
effective number of additional QCD flavors is unchanged. As reported in Eq.
(29), the relation between $f_{a}$ and $f_{\rm DC}$ is also unchanged. The
conditions of confinement and perturbativity, Eqs. (3) and (5), then imply on
$N_{\rm DC}$ the same constraint of Eq. (22) as in the QCD model. As before,
the renormalizable lagrangian includes only kinetic terms: no Yukawa terms
with GUT scalars are allowed by gauge and Lorentz invariance if these fields
do not carry ${\rm U}(1)_{\rm 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
$\Lambda_{DC}$, hence the symmetry breaking pattern of Eq. (12) is enlarged to

${\rm U}(1)_{L}^{6}\times{\rm U}(1)_{R}^{6}\longrightarrow{\rm
U}(1)_{V}^{6}\,,$ (25)

resulting in two additional pseudo NGBs. These correspond to the two SM
singlets contained in ${\rm U}(1)_{\rm D}$-neutral adjoints of ${\rm
SU}(5)_{\rm GUT}$ in the spectrum of NGBs. We will denote them by $\tilde{a}$
and $\tilde{s}$ and indicate with ${\rm U}(1)_{\tilde{a}}$ and ${\rm
U}(1)_{\tilde{s}}$ the associated spontaneously broken symmetries. Both ${\rm
U}(1)_{\tilde{a}}$ and ${\rm U}(1)_{\tilde{s}}$ are explicitly broken by the
GUT dynamics, and a mass for $\tilde{a}$ and $\tilde{s}$ of order

$m_{\tilde{a},\tilde{s}}^{2}\sim\frac{g^{2}_{\rm
GUT}}{16\pi^{2}}\frac{\Lambda_{\rm DC}^{4}}{M^{2}_{\rm GUT}}$ (26)

is generated radiatively by loops of GUT gauge bosons. The same effect arises
in the low energy theory from dimension-6 operators that are generated at
$M_{\rm GUT}$ by integrating out GUT states. We define the generators so that
the ${\rm U}(1)_{\tilde{a}}$ group is anomalous under $\,{\rm
SU}(3)_{c}\times\,{\rm SU}(2)_{EW}\times\,{\rm U}(1)_{Y}$, hence $\tilde{a}$
is a heavy axion-like particle, while ${\rm U}(1)_{\tilde{s}}$ is free from
gauge anomalies. Under $CP$ $\tilde{a}$ and $\tilde{s}$ are both odd, while
under $C^{\prime}$ they have charge $+1$ and $-1$ respectively. Since
$C^{\prime}$ is an exact symmetry in the dark sector, the mass matrix of
$\tilde{s}$ and $\tilde{a}$ is diagonal up to higher-order
$C^{\prime}$-breaking effects. A mixing between $\tilde{s}$ and $\tilde{a}$
can only arise through higher order loops involving SM fermions.

According to Eq. (26), $\tilde{a}$ and $\tilde{s}$ have masses that range from
$\sim 10^{4}\,$GeV to $\sim 10^{8}\,$GeV when $\Lambda_{\rm DC}$ varies from
$10^{11}\,$GeV to $10^{13}\,$GeV. Their only interactions come from non-
renormalizable operators and are too feeble to give appreciable signals in
laboratory experiments. Similarly to the axion, both $\tilde{a}$ and
$\tilde{s}$ can be very long-lived: $\tilde{a}$ can decay to two photons
through the anomaly, with decay width

$\Gamma_{\tilde{a}\to\gamma\gamma}\sim\frac{\alpha_{em}^{2}}{2(4\pi)^{3}}\frac{m^{3}_{\tilde{a}}}{f_{\rm
DC}^{2}}\,,$ (27)

while $\tilde{s}$ has an even smaller decay width that can be induced only by
a mixing with $\tilde{s}$ generated by higher-order $C^{\prime}$-breaking
loops. The metastability of $\tilde{a}$ and $\tilde{s}$ can have an impact on
cosmology, as discussed in Section 6.

Apart from $\tilde{a}$ and $\tilde{s}$, 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. (14) is straightforward,

$Q_{s}={\rm diag}(5,-5,-\mathbb{1},\mathbb{1}),\qquad Q_{a}={\rm
diag}(-5,-5,\mathbb{1},\mathbb{1})\,,$ (28)

and the same properties for $a$ and $s$ are derived as for the QCD model. At
energies of order $\Lambda_{\rm DC}$ the theory can be matched to an effective
Lagrangian for the four SM-singlet NGBs. This is similar to that of Eq. (20),
with $f_{a}$ and $N_{DW}$ given again by

$f_{a}=\dfrac{f_{\rm DC}}{4N_{DC}}\qquad\quad N_{DW}=2N_{DC},$ (29)

but with two additions. First, there is an anomalous coupling of the axion
also to $SU(2)_{EW}$ gauge fields and, as a consequence of the $SU(5)$
invariance, the axion-photon coupling takes the standard value for a model
with complete ${\rm SU}(5)_{\rm GUT}$ representations, i.e. $E/N=8/3$. Second,
four-fermion operators of the same form as in Eq. (24) are generated at
$M_{\rm GUT}$ by the tree-level exchange of heavy gauge bosons; at the scale
$\Lambda_{\rm DC}$, they in turn generate derivative interactions of the form
$c_{j}\partial_{\mu}aj^{\mu}/f_{a}$ with coefficients

$c_{j}\sim\dfrac{g_{\rm GUT}^{2}f_{a}^{2}}{M^{2}_{\rm{GUT}}}\,.$ (30)

For $g_{\rm{GUT}}\lesssim 0.1$ and $\Lambda_{\rm{DC}}<M_{\rm{GUT}}$, the
estimated value of $c_{j}$ is always smaller than $10^{-4}$.

###  4.3 Generic Models with $n_{f}=2$

We now generalize the classification to the case of arbitrary ${\rm G}_{\rm
SM}$ representations $r_{i}=\sum_{\alpha}r_{i}^{(\alpha)}$, with
$r_{i}^{(\alpha)}$ irreducible. To this aim, it is convenient to introduce the
following short-hand notation:

$d_{i}={\rm dim}(r_{i})=\sum_{\alpha}{\rm dim}(r_{i}^{(\alpha)})\,,\qquad
T_{i}=\sum_{\alpha}T(r_{i}^{(\alpha)})\,.$ (31)

The choice of ${\rm U}(1)_{\rm D}$ charges is restricted by the requirement of
gauge anomaly cancellation. Linearly independent, non-trivial constraints come
from the anomalies $[{\rm G}_{\rm SM}]^{2}{\rm U}(1)_{\rm D}$, $[{\rm
SU}(N_{\rm DC})]^{2}{\rm U}(1)_{\rm D}$ and $[{\rm U}(1)_{\rm D}]^{3}$;
following the notation of Eq. (1), the corresponding conditions are:

$\begin{cases}(p_{1}+q_{1})T_{1}+(p_{2}+q_{2})T_{2}=0,\\\
(p_{1}+q_{1})d_{1}+(p_{2}+q_{2})d_{2}=0,\\\
(p_{1}^{3}+q_{1}^{3})d_{1}+(p_{2}^{3}+q_{2}^{3})d_{2}=0\,.\end{cases}$ (32)

Combining the first two equations we obtain

$(p_{1}+q_{1})\left(d_{2}T_{1}-d_{1}T_{2}\right)=0\,.$ (33)

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 $\,{\rm SU}(N)$ with rank 2 or lower, as done in Ref.
Contino:2020god  to maintain perturbativity up to the Planck scale. In that
case, the possible charge assignments are the same as those of Ref.
Contino:2020god , which lead to Type I, II and III theories defined at the
beginning of this section. For reducible $r_{i}$, new solutions of Eqs. (32)
and (33) exist with $d_{1}\neq d_{2}$; in the case of a chiral assignment of
charges and for given arbitrary $q_{1}$ and $p_{1}$, we find

$\begin{split}q_{2}=-\dfrac{1}{2}\dfrac{d_{1}}{d_{2}}(p_{1}+q_{1})\pm\dfrac{\sqrt{12d_{2}^{4}(p_{1}^{2}-p_{1}q_{1}+q_{1}^{2})-3d_{1}^{2}d_{2}^{2}(p_{1}+q_{1})^{2}}}{6d_{2}^{2}},\\\\[5.69046pt]
p_{2}=-\dfrac{1}{2}\dfrac{d_{1}}{d_{2}}(p_{1}+q_{1})\mp\dfrac{\sqrt{12d_{2}^{4}(p_{1}^{2}-p_{1}q_{1}+q_{1}^{2})-3d_{1}^{2}d_{2}^{2}(p_{1}+q_{1})^{2}}}{6d_{2}^{2}}.\end{split}$
(34)

It is always possible to redefine the ${\rm U}(1)_{\rm 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 ${\rm G}_{\rm 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_{\rm 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. (5) we assume a
relation $\Lambda_{\rm 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_{\rm DC}$ and $f_{a}$ (by determining the embedding of ${\rm
U}(1)_{\rm PQ}$ in the global symmetry group) and varying $f_{a}$ in the whole
allowed range. Some models with larger ratios $f_{\rm 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.

${\rm G}_{\rm SM}$ | $r_{1}$ | $r_{2}$ | $d_{1}$ | $d_{2}$ |  | $\Delta^{\rm max}_{\cancel{PQ}}$
---|---|---|---|---|---|---
$\,{\rm SU}(3)_{\rm c}$ | $(\mathbf{3}\oplus m\,\mathbf{1})$ | $(\mathbf{3}\oplus m\,\mathbf{1})$ | $3+m$ | $3+m$ | $1\leq m\leq 23$ | 6
| $(\mathbf{3}\oplus m\,\mathbf{1})$ | $2\,(\mathbf{3}\oplus m\,\mathbf{1})$ | $3+m$ | $2(3+m)$ | $1\leq m\leq 8$ | 9
| $(\mathbf{3}\oplus m\,\mathbf{1})$ | $3\,(\mathbf{3}\oplus m\,\mathbf{1})$ | $3+m$ | $3(3+m)$ | $1\leq m\leq 3$ | 12
| $2\,(\mathbf{3}\oplus m\,\mathbf{1})$ | $2\,(\mathbf{3}\oplus m\,\mathbf{1})$ | $2(3+m)$ | $2(3+m)$ | $1\leq m\leq 3$ | 6
| $(\mathbf{3}\oplus\mathbf{1})$ | $4\,(\mathbf{3}\oplus\mathbf{1})$ | $4$ | $16$ |  | 15
| $2\,(\mathbf{3}\oplus\mathbf{1})$ | $3\,(\mathbf{3}\oplus\mathbf{1})$ | $8$ | $12$ |  | 15
$\rm SU(5)_{GUT}$ | $(\mathbf{5}\oplus m\,\mathbf{1})$ | $(\mathbf{5}\oplus m\,\mathbf{1})$ | $5+m$ | $5+m$ | $1\leq m\leq 21$ | 6
| $(\mathbf{5}\oplus m\,\mathbf{1})$ | $2\,(\mathbf{5}\oplus m\,\mathbf{1})$ | $3+m$ | $2(3+m)$ | $1\leq m\leq 6$ | 9
| $(\mathbf{5}\oplus\mathbf{1})$ | $3\,(\mathbf{5}\oplus\mathbf{1})$ | $6$ | $18$ |  | 12
| $2\,(\mathbf{5}\oplus\mathbf{1})$ | $2\,(\mathbf{5}\oplus\mathbf{1})$ | $12$ | $12$ |  | 6
| $\mathbf{5}$ | $(\mathbf{10}\oplus 5\,\mathbf{1})$ | $5$ | $15$ |  | 12
Table 3: Choices of representations for models with $n_{f}=2$. The full list
of possibilities contains also theories that are obtained by replacing any of
the irreducible representations in this table with its conjugate. In the case
of the $\,{\rm SU}(3)_{\rm c}$ theories listed in the upper panel, electroweak
quantum numbers are not specified but can be assigned freely as long as the
conditions from cancellation of gauge anomalies and perturbativity of gauge
couplings are satisfied.

We still need to assess to what level the PQ symmetry is protected in these
theories. From Eq. (32) it follows that an operator of the form

$(\psi_{1}\chi_{1})^{d_{1}}(\psi_{2}\chi_{2})^{d_{2}}$ (35)

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 ${\rm gcd}(d_{1},d_{2})>1$, where
${\rm 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}\rightarrow d_{1}/{\rm gcd}(d_{1},d_{2})$ and $d_{2}\rightarrow
d_{2}/{\rm gcd}(d_{1},d_{2})$. 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 $\Delta_{\cancel{PQ}}$, as defined in Sec. 3.1, has
an upper bound $\Delta_{\cancel{PQ}}\leq\Delta^{\rm
max}_{\cancel{PQ}}=3(d_{1}+d_{2})/{\rm gcd}(d_{1},d_{2})$. 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 $\Delta^{\rm max}_{\cancel{PQ}}\geq 9$
have instead a chance to solve the quality problem, for suitable values of
$f_{a}$, depending on the actual value of $\Delta_{\cancel{PQ}}$. Determining
the latter requires specifying the ${\rm U}(1)_{\rm D}$ charges and the
embedding of ${\rm U}(1)_{\rm PQ}$ into the global symmetry group. We leave
such analysis for future work Podo:2022gyj .

##  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 ${\rm U}(1)_{\rm 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.

| ${\rm SU}(N_{\rm DC})$ | ${\rm U}(1)_{\rm D}$ | ${\rm G}_{\rm SM}$ | ${\rm U}(1)_{\rm PQ}$
---|---|---|---|---
$\psi_{1}$ | $\square$ | $p_{1}$ | $r_{1}$ | $\alpha$
$\psi_{2}$ | $\square$ | $p_{2}$ | $r_{2}$ | $\beta$
$\psi_{3}$ | $\square$ | $p_{3}$ | $r_{3}$ | $\gamma$
$\chi_{1}$ | $\bar{\square}$ | $q_{1}$ | $\bar{r}_{1}$ | $\alpha$
$\chi_{2}$ | $\bar{\square}$ | $q_{2}$ | $\bar{r}_{2}$ | $\beta$
$\chi_{3}$ | $\bar{\square}$ | $q_{3}$ | $\bar{r}_{3}$ | $\gamma$
Table 4: Composite axion model with $n_{f}=3$. The representations
$r_{1},r_{2},r_{3}$ are assumed to be irreducible, and $G_{\rm SM}$ is either
the SM gauge group or one of its GUT extensions. The ${\rm U}(1)_{\rm D}$
charges are defined so that $p_{i}\not=p_{j}$, $q_{i}\not=q_{j}$ for $i\not=j$
and can be normalized to be all integer numbers.

The conditions for the cancellation of $[{\rm G}_{\rm SM}]^{2}{\rm U}(1)_{\rm
D}$, $[{\rm SU}(N_{\rm DC})]^{2}{\rm U}(1)_{\rm D}$ and $[{\rm U}(1)_{\rm
D}]^{3}$ gauge anomalies are:

$\begin{cases}(p_{1}+q_{1})T_{1}+(p_{2}+q_{2})T_{2}+(p_{3}+q_{3})T_{3}=0,\\\
(p_{1}+q_{1})d_{1}+(p_{2}+q_{2})d_{2}+(p_{3}+q_{3})d_{3}=0,\\\
(p_{1}^{3}+q_{1}^{3})d_{1}+(p_{2}^{3}+q_{2}^{3})d_{2}+(p_{3}^{3}+q_{3}^{3})d_{3}=0.\\\
\end{cases}$ (36)

If the representations $r_{1},r_{2},r_{3}$ are all equal, then the accidental
NGB in the spectrum is not anomalous under SM color and cannot play the role
of the QCD axion. If two representations are equal (up to conjugation), then
the system (36) admits only vector-like solutions. Indeed, in that case from
the first two equations of Eq. (36) it follows
$(p_{i}+q_{i})(T_{i}d_{j}-T_{j}d_{i})=0$, where $r_{i}$ is the representation
which differs from the other two. The relation $T_{i}d_{j}=T_{j}d_{i}$ is not
satisfied by any pair of representations of $\,{\rm SU}(N)$ with rank 2 or
lower, while higher-dimensional representations are excluded by the
perturbativity constraint. This implies $p_{i}=-q_{i}$, hence vector-like
representations. Therefore, we are led to consider models where the three
irreducible representations are all different and inequivalent under complex
conjugation.

The choices of representations are restricted by the conditions of Section
2.1. We enforce the constraint of Eq. (5) by estimating $\Lambda_{\rm DC}=4\pi
f_{\rm DC}/\sqrt{N_{\rm DC}}$ and computing the ratio $f_{\rm DC}/f_{a}$ in
each model. Together with Eq. (3), this leads to a range of allowed values of
$N_{\rm DC}$ and $f_{a}$. The minimum number of dark colors is determined by
Eq. (3), from which it follows the minimum value of $f_{a}$ through Eq. (5).
We require $f_{a}\lesssim 10^{12}\,$GeV, and this in turn implies an upper
bound on $N_{\rm DC}$ again through Eq. (5). 161616We could consider larger
values of $f_{a}$, up to $M_{\rm GUT}$. For $f_{a}>10^{12}\,$GeV, on the other
hand, the axion energy density from misalignment is too large unless one
assumes small initial conditions, and furthermore the axion quality problem
exacerbates. According to these criteria, there are only a few possibile
representations $r_{1},r_{2},r_{3}$ leading to viable models. They are
reported in Table 5, together with the allowed range of dark colors and the
minimum value of $f_{a}$.

${\rm G}_{\rm SM}$ | $r_{1}$ | $r_{2}$ | $r_{3}$ | $\Delta^{\rm max}_{\cancel{PQ}}$ | $N_{\rm DC}$ | $f_{a}^{\rm min}$
---|---|---|---|---|---|---
${\rm SU}(5)_{\rm GUT}$ | $\mathbf{1}$ | $\mathbf{\bar{5}}$ | $\mathbf{10}$ | $12$ | $4,5,\dots,11$ |  $4\cdot 10^{8}\,$GeV
| $\mathbf{1}$ | $\mathbf{\bar{5}}$ | $\mathbf{15}$ | $15$ | $6,7$ |  $10^{11}\,$GeV
| $\mathbf{1}$ | $\mathbf{10}$ | $\mathbf{15}$ | $18$ | $7$ |  $10^{12}\,$GeV
$\,{\rm SU}(3)_{c}$ | $\mathbf{1}$ | $\mathbf{3}$ | $\mathbf{6}$ | $12$ | $3,4,\dots,9$ |  $4\cdot 10^{8}\,$GeV
| $\mathbf{8}$ | $\mathbf{3}$ | $\mathbf{6}$ | $15$ | $5$ |  $5\cdot 10^{11}\,$GeV
| $\mathbf{1}$ | $\mathbf{3}$ | $\mathbf{8}$ | $15$ | $3,4,\dots,7$ |  $10^{9}\,$GeV
| $\mathbf{1}$ | $\mathbf{6}$ | $\mathbf{8}$ | $12$ | $4,5$ |  $10^{11}\,$GeV
Table 5: Choices of the representations defined in Table 4 allowed by the
constraints of Section 2.1. Additional viable models are obtained by
conjugating the representation $r_{3}$ (or equivalently $r_{2}$), except for
models in the last two rows. The last two columns report, respectively, the
allowed values of $N_{\rm DC}$ and the smallest allowed 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 $\Delta_{\cancel{PQ}}$ 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 $\Delta_{\cancel{PQ}}$ by exploiting the conditions
of anomaly cancellation, Eq (36). From the first two equations, in particular,
it follows that operators of the form

$\prod_{i=1}^{3}\left(\psi_{i}\chi_{i}\right)^{\kappa_{i}}\left(\psi_{i}^{*}\chi_{i}^{*}\right)^{\bar{\kappa}_{i}}$
(37)

are neutral under ${\rm U}(1)_{\rm D}$ provided that 171717Here we use the
fact that $p_{i}\not=-q_{i}$ and $T_{i}d_{j}\not=T_{j}d_{i}$ for any $i,j$ as
discussed earlier.

$(\kappa_{1}-\bar{\kappa}_{1})(T_{2}d_{3}-T_{3}d_{2})+(\kappa_{2}-\bar{\kappa}_{2})(T_{3}d_{1}-T_{1}d_{3})+(\kappa_{3}-\bar{\kappa}_{3})(T_{1}d_{2}-T_{2}d_{1})=0\,.$
(38)

These operators are invariant also under ${\rm SU}(N_{\rm DC})\times{\rm
G}_{\rm SM}$, manifestly violate ${\rm U}(1)_{\rm PQ}$, and are of the right
form to generate an axion potential through a single insertion. Their
dimensionality is equal to $3\sum_{i}(\kappa_{i}+\bar{\kappa}_{i})$.
Therefore, for a given model of Tab. 5, there exists an upper bound on the
dimension of the PQ-violating operators given by $\Delta^{\rm
max}_{\cancel{PQ}}=3\kappa_{\rm min}$, where $\kappa_{\rm
min}=\min\\{\sum_{i}(\kappa_{i}+\bar{\kappa}_{i})\\}$ for positive integer
$\kappa_{i}$, $\bar{\kappa}_{i}$ satisfying Eq. (38). The values of
$\Delta^{\rm max}_{\cancel{PQ}}$ are reported in Tab. 5 for the selected
models.

All the models of Tab. 5 can address the quality problem if their value of
$\Delta_{\cancel{PQ}}$ equals the upper value $\Delta_{\cancel{PQ}}^{\rm max}$
computed as above. It is in fact plausible that this equality can be achieved
with a suitable choice of the ${\rm U}(1)_{\rm D}$ charges. In the following
we shall focus on GUT models with
$(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{\bar{5}},\mathbf{10})$, where $f_{a}$
can be as low as $4\cdot 10^{8}\,$GeV. We will present their relevant
properties and show that choices of ${\rm U}(1)_{\rm D}$ charges exist for
which $\Delta_{\cancel{PQ}}=\Delta_{\cancel{PQ}}^{\rm max}=12$. QCD models
with $(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{3},\mathbf{6})$ will be briefly
discussed in Appendix B.

###  5.1 Analysis of $(1,\bar{5},10)$ models

Let us consider models of Table 4 with representations
$(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{\bar{5}},\mathbf{10})$ or
$(\mathbf{1},\mathbf{5},\mathbf{10})$ of $\,{\rm SU}(5)$. In this case the
number of dark flavors is $N_{f}=16$, and the effective number of additional
QCD flavors is $\Delta N_{f}^{\rm(QCD)}=4N_{\rm DC}$. The conditions of Eqs.
(3), (5) restrict the number of dark colors to $4\leq N_{\rm DC}\leq 11$, and
values of $f_{a}$ as small as $4\cdot 10^{8}\,$GeV are allowed. The ${\rm
U}(1)_{\rm D}$ charges must be chosen so as to have chiral representations and
this implies that no gauge invariant mass terms are possible. However,
differently from the $n_{f}=2$ case, this condition alone does not
automatically forbid the existence of other gauge-invariant 2-fermion
operators featuring GUT scalars. These operators explicitly violate the
Peccei-Quinn symmetry and their impact on the axion solution should be
analyzed carefully. 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 ${\rm U}(1)_{\rm D}$
charges that do not allow any operator with two dark fermions. In practice,
this is equivalent to require $q_{i}\not=-p_{j}$ for any $i,j$, and
$p_{i}\not=p_{j}$, $q_{i}\not=q_{j}$ for any $i\not=j$.

From Section 2 we know that the number of NGBs that are singlets of ${\rm
G}_{\rm 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}=\text{diag}(5,\mathbb{1}_{5},-\mathbb{1}_{10})\,.$ (39)

The corresponding current $J^{\mu}_{a}$ has an anomaly in the background of
$\,{\rm SU}(5)$ gauge fields

$\langle\partial_{\mu}J_{a}^{\mu}\rangle=-\dfrac{\alpha_{5}}{4\pi}A_{5}\,G\tilde{G}\,,\qquad
A_{5}=-2N_{\rm DC}\,,$ (40)

where $G_{\mu\nu}$ is the $\,{\rm SU}(5)$ field strength. In this case the
axion decay constant and domain wall numbers are

$f_{a}=\frac{f_{\rm DC}}{4N_{\rm DC}}\qquad\quad N_{DW}=2N_{DC}\,.$ (41)

Besides the axion there are no additional SM-singlet NGBs that play the role
of $s$ in $n_{f}=2$ models. On the other hand, there are three pseudo NGBs,
whose mass can be estimated as in Eq. (26). Their number can be easily deduced
by noticing that the number of irreducible representations increases from
three to six when the GUT group is broken to the SM one. In particular, the
three pseudo NGBs originate from the $\mathbf{24}$ contained in the product
$\mathbf{5}\times\mathbf{\bar{5}}$, and from the $\mathbf{24}$ and
$\mathbf{75}$ contained in the product
$\mathbf{10}\times\mathbf{\overline{10}}$. We shall denote them by
$\varphi_{I}$, with $I=1,2,3$. In an orthonormal basis (i.e. ${\rm
Tr}(T_{I}T_{J})=\delta_{IJ}/2$), the corresponding generators can be written
as

$\begin{split}&T_{1}=\frac{1}{2\sqrt{5}}\,{\rm
diag}(0,3\cdot\mathbb{1}_{2},-2\cdot\mathbb{1}_{3},0,0\cdot\mathbb{1}_{3},0\cdot\mathbb{1}_{6})\\\
&T_{2}=\frac{1}{6}\,{\rm
diag}(0,0\cdot\mathbb{1}_{2},0\cdot\mathbb{1}_{3},0,-2\cdot\mathbb{1}_{3},\mathbb{1}_{6})\\\
&T_{3}=\frac{1}{6\sqrt{5}}\,{\rm
diag}(0,0\cdot\mathbb{1}_{2},0\cdot\mathbb{1}_{3},-9,\mathbb{1}_{3},\mathbb{1}_{6}).\\\
\end{split}$ (42)

Then, the anomalous couplings to the SM gauge fields can be written as

$\mathcal{L}\supset\dfrac{\varphi_{I}}{f_{\rm
DC}}A_{Ii}\dfrac{\alpha_{i}}{8\pi}\,G_{i}\,\tilde{G}^{i}\,,$ (43)

where $i$ runs over the three simple factors of the SM gauge group $\,{\rm
SU}(3)_{c}\times\,{\rm SU}(2)_{EW}\times\,{\rm U}(1)_{Y}$, and the anomaly
matrix has the form

$A_{Ii}=N_{\rm{DC}}\begin{pmatrix}-\frac{1}{\sqrt{15}}&\frac{3}{2\sqrt{15}}&\frac{5}{6\sqrt{15}}\\\
0&\frac{1}{2}&-\frac{5}{6}\\\
\frac{1}{2\sqrt{5}}&\frac{1}{2\sqrt{5}}&-\frac{5}{2\sqrt{5}}\end{pmatrix}.$
(44)

Since the matrix $A_{Ii}$ has rank $2$, one might be tempted to conclude that
there is a particular combination of the $\varphi_{I}$ which does not couple
to any of the anomalies, and is thus stable. However, the mass generated
radiatively by ${\rm SU}(5)_{\rm GUT}$ interactions for the two pseudo NGBs
contained in the $\mathbf{24}$ is different from that of the pseudo NGB coming
from the $\mathbf{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. Contino:2010rs . Therefore, the rotation required to decouple
one of the pseudo NGBs from the anomalies would induce a mixing between the
$\varphi_{I}$’s,181818We thank Luca Vecchi for bringing this point to our
attention. and we can conclude that all three states can decay through the
anomalies. From a qualitative viewpoint, the situation is very similar to that
of the $\tilde{a}$ in the previous section.

Finally, ${\rm U}(1)_{\rm D}$ is spontaneously broken at dark confinement and
the dark photon obtains a mass

$m_{\gamma_{D}}=\frac{4}{\sqrt{5}}|p_{1}+q_{1}|e_{D}f_{\rm{DC}}.$ (45)

To derive this formula we have used the anomaly conditions (36).

###  5.2 Explicit solutions and axion quality for $(1,\bar{5},10)$ models

To fully specify a model, one must assign ${\rm U}(1)_{\rm D}$ charges that
satisfy the conditions of anomaly cancellation in Eq. (36). The choice of
charges determines the structure of higher-dimensional 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 (36) does not appear to be an
easy problem and it will be tackled in a separate work Podo:2022gyj . However,
some explicit solutions are easy to find. A particularly interesting family of
solutions is given by

$(p_{1},p_{2},p_{3})=(5,-3,1)\quad\quad\quad(q_{1},q_{2},q_{3})=-q\,(5,-3,1)\
,$ (46)

where the parameter $q$ is an integer number that must be different than 1 to
have chiral representations. That of Eq. (46) turns out to be the only
solution of the form $q_{i}=-q\,p_{i}$ with $q$ constant. Its existence can be
traced back to the possibility of embedding $\,{\rm SU}(5)\times\,{\rm U}(1)$
into $\,{\rm SO}(10)$ with the branching rule191919This embedding is valid
only if $(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{\bar{5}},\mathbf{10})$, but
since the mixed anomaly cancellation conditions are not sensitive to the
exchange $\mathbf{5}\leftrightarrow{\mathbf{\bar{5}}}$, the solution exists
also for $(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{5},\mathbf{10})$.

$\mathbf{16}=\mathbf{1}_{5}\oplus\mathbf{\bar{5}}_{-3}\oplus\mathbf{10}_{1}\,.$
(47)

With the correct normalization for $\alpha_{D}$, the two subgroups can be
identified with ${\rm SU}(5)_{\rm GUT}$ and ${\rm U}(1)_{\rm D}$, the
$\psi_{i}$ fields can be organized into a single spinorial representation of
$\,{\rm SO}(10)$, and the same is true for the $\chi_{i}$. Since the
$\mathbf{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 $\,{\rm SO}(10)_{\rm L}\times\,{\rm SO}(10)_{\rm
R}$, and a structure that is highly reminiscent of the theories discussed in
Ref. Vecchi:2021shj . The precise dark fermion content is shown in Tab. 6,
from which it is apparent that the doubling of the $\,{\rm SO}(10)$’s is
necessary to ensure a chiral structure.202020A single $\,{\rm SO}(10)$ would
lead to $q=1$ in Eq. (46), hence to vectorlike ${\rm U}(1)_{\rm D}$ charges.

| ${\rm SU}(N_{\rm DC})$ | $\,{\rm SO}(10)_{\rm L}$ | $\,{\rm SO}(10)_{\rm R}$
---|---|---|---
$\psi$ | $\square$ | $\mathbf{16}$ | $\mathbf{1}$
$\chi$ | $\bar{\square}$ | $\mathbf{1}$ | $\mathbf{\overline{16}}$
Table 6: UV completion of the $n_{f}=3$ model with the choice of ${\rm
U}(1)_{\rm D}$ charges in Eq. (46). SM fermions are assumed to be charged
under $\,{\rm SO}(10)_{L}$.

As in Vecchi:2021shj , all SM generations can be taken to transform as
$\mathbf{16}$ of either $\,{\rm SO}(10)$, and without loss of generality we
choose $\,{\rm SO}(10)_{\rm L}$. With a suitable choice of GUT scalars, the
group $\,{\rm SO}(10)_{\rm L}\times\,{\rm SO}(10)_{\rm R}$ can be broken down
to ${\rm SU}(5)_{\rm GUT}\times{\rm U}(1)_{\rm D}$ at an intermediate scale
$M_{\rm GUT}^{\prime}$, where ${\rm SU}(5)_{\rm GUT}$ is the diagonal $\,{\rm
SU}(5)$ subgroup and ${\rm U}(1)_{\rm D}$ is a linear combination of the two
${\rm U(1)}$’s contained in $\,{\rm SO}(10)_{\rm L}\times\,{\rm SO}(10)_{\rm
R}$. One thus recovers the original model. The Peccei-Quinn symmetry is
explicitly violated by the $\,{\rm SO}(10)_{\rm L}\times\,{\rm SO}(10)_{\rm
R}$ gauging, and emerges as accidental at low energies. Unfortunately, the
appeal of this construction is diminished by the existence of the PQ-violating
operators

$\mathcal{O}_{L}=\psi_{1}\psi_{2}\big{(}\psi_{3}^{*}\big{)}^{2}\,,\quad\quad\quad\mathcal{O}_{R}=\chi_{1}\chi_{2}\big{(}\chi_{3}^{*}\big{)}^{2}$
(48)

that can generate a potential for the axion when combined into a double
insertion. Their product indeed has the right form, as discussed in Section
3.1. Since $\mathcal{O}_{L}$, $\mathcal{O}_{R}$ have dimension 6, their double
insertion is equivalent to a single insertion of a dimension-8 operator.
Therefore, in this model $\Delta_{\cancel{PQ}}\leq 8$ and the axion solution
of the strong CP problem is spoiled for any phenomenologically viable value of
the Peccei-Quinn scale.

In general, one can perform a numerical scan of the charges lying within a
finite interval $|p_{i}|,|q_{i}|\leq n_{\rm{max}}$, and identify models with a
high quality axion. An intrinsic limitation is given by the existence of
higher-dimensional operators that can affect the axion potential for any
values of the charges, as a consequence of anomaly cancellation. Solving Eq.
(38) leads to three dimension-12 operators that are gauge-invariant for any
choice of the ${\rm U}(1)_{\rm D}$ charges and have the correct flavor
structure to generate an axion potential at low energy:

$\mathcal{O}_{1}=\psi_{2}\chi_{2}\big{(}\psi_{3}\chi_{3}\big{)}^{3}\,,\quad\mathcal{O}_{2}=\psi_{1}\chi_{1}\psi_{2}\chi_{2}\big{(}\psi_{3}^{*}\chi_{3}^{*}\big{)}^{2}\,,\quad\mathcal{O}_{3}=\psi_{1}\chi_{1}\big{(}\psi_{2}\chi_{2}\big{)}^{2}\psi_{3}\chi_{3}\,.$
(49)

Therefore, $\Delta_{\cancel{PQ}}\leq 12$ in this case, as shown in Tab. 5. We
have thus performed a search for models with $\Delta_{\cancel{PQ}}=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_{\rm 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.

  * •

One 6-fermion and one 4-fermion operator, with no GUT scalars.

In principle these operators might contain SM fermions in flavor-singlet
combinations. In practice, ${\rm SU}(N_{\rm 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.212121For $N_{DC}=5$, one can construct a gauge
invariant operator with 5 SM fermions and one dark fermion. However, it would
need to be paired with another 6-fermion operator to be of the correct form,
resulting in an effective dimension $d_{\rm{eff}}=14$. However, for
representations $r_{1}=\mathbf{1},r_{2}=\mathbf{\bar{5}},r_{3}=\mathbf{10}$,
it is possible to show that the product of four dark fermions can be always
turned into a singlet of $\,{\rm SU}(5)$ by multiplying it with one or two GUT
scalars.222222Given 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 $\,{\rm 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 $\,{\rm SU}(5)$ is
instrumental to embed the Higgs boson. 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 $\,{\rm 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 ${\rm U}(1)_{\rm D}$
charge.232323In practice, robust models are selected by imposing $\,{\rm
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. Besides selecting robust models, we also performed a
search for high-quality models with a minimal scalar sector, i.e. GUT scalars
transforming as $\mathbf{5}+\mathbf{24}$ of $\,{\rm SU}(5)$.

Table 7 illustrates the result of our numerical study of the parameter
space.242424To perform the classification, we have made use of the Mathematica
package LieART Feger:2012bs ; Feger:2019tvk .

|  | $5\leq N_{\rm{DC}}\leq 11$ |  |  | $N_{\rm{DC}}=4$ |
---|---|---|---|---|---|---
$n_{max}$ | $10$ | $15$ | $20$ | $10$ | $15$ | $20$
AC solutions | $77$ | $189$ | $341$ | $77$ | $189$ | $341$
HQ axions, $(\mathbf{1},\mathbf{\bar{5}},\mathbf{10})$ | $22$ | $68$ | $150$ | $9$ | $31$ | $82$
Robust HQ axions, $(\mathbf{1},\mathbf{\bar{5}},\mathbf{10})$ | $4$ | $16$ | $47$ | $0$ | $5$ | $22$
HQ axions, $(\mathbf{1},\mathbf{5},\mathbf{10})$ | $14$ | $44$ | $99$ | $2$ | $12$ | $33$
Robust HQ axions, $(\mathbf{1},\mathbf{5},\mathbf{10})$ | $3$ | $16$ | $36$ | $0$ | $7$ | $21$
No $d\leq 8$ operators | $0$ | $1$ | $10$ | $0$ | $0$ | $0$
Table 7: Number of $(\mathbf{1},\mathbf{\bar{5}},\mathbf{10})$ and
$(\mathbf{1},\mathbf{5},\mathbf{10})$ GUT models with
$\Delta_{\cancel{PQ}}=12$ and integer charges satisfying $|p_{i}|,|q_{i}|\leq
n_{\rm{max}}$. Models are found through a numerical scan of the parameter
space for different values of $n_{\rm{max}}$ and $4\leq N_{\rm{DC}}\leq 11$.
The counting is made identifying solutions that are equal up to an overall
factor or by exchanging $p_{i}\leftrightarrow q_{i}$. The first row indicates
the total number of independent solutions of the Anomaly Cancellation (AC)
conditions of Eq. (36), while the following rows list how many of these
solutions yield High Quality (HQ) models, either with a minimal GUT sector or
for generic GUT scalar sectors (robust models). The last row shows the number
of models with no flavor-violating operators with dimension $d\leq 8$.

Solutions are identified if they only differ by an integer multiplicative
factor, or if they can be obtained from one another through the exchange
$p_{i}\leftrightarrow q_{i}$. The anomaly cancellation conditions alone admit
a large number of solutions even for small $n_{\rm{max}}$, and it is precisely
thanks to this abundance of solutions that we are able to find examples where
the axion quality is protected up to a sufficiently high level. The
classification does not depend on the number of dark colors, except for the
special value $N_{\rm DC}=4$. In the latter case it is possible to make ${\rm
SU}(N_{\rm DC})$ singlets by pairing 4 (anti)fundamentals, hence more flavor-
violating operators are allowed, implying fewer solutions. With
$n_{\rm{max}}=10$, the following examples are robust in the sense discussed
above for $(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{\bar{5}},\mathbf{10})$ and
$5\leq N_{\rm{DC}}\leq 11$:

$\begin{split}&(p_{1},p_{2},p_{3})=(-5,-6,+2)\quad\quad\quad(q_{1},q_{2},q_{3})=(+10,+3,-1)\\\
&(p_{1},p_{2},p_{3})=(+0,-7,+7)\quad\quad\quad(q_{1},q_{2},q_{3})=(+5,+4,-6)\\\
&(p_{1},p_{2},p_{3})=(+0,+4,-6)\quad\quad\quad(q_{1},q_{2},q_{3})=(+5,-7,+7)\\\
&(p_{1},p_{2},p_{3})=(+2,+3,-5)\quad\quad\quad(q_{1},q_{2},q_{3})=(+3,-6,+6)\,,\\\
\end{split}$ (50)

and the last three are also robust for
$(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{5},\mathbf{10})$. The other 18 high-
quality solutions are reported in Appendix C, together with those for
$N_{\rm{DC}}=4$. We stress that we have not found any charge assignment with
$n_{\rm{max}}=20$ where a single insertion of an operator with dimension
$d<12$ can generate a potential for the axion: it is the presence of multiple
insertions which rules out the excluded models. In this sense, the analysis
carried out in Section 3.1 is crucial to correctly identify the high-quality
models.

Finally, it interesting to notice that $n_{f}=3$ also allows models with the
much stronger property of having no PQ-violating operators (of any form) up to
dimension 12. Assuming that invariance under ${\rm SU}(5)_{\rm GUT}$ can
always be obtained with suitable GUT scalars, there are no such models with
$n_{\rm{max}}=20$, and the one with the smallest charges that we have found is

$(p_{1},p_{2},p_{3})=(+30,+39,+87)\quad\quad\quad(q_{1},q_{2},q_{3})=(+55,-90,-70).$
(51)

##  6 Phenomenology

The rich spectrum of composite states predicted by our models in the Peccei-
Quinn sector is rather elusive and out of the reach of collider experiments.
States with mass of order $\Lambda_{\rm DC}$ are too heavy to be produced,
while SM-neutral NGBs have too weak interactions 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. Marsh:2015xka ;
Hui:2021tkt . 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 Contino:2020god ). Moreover, as we show in Appendix A, the domain wall
number is $N_{\rm{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 Zeldovich:1974uw ;
Sikivie:1982qv .

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 accidentally-stable states with mass $\sim\Lambda_{\rm 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 $\tilde{a}$ and $\tilde{s}$; GUT models with $n_{f}=3$ and
irreducible representations have instead only pseudo NGBs, whose cosmology and
phenomenology are similar to those of $\tilde{a}$.

Let us discuss first the massive pseudo NGBs. The decay rate of Eq. (27) is
sufficiently large so that $\tilde{a}$ decays to photons before Big Bang
Nucleosynthesis (BBN) for all the phenomenologically allowed values of the PQ
scale, $f_{a}\gtrsim 4\cdot 10^{8}\,$GeV. On the other hand, $\tilde{s}$ is
more long-lived, since it decays only through higher-order
$C^{\prime}$-breaking effects. The mass of $\tilde{a}$, $\tilde{s}$ is
parametrically suppressed with respect to $f_{a}$ by a factor of order
$16\pi^{2}\sqrt{(\alpha_{\rm GUT}/4\pi)}\left(f_{a}/M_{\rm GUT}\right)$, see
Eq. (26), and there appear to be two broad scenarios: either i) the reheating
temperature is smaller than $m_{\tilde{a},\tilde{s}}$, or ii) it is in the
range $m_{\tilde{a},\tilde{s}}<T_{\rm RH}<f_{a}$. In the first case, the
states $\tilde{a}$ and $\tilde{s}$ are not populated 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, $\tilde{a}$ and $\tilde{s}$ can be populated at reheating. If $T_{\rm
RH}$ is larger than the temperature $T_{\rm fo}\sim m_{\tilde{a}}(M_{\rm
GUT}/M_{\rm Pl})\sqrt{4\pi}/(\alpha_{s}^{3}\alpha_{\rm
GUT}^{1/2})>m_{\tilde{a}}$ at which processes $gg\leftrightarrow g\tilde{a}$
go out of equilibrium, $\tilde{a}$ 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 of
$\tilde{s}$ with the SM bath are more suppressed than those of $\tilde{a}$,
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, $\tilde{a}$ and $\tilde{s}$ 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 of $\tilde{s}$ can be problematic and in conflict with observations.
It would be particularly interesting to explore the possible role of
$\tilde{a}$ and $\tilde{s}$ in cosmic inflation and their imprints as
primordial non-gaussianities in cosmological correlators, especially because
their mass is parametrically separated from $\Lambda_{\rm DC}$ and can be
close to the scale of inflation. We leave a more detailed analysis of the role
of $\tilde{a}$ and $\tilde{s}$ in cosmology to a future work. In the
following, we focus on the first scenario where $T_{\rm RH}$ is smaller than
the mass of $\tilde{a}$ and $\tilde{s}$, 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
diCortona:2015ldu ; Bonati:2015vqz . Using the notation and the results of
Ref. diCortona:2015ldu , 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:

$\begin{split}&m_{a}=5.70(7)\left(\dfrac{10^{12}\,\rm
GeV}{f_{a}}\right)\,\rm\mu eV,\\\\[5.69046pt]
&g_{a\gamma\gamma}=\dfrac{\alpha_{em}}{2\pi
f_{a}}\left(\dfrac{E}{N}-1.92(4)\right),\\\\[5.69046pt] &c_{p}=-0.47(3),\quad
c_{n}=-0.02(3)\,,\end{split}$ (52)

while the coupling with electrons is generated only through radiative
electromagnetic corrections and is suppressed Srednicki:1985xd ; Chang:1993gm
, $c_{e}\lesssim 3\cdot 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.2, but are suppressed for $\Lambda_{\rm
DC}<M_{GUT}$. The low-energy coupling to photons is predicted in the $\,{\rm
SU}(5)$ models where $E/N=8/3$, while in QCD models it depends on the
hypercharge of dark fermions and is effectively a free parameter.

The cosmological abundance of QCD axions given by the misalignment mechanism
can match the observed dark matter abundance for an initial misalignment
$\theta_{0}=a_{0}/f_{a}$ of order one and values of the axion decay constant
$f_{a}\sim 10^{12}\,\rm GeV$. Larger values of the decay constant require
smaller initial misalignment, up to $\theta_{0}\sim 0.01$ for $f_{a}\sim
10^{16}\,\rm GeV$, or a scale of inflation lower than the QCD scale
Graham:2018jyp ; Guth:2018hsa . Values of $f_{a}$ smaller than $\sim
10^{12}\,\rm GeV$ are allowed and can correspond to the axion accounting for
either a subdominant DM component, if $\theta_{0}\simeq 1$, or all of the dark
matter, if its initial misalignment is tuned to $\theta_{0}\simeq\pi$.

In the scenario we are considering, $f_{a}>H_{I}$, quantum fluctuations of the
axion field during inflation give rise to isocurvature perturbations
Linde:1985yf . If the axion accounts for the whole DM density, the current
limits on isocurvature perturbations give an upper bound on the scale of
inflation $H_{I}\lesssim 10^{9}\,\rm GeV$ for $f_{a}\lesssim 10^{16}\,\rm
GeV$, and also a lower bound on the axion decay constant $f_{a}\geq
10^{10}\,\rm GeV$ diCortona:2015ldu . The upper bound on $H_{I}$ implies that
in single-field slow-roll inflation the tensor to scalar ratio is extremely
small $r\lesssim 2\times 10^{-6}$. The observation of primordial gravitational
waves could therefore falsify our models under the assumption of single-field
slow-roll inflation, unless the axion is a subdominant component of the DM
abundance.

The $n_{f}=2$ models of Section 4 predict an additional light NGB, the singlet
$s$. This is an exact NGB up to non-perturbative gravitational effects, and
has extremely suppressed interactions with the other fields. It behaves
effectively as a free massless real scalar field and, as long as its
cosmological density is small enough, it can be a harmless new light particle
with feeble interactions. Gravitational corrections at the Planck scale can
violate the symmetry ${\rm U}(1)_{s}$ and generate a mass $m_{s}$ such that
$\log\left(m_{s}/M_{\rm Pl}\right)\sim-M_{\rm Pl}/f_{a}$. For values of the
mass in the range $10^{-20}\,{\rm eV}<m_{s}<10^{-10}\,{\rm eV}$, the existence
of this particle may be tested through black hole observations Brito:2015oca ;
Arvanitaki:2016qwi . An estimate along the lines of Ref. Hebecker:2018ofv
seems to suggest $f_{a}\gtrsim 5\cdot 10^{16}\,\rm GeV$ for this to be the
case, with $m_{s}$ exponentially smaller for lower $f_{a}$. Since
$f_{a}<1.2\cdot 10^{16}\,\rm GeV$ is required in our models to have a high-
quality 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 Hui:2016ltb .

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 $\Delta N_{\nu}^{\rm eff}$. An
irreducible component for the QCD axion comes from the thermal production of
axions from the SM bath Salvio:2013iaa ; Arias-Aragon:2020shv ; DEramo:2021psx
. Its value depends on the reheating temperature and on $f_{a}$, but it is
always bounded by $\Delta N_{\nu}^{\rm eff}\leq 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 Hooper:2020evu . 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.

##  7 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\,\text{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) unified dynamics at an
intermediate scale, although the existence of such a dynamics is not required.
The SU(5) models that we describe in Sections 4 and 5 can be downgraded to
purely SM models with straightforward changes, without impairing the solution
to the quality problem. Moreover, purely QCD models based on a different
choice of representations have been described in Section 4.1 and in Appendix
B. In all cases, the QCD axion emerges as the pseudo NGB of an accidental
${\rm U}(1)_{\rm PQ}$ global symmetry, and it can be accompanied by other
light particles, depending on the specific field content. For $n_{f}=2$ models
(both GUT and SM), there is an additional scalar singlet, lighter than the QCD
axion and massless up to UV symmetry breaking effects. Moreover, in the GUT
case (for both $n_{f}=2$ and $n_{f}=3$) axion-like particles (ALPs) or
additional pseudo NGBs appear, with an intermediate mass scale of order
$(g_{\rm GUT}/4\pi)(\Lambda^{2}_{\rm DC}/M_{\rm GUT})$.

Since PQ-violating higher-dimensional operators can spoil the solution of the
strong CP problem, it becomes crucial to understand their effects on the low-
energy physics in detail, in particular on the axion potential. For a
confining ${\rm SU}(N_{\rm DC})$ gauge group and ${\rm SU}(N_{\rm
DC})\times{\rm G}_{\rm SM}$ vectorlike dark fermion representations, we
provided in Section 3.1 simple criteria to discern whether a given set of
higher-dimensional operators can generate a potential for the axion (and
possibly other NGBs). These are derived from a spurion analysis within the
effective field theory as well as using the selection rules associated to the
exact global symmetries. Although the form of operators contributing to the
potential with a single insertion was already stated by Ref. Dobrescu:1996jp ,
the inclusion of flavour breaking weak gaugings and the importance of multiple
insertions had not been discussed in the previous literature. In our examples,
these effects turn out to play a crucial role in excluding otherwise viable
models. Another result of our work is a general analysis of the discrete
symmetries associated to $\,{\rm U}(1)_{\rm PQ}$ in composite axion models
based on a QCD-like confining gauge group, detailed in Appendix A.

The first class of models we considered are those with two dark flavours
($n_{f}=2$) and generic reducible representations under the SM. They are
discussed in Section 4. The simplest version has the matter content shown in
Tab. 1, while Tab. 3 reports all the assignments compatible with
perturbativity of SM gauge couplings and dark color confinement. The ${\rm
U}(1)_{\rm D}$ charges must satisfy the anomaly cancellation conditions but
are otherwise left as free parameters in the analysis. The second class of
models, with three irreducible flavours ($n_{f}=3$), is discussed in Section
5; see Tabs. 4 and 5 in particular. For both $n_{f}=2$ and $n_{f}=3$, we have
found several models where the leading PQ-breaking operators have dimension 9
or higher, and that can thus address the axion quality problem.

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. Randall:1992ut , the PQ-violating dimension
depends on the number $m$ of colors of the weakly-gauged non-abelian group:
$\Delta_{\cancel{PQ}}=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 Dobrescu:1996jp . In the case of
simple GUT extensions of the models of Ref. Randall:1992ut , we find that
multiple insertions with effective dimension 10 can generate an axion
potential, so that $\Delta_{\cancel{PQ}}=\min\\{3m,10\\}$, see Appendix D. The
models of Ref. Gavela:2018paw  can have an accidental PQ symmetry up to
dimension 9, while $\Delta_{\cancel{PQ}}=12$ in the models of Refs.
Redi:2016esr ; Vecchi:2021shj . The constructions of Ref. Fukuda:2017ylt  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 ${\rm
U}(1)_{\rm D}$ charges. By performing a detailed scan of $n_{f}=3$ models with
GUT representations
$(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{\bar{5}},\mathbf{10})$ or QCD
representations $(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{3},\mathbf{6})$, we
have found several charge assignments that lead to $\Delta_{\cancel{PQ}}=12$.
These models are perturbative for values of the axion decay constant as low as
$4\cdot 10^{8}\,$GeV and thus give a robust solution to the axion quality
problem. For $f_{a}\sim 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 $\Delta_{\cancel{PQ}}$ 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 Podo:2022gyj .

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_{\rm
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}\sim 10^{12}\,{\rm GeV}$ and the breaking of
PQ symmetry is induced by $d=12$ local operators generated at $M_{\rm 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
Co:2019jts ; Chang:2019tvx  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.

## Acknowledgements

We are grateful to Luca Vecchi for useful exchanges on his recent work on
composite axions and unified theories Vecchi:2021shj , that appeared while
this work was under completion, and for comments on a preliminary version of
this article. We would like to thank Lam Hui for useful discussions on
cosmological bounds on light scalar dark matter, and Joe Conlon for
conversations on axions in String Theory. This research was supported in part
by the MIUR under contract 2017FMJFMW (PRIN2017). F.R. acknowledges support
from the Dalitz Graduate Scholarship, jointly established by the Oxford
University Department of Physics and Wadham College. The work of A.P. is
supported by the Simons Foundation Award No. 658906.

##  Appendix A Discrete symmetries in composite axion models

The global PQ symmetry is broken both explicitly (by the QCD anomaly) and
spontaneously (by the ${\rm SU}(N_{\rm DC})$ condensate). As first pointed out
by Sikivie in Ref. Sikivie:1982qv , 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 $\,{\rm U}(1)_{\rm
PQ}$ in composite axion models based on a QCD-like confining gauge group.

Let us denote by ${\rm G}_{V}$ the exact global flavour symmetry group
preserved by the vacuum. In our models this is a vectorial group of the form
${\rm G}_{V}=\,{\rm U}(1)_{V,1}\times\dots\times\,{\rm U}(1)_{V,n_{irr}}$,
where each factor $\,{\rm U}(1)_{V,i}$ acts with charge
$Q_{V,i}(\psi_{j})=\delta_{ij}$ and $Q_{V,i}(\chi_{j})=-\delta_{ij}$. As in
the main text, we fix the normalization of the $\,{\rm U}(1)_{\rm PQ}$ charges
of the dark fermions $\psi_{i},\chi_{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\pi)$. In general, the $\,{\rm
U}(1)_{\rm PQ}$ PQ group can have a non-vanishing intersection with the
unbroken flavor group, $\mathbb{Z}_{V}=\,{\rm U}(1)_{\rm PQ}\cap{\rm G}_{V}$,
where $\mathbb{Z}_{V}$ is a discrete group. In other words, the ${\rm
SU}(N_{\rm DC})$ condensates induce the spontaneous breaking

$\,{\rm U}(1)_{\rm PQ}\longrightarrow\mathbb{Z}_{V}\,.$ (53)

The local dark color gauge group ${\rm SU}(N_{\rm 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
Strocchi ). 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
$\,{\rm U}(1)_{\rm PQ}$. Since $\,{\rm U}(1)_{\rm PQ}$ commutes with ${\rm
SU}(N_{\rm DC})$, the transformations of such a global ${\rm SU}(N_{\rm DC})$
that are contained in $\,{\rm U}(1)_{\rm PQ}$ have to be part of the (global)
center group $\mathbb{Z}_{N_{\rm 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: $\mathbb{Z}_{N_{\rm DC}}\subset{\rm G}_{V}$. Therefore, there
is no need to consider separately transformations of the center
$\mathbb{Z}_{N_{\rm DC}}$.

Invariance under the symmetry $\mathbb{Z}_{V}$ restricts the physical domain
of the axion. Denoting by $N_{V}$ the order of $\mathbb{Z}_{V}$, one has

$\dfrac{a}{f_{\rm PQ}}\in\Bigg{[}0,\dfrac{2\pi}{N_{V}}\Bigg{)},$ (54)

and the domain of PQ shift transformations on the axion is correspondingly
restricted.

On the other hand, the anomaly breaks explicitly the $\,{\rm U}(1)_{\rm PQ}$
symmetry down to a discrete group $\mathbb{Z}_{A}$. Since the vectorial group
is anomaly free, it follows that $\mathbb{Z}_{A}$ contains $\mathbb{Z}_{V}$ as
a subgroup, and in particular its order $N_{A}$ is an integer multiple of
$N_{V}$. The quotient group

$\mathbb{Z}_{\rm PQ}=\mathbb{Z}_{A}/\mathbb{Z}_{V},$ (55)

has order $N_{\rm 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:

$\dfrac{a}{f_{\rm PQ}}\rightarrow\dfrac{a}{f_{\rm PQ}}+\dfrac{k}{N_{\rm
PQ}}\dfrac{2\pi}{N_{V}},\qquad\text{with }k=1,\dots,N_{\rm PQ}\,.$ (56)

It is straightforward to show, using the effective chiral lagrangian approach
reviewed in diCortona:2015ldu , that (56) is a symmetry of the non-
perturbative axion potential. It follows that the axion potential has $N_{\rm
PQ}$ distinct minima and the so-called domain wall number $N_{\rm DW}$ equals
$N_{\rm PQ}$.

Let us now determine $N_{A}$ and $N_{V}$ in our models. As in the main text,
we use a normalization where ${\rm Tr}(t_{a}t_{b})=\delta_{ab}/2$ for fields
in the fundamental representation and use a Dirac notation. Under an anomalous
PQ transformation with parameter $\alpha$, the QCD theta angle shifts
by252525The $\,{\rm SU}(5)$ GUT case is obtained by replacing the QCD
generators with the corresponding GUT generators, obtaining the corresponding
$A_{5}$. $\theta_{0}\rightarrow\theta_{0}-2A_{3}\alpha$, with
$\delta_{ab}A_{3}=N_{\rm DC}{\rm Tr}\left[\\{t_{a},t_{b}\\}Q_{\rm PQ}\right]$.
It follows that $N_{A}=2|A_{3}|$.

In order to compute $N_{V}$, we need to determine $\mathbb{Z}_{V}=\,{\rm
U}(1)_{\rm PQ}\cap{\rm G}_{V}$. Notice, first, that fermions $\psi_{i}$ and
$\chi_{i}$, with $i=1,\dots,n_{irr}$, have equal PQ charges $Q_{\rm PQ}^{i}$
and opposite $\,{\rm U}(1)_{V,i}$ charges $\pm 1$. Performing a PQ
transformation with parameter $\alpha$ and independent ${\rm G}_{V}$
transformations with parameters $\beta^{i}$, the $\mathbb{Z}_{V}$ is
determined by solving the system of equations

$\begin{split}&Q_{\rm PQ}^{i}\,\alpha=\beta^{i}+2\pi\kappa^{i}\\\ &Q_{\rm
PQ}^{i}\,\alpha=-\beta^{i}+2\pi\kappa^{\prime\,i},\end{split}$ (57)

with $\kappa^{i}$, $\kappa^{\prime\,i}$ arbitrary integers. Summing the
equations it follows that $\alpha\,Q_{\rm
PQ}^{i}=\pi(\kappa^{i}+\kappa^{\prime\,i})$, and since the charges $Q_{\rm
PQ}^{i}$ are coprime integers, the only non trivial solution in the range
$\alpha\in(0,2\pi)$ is $\alpha=\pi$. We find, therefore, that
$\mathbb{Z}_{V}=\mathbb{Z}_{2}$ and $N_{V}=2$.

We arrive at our final result

$N_{\rm DW}=N_{\rm PQ}=\dfrac{N_{A}}{N_{V}}=|A_{3}|,$ (58)

which is always an integer multiple of $N_{\rm DC}$ and, in particular,
greater than $1$.

##  Appendix B $(\mathbf{1},\mathbf{3},\mathbf{6})$ QCD models

In this Appendix we give a brief overview of $n_{f}=3$ QCD models with
$(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{\bar{3}},\mathbf{6})$ or
$(\mathbf{1},\mathbf{3},\mathbf{6})$, and classify their charge assignments
that can solve the quality problem.

The constraints of Sec. 2.1 imply $3\leq N_{\rm DC}\leq 9$ for these models.
The PQ generator takes the form

$Q_{a}=\text{diag}(+18,+8\,\mathbb{1}_{3},-7\,\mathbb{1}_{6})\,.$ (59)

By counting the residual $\rm{U}(1)$ global symmetries, one can see that there
are no additional massless NGBs in the spectrum, so that the low-energy
phenomenology is entirely characterized in terms of the axion. The anomaly of
$\,{\rm U}(1)_{a}$ with respect to $\rm{SU}(3)_{c}$ is given by

$\langle\partial_{\mu}J_{a}^{\mu}\rangle=-\dfrac{\alpha_{3}}{4\pi}A_{3}\,G\tilde{G}\,,\qquad
A_{3}=-27N_{\rm DC}\,,$ (60)

so that

$f_{a}=\frac{f_{\rm{DC}}}{54N_{\rm{DC}}}\qquad\quad N_{DW}=27N_{DC}\,.$ (61)

The anomaly equations (36) imply the existence of the PQ-violating operator

$\mathcal{O}=\psi_{1}\chi_{1}\big{(}\psi_{2}\chi_{2})^{2}\psi_{3}\chi_{3}\,,$
(62)

that has dimension 12 and the right flavour structure to generate an axion
potential through a single insertion. Solving the strong CP problem for
$f_{a}\sim 10^{12}\,$GeV thus requires the theory to be valid up to the Planck
scale. An intermediate scale of Grand Unification in this case would lead to
too large corrections to the axion potential, as estimated in Section 3, and
is therefore not compatible with the axion solution.

High quality models with $\Delta_{\cancel{PQ}}=12$ can be characterized as
done in Section 5.2. The only difference with respect to the GUT case is the
absence of $\rm{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.

|  | $N_{\rm{DC}}=3,5,\dots,9$ |  |  | $N_{\rm{DC}}=4$ |
---|---|---|---|---|---|---
$n_{\rm max}$ | $10$ | $15$ | $20$ | $10$ | $15$ | $20$
AC solutions | $16$ | $40$ | $96$ | $16$ | $40$ | $96$
HQ axions, $(\mathbf{1},\mathbf{\bar{3}},\mathbf{6})$ | $3$ | $14$ | $50$ | $3$ | $11$ | $43$
HQ axions, $(\mathbf{1},\mathbf{3},\mathbf{6})$ | $4$ | $17$ | $50$ | $2$ | $15$ | $43$
No $d\leq 8$ operators | $0$ | $1$ | $8$ | $0$ | $0$ | $1$
Table 8: Number of $(\mathbf{1},\mathbf{\bar{3}},\mathbf{6})$ and
$(\mathbf{1},\mathbf{3},\mathbf{6})$ QCD models with $\Delta_{\cancel{PQ}}=12$
and integer charges satisfying $|p_{i}|,|q_{i}|\leq n_{\rm{max}}$. Models are
found through a numerical scan of the parameter space for different values of
$n_{\rm{max}}$ and $3\leq N_{\rm{DC}}\leq 9$. The counting is made identifying
solutions that are equal up to an overall factor or by exchanging
$p_{i}\leftrightarrow q_{i}$. The first row indicates the total number of
independent solutions of the Anomaly Cancellation (AC) conditions of Eq. (36),
while the following rows list how many of these solutions yield High Quality
(HQ) models. The last row shows the number of models with no operators of
dimension $d\leq 8$ with a non-trivial flavour structure.

The three charge assignments with $n_{\rm max}=10$ that solve the quality
problem with $(r_{1},r_{2},r_{3})=(\mathbf{1},\mathbf{\bar{3}},\mathbf{6})$
and $3\leq N_{\rm DC}\leq 9$ are:

$\begin{aligned} &(-1,-7,+1;+10,+2,+0)\\\ &(+3,-7,+7;+6,+2,-6).\\\
\end{aligned}\qquad\qquad\qquad\begin{aligned} &(-1,+2,+1;+10,-7,+0)\\\ &\\\
\end{aligned}$ (63)

Here and in the following, the ${\rm U}(1)_{\rm 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})=(\mathbf{1},\mathbf{3},\mathbf{6})$ and
$N_{\rm{DC}}=3,5,\dots,9$, the solutions are

$(-1,-7,+1;+10,+2,+0)\,\,\,\quad\quad\quad\quad\quad(-1,+2,+1;+10,-7,+0)$ (64)

and

$(+3,-7,+7;+6,+2,-6)\,\,\,\quad\quad\quad\quad\quad(+4,-4,+0;+5,-1,+1).$ (65)