Window for preferred axion models

We discuss phenomenological criteria for defining"axion windows", namely regions in the parameter space of the axion-photon coupling where realistic models live. Currently, the boundaries of this region depend on somewhat arbitrary criteria, and it would be highly desirable to specify them in terms of precise phenomenological requirements. We first focus on hadronic axion models within post-inflationary scenarios, in which the initial abundance of the new vectorlike quarks $Q$ is thermal. We classify their representations $R_Q$ by requiring that $i)$ the $Q$ are sufficiently short lived to avoid issues with long-lived strongly interacting relics, $ii)$ the theory remains weakly coupled up to the Planck scale. The more general case of multiple $R_Q$ is also studied, and the absolute upper and lower bounds on the axion-photon coupling as a function of the axion mass is identified. Pre-inflationary scenarios in which the axion decay constant remains bounded as $f_a\leq 5\cdot 10^{11}\,$GeV allow for axion-photon couplings only about 20% larger. Realistic Dine-Fischler-Srednicki-Zhitnitsky type of axion models also remain encompassed within the hadronic axion window. Some mechanisms that can allow to enhance the axion-photon coupling to values sizeably above the preferred window are discussed.


I. INTRODUCTION
It is well known that the standard model (SM) of particle physics does not explain some well established experimental facts like dark matter (DM), neutrino masses, and the cosmological baryon asymmetry, and it also contains fundamental parameters with highly unnatural values, like the coefficient µ 2 ∼ O((100 GeV) 2 ) of the quadratic term in the Higgs potential, the Yukawa couplings of the first family fermions h e,u,d ∼ 10 −6 − 10 −5 and the strong CP violating angle |θ| < 10 −10 . This last quantity is somewhat special: its value is stable with respect to higher order corrections [1] (unlike µ 2 ) and (unlike h e,u,d [2]) it evades explanations based on environmental selection [3]. Thus, seeking explanations for the smallness of θ independently of other "small values" problems is theoretically motivated. While most of the problems of the SM can be addressed with a large variety of mechanisms, basically only three types of solutions to the strong CP problem have been put forth so far. The simplest possibility, a massless up-quark, is now ruled out (m u = 0 by 20 standard deviations [4,5]). The so-called Nelson-Barr (NB) type of models [6,7] either require a high degree of fine tuning, often comparable to setting |θ| < ∼ 10 −10 by hand, or additional and rather elaborated theoretical structures to keep θ sufficiently small at all orders [8,9]. The Peccei-Quinn (PQ) solution [10,11] arguably stands on better theoretical grounds, and from the experimental point of view it * luca.di-luzio@durham.ac.uk † mescia@ub.edu ‡ enrico.nardi@lnf.infn.it also has the advantage of predicting an unmistakable signature: the existence of a new light scalar particle, universally known as the axion [12,13]. Therefore, the issue if the PQ solution is the correct one, could be set experimentally by detecting the axion. In contrast, no similar unambiguous signature exists for NB models.
A crucial challenge for axion models is to explain through which mechanism the global U (1) P Q symmetry, on which the solution relies (and that presumably arises as an accident), remains protected from explicit breaking to the required level of accuracy [14][15][16], and it seems fair to state that only constructions that embed such an explanation can be considered theoretically satisfactory. A wide variety of proposals to generate a high quality U (1) P Q have been put forth based, for example, on discrete gauge symmetries [17][18][19][20], supersymmetry [15,21,22], compositeness [23][24][25][26], flavour symmetries [27] or new continuous gauge symmetries [28,29]. Regardless of the details of the different theoretical constructions, many properties of the axion remain remarkably independent from specific model realizations. It is then very important, in order to focus axion searches, to identify as well as possible the region in parameter space where realistic axion models live. The vast majority of axion search techniques are sensitive to the axion-photon coupling g aγγ which is inversely proportional to the axion decay constant f a . Since the axion mass m a has the same dependence, the experimental exclusion limits, as well as the theoretical predictions for specific models, can be conveniently presented in the m a -g aγγ plane (see Fig. 3). The commonly adopted "axion band" cor-arXiv:1705.05370v2 [hep-ph] 26 Sep 2017 responds roughly to with a somewhat arbitrary width chosen to include representative models as e.g. those of Refs. [30][31][32]. Recently, in Ref. [33] we have put forth a definition of a phenomenologically preferred axion window as the region encompassing hadronic axion models which i) do not contain cosmologically dangerous strongly interacting relics; ii) do not induce Landau poles (LP) below a scale Λ LP of the order of the Planck scale. In this paper we will first present a more detailed analysis of the phenomenological constraints on hadronic axion models (to which we will often refer also as Kim-Shifman-Vainshtein-Zakharov (KSVZ) [34,35] type of axion models) on which the study of Ref. [33] was based.
Since the first condition i) is relevant only when the heavy quarks Q have an initial thermal abundance, the validity of the analysis in Ref. [33] is restricted to the case when T reheating m Q . The Q acquire their mass via a Yukawa coupling with the complex axion field so that, for Yukawa couplings not exceeding unity, this translates into T reheating f a (where f a is the axion decay constant) a condition that can be only realized when the PQ symmetry is broken after inflation, and will be referred as post-inflationary scenario. However, astrophysical considerations imply a lower bound f a 10 9 GeV, while the only firm limit on the scale of inflation is provided by big bang nucleosynthesis (BBN) to merely lie above a few MeV. Since this leaves ample space for axion models to be realized in pre-inflationary scenarios, in which the initial Q abundance is completely negligible, it would be interesting to generalize the analysis of [33] by dropping condition i). Such a generalization will be carried out in section VI, subject to the only condition that f a ≤ 5 × 10 11 GeV, which restricts the class of models to those which do not require any ad hoc tuning (or anthropic selection arguments) to justify particularly small initial values of θ. As we will show, in pre-inflationary scenarios the LP condition ii) alone is sufficiently strong that the limits found in [33] get relaxed at most by ≈ 20%. In section VII we extend the analysis to include also the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) axion [36,37], that was not considered in [33], together with several of its variants, to which we will collectively refer as DFSZ-type of models. We will argue that the same window that encompasses preferred hadronic axion models, also includes the majority of realistic DFSZ scenarios.
The layout of the paper is the following: in section II we introduce hadronic axion models with some focus on the issue of the stability of the new heavy quarks Q. Section III is devoted to the cosmological consequences of stable or long-lived Q's: we estimate the present abundances of strongly interacting relics arguing that absolutely stable Q's are likely excluded, and we review the constraints on the lifetimes of metastable Q's. In section IV we put forth a definition of preferred hadronic (KSVZ) axion models on the basis of our two well-defined discriminating criteria. In section V we identify the window in parameter space where preferred axion models live, and we discuss the corresponding maximum and minimum values allowed for the axionphoton coupling. In section VI we address KSVZ models in pre-inflationary scenarios showing that the previous results get only mildly relaxed as long as the requirement f a ≤ 5 × 10 11 GeV on the axion decay constant is maintained. In section VII we address DFSZ-type of axion models showing that the same window also includes realistic models of this type. Finally, in section IX we review the main results and draw the conclusions.

II. HADRONIC AXIONS
The basic ingredient of any axion model is a global U (1) P Q symmetry. The associated Noether current must have a color anomaly and, although not required for solving the strong CP problem, in general it also has an electromagnetic anomaly: where G, F are the color and electromagnetic field strength tensors,G,F their duals (e.g. F ·F ≡ 1 2 µνρσ F µν F ρσ , etc.), and N and E are the color and electromagnetic anomaly coefficients. In a generic axion model of KSVZ type [34,35] the anomaly is induced by pairs of heavy fermions Q L , Q R which must transform non-trivially under SU (3) C and chirally under U (1) P Q . Their mass arises from a Yukawa interaction with a SM singlet scalar field Φ which develops a PQ breaking vacuum expectation value. Therefore their PQ charges X L,R ≡ X (Q L,R ), normalized to X (Φ) = 1, must satisfy In KSVZ models the SM fermions do not contribute to the color or electromagnetic anomalies so that their PQ charges can be set to zero. We denote the (vectorlike) representations of the SM gauge group , so that the anomaly coefficients read: Here Q denotes the sum over all irreducible SU (3) C × U (1) em representations (we allow for the simultaneous presence of more R Q ). The color index is defined by Tr T a Q T b Q = T (C Q )δ ab with T Q the generators in C Q (in particular, T (3) = 1/2, T (6) = 5/2, T (8) = 3, T (15) = 10) and Q Q denotes the U (1) em charge. Different choices for R Q imply different phenomenological consequences, and we will use this fact to identify phenomenologically preferred models. Let us parametrize the scalar field Φ as ρ(x) acquires a mass m ρ ∼ V a while a(x) is the axion field which would remain massless in the absence of explicit U (1) P Q breaking. In invisible axion models, in order to sufficiently suppress the axion couplings that scale as 1/f a ≡ 2N/V a , it is assumed that V a ( √ 2G F ) −1/2 = 247 GeV. More quantitatively, astrophysical constraints hint to a lower limit f a > ∼ 4 · 10 8 GeV [38,39].
The renormalizable Lagrangian for a generic hadronic axion model can be written as: where L SM is the SM Lagrangian, with Q = Q L + Q R and, from the last term, m Q = y Q V a / √ 2. V HΦ contains the new scalar couplings: Finally, L Qq contains possible renormalizable terms coupling Q L,R to the SM quarks q = q L , d R , u R , which can allow Q decays [20]. Note however, that SM gauge invariance allows L Qq = 0 only for few specific R Q and, for example, the original KSVZ assignment R Q = (3, 1, 0) [34,35] implies L Qq = 0 (and it would in fact forbid Q-decays to all orders).

A. Q stability and PQ quality
The issue whether the Q's are exactly stable, metastable, or decay with safely short lifetimes, is of central importance for KSVZ models in postinflationary scenarios, and we will now discuss it in more detail. The gauge invariant kinetic term in symmetry corresponding to independent rephasing of the Q L,R and Φ fields. The PQ Yukawa term (y Q = 0) breaks U (1) 3 → U (1) P Q ×U (1) Q where, in analogy to ordinary baryon number U (1) B for the SM quarks, U (1) Q is the Q-baryon number of the new quarks [34] under which Q L,R → e iβ Q L,R and Φ → Φ. Moreover, U (1) Q being vectorlike is not broken by anomalies. If it were an exact symmetry, U (1) Q would ensure absolute Q stability, a possibility which is preferable to avoid. In the few cases in which L Qq = 0 is allowed by G SM gauge invariance, U (1) Q × U (1) B is further broken to U (1) B , that is a generalized baryon number extended to the Q's, which can then decay into SM quarks with unsuppressed rates. However, whether L Qq is allowed at the renormalizable level, does not depend solely on R Q , but also on the specific PQ charges. For example, independently of R Q , if U (1) P Q were an exact symmetry, the common assignment X L = −X R = 1 2 would forbid PQ invariant decay operators of any dimension. More realistically, both U (1) P Q and U (1) Q are expected to be broken at least by Planckscale effects, inducing PQ violating contributions to the axion potential V d>4 Φ as well as an effective Lagrangian L d>4 Qq . In particular, in order to preserve |θ| < 10 −10 , operators in V d>4 Φ must be of dimension d ≥ 11 [14][15][16], and if L d>4 Qq had to respect U (1) Q to a similar level of accuracy, then the Q's would behave as effectively stable. However, a scenario in which the global U (1) Q symmetry arises as an accident because of specific assignments for the charges of another global symmetry U (1) P Q seems theoretically untenable. It would be instead desirable to enforce on the basis of first principles a situation in which (i) U (1) P Q arises accidentally and is of the required high quality, (ii) U (1) Q is either broken at the renormalizable level, or it can be of a sufficient bad quality to allow for sufficiently fast Q decays. Here we will not commit ourselves to any specific mechanism to realize such a scenario, and we will simply adopt a technical solution to this issue: a discrete (gauge) symmetry Z N under which Φ → ω Φ (with ω ≡ e i2π/N ) which can automatically ensure that the minimum dimension of the PQ breaking operators in V d>4 Φ is N, so that the first condition is satisfied if N ≥ 11. In Ref. [33] it was shown that at the same time suitable transformations for Q L,R under Z N can be found that allow the Q's to decay via operators of much lower dimension d ≤ 5. Although, admittedly, such a solution seems just as an ad hoc construction, it suffices to ensure that it is consistent to assume that a high quality U (1) P Q can live together with a U (1) Q of sufficiently bad quality.

III. HEAVY QUARKS COSMOLOGY
We start by assuming a post-inflationary scenario (U (1) P Q broken after inflation). In this case, requiring that the axion energy density from vacuum realignment does not exceed Ω DM implies f a 5·10 11 GeV [40][41][42]. We further assume m Q < T reheating so that via gauge interactions the Q's will attain a thermal distribution, providing the initial conditions for their cosmological history, which will then depend only on their mass m Q and representation R Q .
For some R Q the heavy quark can only hadronize into fractionally charged hadrons, and in this case, as detailed in Appendix A, decays into SM particles are forbidden. These Q-hadrons must then exist today as stable relics. Searches for fractionally charged particles limit their abundance with respect to ordinary nucleons to n Q /n b < ∼ 10 −20 [43]. This is orders of magnitude below any reasonable estimate of their relic abundance and of their concentrations in bulk matter. This restricts the possible R Q 's to the much smaller subset which allows for integrally charged or neutral color singlet Q-hadrons, in which case the limits on cosmologically stable heavy relics are less tight. However, for each R Q belonging to this subset it is always possible to construct gauge invariant operators that can allow the Q to decay into SM particles (see Appendix A). Let us start by discussing the case of lifetimes τ Q shorter than the age of the Universe, so that no heavy relics are left around during the present era. Cosmological observations severely constrains the allowed values for τ Q . For τ Q ∼ (10 −2 − 10 12 ) s the decays of superheavy quarks with m Q 1 TeV would affect BBN [44][45][46][47]. Early energy release from heavy particles decays with lifetimes ∼ (10 6 − 10 12 ) s is strongly constrained also by limits on CMB spectral distortions [48][49][50], while Q's decaying around the recombination era (t rec ∼ 10 13 s) are tightly constrained by measurements of CMB anisotropies. Decays after recombination would give rise to free-streaming photons visible in the diffuse gamma ray background [51], and tight constraints from Fermi LAT [52] allow to exclude τ Q ∼ (10 13 − 10 26 ) s. Note that these last constraints are also able to exclude lifetimes that are several order of magnitude larger than the age of the Universe, t U ∼ 4 · 10 17 s.

A. Abundance of strongly interacting relics
Cosmologically stable Q's are severely constrained by the requirement that their present energy density does not exceed that of the DM Ω Q ≤ Ω DM ∼ 0.12 h −2 . Obtaining reliable estimates of Ω Q is a non-trivial task though. Some controversy in the results exists, mainly related to possible large enhancements of the annihilation rate with respect to free Q's annihilation, which can occur after the Q's get confined into hadrons. We now review the state of the art, and we formulate a motivated guess about the most reasonable range of values for Ω Q .
At temperatures above the QCD phase transition T C ∼ 180 MeV the Q's annihilate as free quarks. One generally assumes a symmetric scenario n Q = n Q since any asymmetry would eventually quench QQ annihilation resulting in stronger bounds. Perturbative computations in this regime are reliable, and give: where n f is the number of quark flavors into which Q can annihilate, and (c f , c g ) = ( 2 9 , 220 27 ) for triplets [53], ( 3 2 , 27 4 ) for octets [54], etc. Freezeout of free Q annihilation occurs around T f o ∼ m Q /25 and, for m Q > few TeV, at T f o there are g * = 106.75 effective degrees of freedom in thermal equilibrium. Together with Eq. (10) this gives: where the second equality holds for color triplets and for reference values of the relevant parameters. The upper lines in Fig. 1 give Ω Q h 2 Free as a function of m Q for SU (3) C triplets (dotted) and octets (dashed), including the running of α s (µ = m Q ) computed at two-loops. We see that only in a narrow interval at low m Q the limit Ω Q ≤ Ω DM ∼ 0.12 h −2 is respected, and an improvement of the lower limit on m Q by a factor of a few would exclude also this region.
For T < T C the Q's get confined in color singlet hadrons, that can be pictured as a heavy parton surrounded by a QCD cloud ("brown muck") of light degrees of freedom. As the temperature decreases below T C , the presence of a baryon asymmetry for the SM quarks implies that the brown muck is preferentially constituted by light quarks q (and eventually gluons g) rather than by antiquarks q. For example, for a color triplet, the heavy meson Qq will readily scatter with ordinary nucleons which are relatively much more abundant, giving rise to a heavy baryon: Qq + qqq → Qqq + qq. For Q's belonging to different SU (3) C representations different types of color singlet baryons and mesons, including exotic and hybrid states, will eventually form, for example: Qq (Q ∈ 3), Qqq (Q ∈ 3, 6), Qqqq (Q ∈ 8, 10), Qgg (Q ∈ 10, 10), etc. However, the most important feature concurring to determine possible large enhancements of the annihilation cross section after confinement is largely independent of many fine details, and is essentially related to the fact that the QCD cloud of light quarks and gluons surrounding the heavy partons results in composite states of typical hadronic size R h ∼ fm. Then, because of finite size effects of the hosting composite state, Q annihilation can restart below T C , and the relic abundance of Q-hadrons can get further depleted until a new freezeout temperature is reached. Clearly, if this picture is correct, annihilation in the perturbative regime will be to a large extent irrelevant, and the relic density of Q-hadrons will be essentially fixed by non perturbative processes occurring at T < T C . Presently, agreement on quantitative estimates of the annihilation rate for massive colored particles confined into hadrons has not been reached, and published estimates for the relic density span over several orders of magnitude. This issue is of central importance for the present study so that, after reviewing the relevant literature, we will attempt to pin down some general conclusion.
In Ref. [55] the relic density of confined heavy stable color sextet quarks was estimated by assuming an annihilation cross section of typical hadronic size σ ann ∼ (m 2 π v) −1 ∼ 30v −1 mb. This resulted in n Q /n b ∼ 10 −11 , where n b is the present abundance of baryons. In [56] it was remarked that the quoted value of σ ann was likely overestimating the annihilation cross section by a few orders of magnitude. This is because it corresponds to a typical inclusive hadronic scattering cross section, and it was argued that the relevant exclusive annihilation channel (not containing the heavy Q quarks in the final state) would not exceed a fraction f < 1 of the geometrical cross section σ r ∼ πR 2 h . In this case, even assuming f ∼ 1, a much larger abundance n Q /n b ∼ 10 −9 (m Q /10 TeV) −1 would result. The relic density of heavy stable gluinos (that is Q's in the adjoint of SU (3) C ) was studied in [54] considering several different possibilities for σ ann ranging from perturbative, perturbative dressed with Sommerfeld enhancements, and various non-perturbative possibilities. For the reference value m Q = 10 TeV they found results ranging from overclosure Ω Q h 2 ∼ 1 down to Ω Q h 2 ∼ 5 · 10 −10 which corresponds to Cosmologically long-lived gluinos were reconsidered in [57]. The authors correctly identify the relevant cross section as the one characterizing, after the QCD phase transition, the annihilation of Qhadrons. Similarly to Ref. [56] they argue that the expected cross section is not of hadronic size, but it should rather be characterized by the size of the heavy parton localized inside the hadron core, i.e. σ ∝ 1/m 2 Q , a conclusion apparently supported by the partial wave unitarity limit for the inelastic cross section [58] . However, the very large mass m Q TeV relative to the typical energy E Q < ∼ T C implies that the corresponding momentum is large, and many partial waves are involved in the collision. If all partial waves up to [56] would be recovered. Clearly, σ inel is not by itself the QQ annihilation cross section, but includes all inelastic processes as e.g. the formation of bound states out of the collision of two Q-hadrons, apparently supporting the conclusion that f < 1. However, if the bound states can efficiently radiate off large amounts of angular momentum, the Q and Q wavefunctions will eventually be brought to overlap so that also bound state formation would contribute to annihilation. Collapse to states of low angular momentum must however occur in a relatively short time, so that annihilations will occur well before the BBN era. Ref. [57] estimated the rate for angular momentum radiation via π emission, and found it to be very small, concluding that bounded Q's would remain incapable of annihilating on a sufficiently short time scale. The most conservative limits quoted in that paper are then obtained under the assumption that the annihilation cross section saturates the s and p wave unitarity limits, which for m Q = 10 TeV yields n Q /n b > ∼ 10 −4 , close to saturation of the cosmological limit Ω Q < Ω DM .
Annihilation via bound state formation was reconsidered in [59] and, as regards the radiation time to collapse down to low angular momentum states, opposite conclusions with respect to [57] were reached: they find that for charged Q-hadrons, photon radiation can collapse the bound state with a time scale τ rad < ∼ 1 s for all masses m Q < ∼ 10 11 GeV (for neutral partons and neutral host Q-hadrons however, this reduces to m Q < ∼ 2.5 TeV). Their conclusion is that for charged states the relic density of Q's gets sizeably depleted by the second stage of annihilation after hadronization. For the Q contribution to the present energy density the results of [59] imply: The mechanism of annihilation via bound state formation was put under closer scrutiny in [60], where previously neglected effects of the large number of thermal bath pions (n π n Q ) on the bound states were considered. The two relevant effects are breakup of the bound states due to collisions with π's with energy larger than the typical binding energy E B ∼ few 100 MeV, and de-excitation processes through which the colliding π carries away two units of angular momentum. These two processes work in the opposite directions of delaying and speeding up QQ annihilation, and it was estimated that eventually they would roughly equilibrate each other, yielding a result not far from the estimates in [59].
A more quantitative study of this mechanism was carried out in Ref. [61]. The conclusion was that Eq. (12) represents a conservative lower limit on Ω Q , but that much larger values are also possible. In fact [59,60] did not consider the possible formation of (QQ...) bound states which, opposite to QQ, would hinder annihilation rather than catalyze it. This possibility was discussed in Ref. [61] but was not included in their quantitative analysis. However, doubly and triply heavy baryons, like Ω ccq , Ω bbb (see [62,63] for compilations of recent results) are a firm prediction of the quark model, also supported by the recent discovery of the doubly heavy hadron Ξ + cc by the LHCb collaboration [64]. Clearly, the size of bound states solely composed by Q's or by Q's would be much smaller than hadronic, quenching all enhancements of the annihilation, and if a relevant fraction of Q's ends up in multi-Q bound states, the final relic density would be better approximated by the free quark result Eq. (11) rather than by Eq. (12).
By evaluating Eq. (12) for reference values of the relevant parameters, we obtain the continuous line in Fig. 1 which, according to the discussion above, should be understood as a rather conservative lower limit on Ω Q h 2 . However, even assuming that the relic abundance approaches this lower limit, the relative concentration of Q-hadrons n Q /n b ∼ 10 −8 (m Q /TeV) 1/2 would still be quite large. If the Q's accumulate with similar concentrations within the galactic disk, existing limits from searches of anomalously heavy isotopes in terrestrial, lunar, and meteoritic materials [65] would be able to exclude their existence for most of the range of masses allowed by the Ω Q < Ω DM constraint. Many other arguments have been put forth disfavoring the possibility of heavy stable Q's: their capture in neutron stars would form black holes on a time scale of a few years [66] and, more generically, they could endanger stellar stability [67], while their annihilation in the Earth interior would result in an anomalously large heat flow [68]. In conclusion, unless an extremely efficient mechanism exists that keeps Q-matter completely separated from ordinary matter, the possibility of stable Q-hadrons that were once in thermal equilibrium is ruled out.

B. Q Lifetimes
We have seen in the previous sections that cosmologically stable heavy Q's with m Q < T reheating are strongly disfavored, and that in case they are unstable, only lifetimes τ Q < ∼ 10 −2 s are safe with respect to cosmological issues. For the post-inflationary case, we will then consider as phenomenologically preferred only scenarios in which this condition can be satisfied. The order of magnitude of τ Q crucially depends on the dimension d of the operators responsible for the decays. Below we derive quantitative estimates for τ Q as a function of m Q and d, and we argue that only for d = 4 and d = 5 the constraint τ Q < ∼ 10 −2 s can be satisfied in a natural way.
Depending on their gauge quantum numbers, the Q's can couple directly to SM quarks via renormalizable operators. All the representations that allow for this possibility are listed in Table I. We have basically two different cases: i) d = 3 super-renormalizable operators like, for example, µ Qq Q L d R as in the first row in Table I, or d = 4 operators involving Φ which generate effective d = 3 mixing operators after PQ symmetry breaking, like for example λ Qq V a Q L d R as in the third row in Table I; ii) genuine d = 4 operators, like for example λ QqH q L Q R H as in the second row in Table I. For m Q > ∼ 10 TeV, unless the relevant couplings have exceedingly small values (µ Qq , λ Qq V a 1 keV, λ QqH 10 −12 ), τ Q < 10 −2 s is always ensured. For some R Q 's, even if renormalizable decay operators are forbidden by gauge invariance, effective operators of dimension d > 4 can still be allowed. We assume conservatively that higher dimensional operators are suppressed by powers of the Planck mass m P = 1.2 × 10 19 GeV, and we write them as: For constant matrix elements and massless final states, the phase space factor for Q decays into n f final states can be integrated analytically (see e.g. [69]), yielding the decay rate Q-decay operators of d = 5, 6, 7 involve at least n f = 2, 3, 4 particles in the final state, thus we obtain: where we have normalized m Q to its (presumably) largest value V a < ∼ 5 · 10 11 GeV compatible, in postinflationary PQ breaking scenarios, with an axion energy density not exceeding Ω DM [40][41][42]. 1 Our results for τ d=5,6,7 are plotted in Fig. 2. We see that for d = 5, decays can occur with lifetimes shorter than 10 −2 s as long as m Q > ∼ 800 TeV. For d = 6, even when m Q ∼ V a decays occur dangerously close to the BBN era. Finally, decays via d = 7 operators are always excluded. We can then conclude that only d ≤ 4 and d = 5 operators naturally yield sufficiently fast decays, so that only the R Q listed in tables Table I and Table II are safe from cosmological issues.

IV. SELECTION CRITERIA
In this section we proceed to select KSVZ-type (or hadronic) axion models which satisfy the follow- ing two criteria: i) cosmologically safe Q lifetimes, and ii) the absence of LP in the SM gauge couplings at sub-Planckian scales. We will define as phenomenologically preferred those post-inflationary models which satisfy these two criteria. We will also briefly comment on two other possible criteria, namely the absence of the domain wall (DW) problem, and R Q -assisted improved gauge coupling unification. However, we will eventually conclude that these two additional conditions do not match a sufficient level of generality to represent reliable selection criteria, and should be better considered just as desirable features of specific models. After discussing the case of hadronic (KSVZ) axions in postinflationary scenarios, in section VI we will generalize the analysis to include pre-inflationary models.
A. First criterium: Q lifetimes As a first discriminating criterium we assume that: Models that allow for sufficiently short lifetimes (τ Q < ∼ 10 −2 s) are phenomenologically preferred with respect to models containing long-lived or cosmologically stable Q's.
According to the analysis in section III B, only d ≤ 4 and d = 5 operators are safe from cosmological issues. The quantum number assignments that allow for d ≤ 4 decay operators (L Qq = 0) are collected in Table I. These operators induce 2body decays either directly or via Q-q mass mixing, allowing in both cases for fast Q decays. Out of the seven possibilities listed in Table I, the ones in the third and fifth row were already identified in Ref. [20] (see also [70][71][72]). They coincide in fact with models KSVZ-II and KSVZ-III of [20] modulo a redefinition of the PQ charges by a shift proportional to baryon number X → X + B , respectively with B = − 3 2 and B = + 1 2 . As long as only B conserving operators are considered, this gives no difference in the physics.
In Table II we give the list of R Q for which d = 5 operators involving a single Q-insertion are allowed. In this case Q decays occur with sufficiently short lifetimes only for m Q > ∼ 800 TeV. Large representations can often induce LP in the hypercharge, weak, or strong gauge couplings g 1 , g 2 , g 3 at some uncomfortably low-energy scale Λ LP < m P . At energy scales approaching m P , gravitational corrections to the running of the gauge couplings can become relevant, and explicit computations show that they go in the direction of delaying the emergence of LP [73]. Therefore, to be conservative, we choose a value of Λ LP for which gravitational corrections are presumably negligible. As a second discriminating criterium we then assume that: Models in which LP in the SM gauge couplings appear only above Λ LP ∼ 10 18 GeV are phenomenologically preferred. We evaluate the scale at which the LP arise by setting conservatively the threshold for the R Q representations at m Q = 5 · 10 11 GeV. In first approximation the scaling of the LP is linear with m Q , though sizable deviations from linearity are expected in case of several decades of running. We employ two-loop beta functions to avoid possible accidental cancellations which can arise for some representations in the one-loop coefficients [69,74]. Among the R Q which allow for d = 5 decay operators, the five highlighted in light red in Table II lead to LP below 10 18 GeV and we consider them as theoretically disfavored.

C. Domain walls
The axion field a, being an angular variable, takes values in the interval [0, 2πV a ). The QCD induced axion potential is periodic in a with period ∆a = 2πV a /(2N ), and is thus characterized by an exact Z 2N discrete symmetry. Once at T ∼ Λ QCD the explicit U (1) P Q breaking from non-perturbative QCD effects starts lifting the axion potential, within the domain of definition of a, N DW = 2N degenerate vacua appear, and if the initial value of a is different in different patches of the Universe out of causal contact (as is the case in post-inflationary scenarios), in each of these patches the axion field will flow towards a different minimum, breaking spontaneously Z N DW with N DW different vacuum values of a. Then, at T < ∼ Λ QCD DWs will form at the boundaries between regions of different vacua. This leads to the so-called cosmological DW problem [75] which consists in the fact that the energy density of the DWs will largely overshoot the critical density of the Universe.
In pre-inflationary scenarios the problem is avoided at once since the whole observable Universe corresponds to an initial patch characterized by a unique value of a, and which gets exponentially inflated to super-horizon scales. For post-inflationary scenarios some solutions also exist [76,77]. The DW problem is avoided if N DW = 1 [78,79] so we might want to consider this specific value as an additional desirable feature for axion models. In the last column in Table III we list the number of DW for each R Q , and we see that only the first two cases have N DW = 1. However, if multiple R Q are considered, one can conceive new solutions. For instance, given the color indices T (8) = 3 and T (6) = 5/2, N DW = 1 models can be constructed by combining one (8, 1, Y ) with a (6, 1, Y) with opposite PQ charge difference. With three R Q it is also possible to have N DW = 1 in a trivial way: by canceling the DW contribution of two R Q and leaving a third one with N DW = 1.
Models with N DW > 1 can also remain viable in post-inflationary scenarios, but additional assumptions are needed. The DW problem can be disposed of in a simple way by introducing an explicit small soft breaking of the PQ symmetry [75]. On one hand, the size of this breaking should be large enough so that a single (true) vacuum can take over before the DWs start dominating the energy density. On the other hand, it should be sufficiently small to ensure that the PQ solution does not get spoiled, as it would occur if θ gets shifted away from zero by more than ∼ 10 −10 . Since there is a sizable region in parameter space where these two conditions can be simultaneously matched [39], we can conclude that the DW problem can be solved also in N DW > 1 models and, accordingly, we prefer not to consider N DW = 1 as sufficiently motivated condition for selecting post-inflationary preferred axion models.

D. Q-assisted unification
Fixing the threshold for the new quark representations R Q at some suitable intermediate scale could improve on SM gauge coupling unification (see Ref. [80] for a dedicated analysis). Out of all the representations listed in Table III, we find that only Q ∼ (3, 2, 1/6) can considerably improve unification with respect to the SM while, at the same time, keeping the unification point at a reasonably high scale Λ GUT ∼ 10 15 GeV. This possibility was already pointed out in [80], the only difference is that in our two-loop analysis the value of the optimal threshold m Q = 2 × 10 7 GeV is about a factor of twenty larger than what found in [80].
While gauge coupling unification is a desirable feature for any particle physics model, envisaging a GUT completion of KSVZ axion models featuring a hierarchy V a Λ GUT in which only the fragment R Q of a complete GUT multiplet receives a mass m Q < ∼ V a Λ GUT , while all the other fragments acquire GUT-scale masses, is not straightforward. This appears especially challenging in all the cases in which U (1) P Q commutes with the GUT gauge group. Besides these theoretical considerations, we must also consider the possibility that improved gauge coupling unification might simply occur as an accident because of the many representations we have considered. We then conclude that also improved unification is not a sufficiently well motivated condition to be chosen as a selection criterium for preferred axion models.

E. Summary
The fifteen R Q 's that satisfy our two criteria are collected in Table III. In this table we also give, for each R Q , in the third column the energy scale Λ LP at which the first LP occurs, together with the corresponding gauge coupling, in the fourth column the value of E/N which determines the strength of the axion-photon coupling, and in the last column the number of DW. It should be clear that the two criteria on which our selection has been based should not be understood as strict no-goes, since under specific conditions models that do not satisfy these conditions can also be viable. For example, the first requirement of sufficiently fast decays for the strongly interacting relics applies only to scenarios for which m Q < T reheating , and there is no similar issue in pre-inflationary scenarios. However, in section VI we will show that, as long as the threshold for integrating out the heavy Q is kept at, or below, m Q ∼ 5 × 10 11 GeV, the LP condition alone is sufficiently constraining that the previous results get only mildly relaxed. If the PQ breaking scale is allowed to lie sizeably above 5 × 10 11 GeV, as it is possible in pre-inflationary models, then also the LP condition will progressively diminish its strength due to the possibility of increasing correspondingly the heavy quark threshold. This however, can be done only at the cost of an increased fine tuning in the initial value of θ, something that can well be considered as theoretically unpleasant.

V. AXION COUPLING TO PHOTONS
From the experimental point of view, the most promising way to unveil the axion is via its interaction with photons, which is described by the effective term L aγγ = −(1/4)g aγγ a F ·F . The axion-photon coupling is given in terms of the anomaly coefficients in Eq. (2) by [30,92]: where the uncertainty is evaluated with the NLO chiral Lagrangian [93]. The strongest coupling is obtained for R s Q = (3, 3, −4/3) that gives E s /N s − 1.92 ∼ 12.75, almost twice the usually adopted value of 7.0 [53], while the weakest coupling is obtained for R w Q = (3, 2, 1/6) for which E w /N w − 1.92 ∼ −0.25 is about 3.5 times larger than the usual lower value of 0.07. The corresponding couplings are depicted in Fig. 3 with the two oblique black lines labeled E/N = 44/3 and E/N = 5/3. According to our two selection criteria all preferred hadronic axion models containing a single R Q fall within the light green strip enclosed by these two lines.
Let us now study to which extent the previous results can be changed by the presence of more R Q 's. It would be quite interesting if, for example, g aγγ could receive sizeable enhancements. However, we can easily see that, as long as the sign of ∆X = X L − X R is the same for all R Q 's, this cannot occur. Let us write the combined anomaly factor for R Q + R s Q : Since by construction the anomaly coefficients of any R Q in our preferred set satisfy E/N ≤ E s /N s , the factor in parenthesis cannot be larger than unity, implying E c /N c < E s /N s . This is not so, however, if we allow for opposite signs in the PQ charge differences: ∆X = −∆X s . In this case E/E s and N/N s become negative and g aγγ can get enhanced. The largest enhancement attainable with two R Q 's is obtained with R s Q ⊕R w Q . This still respects the LP selection criterium and yields E c /N c = 122/3. Further enhancements are possible with three or even more R Q 's, but adding multiple R Q 's will eventually lead to sub-Planckian LP, so that there is in fact an absolute upper bound on E c /N c compatible with the requirement of no LP below 10 18 GeV. We have searched for this maximum coupling, and we have found that it corresponds to the combination R 8 ⊕ R 6 R 9 , where the meaning of the symbols ⊕ and is that the representation has to be taken with the same or with the opposite sign of the PQ charge difference ∆X . The resulting maximum axion-photon coupling corresponds to E/N = 170/3 and is depicted in Fig. 3 with the uppermost dotdashed oblique line.
limit E/N = 5/3 for a single R Q , and even yield complete axion-photon decoupling (within theoretical errors), a possibility that requires an ad hoc choice of R Q 's, but no numerical fine tuning. With two R Q 's there are three such cases: R 6 ⊕ R 9 ; R 10 ⊕R 12 and R 4 ⊕R 13 giving respectively E c /N c = (23/12, 64/33, 41/21) ≈ (1.92, 1.94, 1.95). In all these cases the axion could be detected more easily via its coupling to nucleons, providing additional motivations for axion searches which do not rely on the axion coupling to photons [94][95][96]. 2 Finally, let us mention that in the cases in which g aγγ is strongly suppressed, the limits from stellar evolution are accordingly weakened. However, the region m a 1 eV is still excluded due to the hot DM bound [53].

VI. KSVZ MODELS IN PRE-INFLATIONARY SCENARIOS
The discussion of the KSVZ models in the previous sections pertained to post-inflationary scenarios, with m Q < T reheating ensuring, as initial condition, a thermal abundance for the Q. However, the scale of inflation is firmly bounded from below only by BBN considerations, which imply a loose limit of just a few MeV, and thus m Q T reheating is certainly not an unlikely possibility. It is then mandatory to explore to which extent our results can represent acceptable estimates of the preferred axion window also in pre-inflationary scenarios, when the first condition of forbidding long lived or stable strongly interacting relics must be dropped.
The requirement that the contribution to the cosmological energy density from axion misalignment does not exceed the energy density of DM implies, in post-inflationary scenarios, f a < ∼ 5·10 11 GeV. In preinflationary scenario this condition can be avoided by assuming that in the original patch corresponding to the present observable Universe the initial value of θ is sufficiently close to the minimum of the zero temperature axion potential. However, values largely in excess of f a ∼ 5 · 10 11 GeV require correspondingly large fine tunings in the initial conditions, or invoking anthropic selection arguments to justify a sufficiently small initial value of θ. This might well be considered an unwanted feature of preferred axion models, and therefore we will restrict our study of the pre-inflationary case by keep-ing the condition f a < ∼ 5 · 10 11 GeV. Taking also in this case a maximum value for the Yukawa couplings y Q ≤ 1, we can again apply the LP criterium with m Q = 5 · 10 11 GeV as the threshold for integrating out the heavy quarks.
To see which constraints can be implied by the LP condition alone, let us start by considering a single representation R Q = (C Q , I Q , Y Q ). The E/N factor can be conveniently written as where C Q is the dimension of the color representation, and T (C Q ) is the color Dynkin index previously introduced. In terms of the Dynkin indices labeling the representation C Q = (α 1 , α 2 ) the index can be written as [97]: P (α 1 , α 2 ) = α 2 1 + 3α 1 + α 1 α 2 + 3α 2 + α 2 2 . (22) The polynomial P (α 1 , α 2 ) which appears in the denominator of E/N has its minimum value for the fundamental representation (α 1 , α 2 ) = (1, 0), so we learn that the largest values of E/N are obtained for a color triplet.
In order to study the values of E/N in (I Q , Y Q ) representation space, we start by establishing the 'corners' that saturate the LP condition. Respectively for hypercharge and weak isospin we find: Any larger value of hypercharge or isospin would induce a sub-Planckian LP in the corresponding coupling (although in our search we allow for continuous values of Y Q , for convenience we round the result to a close fractional value). We can now use eq. (20) to find the maximum of E/N subject to the condition Y Q I Q ≤ 5 2 , which is implied by the maximum allowed coefficient for the hypercharge coupling βfunction. This value is given by the value of I Q that maximizes the function In the case of more representations, the maximum value of E/N can be found in correspondence of the 'corners' of the (I Q , Y Q ) representation space (as well as for sets of representations that are 'equivalent' in the sense specified below). The combination (3, 1, 5/2) ⊕ (3, 4, 0) which maximizes the numerator of E/N requires the addition of SU (2) L × U (1) Y singlet and color non-trivial representations, in order to minimize the denominator to the minimum possible value of ±1/2. With three representations, adding (8, 1, 0) (with a negative sign of the PQ charge difference) allows to arrange for this given that (1 + 4) · T (3) − T (8) = −1/2. For this combination of three R Q we then obtain E/N = −135/2. With four representations, the +1/2 in the denominator can be obtained by including ⊕ (3, 1, 0) (6, 1, 0). This gives E/N = +135/2 which results in a slightly smaller axion-photon coupling, due to the negative sign of the chiral perturbation contribution, see eq. (18). Equivalent representations can be obtained for example by the replace- , as well as in other similar ways. To check that the result obtained relying on the previous simple argument indeed corresponds to the maximum (E/N ), we have carried out a thorough computer search exploring (I Q , Y Q ) representation space, which confirmed that |E/N | = 135/2 gives the maximum axion-photon coupling compatible with the LP condition. The corresponding upper limit on g aγγ is only about 20% larger than the maximum coupling labeled E/N = 170/3 obtained in post-inflationary scenarios, and well below the limit on DFSZ-IV models represented by the uppermost red line. Not to clutter too much the plot, we have not included in Fig. 3 the corresponding line.
In conclusion, for KSVZ models the preferred axion window for the different cases considered is well represented by the black oblique lines in Fig. 3, restricted to the region on the right hand side of the violet vertical line labeled f a > 5 × 10 11 GeV. On the left of this line only pre-inflationary models with progressively larger values of f a are allowed. In this case the heavy quark threshold can be correspondingly increased, thus weakening the constraints from the LP condition. Therefore for KSVZ models larger values of the axion-photon coupling become allowed within this region. However, this goes at the expense of a progressively larger amount of fine tuning in the initial value of θ, which might well be considered as an unwanted feature in phenomenologically preferred axion models.

VII. DFSZ-TYPE OF AXION MODELS
In DFSZ-type of models [36,37] two or more Higgs doublets H i , carrying PQ charges, together with the SM singlet axion field Φ are introduced. The SM fermion content is not enlarged, but in general both quarks and leptons carry PQ charges. The electromagnetic and color U (1) P Q anomalies then de-pend on the known fermions assignments under the SM gauge group, but also on their model dependent PQ charge assignments. Hence, several variants of DFSZ axion models are possible, some of which have been discussed, for instance, in Refs. [31,32]. Here we argue that for most of these variants the axionphoton coupling falls within the regions highlighted in Fig. 3. Only in some specific cases the KSVZ upper limit E/N = 170/3 can be exceeded. We will point out under which conditions this can occur.
Let us start with some general considerations: we assume n H ≥ 2 Higgs doublets H i which are coupled to quarks and leptons via Yukawa interactions, and to the axion field Φ through scalar potential terms. The kinetic term for the scalars carries a U (1) n H +1 rephasing symmetry that must be explicitly broken to U (1) P Q × U (1) Y in order that the PQ current in Eq. (2) is unambiguously defined, and to avoid additional Goldstone bosons with couplings only suppressed as the inverse of the electroweak scale. By considering from the start only gauge invariant operators, the relevant explicit breaking U (1) n H +1 → U (1) P Q must be provided by non-Hermitian renormalizable terms in the scalar potential involving H i and Φ. This implies that the PQ charges of all the fermions and Higgs doublets are interrelated and cannot be chosen arbitrarily. In the most general scenario, each SM fermion field carries a specific PQ charge. However, given that the anomalies of the PQ current depend on the difference between the PQ charges of L-and R-handed fermions, without loss of generality we can set the PQ charges of the L-handed fermions to zero, and only consider the charges of the R-handed fermions X uj , X dj , X ej , where j is a generation index. The ratio of anomaly coefficients E/N reads and it is particularly convenient to write it as in the second equality. Note that in order to have a nonvanishing PQ-color anomaly, the denominator must be non-vanishing. The original DFSZ model [36,37] includes two Higgs doublets, H u,d , coupled to the singlet scalar field via the quartic term H u H d Φ 2 , and family independent PQ charges for the SM fermions. Then the factor E/N is fixed up to the two-fold possibility of coupling the leptons either to H d or to H * u . Eq. (26) shows that these two cases yield, respectively DFSZ-I : X e = X d , E/N = 8/3 , DFSZ-II : which in both cases give axion-photon couplings that fall inside the KSVZ band in Fig. 3.
Let us now consider the so called DFSZ-III variant [31] in which the scalar sector is enlarged to con-tain n H = 3 Higgs doublets H e,d,u coupled respectively to leptons, down-type and up-type quarks. Although here we have some more freedom in choosing the values of the charges X e , in order to enforce the breaking U (1) 4  It is easy to verify that the bilinear terms alone yield the same two possibilities listed in Eq. (27). Let us then consider quadrilinear couplings. Since Φ 2 has the same PQ charge than (H u H d ) † , the four cases below exhaust all the possible relations between X e and the other PQ charges: These four possibilities yield, respectively: all of which give axion-photon couplings that fall within the N Q = 1 band in Fig. 3. 3 Many more possibilities in choosing the PQ charges become possible if we allow for generation dependent assignments, as was done for example in Ref. [98]. The maximum freedom corresponds to the case in which there are three Higgs doublets for each fermion species (H e1 , H e2 , H e3 , etc.) so that the scalar rephasing symmetry is U (1) n H +1 with n H = 9. Here, we refer to this scenario as DFSZ-IV. Although such a model might be plagued by various phenomenological issues, bounding from above the maximum possible E/N in DFSZ-IV is useful, since it provides an upper bound to E/N for all cases with generation dependent PQ charges, and with n H ≤ 9 Higgs doublets coupled to the SM fermions.
From Eq. (26) we see that in order to maximize E/N we have to find the maximum possible value of j (X uj + X ej ) (namely the largest possible PQ charges for the up-type quarks and leptons, all with the same sign) together with the minimum value of the denominator j (X uj + X dj ) compatible with a nonzero QCD anomaly, which is 2X Φ . This second condition is realized, without loss of generality, by choosing The last equality is satisfied by where the values of y and z are arbitrary. The scalar terms allowed by this choice break the U (1) 4 symmetry in the second and third generations to U (1) y ×U (1) z , which in turn must be broken by couplings between these scalars and scalars of the first generation. Starting with the second generation, the term H d2 H u1 (Φ * ) 2 yields the largest possible charge X u2 = y = 6X Φ (and X d2 = −4X Φ ). Note that the term H u2 H * u1 Φ 2 would instead only yield y = 4X Φ . The relatively large charge X u2 allows to get an even larger charge X u3 via the term (H u3 H * u2 )(H * d1 H * u2 ) giving X u3 = z = 12X Φ and X d3 = −14X Φ . This accomplishes the breaking of all the redundant symmetries in the quark sector. Regarding the breaking of the U (1) 3 symmetries in lepton sector, we need to couple at least one lepton scalar (H e1 without loss of generality) to the scalars of the quark sector. The possible bilinears are either of the form (H e1 H uj ) or (H e1 H * dj ). The most favorable possibility to get a large charge X e1 is to start with (H e1 H * d1 ), since H d1 has the only non-negative charge X d1 = 0, and next to couple this term to the bilinear with the largest possible positive charge, which is (H d3 H * d1 ). This yields X e1 = 14X Φ , which can be used to push up X e2 and X e3 to even larger values, via the following sequence of couplings: , yielding X e2 = 92. The values of the PQ charges derived in this way give the maximum possible axion-photon coupling in DFSZ-IV models, which corresponds to DFSZ-IV : (E/N ) max = 524/3 .
In this class of models it is also easy to obtain axion-photon decoupling ensuring at the same time a correct breaking of the U (1) 9+1 global symmetries down to U (1) P Q . An example is given by: which yields E/N = 23/12 ≈ 1.92. In conclusion, although the value in Eq. (32) exceeds by a factor of three the maximum KSVZ value E/N = 170/3, the construction through which (E/N ) max has been obtain is sufficiently cumbersome to suggest that the N Q > 1 region in Fig. 3 can be considered as representative also of most of DFSZ-IV models. The values of E/N associated to the maximum and minimum of g aγγ for different classes of models are summarized in Table IV. Note that differently from the KSVZ models analyzed previously, the limits on the axion-photon coupling in DFSZ models do not depend on details of the Universe cosmological evolution, and therefore hold also within the region on the left of the violet vertical line in Fig. 3 labeled f a > 5 × 10 11 GeV. As we have seen, the criterium of the absence of LP plus the requirement of no cosmological issues in post-inflationary scenarios, or of natural initial values for θ in pre-inflationary scenarios, imply well defined limits for the axion-photon coupling in all the type of models we have considered so far. However, it is also possible to envisage scenarios in which our selection criteria are satisfied, and still the axion-photon coupling can lie well above the preferred window. The models we have considered in our analysis are characterized by a specific sector of scalar fields carrying PQ charges. In KSVZtype of models we have included only one SM singlet scalar Φ. In DFSZ-type of models we have allowed for up to one scalar doublet for each type of SM fermion, for a total of nine electroweak doublets carrying PQ charges, in addition to the singlet Φ. However, more PQ charged singlets Φ k could be introduced without conflicting with phenomenological constraints, and up to about fifty electroweak scalar doublets H k could be also added before violating the LP condition. By adding scalar doublets that do not couple directly to the fermions, it is possible to obtain very large PQ charges for the leptons, and huge enhancements in the value of E/N . To see how this can work let us start from the quadrilinear scalar coupling H u H d Φ 2 and the PQ charges X (H u ) = −2X Φ and X (H d ) = 0. Let us define H 1 = H u and next let us add a whole set of scalar doublets H k (k = 2, 3, . . . , n) with charges . Finally, let us couple the lepton Higgs doublet as (H e H n )(H n H d ). We obtain X (H e ) = 2 n+1 X Φ . As mentioned above, n can be as large as fifty before a LP is hit below the Planck scale, so that exponentially large lepton charges |X e | ∼ 2 50 are possible (in this construction the axion-electron coupling gets exponentially enhanced as well). An analogous construction is possible also in KSVZ models by adding more scalar singlets Φ k . This possibility was put forth in [99] to which we refer for further details.

IX. CONCLUSIONS
Axions are well-motivated candidates for physics beyond the SM. Axion models solve the strong CP problem and provide an excellent DM candidate. Most importantly, experiments are starting now to probe the parameter space region for the axionphoton coupling predicted by realistic axion models, and an outburst of new experimental proposals for axion searches have been put forth in last few years (see e.g. [91,[94][95][96][100][101][102]). It is then very important to define precisely the region in parameter space where axion searches should focus, possibly assessing which are the desirable properties common to axion models that fall within this region. The commonly adopted axion window considered so far [53] corresponds to a selection of realistic KSVZ and DFSZ axion models, e.g. those of Refs. [30][31][32]; but, lacking of a well-defined guiding principle, it unavoidably contains some degree of arbitrariness.
In this work we have put forth a recipe for defining a window for preferred axion models on the basis of well defined sets of selection criteria. We have considered first KSVZ-type of axion models, for which all the particles carrying PQ charges are new, beyond the SM states. Starting with post-inflationary scenarios we have classified the representations R Q of the new quarks Q on the basis of the following two criteria: (i) Q decays must be fast enough in order not to bring in cosmological issues; (ii) the new representations R Q should not generate Landau poles below 10 18 GeV. Only fifteen representations which we have collected in Table III satisfy these two selection criteria. The ratio of their anomaly coefficients E/N can then be used to define a first window for preferred axion models, which is displayed in Fig. 3 in the m a -g aγγ plane. We have then shown that models containing multiple R Q representations, but which still satisfy the two criteria, allow to enlarge sizeably this window, and that at fixed values of m a both stronger and weaker couplings are possible. While the weakening of the coupling can reach the limit of complete axion-photon decoupling (within current theoretical errors), the size of the possible enhancements is still bounded by an upper limit on E/N , which is set by the LP condition. We have also discussed the possibility of considering additional criteria, like requiring the absence of cosmologically dangerous domain walls, or considering the possible improvements in SM gauge coupling unification induced by R Q , but we have concluded that these conditions are not sufficiently strong and general to represent valid selection criteria, and should just be considered as desirable features of specific models.
We have then extended the analysis to include also KSVZ models in pre-inflationary scenarios, in which case the first condition (i) must be dropped. The second condition (ii) on the absence of sub-Planckian LP still holds, and maintains in full its constraining power under the assumption that the threshold for integrating out the new heavy quarks remains at, or below, m Q ∼ 5 · 10 11 GeV. This corresponds to axion models in which the QCD θ parameter can assume natural initial values, of order unity, without generating an overabundance of DM. We have shown that, with respect to the case when also condition (i) can be consistently applied, the upper limits on the axion-photon couplings in preinflationary scenarios are relaxed by a factor of 2.5 in the case of a single R Q , and only by 20% in the case of multiple R Q .
Finally, we have argued that the definition of a preferred axion window based on the analysis of post-and pre-inflationary KSVZ-type of axion models, is also representative of the vast majority of realistic DFSZ-type of models. The minimal DFSZ realization contains two Higgs doublets, and the axionphoton coupling is fixed up to a two-fold choice. The next-to-minimal realization includes one additional Higgs doublet coupled to the leptons. Since leptons contribute only to the anomaly coefficient E, this allows for larger values of E/N . Nevertheless, for all these cases the axion-photon coupling falls within the window for post-inflationary KSVZ models with a single R Q . A much more general (and probably unrealistic) possibility includes nine Higgs doublets, each coupled to a different SM fermion. The maximum E/N allowed in this extreme case still exceeds the limit for post-inflationary hadronic axion models with multiple R Q representations by only a factor of three. Sizeable enhancements of the axion-photon couplings seem to become possible only in rather cumbersome constructions which introduce a large number of scalars carrying PQ charges (electroweak singlets for KSVZ models [99] or doublets for DFSZ models) which do not couple to the SM fermions or to the Q, but only among themselves through a quite specific pattern of mixed operators.
(e.g. Q = 1/5, π, etc.). They cannot decay into SM particles in force of electric charge conservation. Since quarks of this type do not get confined into hadrons of integer charge, they also cannot get bounded into neutral hadrons, atoms or molecules. In other words, their fractional charge must remain naked.
In case the Q's have less exotic (for example integer) electric charges, the connection with their absolute stability is less obvious. However, also in this case it is possible to reach the same conclusion. Let us first consider the fundamental particles carrying color of the SM. Let us assign SU (3) C triality τ = +1 to the fundamental representation q j ∈ 3 for the quarks. Then the reducible representation 3 n = 3 × 3 × 3 . . . , as well as all its irreducible fragments, have triality τ = n [mod 3]. For example the antisymmetric 2-index representation for the antiquarks q jk ∈ 3 has triality τ = 2 since it is an irreducible fragment of 3 × 3 = 3 + 6 (also containing the 2-index symmetric 6), while the 3-index q jkl ∈ 3×3×3 containing the totally antisymmetric singlet, the totally symmetric 10, and two adjoints 8 where the gluons sit, has triality τ = 3 = 0 [mod 3]. SM hadrons are color singlets and are built by contracting SU (3) C indices with the totally antisymmetric tensor ijk into invariant index-less tensors. We can then build hadrons only from combinations of quarks and gluons of 0-indices [mod 3], e.g. q i q j q k , q i q jk , q i q j q k q lmn etc.. The SM fundamental colored particles have charge Q(q j ) = −1/3 + n (n = 0, 1) for 1-index (quarks), Q(q ij ) = −2/3 + n for 2-index (antiquarks) and Q(q ijk ) = 0 for 3-index (gluons). Therefore, any combination of a number of indices multiple of three results in an integrally charged or neutral state. So, as is experimentally well known, all SM hadrons, being color singlets, are integrally charged. Electric charge conservation then precludes the possibility that fractionally charged Q-hadrons of any type could decay into lighter SM states.
It is in fact possible to prove a slightly stronger statement: gauge invariant operators inducing decays of exotic heavy quarks Q are allowed if and only if all color singlet Q-hadrons are integrally charged or neutral. The necessary condition in this statement is equivalent to the previous result (which can also be proven in a more direct way and without appealing to electric charge conservation). As regards the sufficient condition, it is not very useful since the decay operators which will mandatorily appear could be so suppressed to render the Q's effectively stable with respect to all phenomenological consequences. This can happen for example if we choose R Q 's with particularly large isospin/hypercharge values, since in this case gauge invariant decay operators would arise at rather high dimensions. The proof of the sufficient condition is then uninteresting and can be omitted.