Anomaly Ratio Distributions of Hadronic Axion Models with Multiple Heavy Quarks

We consider hadronic axion models that extend the Standard Model by one complex scalar field and one or more new heavy quarks, i.e. $N_\mathcal{Q} \geq 1$. We review previously suggested selection criteria as well as categorize and catalog all possible models for $N_\mathcal{Q} \leq 9$. In particular, allowing for $N_\mathcal{Q}>1$ can introduce models that spoil the axion solution of the strong CP problem. Demanding that Landau poles do not appear below some energy scale limits the number of preferred models to a finite number. For our choice of criteria, we find that $N_\mathcal{Q} \leq 28$ and only 820 different anomaly ratios $E/N$ exist (443 when considering additive representations, 12 when all new quarks transform under the same representation). We analyze the ensuing $E/N$ distributions, which can be used to construct informative priors on the axion-photon coupling. The hadronic axion model band may be defined as the central region of one of these distributions, and we show how the band for equally probable, preferred models compares to present and future experimental constraints.


I. INTRODUCTION
QCD axions [1,2], initially proposed as a solution to the strong CP problem [3,4], are excellent cold dark matter (DM) candidates [5][6][7][8][9]. Numerous experimental searches are currently underway to find such particles [10]. One major challenge of axion detection is that the axion mass is set by an unknown parameter, the axion decay constant f a , which can range across many orders of magnitude. Moreover, the axion's interactions with the Standard Model (SM) are usually model-dependent, and a UV axion model has to be constructed in order to determine the exact relationship of f a and the axion couplings.
One class of such UV models are hadronic (also called KSVZ-type) axion models [11,12], which extend the SM by a new complex scalar field and N Q ≥ 1 heavy, exotic quarks. For a given value of N Q there exist multiple, discrete models, which trace out lines in the axion mass and axion-photon coupling parameter space. The locations of these lines are determined by the anomaly ratio E/N and a model-independent contribution from axionmeson mixing.
To map and restrict the resulting landscape of axion models, it has been suggested that phenomenological selection criteria can be used to single out preferred models [13,14]. This allows us to restrict the parameter space and helps experiments to assess their sensitivity requirements. However, so far only the case of N Q = 1 has been fully cataloged, which is why we want to study models with N Q > 1 as far as this is feasible. First, we summarize the construction of KSVZ-type axion models and phenomenological selection criteria in Secs. II and III. Subsequently, a catalog of all possible models with N Q ≤ 9 is presented and the resulting E/N distributions are discussed. We catalog all preferred models, * vaisakh.plakkot@stud.uni-goettingen.de † hoof@uni-goettingen.de for which we find that the maximum possible number of Qs is N Q = 28. In Sec. V we outline how the catalog of models can be used to construct informative prior distributions on E/N . These can be used to define the KSVZ axion model band and we show how it compares to current and future experimental constraints. Finally, we summarize our work and end with some closing remarks. Model catalogs and further supplementary material are available on Zenodo [15].

II. HADRONIC AXION MODELS
Let us denote a representation of a particle as (C, I, Y), where C and I are the SU(3) C color and SU(2) I isospin representations, respectively, while Y denotes the particle's U(1) Y hypercharge.
For example, the traditional KSVZ axion model contains a heavy chiral quark Q = Q L + Q R ∼ (3, 1, 0), charged under the U(1) PQ Peccei-Quinn (PQ) symmetry with charge X = X L − X R = ±1, and the complex scalar field Φ ∼ (1, 1, 0) with PQ charge normalized to X Φ = 1. All SM fields are uncharged under the PQ symmetry in the KSVZ model, and the relevant part of the Lagrangian is where y Q is the Yukawa coupling constant and the last term is a potential for the complex scalar field with order parameter v a . The Lagrangian is invariant under a chiral U(1) PQ transformation Φ → e iα Φ, Q L/R → e ±iα/2 Q L/R . The field Φ attains a non-zero value at the minimum of the potential, resulting in a spontaneously broken PQ symmetry. Expanding Φ around its vacuum expectation value gives the axion as the corresponding angular degree of freedom, with value in the interval [0, 2πv a ). The mass of Q is then m Q = y Q v a / √ 2.
Performing a chiral U(1) transformation such that Q L/R → e ±ia/(2va) Q L/R , the mass term for Q can be made independent of the axion field phase. This transformation adds an anomalous G G term to Eq. (1) as well as an F F term, where G and F are the gluon and photon field strength tensors, respectively, and the tilde denotes their duals. With the electromagnetic (EM) and color anomaly contributions due to the U(1) PQ charged quarks labeled E and N respectively, the coupling terms become where f a = v a /(2N ). The axion-photon coupling is thus parameterized by the anomaly ratio E/N alone. More precisely, the mass and coupling to photons for QCD axion models are given by [16,17] For some representation r under which the heavy quark Q in the KSVZ axion model transforms, the EM and color anomalies can be calculated as where d(·) denotes the dimension of a representation, q = I (3) − Y is the EM charge of Q, and T (C) is the SU(3) C Dynkin index (see Ref. [18]). In KSVZ-type models, only Q is charged under the PQ symmetry (apart from Φ) and e.g. for Q ∼ (3, 1, 0) we have N = X /2 and E = 3X tr(q 2 ), using that T (3) = 1/2. In general, one finds for a single Q that where the last equality holds only when Q is a singlet under SU(2) I . This e.g. leads to the well-known result that the original KSVZ model has E/N = 0. When considering models with multiple Q i , which have representations r i and anomaly coefficients E i and N i given by Eqs. (5a) and (5b), respectively, the overall anomaly ratio is simply where the index i runs over the different quarks, labeled i = 1, . . . , n. Note that, when labeling a tuple of Qs in a model, there exists a "relabeling symmetry." For example, assume that two Qs with the same U(1) PQ charge respectively transform under representations r 1 and r 2 , denoted by r 1 ⊕ r 2 . Then there is an equivalency relation such that r 1 ⊕ r 2 ∼ r 2 ⊕ r 1 , in the sense that they trivially give the same anomaly ratio E/N . Similarly, we can also consider combinations of representations with " ", the symbol we use to denote Qs with opposite U(1) PQ charges such that r i r j ⇒ X i = −X j . Here we have e.g.
, as all three models trivially give the same overall anomaly ratio.
The relabeling symmetry allows us to simplify the presentation of the catalog, and we refer to a list of models where this symmetry has been accounted for as "nonequivalent." It may also play a role in the statistical interpretation of the catalog: if not all Qs are indistinguishable, the multiplicity arising from the equivalency relation must be taken into account. We comment on this further in Section V A.

III. PHENOMENOLOGICAL SELECTION CRITERIA
Let us now review the various selection criteria for preferred axion models, most of which have already been proposed and discussed extensively in Refs. [13,14].
Here, we focus on the applicability in the pre-and postinflationary PQ symmetry breaking scenarios and observe that N Q > 1 allows for the existence of a new criterion related to the axion's ability to solve the strong CP problem.

A. Dark matter constraints
A natural requirement is to demand that axions do not produce more DM than the observed amount, Ω c h 2 0.12 [19]. For QCD axions this results in an upper bound on f a and previous studies of preferred axion models used f a < 5 × 10 11 GeV [13,14], assuming a post-inflationary cosmology with realignment axion production. Let us extend this discussion and make a few comments regarding the different cosmological scenarios and their impact on the f a bound.
First, in the pre-inflationary PQ symmetry breaking scenario, the initial misalignment angle of the axion field, denoted by θ i , is a random variable. Since any topological defects are inflated away, realignment production is the only relevant contribution and the limit on f a depends on its "naturalness." While this is not a uniquely defined concept, using the usual assumption of uniformly distributed angles, θ i ∼ U(−π, π) the code developed in Ref. [20] finds f a < 4 × 10 12 GeV for the 95% credible region of posterior density. 1 This limit on f a effectively relies on the naturalness being encoded automatically in prior on θ i . Second, when topological defects can be neglected in the post-inflationary symmetry breaking, the relic axion density is determined by an average of misalignment angles over many causally-disconnected patches. This corresponds to the benchmark scenario of Refs. [13,14]. Again using the code developed in Ref. [20], we obtain f a < 2 × 10 11 GeV (at the 95% CL).
The third and last case is the post-inflationary scenario including a significant contribution from topological defects i.e. cosmic strings and domain walls (DWs). In fact, recent studies indicate that the production of axions via topological defects dominates the vacuum realignment production [21,22]. For models with domain wall number N DW ≡ 2N = 1 (cf. Section III E), the authors find that f a 10 10 GeV, while models with N DW > 1 reduce the value of f a by a factor O(N DW ) [22]. For the preferred models considered in this work, N DW ≤ 28 such that the bound might be loosened to about f a 3 × 10 8 GeV. It should be noted that these results rely on extrapolating the outcome of numerical simulations more than 60 orders of magnitude, and they hence are potentially subject to large systematic uncertainties.
In summary, the upper limit on f a , and hence the results presented in what follows, very much depend on the cosmological scenario at hand. To simplify the discussion, to avoid the potentially large uncertainties mentioned above, and to better compare with previous work of Ref. [14], we also adopt f a < 5 × 10 11 GeV.
However, we stress again that a different choice of f a will affect the number of preferred models, as f a is one of the factors that determines the value of m Q . This is because provides an upper bound on m Q . Moreover, a universal bound on the m Q (up to the Yukawa couplings) requires that all Qs are coupled to the Φ field in the same way to get a single v a parameter. So long as the coupling y Q ∼ O(1) or lower, the upper bound on f a = v a /N DW is indeed an upper limit to m Q . Larger values of the coupling require fine-tuning of parameters, and are hence deemed undesirable from a theoretical viewpoint. In what follows, we choose m Q = 5 × 10 11 GeV as a conservative value for all Q masses (see Sec. III D for more details on the influence on Landau pole constraints). Finally, note that the Qs themselves contribute to the matter content in the Universe, and we need to consider the possibility that their abundance exceeds Ω c h 2 . Since this issue can be avoided if the lifetime of the quarks is short enough, we discuss this in the next section.

B. Lifetimes
Other than the possibility that the Qs' abundances exceed Ω c h 2 , there also exist additional experimental and observational constraints, which have already been discussed before [13,14].
To avoid the DM constraints, we require the Qs to decay into SM particles with a reasonably low lifetime. Heavy quarks with m Q 1 TeV and lifetimes 0.01 s < τ Q < 10 12 s are severely constrained, as they would also affect Big Bang Nucleosynthesis and observations of the Cosmic Microwave Background [23,24]. Fermi-LAT excludes 10 13 s < τ Q < 10 26 s, thus excluding lifetimes greater than even the age of the Universe (∼ 10 17 s) [25]. As a result, for heavy quarks (m Q 1 TeV), only representations with τ Q < 10 −2 s are considered to be a part of the preferred window. Lighter relics would be excluded from experimental bounds e.g. at the LHC [26].
Such a constraint on the Q lifetime, when applied to the heavy quark decay rate, translates to restrictions on the dimensionality of the possible Q to SM fermion decay operators. With m Q 5 × 10 11 GeV, the lifetime constraints in turn constrain operators to have dimensions d ≤ 5 [13,14]. This implies a total of 20 possible representations for Q, all charged under SU(3) C and U(1) Y . The lifetime constraint has no further consequence on cases with N Q > 1 under the assumption that the different Q i do not interact among themselves or decay into particles other than SM fermions.
As noted before [14], the lifetime constraints are typically not required in the pre-inflationary PQ symmetry breaking scenario. This is because the Qs can get diluted by inflation, which prevents them from becoming cosmologically dangerous relics after they freeze out. Without these constraints, many more models with even higherdimensional operators can exist, and restricting ourselves to at most five-dimensional operators therefore only becomes an assumption in this case.

C. Failure to solve the strong CP problem
This criterion is specific to models with N Q > 1 that allow the Qs to have opposite U(1) PQ charges. It is clear from Eq. (5b) that the addition of multiple heavy quarks can lead to a smaller overall N than the individual N i , but only when one or more of the quarks have a (relative) negative U(1) PQ charge. In some cases a total cancellation of the N i terms occurs (N = 0). While these models give rise to massless axion-like particles with a coupling to photons governed by E, they do not solve the strong CP problem: as can be seen from Eq. (2), N = 0 means that there is no G G contribution in the Lagrangian. Considering that the primary objective of QCD axion models is to solve the strong CP problem, we propose that only models with N = 0 should be considered preferred.

D. Landau poles
The single most powerful criterion amongst the ones proposed by Refs. [13,14] in the context of this work comes from the observation that representations with large C, I, or Y can induce Landau poles (LPs) at energies well below the Planck mass. At an LP, the value of a coupling mathematically tends to infinity, signaling a breakdown of the theory. Since quantum gravity effects are only expected to appear at energies near the Planck mass, a breakdown of the theory before that point can be regarded as problematic or undesirable.
It has thus been proposed that preferred models have LPs at energy scales Λ LP 10 18 GeV. From the 20 representations mentioned previously, only 15 fulfil this criterion [13,14]; we refer to these as "LP-allowed" models and label them r 1 to r 15 (as per Table II in Ref. [13]).
The running of the couplings are computed at two-loop level with the renormalization group equation [27,28] where with i, j ∈ {1, 2, 3} for the three gauge groups, α i = for energy scale µ and Z boson mass m Z , while a i and b i are the one-and two-loop beta functions. C 2 and T are the quadratic Casimir and Dynkin indices of the corresponding gauge group, respectively, and F and S denote fermionic and scalar fields.
G i denotes the adjoint representation of the gauge group, and κ = 1 2 , 1 for Weyl and Dirac fermions, while η = 1 for complex scalars. 2 Adding multiple Qs to the theory increases the coefficients of beta functions through the fermionic terms. As a consequence, the couplings diverge faster i.e. induce LPs at lower energy scales, as has been anticipated before [14].
Since the addition of more particles with a given representation into a gauge theory only worsens the running of the corresponding gauge coupling, it is possible to find the number of copies of a particle that can be included in the theory before it induces an LP below 10 18 GeV. This drastically reduces the number of LPallowed combinations possible. Integrating all Q i in at m Q = 5 × 10 11 GeV, we find that there are 59,066 nonequivalent combinations of Q i from the representations r 1 , r 2 , . . . , r 15 that do not induce LPs below 10 18 GeV.
As the Q i contribute to the beta functions above the 2 The case of η = 1 2 for real scalars is not relevant for the present study. Also note that the expression for b ij in Ref. [28] is slightly erroneous since the second term applies only to the non-diagonal elements of b ij , as found when comparing with the SM beta functions in Ref. [27]. energy scale m Q , the running of the gauge coupling begins to deviate from the SM only at this scale. Different values of m Q are bound to produce different results for the LPs; the lower m Q is, the earlier an LP appears. As an example, consider m Q = 10 10 GeV: for N Q = 3, we find that 888 models are preferred (they have Λ LP > 10 18 GeV and N = 0 as per the discussion in III C), compared to 1,442 models when m Q = 5 × 10 11 GeV. Furthermore, we use the same mass for all Q i in the models, which may not be the case in reality (due to different y Qi or e.g. in multi-axion models). However, setting the masses to the highest possible value in the preferred window allows us to keep the number of disfavored models to a minimum. Without further information on the values of f a and individual m Qi , excluding fewer models may be advantageous in the sense of presenting a more inclusive E/N catalog.

E. Other interesting model properties
Let us summarize a few other possible model properties, already discussed in Refs. [13,14].
Since the axion is the angular degree of freedom of the PQ scalar field, it has a periodic potential and several degenerate vacua, given by the domain wall number N DW = 2N . During PQ symmetry breaking, the axion field can settle into any of these degenerate minima in different Hubble patches, giving rise to domain walls. The energy density contained in such topological defects can far exceed the energy density of the Universe [29] in the post-inflationary PQ breaking scenario. However, in models with N DW = 1, the string-domain wall configuration would be unstable [30], which presents a possible solution and makes N DW = 1 a desirable property of such models.
However, the DW problem can be avoided by allowing for a soft breaking of the PQ symmetry [29]. Moreover, in a pre-inflationary PQ symmetry breaking scenario, the patches and the topological defects are inflated away [31]. In line with Refs. [13,14], we therefore do not impose this criterion.
Among the 15 LP-allowed representations, only two have N DW = 1. When all Q i have the same U(1) PQ charges, such a restriction would forbid any models with multiple heavy quarks. With this in mind, a constraint on N DW is not used to exclude N Q > 1 models. In cases where the Q i are permitted to have opposite U(1) PQ charges, more complicated models with N DW = 1 can be built by choosing the Q i such that i N i = 1/2. Even then, the number of such models is few in comparison to the whole set of LP-allowed models.
Another intriguing property is the unification of the gauge couplings due to the presence of the Qs. The authors of Refs. [13,14] note that one of the 15 LP-allowed representations induces a significant improvement in unification. While we do not investigate this further, we expect to find more models that improve unification for higher N Q , which might be an interesting topic for a future study.

IV. MODEL CATALOG AND ANOMALY RATIO DISTRIBUTIONS
Let us discuss a few key findings and properties of the model catalog created in this work, which we summarize in Table I and Fig. 1.
To structure the discussion, we single out two subsets of the total model space: one where all Q i transform under the same representation and one where the representations are arbitrary but the U(1) PQ charges of the quarks have the same sign (we call these "additive models").

A. Subset I. Identical representations
First, consider the case where only representations of the form N Q j=1 r i with fixed i ∈ [1,20] are allowed. The number of possible models for a given N Q is then simply N r = 20, such that the total number of models up to and including some N Q is N tot = N r N Q .
Given that all quarks in such models have the same representation and U(1) PQ charge, only twelve discrete values of E/N are allowed when the LP criterion is taken into account [13]. However, the relative distribution is determined by the effect of each representation on the gauge group beta functions. We find that only LP-allowed model for N Q = 28 and that there are in total 79 preferred models in this subset.

B. Subset II. Allowing different additive representations
Next, consider the case where we can have arbitrary additive representations, written in such a way that they respect the relabeling symmetry: We find that, after applying the selection criteria, there are 59,066 preferred models for N Q ≤ 28. In particular, for N Q = 28, there are only nine LP-allowed models, none of which can be extended by another quark while preserving the criterion. The highest freedom in this subset is found for N Q = 10, where 5,481 models fall in the preferred region.
Among these models, the smallest and largest anomaly ratios are 1/6 and 44/3 respectively, both of which come from N Q = 1 models. The median of the distribution of this set of models is med(E/N ) ≈ 1.87, indicating TABLE I. Selected statistics for the complete set of models with NQ ≤ 9. We include information about the E/N ratios that give rise to the largest axion-photon coupling i.e. E/N ≡ argmax E/N (|E/N − 1.92|), photophobic models (|E/N − 1.92| < 0.04), and preferred (LP-allowed and N = 0) models.

NQ
Total #models E/N LP-allowed aγγ . We define models as "photophobic" if their E/N ratio is within one standard deviation of the nominal C (0) aγγ value: We find that 3,255 models (≈ 5.5%) among the 59,066 non-equivalent models are photophobic. Considering all preferred additive models up to N Q ≤ 28, there are 443 different E/N values. Out of these, 28 are unique in the sense that they are uniquely identifiable since their anomaly ratio E/N is different from any other nonequivalent model.

C. Complete set
Finally, let us comment on the complete set of possible models where we may also subtract representations, denoted by " ." Allowing U(1) PQ charges to have one of the two possible values for each Q i , we open the window to a much wider range of possible E/N values. In particular, the anomaly ratio, and thus the axion-photon coupling, can become negative (see Fig. 2) and, as mentioned before, the solution to the strong CP problem can be spoilt in models with N = 0.
For n ⊕ + n = N Q , where n ⊕ and n are the number of Qs with "positive" and "negative" U(1) PQ charges, 3 respectively, the number of models with n ⊕ > n is simply In the case where n ≡ n ⊕ = n , accounting for the fact that the anomaly ratio depends on the relative U(1) PQ charges of the Qs such that we have an equivalence of the type (r i ⊕ r j ) (r k ⊕ r l ) ∼ (r k ⊕ r l ) (r i ⊕ r j ), we also need to take care not to double-count models exhibiting this symmetry, giving With this, we find that the number of models grows very fast as N Q increases. This also makes it computationally difficult to compute and store all of the different combinations -let alone check the criteria for preferred models. We therefore restrict the complete analysis in this case to N Q ≤ 9.
The anomaly ratio distribution in the complete set exhibits a peak near zero, and we expect the trend to continue even for larger N Q . However, in general care should be taken when interpreting the "trends" visible in Fig. 1. For example, the number of LP-allowed models will eventually go down again as we move towards N Q = 28, despite the quickly growing total number of possible models. One may speculate that the number of uniquely identifiable E/N ratios could exhibit a similar behavior as the number of LP-allowed models, while the number of different E/N might eventually saturate.
Allowing for opposite U(1) PQ charges gives rise to models with large axion-photon coupling; the largest and smallest values of E/N found, 170/3 and −166/3 respectively, give larger |C aγγ | than what is possible in the previously discussed subsets. Note that the N Q = 8 model for E/N = −166/3 (r 2 ⊕ r 2 ⊕ r 5 ⊕ r 6 ⊕ r 7 r 1 r 9 r 9 ) was not reported in Refs. [13,14] as giving the highest possible |C aγγ |; instead the authors indicated that E/N = 170/3 led to the largest absolute value of the coupling. We find that among the complete set of 5,753,012 preferred models, there are 81,502 photophobic models and 820 different anomaly ratios, with 79 out of those also being from uniquely identifiable models.

V. IMPACT ON AXION SEARCHES
In this section, we discuss possible statistical interpretations of the hadronic axion model catalog and show the impact of these on the mass-coupling parameter space.

A. On constructing E/N prior distributions
The catalog of KSVZ models -even after applying the selection criteria -is but a list of possible models. It does not inherently contain information about how probable each model is. The model with E/N = −166/3 gives the largest |C aγγ | ≈ 57, which will place an upper bound on the axion-photon coupling and delimit the upper end of the KSVZ axion band. On the other end, complete decoupling with photons (C aγγ ≈ 0) is also possible within the theoretical errors. Since any of the models might be realized in Nature, perhaps due to a deeper underlying reason that is not obvious at present, one might be satisfied with this picture.
However, the boundaries of the band are extreme cases and do not take into account where the bulk of possible models can be found. For example, defining a desired target sensitivity for an experiment becomes non-trivial in the face of C aγγ potentially being extremely close to zero. We propose instead that covering a certain fraction of all possible models or constructing a prior volume might be more meaningful ways to define such a target.
To directly interpret an E/N histogram as a distribution implicitly makes the assumption that each model is equally likely to be realized in Nature. While this interpretation might be considered "fair," one could argue that models with many Qs are more "contrived" and consequently introduce a weighting factor that penalizes models with N Q 1. This could be achieved with e.g. exponential suppression via a weighting factor ∝ e −N Q , or ∝ 2 −N Q . Another option could be to choose models that are minimal extensions (N Q = 1) or similar to the family structure of the SM (N Q = 3 or e.g. a weighting Such consideration are aligned with the Bayesian interpretation of statistics, and will probably meet criticism for this reason. However, as pointed out in Ref. [20], at least in the pre-inflationary PQ symmetry breaking scenario, which is fundamentally probabilistic in nature, the Bayesian approach is well motivated. Furthermore, Ref. [20] also proposed that the discrete nature of KSVZ models should be reflected in the prior choice of E/N . Such a physically-motivated prior should further reflect the combinatorics of KSVZ model building by including the multiplicity of E/N ratios. As mentioned at the end of Section II, this multiplicity also depends on whether or not the Q i are distinguishable by e.g. having different masses.
With this in mind, we show different statistical interpretations of the anomaly ratio in Fig. 3. For visualization purposes, we show kernel density estimates of the distributions for different weighting factors mentioned above, while reminding the reader that the underlying histograms and distributions are actually discrete and not continuous. From Fig. 3 it becomes clear that the different weightings can change the width of the distribution, introducing a prior dependence in an analysis. However, the modes of the distributions remain around E/N ∼ 2, which means that a partial cancellation of the axion-photon coupling C aγγ is typically possible, as already observed in Fig. 2.

B. Experimental constraints on preferred KSVZ axion models
Of course, a possible partial cancellation of the axionphoton coupling has consequences on the various astrophysical, cosmological, and laboratory searches (see e.g. Ref. [10]) for axions. The most powerful analyses combine the results of different experiments to place joint limits on the properties of different types of axions (e.g. Refs. [20,35,36]).
To investigate this further, consider e.g. a prior on E/N where all preferred (LP-allowed models with d ≤ 5 operators and N = 0), non-equivalent KSVZ models are considered equally probable. 4 We can then generate samples for C aγγ = E/N − C  Fig. 4, and the bulk of these models can be constrained by present and future experiments. In fact, while complete cancellation of C aγγ is possible within the theoretical uncertainty for some E/N values, we find that the bulk of models is at worst somewhat suppressed. This is very encouraging for experimental searches.
Had we only considered additive models, the 95% region would be |C aγγ | ∈ [0.02, 1.67], such that the upper end of the band would be lower than the traditional KSVZ model with E/N = 0. This can be readily understood from the E/N distributions in Fig. 3, whose mode is typically close to C (0) aγγ such that the value of |g aγγ | is lower than what would be expected for |C aγγ | ∼ O(1). In this case, the planned future experiments would not be able to probe large parts of the band, indicating that the choice of prior -even if physically-motivated -can induce a noticeable impact on the results.

VI. SUMMARY AND CONCLUSIONS
We provide a catalog of all hadronic, or KSVZ, axion models with N Q ≤ 9, featuring 1,027,233,129 nonequivalent models in total. When we apply the selection criteria for preferred models, we find a limit of N Q ≤ 28 and that only 5,753,012 non-equivalent models with 820 different E/N values exist (59,066 non-equivalent models with 443 different E/N values for additive representations). While relaxing existing or adding new criteria can increase or reduce these numbers, we generically expect that the Landau pole (LP) criterion will be a powerful tool to limit the number of possible models -even with modified constraints or in other axion models. This is similar to the Standard Model, where the number of families can also be restricted by demanding that LPs do not appear below some energy scale. We further propose that only models with QCD anomaly N = 0 be considered preferred.
Our model catalog can be a useful, searchable database for researchers wishing to study the KSVZ axion model space. It allows to e.g. make statements about what fraction of possible models a given experiment is sensitive to. We made catalogs, histograms, and example Python scripts available on the Zenodo platform for this purpose [15].
Some models in the catalog might be considered "contrived" as they add many new particles to the theory. Of course, in case of a discovery or if any other appealing reason for a seemingly more complicated models is put forward, this perception might change. In absence of such reasons, the E/N values may be interpreted as statistical distributions, which encode assumptions about the probability of the different models. We generally outlined how prior distributions can be constructed from the catalog and gave concrete examples of such choices.
For the specific choice of equally probable preferred models, we consider the consequences for axion searches and the definition of the KSVZ axion band. Here we suggest that the latter may be defined as the central 95% region of all models, taking into account uncertainties from the model-independent contribution to the axionphoton coupling. If only "additive models" are considered, the bulk of the preferred models can unfortunately not be probed by current or future experiments since the anomaly ratio distributions in this case tend to peak around E/N ∼ 2. In general, using the discrete E/N distributions improves on unphysical prior choices considered in the past (e.g. Ref. [20]).
Even when ignoring the statistical perspective, it is useful for axion searches to know that the preferred models only admit 820 different E/N values. In case of an axion detection, one may therefore test these discrete models against each other to see which models are most compatible with the detected signal. One could further test them against a generic axion-like particle or other QCD axion models. In an ideal scenario, this might even allow an experiment to infer the underlying high-energy structure of a model, which highlights the known property of axion models to connect high-energy physics to low-energy observables.
In summary, the powerful LP criterion restricts the number of KSVZ models to a finite value. In that sense, the catalog presented here is a complete list of all preferred KSVZ models, which may be used as input for axion searches and forecasts. Since KSVZ models could e.g. be extended by also considering multiple complex scalar fields or feature more complex couplings to the SM, and . For context, we show various present (shaded regions) and future (dashed lines) haloscope (blue) and helioscope (red) limits and forecasts [32] as well as bounds from hot dark matter [33], energy loss in SN1987A [34], and recent string simulations [22].
since there are other kinds of QCD axion models such as the DFSZ-type models, this work presents another step forward in mapping the landscape of all phenomenologically interesting axion models.