The Quest for an Intermediate-Scale Accidental Axion and Further ALPs

The recent detection of the cosmic microwave background polarimeter experiment BICEP2 of tensor fluctuations in the B-mode power spectrum basically excludes all plausible axion models where its decay constant is above $10^{13}$ GeV. Moreover, there are strong theoretical, astrophysical, and cosmological motivations for models involving, in addition to the axion, also axion-like particles (ALPs), with decay constants in the intermediate scale range, between $10^9$ GeV and $10^{13}$ GeV. Here, we present a general analysis of models with an axion and further ALPs and derive bounds on the relative size of the axion and ALP photon (and electron) coupling. We discuss what we can learn from measurements of the axion and ALP photon couplings about the fundamental parameters of the underlying ultraviolet completion of the theory. For the latter we consider extensions of the Standard Model in which the axion and the ALP(s) appear as pseudo Nambu-Goldstone bosons from the breaking of global chiral $U(1)$ (Peccei-Quinn (PQ)) symmetries, occuring accidentally as low energy remnants from exact discrete symmetries. In such models, the axion and the further ALP are protected from disastrous explicit symmetry breaking effects due to Planck-scale suppressed operators. The scenarios considered exploit heavy right handed neutrinos getting their mass via PQ symmetry breaking and thus explain the small mass of the active neutrinos via a seesaw relation between the electroweak and an intermediate PQ symmetry breaking scale. We show some models that can accommodate simultaneously an axion dark matter candidate, an ALP explaining the anomalous transparency of the universe for $\gamma$-rays, and an ALP explaining the recently reported 3.55 keV gamma line from galaxies and clusters of galaxies, if the respective decay constants are of intermediate scale.

IV. Intermediate-Scale accidental axion and further ALPs from Field Theoretic Bottom-Up SM Extensions 15 A. Z 13 ⊗ Z 5 ⊗ Z 5 model with axion-ALP mixing 16 B. Z n ⊗ Z m models with a photophilic ALP 22 1. Z 11 ⊗ Z 9 model with extended Higgs sector 23 2. Z 11 ⊗ Z 9 model with minimal SM Higgs sector 25 3. Z 11 ⊗ Z 7 model with minimal SM Higgs sector 27 C. Z 11 ⊗Z 9 ⊗Z 7 model with two photophilic ALPs 27 V. Summary and Outlook 29

Acknowledgments 31
A. Axion and ALP Couplings to Gluons, Photons, and Electrons 31 B. Effects from Gravity through Planck-Scale Suppressed Operators 33 C. Domain-Wall Problem in the Models 36 References 38

I. INTRODUCTION
The axion A is a strongly motivated very weakly interacting ultralight particle beyond the Standard Model (SM). Its existence is predicted in the course of an elegant solution to the strong CP problem [1][2][3], that is the non-observation of flavor conserving CP violation originating from the theta-angle term in the Lagrangian of QCD, involving the gluon field strength G a µν and its dual,G a,µν ≡ µνλρ G a λρ /2, with ε 0123 = 1. This solution consists in adding to the Standard Model a scalar field theory describing a (pseudo) Nambu-Goldstone boson arising from the breaking of a global chiral (Peccei-Quinn (PQ)) U (1) PQ symmetry: that is, the corresponding scalar field A(x) satisfies a shift symmetry A(x) → A(x) + const. which is only violated by the anomalous coupling to gluons, In fact, theθ-term can then be eliminated by absorbing it into the axion field, A +θf A → A. Moreover, nonperturbative QCD effects 1 provide for a non-trivial potential for the shifted axion field A -minimized at zero expectation value, A = 0, thus wiping out strong CP violation -and predict a mass for the particle excitation around this minimum, the axion A, in terms of the pion mass, m π = 135 MeV, its decay constant, f π ≈ 92 MeV, the ratio z = m u /m d ≈ 0.56 of up and down quark masses, and the decay constant f A . For large f A , the axion is an ultralight particle with very weak interactions with the Standard Model [4][5][6][7]. Its universal and phenomenologically most important interaction with photons [8][9][10][11], is tiny (see the yellow band in Fig. 1), since observations in astrophysics -in particular the observed duration of the neutrino signal from supernova SN 1987A -require a large decay constant f A (small mass m A ) [18], A further strong motivation for the axion comes from the fact that, for sufficiently large decay constant f A , 10 9 GeV f A 10 13 GeV ⇒ 10 −7 eV m A 1 meV, (6) it is a cold dark matter candidate, being non-thermally produced in the early universe by the vacuum-realignment mechanism and, in some models and under certain circumstances, also via the decay of topological defects such as axion strings and domain walls [19][20][21][22][23][24][25][26][27].
The upper bound on the decay constant in Eq. (6) follows from the recent discovery of the cosmic microwave background polarimeter experiment BICEP2 of tensor fluctuations in the B-mode power spectrum [28]. This implies a high value for the Hubble scale during inflation, which, together with constraints on isocurvature fluctuations [29][30][31][32][33], rule out plausible scenarios where inflation occurs after PQ symmetry breaking [34][35][36], that is where with the Gibbons-Hawking temperature during inflation and the maximum thermalization temperature after reheating. Here M Pl = 1.22 × 10 19 GeV is the Planck mass and 0 < eff < 1 is an efficiency parameter. Ways out of this conclusion have been put forward in Refs. [34,35]. The pre-inflationary PQ symmetry breaking scenario would have allowed, in principle, much higher values of the decay constant without overshooting the dark matter abundance, by invoking small values of the initial misalignment angle. After BICEP2, the plausible axion possibilities have narrowed down to scenarios with post-inflationary PQ symmetry breaking. A conservative upper bound on the axion decay constant is then A more stringent upper bound can be obtained by requiring that the axion dark matter abundance generated via the vacuum-realignment mechanism should not exceed the observed dark matter abundance, leading to [23-25, 31, 36] f A (1 ÷ 10) × 10 11 GeV, m A (6 ÷ 60) µeV.
There is also a theoretical motivation for the existence of additional axion-like particles (ALPs) (for reviews, see [50][51][52]), emerging as pseudo Nambu-Goldstone bosons from the breaking of other global symmetries, such as lepton number [53,54] or family symmetries [55][56][57]. Most notably, the low-energy effective field theories emerging from string theory generically contain candidates for the axion plus possibly a multitude of additional ALPs. Indeed, when compactifying the six extra spatial dimensions of string theory, Kaluza-Klein zero modes of antisymmetric form fields -the latter belonging to the massless spectrum of the bosonic string propagating in ten dimensions -appear as axion and further ALP candidates in the low-energy effective action [58][59][60][61][62][63], their number being determined by the topology of the internal compactified manifold. Moreover, the axion and a multitude of additional ALP candidates may also FIG. 1: Prediction of the axion-photon coupling versus its mass (yellow band). Also shown are excluded regions arising from the non-observation of an anomalous energy loss of massive stars due to axion or ALP emission [12], of a γ-ray burst (in coincidence with the observed neutrino burst) from SN 1987A due to conversion of an initial ALP burst in the galactic magnetic field, of changes in quasar polarizations due to photon-ALP oscillations, and of dark matter axions or ALPs converted into photons in microwave cavities placed in magnetic fields [13][14][15][16][17]. Axions and ALPs with parameters in the regions surrounded by the red lines may constitute all or a part of cold dark matter (CDM), explain the cosmic γ-ray transparency, and the soft X-ray excess from Coma. The green regions are the projected sensitivities of the light-shining-through-wall experiment ALPS-II, of the helioscope IAXO, of the haloscopes ADMX and ADMX-HF, and of the PIXIE or PRISM cosmic microwave background observatories.
arise as pseudo Nambu-Goldstone bosons from the breaking of accidental global U (1) symmetries that appear as low energy remnants of exact discrete symmetries occurring in string compactifications [64][65][66]. Intriguingly, very light and weakly coupled generic ALPs a i , with decay constants in the same range as Eq. (6), share with the axion the property of being cold dark matter candidates since they are also produced via the vacuum-realignment mechanism [67][68][69]. In fact, their relic abundance, in the now strongly favored post-inflationary symmetry breaking scenario, is roughly given by 2 Therefore, there is a strong motivation to search not only for the axion A, but also for additional light ALPs a i , for which the low-energy effective Lagrangian describing their interactions with (the lightest Standard Model particles) photons can generically be written as However, unlike the axion, which inherits many of its properties (m A , g Aγ ), from non-perturbative QCD effects associated with chiral symmetry breaking, the masses m ai and the photon couplings g aiγ of the additional ALPs are model dependent. Thus, there exists the possibility that ALPs are hierarchically more strongly coupled to photons than axions with the same mass and therefore easier to detect. Interestingly, there are indications from gamma-ray astronomy, which may point to the existence of at least one ALP beyond the axion. Gamma-ray spectra from distant active galactic nuclei (AGN) should show an energy and redshiftdependent exponential attenuation, exp(−τ (E, z)), due to e + e − pair production off the extragalactic background light (EBL) -the stellar and dust-reprocessed light accumulated during the cosmological evolution following the era of re-ionization. The recent detection of this effect by Fermi [70] and H.E.S.S. [71] has allowed to constrain the EBL density. At large optical depth, τ 2, however, there are hints that the Universe is anomalously transparent to gamma-rays [72][73][74][75][76]. This may be explained by photon ↔ ALP oscillations: the conversion of gamma rays into ALPs in the magnetic fields around AGNs or in the intergalactic medium, followed by their unimpeded travel towards our galaxy and the consequent reconversion into photons in the galactic/intergalactic magnetic fields [77][78][79][80][81]. This explanation requires a very light ALP, which couples to two photons with strength [82], Note that the entire parameter region (15) has no overlap with the universal prediction of the axion, see Fig. 1. In fact, an axion with m A 10 −7 eV would have a photon coupling |g Aγ | 10 −16 GeV −1 . Therefore, if this hint is taken seriously, it points to the existence of an ALP beyond the axion.
Intriguingly, an observed soft X-ray excess in the Coma cluster may also be explained by the conversion of a cosmic ALP background radiation (CABR) -produced by heavy moduli decay and corresponding to an effective number N eff ∼ 0.5 of extra neutrinos species [83,84] -into photons in the cluster magnetic field [85,86]. This explanation requires that the spectrum of the CABR is peaked in the soft-keV region and that the ALP coupling and mass satisfy respectively, overlapping with the parameter range (15) preferred by the ALP solution of the gamma-ray transparency puzzle, as is apparent in Fig. 1. Astrophysical bounds on ALPs arising from magnetised white dwarfs [87] and from the non-observation of a γ-ray burst in coincidence with neutrinos from the supernova SN 1987A [88,89] provide limits close to |g aγ | 10 −11 GeV −1 , for masses m a 10 −7 eV and m a 10 −9 eV, respectively, and thus cut into the parameter range (15) preferred by the cosmic γ-ray transparency anomaly. Even stronger limits, |g Aγ | 2 × 10 −13 GeV −1 , for m A 10 −14 eV, have been obtained by exploiting high-precision measurements of quasar polarizations [90,91]. But there remains still a sizeable region in ALP parameter space motivated by the above anomalies and at the same time consistent with all astrophysical constraints, as can be seen in Fig. 1.
At small masses below 10 −14 eV, a part of the region of interest in ALP parameter space will be probed indirectly by the next generation of cosmic microwave background (CMB) observatories such as PIXIE [92] and PRISM [93] (see the region labelled "PIXIE/PRISM" in Fig. 1, based on the assumption of an extragalactic magnetic field B of nG size; the projected sensitivity scales with the magnetic field as B −1 ), because resonant photon-ALP oscillations in primordial magnetic fields may lead to observable spectral distortions of the CMB [94][95][96].
A complementary part of the region of interest in parameter space will soon be probed by a pure laboratory experiment: the light-shining-through-a-wall (LSW) [97] experiment Any-Light-Particle-Search II (ALPS-II) is designed to detect photon-ALP-photon oscillations in the range [98] (see the green region labelled "ALPS-II" in Fig. 1) Further experimental opportunities covering this region in ALP parameter space will open if the International Axion Observatory (IAXO), a helioscope searching for solar axions and ALPs, is realized [99]. Its projected sensitivity is The latter instrument has also the possibility to probe the possible coupling of the axion or further ALPs to electrons, via their solar production by atomic axio-recombination, axio-deexcitation, axio-Bremsstrahlung in electron-ion or electron-electron collisions, and Compton scattering [100]. This is of considerable interest because of hints of an extra stellar cooling mechanism not accounted by the Standard Model. In fact, the white dwarf (WD) luminosity function seems to require a new energy-loss channel that can be interpreted in terms of losses due to sub-keV mass axions or ALPs, with Yukawa couplings [101], |g Ae | ≡ |C Ae | m e /f A ∼ 10 −13 and/or |g ae | ≡ | C ae |m e /f a ∼ 10 −13 , which are well in the range expected for an intermediate scale axion or further ALPs (if the model-dependent couplings C Ae and C ae are of order unity). The same parameter range is preferred to explain the anomalous size of the observed period decrease of the pulsating WDs G117-B15A and R548 by additional axion/ALP losses [102,103]. A third independent hint of anomalous stellar losses has recently been found in the red-giant branch of the globular cluster M5, which seems to be extended to larger brightness than expected within the Standard Model. A possible explanation of this observation is that the helium cores of red giants lose energy in axions or further ALPs with electron couplings of the same order as in Eq. (20) [104]. Very recently, two groups have found an unidentified X-ray line signal at 3.55 keV in the stacked spectrum of a number of galaxy clusters [105] and the Andromeda galaxy [106]. It is tempting to identify this line with the expected signal from two photon decay of 7.1 keV mass ALP dark matter [107][108][109][110]. To match the observed X-ray flux, but allowing for the likely possibility, that the ALP dark matter makes only a fraction x a ≡ ρ a /ρ DM of the total density of dark matter, the required lifetime and thus coupling of the ALP is [108] Therefore, it is timely to have a close look onto ultraviolet extensions of the Standard Model featuring, apart from an intermediate scale axion, also further ALPs and to investigate possible correlations between the low-energy axion and ALP couplings. In fact, as we will show in Sec. II, such correlations inevitably occur if there are originally two (or more) Nambu-Goldstone fields coupled to GG. The crucial determination of the decay constants and couplings of axion-like fields from original high-scale theories to gluons, photons, and electrons will be done in Secs. III and IV. In the latter section, we construct particularly well-motivated ultraviolet completions of the Standard model featuring accidental Peccei-Quinn symmetries and deduce their low-energy phenomenology. Finally, we summarize, conclude and give an outlook for further investigations in Sec. V.

II. CORRELATIONS BETWEEN AXION AND ALPS COUPLINGS IN MULTI-AXION MODELS
Up to now, most phenomenological studies have considered just one particle type at a time: either the axion or an ALP different from the axion, without taking into account possible relations between their low-energy parameters (see, however, Refs. [65,66,[111][112][113]). In this section, we will study, in multi-axion models, the phenomenologically most important axion and ALP couplings to photons and electrons in depth. We will find that there are possibly strong correlations in these couplings.
The most general low-energy effective Lagrangian of a model with n ax axion-like fields enjoying shift symmetries a i → a i + const., i.e. realizing a non-linear representation of n ax U (1) PQ i symmetries, coupling to gluons and photons, with field strengths G and F , respectively, and to electrons, reads [63], where f a i are the decay constants of the axion-like fields a i . The anomaly coefficients, C ig and C iγ , and the electron coupling, C ie , are typically of order unity, see Secs. III and IV 3 . The proper axion field A, that is the field which solves the strong CP problem, is the linear combination of the axion-like fields in front of the GG term in Eq. (22) [58], Its particle excitation, the axion A, mixes with the pion, rendering it a pseudo Nambu-Goldstone boson, whose mass has been given in Eq. (3). To this end, f A has to be chosen such that the kinetic term of A is normalized canonically. The remaining n ax − 1 axion-like fields a i , orthogonal to the axion A, are still massless Nambu-Goldstone bosons.
To be more explicit, let us specialize first to the phenomenologically well motivated two-axion model, n ax = 2 (see the appendix of Ref. [63] for a further exposition of the multi-axion case) and defer the discussion of the more general case to the end of this section. The properly normalized axion A and the additional ALP a are related to the original axion-like fields a i by with normalization [66] 1 The low-energy Lagrangian of this model, below the chiral symmetry breaking scale, is then where As in a single axion model [2,[8][9][10][11], the axion A has a universal, model-independent contribution to its coupling to the photon (the last term in Eq. (27)) which arises from its mixing with the pion (see, e.g., appendix of Ref. [63]). The parameters in the above expressions for the couplings are however redundant, because f A is constrained by Eq. (25). The physics is more obvious if one expresses the transformations in Eq. (24) in terms of mixing angles [65], A = a 1 cos δ + a 2 sin δ, a = − a 1 sin δ + a 2 cos δ, with In terms of this parametrization, the couplings can then be written as Effectively, the couplings depend on the dimensionful parameter f A , on the dimensionless ratios C 1γ /C 1g , C 2γ /C 2g , C 1e /C 1g , C 2e /C 2g , and on the mixing angle δ.
Apart from the model-independent contribution due to the mixing with the pion, the photon coupling of the axion has a slight, order one, model dependence due to the ratio of the anomaly coefficients C 1γ /C 1g , C 2γ /C 2g , and the mixing angle, cf. Eq. (33). On the other hand, the electron coupling of the axion A is very much model dependent: a non-vanishing g Ae requires that at least one of the original axion-like fields has a non-zero coupling to the electron, C ie = 0 for i = 1 and/or 2, cf. Eq. (34).
This strong model dependence is shared by the photon and electron couplings of the ALP a: here one needs that at least one of the original axion-like fields has a non-zero coupling to the photon (electron), C iγ = 0 for i = 1 and/or 2 (C ie = 0 for i = 1 and/or 2), otherwise g aγ (g ae ) vanishes, cf. Eq. (35) (Eq. (36)). In this context it is also important to mention that eventual explicit mass terms, m 2 a i a 2 i /2, for axion-like fields arising from breaking of the Peccei-Quinn symmetries, e.g. by higher dimensional operators suppressed by some high scale (see next section), will have negligible effects on the photon and electron couplings of the axion A and the ALP a, as long as the masses are much smaller than the dynamically generated axion mass, m a i m A (cf., e.g., appendix of Ref. [63]). We conclude that in the case that the axion A and the ALP a consist of an appreciable mixture of the original axion-like fields, i.e. as long as | tan δ| is of order unity, the couplings are all determined by f A and are therefore expected to be approximately of the same size for for the axion A and the ALP a, up to order one factors. A hierarchical difference between the axion and the ALP coupling can possibly only arrive in the situation where the axion (and correspondingly also the ALP) originates essentially only from one axion-like field, which occurs for δ ≈ 0 (π/2), meaning that cos δ ≈ 1 (sin δ = 1), i.e. that a 1 (a 2 ) constitutes the axion. In fact, concentrating without loss of generality on the former case, in which the sizes of the axion and ALP couplings appear to decouple from each other, the former (latter) being determined by f a 1 (f a 2 ). However, as long as C 2g = 0, the ALP couplings still can not be hierarchically larger than the axion couplings, since the above relations imply an upper bound on the ratio A true hierarchy in the couplings, in particular an ALP-photon coupling much larger than an axion-photon coupling, is only possible if there is no mixing at all, C 2g ≡ 0, but still C 2γ = 0, i.e. if the ALP a 2 is photophilic 4 . In this case we get such that e.g. |g aγ | can very well be in the range suggested by the ALP explanation of the cosmic γ-ray transparency puzzle and be accessible to the next generation of laboratory experiments, 10 −10 GeV −1 |g aγ | 10 −12 GeV −1 ; m a 10 −7 eV, (47) provided that f a 2 /|C 2g is in the range, while at the same time |g Aγ | can be right in the cosmic axion window, provided that f A = f a 1 /|C 1g | is in the range We therefore conclude that in models with multiple axion-like fields it is of very high phenomenological relevance to know whether several of them couple to GG simultaneously. This question will be considered in the following sections.
Let us end this section by noting that even further insight into the general relation between the axion and ALP couplings can be obtained by purely geometrical means. To this end, we rewrite Eqs. (27) and (29) as Aγ . where and Thusŷ points into the direction of the axion, A(x) =ŷ i a i (x), whereasẑ points into the direction of an ALP that couples to photons (the orthogonal direction decouples from photons). Plane geometry ensures us that we can write then Aγ , g aγ = g γmax sin Θ , where and Θ is the angle betweenŷ,ẑ, fromŷ toẑ, This implies see Fig. 2. In particular, the ALP-photon coupling is constrained by while the axion-photon coupling satisfies If we take the hierarchy f a 1 f a 2 , with order one coefficients C ig , C iγ = 0, we realize that the vectors in Eqs. (51) are both approximately aligned along (0, ±1). This situation leads to Θ = 0 or π and hence to Then f A ∼ f γ and |g aγ | |g Aγ |, barring accidental cancellations. The opposite hierarchy, |g aγ | |g Aγ |, is more difficult to be obtained without fine tuning and requires g (0) Aγ ≤ g γmax , or, equivalently, f γ ≤ f A /C 0 f A /1.95. Analyzing the various possibilities, the only way we can achieve such a hierarchy from hierarchical scales f a i is to require one of the axion-like fields to be photophilic.
Notice that Eq. (56) can be directly generalized for the case of n ALP > 1 ALPs, besides the axion. The constraint generalizes to where Eq. (54) still defines g γmax with the generalization where n ax = n ALP + 1 is the number of axion-like fields. The couplings in this case are given by where u i are orthonormal basis vectors spanning the space orthogonal toŷ. An analogous constraint for the electron couplings can be obtained by introducing the angle Θ e between the vector y in (51), pointing in the axion direction, and the vector pointing in the direction of an ALP that couples to electrons. With the help of this, we can rewrite Eqs. (28) and (30) in the form g Ae = g emax cos Θ e , g ae = g emax sin Θ e , where Therefore, the ALP-electron couplings are constrained by This constraint can be generalized also straightforwardly to the case of more than two axion-like fields.

III. AXION AND FURTHER ALPS IN AD-HOC PECCEI-QUINN SM EXTENSIONS
In this section we consider a number of field theoretic extensions of the SM where an axion or further ALPs occur as Nambu-Goldstone bosons from the breaking of ad-hoc U (1) Peccei-Quinn (PQ) symmetries. We will derive the low-energy couplings in terms of the fundamental parameters of the underlying high-scale theories, with the aim of elucidating, in the following section, to which extent current and upcoming axion and ALP experiments can probe them.

A. KSVZ-type models
We start with the so called KSVZ model [4,5]. In this construction, a color triplet, but SU (2) L singlet heavy vector-like fermion, Q = (Q L , Q R ), is added to the SM. In addition, a hidden (i.e. SM singlet) complex scalar field σ is introduced. The SM particles are assumed to be uncharged under the PQ symmetry, while the hidden scalar field and the exotic color triplet are supposed to transform as respectively. Correspondingly, the most general Yukawa interactions involving the PQ charge fields and the most general renormalizable scalar potential involving also the SM Higgs field H read respectively, both being invariant under U (1) PQ . The self couplings in the latter are assumed to satisfy λ H , λ σ > 0 and λ 2 Hσ < λ H λ σ , to ensure that the minimum of the scalar potential is attained at the vacuum expectation values (vevs) where v = 246 GeV. The real scalar field a parameterizing the phase of the hidden scalar field σ in the expansion around the vev, then corresponds to the Nambu-Goldstone boson arising from the breaking of the global PQ symmetry. Assuming that the vev v PQ is much larger than the electroweak scale, one may integrate out the heavy fields ρ and Q. The low-energy effective Lagrangian of a matches Eq. (22), with n ax = 1, a 1 ≡ a , and The fermionic current associated with the PQ symmetry is indeed anomalous, leading, as follows from the general formulae (A10a) in Appendix A, to axionic couplings to gauge bosons in the low-energy effective Lagrangian (22), with C (Q) emi denoting a possible electric charge of the exotic color triplet. The coupling to the electron, on the other hand, vanishes in the KSVZ model, as follows from the general formula in Appendix A, B. DFSZ-type models Along the same lines one can also construct DFSZ-type ad-hoc models [6,7], which require, apart from the generic hidden complex scalar σ, at least two Higgs doublets, H u and H d , giving masses to up-type and down-type quarks, respectively, and in which the SM quarks, instead of an exotic vector-like color triplet, carry PQ charges. The fields are supposed to transform as follows under the PQ symmetry, leaving the Yukawa interactions, with H u = H * u , as well as the most general renormalizable scalar potential, invariant under the PQ symmetry, provided that After PQ and EW symmetry breaking at the scales v PQ v = v 2 u + v 2 d = 246 GeV, the Higgs and singlet scalar fields can be parametrized as The fields in the radial components, h d , h u , and ρ, correspond to the two physical neutral Higgs bosons and the physical scalar singlet, respectively. The phases of the fields involve the Nambu-Goldstone boson ζ which is eaten by the Z 0 to generate its mass, and the Nambu-Goldstone boson a which is an axion-like field. Orthogonality of ζ and a requires The condition (79) then determines that The low-energy effective Lagrangian of a matches Eq. (22), with n ax = 1, a 1 ≡ a , and the latter arising from the requirement of canonical field normalization of a using f PQ = v PQ . The fermionic content of the theory is anomalous with respect to U (1) PQ , and the couplings to SM gauge bosons can be obtained from the general formulae in Appendix A. The coupling to gluons is given by (cf. Eq. (A10a) in Appendix A) where N g is the number of generations. The couplings to photons (see for example [114]) and electrons (see Appendix A for the general formulae), on the other hand, depend also on the PQ charge assignment for the leptons. In fact, one finds for the coupling to the gluon. Furthermore, for the coupling to the electron, we get On account of Eq. (81), the physically most relevant ratio of the electron coupling relative to the gluon coupling, cf. Eq. (36), are then finally obtained as if only a third H l couples to l eR .
C. Models relating the PQ scale with the seesaw neutrino scale The Peccei-Quinn breaking scale in these models is however still ad-hoc. In fact, the strong CP problem is solved for any value of the Peccei-Quinn scale and the preferred range for it, 10 9 GeV f A 10 12 GeV, arises just from astrophysical (stellar bounds) and cosmological (dark matter) constraints. Intriguingly, however, this range overlaps with the preferred one for the breaking scale of lepton number in seesaw explanations of the smallness of the active neutrino masses. These scales can indeed be identified in an extension of the KSVZ-type or DFSZ-type models by the further inclusion of three right handed Majorana neutrinos N iR , i = 1, 2, 3 (for similar considerations, see Ref. [115,116]). Considering for example the KSVZ-type model and assigning to the right handed neutrinos, as well as to the left handed SM lepton doublets L and to the right handed SM lepton singlets l R , the same Peccei-Quinn charges as to the right-handed exotic color triplet (cf. Eq. (68)), the most general PQ invariant Yukawa interactions involving the fields charged under the Peccei-Quinn symmetry generalize to Here H denotes the Higgs doublet, H = H * , and G ij , F ij , Y ij , are arbitrary complex 3 × 3 matrices, while y ij is a symmetric 3 × 3 matrix.
In this case, the Peccei-Quinn symmetry is an extension of the global lepton number U (1) L in the SM (including now right-handed neutrinos), which is spontaneously broken when σ acquires a vev. Therefore, in this model, the Nambu-Goldstone boson a is in fact what is usually known as a majoron [53,54,115,116]. After integrating out the heavy fields, the low-energy couplings of a are still given by , and C a e = 0. (90) However, in this case, v PQ is at the same time the seesaw scale. In fact, below EWSB, the last two terms of Eq. (89) give rise to a neutrino mass matrix of the form realizing the seesaw mechanism [117][118][119], i.e. explaining the smallness of the masses of the left-handed SM active neutrinos by the large mass of the right-handed SM singlet neutrinos, Remarkably, for an intermediate scale symmetry breaking scale, f a = v PQ ∼ 10 11÷12 GeV, such a simple single KSVZ-like extension of the SM can explain simultaneously • the non-observation of strong CP violation (due to the coupling of a to GG), • the smallness of the active neutrino masses (seesaw with natural sizes of Yukawa couplings), • the nature of cold dark matter (axions), • the baryon asymmetry of the universe (due to thermal leptogenesis from N iR decay [120], requiring M N iR 10 9 GeV [121]), and • the stabilization of the electroweak vacuum (due to a threshold effect associated with the Higgs portal term proportional to λ Hσ in Eq. (70) [122]).
Clearly, these single axion models can be easily extended to have two (or more) Peccei-Quinn symmetries, realizing KSVZ × KSVZ, DFSZ × DFSZ, and mixed KSVZ × DFSZ models. We will present a particularly well motivated KSVZ × DFSZ type model in Sec. IV A.

D. Photophilic models
We have, however, not yet exhausted our model building possibilities. Both in KSVZ-type as well as in DFSZtype models, there are always colored fermions carrying PQ charges in a chiral manner -the exotic color triplets and ordinary SM quarks, respectively. Therefore, in these cases the axion-like fields always have a non-zero axionic coupling to GG, C a g = 0. In view of the astrophysical hints for the existence of an ALP different from the axion it is of considerable phenomenological interest to construct also models in which the axion-like field is photophilic, i.e. C a g = 0, but C a γ = 0, as has been discussed in Sec. II.
The simplest model predicting a photophilic ALP is a variant of the KSVZ model where the exotic fermion is a colorless but electrically charged particle, which we denote by E = (E L , E R ). The hidden complex singlet σ is introduced as usual and the PQ symmetry 5 acts as The Yukawa Lagrangian is analogous to the KSVZ model and reads The parametrization of the axion-like field a remains as in Eq. (72), with Eq. (73) setting the ALP decay constant f a . The couplings of a to gluons, photons, and electrons are given by (cf. Appendix A) where C (E) em denotes the electric charge of E. If, however, the interaction in Eq. (94) is the only source of interaction for E at the PQ scale or below, the E-number would be conserved and this exotic lepton would constitute a stable charged particle which is cosmologically problematic, unless its mass m E = y v PQ TeV (for a review, see Ref. [123]), requiring then an unnaturally small Yukawa coupling y 10 −6 if v PQ 10 9 GeV. One possible remedy for this situation is to set C (E) em = ±1, so that either E R or E c L has the SM gauge quantum numbers of l iR and then E is allowed to mix with the SM charged leptons and to decay to SM particles. A similar mixing mechanism has been proposed to allow for the decay of the exotic superheavy quark Q in the KSVZ model which shares with the superheavy exotic fermion E the cosmological problem [124][125][126]. We present models of this type with both an axion and an additional photophilic axion-like field in Secs. IV B 2 and IV B 3.
Another simple alternative photophilic model exploits, apart from the usual hidden complex scalar σ and similar to the DFSZ model, a second Higgs doublet, H l , which gives rise to the masses of the charged leptons, cf.
This extra Higgs doublet is necessary for the restriction of the U (1) PQ symmetry to the lepton sector. The fields are supposed to transform as σ → e iα σ , under the U (1) PQ symmetry, such that both the above Yukawa Lagrangian and the scalar potential, where where v l is the vev of H l , and the couplings to gluons, photons, and electrons are (cf. Appendix A) We present a well-motivated model of this type with both an axion and an additional photophilic axion-like field in Sec. IV B 1. These types of model have the characteristic feature that the presence of a photophilic axion-like field is always accompanied by possible signals at the electroweak scale through H l .

IV. INTERMEDIATE-SCALE ACCIDENTAL AXION AND FURTHER ALPS FROM FIELD THEORETIC BOTTOM-UP SM EXTENSIONS
The ad-hoc Peccei-Quinn extensions of the Standard Model discussed above give some insight into the origin of the decay constants and the sizes of the anomaly coefficients, but still have several drawbacks. Most importantly, new fields (hidden complex scalar fields, heavy vector-like quarks, extra SM Higgses) and new PQ symmetries have been introduced and imposed by hand. Moreover, these symmetries are not protected from explicit symmetry breaking effects by Planck-scale suppressed operators -e.g.
for a model with two singlets -expected to appear generically in the low-energy effective Lagrangian [127][128][129][130][131][132]. These operators modify the axion potential, eventually shifting its minimum away from zero, thereby destroying the solution of the strong CP problem, cf. Appendix B. Moreover, they induce masses, e.g.
lifting also the additional ALPs to pseudo Nambu-Goldstone bosons. Both drawbacks can be absent in models where the Peccei-Quinn symmetries are not ad-hoc, but instead automatic or accidental symmetries [127]. In fact, Peccei-Quinn symmetries could be accidental consequences of exact discrete Z N symmetries, which in addition, if N = D + 1 is large enough, i.e.
We now present a SM extension with two axion-like fields corresponding to (pseudo-) Nambu-Goldstone bosons from the breaking of two continuous Peccei-Quinn symmetries U (1) PQ1 × U (1) PQ2 , which themselves are accidental consequences of exact discrete Z N symmetries. The model is a further extension of an n ax = 1 model proposed in Ref. [133], where the SM was extended by several fields: not only by a SM complex scalar singlet σ, but also by several SM Higgs multiplets and three right-handed SM singlet neutrinos. These extra fields were primarily included to enable the formulation of a large enough discrete symmetry ensuring the protection of the QCD axion against de-stabilization by Planck-scale suppressed explicit symmetry-breaking effects. Among the attractive features of the model are • the close connection between the seesaw scale to explain the smallness of the active neutrino masses and the Peccei-Quinn scale (see also [56,115,116,[139][140][141][142][143][144][145][146][147]), • the stabilization of the accions and the proton, as well as the non-occurence of flavor changing neutral currents, due to the large discrete symmetries, and The extended n ax = 2 model has four SM Higgs doublets H u , H d , H l , H N , which give Dirac mass terms for up-type quarks, down-type quarks, charged leptons and neutrinos, respectively, see Eq. (104) below; an SU (2) L triplet T with hypercharge Y = 2; and two SM singlet complex scalar fields σ 1 and σ 2 , which carry only Peccei-Quinn charges. The SM fermionic content is augmented by the addition of a vector-like color triplet, (Q L , Q R ), as in the KSVZ model, and three right-handed neutrinos, N iR , i = 1, 2, 3.
We impose an exact discrete symmetry Z 13 ⊗Z 5 ⊗Z 5 , with the fields transforming according to Table I. As mentioned above, the reason for such a large discrete symmetry is to suppress operators up to some mass dimension, thereby ensuring the effectiveness of the Peccei-Quinn mechanism. The order of each group factor is chosen to be a prime number because it forbids dangerous operators formed by invariant powers of the singlets σ k 1,2 , with k less than the group order, as dictated by the Lagrange theorem.  The discrete symmetry in Table I allows for the following Yukawa interaction terms: Here q iL , i = 1, 2, 3, denote the left-handed SU (2) L quark doublets, H = H * , and Y ij , Γ ij , F ij and G ij are complex arbitrary 3 × 3 matrices, and y ij is a symmetric 3 × 3 matrix. Of phenomenological great importance is the fact that, despite the occurence of multiple Higgs doublets, no flavor changing neutral currents emerge, since there is only one Higgs doublet for each type of fermion.
The scalar potential with all renormalizable terms allowed by the discrete symmetry is V = V H + V N H , where the hermitian terms are where k, k = u, d, l, N , a = 1, 2, while the nonhermitian terms read Note, that if we did not impose the additional Z 5 ⊗ Z 5 symmetry, dangerous terms of the type ( , and σ 2 1 σ 2 , would still be allowed, breaking the accidental global symmetries discussed below, in particular the Peccei-Quinn symmetries. Accidentally, the model (104) -(106) possesses five U (1) (i.e. rephasing) symmetries which are constrained by the non-hermitean terms in (104) and (106): the imposed U (1) Y invariance and four additional accidental global U (1) symmetries. Among the latter, we can identify two non-chiral ones, namely U (1) B , ensuring baryon number conservation for usual quarks and U (1) Q , ensuring Q-number conservation for the exotic quark. The two remaining accidental symmetries, the Peccei-Quinn symmetries U (1) PQ1 and U (1) PQ2 , involve the singlets σ 1 and σ 2 , and are chiral for quarks.
The latter accidental symmetries can be chosen in a way that σ 1 is only charged under U (1) PQ1 , while σ 2 is only charged under U (1) PQ2 . For example, the potential in (106) is independent of σ 1 and so a change of the global phase of this field can be compensated by changing the phases of Q L,R , N iR , L i , and l iR . We identify this symmetry with U (1) PQ1 .
We denote the charges of a field ψ under the Peccei-Quinn symmetries U (1) PQ1 and U (1) PQ2 as X ψ and K ψ , respectively, If we write X σ1 = (X u + X d ), X u ≡ −X u R , X d ≡ −X d R , and choose the normalization then the transformation of the singlets reads and the theory (104) -(106) is obviously invariant under the following U (1) PQ1 transformation (displaying only the non-trivially transforming fields) This symmetry is an extension of the SM lepton number, including right-handed neutrinos, and is spontaneously broken when σ 1 acquires a vev, cf. Sec. III C. We show the U (1) PQ1 charges also in Table II.
The U (1) PQ2 charges are also easily extracted by requiring invariance of (104) -(106) under U (1) PQ2 , and we find It should be noted that we have written some charges in terms of X u , X d , X L , X l but there is only a freedom to choose X u or X d , which is subsequently fixed after electroweak symmetry breaking. The rest of the charges can be found in the second line of Table II and they are equivalent to the ones in Ref. [133] when the nonhermitian terms in the potential are the same. When defining the U (1) PQ2 charges in Eq. (111), we have chosen σ 1 and q L to be uncharged under U (1) PQ2 , i.e., X σ1 = X q L = 0, and X Q L = −X Q R , for convenience. The first choice is what turns N iR and Q L,R uncharged under U (1) PQ2 . A different choice would lead to a different value for the coefficients C ig and C iγ , which are dependent on the choice of U (1) PQ1 and U (1) PQ2 . The couplings to the physical axion and ALP, however, are fixed for a given model and thus independent of the choice of U (1) PQ1 and U (1) PQ2 .
The parameters in the scalar potential are assumed to be such that the SM doublet Higgs fields get non-zero vevs GeV, while the SM singlet fields are assumed to attain their vevs, σ a = v a / √ 2, a = 1, 2, at much higher scales. Therefore, the right-handed neutrinos N iR and the extra colored triplet, which get their masses via Yukawa couplings with σ 1 in Eq. (104), are heavy, decoupling from the electroweak scale. The first two terms in the second line of Eq. (104) give rise to a neutrino mass seesaw relation, which now involves the vev v N of H N , which might be much smaller than v, At low energies, two Nambu-Goldstone bosons a a with decay constants arise from the breaking of the two accidental accidental Peccei-Quinn symmetries, describing the phases of the singlet fields σ b in the expansion around the vevs. The Nambu-Goldstone boson a 1 is thus also a majoron-like particle [53,54]. The interaction term y Q Q L σ 1 Q R implies that U (1) PQ1 is in fact a chiral symmetry for the exotic quark Q. As noted previously, the number of this quark is separately conserved.  The ALP-gluon and ALP-photon couplings can be extracted from the charges in Table II as (cf. Appendix A) leading to The ALP-electron couplings can be extracted in the same way (cf. Appendix A), so that where tan β ≡ v u /v d . Therefore, the considered model automatically leads to a low energy theory with two axion-like fields a b , both coupling to GG and FF . Successful implementation of the neutrino seesaw mechanism with not too much fine-tuning in the Yukawa couplings and thermal leptogenesis prefer the symmetry breaking scale v 1 and thus the decay constant f a 1 of the majoron-like particle a 1 to be in the range 10 9 GeV v 1 f a 1 10 13 GeV.
The symmetry breaking scale v 2 and thus the decay constant f a 2 of the axion-like field a 2 , however, is not fixed by bottom up considerations.
Remarkably, the upper end of the above parameter range is distinguished by another observation, namely gauge coupling constant unification. In fact, if we take both Peccei-Quinn symmetry breaking scales see Fig. 3. This value is surprisingly close to the PQ scales and indicates that the model can be embedded in a grand unification group, so pointing to a connection between such scales [149]. In this connection it is important to note that dangerous operators leading to baryon number violation processes [150] are forbidden by the discrete symmetry Z 13 ⊗ Z 5 ⊗ Z 5 to very high mass dimensions. Therefore, in this model proton decay is not observable, despite the low unification scale. Importantly, as detailed in Appendix B, the discrete symmetry Z 13 ⊗ Z 5 ⊗ Z 5 is large enough to protect the axion from destabilization by Planck scale suppressed explicit PQ symmetry breaking operators and that the corresponding induced masses are negligible, as long as This is based on the fact, that the operators of lowest mass-dimension that are invariant under the discrete symmetries of Table I, but not invariant under the Peccei-Quinn symmetries of Table II in the low-energy effective Lagrangian involve high powers of the scalar singlet fields, The theoretically most favored value for the decay constant of the axion in this model is then (cf. Sec. II) leading to the following mass predictions for the axion A and the ALP a (cf. Sec. II and Appendix B), Therefore, this model favors axion dark matter in the mass range probed by ADMX [37]. The ultralight ALP a in this model, however, will escape detection in current and near future laboratory experiments, such as ALPS-II [98] and IAXO [99], and can not explain the astrophysical puzzles such as the anomalous cosmic γ-ray transparency, the cooling excess of stars and the soft X-ray excess from Coma (cf. Sec. I). In fact, the axion and the ALP consist here of a non-negligible mixture of the original axion-like fields, with mixing cos δ 0.083, sin δ 0.996, tan δ 12.
Correspondingly, the ALP couplings to photons and electrons are of quite similar size as the respective axion couplings (for definiteness, we chose here C (Q) em = 1), and thus way too small to explain the cosmic γ-ray transparency (the latter requiring |g aγ | 10 −12 GeV −1 , cf. Eq. (15) and Fig. 1) and the excess star cooling (the latter requiring |g Ae | 10 −13 and/or |g ae | 10 −13 , cf. Eq. (20)). Moreover, the expected cosmic mass fraction in ALP dark matter is negligible [67], making also the direct laboratory detection of dark matter comprised of the ALP a impossible. Nevertheless, the indirect detection of an ultra-light ALP which such a small coupling to photons may still be possible with the next generation of cosmic microwave background (CMB) observatories such as PIXIE and PRISM (cf. Fig. 4) because resonant photon-ALP oscillations in primordial magnetic fields may lead to observable spectral distortions of the CMB [94][95][96].  On the other hand, if we give up unification and allow for Yukawa couplings of order 10 −4 to obtain the small active neutrino masses even in the case of a lower seesaw scale, we can get an ALP in the parameter range favored by the explanation of the soft X-ray excess from the Coma cluster. Choosing for example v 1 f a 1 = 1 × 10 10 GeV, v 2 f a 2 = 7.5 × 10 10 GeV, leads to an axion A with a quite small decay constant and thus rather large mass, constituted mainly by the axion-like field a 1 (small mixing δ), cos δ 0.96, sin δ 0.29, tan δ 0.3, with a quite sizeable coupling to the photon, At the same time, the ALP a, with mass m a 1.8 × 10 −33 eV, is constituted mainly by the axion-like field a 2 , and has also a sizeable coupling to the photon, right in the range of phenomenological interest. The couplings to the electron, however, are in this case still too small, in order to explain the anomalous cooling of stars. For illustration, let us take also a third pair of benchmark values for the PQ scales which leads to a case with very tiny mixing between the axion an the ALP. For definiteness, we chose v 1 f a 1 = 10 13 GeV, v 2 f a 2 = 10 10 GeV, corresponding to f A 3 × 10 9 GeV, m A 2 meV, m a 10 −28 eV, cos δ 3 × 10 −4 , sin δ 1, tan δ 3 × 10 3 , (140) with couplings of size In summary, cf. Table III and Fig. 4, we have now verified in an explicit example what has been anticipated in the general discussion of Sec. II: in the case that the axion A and the ALP a consist of an appreciable mixture of the original axion-like fields, i.e. as long as | tan δ| is of order unity, the couplings are all determined by f A and are therefore expected to be approximately of the same size for for the axion A and the ALP a, up to order one factors. This is clearly the case for models A1 and A2 in Table III. A hierarchical difference between the axion and the ALP couplings can possibly only arrive in the situation where the axion (and correspondingly also the ALP) originates essentially only from one axion-like field, which requires | tan δ| ≈ 0 or | tan δ| 1. In this case, the (photon and electron) couplings of the ALP are suppressed by the factor tan δ/(1 + tan 2 δ) in comparison to the couplings of the axion and thus hierarchically smaller than the latter.

B. Zn ⊗ Zm models with a photophilic ALP
As has been emphasized in Sec. II, in order to arrive at the same time at the prediction of a sizeable amount of axion dark matter and at the explanation of the astrophysical hints by an additional ALP one needs • that only one of the original axion-like fields couples to GG, while the other is photophilic, and • that there is a hierarchy between the vevs ∼ decay constants of the two PQ symmetries.
Therefore, in this subsection we build models featuring one accidental KSVZ-type axion and another accidental photophilic ALP, which does not couple to GG through a U (1) PQ2 ×SU (3) C ×SU (3) C chiral anomaly, but, nevertheless, couples to FF due to a U (1) PQ2 × U (1) em × U (1) em anomaly. The latter will be introduced as a photophilic analog of the DFSZ or KSVZ axion, cf. Sec. III D.

Z11 ⊗ Z9 model with extended Higgs sector
Our first model has two SM Higgs doublets, H q and H l , where the first one gives Dirac mass terms for quarks and neutrinos, while the second one gives them for charged leptons, see Eq. (145) below. Furthermore, it features two SM singlet complex scalar fields σ 1 and σ 2 , which carry only Peccei-Quinn charges. The phase of σ 1 will carry a KSVZ-like axion, while the one of σ 2 will carry a DFSZ-like photophilic ALP. The SM fermionic content is augmented by the addition of a vector-like color triplet, (Q L , Q R ), as in the KSVZ model, and three right-handed neutrinos, N iR , i = 1, 2, 3.
qL uR dR L lR NR Hq H l QL QR σ1 σ2  Given this field content, we impose an exact discrete Z 11 ⊗ Z 9 symmetry, with the fields transforming according to Table IV. Then the most general Yukawa interactions are given by The scalar potential V = V H + V N H is supposed to be given by the sum of all possible renormalizable hermitian terms involving the scalar fields and respecting the discrete symmetry in Table IV, V H (H q , H l , σ 1 , σ 2 ), and the most general renormalizable non-hermitian term ("DFSZ-type" term),  The Yukawa interactions (145) feature an accidental U (1) PQ1 ⊗ U (1) PQ2 symmetry which is summarized in Table  V. In fact, one PQ symmetry, U (1) PQ1 , may be defined by the following action on the fields: while the rest of the fields transform trivially. This symmetry is anomalous with respect to QCD and gives thus rise to the axion, whose corresponding field appears in the phase of σ 1 in the expansion around its vev v 1 , with The other accidental PQ symmetry U (1) PQ2 acts on the fields as The PQ charges X l and X q are restricted by (146) to obey This symmetry is non-anomalous with respect to QCD, but anomalous with respect to QED. Therefore, strictly speaking, U (1) PQ2 is a PQ-like symmetry, rather than a proper PQ symmetry. Its breaking will lead to a photophilic axion-like field a 2 appearing in the phase σ 2 in the expansion around its vev v 2 , where The low-energy couplings to the gluon are found to be (cf. Appendix A) Therefore, a 1 can be identified with the axion field A(x), while a 2 is a photophilic axion-like field,. The couplings to the photon, on the other hand, read where C (Q) em is the the electric charge of Q L , Q R . The coupling with electrons can be also extracted from Appendix A as so that where here tan β = x = v q /v l . As in our GUT example from Sec. IV A, U (1) PQ1 is in fact a lepton number symmetry and its breaking will thus give rise to a seesaw neutrino mass relation, where v q = √ 2 H k is naturally of order the electroweak scale, since v = v 2 q + v 2 l = 246 GeV. For definiteness and in order to get an axion well in the mass range favored after BICEP2 (cf. Sec. I), we allow for a modest fine-tuning in the Yukawa's and chose the vev v 1 to be v 1 = 3 × 10 11 GeV, (159) leading to an axion decay constant, mass, and photon coupling of order (we take again C (Q) em = 1) right in the ballpark to explain and detect dark matter in terms of axions. The ALP couplings, on the other hand, can be simultaneously in the range of highest phenomenological interest, |g aγ | ∼ 10 −12 GeV −1 , |g ae | ∼ 10 −13 , if we chose the vev v 2 of the second PQ symmetry much smaller, v 2 ∼ 10 9 GeV. (161) In fact, then In this context it is also important to note that the photophilic ALP has a quite sizeable mass, for the chosen discrete symmetry and PQ symmetry breaking scale. In fact, the lowest mass dimension operator violating U (1) PQ2 , but respecting the discrete symmetry, occurs already at mass dimension D = 9, and therefore induces a mass which does not suffer from too much Planck-suppression (cf. Appendix B), Finally, the axion in this model is protected from destabilization by Planck suppressed operators. In fact, the lowest mass dimension operator violating U (1) PQ1 , but respecting the discrete symmetry, has no destructive effect on the axion properties, as long as (v 11 1 v 2 ) 1/12 10 12 GeV (cf. Appendix B), which is the case for the favored choice of parameters.
We have summarized the values predictions of masses and couplings of this model for the chosen input values of the PQ breaking scales in Table III and Fig. 4, labelled as "Model B1.1". We see that they are well in a region to explain the astrophysical puzzles and to be probed in the next generation of experiments.
Interestingly, for another set of initial input values, which we label as "Model B1.2" in Table III cf. Fig. 4. In this case, it may in fact explain the recently reported unidentified 3.55 keV line from galaxies and galaxy clusters [105,106], mentioned in Sec. I, in terms of two photon decay [107][108][109][110], if its fractional matter density x a ≡ ρ a /ρ CDM is of order x a 10 −8 − 10 −7 .
2. Z11 ⊗ Z9 model with minimal SM Higgs sector As a further example, we consider a model which features, like the KSVZ model, a minimal SM Higgs sector and where all fields beyond the SM are SU (2) L singlets. These are: two scalar singlet fields σ 1 , σ 2 , an exotic quark singlet field Q, a colorless and electrically charged fermionic field E, and three types of right-handed neutrinos, N iR . The imposed discrete Z 11 ⊗ Z 9 symmetry transformations are given in Table VI. qL H uR dR L lR NR σ1 QL QR σ2 EL ER  The Z 11 ⊗ Z 9 symmetry allows the following Yukawa interaction terms In the scalar potential all interaction terms are hermitian. Eq. (169) has an automatic U (1) PQ1 ⊗ U (1) PQ2 symmetry in which the fields carry charges as shown in Table VII. Also, there is a seesaw mechanism for the neutrinos working like in the previous models. The axion-like fields in this case appear only in the singlets fields parametrization of Eq. (72), such that The low-energy couplings to the gluon they are found to be (cf. Appendix A) rendering a 2 a photophilic axion-like field, while a 1 can be identified with the axion field A(x). The couplings to the photon, on the other hand, read where C (Q) em = −1/3 and C (E) em = −1 are the electric charge of Q and E, respectively. The couplings to the electron can be also extracted from Appendix A as such that Again, the active neutrino masses favor v 1 10 12÷14 GeV, while the explanation of the astrophysical anomalies favor v 2 ∼ 10 8÷10 GeV. As in the previous model, the lowest mass dimension operator violating U (1) PQ2 , but respecting the discrete symmetry, is of the form (183), inducing an ALP mass of size m a 1 × 10 −7 eV |g| 1/2 v 2 10 9 GeV Moreover, as in the previous model, the axion is protected from destabilization by Planck suppressed operators (cf. Appendix B). The masses of E and Q are given by We see in this scenario E could naturally have a mass around the TeV scale or below. By means of a mixing proportional to v 2 /M Pl ≈ 10 −10 we may have a decay of E into SM charged leptons, l i , plus the Higgs boson, h, according to E → l i + h. This could leave a signal in Drell-Yan E pair production at a hadron collider such as LHC. For Q, its mass would be out of the present direct searches unless there is a fine tuning of order y Q ≈ 10 −4 .
Assuming the latter, a QCD pair production process of Q, followed by decays into SM d-type quarks plus a Higgs boson, Q → d i + h, is possible at the LHC. 3. Z11 ⊗ Z7 model with minimal SM Higgs sector Motivated by the possibility to explain the 3.55 keV line by decaying ALP dark matter, we introduce yet another model with a somewhat smaller discrete symmetry to get a keV scale ALP mass with a PQ symmetry breaking scale of order 10 9 GeV, instead of the ∼ 10 12 GeV required in the model of Sec. IV B 1. It is based on the fact that, with the same field content as in the model from Sec. IV B 2, we can also realize a Z 11 ⊗ Z 7 symmetry as given in Table  VIII. qL H uR dR L lR NR σ1 QL QR σ2 EL ER  In this case the lowest dimensional operator breaking U (1) PQ2 and involving just σ 2 is and this induces the following mass for the ALP (cf. Appendix B), Therefore, an explanation for the 3.55 keV line from galaxies can be found with this setup of fields and the Z 11 ⊗ Z 7 symmetry, cf. Table III and Fig. 4, labelled "Model B3". The lowest dimensional operator involving σ 1 and breaking U (1) PQ2 is the following one Therefore, there are no dangerous Planck suppressed operators destabilizing the axion for v 1 = 10 13 GeV. Also, the mass generation mechanisms for neutrinos, Q and E keep the same as before.

C. Z11⊗Z9⊗Z7 model with two photophilic ALPs
Finally, we construct a model with two photophilic ALPs: one very light, to explain the astrophysical hints such as the cosmic gamma ray transparency, and one with a mass of 7.1 keV, to explain the recently found 3.55 keV line from Andromeda and galaxy clusters by its decay in two photons. Clearly, such a model will involve three complex scalar singlet fields σ i , i = 1, 2, 3. Also, besides an exotic quark singlet field Q, two colorless and electrically charged fermionic fields, E, E , and three right-handed neutrinos, N iR , are introduced.  The field content of our model is specified in Table IX, which displays also the field transformations which we impose to ensure an exact discrete Z 11 ⊗Z 9 ⊗Z 7 symmetry. Based on this, the most general interactions involving scalars and fermionic fields are given by This interaction Lagrangian features an accidental U (1) PQ1 ⊗ U (1) PQ2 ⊗ U (1) PQ3 symmetry shown in Table X. The lowest mass dimension operator violating U (1) PQ2 , but respecting the discrete symmetry, occurs already at mass dimension D = 9, Thus, the mass of the respective ALP a 2 will be the same as in Eq. (165), The lowest mass dimension operator violating U (1) PQ3 , but respecting the discrete symmetry, occurs even at lower mass dimension, D = 7, leading to a mass (cf. Eq. (178)) m a3 7.1 keV |g | 1/2 v 3 1.84 × 10 9 GeV The mass of Q is the same as in Eq. (176), and the mass matrix involving E and E on the basis {E L , E L } is whose eigenvalues are where These masses can be naturally in the 10 2 -10 3 GeV scale without severe fine tuning. As an example, taking k 11 = 0.5, k 12 = 0.8, k 13 = 0.8, k 23 = 0.5, and v 1 = 10 13 GeV, v 2 = 10 9 GeV, v 3 = 1.8 × 10 9 GeV, we obtain M + ≈ 932 GeV and M − ≈ 188 GeV. We see that, as in the model of Sec. IV B 2, the exotic fermions in this model also have masses eventually in reach of colliders.
The axion in this model is protected from destabilization by Planck suppressed operators. In fact, the lowest mass dimension operator violating U (1) PQ1 , but respecting the discrete symmetry, has no destructive effect on the axion properties, as long as (v 11 1 v 2 v 3 3 ) 1/15 10 12 GeV (cf. Appendix B), which is marginally the case for the favored choice of parameters.  The low-energy couplings to the gluon they are found to be (cf. Appendix A) The couplings to the photon, on the other hand, read where C For v 1 = 9 × 10 11 GeV , v 2 = 4 × 10 9 GeV , v 3 = 3 × 10 9 GeV , the axion A = a 1 has a mass in the regime where it could be the dominant part of dark matter, and the two photophilic ALPs, a 2 and a 3 have masses and couplings in the regime of interest of the astrophysical hints (gamma ray transparency; X-ray line), see Table XI. In fact, as can be seen in Fig. 4, it is easy to accommodate all the hints in our model.

V. SUMMARY AND OUTLOOK
Motivated by theory, cosmology, astrophysics, and upcoming terrestrial experiments, we have undertaken a general analysis of models featuring the axion A plus a further ALP a. In these models, the former may constitute dark matter, while the latter may explain simultaneously the anomalous cosmic γ-ray transparency, the anomalous cooling of white dwarfs and red giants, and the soft X-ray excess from the Coma cluster (in case that there is a cosmic ALP background radiation peaked in the soft-keV regime). To this end, the ALP-photon coupling is required to be of order |g aγ | 10 −12 GeV −1 , for sufficiently small ALP mass, cf. Sec. I and Fig. 1. Fortunately, the next generation of light-shining-through-a-wall experiments (ALPS-II) and helioscopes (IAXO) will be sensitive in this range of ALP parameter space. Intriguingly, a further ALP, decaying into two photons with a similar coupling strength, but with a much larger mass of 7.1 keV, may explain the recently reported 3.55 keV line from Andromeda and galaxy clusters if it constitutes a tiny part of cold dark matter, cf. Sec. I and Fig. 1.
In models with several axion-like fields, a i , i = 1, .., n ax , coupling to the topological charge density of QCD, the field A corresponding to the proper QCD axion is a mixture, cf. Sec. II, In general, if the axion-like field associated to the lowest PQ scale couples to all gluons, photons and electrons, the ALP couplings to photons and electrons can not be hierarchically larger than the respective axion couplings (although they can be hierarchically smaller). In fact, the couplings satisfy the constraints, , with 1 which generalize Eqs. (56) and (67) in Sec. II to more than two axion-like fields. Barring accidental cancellations and considering order one coefficients, we have typically f A ∼ f γ ∼ f e , |g aiγ | |g Aγ | ∼ α/(2πf A ), and |g aie | |g Ae | ∼ m e /f A . This has important consequences for phenomenology. In particular, in models featuring mixing, the discovery of an ALP in the next generation of laboratory experiments would imply f A ∼ α/(2π|g aγ |) ∼ 10 9 GeV and thus an axion mass in the m A ∼ meV range -unfavorably large for axion dark matter, cf. Fig. 1. Conversely, the detection of axion dark matter with a mass in the favored (0.1 -100) µeV range would imply, in these models, a tiny ALP-photon coupling, making the latter unaccessible to the next generation laboratory experiments. Micro-eV mass axion dark matter and a large ALP-photon coupling, |g aγ | 10 −12 GeV −1 , can occur simultaneously only in models where the axion-like field dominantly coupled to photons (lowest scale) does not couple to gluons, i.e., it is photophilic (e.g. C 1g = 0, C 2g = 0, C 2γ = 0).
We have reviewed in Sec. III plausible ultraviolet completions of the Standard Model in which axions and more general ALPs occur as Nambu-Goldstone bosons from the breaking of U (1) Peccei-Quinn (PQ) symmetries. Parameters occurring in the low-energy effective Lagrangian (decay constants f a and dimensionless couplings C ij ) were then determined in terms of the fundamental parameters of the underlying high-scale theories (PQ symmetry breaking scales v PQ and PQ charges). We have emphasized two especially well motivated classes of models: i) Models, in which the PQ symmetry coincides with the lepton number symmetry and thus the PQ breaking scale with the seesaw scale, favoring therefore a decay constant of order f a ∼ v PQ ∼ 10 11÷15 GeV, to explain the active neutrino masses, m ν ∼ v 2 /v PQ ∼ 0.01 eV (the spread arises from our ignorance of the relevant Yukawa couplings). Intriguingly, such a simple, minimalistic extension of the SM can explain at the same time and with the same PQ breaking scale also the nature of dark matter (axions), the baryon asymmetry of the universe (through leptogenesis) and the stability of the electroweak vacuum, as has been outlined in Sec. III C. ii) Models, which feature a photophilic axion-like field which then possibly has a much larger coupling to photons than the axion.
The high-scale models introduced in Sec. III can be combined to construct low-energy theories featuring simultaneously an axion and one or more further ALPs. However, such ad-hoc models in which the new PQ symmetries are introduced and imposed by hand, have the drawback that these symmetries are not protected from explicit symmetry breaking by Planck-scale suppressed operators, which possibly spoil the solution of the strong CP problem. Therefore, in Sec. IV, we have considered several classes of models in which the Peccei-Quinn symmetries are not ad-hoc, but instead accidental symmetries originating from large exact discrete symmetries: Z 13 ⊗ Z 5 ⊗ Z 5 , Z 11 ⊗ Z 9 , Z 11 ⊗ Z 7 , and Z 11 ⊗ Z 9 ⊗ Z 7 , respectively. In the first three classes, the SM field content was enlarged by two hidden complex scalars, by extra Higgs doublets, and an exotic colored fermion; in the fourth class there was an additional third hidden complex scalar introduced. In the considered models, both the right amount of axion dark matter, viable neutrino masses, and the baryon asymmetry of the universe can be obtained by properly choosing the vev v PQ 1 f a 1 ∼ 10 11÷13 GeV of one of the hidden scalars. Only in the first, Z 13 ⊗ Z 5 ⊗ Z 5 , model the two axion-like fields turn out to mix (both C 1g and C 2g are non-vanishing), while in the other model classes, the remaining axion-like fields are photophilic (C ig = 0, C iγ = 0, for i = 2, 3). Remarkably, this model predicts the unification of the SM gauge couplings at M U ≈ 2.8 × 10 13 GeV, while proton decay is very heavily suppressed due to the large discrete symmetry. For all classes of models we have given benchmark parameter values for the vevs of the hidden complex scalars and determined the corresponding low-energy parameters, exploiting the general expressions obtained in Sec. II. In particular, we have determined the ALP mass(es) induced by the explicit symmetry breaking operators. In the Z 13 ⊗ Z 5 ⊗ Z 5 model, it is very small because it arises from a term in the Lagrangian suffering a suppression by a rather high power of the Planck mass, while, in the Z 11 ⊗ Z 9 and Z 11 ⊗ Z 7 models, it is much larger, and m a ∼ 7 respectively, even for smaller PQ vevs. The Z 11 ⊗ Z 9 ⊗ Z 7 model has two ALPs: an ultralight ALP whose mass is given by (197) and which may explain e.g. the cosmic γ-transparency, and a much heavier ALP with mass according to Eq. (198) which may explain the 3.55 keV line. These results have been summarized in Tables III and XI and in Fig. 4. For the chosen benchmark values of the PQ breaking scales, these models indeed predict an axion and ALP(s) in the astrophysically, cosmologically, and experimentally favored parameter regions. We leave the study of possible tests of this class of models by looking for signatures of the extra -from the assumed field content beyond the SM -particles at the LHC for future work. A detailed study of the early cosmology of these models, in particular for the case where the reheating temperature of inflation is above the PQ phase transition, was beyond the scope of the present investigation, but is certainly of high interest. In this context, it is important to note that the models considered in Sec. IV do not suffer from the cosmological domain wall problem, as is shown in Appendix C. This is crucial in view of the fact that the recent discovery of the tensor modes in the B-mode power spectrum by BICEP2 strongly disfavor scenarios where the PQ phase transition occurs before inflation, depriving the universe from the possibility to wipe out topological defects created at the phase transition, such as domain walls, by inflation.
There are still ad-hoc features in these scenarios, however. In particular, the extension of the SM field content and the choice of the discrete symmetry group are mostly motivated from low-energy phenomenology. Intriguingly, it appears that large discrete symmetries able to give rise to multiple accidental Peccei-Quinn symmetries are a generic feature of top-down motivated heterotic string scenarios. In fact, orbifold compactifications of the heterotic string are known to predict a plenitude of hidden complex scalars and vector-like exotic particles and a multitude of discrete symmetries (R symmetries from the broken SO(6) symmetry of the compactified space and stringy symmetries from the joining and splitting of strings) which are exact at the perturbative level. We plan to scan the mini-landscape of certain heterotic orbifolds for models with accidental axions and further ALPs along the lines pioneered in Ref. [65,66]. Here, the phenomenologically most pressing question is whether there exist models where the vevs of the accidental PQ symmetries are naturally in the range of the intermediate scale, ∼ 10 10÷13 GeV, rather than in the range of the heterotic string scale, M s ∼ M Pl .
In models involving SM quarks, leptons, and right handed singlet neutrinos, the axion-like fields appear as phases of SM Higgs doublets H k in the Yukawa terms Multiple PQ symmetries can be treated analogously. The PQ charges X H k and X σ = 1 of the scalar fields fix the parametrization of the axion-like field a as This parametrization induces couplings between axion-like fields and fermions in Eqs. (A2) and (A3), which can be removed by the chiral rotations for α = a (x)/f a . After these chiral rotations we obtain the second and third term of Eq. (22) through Eq. (A1) with coefficients Here, the first terms on the right hand sides are due to the exotic quark Q, while the second terms are due to the SM quarks and leptons. The coupling of a with electrons come from the kinetic terml R iγ µ ∂ µ l R , after the chiral rotation (A6), which induces We have discarded a total derivative in the last equivalence. Therefore, we obtain the coefficient in Eq. (22). For example, for the DFSZ axion, we have X H d = X d , X Hu = −X u and X H l = X d , −X u , 0 for the three types of model in Eq. (85), respectively. All the supplied formulae are written in terms of the PQ charges of the scalars. We can rewrite them in terms of the PQ charges of fermions as follows: Formula (A10a) assumes the summation of all colored fermions ψ i which are in the fundamental representation of SU (3) c . Analogously, the summation in Eq. (A10b) goes over all electrically charged fermions such that quarks should be counted for each color. The formulae in (A10) have the advantage that they are independent of the specific implementation of the model (e.g. Yukawa terms) but depend only on the PQ charges of the fermions. We emphasize that the PQ charges are defined by relations such as (A4) and we are assuming the parametrization (A5). The decay constant f a is fixed by canonical normalization of the kinetic term (∂a ) 2 as in Eq. (22) so that rescaling of PQ charges also leads to rescaling of f a , keeping the terms in Eq. (22) the same. We have normalized the PQ charges such that the smallest of the PQ charges of the singlet scalars is set to X σ = 1, which usually leads to f a = √ 2 | σ | = v σ . This normalization is also relevant for the domain wall problem; see Appendix C.
In more general situations where we have colored fermions in representations other than the fundamental one, Eq. (A10a) generalizes to where b is a fixed index without summation. The SU(3) c generators should be normalized such that Tr[T a T b ] = I(R) δ ab , with I(R) = 1 2 for the fundamental representation.

Appendix B: Effects from Gravity through Planck-Scale Suppressed Operators
At this point we want to estimate the impact of high dimensional operators suppressed by the Planck scale, M Pl , on the Peccei-Quinn solution for the strong CP problem. We generalize the arguments in Refs. [128][129][130][131][132] for the case involving more Higgs fields, like the models we worked out here. We consider a generic multi-Higgs model with two Peccei-Quinn symmetries as in Sec. IV. We assume that the triplet T does not have a vev so that it will not contribute to the axion and ALP low-energy effective potential. A vev for T which is much lower than those for the doublets, i.e., T could also be considered without impacting significantly such effective potential, but we will not take it into account. Thus, we will only consider Planck-scale suppressed interactions involving Higgs doublets and singlets. Let us assume then that the operator with the lowest dimension that is suppressed by M Pl and breaks the Peccei-Quinn symmetries gives a contribution to the high-scale effective Lagrangian according to where b, c = u, d, l, N . The integer D = 2m + n + k > 4 is the operator dimension, and g = |g|e i∆ is the effective dimensionless coupling supposedly produced by gravitational interactions, which may violate CP (∆ = 0). The operator in Eq. (B1) contributes to the effective potential for the axion and the ALP as follows where Λ 4 (D) = |g| and The minus sign in ±n (±k) denotes σ * 1 (σ * 2 ) instead of σ 1 (σ 2 ) in Eq. (B1), K b,c (X b,c ) is the PQ charge of H b,c associated to the PQ symmetry for which σ 1 (σ 2 ) is charged and δ is the mixing angle introduced in Eq. (31). The first term in Eq. (B2) is the QCD instanton potential whereas the second terms arises from the dominant nonrenormalizable operator, Eq. (B1), from gravity effects breaking the Peccei-Quinn symmetries. We should check if these effects do not spoil the solution to the strong CP problem by destabilizing the effective potential of Eq. (B2). Experimental bounds on the strong CP violation parameter require A /f A < 10 −10 . As the minimum of V eff in the ALP direction requires this implies that the minimum of V eff still occurs for A = 0. We should note that if there is more than one operator below a "critical" dimensionality, the above observation may not be valid anymore. In this case there is not a single condition like Eq. (B5) and some symmetry suppression must take place. The critical dimensionality is obtained by considering that there is no significant destabilization of V eff due to operators like Eq. (B1) that break the Peccei-Quinn symmetries. This is ensured, if we require where we take sin C m,n,k to define the lowest dimensional operators that could be allowed. The vevs of the Higgs doublets can be naturally assumed as v b ≈ 10 2 GeV. So, taking We have verified that the model we have presented in Sec. IV A is safe from gravitational effects since the dangerous operators defined above are all forbidden by the discrete symmetries. The operators of the form (B1) with lowest dimension and m = 1, allowed by the symmetries of Table I, are except for the term H † d H u σ 2 in (106) which is invariant under U (1) PQ2 (and U (1) PQ1 ). We then seek all operators of dimension equal or lower than 15 with larger m. For m = 2, 4, 5 there are no operators besides (H † d H u σ 2 ) n , which we discard. For m = 3, 6, we can find the operators: We have again dropped the U (1) PQ1 invariant terms (H † d H u σ 2 ) n . We can see that the operators in (B10) are harmless because they are also invariant under the Peccei-Quinn symmetries. Therefore, the operators in (B9) are the lowest order operators that are invariant under the discrete symmetries of Table I but not invariant under the Peccei-Quinn  symmetries. An estimate for the mass corrections induced by gravity for both the axion and the ALP can be obtained from V eff in Eq. (B2). The leading operator breaking the Peccei-Quinn symmetries is the first one, O 14 , in Eq. (B9). Such operator leads to the following contribution to the mass matrix of the axion A and ALP a: where C A,a ≡ C 1,7,5 A,a , defined in Eq. (B3), are not expected to be much greater than unity. Thus, the correction δm A for the axion mass, and the generated ALP mass, m a , will be of order for the same set of values for the scales used above. Hence, no significant mass correction arises for both the axion and the ALP. The mixing is also quite small and can be neglected. Therefore, although the ALP gain a mass from gravitational interactions, turning it into a pseudo-Goldstone boson, its feeble mass highly suppresses any physical effect such as the ALP decay into photons. For models containing a photophilic ALP, the contribution from operators of the form (B1) to the effective potential (B2) should be analyzed differently and needs to be separated into two parts: (a) the axion potential and (b) the photophilic ALP mass. The contribution for (a) should be still negligible and we seek here a photophilic ALP mass satisfying the constraints (15). Note that in this kind of models the two scales associated to the axion and to the photophilic ALP are not related in general.
Let us assume σ 1 carries the axion, A = a 1 , whereas σ 2 carries the photophilic ALP, a = a 2 . To avoid a substantial A-a mixing, we need to ensure that the mixing terms are negligible compared to the dominant contributions giving mass to the axion (QCD potential) and to the photophilic ALP. The latter should be given by an operator (B1) with n = 0. Taking also m = 0, for simplicity, the operator has the form We can see their contribution to the QCD potential or to the photophilic ALP mass are negligible. For example, the first two operators in the first line of Eq. (B20) give, respectively, δm a ∼ |g| 1/2 10 −14 eV f a 2 10 9 GeV 5/2 , δm A ∼ |g| 1/2 10 −20.5 eV f a 1 10 12 GeV Other operators give even smaller contributions.

Appendix C: Domain-Wall Problem in the Models
Although the PQ symmetry is broken by QCD instantons, the vacuum state may still be symmetric under a discrete group Z N ⊂ U (1) PQ . As a result a vacuum degeneracy leading to topological defects such as domain wall might occur causing cosmological problems for certain axion models [153].
It was observed long ago that in some theories with spontaneous symmetry breakdown a domain structure for the vacuum is expected [154]. The reason is that when there are degenerated vacuum states the fields might settle down in different low energy configurations in causally disconnected regions, as the Universe temperature decreases below the energy scale where a symmetry is broken. In the case of axion models this symmetry is the U (1) P Q . Thus, domains with distinct vacuum configurations can be form in different regions of the Universe. The boundary between two domains is a domain wall, and it is a solution of minimal energy field configurations interpolating the neighboring vacua (for a review see [155,156]). Stable domain walls bring a cosmological problem because their contribution to the total energy density would be too large for a flat, homogeneous isotropic Universe [154], [155], [22]. Many studies and solutions of the domain wall problem in axion models were performed [153,[157][158][159][160][161]. One of the simplest solution is to ensure that the PQ phase transition occurs before inflation. In this case -which seems however to be strongly disfavored by the recent discovery of the tensor modes in the B-mode power spectrum by BICEP2, as has been discussed in Sec. I -topological defects like domain walls created at the phase transition would be wiped out by inflation.
Before assessing the domain wall problem in the models in Sec. IV we highlight the vacuum degeneracy which occurs generically in axion models. We take into account at first just one U (1) PQ symmetry and assume its only explicit break is due to the QCD instantons. Thus, there is a potential V (θ) as a result of the effective interaction L eff = α s 8π C a g θ G·G , where θ(x) = a (x)/f a , and C a g an integer number obtained from Eq. (A11). The potential is not invariant under arbitrary axion field shift transformation a → a + α f a realizing the U (1) PQ symmetry. For the purposes of exposing the domain wall problem we use the parametrization of the singlet field as (omitting the heavy fields) with the PQ charge X σ of σ normalized so that all fields have integer PQ charges. We assume the periodicity θ ≡ θ + 2π, or, equivalently, the periodicity of a field is 2πf a . It means that changing θ → θ + 2π all fields that depend on θ return to the same value. Although the shift symmetry θ → θ + α is broken for arbitrary values of α, a shift with any discrete value in the set α k = 2πk/|C a g |, with k = 0 , ..., |C a g | − 1, remains as a symmetry. It means that if C a g = 1 there are different vacuum expectation values configurations which are all degenerated with respect to the minimum of the potential, i. e., Therefore, the vacuum is invariant under a discrete symmetry subgroup of U (1) PQ where Q is the U (1) P Q charge operator. The order of the discrete group, N , is in fact here the domain wall number N DW and is given by N DW = |C a g |, if there is no subgroup of Z N which acts trivially on the vacuum. When there is a subgroup Z M acting trivially on the vacuum its symmetry group is Z N/M ≡ Z N /Z M , and the domain wall number equal to N DW /M . This is the essence of the domain wall problem in axion models with a single U (1) PQ symmetry. The analysis of this problem for the models in Sec. IV requires additional observations due the fact they have two global chiral symmetries. For these cases the domain wall number is better defined as the disconnected degenerated vacuum configurations. A treatment for the domain wall problem in theories with more than one global chiral symmetry is developed in [159], [160], [161]. In what follows we disregard the mixing and mass corrections due nonrenormalizable effective operators generated by gravitational interactions.
We work on the same basis of the global chiral symmetries for the models in Sec. IV, with K σ1 being the U (1) PQ1 charge of σ 1 , and X σ2 being the U (1) PQ2 charge of σ 2 . But differently of what we have taken before we choose here K σ1 and X σ2 not equal to one.
For the model in Sec. IV A we take K σ1 = 2 and all U (1) PQ1 charges are integers as shown Table XII. This is the same as multiplying by two all U (1) PQ1 charges in Table II. Concerning the U (1) PQ2 charges, starting with the ones defined in Table II we can use the hypercharge gauge transformation to redefine the PQ charges of all fields according to X ψ → X ψ + Y ψ X u , with Y ψ the hypercharge of the field ψ. Using the baryon number symmetry we can still redefine the quark PQ charges by a shift X q → X q − 1 3 X u . Finally, using the normalization X u + X d = X σ2 = 3, all the U (1) PQ2 are then integers as shown in Table XII.  The sets of charges in Table XII are equivalent to the previous ones in Table II in the sense that although the previous anomaly coefficients C ig obtained in Sec. IV A change now to C 1g = K σ1 C 1g and C 2g = X σ2 C 2g the scales also change according to f a 1 = f a 1 K σ1 an f a 2 = f a 2 X σ2 , so that the ratios C ig /f a i do not change. The axion field A is defined as, according Eq. (24), with C 1g = K σ1 = 2, C 2g = N g X σ2 = 9, and θ i (x) = a i (x)/f a i . It can be seen from the parametrization of the scalar fields that both θ 1 and θ 2 have period 2π, i. e., all scalar fields return to the same value when shifting θ 1,2 → θ 1,2 + 2π. The U (1) PQ1 × U (1) PQ2 symmetry transformations leads to θ 1,2 → θ 1,2 + α 1,2 , with α 1,2 being the group parameters. These general transformations are anomalous but, due the fact that C 2g = 9, there is the following subset of discrete transformations leaving de vacuum invariant with k = 0, ..., 8, and Q 2 the U (1) PQ2 PQ charge operator. Although C 1g = 2 the group Z 2 ≡ e iπk Q1 , with k = 0, 1, acts trivially on the scalar and quark condensates so that it is not relevant here. Therefore, there are nine vacuum expectation values configurations differing from one to another by action of Eq. (C6), being all degenerated with respect to the vacuum state. Denoting collectively the vacuum expectation values of scalar fields and quark condensates as the elements { ϕ k , q L q R k } in the set o degenerated vacua, the elements are such that It happens that all elements in the set of vacua are connected. To see this we follow with the arguments in Ref. [161] observing first that since θ i ≡ θ i + 2π the space of possible vacua is a torus. To reach a degenerated point one has to go to a distance α 1 = 2π in the θ 1 direction and α 2 = 2π/9 in the θ 2 direction. Also, there is a group U (1) β whose generator is defined in terms of the generators Q 1,2 of U (1) PQ1,2 according to where the integers are (β 1 , β 2 ) = (−4, 1). These numbers are chosen in order that the U (1) β anomaly is equal to one, To perform a 2π rotation in U (1) β one needs a −2π × 4 in the group U (1) PQ1 and a 2π in the group U (1) PQ2 , but these two rotation brings back to the same point once θ i ≡ θ i + 2πβ i . Thus, all elements in the set of vacuum states are connected, and the degeneracy of the vacuum is just one by the reason that N β = 1, resulting that the domain wall number is N DW = 1. The condition N β = 1 can always be reached when C 1g and C 2g are relative prime [161], as is the case in the present model. Another way to see that the elements { ϕ k , q L q R k } are connected is to note that the symmetry is free from QCD anomaly, and related to the Goldstone boson a. This symmetry can be used to connect all elements transforming the scalar fields and ordinary quarks, with γ = π n 9 , with n = 0, ... 8, so that there is just one vacuum (e iπnQ1 acts trivially on the vacuum).
The analysis above requires that the vacuum expectation values σ 1 and σ 2 are the same. If T 1 is the temperature below which σ 1 get non-vanishing, and T 2 is the temperature below which σ 2 get non-vanishing, it is shown in Ref. [161] that axion strings will form for T 1 > T 2 , (or T 1 < T 2 ). These "σ 1 axion strings" and "σ 2 axion strings" have n 1 = C 1g /K σ1 = 1 and n 2 = C 2g /X σ2 = 3 axion domain walls, respectively. Being n 1 and n 2 not equal to ±1 is still a domain wall problem even if N β = 1 [161]. But a solution may still be found with the inclusion of extra quarks with specific PQ charges in order to furnish n 2 = 1 [158].
For the model in Sec. IV B 1 only U (1) PQ1 has the QCD anomaly. Therefore, the domain wall problem analysis goes as exposed in the beginning of this section (the same reasoning applies also to the model in Sec. IV C). The normalization leaving all the U (1) PQ1 charges integers is such that for the model in Sec. IV B 1 K σ1 = 2, resulting in the anomaly number C 1g = 2. But in this model the discrete group with k = 0, 1, acts trivially on the vacuum defined by the condensates σ 1 and q L q R . Therefore, there is no domain wall problem in the model. For the model in Sec. IV B 2 we also have that only U (1) PQ1 has the QCD anomaly. But in this case the normalization K σ1 = 2 shows that C 1g = 4. The vacuum is invariant under Z 2 ≡ Z 4 /Z 2 . Therefore, the domain wall number is N DW = 2. In order to avoid a domain wall in this model we could still add a second quark singlet Q 2 L,R with the interaction σ * 1 Q 2L Q 2R . This additional quark would lead to C 1g = 2, so that N DW = 1, eliminating the domain wall formation.
In our analysis here it was considered that the only explicit breaking of U (1) PQ1 and U (1) PQ2 is through QCD anomalies. Taking into account the gravitational interactions through Planck scale suppressed nonrenormalizable operators violating global symmetries, degenerated vacuum configurations for values within the periods of a 1,2 may not arise due the corrections for the axion potential like the second term in Eq. (B2). In fact, a solution for the domain wall problem in axion models is to allow for some explicitly breakdown of the PQ symmetry [153]. In our case the large discrete symmetries make a small explicit breaking of U (1) PQ1 and U (1) PQ2 . On the other hand, it is argued in [162] that with only highly suppressed Planck scale suppressed operators long lived domain walls may still be present, posing problems for the standard cosmology in what concerns the Universe evolution.