High-quality axions in solutions to the $\mu$ problem

Solutions to the $\mu$ problem in supersymmetry based on the Kim-Nilles mechanism naturally feature a Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) axion with decay constant of order the geometric mean of the Planck and TeV scales, consistent with astrophysical limits. We investigate minimal models of this type with two gauge-singlet fields that break a Peccei-Quinn symmetry, and extensions with extra vectorlike quark and lepton supermultiplets consistent with gauge coupling unification. We show that there are many anomaly-free discrete symmetries, depending on the vectorlike matter content, that protect the Peccei-Quinn symmetry to sufficiently high order to solve the strong CP problem. We study the axion couplings in this class of models. Models of this type that are automatically free of the domain wall problem require at least one pair of strongly interacting vectorlike multiplets with mass at the intermediate scale, and predict axion couplings that are greatly enhanced compared to the minimal supersymmetric DFSZ models, putting them within reach of proposed axion searches.


I. INTRODUCTION
Extensions of the Standard Model (SM) of particle physics are plagued by several apparent hierarchy problems, which can be viewed as hints towards the ultimate completion of the theory. The vacuum energy, as manifested in the observed expansion rate of the universe due to the cosmological constant, is approximately 120 orders of magnitude smaller than its naive dimensional analysis estimate M 4 P , where M P = 2.4 × 10 18 GeV is the reduced Planck mass scale. The "big hierarchy problem" is that the electroweak scale set by the Higgs field squared mass parameter is 32 orders of magnitude smaller than M 2 P . The strong CP problem is that the CP-odd angle θ in the QCD Lagrangian could be of order 1, but is constrained to be smaller than 9 × 10 −11 by the measured value [1] of the electric dipole moment of the neutron.
In this paper, we will be concerned with a class of models that simultaneously address the latter two problems. The big hierarchy problem of the Higgs field is addressed by supersymmetry, with supersymmetry breaking terms characterized by the TeV scale, † while the strong CP problem is addressed by including a global Peccei-Quinn (PQ) U (1) symmetry [2,3] that is explicitly broken by a QCD anomaly but also spontaneously broken, minimizing the effective value of |θ| and giving rise to a light, very weakly coupled axion [4]- [12]. For reviews of the axion solution to the strong CP problem from various points of view, see refs. [13]- [22].
The "µ problem" of supersymmetry relates the issues of supersymmetry breaking and the PQ symmetry breaking. The superpotential of the Minimal Supersymmetric Standard Model (MSSM, for a review see [23]) contains terms of the form (1.1) where H u and H d are the Higgs doublet chiral superfields. The q, u, d, , and e are quark and lepton chiral superfields, for which we suppress flavor and gauge indices and use lowercase letters to distinguish them from additional vectorlike quark and lepton superfields to be considered below. Besides the Yukawa coupling matrices y u , y d , and y e , this superpotential contains a single dimensionful parameter µ, which preserves supersymmetry but nevertheless should be roughly of the same order as the supersymmetry breaking mass scale in order to allow for electroweak symmetry breaking. To accomplish this, the Kim-Nilles mechanism [24] takes the µ term to be absent from eq. (1.1) in the ultraviolet theory, and in its place introduces non-renormalizable terms, for example of the form proposed in ref. [25]: 2) † In this paper, we will loosely refer to the mass scales associated with superpartner masses and supersymmetry breaking as the TeV scale. Given the negative search results so far at the Large Hadron Collider (LHC), the masses of the strongly interacting superpartners are evidently somewhat larger than 1 TeV. The fact that their masses exceed the Higgs vacuum expectation value by an order of magnitude is the "little hierarchy problem", which we do not address in this paper; it is obviously much less severe than the big hierarchy problem, for which the LHC results have not favored any competing hypotheses.
where X, Y are gauge-singlet chiral superfields, and λ µ and λ are dimensionless couplings. The term λ in eq. (1.2) is included to stabilize the potential at large X, Y scalar field strengths. Including the effects of supersymmetry breaking, there are also contributions to the Lagrangian (in terms of scalar fields, for which we use the same symbol as the corresponding superfield): where a µ , a are of order m soft and m 2 X , m 2 Y are of order m 2 soft , with m soft roughly at the TeV scale. If m 2 X and m 2 Y are negative, or just sufficiently small compared to |a/λ| 2 , then the resulting scalar potential has a local minimum for vacuum expectation values (VEVs) parametrically of order ‡ (1. 4) a scale intermediate between the Planck and TeV scales. This can always occur if m 2 X and m 2 Y are negative, but this is not a necessary condition. For example, if m 2 X = m 2 Y = m 2 , then there is a non-trivial local minimum if |a/λ| 2 − 12m 2 > 0, and it will be a global minimum if |a/λ| 2 − 16m 2 > 0. Symmetry breaking along such flat directions is stable against radiative corrections [26]. In the low energy theory, the µ and b terms of the MSSM Lagrangian are (1. 6) solving the µ problem.
If m soft is of order the TeV scale, then 10 9 GeV ∼ < M int ∼ < 10 12 GeV, (1.7) depending on the dimensionless parameters involved. Equation (1.7) roughly coincides with the preferred range for the decay constant of the axion from astrophysical constraints on the low end and from cosmological dark matter density on the upper end. It is therefore notable that X and Y carry non-zero charges for a PQ symmetry, which the VEVs spontaneously break, so that the axion could be a linear combination of the pseudo-scalar components of the fields X and Y (with a very small admixture of H u and H d ). Since the MSSM quarks also carry non-zero PQ charge, this is an example of a Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) [8,9] axion model. In this way, the solution of the µ problem in supersymmetry can also be the solution of the strong CP problem. ‡ Note that it is important to include the contribution of the allowed holomorphic supersymmetry-breaking coupling a here. Omitting it would seemingly lead to Y X ∼ M int . Base model Superpotential terms PQ charges of (X, Y ) Besides the axion, the other components of the supermultiplets X and Y all get TeV scale masses. These include a scalar saxion, and gauge-singlet axino fermions, one of which could be the lightest supersymmetric particle. The lightest MSSM superpartner could decay to the axino with a macroscopic proper decay length, which can be much larger than the size of a collider detector. However some of the decays can occur within the detector, providing a striking search signal [27] for the Large Hadron Collider.
In addition to eq. (1.2), there are three other similar but distinct superpotential structures involving two PQ-breaking superfields X and Y . (For simplicity and economy, we restrict our attention to only two such fields. Although it is possible to have more than two, this would seem to make it harder to find solutions to the axion quality problem to be discussed shortly.) They are: (1. 10) with corresponding holomorphic soft supersymmetry-breaking parameters a µ and a in each case, in the obvious way. The structure in eq. (1.8) was proposed in ref. [28], and those in eqs. (1.9) and (1.10) in ref. [27]. For a review and further elucidation, see [29]. In the following, we will refer to the models defined by eqs. (1.2), (1.8), (1.9), and (1.10) as base models B I , B II , B III , and B IV , respectively, since we will be interested in extensions of them. Each of them implies a different assignment of PQ charges, which are summarized in Table  1.1. We choose the normalization of the PQ charges so that H u H d has charge −2.
The axion quality problem [30]- [34] results from the possible presence of higher dimensional contributions to the Lagrangian that explicitly violate the U (1) PQ symmetry, since these can displace the QCD θ parameter away from 0 at the minimum of the scalar potential, spoiling the solution to the strong CP problem. Such contributions are expected to be allowed if they are not forbidden, because ungauged symmetries are not respected by quantum gravitational effects, at least in thought experiments such as black hole evaporation processes. In supersymmetry, we therefore consider superpotential operators of the form (1. 11) with dimensionless κ. Together with the λ term in the superpotential, this gives rise to scalar potential terms (from |F X | 2 and |F Y | 2 ) that have p+2 powers of X and Y and are suppressed by 1/M p−2 P . These are also accompanied by holomorphic supersymmetry-breaking terms in the scalar potential where a κ should be of order the TeV scale. For generic phases of κ and a κ , both of these types of contributions result in tadpoles for the axion field A (a linear combination of the imaginary parts of X and Y expanded around the CP-conserving vacuum) that are parametrically of the same form where we have identified the axion decay constant f A with the intermediate scale VEVs M int , and the dimensionless quantity δ depends on the parameters (including the magnitudes and phases of λ, κ, and a κ and the integer k) in a complicated way. The axion potential also includes the usual squared mass term, approximately for small A, The combination of eqs. (1.13) and (1.14) gives rise to (1. 16) Now by requiring that the θ eff parameter at the minimum of the potential is less than 10 −10 from the experimental bound on the neutron electric dipole moment, one therefore finds [32][33][34] that one should have p + 2 > 88 + log 10 (δ) 9.4 − log 10 (f A /10 9 GeV) . (1.17) Note that the bound on p becomes weaker for smaller f A . With the naive δ ≈ 1, if PQviolating superpotential terms with p = (8, 9, 10, 11, or 12) are present one should have f A ∼ < (4×10 9 , 3×10 10 , 10 11 , 4×10 11 , or 10 12 ) GeV, respectively. However, like all naturalness criteria, this one is inherently fuzzy, as δ could be significantly less than 1, for example because the corresponding coupling(s) happen to have a small magnitude and/or a phase alignment with the bare θ. We will therefore not commit to a specific requirement for p, with the understanding that larger p is safer in some sense. Previous works exploring solutions to the axion quality problem have invoked composite axions models [35]- [45], additional continuous gauge symmetries [32,34], [46]- [53], and discrete gauged symmetries in non-supersymmetric [54]- [58] and supersymmetric [59]- [66] models. Here, we will be interested in the latter type of idea, in which the U (1) PQ symmetry arises as an approximate accidental consequence of a discrete symmetry imposed on superfields. We will consider models in which discrete symmetries forbid PQ-violating Lagrangian terms up to some mass dimension, and refer to the smallest exponent allowed in a particular model for the PQ-violating terms of the types in eqs. (1.11) and (1.12) as the PQ-violation suppression p. Note that if the suppression p for a particular model is odd, then it can be bumped up to the next integer by simply imposing another Z 2 symmetry under which both X and Y are odd, since this forbids all terms of the forms eq. (1.11) and (1.12) with odd p. We therefore consider as potentially viable any models whose discrete symmetries predict that p should be 7 or more, and consider models with p ≥ 12 for all PQ-violating terms involving only X and Y as presumptively high-quality.
Since the Kim-Nilles mechanism provides a µ term for the H u and H d fields at the TeV scale, it is reasonable to suppose that the same mechanism can give masses to other vectorlike pairs of chiral superfields as well, some of which could therefore be at the TeV scale just like the MSSM Higgs and higgsino particles. For each additional pair Φ + Φ of chiral superfields, supersymmetric mass terms near the TeV scale can arise in three possible ways, due to non-renormalizable superpotential terms of the forms: assuming λ Φ is of order one. Alternatively, masses at the intermediate scale can be achieved by renormalizable superpotential terms of the forms: (1.19) for vector-like chiral superfield pairs, in the extensions of the four base models, depending on the possible superpotential terms that provide masses due to the scalar components of X and Y obtaining VEVs of order M int . The first three mass terms provide for a TeV scale mass, and the last two provide for an intermediate scale mass.

Mass terms B
(1, 2, −1/2) + (1, 2, 1/2) D 6 + D 6 (6, 1, 1/3) + (6, 1, −1/3) The net PQ charges of ΦΦ are important for understanding the low-energy axion couplings, as discussed in the next section. Table 1.2 shows the possible values of PQ charges Q ΦΦ ≡ (Q Φ + Q Φ ) in the extensions of the four base models, for the various possible superpotential mass terms. We will consider pairs of additional vectorlike quark or lepton superfields Φ + Φ chosen from among those in Table 1.3. The first 5 pairs each include one field with the same color and electroweak quantum numbers as the MSSM chiral superfields, from which they are distinguished by the use of capital letters. While any of the possible mass terms could be used independently for each of types of fields, it is well-known that nearly degenerate sets of chiral superfields in 5 + 5 = D + D + L + L or 10 + 10 = Q + Q + U + U + E + E representations of the SU (5) grand unified theory [67] will preserve the apparent unification of gauge couplings observed in the MSSM, as illustrated for example in the left panel of Figure 1.1. Note that we do not assume that SU (5) is actually the unbroken gauge group in the ultraviolet, and we allow for different components of the 5 + 5 and/or 10 + 10 to have different mass source terms and therefore different PQ charges. In Table 1.3 we have also included a more exotic pair, an SU (3) c color sextet (quix) and its conjugate with electric charges ±1/3, denoted D 6 + D 6 . If these have masses at an intermediate scale M int ≈ 10 11 GeV and there are two L + L pairs near the TeV scale, then one can again have gauge coupling unification at a scale somewhat higher than in the MSSM, as shown in the right panel of Figure 1.1. Although this combination might seem somewhat of an ad hoc choice, it is of interest because it provides an example of an even more enhanced coupling of the axion to photons, as we will see. There are other classes of quixotic models that can have gauge coupling unification, for example D 6 + D 6 or U 6 + U 6 = (6, 1, −2/3) + (6, 1, 2/3) quixes at the intermediate scale that have one or two weak isotriplets T ∼ (1, 3, 0) at the same scale. Such models do not predict any new particles at the TeV scale, and therefore are omitted from our discussion.
The rest of this paper is organized as follows. In section II, we review the general properties of axion models, in particular the relation between PQ charge assignments and axion couplings, the domain wall problem, and the axion quality problem. The relevant relations are obtained for supersymmetric DFSZ axion models of the type described above. In section III, we explore possible discrete symmetries that can solve the axion quality problem, and establish that there are many available for all of the base models and their extensions. In section IV, we discuss axion signals and detection prospects for our models, including both present constraints and future detection prospects. Section V contains some summarizing remarks.

A. Axion properties in general
In this section we review the general properties of QCD axion models, following the discussion in [20]. The anomalous divergence of the Peccei-Quinn current in a model with each left-handed Weyl fermion ψ f transforming under SU (3) c in the representation R f with a PQ charge Q f (such that ψ f → e iQ f α ψ f ) is given by: respectively. Here, the trace in N is over all left-handed fermion representations with SU (3) c index T (R f ), while the trace in E is over all left-handed fermion fields with electromagnetic (EM) charge q f . The PQ current can be explicitly written as: where ϕ s are complex scalar fields, each with a PQ charge Q s (such that ϕ s → e iQsα ϕ s ). The PQ symmetry can be spontaneously broken by non-zero vacuum expectation values (VEVs) of some scalars ϕ s ≡ v s (without breaking color or EM) that can be written as: Here, ρ s and a s are the canonically normalized radial and angular fields, respectively. In our case, the scalars with VEVs will be either doublets or singlets under SU (2) L . The PQ current becomes: in the decoupling limit, i.e. ignoring the terms with heavy radial fields ρ s , which are PQ invariant. Now, defining the pseudo-Nambu-Goldstone boson axion field A is given by: This results in We also note that under a U (1) PQ transformation parameterized by α, a s → a s + αQ s v s , (2.10) Assuming the PQ symmetry is not explicitly broken, the PQ charges in a theory are determined by the terms in the Lagrangian. However, there is an ambiguity in the charges due to the freedom to add multiples of the weak hypercharge and any other non-anomalous U (1) symmetries that might be present in the theory. The "physical" PQ charges can be obtained by imposing an orthogonality condition: 12) to ensure that the axion and the Z boson do not mix. Here, Y s are the U (1) Y weak hypercharges of the scalars that get VEVs. The above condition can be obtained by requiring the axion A to be invariant under U (1) Y gauge transformations parameterized by α Y : (2.14) The would-be Nambu-Goldstone boson G that becomes the longitudinal component of the Z-boson is: The anomalous divergence of the PQ current in eq. (2.1), with j µ PQ in eq. (2.9), now implies: This equation can be obtained as the Euler-Lagrange equation of the following effective Lagrangian: where, provided that the PQ symmetry has a non-zero QCD anomaly N , the axion decay constant is defined by As a convention, we choose to define f A as always positive. In some models N is negative, in which case v A is negative in our convention.
In terms of four-component Dirac fermions Ψ f with left-handed and conjugate of righthanded components carrying PQ charges Q f and Q f , respectively, One can now express the axion-fermion couplings in terms of axial currents only, eliminating the vector currents. This can be achieved by redefining Dirac fermions Ψ f by the following U (1) vector transformation (as opposed to the axial transformation that contributes to the PQ anomaly): under which the kinetic terms of the fermions give an additional contribution that cancels the axion-fermion couplings involving the vector currents, After evaluating the coefficients of axion-photon and axion-fermion couplings given in eqs. (2.20) and (2.21), respectively, the resulting low-energy axion interaction Lagrangian can be obtained as where the axion mass is given in terms of the axion decay constant f A by [68] m A = 5.691(51) 10 12 GeV f A µeV, (2.26) and [69] C Aγ = c γ − 1.92(4), (2.27) The tree-level couplings to fermions above do not include the effects of renormalization group running [70]- [75], which are usually small compared to experimental uncertainties. A notable exception is the coupling C Ae when c e vanishes or is very small in hadronic axion models, which we will include as non-supersymmetric benchmarks below. In those cases we include the leading logarithmic contribution [70,71] (2.32) Finally, the effective Lagrangian which provides the axion coupling to photons and on-shell electrons, protons, and neutrons can be written as with the definitions Here we note that the contribution from the CP-violating terms where axion derivatively couples to vector fermion current in eq. (2.19) vanish after integrating by parts and using the Dirac equation (without having to redefine the fermion fields as done in eq. (2.23)). After the PQ symmetry breaking, a discrete subgroup Z N DW (= e 2kπi/N DW , k = 0, 1, . . . , N DW −1) is left unbroken. N DW is the Domain Wall (DW) number that corresponds to the number of discrete set of inequivalent degenerate minima of the axion potential [76]. With the above definition of PQ-QCD-QCD anomaly coefficient N , the DW number can be computed as [77]: where n s ∈ Z. The domain wall number must be invariant under the rescalings of the PQ charges, as is reflected in the above formula. Models with N DW > 1 may have non-trivial cosmological implications. In particular, formation of topological defects such as stable domain walls, due to degenerate vacua with different possible phases of the axion, can dominate the early universe [76]. However, the cosmological domain wall problem may not arise if PQ symmetry is broken in the preinflationary era and is not restored after inflation ends, so that the observable universe today consists of a single patch that initially had a common value of θ. Another possibility [78] is that the apparent Z N DW discrete symmetry relating different vacua are actually embedded within a continuous gauged symmetry that is spontaneously broken at high energies. The N DW apparently distinct vacua are then actually the same, being connected by the continuous gauge symmetry. In this paper we will concentrate instead on the possibility that safety is achieved by N DW = 1, in which case any domain walls that form are bounded by strings on their edges, and are highly unstable [79]. This can be achieved using a strategy similar to that in [80], by introducing heavy fermions charged under the PQ symmetry, in our case components of chiral supermultiplets that are vectorlike under the Standard Model gauge group.

B. Application to DFSZ axions in supersymmetry
In the following, we consider supersymmetric DFSZ-type axion models with two PQbreaking gauge-singlet fields X, Y as discussed in the Introduction, consisting of the four base models B I , B II , B III , and B IV summarized in Table 1.1, with possible extensions by additional vectorlike superfields that maintain approximate gauge coupling unification. The additional chiral superfield content in the models considered below are listed in Table 1 From our choice of normalization in Table 1.1, we get the following condition on the PQ charges of the MSSM Higgses: (2.36) Requiring that the Yukawa terms in the MSSM superpotential are U (1) PQ invariant, we can express the following PQ charge combinations−the only combinations that enter the low-energy axion-fermion couplings−in terms of Q Hu : We have used eq. (2.36) in the last two equations. Using the above constraints, the PQ charges of all the MSSM fields can be expressed in terms of only three numbers Q Hu , Q q , and Q , which at this point are still unconstrained and therefore can be treated as free parameters. This corresponds to the freedom to add multiples of other U (1) charges to the PQ symmetry charges. In order to obtain the low-energy axion couplings and the DW number, the scalars (represented by the same symbols as their respective supermultiplets) that acquire a VEV are parameterized as: where H 0 u and H 0 d are the neutral MSSM Higgs scalars, a s = {a x , a y , a u , a d } are the pseudoscalar bosons that contribute to the axion, and v x , v y v u , v d in an invisible axion model. Imposing the orthogonality condition eq. (2.12) in each of the base models yields: The orthogonality condition has essentially used the freedom to add arbitrary multiples of the U (1) Y charges to the fields. The PQ charges for the MSSM Higgs supermultiplets using eq. (2.36), consistent with the Higgs VEVs and our normalization convention Q Hu + Q H d = −2, are: and we have, from eq. (2.7), 43) with N to be given soon below. The contribution 4s 2 β c 2 β v 2 is numerically negligible. The superfield charges of the Peccei-Quinn symmetry and two linearly independent anomaly-free U (1) symmetries of the dimensionless part of the MSSM superpotential.
Here, ν represents either a gauge-singlet neutrino or the combination −H u appearing in the non-renormalizable Weinberg operator for neutrino masses. The physical Peccei-Quinn charges given here are obtained by imposing the orthogonality condition eq. (2.12) and the normalization condition Q Hu + Q H d = −2, and are given in terms of two free parameters Q q , Q , and the ratio of Higgs expectation values tan β = s β /c β . If a seesaw neutrino mass term S ν ν is included in the superpotential for some gauge-singlet field S whose scalar component gets a VEV, then 2T 3 Using eqs. (2.37)-(2.39) the PQ charges of the rest of the MSSM chiral superfields can be fixed in terms of Q q and Q . We can also require that neutrino masses are provided by the superpotential version of the Weinberg operator (H u )(H u ). Equivalently from the point of view of anomalies, we can introduce a gauge-singlet neutrino superfield ν which has superpotential terms for example either a bare mass term or X or Y . This then gives us the further constraint where Q S = 0 in the case of a bare seesaw mass (λ ν S → M ν ). The resulting PQ charge assignments for the MSSM chiral superfields are shown in Table 2.1. These physical PQ charge assignments also hold for any extensions of the base models, as long as there are no additional scalars with non-zero VEVs that can potentially feed into the orthogonality condition in eq. (2.12). Also included for future reference are the charges for two non-anomalous family-independent U (1) symmetries; one is 6Y , the weak hypercharge normalized to integer charges for all fields, and the second is 2T 3 R = 2Y −(B−L), also with integer charges.
Using eqs. (2.2) and (2.3), the PQ anomaly coefficients in the base models and their extensions to include vectorlike fields are: where the contributions Q ΦΦ can be read from Note that all of N , E, c γ , c u , c d , and c e , and therefore the low-energy effective Lagrangian parameters eqs. (2.27)-(2.31) and (2.34), are independent of the choices of Q q and Q , and we arrive at for the central value with parenthetical uncertainty estimates in our base models and extensions, parameterized only by N and E. Note that the model-independent (N -independent) part of the neutron coupling g An suffers from an accidental cancellation [69]. We now compute the domain wall number (i.e. the number of independent degenerate minima of the axion potential) for each of the four base models listed in Table 1.1 using the formula given in eq. (2.35).
where n x , n y , n u , and n d are integers accompanying the scalars X, Y , H u , and H d , respectively. In the above expression, our goal is to arrange each n s such that eq. (2.53) is an integer. This is obtained provided that n y = −3n x and n u +n d = 2n x . Then from eqs. (2.35) and (2.53) we find: where n x is a non-zero integer. Therefore, since N = 3 for the MSSM field content, the domain wall number of the base model B I is N DW = 6.
Model B II : Following a similar procedure as above, we have: which is again an integer if n y = −3n x but this time with n u + n d = −2n x . The domain wall number is therefore given by: where again n x is a non-zero integer. For the base model B II with N = 3, this again amounts to N DW = 6.
Model B III : Proceeding as above, we begin with: (2.57) This will be an integer if n y = −3n x and n u + n d = 6n x , resulting in for n x a non-zero integer. For the base model B III with N = 3, we obtain N DW = 18 .

Model B IV :
To calculate the domain wall number for this model, we start with: This will be an integer if n y = −n x and n u + n d = −2n x , with the result where again n x is a non-zero integer. This gives N DW = 6 for the base model B IV .
The formulas for N DW in eqs. (2.54), (2.56), (2.58), and (2.60) in terms of N are also applicable to any extensions of the models B I , B II , B III , and B IV , respectively, provided that there are no additional scalars getting VEVs. In each case, the value of the anomaly N should be computed by including contributions from the other strongly interacting chiral superfields present in the model extension, as in eqs. (2.46). Then N DW is the smallest integer |2N n x | (for extensions of base models B I , B II , and B IV ), or the smallest integer |6N n x | (for extensions of base model B III ), where n x is a non-zero integer.
Although N DW = 1 in the base models and many of their extensions, it follows from the preceding that it is possible to achieve N DW = 1 models and avoid the cosmological domain wall problem in certain extensions. This requires: It should also be noted that the new heavy vectorlike particles in the extended models must be allowed to decay to Standard Model particles in order to avoid dangerous cosmological charged relics. This is easy to arrange, as the decays can be mediated by couplings to the Standard Model quark and lepton superfields via Yukawa terms in the superpotential, in several different ways. For example, in models with an extra 5 + 5, the vectorlike weak isosinglet down-type quark can decay through any of the dimensionless couplings H d qD or qD or udD, while the vectorlike weak isodoublet lepton can decay through any of the superpotential couplings H d eL or eL or qdL. It is easy to construct similar couplings that allow Q + Q, U + U , and E + E to decay in models that have an extra 10 + 10. In each case, in order to not overconstrain the system of equations determining the PQ charges, or the discrete symmetry charges discussed in the next section, one may select only one of the couplings for a given Φ + Φ pair. (This also avoids tree-level violation of baryon number and/or lepton number, which could lead to proton decay.) In models with quixes (the last two rows of Table 2.2), they can decay via a superpotential term of the form D 6 ud. This fixes the PQ charge to be Q D 6 = 2Q q − 2. The important point for the low-energy axion phenomenology is that the existence of all such Yukawa couplings maintains the freedom to choose the net PQ charges Q ΦΦ consistently with the vectorlike mass terms, and therefore does not affect the PQ anomalies.

Model extension
Base Mass terms

A. Conventions and assumptions for discrete symmetries
In this section we consider discrete symmetries that can protect the Peccei-Quinn U (1) symmetry up to dimension p in the X and Y superfields in the superpotential. In equivalent language, the Peccei-Quinn symmetry is an accidental consequence of imposing the discrete symmetry. We consider as possibilities an Abelian discrete symmetry Z n or a discrete Rsymmetry Z R n . In both cases, each chiral superfield Φ has integer charge z Φ = 0, 1, . . . n − 1 (mod n). It will be convenient to treat both cases in a unified framework. We therefore take the gauginos and the anti-commuting coordinates θ α to have charge r, and each superpotential term is required to have total Z R n charge 2r (mod n). For the case of ordinary (non-R) discrete symmetries we take r = 0. In both cases, all Lagrangian terms have total discrete charge 0. [For continuous R-symmetries it is customary to take r = 1 by a choice of normalization, but for discrete R-symmetries this is not always possible, since we require that all Z R n charges for all fields are integers (mod n).] In Table 3.1, we list the allowed † values of n that can provide for a maximum possible suppression p, for the base models and their extensions. This means that the lowest dimension superpotential term that is allowed by the discrete symmetry but violates the PQ symmetry has dimension p, and so has the form X j Y p−j /M p−3 P . From these results, we can conclude that n must be at least 13 (for models B I,II,III ) or 12 (for model B IV ) even if the axion decay constant f A is as low as 10 3.1: Allowed n, for discrete symmetries Z n and discrete R-symmetries Z R n , that can provide a maximum possible suppression up to p for the base models and their extensions. The meaning of p is that the lowest dimension PQ-violating superpotential term(s) allowed by the discrete symmetry and involving only X and Y are of the form X j Y p−j . No anomaly cancellation constraints are imposed, yet. † Note that n must be even for non-R (r = 0) symmetries Z n for model B IV and its extensions, in order to allow the term X 2 Y 2 but forbid the term XY . required p, the smallest possible values of the order n of the discrete symmetry group can be read from the table. Of course, not every discrete symmetry at order n will provide the suppression listed.
We now consider the conditions imposed by anomaly cancellation on the discrete symmetry group, which depend on the discrete symmetry charges of the superfields that are charged under the Standard Model gauge group. For simplicity, we assume that the MSSM quark and lepton chiral superfields have charges that are generation-independent. Using the symbol z Φ for the additive Z R n charge of the chiral superfield Φ, the Z R n × G × G anomalies for G = SU (3) c , SU (2) L , and U (1) Y are respectively:

2)
Here, we have adopted a convenient but somewhat unusual normalization, by taking the index for the fundamental representation of the non-Abelian groups to be 1 rather than 1/2, so that A 2 and A 3 are integers defined (mod n). Likewise, the anomaly A 1 has been given in a normalization so that it is always an integer (assuming the weak hypercharge quantization of the Standard Model), and again defined (mod n). The first term in each of A 3 and A 2 is the gaugino contribution. The number of chiral quark and lepton generations is n g = 3. Finally, the contributions ∆ ΦΦ are equal to the sum of the fermion discrete symmetry charges z Φ − r and z Φ − r for the indicated vectorlike chiral superfields, and are also equal to the negative of the sum of z X and z Y contributions for the terms in eqs. (1. 18) and (1.19) that produce their masses. With our normalization of the A i , the conditions [81][82][83][84], [61] for the discrete symmetry to be anomaly free can be obtained by treating Z R n as a subgroup of an anomaly-free continuous U (1) symmetry that is spontaneously broken by a scalar field VEV of charge n [85], in the presence of other possible very heavy fermions. The anomaly free condition is that one must have where ρ GS is a constant that can arise from the Green-Schwarz (GS) mechanism [86] that may be used to cancel the anomalies for the U (1) group that contains Z R n , and m 1,2,3 are integers, and k 2,3 are positive integer Kac-Moody levels, while k 1 in general could be arbitrary. For simplicity, we will therefore consider as a weak assumption that k 2 = k 3 , but also consider as a stronger assumption motivated by gauge coupling unification that k 3 = k 2 = k 1 = 1, corresponding to the normalization in which SU (3) c × SU (2) L × U (1) Y could be embedded in a simple GUT group, SU (5) or SO(10) or E 6 . We then have: for the weaker condition in which the U (1) Y normalization is considered arbitrary, with the additional stronger condition if the gauge couplings are required to unify in the usual way. The value ρ GS = 0 corresponds to the notable special case in which the Green-Schwarz mechanism does not play a role, in which case From the requirement that the MSSM superpotential couplings are allowed, we have z u = −z Hu − z q + 2r, and z d = −z H d − z q + 2r, and z e = −z H d − z + 2r. Furthermore, z H d is determined either by z H d = −z Hu − z X − z Y + 2r for base model B I and its extensions, or z H d = −z Hu − 2z X + 2r for base models B II and B IV and their extensions, or z H d = −z Hu − 2z Y + 2r for base model B III and its extensions. We also require that neutrino masses are provided by the superpotential version of the Weinberg operator (H u )(H u ). Equivalently from the point of view of anomalies and low-energy phenomenology, we can introduce three gauge-singlet neutrino superfields ν which have superpotential terms and then we have that z ν = 2r − z − z Hu must either be r, or perhaps r + n/2 if n is even. If we temporarily neglect the M ν term, the superpotential y u H u qu − y d H d qd − y e H d le + y ν H u ν is invariant under two anomaly-free U (1) symmetries with integer charges. One is U (1) 6Y where Y is the weak hypercharge, and the other can be taken to be either 3(B − L) or 2T 3 R = 2Y − (B − L). The charges were listed in Table 2.1. It is apparent that we can always redefine the Z R n charges by adding a multiple of 6Y in order to make q = 0, as a convention without loss of generality. Now, if n is even, we can add (n/2)2T 3 R to every z Φ to obtain another discrete symmetry which has the same anomalies A 1,2,3 , and differs from the original discrete symmetry only by adding charges n/2 (mod n) for the fields H u , H d , u, d, e, and ν. This amounts to simply toggling whether matter-parity-violating terms are allowed or forbidden by Z R n . We therefore always take z ν = r, rather than the alternative z ν = r + n/2 in the case that n is even, with the understanding that there is always a corresponding discrete symmetry that can be obtained in this way. In the case of odd n, one can also always choose to impose matter parity, or not.
All of the charges and r are defined (mod n). If n is even, then any Z R n symmetry with r = n/2 is equivalent to the corresponding Z n non-R symmetry from the low-energy point of view, because the two symmetries forbid and allow precisely the same Lagrangian terms. Also, if r > n/2, then one can replace (z Φ , r) → (n−z Φ , n−r) to obtain an equivalent version of the discrete symmetry. Therefore, for the R-symmetries, one can take 0 < r < n/2 without loss of generality.
For any given discrete symmetry, equivalent versions of it can be obtained by multiplying Charges for discrete symmetry (Z n for r = 0, or Z R n for non-zero integer 0 < r < n/2), for chiral superfields in the base models, in terms of two integers h and x, in the convention adopted here. In the case of B IV , m = 0 for odd n, while m = 0, 1 for even n and r = 0, and m = 1 for even n if r = 0. If vectorlike pairs are added to the model, their net discrete symmetry charges can be obtained from those for X and/or Y , given the mass term as specified in eq. (1. 18 all of the charges (including r) by a common integer relatively prime to n, and then taking the results (mod n). We use this to eliminate redundant low-energy descriptions of a given discrete symmetry. We also reduce every discrete symmetry to the smallest n that describes it, by removing any common factors from the list of charges z Φ and r. It follows that the Z n or Z R n charges of the base model chiral fields can always be written, by a conventional choice, in terms of only two independent charges x and h (and one binary choice m = 0, 1 for model B IV when n is even and r = 0), as summarized in Table 3.2. Therefore, for the sake of brevity, we will not list discrete symmetry charges for H d and the MSSM quark and lepton fields below, since they can be obtained from Table 3.2. For extensions of the base models, the total discrete symmetry charges for each vectorlike pair can be obtained from those for X and/or Y , given the mass term as specified in eq. (1. 18 Gaugino masses, and holomorphic soft-supersymmetry breaking terms corresponding to terms in the superpotential, are forbidden by unbroken discrete R-symmetries. Since they are a phenomenological necessity, a discrete R-symmetry must be spontaneously broken. This can be done with a spurion with Z R n charge 0, whose F -term component (with Z R n charge equal to −2r) obtains a VEV. In the low-energy theory, this will give rise to the usual MSSM supersymmetry breaking terms including gaugino masses, as well as PQ-violating terms of the form in eq. (1.12). As already noted in the Introduction, the resulting axion tadpole contribution is parametrically of the same order as the corresponding superpotential terms in eq. (1.11), and therefore has the same order-of-magnitude effect on θ eff .

B. Discrete non-R symmetries Z n
In Table 3.3, we give a complete list of the possible distinct non-R (r = 0) Z n symmetries for the base models, consistent with a bare Weinberg operator (H u )(H u ) (or equivalently ν ν) for neutrino masses and the anomaly constraint A 2 = A 3 (mod n). In the second column of the table, we show the resulting suppression factor p. We see that there is a large selection of Z n symmetries that satisfy the weaker anomaly-free constraint. These include, for example, a Z 22 symmetry with p = 11 for base model B IV , equivalent to one previously proposed and studied in ref. [60]. In our Table 3.3, this symmetry has (X, Y, H u ) charges (2, 6k + 3, 4k − 2) with k = 1; the charges listed in ref. [60] are different but the symmetry is equivalent in the sense discussed above. There are many other discrete symmetries with equal or greater suppression p. For example, with the same base model B IV , we see that there are also Z 16 and Z 20 symmetries that provides protection up to p = 10 and 12, respectively. Despite this infinite number of possible Z n symmetries, the only base model cases that satisfy the additional condition A 1 = 5A 3 (mod n) and have adequate suppression (p ≥ 7) of PQ-violating terms in the superpotential are Z 36 for model B III with p = 12 (with three distinct possibilities for the H u charge, labeled in the table by m = 0, 1, 2), and the Z 36 for model B IV with p = 8. In both cases, the Green-Schwarz mechanism is needed, as ρ GS = A 3 = 18 (mod 36).
If there are extra vectorlike chiral supermultiplets that owe their masses to X and Y , then there are so many possible discrete symmetries that we cannot provide an exhaustive list. Some examples are given in Tables 3.4 and 3.5, in which we have limited ourselves to cases with at most one 5 + 5 or 10 + 10 of SU (5) at the TeV scale, and in which both anomaly-free constraints A 2 = A 3 (mod n) and A 1 = 5A 3 (mod n) are required. The order n, the Z n MSSM charges, and the suppression p for these symmetries are the same as found in Table 3.3, but the anomalies N and E that enter into the low-energy axion phenomenology are different because of the contributions of the vectorlike fields. We also note that there are several examples that have ρ GS = 0 and therefore do not require the Green-Schwarz mechanism. The examples with ρ GS = 0 and N = 0 that have the smallest order n include Z 20 for extensions of B II and B III , and Z 12 and Z 36 for extensions of B IV . Each of these has p = 8, which provides sufficient quality for the PQ symmetry for f A ∼ < 4 × 10 9 GeV if the PQ-violating terms have order one couplings with generic phases, and for larger f A otherwise. We note that when n is not too large, ρ GS = 0 tends to allow for higher p.
Also shown in Tables 3.4 and 3.5 are some examples that have ρ GS = 0 but also N = 0, so that there is no PQ-QCD-QCD anomaly. These examples will have a light axion-like particle (ALP), but will not provide a solution to the strong CP problem. One peculiar case is that of base model B IV extended by X 2 QQ + X 2 U U + XY EE. This particular vectorlike content implies a completely anomaly-free PQ symmetry with E = N = 0, resulting in an infinite number of different Z n symmetries that protect it; only the first one (with n = 14) is shown. However, all cases that give N = 0 are moot for our main purpose in this paper, as there is no reason to forbid high-dimension PQ-violating terms if the strong CP problem cannot be solved anyway.
The examples shown in Tables 3.4 and 3.5 are only some of the many possibilities. There are similar Z n symmetries available for models that have other combinations of vectorlike supermultiplets both at the TeV and intermediate scales. This is also true of the models with quixes, as for example in the last two rows of Table 2.2.  TABLE 3.3: Exhaustive list of distinct Z n symmetries that satisfy the discrete anomaly cancellation constraint A 2 = A 3 (mod n) for the MSSM base models including the bare Weinberg operator (H u )(H u ) for neutrino masses. The meaning of p is that the lowest dimension PQ-violating superpotential term(s) allowed by the discrete symmetry and involving only X and Y are of the form X j Y p−j . In this table, k is any integer, and m = 0, 1, 2. The middle three columns are the Z n charges for X, Y , and H u . The H d and quark and lepton superfield charges are then determined as in Table 3.2, with the charges of q and H u set equal to 0 without loss of generality, as discussed in the text. Despite this infinite number of possible Z n symmetries, the only cases that satisfy the additional condition A 1 = 5A 3 (mod n) and have adequate suppression of PQ-violating terms (p ≥ 7) are Z 36 for model B III with p = k = 12, and the Z 36 for model B IV with p = 8 and k = 1.  : Some examples of non-R discrete symmetries Z n satisfying the anomaly cancellation conditions A 2 = A 3 (mod n) and A 1 = 5A 3 (mod n), obtained from base models B I , B II , or B III by adding up to one 5 + 5 or 10 + 10 of SU (5) at the TeV scale. The suppression of PQ violation p is defined so that the lowest dimension PQ-violating superpotential term(s) allowed by the discrete symmetry and involving only X and Y and are of the form X j Y p−j . In all cases in this table, the X and Y charges are respectively 1 and −3, and the H u charges are as listed with m = 0, 1, 2. The H d and quark and lepton superfield charges are then determined as in Table 3.2. The last three columns give the Green-Schwarz contribution ρ GS to the Z n anomalies and the PQ-QCD-QCD and PQ-EM-EM anomalies N and E. The five cases with ρ GS = 0 do not require a Green-Schwarz mechanism. However, of those, the two cases with N = 0 have no PQ-QCD-QCD anomaly and therefore have an axion-like particle but do not provide a solution to the strong CP problem.  : Some examples of non-R discrete symmetries Z n satisfying the anomaly cancellation conditions A 2 = A 3 (mod n) and A 1 = 5A 3 (mod n), obtained from base model B IV by adding up to one 5 + 5 or 10 + 10 of SU (5) at the TeV scale. The meaning of p is that the lowest dimension PQ-violating superpotential term(s) allowed by the discrete symmetry and involving only X and Y and are of the form X j Y p−j . The next three columns list the charges of X, Y , and H u , with m = 0, 1, 2. The remaining MSSM charges are determined in each case as in Table 3.2. The last three columns give the Green-Schwarz contribution ρ GS to the Z n anomalies and the PQ-QCD-QCD and PQ-EM-EM anomalies N and E. The cases with ρ GS = 0 do not require a Green-Schwarz mechanism. However, the last two examples shown, with N = 0, have no PQ-QCD-QCD anomaly and therefore have an axion-like particle but do not provide a solution to the strong CP problem.

Base Extension
Discrete R-symmetries provide even more possibilities for an accidental Peccei-Quinn symmetry protected to a high power p, including a larger number of cases that satisfy both A 2 = A 3 (mod n) and A 1 = 5A 3 (mod n). A non-exhaustive selection of such symmetries for the base models is shown in Table 3.6. One of the discrete Z R 24 symmetries had been previously proposed and studied for the MSSM in ‡ ref. [61], and was found in ref. [29,64] to extend to base models B II and B III with suppression p = 10. There are actually several inequivalent Z R 24 symmetries; the one found in [61] and applied to B II and B III in [29,64] is equivalent (after shifting by a multiple of weak hypercharge) to what we have listed in Table 3.6 in the second row under B II and the third row under B III , with k = m = 0 and ρ GS = 18 in each case. It was also pointed out in [29,64] that when this Z R 24 symmetry is realized in base models B I and B IV , it only provides suppression p = 7; we have chosen not to list these among the examples, although suppression p = 7 may be good enough if f A is up to 4 × 10 9 GeV, since it can be promoted to p = 8 by imposing an extra Z 2 symmetry that acts only on X and Y .
For model B IV , there are two distinct Z R 12 symmetries that give p = 7 (or p = 8 by the trick just mentioned), one of which requires use of the Green-Schwarz mechanism, and one of which does not because ρ GS = A 1 = A 2 = A 3 = 0 (mod 12). For ρ GS = 0 with protection up to p = 8, one also has symmetries with smallest order n = 54 for models B I , B II and B IV , and with protection up to p = 10 for n = 54 for models B III and n = 108 for models B I and B IV . There is a plethora of other possibilities beyond those shown in the table.
Adding vectorlike supermultiplets adds to the possibilities for discrete R symmetries. A few examples are shown in table 3.7, limited for reasons of brevity to only a subset of the cases that have ρ GS = 0, so no Green-Schwarz mechanism necessary, and that have at most one SU (5) multiplet at the TeV scale.
We have checked that for every available pair of anomaly coefficients N and E that determine the low-energy axion phenomenology as discussed in the next section, one can find a variety of corresponding discrete symmetries. From the point of view of low-energy phenomenology, the exact identity of the discrete symmetry may be of limited interest, since it is not possible to determine whether the discrete symmetry is an R-symmetry, or its order n, or whether it may have ρ GS = 0. The more important point seems to be the existence proof that it is always possible to find such discrete symmetries, so that the global PQ symmetry is consistent. The correlation of N , E, and the possible presence of TeV-scale vectorlike quarks or leptons could eventually point the way to specific ultraviolet completions.

D. Impact on baryon number and lepton number violation
Recall that in the MSSM there are renormalizable superpotential terms that violate lepton number L and baryon number B, schematically: (3.9) ‡ Ref. [61] imposed a requirement that the discrete symmetry charges for MSSM quark and lepton superfields respect SU (5) invariance, which limits the possibilities to only a few. We do not impose this requirement. We note that all of the symmetries we find can be made consistent with a partial unification with a Pati-Salam SU (4) PS × SU (2) L × U (1) R [87] embedding, although this is not always immediately obvious with our discrete charge conventions. : Some examples of discrete R-symmetries Z R n that satisfy both discrete anomaly cancellation constraints A 2 = A 3 (mod n) and A 1 = 5A 3 (mod n) for the base models. The superpotential has charge 2r and gauginos have charge r. The meaning of p is that the lowest dimension PQ-violating superpotential term(s) allowed by the discrete symmetry and involving only X and Y and are of the form X j Y p−j . The next three columns list the charges of X, Y , and H u in terms of k, k = 0, 1, and m = 0, 1, 2. The remaining MSSM charges are determined in each case as in Table 3   : Some examples of discrete R-symmetries Z R n that satisfy the anomaly cancellation conditions A 1 = A 2 = A 3 = 0 (mod n) without the Green-Schwarz mechanism (so ρ GS = 0). The superpotential has charge 2r and gauginos have charge r. The meaning of p is that the lowest dimension PQ-violating superpotential term(s) allowed by the discrete symmetry and involving only X and Y and are of the form X j Y p−j . The next three columns list the charges of X, Y, H u , in terms of k = 0, 1, and m = 0, 1, 2, and s = 0, 1, 2, 3. The remaining MSSM charges are determined in each case as in Table 3 Taken together, these would predict very rapid proton decay. The most common way of avoiding this is to impose the Z 2 matter parity = (−1) 3(B−L) discrete symmetry (or equivalently R-parity) [88]- [90]. There are also non-renormalizable operators that suppressed by the cutoff scale but violate both B and L, and so could directly mediate proton decay in violation of current bounds: As was noted in [60] for the Z 22 discrete symmetry example found there, and in [61] for their Z R 24 symmetry, the discrete symmetry that protects the PQ symmetry can also help by forbidding these dangerous baryon number and lepton number violating operators.
In Table 3.8, we show the discrete symmetry charges z O − 2r of the above superpotential operators O. The operator is allowed if the entry in the table vanishes. Note that the lepton-number violating operator H u is always forbidden if the discrete symmetry is an R symmetry, but is always allowed if it is a non-R symmetry (with r = 0). This is because we required that the square of this operator is allowed to provide neutrinos masses. Also, the discrete symmetry always forbids the operators q d and e. This is because if they are allowed, one can check that the term Y 4 would also be allowed, and would violate the PQ symmetry, against our defining requirement for the discrete symmetry. Therefore, the PQ-protecting Z R n discrete symmetry never allows renormalizable L violation, and Z n only allows soft L violation. Regarding the other B and L violating operators, the charges in the table have no particular reason to vanish, and we find that in the vast majority of discrete symmetry cases they do not. A catalog of specific cases will not be attempted here, but the general lesson is that quite often the discrete symmetries are powerful enough to automatically suppress proton decay sufficiently to satisfy present bounds, in addition to providing the high-quality PQ symmetry. One can also always supplement the PQ-protecting discrete symmetry with either matter parity, or the Z 3 baryon triality [82], which forbids all proton decay.

IV. AXION SIGNALS AND DETECTION PROSPECTS
In this section, we discuss the low-energy axion couplings to photons, electrons, and nucleons, and their impact on limits and detection prospects, for the supersymmetric base models and their extensions discussed above. For comparison, we will also give results for the standard non-supersymmetric benchmark QCD axion models of Kim-Shifman-Vainshtein-Zakharov (KSVZ) [6,7], and DFSZ [8,9] of types I and II.
In the non-supersymmetric benchmark DFSZ-I model, the Standard Model leptons and down-type quarks get their mass from one Higgs doublet, and the up-type quarks couple to the other Higgs doublet, as is the case with the MSSM. In DFSZ-II, the Standard Model leptons instead couple to the Higgs doublet that provides for the mass of the up-type quarks. In the KSVZ models, the Standard Model quarks and leptons are not charged under U (1) PQ , unlike the DFSZ models. We consider several KSVZ constructions in which the heavy vectorlike quarks mediating the PQ anomaly have different electroweak quantum numbers. The original Kim model [6] in which the vectorlike quark transforms under the Standard Model gauge group as (3, 1, 0) + (3, 1, 0), will be called KSVZ 0 , while KSVZ D , KSVZ U , KSVZ Q have vectorlike quarks with gauge quantum numbers as that of D + D, U + U , Q + Q, respectively. Table 4.1 summarizes these benchmarks, with their PQ-QCD-QCD anomaly coefficients N , and the coupling coefficients defined in eqs. (2.20) and (2.21) and appearing in eqs. (2.27)-(2.32) and (2.34). Among these benchmarks models, only KSVZ 0 , KSVZ D , and KSVZ U models have domain wall number N DW = 1.
As noted in ref. [91], all four supersymmetric DFSZ base models B I , B II , B III , and B IV have E/N = 2, which is close to cancelling the model-independent contribution −1.92 ± 0.04 in eqs. (2.27) and (2.49), and so have suppressed axion-photon couplings compared to the standard non-supersymmetric benchmark cases. However, as we will see below, adding extra vectorlike supermultiplets that get their masses from the VEVs of the X and Y scalar fields can enhance the axion couplings. This is because they have non-zero net PQ charges contributing to the anomaly coefficients N and E. Interestingly, the extensions with N DW = 1 that avoid the cosmological domain wall problem have the smallest |N | and therefore can give rise to the most enhanced low-energy axion couplings. As discussed in Section III, for every (N, E) found in the base models and their extensions, the PQ symmetry can be protected to a high degree (generally in a variety of different ways). Consequently, high-quality QCD axions in many of these extensions, including some N DW > 1 cases, should be accessible to future axion searches. We begin with the axion coupling to electrons, because this gives the most robust bound on f A from present astrophysical data. Figure 4.1 shows the axion-electron coupling |g Ae | as a function of the axion mass m A [and equivalently f A , through eq. (2.26)], for large tan β. For the DFSZ-I benchmark model and all supersymmetric models, the results are not very sensitive to the value of tan β as long as it is large, and |g Ae | = m e /|N |f A in the large tan β limit, so the lines can be simply labeled by |N |. The line |N | = 3 is thus the common result for DFSZ-I and for all four base models. It is surrounded by a green shaded region bounded by the labels |N | = 3/2 and |N | = 6, which is the range of possibilities for all extensions of the base models with only one 5 + 5 at either the intermediate or TeV scale. The larger blue shaded band between |N | = 1/2 and |N | = 33/2 is the range of results for the more general class of models obtained by extending the base models to include a 5 + 5 or 10 + 10 at the TeV scale and/or 5 + 5 or 10 + 10 at M int , or 2 × (L + L) at the TeV scale and an exotic quix pair D 6 + D 6 at M int . The largest possible couplings g Ae for a given f A are obtained for |N | = 1/6 for extensions of B III , and |N | = 1/2 for extensions of B I,II,IV . There are no extensions that have |N | in between the allowed values 1/2 and 1/6, so we did not include that range in the blue shaded band. We also show the results for KSVZ 0 and DFSZ-II with tan β = 10. Note that |g Ae | is much smaller, and more sensitive to large tan β, in DFSZ-II than in DFSZ-I. In these cases, since the tree-level g Ae is zero and very small respectively, we include the leading log renormalization contribution from eq. (2.32).
The models with N DW = 1 saturate the upper limit on the axion-electron coupling as shown in the figure. These include extensions of base models B I , B II , and B IV with |N | = 1/2 and B III with |N | = 1/6 shown in Table 2.2. There are also extensions of base model B III that have |N | = 1/2 but with N DW = 1. The lower limits in all extensions correspond to the models with largest |N | that can occur. In all base model extensions with SU (5) pairs, or D 6 + D 6 quixes, the largest |N | = 33/2 occurs in B II extended with a 10 + 10 at the TeV scale and a 10 + 10 at M int . And, in the extensions with exactly one 5 + 5, the largest |N | = 6 occurs in the extension of B II with a 5 + 5 at the TeV scale.
Current experimental limits from the LUX experiment [92], and much stronger bounds from the brightness of the tip of the red-giant branch [93,94] are shown as shaded regions in Figure 4.1. There are also bounds [95]- [97] from the cooling of white dwarfs, which are somewhat less strong than the red-giant bounds, and some hints of stellar cooling that might be ascribed to axions [98]. The red giant bound on the axion-electron coupling g Ae < 1.3 × 10 −13 sets the most stringent astrophysical constraint throughout our supersymmetric  The light blue region is the allowed range for extensions of base models that include a 5 + 5 or 10 + 10 at the TeV scale and/or 5 + 5 or 10 + 10 at M int , or 2 × (L + L) at the TeV scale and an exotic quix pair D 6 + D 6 at M int . The cases with N DW = 1 have either |N | = 1/2 (extensions of B I,II,IV , heavier solid blue line) or 1/6 (extensions of B III , solid purple line), and saturate the upper limit on |g Ae | in these extensions. The shading for allowed supersymmetric models excludes 1/6 < |N | < 1/2, where there are no cases. The much smaller predictions of the benchmarks DFSZ-II with tan β = 10 and KSVZ 0 are also shown for comparison. The current experimental limits on |g Ae | from the LUX experiment and from stellar cooling are also shown on the plot, as labeled.
DFSZ axion model space: We now turn to the axion-photon coupling and direct detection prospects for axions. In Figure 4.2, we show the axion-photon coupling |g Aγ | as a function of the axion mass m A and the axion decay constant f A . The axion-photon coupling depends only on the ratio E/N and f A , so each line in this plot corresponds to a value of |E/N − 1.92(4)|. The results for the base models (which have E/N = 2), including the uncertainty on the modelindependent contribution, are shown as the green shaded band. The blue shaded region shows the range of much larger couplings obtained for extended models made from SU (5) multiplet combinations that give N DW = 1. The E/N values that occur in these extensions were listed in Table 2.2. This blue shaded region is bounded from above by a solid blue line that corresponds to the extensions that have E/N = 68/3, and the lower bound for these models is E/N = 8/3 (which is the same as DFSZ-I). This shows that requiring N DW = 1 guarantees that the axion-photon coupling is at least as large as for the DFSZ-I model, and is usually considerably larger, for a given axion mass. The extensions of base models with 2 × (L + L) at the TeV scale and a quix pair D 6 + D 6 at M int can give an even larger axion-photon coupling compared to the extensions with SU (5) multiplets. This is shown by the dashed black line that corresponds to E/N = 104/3, which occurs in the extension of B II with N DW = 1.
Also shown in Figure 4.2 is a light gray shaded region obtained for models with N DW > 1, obtained by extending the base models to include a 5 + 5 or 10 + 10 at the TeV scale and/or 5 + 5 or 10 + 10 at M int . The upper limit in these models with N DW > 1, shown with a solid brown line, is realized by an extension with a 10 + 10 at the TeV scale and another 10 + 10 at M int with E/N = −20/3, N DW = 3. The base models extended with only a 5 + 5 at the TeV scale or M int cannot have N DM = 1, and have an upper limit for |g Aγ | shown as the solid red line in Figure 4.2. This corresponds to E/N = 17/3, realized in an extension of the model B II with a 5 + 5 at the TeV scale. The axion-photon coupling can vanish within errors in some of the extensions (including quixotic extensions) of base models when E/N happens to be close to the imperfectly known model-independent contribution −1.92 ± 0.04, as will be further illustrated below in the bottom panel of Figure 4.3. This is the reason for shading the entire region in light gray below the solid brown line. Also shown in Figure 4.2 are the lines for the non-supersymmetric DFSZ-I, DFSZ-II, and KSVZ 0 benchmark models. Note that, as can be seen in Table 4.1, DFSZ-I and DFSZ-II have the same axion-photon couplings as KSVZ U and KSVZ D , respectively.
The current experimental bounds on g Aγ as a function of m A are shown as shaded regions in the top panel of Figure 4.2. These include limits from the helioscope CAST [99,100] and from evolution of Horizontal Branch (HB) stars [101]. By coincidence these happen to give almost the same limit g Aγ < 6.5 × 10 −11 GeV −1 over a very wide range of f A , although the CAST bound becomes weaker for f A ∼ < 4 × 10 8 GeV. Also shown are the results of searches over much narrower bands in f A by haloscopes, under the assumption that axions are the dark matter: ADMX [102][103][104][105], CAPP [106], RBF [107], UF [108], HAYSTAC [109,110], QUAX [111], and ORGAN [112]. (In making Figures 4.2 and 4.4, we made substantial use of the axion limit data collected at ref. [113].) This shows that some limited ranges of f A for extended supersymmetric models with N DW = 1 have already been probed. In the lower panel of Figure 4.2, the shaded regions with dashed borders show the projections for future sensitivity from the helioscope IAXO [114], and haloscopes (ADMX [115], KLASH [116], MADMAX [117], Plasma haloscope [118], and TOORAD [119]). Except for a narrow range in the case of TOORAD, these do probe the supersymmetric DFSZ base models, but it is encouraging that the haloscopes do cover many of the possibilities for the model extensions, including essentially all of the allowed range predicted by our models with N DW = 1. However, it is again important to note that the haloscope search projections assume that axions are the main component of the dark matter.
The axion couplings to electrons and photons both scale with 1/f A . It is therefore interesting to compare the ratios g Ae /g DFSZ−I,tan β=10 Ae and g Aγ /g DFSZ−I Aγ , in which the dependence on the scale f A (and the axion mass) very nearly cancels. We choose the normalizing denominators to be the results for the non-supersymmetric benchmark DFSZ-I (with tan β = 10). Figure 4.3 shows a scatterplot of these ratios for the base models and their extensions along with the benchmarks DFSZ-I, DFSZ-II, and KSVZ Q,U,D,0 , as labeled. In the plot, the extensions of base models are categorized based on the additional particle content at the TeV scale, which could eventually be discovered in collider experiments. First, there are models with a 5 + 5 or 10 + 10 at the TeV scale with or without a 5 + 5 or 10 + 10 at M int , then there are models with 2 pairs of vectorlike leptons L + L at the TeV scale and D 6 + D 6 quixes at M int , and lastly there are models with no new particle content at the TeV scale which includes the extensions with a 5 + 5 or 10 + 10 at M int . In these extensions, different combinations of the superpotential mass terms involving the X and Y fields for the additional vectorlike supermultiplets give rise to different points in the figure. In the top panel of Figure 4.3, we only show cases with N DW = 1, whereas in the bottom panel we show all possibilities (i.e. N DW ≥ 1) that occur in these extensions. Also, in the top panel the uncertainty bars are only shown for the base models, for which the uncertainties are significant due to the proximity of E/N = 2 to 1.92 (4). In the bottom panel, the g Aγ uncertainty bars are shown for all points.
Except for the f A dependence that enters the ratio g Ae /g DFSZ−I,tan β=10 Ae through ∆C Ae for models in which C Ae vanishes at tree level, the plots in base model B III with |N | = 1/2, N DW = 3. Also, in the figure, the DFSZ-I and DFSZ-II benchmark models are shown for various values of tan β. We can therefore understand the impact of tan β on the points for base models and its extensions, as they have a similar tan β dependence as that of DFSZ-I. Finally, we turn to the axion-nucleon couplings, for which the neutron coupling is most accessible to constraint. Figure 4.4 shows the axion-neutron coupling as a function of the axion mass m A and the decay constant f A , choosing tan β = 10 for plotting purposes. Due to significant uncertainties in the axion coupling predictions, each case in the figure is a shaded band instead of a line. First, bands are shown for the N DW = 1 cases that occur in N = ±1/2 and ±1/6 extensions of the base models. The axion-nucleon coupling, in general, depends on the sign of N , but in the case of the neutron the shaded bands for N = 1/2 and N = 1/6 very nearly overlap with the shaded bands for N = −1/2 and N = −1/6, respectively, so we have combined them. Also shown is the band of maximum allowed axion-neutron coupling in the extensions with only a 5 + 5, which occurs for N = 3/2 in base model B I extended with a 5 + 5 at M int . In addition, ranges are shown for the supersymmetric base models and the DFSZ and KSVZ models, as labeled. Note that for a fixed tan β and N , the base models and DFSZ-I and DFSZ-II all have the same axion-nucleon couplings. Within uncertainties, the KSVZ axion can accidentally decouple from the neutron, so the shading in that case extends to arbitrarily low g An .
The existence of the neutrino signal from Supernova SN 1987A has been interpreted by ref. [120] to put the following constraint on axion-nucleon couplings: However, it has been suggested e.g. in ref. [121] that such bounds should be taken as a guide rather than a sharp bound, given the uncertainties involved. Taken at face value the bound implies, for large tan β, f A ∼ > 0.15 + 0.66/N 2 (10 9 GeV), (4.4) for our models, by using eqs. (2.51) and (2.52). (There is also a term proportional to 1/N under the square root, but its coefficient is consistent with 0 within uncertainties.) For the critical case of small |N |, this is not quite as strong as eq. (4.1) obtained above from the constraints on g Ae from red-giant and white dwarf cooling. In Figure 4.4, we show the bound as it applies to the supersymmetric DFSZ models as parameterized by eqs. (2.51) and (2.52) with large tan β, recognizing that it is not a model-independent bound on g An (because other models can have different correlations between g An and g Ap ) and in any case may not be robust in detail. There is also an even stronger candidate bound [122] of |g An | < 2.8 × 10 −10 from avoiding too much cooling of the hot neutron star HESS J1731-347. However, since it is difficult to account for the temperature of this particular object even without axionic cooling [123], we do not include this bound in Figure 4 The axion-neutron coupling as a function of the axion mass and the axion decay constant. The model predictions for large tan β are diagonal bands including uncertainties in the coupling. As labeled, they correspond to the supersymmetric base models (same as nonsupersymmetric DFSZ-I and DFSZ-II benchmarks), the extensions with N DW = 1 that can have enhanced axion-nucleon couplings, and the case that gives rise to the largest axion-nucleon coupling in the extension with a single 5 + 5. Within uncertainties, the neutron may decouple from KSVZ axions, as emphasized by the orange shading extending to very low couplings. The experimental bound from supernova SN 1987A is also shown as it applies to the supersymmetric DFSZ models as parameterized by eqs. (2.51) and (2.52) with large tan β, taking into account the correlation with g Ap . Also shown are future projections for CASPEr wind (phase 2) and ARIADNE. The projection for the ARIADNE experiment is based on an optimistic choice regarding the CP-violating coupling of axions to neutrons, as discussed in the text.
Phase-II of the proposed CASPEr wind experiment [124], and a projection for the ARIADNE experiment [125,126]. ARIADNE would be sensitive to the product of g An and the CP-odd scalar coupling g s An (as in L CP−odd ⊃ −Ag s An Ψ n Ψ n ), and the region shown in Figure 4.4 is based on an optimistic choice of g s An = 10 −12 GeV/f A [20,125,127]. The coupling g s An could arise from a small non-zero CP violation in θ eff , and so cannot be uniquely predicted by the other relevant model parameters N and tan β. We note that, at least with this optimistic assumption, ARIADNE can probe g An for f A up to about 2.5 × 10 10 GeV for the supersymmetric DFSZ base models, which are invisible to the projected searches for g Aγ .

V. CONCLUSION
In this paper, we considered supersymmetric DFSZ type axion models with the field content of the MSSM plus two gauge-singlet fields X, Y that spontaneously break the PQ symmetry, and some extra vectorlike quark and lepton supermultiplets. These models simultaneously give a solution to the µ problem and the strong CP problem. The extra vectorlike content is chosen such that the perturbative gauge coupling unification is maintained. The PQ symmetry is spontaneously broken by the scalar components of X and Y , which acquire intermediate scale vacuum expectation values, simultaneously giving rise to a high-quality nearly invisible axion with a decay constant within the current astrophysical limits and a µ term around the TeV scale. The axino and saxion also have masses of order the TeV scale. Different combinations of mass terms for the additional vectorlike supermultiplets involving X and Y result in different combinations of the PQ anomaly coefficients (N, E), which in turn determines the low-energy axion couplings.
For the base models and their extensions, we studied how to obtain the PQ symmetry as an accidental symmetry that is guarded against the dangerous PQ violating superpotential terms of the form X j Y p−j /M p−3 P to an extent that is compatible with the experimental constraint on the QCD θ parameter, by imposing anomaly-free discrete non-R or R Z n symmetries with or without the Green-Schwarz mechanism. If the axion decay constant f A is as low as 10 9 GeV, we can allow p = 7, and for f A 10 12 GeV, we may instead need a higher suppression of up to p = 12. In order for the Z R n symmetries (which reduces to non-R Z n symmetry if the Z R n charge of the gauginos r = 0) to be anomaly-free, we impose a weaker constraint along with the additional stronger constraint on the Z R n ×G×G anomalies for G = SU (3) c , SU (2) L , and U (1) Y as discussed in Section III. Out of all possible non-R Z n symmetries for the base models, the only cases with adequate suppression (p ≥ 7) that satisfy the stronger set of anomaly constraints are Z 36 for base model B III with p = 12, and a Z 36 for B IV with p = 8. Both of these cases require the Green-Schwarz mechanism. On the other hand, there are a lot more anomaly-free Z n R-symmetries for each base model, some of which do not require the Green-Schwarz mechanism. For example, there is a Z R 54 for base model B III with p = 10, Z R 12 for base model B IV with only p = 7, both of which do not require the Green-Schwarz mechanism. With the Green-Schwarz mechanism, there are many more Z R n symmetries that have lower n and higher suppression p. In our approach, the discrete symmetry is fundamental and exact (being anomaly free) and the PQ symmetry is not fundamental and approximate (being an accidental consequence of the discrete symmetry).
Adding vectorlike supermultiplets adds greatly to the possibilities for both discrete non-R and R symmetries. For each possible pair of PQ anomaly coefficients (N, E) in the extensions, we find that there are always anomaly-free discrete symmetries that protect the PQ symmetry to a high degree of accuracy. Not only do these discrete symmetries provide for an accidental U (1) PQ symmetry, but they also can forbid dangerous baryon number and lepton number violating operators, in most cases, that could mediate dangerous proton decay. There are so many available discrete symmetries that we did not attempt a complete categorization. The additional fields and the discrete symmetry of the models we have proposed increase the complexity of the MSSM. The discrete symmetry could arise from an additional U(1) symmetry in the ultraviolet, but other than that we do not have any insight to offer as to the reasons for this additional complexity, other than the problems that it solves.
In our models, the vectorlike supermultiplets are coupled to Standard Model quark and lepton superfields to avoid cosmological problems. These couplings could be very small in magnitude while still allowing the vectorlike particles to decay promptly. Thus, while they would in general include additional CP-violating phases, the resulting CP-violating effects on the Standard Model could easily be negligible if the magnitudes of the couplings are sufficiently small. We also note that in the MSSM there are potential CP-violating effects from gaugino masses and the µ parameter. Our framework has nothing to say about these effects, except that they are not worse than in the usual MSSM, and as usual they are ameliorated for superpartners in the TeV range.
Due to the soft supersymmetry-breaking terms, the scalar components of the vectorlike chiral supermultiplets tend to be heavier than their corresponding fermions, as studied for example in ref. [128]. Therefore, the constraints on the masses of the extra vectorlike supermultiplets come from the searches for pair-production of vectorlike quarks and leptons at the particle colliders. From the most recent LHC searches by the ATLAS and CMS collaborations for vectorlike quarks [129,130], the up-type (down-type) vectorlike quarks that are assumed to decay to the top (bottom) quarks are excluded at 95% confidence level up to 1310 GeV -1600 GeV (1200 GeV -1570 GeV), depending on its branching ratios. And, the weak isodoublet vectorlike leptons that decay to the tau leptons are excluded at 95% confidence level from 120 GeV to 790 GeV by the CMS collaboration in ref. [131]. Lastly, as it was pointed out in the refs. [132,133], the weak isosinglet charged vectorlike leptons that mix with the tau lepton have essentially no reach prospects at the current and the future proton-proton colliders due to low pair-production cross-sections and unfavorable branching ratios.
After the PQ breaking the axion potential typically acquires more than one inequivalent degenerate minimum, given by the domain wall number N DW , leading to a cosmological domain wall problem if the symmetry is broken in the post-inflationary era. Models with N DW = 1 evade this problem. The four base models without extra vectorlike supermultiplets have N DW = 1, and therefore may suffer from the domain wall problem. They also have suppressed axion-photon couplings compared to the standard benchmark QCD axion models, and so may be invisible to future proposed direct axion searches. However, in the extensions of base models that include extra vectorlike supermultiplets, we obtained a wide variety of larger axion couplings. It is notable that the extensions with N DW = 1, which always have at least one strongly interacting vectorlike supermultiplet at the intermediate scale, have the smallest |N | and therefore give rise to enhanced axion couplings, which are likely to be within reach of future axion searches.