Axion Predictions in $SO(10)\times U(1)_{\rm PQ}$ Models

Non-supersymmetric Grand Unified $SO(10)\times U(1)_{\rm PQ}$ models have all the ingredients to solve several fundamental problems of particle physics and cosmology -- neutrino masses and mixing, baryogenesis, the non-observation of strong CP violation, dark matter, inflation -- in one stroke. The axion - the pseudo Nambu-Goldstone boson arising from the spontaneous breaking of the $U(1)_{\rm PQ}$ Peccei-Quinn symmetry - is the prime dark matter candidate in this setup. We determine the axion mass and the low energy couplings of the axion to the Standard Model particles, in terms of the relevant gauge symmetry breaking scales. We work out the constraints imposed on the latter by gauge coupling unification. We discuss the cosmological and phenomenological implications.

Employing a bottom-up approach, it has been shown recently that a minimal extension model of the SM -dubbed Standard Model -Axion -Seesaw -Higgs portal inflation model or SMASH -may explain these problems in one stroke [1,2]. It consists in extending the SM by three right-handed extra neutrinos, by an extra quark, and a complex scalar field, which are charged under a new global U (1) PQ Peccei-Quinn (PQ) symmetry. The latter is assumed to be broken spontaneously at an intermediate scale v σ ∼ 10 11 GeV given by the vacuum expectation value of the scalar field. The neutrino flavour oscillation puzzle is solved by the well-known type-I seesaw mechanism [3][4][5][6]: the neutrino mass eigenstates split into a heavy set comprising three states with masses proportional to v σ , composed of mixtures of the new right-handed neutrinos, and a light set of three states with masses inversely proportional to v σ , composed of mixtures of the SM left-handed neutrinos. The extra quark and the excitation of the modulus of the new scalar field also get large masses proportional to v σ , while the excitation of the phase of the new scalar field stays very light, its mass being inversely proportional to v σ . Crucially, this phase field acquires a linear coupling to the gluonic topological charge density from a loop correction due to the extra quark. Correspondingly, it replaces the θ angle of QCD by a dynamical field and thus solves the strong CP problem [7]. Its particle excitation can therefore be identified with the axion [8,9]. Loop effects involving gravitons induce non-minimal gravitational couplings of the Higgs boson and of the new scalar field. These couplings make the scalar potential energy in the Einstein frame convex and asymptotically flat at very large field values.
Correspondingly, the modulus of the new scalar field or a mixture of it with the modulus of the Higgs field can play the role of the inflaton. The baryon asymmetry is produced via thermal leptogenesis [10]. Soon after reheating of the universe and breaking of the U (1) PQ symmetry, the decays of the right-handed neutrinos produce a lepton asymmetry, which is partly converted into a baryon asymmetry by sphaleron interactions before the breaking of the electroweak symmetry. Finally, at temperatures around the QCD transition between a quark-gluon and a hadron plasma phase, dark matter is produced in the form of a condensate of extremely nonrelativistic axions [11][12][13]. To account for all of the cold dark matter in the universe, the U (1) PQ symmetry breaking scale is required to be around v σ ∼ 2 × 10 11 GeV, corresponding an axion mass around 30 µeV [1,2,14,15]. Adding a cosmological constant to account for the present acceleration of the universe, SMASH offers a self-contained description of particle physics, from the electroweak scale to the Planck scale, and of cosmology, from inflation until today.
It is an interesting question whether one can get a similar self-contained description by exploiting a grand unified theory (GUT) extension of the SM: a GUT SMASH variant 1 . In fact, it is well established that GUTs based on the gauge group SO(10) [17,18] may solve the fundamental problems i)-iv) discussed above exploiting the same mechanisms as our bottom-up variant of SMASH [19][20][21][22][23][24]. In fact, the right-handed neutrinos and thus the seesaw mechanism and the possibility of baryogenesis via thermal leptogenesis [25] occur automatically in these models. Moreover, an axion suitable to solve the strong CP problem and to account for the observed amount of dark matter can arise from the rich Higgs sector of these models. An intermediate scale M I ∼ 10 11 GeV (as required in order to have axion dark matter) between the GUT scale and the electroweak scale may arise naturally from the necessity of one or more intermediate gauge groups between the SM gauge group and SO (10) to get gauge coupling unification without invoking TeV scale supersymmetry [26][27][28]. Finally, the Higgs sector of these models necessary for the breaking of SO (10) and its intermediate scale subgroups also provides candidates for the inflaton if they are non-minimally coupled to gravity [29].
As a first step towards an SO(10) GUT SMASH model, in this paper we identify the physical axion field and determine its low energy effective Lagrangian in a set of wellmotivated non-supersymmetric SO(10) × U (1) PQ models -a missing piece in the existing literature. In our treatment we bridge the gap between GUT and low scales: we pay particular attention to low energy constraints, ensuring orthogonality of the physical axion with respect to the gauge bosons of all the broken gauge groups, and we are able to identify the global symmetry associated with the physical axion, given by a combination of the original PQ symmetry and transformations in the Cartan subalgebra of SO (10). We provide calculations of the domain-wall number for the physical axion -which match the expectations from the simple UV symmetries-and we also explore how gauge coupling unification, proton decay, B-L, black hole superradiance and stellar cooling constraints affect the allowed window of axion masses.
The paper is structured as follows. In Section 2 we revisit how a PQ symmetry can be motivated in non-supersymmetric SO(10) models independently of the strong CP problem: it forbids some terms in the Yukawa interactions responsible for the fermion masses and their mixing, thereby crucially improving the economy and predictivity of the models. The general construction of the low energy effective Lagrangian of the axion in theories with multiple scalar fields is reviewed in Section 3. This formalism is then used in Section 4 to work out the axion predictions for a number of promising SO(10) × U (1) PQ models. The constraints imposed by gauge coupling unification are explored in Section 5. We summarise and discuss the cosmological and phenomenological implications of our results in Section 6.

The case for a Peccei-Quinn symmetry
The SM matter content nicely fits in three generations of a 16-dimensional spinorial representation 16 F of SO(10), cf. Table 1. On the other hand, there are many possible Higgs representations corresponding to various possible symmetry stages between SO(10) and SU (3) C ⊗ U (1) em , cf. Table 2. Group theory requires at least the following representations  in order to achieve a full breaking of the rank five group SO (10) down to the rank 4 SM group SU (3) × SU (2) × U (1): • 16 H or 126 H : they reduce the rank by at least one unit, either leaving a rank four SU (5) little group unbroken, or else breaking the SM group.
• 45 H or 54 H or 210 H : they admit for rank five little groups, either SU (5) ⊗ U (1) or different ones, like the Pati-Salam (PS) group SU (4) ⊗ SU (2) ⊗ SU (2) [30]. In the latter case, the intersection of the little group with the SU (5) preserved by a 16 H or 126 H can give the SM gauge group.
We will exploit in our explicit models the 126 H and the 210 H representations. Since 16 F × 16 F = 10 H + 120 H + 126 H , the most general Yukawa couplings involve at most three possible Higgs representations, L Y = 16 F Y 10 10 H + Y 120 120 H + Y 126 126 H 16 F + h.c., (2.1)  Table 2. Decomposition of the scalar multiplets according to the various subgroups in our breaking chains. We only display the multiplets which get nonzero vacuum expectation values (VEVs) in the different models considered in the paper.
where Y 10 and Y 126 are complex symmetric matrices, while Y 120 is complex antisymmetric. It is then natural to ask: what is the minimal Higgs sector to reproduce the observed fermion masses and mixings? Clearly, in order to get fermion mixing at all, one needs at least two distinctive Higgs representations 2 . Out of the six remaining combinations, however, only three turn out to give realistic fermion mass and mixing patterns: 10 H +126 H , 120 H +126 H , and 10 H + 120 H (see for example [31,32] and references therein). From these combinations, the first two are phenomenologically preferred since the 126 H is required for neutrino mass generation via the seesaw mechanism. The first one is the most studied, in particular because it is the one occurring in the minimal supersymmetric version of SO (10). We will also exploit it in our PQ extensions of SO(10), as elaborated next. First of all, it is important to note that the components of 10 H can be chosen to be either real or complex. In the non-supersymmetric case it is natural to assume a real 10 H representation. However, as pointed out in [22,33], this is phenomenologically unacceptable, because it predicts m t ∼ m b . In the alternative case in which the complex conjugate fields differ from the original ones by some extra charge, 10 H = 10 * H , both components are allowed in the Yukawa Lagrangian, L Y = 16 F Y 10 10 H +Ỹ 10   (Throughout the paper we will consider decompositions of representations under the PS gauge group by default). From the above it follows that the fields which can develop a VEV in which the SM subgroup SU (3) C ⊗SU (2) L ⊗U (1) Y is only broken by SU (2) L doublets, as in the standard Higgs mechanism, are (1, 2, 2), (10, 3, 1), (10,1,3), and (15,2,2): as seen in with U (1) R being the usual T 3 generator within the Lie algebra of SU (2) R . Given this embedding of the SM fermion families into PS representations, we can express the fermion mass matrices arising from the interactions in (2.2) after electroweak symmetry breaking as Here, M D , M R and M L enter the neutrino mass matrix defined on the symmetric basis (ν, n) 3 , The three different Yukawa coupling matrices in (2.6) weaken the predictive power of the model. This motivated the authors of Ref. [22] to impose a PQ symmetry [7], under which the fields transform as

8)
126 H → 126 H e −2iα , 3 In the notation of table 1, ν denotes the left-handed neutrinos included in the lepton doublets l, and n designates the right-handed neutrinos.
which forbids the couplingỸ 10 in (2.6) (see also Ref. [33]). As mentioned above, the 126 H alone breaks SO(10) to the experimentally disfavoured SU (5) -or else it would also break the SM group-so that we have to introduce a third Higgs representation to achieve a symmetry breaking pattern that arrives at the SM gauge group at a scale above that of electroweak symmetry breaking. We exploit in this paper the 210 H representation, which has the following PS decomposition: The former allows for a VEV that preserves the SM gauge group, v 210 ≡ (1, 1, 1) 210 . (2.10) We will further assume v L = 0 (see equation (2.4)), which implies M L = 0 in the mass matrices in equations (2.6) and (2.7), thus giving a type-I seesaw, and yielding the following two-step breaking chain: The symmetry breaking VEVs are constrained by the requirement of gauge coupling unification and can be calculated from the renormalisation group running of the coupling constants, see Section 5. v R and v 210 are further constrained by proton decay and leptonnumber violation bounds, but the former still allow for excellent fits to the fermion masses and mixings, as was seen in [23,34,35] and references therein. For a recent analysis of unification with intermediate left-right groups, see [36].

Axion generalities
As argued before, within the SO(10) framework one can use predictivity to motivate a global PQ symmetry under which the fermions are charged, see (2.8). The chiral fermionic content of the theory ensures that the symmetry is anomalous under the GUT group, and by extension under the subgroups that survive at low energies, such as SU (3) C . This allows to embed the axion solution of the strong CP problem in the GUT theory, as such solution requires a spontaneously broken global U (1) symmetry with an SU (3) C anomaly. Moreover, the resulting axion excitation can play the role of dark matter. In this section we will review generalities of axion fields in models with multiple scalar fields. We will first introduce the strong CP problem and its axionic solution, followed by a review on how the axion excitation is identified in terms of the VEVs and PQ charges of the fields, and how its effective Lagrangian is determined. Then we will elaborate on the orthogonality conditions of the physical axion -which imply that the global U (1) symmetry of the axion is not simply given by the PQ symmetry in (2.8) -and on the axion domain-wall number.

The guts of the strong CP problem
Gauge theories with field strength where µνρσ is the Levi-Civita antisymmetric tensor with 1234 = 1, andTr denotes a normalised trace over an arbitrary representation ρ of the Lie algebra. Denoting the Dynkin index of ρ as S(ρ) -defined from the identity Tr ρ ρ(T m )ρ(T n ) = S(ρ)δ mn -one has The theta term L θ can be seen to be total derivative, so that its contribution to the action only picks a surface term and becomes topological. In fact this contribution is proportional to the integer topological charge n top of a given gauge field configuration 4 : The coupling θ is known as the "θ angle", because physics is invariant under shifts θ → θ + 2nπ, n ∈ Z. This follows simply from the fact that the partition function of the theory involves the functional integral (after rotation to Euclidean time): (where ϕ is a shorthand for all the fields in the theory), which is invariant under the above shifts of θ. Crucially, L θ is of the same form as the chiral anomaly. Indeed, let's assume that the gauge theory has Weyl fermions ψ a in representations T a of a non-Abelian group. Let's further assume that the fermions have mass-terms (which could be field-dependent) L ⊃ −1/2M ab ψ a ψ b + c.c. . Then under chiral transformations the associated chiral current J µ is anomalous [37][38][39], where S(ρ a ) is the Dynkin index of the representation ρ a associated with the Weyl fermion ψ a .
Under a chiral transformation the quantum effective action Γ transforms as 5  9) where N is the number of Weyl fermions in nontrivial representations, and N is given in (3.8), is a chiral invariant that will show-up in CP-violating observables.
In the SM one could in principle have θ angles for all gauge groups. The U (1) Y angle θ 1 cannot have any effect, as there are no finite-action hypercharge field configurations with nonzero topological charge. Since the Weyl fermions charged under SU (2) L only get masses by coupling to the right-handed SU (2) L singlets, the arguments given before equation (3.9) imply that the corresponding θ 2 angle is unobservable. 6 However, for the strong interactions there is a physical θ angle, which would contribute to flavor conserving CP-violating observables such as the neutron dipole moment. Current experimental bounds, however, imply |θ phys | < 10 −10 [43], from which the strong CP problem follows: why is θ phys so small?
In a GUT completion of the SM, the SU (3) C θ term can only come from a GUT θ term as in equation (3.1), modulo rotations of fermion fields. The SM θ SM term can be related 5 This can be derived using path integral methods when accounting for the non-invariance of the fermionic measure [40], or, within dimensional regularisation, when accounting for the non-anticommuting character of the regularised chirality operator [41]. 6 The fact that θ2 can be driven to zero without affecting fermion masses -but changing the unobservable θ1-can be also understood in terms of B and L symmetries [42], as they are also associated with opposite rephasings of left and right-handed Weyl spinors.
to the GUT θ GUT by matching physical invariants in the GUT and its low-energy effective description (SM), so that refer to the SM Weyl fermions charged under SU (3) C . Since Arg det M is a sum of phases over the different eigenvalues, it is clear that θ SM is equal to θ GUT plus a combination of phases associated with the heavy fermion states. In any case, the strong CP problem has again a reflection in the GUT theory, as again the physical invariant combination to the right of (3.10) appears to be tuned to a small value. Solving the CP problem at the level of the GUT by explaining the smallness of θ phys guarantees a solution of the strong CP problem, given that θ phys has to match in both the UV and IR descriptions.

The axion solution
The axion solution relies on replacing the chiral invariant θ phys with a combination involving a dynamical field. Once this is achieved, the CP-problem is solved if it can be shown that the now dynamical θ phys (x) gets a potential with a minimum at θ phys (x) = 0. Given (3.9), having a dynamical θ phys is suggestive of a field-dependent Arg det M , which can be achieved if the mass-matrices depend on some complex scalar fields with dynamical phases. With part of the mass-matrices becoming field-dependent, one can combine chiral rotations of the fermions with a rephasing of the complex scalars to define a classical U (1) symmetry of the theory, which, being chiral, becomes anomalous at the quantum level. In order to avoid massless fermions and have a well-defined Arg det M , the complex scalars must get VEVs. Thus, as anticipated, the axion solution involves an anomalous chiral U (1) symmetry acting on fermions coupled to scalars -the PQ symmetry-which must be spontaneously broken. Then θ phys will involve a combination of phases of the scalars which contains the pseudo-Goldstone excitation A(x) corresponding to the spontaneously broken U (1) symmetry, i.e. the axion. Indeed, at low energies, when the massive excitations of the phases are decoupled, one can recover -as will be shown next-an effective Lagrangian involving the axion in which θ always enters in a combination for some dimensionful scale f PQ and constantN , so that θ phys = θ + N /(2N )Arg det M + N A/f PQ , where det M designates the determinant restricted to the fermions which are not charged under the global U (1) symmetry. The last ingredient of the solution to the strong CP problem is the fact that nonperturbative QCD effects generate a nonzero mass for θ phys . As said before, at low energies one predicts that the QCD θ term will involve the effective combinationθ = θ +N A/f PQ , which coincides with θ phys if the quark masses are chosen to be real. Since the partition function of a theory is related to the vacuum energy density, one can obtain an effective potential forθ from the Euclidean partition function of QCD with real masses, supplemented with the topological term Lθ (see (3.1)): where V designates the four-dimensional Euclidean volume. The effective potential can be shown to have a minimum atθ = 0, because for vector-like fermions Z can be written as a path integral involving positive functions times aθ-dependent phase factor [44]. Then the axion solves indeed the CP problem, and the axion mass is given by Above, we defined the axion decay constant f A as withN appearing in the combination θ in equation (3.11) that enters the axion effective Lagrangian. In (3.13), χ = d 2 V eff (θ)/dθ 2 |θ =0 is the topological susceptibility -the variance of the topological charge distribution divided by the four-dimensional Euclidean volume, χ = n 2 top /V . It can be calculated from chiral perturbation theory or with lattice techniques, with recent agreement [14,45]. Within the error of the NLO calculation of [45], one has m A = 57.0 (7) 10 11 GeV f A µeV. (3.15)

Constructing the axion and its effective Lagrangian
Next we may explicitly construct the axion field and its interactions in a theory with Weyl fermions ψ a and complex scalars φ i , and N g gauge groups. We assume a PQ symmetry under which the fermions and scalars have charges q a and q i , respectively, and which is broken spontaneously by VEVs In the broken phase, we may parameterise scalar excitations as The spontaneous breaking of the global PQ symmetry implies the existence of a Goldstone state A [46], the axion, which corresponds to the following excitation of the phases: where f PQ is a dimensionful scale. Canonical normalisation of A -whose kinetic term follows from applying (3.17) to the sum of kinetic terms of the complex scalars-implies From (3.17) one may then derive The effective Lagrangian for the axion can be obtained from the anomalous conservation of the PQ current [47]. The latter is given by (3.20) and satisfies the anomaly equation Using (3.17) in (3.20), one has that at low-energies -when the heavier excitations of the A i are decoupled and can be ignored on the r.h.s of (3.17)-the anomaly equation (3.21) becomes The latter is equivalent to the Euler-Lagrange equations of the following effective interaction Lagrangian L[A] eff [47], with f PQ andN k determined from the VEVs and PQ charges as in equations (3.18) and (3.21). As anticipated earlier, the θ k parameter for a given group (see (3.1)) and the axion enter the low-energy Lagrangian in the combination θ k + 1/f A,k A. The above effective interactions can also be recovered by rotating the scalar phases away from the Yukawa couplings [48]. The PQ-invariant Yukawa couplings induce contributions to the Lagrangian of the form where we used (3.17) and the fact that the PQ invariance of the above Yukawa coupling demands q i + q a + q b = 0; for simplicity, we suppressed the internal indices for the different representations of the gauge groups, and assumed the appropriate gauge-invariant contractions. The phase factors in (3.25) can be removed by field-dependent chiral rotations of the fermions, After the previous rotations, the kinetic terms of the fermions pick extra contributions, given -up to a minus sign-by the axion-fermion interactions in (3.23). Accounting for the fact that the effective action picks up an anomalous term after chiral rotations, one gets the axion-gauge boson interactions in (3.23), again up to a minus sign. The difference in sign in the terms linear in the axion can be removed with a physically irrelevant redefinition A → −A, so that the results are equivalent.
It should be noted that the fermion rotations in (3.26) are not the only possibility to eliminate the axion dependence in the Yukawas as in (3.25). One may make different choices of fermionic rephasings that will give rise to different effective actions. However, as these just differ by redefinitions of the phases of the fermion fields, they will be physically equivalent. In this respect, one may wonder whether these alternative Lagrangians give different bosonic interactions than the ones following from equations (3.23), (3.24), (3.18) and (3.21); in particular, if the resulting value of f A,3 C for the SU (3) C gauge group were to be sensitive to the chosen rephasings, one would predict different axion masses (see (3.15)), in contradiction with the expectation that physical quantities should remain invariant under field redefinitions. Fortunately this is not the case, as is explicitly shown in Appendix A. An important consequence from this is that f A,3 C is fixed by the scalar PQ charges, as is clear from the fact that one can always remove the phases in (3.25) by rotating a single fermion per Yukawa interaction, with a phase fixed by the phase of the scalar field entering the Yukawa coupling. An explicit formula for f A,3 C in terms of the scalar PQ charges is given in equation (A.5) of Appendix A.
All the SO(10) models considered here have the same couplings to neutral gauge bosons, aside from variations in the scale f A ≡ f A,3 C . This is not surprising, as the GUT symmetry relates the different SM gauge groups [47]. Starting from the anomalous Ward identity of PQ phys in the full GUT theory, it follows that the effective Lagrangian of the axion must contain axion-gauge boson interactions as in (3.23), but with a single gauge group SO(10): One has S(10) = 1, and using the orthogonality of the GUT generators belonging to different subgroups, we may writē (3.29) Using the decomposition in (3.28), and using Tr T a T b = 1/2δ ab in the (anti)fundamental representations of SU (N ), one gets (3.30) At low energies α GUT renormalises differently in each term, and identifying α GUT → α s in the QCD term, and α GUT sin 2 θ W → α in the electromagnetic term, one recovers the result in (A.6) -derived with the method of fermion rephasings-but without any assumption on the matter representations. The ratio of the couplings of the axion to gauge bosons is thus predicted to be the same for any SO(10) theory regardless of the matter content [47].
Moreover, the f A scale of the QCD interactions corresponds to that in the grand unified theory, f A,GUT .

Axial basis
It is customary to write the axion-SM fermion couplings in terms of chiral currents of the massive SM fermions: where Ψ = {ψ α ,ψ †,α } are Dirac fermions constructed from the Weyl spinors paired by mass terms. 7 The axial basis is particularly useful when accounting for nonperturbative QCD effects in the axion-nucleon interactions, either because one may use current algebra techniques [8,47], or because the matching between the UV and nucleon theory is simplified when using axial currents [45]. Moreover, as will be seen, the coefficients of the fermion-axion interactions in the axial basis depend only on the scalar PQ charges. One has Ψγ µ γ 5 Ψ = −ψ † σ µ ψ −ψ † σ µψ , from which it follows that the axion-fermion interactions in the general formula (3.23) can be recasted in terms of chiral currents as in (3.31) if the Weyl fermions connected by mass terms have equal PQ charges. This won't be the case in the GUT models considered here, for which the global symmetry associated with the physical axion enforces different charges for the fermions interacting through Yukawas. However, one can always redefine the fermion fields with axion-dependent phases in such a way that one recovers interactions of the form in (3.31), without affecting the axion coupling to neutral gauge bosons or the Yukawa couplings. Consider for example two SM Weyl spinors ψ,ψ with PQ charges q 1 , q 2 , and which can be grouped into a massive Dirac fermion after electroweak symmetry breaking -e.g. {u L , u}, where u L is the upper component of a q doublet, in the notation of table 1. One can always redefine Under such redefinition with opposite phases, the axion couplings to neutral gauge bosons remain invariant, as the redefinition is a non-anomalous vector transformation, rather than a chiral one. On the other hand, the axion couplings to ψ,ψ change to The combination of charges q 1 + q 2 above can be related to the PQ charge of the Higgs that gives a mass to the SM fermion in question, because the Yukawas have to be invariant under the PQ symmetry. From (2.6) it is clear that in the presence of a PQ symmetry (enforcingỸ = 0 in (2.6)) the up quarks receive their masses from the scalars H u , Σ u , and the down quarks and charged leptons from H d , Σ d . Then in the axial basis the axion interaction with u, d quarks and the electron can be expressed in terms of axial currents involving the corresponding Dirac fields U, D, E as where q Hu and q H d are the PQ charges of the Higgses. In regards to the neutrinos, in the models considered here the Weyl spinors ν ⊃ l and and the SM singlet n (see table 1) have nontrivial charges under the physical PQ symmetry. For a high seesaw scale v R (see equations (2.6) and (2.7)), the light physical states in the neutrino sector will be mostly aligned with the ν i . One may always do an axion-dependent phase rotation such that the ν i end up carrying no PQ charge, and which again does not affect the axion couplings to neutral bosons because the neutrinos are singlets under the strong and electromagnetic groups. The physical light neutrinos can then be described with Majorana spinors constructed from the ν i and which do not couple to the axion, in contrast to the other fermion fields in (3.34). The previous arguments seem to imply that the PQ charges of the fermions might be unobservable, since one can choose a basis of fields in which the axion coupling to charged fermions only depends on the scalar PQ charges, and in which the light neutrinos do not couple to the axion. However, an important subtlety is that the axion-dependent rephasings needed to go to this special basis of fermion fields, although they do not affect the couplings of the axion to the neutral gauge bosons, remain anomalous under SU (2) L and thus alter the couplings of the axion to the weak gauge bosons. Thus the latter retain the lost information of the fermionic PQ charges, and if the corresponding couplings were measurable one could potentially distinguish theories differing in the fermionic PQ charges, as for example a GUT model in which the fermionic PQ charges are related to the global symmetries of the GUT theory, versus a DFSZ model in which the fermion charges are not constrained by the GUT symmetry.

Mixing with mesons and the axion-photon coupling
Although the effective Lagrangian in (3.23) includes couplings of the axion to the photon, such interaction is further modified by QCD effects. The reason is essentially that QCD induces a mixing between the axion and the neutral mesons, which in turn couple to photons through the chiral anomaly, involving the sameF F interaction that appears in the axionphoton coupling. The QCD-induced shift of the axion to photons can be computed with current algebra techniques [8,47] or in chiral perturbation theory, with next-to-leading order results provided in [45]. At lowest order, the modification of the coupling can be recovered by noting that the mixing between the axion and neutral mesons can be removed with an appropriate axion-dependent rotation of the meson fields, which however induces an anomalous shift of the action which is precisely the QCD-induced axion-to-photon coupling. This shift is universal and does not depend on the PQ charges of the quarks, and is given by (3.35)

Axion-nucleon interactions
In regards to axion-nucleon interactions, they can also be obtained by current algebra methods [47,49], or alternatively using a nonrelativistic effective theory for nucleons, with couplings determined from lattice data [45]. The axion-nucleon interactions are not universal, and are given in [45] in terms of the coefficients of the UV axion-fermion effective Lagrangian in the axial basis, i.e. with fermion interactions as in (3.31), (3.34), with coefficients fixed by the scalar PQ charges. Equation (3.34) shows that the UV coefficients are simply determined by the scalar charges of the Higgses H u , H d . Then the results in [45] imply the following axion-nucleon interactions in the chiral basis: with (3.37)

The physical axion: orthogonality conditions
In the presence of both gauge and global symmetries, identifying the axion becomes a bit subtle, as the PQ symmetry is not uniquely defined. This is due to the fact that the gauge symmetries themselves are associated with global symmetries, so that any combination of the PQ symmetry plus a global U (1) symmetry associated with one of the gauge groups defines a new global U (1) symmetry. This arbitrariness implies that one cannot readily identify the PQ charges q i that define the axion as in equation (3.19), as well as determine the ensuing axion interactions and domain wall number, all of which depend on the q i . Nevertheless, there is an important physical constraint that allows to uniquely single out a global PQ symmetry PQ phys : its associated axion must correspond to a physical, massless excitation, and thus it must remain orthogonal to the Goldstone bosons of the broken gauge symmetries. This allows to identify the combination of phases that defines the axion, from which one can reconstruct the scalar charges of PQ phys . This will be the procedure used in Section 4 when studying the properties of the axion in concrete SO(10) models. First, the combination of phases should be massless. Suppose the Lagrangian generates a quadratic interaction for a combination of the phase fields, for some coefficients d m . Then one can simply use (3.17) and demand that the term becomes zero, which gives Writing the axion as which can be interpreted as an orthogonality condition between the mass eigenstate m d m A m (see (3.38)) and the axion m c m A m . Another physical constraint on the axion is the fact that it should not mix with the massive gauge bosons. When the complex scalars are charged under gauge groups, which are then broken by the VEVs, there are additional Goldstone bosons associated with the broken gauge generators. These Goldstones are related to the longitudinal polarisations of massive gauge bosons, and should be orthogonal to the axion. The interaction of the scalars with a gauge field A µ contains the following terms: Using equation (3.17), the cancellation of the axion-gauge boson interaction requires the following constraint on the PQ charges: mn v m T a mn q n v n = 0. (note thatq m and q m represent the gauge and PQ charges, respectively). Again, the avoidance of axion-gauge boson mixing can be also interpreted as an orthogonality condition between the axion combination A = i c i A i and the GoldstoneG ≡ i d i A i of the U (1) gauge group. Indeed, the same reasoning behind equation (3.19) implies Under some conditions satisfied in the theories studied in this paper, it will suffice to consider orthogonality of the axion with respect to the Goldstones associated with diagonal U (1) generators in the Cartan subalgebra of the gauge group. This algebra, of dimension equal to the rank of the Lie group (5 for SO (10)) is spanned by the mutually commuting generators H i , i = 1, . . . 5 of the Lie algebra. The commutation property allows to choose representations of matter fields with well defined quantum numbers, or weights, under the elements of the Cartan algebra. In other words, the Cartan generators will be diagonal. Their associated U (1) symmetries correspond to rephasings of the fields, and thus the orthogonality requirements for the axion along the Cartan generators are of the form of equation (3.44). In regards to the orthogonality conditions for the non-diagonal generators outside the Cartan subalgebra, they can be simplified in terms of the weights of the fields and the roots of the algebra (that is, the weights in the adjoint representation). As shown in Appendix B, orthogonality conditions for non-diagonal generators are satisfied if, within each SO(10) multiplet, the field components that get nonzero VEVs have SO(10) weights whose difference is not a root of the Lie algebra. This is the case in the models considered in this paper, with the fields getting nonzero VEVs indicated in table 2; more details are given in Appendix B.
An important consequence following from the constraint of equation (3.44) is that if a field φ n is charged under a given diagonal generator, then if the axion decay constant f A is to involve the VEV φ n = v n / √ 2, then there has to be at least another scalar charged under the same generator as φ n [50]. This simply follows from the fact that, for a VEV v n to contribute to f A , the associated PQ charge q n has to be nonzero (see equations (3.18) and (3.14)). Then in order to have a solution to (3.44) with nonzero chargeq n one needs at least another scalar charged under both PQ and the diagonal generator. This is a consequence of the fact that if only one scalar charged under PQ and the diagonal gauge generator develops a VEV, a physical global unbroken symmetry survives the breaking. Other fields are needed in order to break this surviving symmetry and give rise to an axion. Even when there are several scalars with nonzero VEVs and charged under both PQ and a diagonal generator, then if one expectation value is much larger than the rest, it follows that f A is bound to be of the order of the smaller VEVs. The orthogonality condition (3.44) implies that the PQ charge q h of the heavy field goes as where the superscripts h, l denote the heavy field and the light fields, respectively. Plugging this into (3.14) and (3.18), one gets which shows explicitly that f A is determined by the light VEVs v l m . This can be interpreted in an effective theory framework as follows: as discussed before, the single large VEV v h leaves a global symmetry unbroken, so that the theory with the heavy field integrated out has a new PQ symmetry that can only be broken by the VEVs of the light fields, which will determine the scale f A .
Finally, once the axion has been identified by starting from a general linear combination as in (3.40) and imposing the orthogonality and masslessness constraints, the effective Lagrangian can be determined in terms of the coefficients c i , which encode the charges of the physical PQ symmetry. Indeed, comparing (3.40) with (3.19), one has that The charges q i correspond to the physical global symmetry PQ phys connected to the axion. This symmetry must be a combination of the original global symmetries in the Lagrangian, and one may find the corresponding coefficients by solving a system of linear equations.
Since we expect PQ phys to act as a rephasing of fields, and since all fields have well-defined quantum numbers (weights) under the generators of the Cartan subalgebra of SO(10), it is natural to expect PQ phys to be a combination of the original PQ symmetry and the transformations in the Cartan subalgebra. The latter includes in particular the charges R and B − L of tables 1 and 2. Once the relevant combination of global symmetries has been identified, then one can immediately obtain the ratios q a /f PQ for the fermion fields. This provides all the necessary information to construct the interaction Lagrangian (3.23), together with the QCD induced photon corrections (3.35) and the axion to nucleon interactions in (3.36),(3.37), which only depend on the former ratios. This is clear for the axion-fermion interactions, while for the axion-gauge boson interactions it follows from the fact that (3.14), (3.18) and (3.21) imply Note how the q i /f PQ , q a /f PQ , and f A are invariant under rescalings of the PQ charges, because under q i → c q i , q a → cq a , one also has f PQ → cf PQ , as follows from (3.18). Thus, as expected, the axion effective Lagrangian (3.23) does not depend on the overall normalisation of the PQ symmetry. The same applies to the axion mass (3.15).
We note that, since the PQ phys symmetry of which the axion is a pseudo-Goldstone boson is a combination of global symmetries of the GUT theory, the arguments leading to (3.30) are still valid when one considers the axion of PQ phys , rather than the original PQ symmetry.

Remnant symmetry and domain-wall number
Under a PQ symmetry, which we may assume to be orthogonal to gauge transformations as discussed in the previous section, the scalar phases transform as δ α A i = q i v i . Together with (3.18), this implies that the axion (3.19) transforms as (3.50) The effective Lagrangian accounting for the PQ anomaly, given in (3.23), breaks the continuous PQ symmetry to a discrete subset Like the periodicity of θ discussed around (3.4), this follows from the invariance of the partition function, once the contribution d 4 xL eff in (3.23) is included.
Within the previous transformations, not all of them are necessarily nontrivial, as some may correspond to rephasings of the original scalar phases A i by 2πv i , which leaves all complex scalar fields unchanged. According to (3.17), these trivial symmetries generate the following group of transformations for the axion: Thus the physical symmetry group left after the anomaly is the quotient S phys = S/P . If S phys is a finite group, then any potential generated for the axion will have a finite number of minima, and there will be domain walls. This is because the potential has to be invariant under S phys , so that its transformations relate minima with other degenerate minima. The number of vacua must be then an integer times the dimension of the finite group. The only truly protected degeneracy is that induced by the finite group, and so we expect as many minima as the dimension of the finite group (as happens for the potential of the axion generated by QCD effects). The domain wall number N DW corresponds then to the dimension of the finite group, dim(S/P ): Next we elaborate on a procedure to determine N DW in terms of f PQ , the PQ charges q i and the VEVs v i . Suppose that n min is the minimum number n for which one transformation in S (eq. (3.51)) can be undone with a transformation in P (eq. (3.52)). This implies for some values n i . Then any transformation S(kn min ), k ∈ Z, can also be undone with an element of P , as is clear by doing n min → kn min , n i → kn i in (3.54). This means that any element in S(n) with kn min ≤ n ≤ (k + 1)n min is equivalent, up to a P transformation, to an element in {S(n), 0 ≤ n ≤ n min }. For the extrema of the interval, this follows from our previous arguments showing that all the S(kn min ), k ∈ Z are equivalent to the trivial transformation S(0). For the transformations inside the interval (kn min , (k + 1)n min ) we have The part involving 2πkn min f PQ /N is by hypothesis equivalent to a transformation in P , and the part involving 2πδf PQ /N is a transformation in {S(n), 0 < n < n min }. This proves that all S(n) are equivalent under P to S(n), n ≤ n min . Thus If there is a finite solution for N DW , since S(N DW ) ∼ S(0) (equivalence up to a P transformation), then one has in fact which is the usual finite symmetry associated with domain walls. We may write N DW in terms of the coefficients c i of the axion combination (3.40). Using (3.48) and (3.18) it follows that Again, N DW is invariant under common rescalings of the PQ charges, as these leave c i and f A invariant. A simple case is that in whichN is an integer and the scalar q i charges have at least one common divisor, which could be unity. Let k denote the maximal common divisor. In this case the domain-wall number is simplyN /k. Indeed, writing then the term in brackets in equation (3.56) reaches a minimum integer value when taking N /k is an integer becauseN is a sum of terms proportional to the charges (see (3.21)). The latter have k as their maximal common divisor, and so k is a maximal divisor ofN . It should be stressed that the domain-wall number for the axion corresponding to PQ phys , computed after imposing orthogonality with respect to the gauge bosons, is not the same as the domain-wall number calculated using the above formulae but with the charges of the original PQ symmetry. The reason is as follows. Starting from the original PQ symmetry, without imposing orthogonality conditions, one has a group of discrete transformations S as in (3.51), but defined in terms of the original PQ charges. Similarly, one can define P transformations as in (3.52). When identifying the physically relevant transformations within S, then one has to remove not only the trivial rephasings in P , but also the discrete transformations in the center Z of the gauge group. Thus we may rewrite equation (3.53) more precisely, emphasising the fact that it has assumed orthogonality with respect to gauge transformations, as follows: For SO(10) the center of the group is Z = Z 2 , so that the naive domain wall number computed from the original PQ symmetry (e.g. (2.8)) using equations (3.56) or (3.58) will be two times larger than the actual physical domain-wall number.
The domain-wall number N DW is relevant because the existence of N DW inequivalent, degenerate vacua implies that, once the PQ symmetry is broken and QCD effects generate a nonzero axion mass, the universe becomes populated with patches in which the axion falls into one of the N DW vacua. These patches are separated by domain walls that meet at axion strings, with each string attached to N DW domain walls. Within a domain wall, the axion field has nonzero gradients, so that the walls store a large amount of energy which may in fact overclose the universe, unless the system of domain walls and strings is diluted by inflation -as is the case if the PQ symmetry is broken before the end of inflation and not restored afterwards-or is unstable [51]. The latter can happen if string-wall systems can reconnect in finite size configurations which may shrink to zero size by the emission of relativistic axions or gravitational waves. This is allowed for N DW = 1 [52], when for example a string loop becomes the boundary of a single membrane; however, for N DW > 1 the loops become boundaries of multiple membranes in between which the axion field takes different values, and such configurations cannot be shrunk continuously to a point, which prevents their decay.

Axion properties in various SO(10) × U (1) PQ models
After motivating the PQ symmetry in predictive SO(10) constructions and reviewing its connection to the axion solution to the CP problem, next we study the properties of the axion in SO(10) models, using the results of the previous section. First, we will show that the Peccei-Quinn symmetry (2.8) -postulated to get a predictive scenario for fermion masses and mixing-is phenomenologically unacceptable unless other scalar fields with nonzero PQ charges are introduced. This is because the model with the 10 H and 126 H scalars predicts an axion decay constant at the electroweak scale, which has been ruled out experimentally (for a review, see [53]). Then we will move on to consider models in which the axion decay constant lies at either the unification scale or in between the latter and the electroweak scale. Table 3. Decomposition of the scalar multiplets according to the various subgroups in our breaking chains. We only display the multiplets which get nonzero vacuum expectation values (VEVs) in the different models considered in the paper. "Scale" refers to the contribution to gauge boson masses induced by the VEV of a multiplet, rather than to the mass of the multiplet itself. According to the extended survival hypothesis, we only keep the multiplets which acquire a VEV at lower scales, (with the exception of Σ u , Σ d , which decouple at M BL in order to give rise to a low-energy 2HDM limit). All submultiplets not in the list are assumed to be at the unification scale M U . In all cases,

Models with an axion decay constant at the electroweak scale
Here we consider the minimal scalar content motivated in Section 2, i.e. a 210 H , a 10 H and a 126 H , with the latter two charged under the PQ symmetry in accordance to equation (2.8). The scale f A will be a combination of the VEVs of the fields charged under PQ, i.e. 10 H , 126 H . These VEVs determine the fermion masses, which include the SM fermions -whose masses and associated VEVs must lie below the electroweak scale-and the righthanded neutrinos, which are allowed to be heavy. The mass of the latter is set only by the VEV v R = (10, 1, 3) 126 within 126 H , as follows from equations (2.4), (2.6), (2.7). The U (1) B−L ⊃ SU (4) C symmetry is broken only by the VEVs v R = (10, 1, 3) 126 and v L = (10, 3, 1) 126 . The latter breaks the electroweak symmetry and contributes to light neutrino masses and low-energy lepton number violation, so that v R v L . Then we are at the situation commented at the end of the previous section, in which a gauge symmetry is broken by several VEVs, with a single dominant one. It follows that f A is of the order of the light VEVs, i.e. v L , v 10 u,d , v 126 u,d , which are at the electroweak scale or below. For an overview of the various VEVs and mass scales, see table 3. The corresponding mass scales of the fermions are presented in table 4.
Despite the lack of viability of the model, it is instructive to construct the axion explicitly using the techniques outlined in Section 3; this will serve as a simple example that will pave the way to the computations in viable models. Again, the axion involves the fields charged under PQ and getting nonzero VEVs, which are contained in the 126 H To simplify the notation as much as possible, we will denote the VEVs with v and the phases with A, with appropriate subindices, as in equation (3.16). We define For simplicity, and as was anticipated in Section 2, we consider a zero VEV v L for the (10, 3, 1) 126 multiplet, in order to avoid B − L violation at low energies, and to realise the simplest version of the seesaw mechanism (see equations (2.4) through (2.7)). We will also denote v BL ≡ v R as this is now the only B-L breaking VEV. A general parameterisation of the axion, without knowledge of the PQ charges, can now be written as in (3.40). As detailed in Section 3, we may constrain the previous coefficients by imposing orthogonality with respect to the Goldstone bosons of the broken gauge symmetries, as well as perturbative masslessness. For the gauge constraints, the choice of nonzero VEVs is such that, as commented in 3.4 and and shown in Appendix A the only nontrivial orthogonality conditions are those with respect to the Goldstones associated with the generators in the five-dimensional Cartan subalgebra of the gauge group. Since all the VEVs corresponds to colour singlets, they carry no weights under the two generators of the Cartan subalgebra of SU (3) C ⊃ SU (4) C . By assumption, the fields also carry no electric charge, which eliminates another combination of Cartan generators (see equation (B.10) for the relation between the electric charge and the weights corresponding to the Cartan generators of the group SU (4) C × SU (2) L × SU (2) R ⊃ SO (10)). This leaves two independent Cartan generators giving rise to two nontrivial orthogonality constraints. We may use the Moving on to impose perturbative masslessness, we note that in the scalar potential the term is allowed by both the gauge and PQ-symmetries. After symmetry breaking, these terms induce masses for some combinations of phase fields. Denoting gaugeinvariant contractions by "inv", we have: The orthogonality conditions as in equations (3.38), More massive combinations can be found under closer inspection of the scalar potential, but they cannot give additional constraints on the axion as we already identified four constraints which, together with the requirement for a canonical normalisation of the axion, fix the five independent coefficients c i . Proceeding in this way we can finally conclude that the axion is given, up to a minus sign 8 , by We remind the reader that the above parameters v i , A i are defined by equations (3.16) and (4.1). The axion couplings to matter can be calculated using the results of section 3.3. Equations (3.23) and (3.49) imply that the effective Lagrangian can be simply derived from the ratios q a /f PQ corresponding to the fermions. The ones corresponding to the scalars can be obtained from the scalar ratios q i /f PQ , which can be immediately derived from the axion coefficients c i by using (3.48). Applying the latter identity to the axion combination (4.4), it follows that , From these we may derive the PQ charges q a /f PQ of the Weyl fermions by identifying the appropriate combination of global symmetries in the Lagrangian that gives rise to the charges in (4.5). The physical symmetry PQ phys can be expressed as a combination of the global PQ, U (1) R and U (1) B−L -as anticipated in 3.4, the modification of PQ involves the symmetries within the Cartan algebra of the group: From the conventions in (4.1), the PQ charges in (2.8) and the U (1) R , U (1) B−L charges given in table 2 one deduces: Finally, the values of q a /f PQ for the fermions follow from (4.6) and (4.7), and the charge assignments in table 4: , q n f PQ = 0. (4.8) From the above we may obtain the f A,k using (3.49): As explained in 3.3, the value of f A,3 C only depends on the scalar PQ charges, and can be also obtained from equation (A.5). The simple relations above reproduce exactly the result of equation (3.30), which was derived in the grand unified theory. Calling f A ≡ f A,3 we obtain the effective Lagrangian for the axion where the q f /f PQ factors (which are the same across generations) are given in equation (4.8).
As stated at the end of section 3.3, one may obtain a physically equivalent effective Lagrangian L int by starting from the usual fermion kinetic terms and Yukawa interactions and perform different phase rotations (A.3) that remove the scalar phases in the Yukawa terms. This does not affect the coupling of the axion to the photon; see Appendix A for more details. At low energy, incorporating the QCD effects from the axion-meson mixing in equation (3.35) and the nucleon interactions in (3.36), (3.37), and expressing the electron interactions in the axial basis as in eq. (3.34), the Lagrangian involving the axion, the photon, nucleons and electrons is: (4), where we defined The couplings to fermions coincide with those in the usual DFSZ model [45,54,55], although the relation between the parameter β and the scalar VEVs now involves additional fields. As commented in 3.3.1, potential differences with respect to DFSZ models could come from the axion interactions with the weak bosons, which in the axial basis leading to (4.11) will contain the information of the PQ phys charges of the fermions. The domain-wall number of the model can be calculated from (3.58). We may first consider the "naive" domain wall number obtained by using the PQ charges of equation (2.8), without imposing orthogonality conditions. In this caseN = 12 (see (3.21)) is an integer -note that the value ofN is the same for the GUT group and all its non-Abelian subgroups, as follows from the fact thatN in (3.21) can be expressed as a single trace over all fermions, which fall into GUT representations. On the other hand, the scalar charges have k = 2 as a maximum common divisor. In this situation, as discussed in section 3.5, the domain wall number would beN /k = 6, as corresponds to a DFSZ axion model. On the other hand, using the physical PQ charges in (4.5), the calculation is a bit more involved. Starting from equation (3.58), the quantity in brackets is a rational function of the v i . In order to have an integer result, we must demand that the numerator is proportional to the denominator. This gives a system of equations, as many as there are independent monomials in the denominator. Denoting the minimum integer as n min (which will be the domain-wall number) one has: n min + 3(n 1 + n 2 ) = 0, n min + 3(n 2 + n 3 ) = 0, n min + 3(n 1 + n 4 ) = 0, n min + 3(n 3 + n 4 ) = 0. (4.13) Since the n i are integers, clearly one has N DW = n min = 3. That is, the domain wall number is half of the naive estimate with the unphysical PQ symmetry in (2.8). As discussed around equation (3.61), this is due to the fact that the naive estimate is not taking into account the need to quotient the remnant discrete symmetry by the center Z 2 of the gauge group .
As anticipated at the beginning of this section, all VEVs appearing in f A ≡ f A,3 C in equation (4.9) have to be at the electroweak scale, as follows from the conventions in (4.1) and equation (2.6). Hence, the axion described in this model is visible, being just a GUT-embedded variant of the original Peccei-Quinn-Weinberg-Wilczek model which is phenomenologically unacceptable (for a review, see [53]).
There are several ways to lift the axion decay constant to higher values, which, as follows from the discussion in 3.4, must involve additional scalars with PQ charges. The simplest way is to give a PQ charge to the Higgs field responsible for the GUT symmetry breaking [20]. An alternative way is to introduce a new scalar multiplet, e.g. a 45 H , also charged under the PQ symmetry [21,23]. A third way is to introduce an SO(10) singlet complex scalar field responsible for the U (1) PQ symmetry breaking [24]. We will consider in this paper benchmark models from all these three categories.

Models with axion decay constants at the unification scale
As follows from the arguments in Section 3.4, in order to have a heavy axion one needs at least two fields charged under PQ and getting large VEVs. In the model of the previous section, the scalar 210 H , which was needed to ensure the breaking of the GUT group, was not charged under PQ. Thus the most minimal way to decouple the axion decay constant from the electroweak scale is to extend the PQ symmetry (2.8)  (4.14) The PQ charge of 210 H follows from the requirement of allowing gauge invariant cubic interactions between the 210 H multiplet the other scalars. The only possibility is 210 H 126 H 10 H , which fixes the above PQ charge.
The construction of the axion field in this model goes along the same lines as in the previous section, yet with an added extra phase associated with the φ = (1, 1, 1 see table  2). We now define The general parameterisation of the axion is now where "inv" denotes a projection into gauge-invariant contractions. Masslessness of the axion requires then (4.16) In addition, since φ is a singlet under U (1) R and U (1) B−L , we still have the same constraints from orthogonality as before, equations (4.2) and (4.3). Solving the linear system of equations and normalising, we construct the axion for this model: Note that in the limit M Z M U the axion is just A = A 6 . This follows from the field assignments in (4.25) and the scales in table 2. Therefore, the dominant contribution to the axion field comes from the 210 H . 9 The PQ phys charges of the scalars are now: (4.18) Once more, the global symmetry PQ phys can be expressed as a combination of PQ and the Cartan generators U (1) R and U (1) B−L , as in equation (4.7), but with coefficients s i that now take the values (4.19) 9 Note that the 210H is not charged under U (1)R and U (1)B−L -see table 3-so that the axion can be aligned with the phase of a field getting a large VEV, like φ ⊃ 210H , without violating any orthogonality condition. This is in contrast to the model 4.1.
The PQ phys charges of the fermions can be obtained from (4.6) and (4.19), and the charges of table 4: , q n f PQ = 0. (4.20) The f A,k , following from (3.49), satisfy again the GUT relations in (3.30), and are given by: , so that the axion decay constant is dominated by the GUT-breaking scale. The effective Lagrangian for the axion is as in equation (4.10), with the PQ phys charges in (4.20). At low energies, incorporating QCD effects and going into the axial basis, one gets the DFSZ-like interactions in (4.11), with the the parameter β in (4.12).
As in the previous model, the original PQ symmetry in (4.14) involves scalar charges with a maximum common divisor k = 2, and one has integerN = 12 (common to the GUT group and its non-Abelian subgroups). Thus the naive domain-wall number -without imposing orthogonality of the axion with respect to the gauge fields-is againN /k = 6. To get the physical domain-wall number we may use (3.58). As was done for the previous model, (3.58) can be converted into a system of equations involving n min and integer n i : n min + 3(n 1 + n 2 ) = 0, n min + 3(n 2 + n 3 ) = 0, n min + 3(n 1 + n 4 ) = 0, n min + (n 3 + n 4 ) = 0, n min − 3n 6 = 0. The only other option would be the (15, 1, 1) of the 210 H , which, as mentioned above does not help in lowering f A , as the multiplet contains a GUT-scale VEV. One could consider an additional 210 H , independent of the GUT breaking, but minimality favours a smaller multiplet like the 45 H . We will adopt this choice, and to comply with existing literature [21,23], we choose to use the field σ = (1, 1, 3  The construction of the axion goes analogous to Section 4.2. The VEV v PQ of the 45 H now plays the role of the VEV of the 210 H , so we define: Once more, the initial PQ symmetry in (4.24) has scalar charges with a maximum common divisor k = 2; on the other hand, for the GUT group and its non-Abelian subgroups one has integerN = 12, giving a naive domain-wall number of 6. The physical domain-wall number follows from (3.58), which is equivalent to the following system of equations: n min + 3(n 1 + n 2 ) = 0, n min + 3(n 2 + n 3 ) = 0, n min + 3(n 1 + n 4 ) = 0, n min + 3(n 3 + n 4 ) = 0, n min − 3n 6 = 0. Once again, one has N DW = n min = 3, half of the naive estimate. The effective Lagrangian for the axion is as in (4.10), with the values of f A and q i /f PQ given in equations (4.18) and (4.21). Accounting for QCD effects in the axial basis, one recovers again the DFSZ-like interactions in (4.11), (4.12).
In [23] v PQ was chosen to lie at the same scale as the VEV of the 126 H . In principle, there is no reason for this equality, so we will not use it in our analysis. Generically, as mentioned before we have two physical scales M BL and M PQ associated with the VEVs of 126 H and 45 H . We will now distinguish between two cases: Case A: M PQ > M BL . If the 45 H acquires its VEV before the 126 H , it takes part in the gauge symmetry breaking. This is because, as said before, v PQ breaks the Pati-Salam group to SU (4) C × SU (2) L × U (1) R (see table 2). We are therefore confronted with the following three-step symmetry breaking chain: Both v PQ , related to M PQ , and the VEV v BL , related to M BL , have to be compatible with gauge coupling unification at M U . Since v PQ ∼ 3f A (see (4.26)), this constrains the axion decay constant. Such constraints will be analysed in Section 5.  The axion is given by the same combination of phases as in section 4.3.1, as the construction only depends on the scalar PQ charges. The axion decay constant can be obtained from equation (3.49) -which, applied to models with N 10 extra fermion multiplets in the 10 F , gives (A.5)-substituting the values of q i /f PQ in (4.18), and using N 10 = 2. This gives Due to the extra fermions, the PQ symmetry in (4.29) has nowN = 4, as opposed to the previous value of 12. Again, the scalar PQ charges have a maximum common divisor of 2, so that the naive domain wall number is 2. When taking the quotient of the discrete symmetry group with respect to the center Z 2 of the gauge group, one expects then a physical domain wall number N DW = 1, which gets rid of the domain-wall problem. This can be explicitly checked using equation (3.58), which now implies the following system of equations for the n i and the minimum integer n min that gives the domain-wall number: n min + (n 1 + n 2 ) = 0, n min + (n 2 + n 3 ) = 0, n min + (n 1 + n 4 ) = 0, n min + (n 3 + n 4 ) = 0, n min − n 6 = 0. As expected, we have N DW = n min = 1. In order to obtain an axion effective Lagrangian we need the PQ phys charges of the extra fermionsD, D,L, L in the 10 F (see table 4). As the scalar content of the theory is as in the previous section, we have that, as before, PQ phys is given by (4.6) with the s i in (4.19). Using the 1 R and 1 B−L assignments in table 4, the charges of the extra fermions are: (4.32) The axion interactions are then given by with the values of f A and q i /f PQ given in equations (4.21), (4.18), and (4.32). Including QCD effects and going to the axial basis, the Lagrangian for axion, photon, nucleon and electrons now has a different relative weight between axion and fermion couplings than in the previous DFSZ-like result of (4.11). In the current domain-wall-free model one has now an extra factor of three in the fermionic couplings:

Models with decay constants independent of the gauge symmetry breaking
A third way to lift the PQ-breaking scale from the electroweak scale, also relying on a new scalar charged under the PQ symmetry and getting a large VEV, relies in the introduction of a complex gauge singlet scalar S and exploits thus an even more minimal scalar sector than the previous intermediate scale axion models exploiting an additional 45 H . However, this choice lacks the predictivity of the previous approaches since the singlet S does not participate in the gauge symmetry breaking. As in the models containing a 45 H , the minimal model has a domain wall problem which can be avoided by introduction of two generations of heavy fermions. Under the Peccei Quinn symmetry, the scalar fields transform as follows for the two models:

Constraints on axion properties from gauge coupling unification
In this section we analyse the constraints put on the axion mass by the requirement of gauge coupling unification in the models introduced in the previous section. In each case, we take into account the running of the gauge couplings at two-loop order. A consistent analysis to this order needs to take into account one-loop threshold corrections. The general and modelspecific β-functions and matching conditions are given in Appendix C. 11 Since the scalar masses depend strongly on the parameters of the scalar potential, which are not known a priori, the scalar threshold corrections due to the Higgs scalars cannot be calculated exactly. Instead, we have assumed that the scalar masses are distributed randomly in the interval [ 1 10 M T , 10M T ], where M T (T ∈ {U, BL, PQ}) is the threshold at which these particles acquire their masses. For the RG running, we have employed a modified version of the extended survival hypothesis [60]. According to the latter, scalars get masses of the order of their VEVs, so that the scalars remaining active in the RG at a given scale are those whose VEVs lie below that scale. We consider however an exception [24,33]: in order to have a 2HDM at low energies, we will assume that Σ u and Σ d , the SM doublets in the (15, 2, 2) P S of the 126 H , have masses of the order of M BL . Such a choice is not arbitrary. First, as commented in Section 2, realistic fermion masses require VEVs v 126 u,d of 11 Our analysis has been performed using Mathematica [58]. In the calculation of the beta functions, we have employed the LieART package [59], as well as our own code. The appearance of our plots was enhanced using Szabolcs Horvat's MaTeX package.  (15,2,2) is the mass of the (15, 2, 2) PS multiplet [24,33]. If the mass is of the order of v BL , the desired electroweak-scale VEVs are achieved. Taking into account the scalar content of each model, the surviving multiplets can be read from table 3. The RG equations also take into account the fermionic representations -three generations of fermions in the 16 F representations, and two additional generations in the 10 F in the models defined in equations (4.29) and (4.36).
To sharpen the predictions of our models, we take into account constraints from the non-observation of proton decay, bounds on the B-L breaking scale obtained from fits to fermion masses, as well as black hole superradiance and stellar cooling constraints.
In regards to proton decay, we use a naive estimate for its lifetime, considering only the decay mediated by superheavy gauge bosons [32]. We approximate the lifetime of the proton by τ ∼ [61] τ (p → π 0 e + ) > 1.6 × 10 34 y. In subsequent plots, constraints imposed by current limits from proton decay will be shown in blue.
The constraints on the B-L scale in SO(10) models can be obtained by fitting the observed values of fermion masses and mixing angles to the relationships implied by the gauge symmetry (eq. (2.6)). Such fits have been performed for example in [35] and [34]. In the former the fit was performed at the weak scale, while in the latter it was done at the GUT scale. As in the models in our analysis, [34]  where we have defined Since the fits only determine the ratios v 126 Finally, black hole superradiance constraints arise from the fact that axion condensates around black holes can affect their rotational dynamics and the emission of gravitational waves [62][63][64]. We will show the associated constraints in black. Bounds from stellar cooling arise from taking into account the loss of energy by axion emission due to photon axion-conversion in helium-burning horizontal branch stars in globular clusters [65]. Such constraints will be shown in gray.

Running with one intermediate scale
Let us first consider Model 1 described by (4.14) with PQ charged scalars in the 210 H , 126 H and 10 H representations. Figure 1 shows the predicted running of the gauge couplings for    The unification scale is well above constraints from proton decay. Exploiting the relation The width of the region in Model 1 is exaggerated to make the bar visible. Note that for the Models 2.1, 2.2 and 3 the exclusion of the higher B-L breaking scales comes from an interplay of the proton stability constraint and the limit on v BL . In all these models, a higher B-L breaking scale corresponds to a lower GUT unification scale, which leads to an instability of the proton and is therefore excluded. (5.4) between the mass of the superheavy gauge bosons and the VEV v U and the relation (4.21) between the axion decay constant and the VEVs, we obtain This result is illustrated in the first three lines of figure 3, which summarises our results for the case of vanishing threshold corrections.
As illustrated in figure 2 and as already pointed out in [66], taking into account the possibility of scalar threshold corrections induces large uncertainties in the prediction of the GUT scale, which result in corresponding large uncertainties in the prediction of the axion mass. Including constraints from proton decay limits and the non-observation of black hole superradiance, the allowed range is 2.6 × 10 15 GeV < f A < 3.0 × 10 17 GeV, 1.9 × 10 −11 eV < m A < 2.2 × 10 −9 eV.

An extra multiplet
In the Model 2.1 described in (4.24), the requirement of gauge coupling unification does not sufficiently constrain the system of differential equations to uniquely fix both intermediate scales -we can only infer a relationship between the three unification scales M PQ , M BL and M U . We have calculated this relationship, and also imposed the aforementioned limits on the unification scale and on the B-L breaking scale. Depending on which VEV is bigger, the RG running is different. In case B, i.e. M PQ <M BL , the Peccei-Quinn breaking VEV does not break any gauge symmetries, and thus it is unconstrained by the evolution. In this case, we have essentially a two-step breaking model, in which the two symmetry breaking scales M U and M BL are fixed. In case A, with M PQ >M BL , v PQ breaks the SU (2) R gauge symmetry and is therefore constrained by the evolution. Both cases are indicated in figure 5.
In The upper limit (5.9) on M PQ -derived by proton decay constraints in the case of minimal threshold corrections -can be turned into an upper limit on the axion decay constant and a corresponding lower limit on the axion mass as follows. Since M PQ is the mass of the gauge bosons that become heavy by the SU (2) R → U (1) R breaking, we get (5.11) The corresponding limit on the axion mass follows straightforwardly from This limit is illustrated in lines 4 to 6 of figure 3. This constraint however is still subject to potentially large corrections from scalar threshold effects. In fact, in figure 7 we display the relation between the different unification scales -as in figure 5-but now for randomised scalar threshold corrections for different ranges of v BL . Obviously, the threshold corrections can increase the bound on M PQ and thus on f A . For 10 9 (10 11 ) GeV < v BL < 10 15 GeV, we 12 See also [23], which considered the same model, however only taking into account only one-loop running and a single Higgs doublet at the weak scale. They find in this case MPQ = 1.3 × 10 11 GeV and MU = 1.9 × 10 16 GeV. We have considered three different ranges of allowed B-L breaking scales. The threshold corrections are randomised in the following way: All scalars take masses among the values { 1 10 , 1, 10} times the corresponding threshold scale, where we have taken care not to make proton decay mediating scalars contained in the (6, 2, 2) (Pati-Salam) multiplet light. We have chosen this discrete set of masses in order to focus on the largest possible corrections coming from the mass degeneracies. The scan was performed using 240 different sets of threshold corrections. Allowing the scalars to take masses in the whole interval [ 1 10 , 10] times the threshold scale, one could "fill the gaps" and find even more compatible solutions. These would however not significantly increase the allowed region of M PQ , whose upper and lower limits we are interested in.
get f A < 6.7 × 10 12 GeV and m A > 8.5 × 10 −7 eV, while for 10 13 GeV < v BL < 10 15 GeV no allowed range of f A remains -in this case, the model is excluded. These findings are summarised in lines 4 to 6 of figure 4.  plots are shown in figure 8. After constraining the B-L breaking scale we obtain an upper limit on the axion mass and decay constant in this model. The minimal threshold case is only allowed if M BL can be as low as 10 9 GeV, in this case the maximal allowed f A is 4.2 × 10 11 GeV -this is also shown in lines 7 to 9 of figure 3.

An extra multiplet, and additional fermions
Including all threshold corrections in our random sample, for 10 9 (10 11 ) GeV < v BL < 10 15 GeV, f A is constrained to be smaller than 8.6 × 10 12 GeV. For 10 13 GeV < v BL < 10 15 GeV, no viable solutions were found in the sample -the model is strongly disfavoured in this case. The results on the axion decay constant are summarised for this model in lines 7 to 9 of figure 4.

Models with a scalar singlet
In the simplest Model 3.1, described in (4.35), the Peccei-Quinn breaking is driven by a scalar singlet and the axion mass is unconstrained, cf. line 12 in figure 4. There is no relation between the PQ breaking scale and the two other scales. The possible ranges for M U and M BL however can be read from figure 2 -the extra scalar singlet in this model does not change the running. Also in this model a lower B-L breaking scale is preferred, and the model is excluded if v BL > 10 13 GeV is imposed.
If, however, we employ the mechanism of reference [57] to reduce the domain wall number and introduce additional heavy fermions (Model 3.2, (4.36)), one has to account for how the latter change the running of the gauge couplings above the scale M PQ at which they acquire their masses, if M PQ < M U . Correspondingly we obtain a relation between the scales M U , M PQ and M BL also in this model. However, this dependence -plotted in figure  9 for different sets of threshold corrections-is very weak, and the additional fermions do not change the beta functions enough to introduce additional constraints. We have verifed that the model is still allowed in the entire parameter space of v PQ .
For M PQ > M U , the extra fermions are integrated out above the GUT scale and do not change the running of the three gauge couplings. This case is always allowed, as long as the model without the additional fermions is not ruled out. The only constraint on both Models 3.1 and 3.2 -which we summarise as Model 3-comes then from the B-L breaking scale. For degenerate scalars at the thresholds, we need to allow for v BL as low as 10 9 GeV, as indicated in lines 10 to 12 of figure 3. If we allow variations in the masses of the heavy scalars, values of v BL of order 10 11 − 10 12 GeV are still allowed. For v BL > 10 13 GeV, Model 3 is excluded. This is illustrated in in lines 10 to 12 of figure 4.

Dependence on the proton lifetime
In our analysis we are dealing with three different scales, none of which have been observed so far. Apart from the axion mass, one could also hope to constrain the unification scale in the future by detecting proton decay. The projected sensitivity of the Hyper-Kamiokande Cherenkov detector to the channel p → e + π 0 after 10 years of measurement is 1.3 × 10 35 yr at 90%CL [67]. Assuming that proton decay is observed during the first decade of Hyper-Kamiokande, we can place further (hypothetical) constraints on the axion decay constant in each of our models. In figure 10 we illustrate how an upper bound on the proton decay scale constrains the scale of Peccei-Quinn breaking and ergo the axion decay constant and mass. We make a naive analysis and assume that proton decay is only mediated by the heavy gauge bosons. As a lower bound for the proton lifetime we use the present limit 1.6 × 10 34 yr [61]. As shown in figure 11, an observation of proton decay is very constraining only for our Model 1 -here we obtain 2.6 × 10 15 GeV < f A < 4.0 × 10 15 GeV-, while in the other models the allowed ranges of f A are still rather large. Figure 11. Viable ranges of f A in the hypothetical case of a known proton lifetime between 1.6 × 10 34 and 1.3 × 10 35 years. Allowed regions are plotted in orange. Regions in black are excluded from black hole superradiance constraints, regions in gray from the non-observation of excessive stellar cooling.

Summary and discussion
We have analysed various non-supersymmetric Grand Unified SO(10) × U (1) PQ models for their predictions on the axion mass, the domain wall number, and the low-energy couplings to SM particles. The basic field content of all the considered models consisted of three spinorial 16 F representations of SO(10) representing the fermionic matter content and three Higgs representations: 210 H , 126 H and 10 H , see table 5. The latter have been assumed to take VEVs in such a way that SO(10) is broken along the symmetry breaking chain SO(10) In some of the models, this basic field content was extended by further scalar and fermion representations. This includes an additional scalar in the 45 H , in which case we have Table 5. Field content, PQ charge assignments, and resulting domain wall number N DW in the various SO(10) × U (1) PQ models considered in this paper.
considered a further symmetry breaking chain, In all models one can choose a basis of fermion fields for which the phenomenologically most important couplings to photons (γ), electrons (f = e), protons (f = p) and neutrons (f = n) read, at energies lower than Λ QCD , where , and N DW is the domain-wall number, which in the models considered is either 3 or 1. For N DW = 3 one recovers the results for the DFSZ axion [45,47,54,55], although the microscopic origin of the parameter β differs (as it is determined by the VEVs of four Higgses, as opposed to two in DFSZ models). The fermion fields for which the above interactions are valid are obtained after special axiondependent rotations of the fermion fields that carry charges under the global symmetry PQ phys compatible with the GUT symmetry. These fermion rotations do not act in the same way over all the components of the 16 F multiplets, and thus will "hide" the GUT symmetry, and moreover modify the axion couplings to the weak gauge bosons. A possible measurement of the latter couplings would open up the possibility of recovering the GUT Figure 12. Possible ranges of the axion mass and decay constant consistent with gauge coupling unification in our four models. Regions in black are excluded by constraints from black hole superradiance, regions in dark blue by proton stability constraints. Regions in gray are excluded by stellar cooling constraints from horizontal branch stars in globular clusters [65]. For comparison, we show also the mass regions preferred by axion dark matter (DM) (lines 5 to 7), cf. [68]. Here, the dark regions indicate the ranges where the axion can make up the main part of the observed DM, with the possibility of fine tuning the initial misalignment angle in the scenario where the PQ symmetry is broken before the end of inflation and not restored thereafter (pre-inflationary PQ symmetry breaking scenario). In the light regions, axions could still be DM, but not the dominant part. The remaining regions are not allowed -axions in this mass range would be overabundant. Note that the region in the N DW = 3 case has been derived under the assumption that the PQ symmetry is protected by a discrete symmetry, so that Planck scale suppressed PQ violating operators are allowed at dimension 10 or higher [69]. In the last two lines the projected sensitivities of various experiments are indicated [70][71][72][73][74][75][76][77].
compatible charges under PQ phys , and discriminate these models from e.g. simpler DFSZ scenarios.
The overall phenomenologically viable range in the axion decay constant of these models spans a very wide range, 10 7 GeV f A 10 17 GeV, corresponding to an axion mass m A between 10 −10 eV and 10 −1 eV (see the orange regions in figure 12). These predictions survive constraints from gauge coupling unification, from black hole superradiance (black), from proton decay (blue), and stellar cooling 13 (gray). The features of the different models are summarised next.
Model 1 -employing just the basic field content mentioned above and assuming that all these fields are charged under the PQ symmetry, cf. equation (4.14) and table 5appears to be most predictive. In fact, we infer from the first line in figure 12 that the aforementioned phenomenological constraints result in an axion parameter region 2.6 × 10 15 GeV < f A < 3.0 × 10 17 GeV, 1.9 × 10 −11 eV < m A < 2.2 × 10 −9 eV, (6.4) if we allow the seesaw scale to get as low as v BL 10 9 GeV. The allowed axion mass range moves towards the upper end, m A 2.2 × 10 −9 eV, if we restrict the seesaw scale to higher values, v BL 10 13 GeV, cf. figure 4 (first and second line).
The small axion mass predicted in Model 1 implies that the PQ symmetry has to be broken before and during inflation and must not be restored thereafter [14] (pre-inflationary PQ symmetry breaking scenario). In fact, in the opposite case (post-inflationary PQ symmetry breaking scenario), the axion mass is bounded from below by about 23 µeV [14,15], cf. figure 12. In order for axion cold dark matter not to become overabundant, the initial value of the axion field in the causally connected patch which contains the present universe had to be small, 10 −3 |θ i | = |A(t i )/f A | 10 −2 [14]. Furthermore, the Hubble expansion rate during inflation must have been small, H I 10 9 GeV, to avoid constraints from the non-observation of traces of isocurvature fluctuations in the cosmic microwave background radiation [78][79][80][81]. The latter constraint disqualifies Model 1 as a possible GUT SMASH candidate, since in Higgs portal inflation the Hubble expansion rate during inflation is of order 10 13 GeV H I 10 14 GeV [2].
Remarkably, the predicted axion mass range of Model 1 will be probed in the upcoming decade by the CASPEr experiment [70], cf. figure 12, which aims to detect directly axion dark matter by precision nuclear magnetic resonance techniques. If successful and interpreted in terms of Model 1, one may translate, via (5.5), the measurement of the axion mass into an indirect determination of the mass of the superheavy gauge bosons, i.e. the unification scale, where χ is the topological susceptibility in QCD, χ = [75.6(1.8)(0.9)MeV] 4 [14,45]. The unification scale can be probed complementarily by the next generation of experiments searching for signatures of proton decay, such as Hyper-Kamiokande [82] or DUNE [83]. Models featuring an axion with a larger mass, m A 23 µeV, compatible with the postinflationary PQ symmetry breaking scenario, can only be obtained if the 210 H responsible for the first gauge symmetry breaking at M U has no PQ charge. But in addition -as emphasised already in references [21,22] -the scalar sector has to be extended by yet another PQ charged field obtaining a VEV. Otherwise the axion decay constant is predicted to be as low as the electroweak scale, which is firmly excluded by laboratory experiments and stellar astrophysics. Correspondingly, we considered also models which have additional scalar fields such as a PQ charged 45 H (Models 2.1 and 2.2) or a PQ charged SO(10) singlet complex scalar field S (Models 3.1 and 3.2), cf. table 5. Crucially, in these extended models, the PQ symmetry breaking scale v PQ and the seesaw scale v BL are independent parameters -in fact, the axion field has to be orthogonal to the Goldstone field which is eaten by the gauge bosons getting mass by B-L breaking. This is different in the original SMASH model [1,2], where the PQ symmetry is at the same time also a B-L symmetry.
Model 2.1 features PQ charges for the 16 F , 126 H and 10 H multiplets and an additional PQ charged 45 H , cf. (4.24) and table 5. The range in the axion decay-constant/mass is predicted to be quite distinct from the one of Model 1, see figure 12. Just accounting for gauge coupling unification with scalar threshold corrections, we have found an upper bound on the decay-constant f A < 4.0 × 10 14 GeV, and a corresponding lower bound on the axion mass, m A > 1.4 × 10 −8 eV. Imposing in addition constraints from proton decay, the upper limit on f A comes further down, while constraints on the photon coupling from stellar cooling introduce also a lower limit on f A , 1.3 × 10 7 GeV < f A < 6.7 × 10 12 GeV, 8.5 × 10 −7 eV < m A < 0.5 eV. In this model, both the pre-inflationary as well as the post-inflationary PQ symmetry breaking scenarios are possible, see figure 12. As far as the latter case is concerned, it is important to note that the possibly inherent domain wall problem can be circumvented if the PQ symmetry is only a low energy remnant of an exact discrete symmetry, so that Planck scale suppressed PQ violating operators -which have been argued to be induced in any case by quantum gravity effects [84][85][86][87], and which render string-wall systems with N DW > 1 unstable -occur at dimension ten 14 [69]. In this case, the preferred mass range for dark  matter is a bit higher than the one singled out in the post-inflationary N DW = 1 scenario. In particular, for N DW = 3, as in Model 2.1, it is given by 0.18 meV m A 2.0 meV, cf. figure 12. Intriguingly, for m A = O(10) meV and tan β 0.3, the axion in this model may explain at the same time recently observed stellar cooling excesses observed from helium burning stars, red giants and white dwarfs [77]. Fortunately, there are a number of running (ADMX [71], HAYSTAC [88], ORGAN [89]), presently being assembled (CULTASK [72], QUAX [90]), or planned (ABRACADABRA [73], KLASH [91], MADMAX [74], ORPHEUS [92]) axion dark matter experiments, which cover a large portion of the mass range of interest for axion dark matter in the preinflationary PQ symmetry scenario in Model 2.1, see figure 12. Furthermore, in the meV mass range of interest for the post-inflationary PQ symmetry breaking scenario, the model can be probed by the presently being build fifth force experiment ARIADNE [75] and the proposed helioscope IAXO [76], cf. figure 12.
Model 2.2 adds to the field content of Model 2.1 two PQ charged vector representations 10 F of SO(10), getting their mass from the 45 H and cancelling two units of the domain wall number. Correspondingly, Model 2.2 has no domain wall problem whatsoever. Allowing v BL to be as small as 10 9 GeV, the allowed mass range in this model is very similar to the one of Model 2.1. Taking into account additional limits from gauge coupling unification, proton decay and the limit from stellar cooling, the preferred ranges in this model are 1.3 × 10 7 GeV < f A < 1.5 × 10 13 GeV, 3.8 × 10 −7 eV < m A < 0.5 eV. (6.7) This model can be probed by the same experiments as Model 2.1. Similar to Model 2.1, this model allows for axions in the post-inflationary DM window. Remarkably, the model features a second potential DM candidate: the lightest stable combination of the additional fermions [93][94][95][96]. Therefore we do not need to insist on the axion being 100% of the observed dark matter and can allow for bigger axion masses (cf. the region labelled "subdominant" in the N DW = 1 bar of figure 12). Finally, we return back to the question posed in the introduction, concerning whether the field content of the considered models may allow for a self-contained description of particle physics and cosmology on the same level as the minimal SMASH model. For Model 3, this question can almost certainly be answered in the affirmative, since nearly all the considerations and calculations done in minimal SMASH can be translated to Model 3 if one exploits the (in general non-minimally gravitationally coupled) complex singlet S (or a linear combination with the Higgs) as the inflaton. Note that, although f A in Model 3 is poorly constrained for v BL < 10 13 GeV, there can be further bounds coming from demanding stability in the Higgs direction, as analyzed for the SMASH theory in refs. [1,2]. Such study was beyond the scope of this paper. Although less minimal, Model 2 may also be a good candidate for a GUT SMASH model, and features a more constrained axion mass.

A Invariance of axion-neutral gauge boson couplings under fermionic rephasings
As mentioned in Section 3.3, one may obtain the axion effective Lagrangian by performing redefinitions of the fermionic phases which eliminate the dependence of the Yukawa interactions on the axion field. Although rephasings fixed by the fermionic PQ charges suffice, one may choose alternative redefinitions -all cancelling the axion dependence coming from the Yukawas-which will give rise to different effective actions. These are physically equivalent, as they only differ by field redefinitions whose effects vanish on-shell. In this appendix we show explicitly that, in keeping with these expectations, the SU (3) C axion decay constant f A,3 C -and with it the axion mass-as well as the coupling of the axion to photons are not sensitive to rephasings of the fermion fields. Although the interactions of the axion with fermions and massive gauge boson remain sensitive to the choice of fermionic PQ charges, the different effective Lagrangians should yield identical on-shell results.
To prove the assertions about the couplings of the axion to the massless bosons, we will consider the Yukawa interactions for the Weyl fermion fields q, u, d, l, e, n (see table 1) generated by the SO(10) invariant couplings in (2.2), withỸ 10 = 0 due to the assumed PQ symmetry (2.8). For completeness, we will also add the Yukawas of additional N 10 fermion multiplets in the 10 F of SO(10) coupled to a scalar in the 45 H , as needed in some of the models of section 4: Since they couple to the same fermion fields, H u and Σ u must have identical P Q phys charges q 1 /f PQ ; similarly, H d and Σ d must have a common charge q 2 /f PQ . We also allow a charge q 6 /f PQ for the field σ. 15 Then we may remove the axion contributions to the Yukawa couplings by performing any of the following fermion rotations, parameterised by arbitrarŷ q q ,q l ,q D ,q L :  Under such anomalous transformations, after redefining A → −A, the axion-gauge boson interactions become: where α s = g 2 s /(4π), α = e 2 /(4π) = g 2 2 g 2 1 /(g 2 1 + g 2 2 )/(4π) are the strong and electromagnetic fine-structure constants, and G aµν , F µν denote the components of the strong and electromagnetic field strengths, respectively. The previous result shows explicitly that, although the fermionic PQ charges appear in the effective Lagrangian, the interactions between the axion and the massless bosons only depend on the PQ charges of the scalars, and thus are independent of possible rephasings of the fermions in the low-energy theory. From our results we may obtain an expression for f A ≡ f A,3 C in terms of the scalar PQ charges: Including the axion-fermion interactions arising from the fermion rephasings, we may finally write the effective Lagrangian defined for the general fermion rotations in (A.3) and the above f A value: B Roots and weights of SO(10) and 4 C × 2 L × 2 R In this appendix we review some basic properties of the roots and weights of SO(10)for more details, see [97]-and prove that for our choice of VEVs in table 2 the axion is automatically orthogonal to the Goldstone bosons associated with the broken, non-diagonal generators of the gauge group. The Lie algebra of SO(10) has rank five, as it contains a subspace of mutually commuting generators -the Cartan subalgebra-of dimension five. The latter is spanned by generators H i , i = 1 . . . 5. As the H i are hermitian and commute with each other, in any representation of the algebra one can find a basis of orthonormal eigenstates of the H i , with real eigenvalues λ called "weights". This applies in particular to the adjoint representation, whose weights are called the roots of the Lie algebra. The roots are defined in a basis of generators in which the adjoint action of the H i is diagonal. Aside from the H i , the Lie algebra is then spanned by generators E α satisfying The E α are not hermitian -in fact (B.1) implies E † α = E −α -but one can always recover hermitian generators from the combinations 1/2(E α −E −α ), 1/2i(E α −E −α ). In an arbitrary representation, one can label the states with well-defined weights as |λ , which satisfy The commutation relations (B.1) imply for some normalisation constant N α,λ . The roots α = {α i } and weights λ = {λ i } can be seen as vectors in an Euclidean space of dimension equal to the rank of the group. There is a convenient choice of non-orthonormal basis called the Dynkin basis in which the α i and λ i can be represented with integer numbers. One starts by identifying a set of linearly independent roots α p = {α p i } called "simple roots", such that any root can be expressed as a linear combination of the simple roots, with coefficients that are all positive or zero, or alternatively all negative or zero. The scalar products of the simple roots -defined as (α p , α q ) = i α p i α q i -allow to define a Cartan matrix A ij Then the Dynkin basis is spanned by the following basis of rootsα p = {α p i } In the Dynkin basis, any root or weight can be expressed as a linear combination of theα p with integer coefficients: The five simple roots α p of SO (10) are those at the bottom of the previous list.
The maximal subgroup 4 C 2 L 2 R ⊃ SO(10) has rank 5, and thus its Cartan generators contain those of SO (10), and the two bases of Cartan generators for each group are related by a linear transformation. Then one can map weights of representations of 4 C 2 L 2 R into weights of SO (10). The relation is given by an invertible projection matrix P such that (B.8) With the former choice of matrix, the weights under 4 C 2 L 2 R are of the form {w i }, i = 1, . . . , 5, where w 1 is the weight corresponding to the generator T 3 of 2 L , w 2 is the weight of T 3 for the group 2 R (or, as denoted in tables 2, 1, the 1 R charge), and w 3 , w 4 , w 5 are the three weights of the Cartan algebra of SU (4), with w 3 , w 4 the weights of the Cartan generators T 3 , T 8 of SU (3). The former assignments can be checked by starting from the SO(10) weights of the 16 of SO(10), computing the 4 C 2 L 2 R weights with (B.8), and identifying the quantum numbers of the SM fermions and the RH neutrino, as in table 1.
One can also identify the charge B − L as a combination of SU (4) weights The electric charge is With the previous tools we may prove now that with the choice of fields getting VEVs in table 2, the orthogonality conditions of the axion with respect to Goldstone bosons associated with the off-diagonal gauge generators are always satisfied. The non-diagonal generators of the Lie algebra in a given representation are spanned by the E α . Let's assume a representation of scalar fields in which the nonzero VEVs v i correspond to states |λ(i) . Then the orthogonality constraints (3.43) from off-diagonal generators can be satisfied with the following sufficient conditions: The previous conditions have to be verified within each SO(10) irreducible representation, as the generators only link field components within them. One has This means that (E α ) mn will be zero -and the orthogonality condition with all the E α (and with them the non-diagonal generators) automatically satisfied-if the difference of the weights associated with the v m = 0 is not a nonzero root of the Lie algebra, that is λ(m) = λ(n)+α for all roots α = 0. This will always be the case if only one component in a given representation has a nonzero VEV, but has to be checked for more general situations. If the property holds, then the only nontrivial orthogonality conditions are those arising from (3.44) , see table 2) gets a VEV, so that orthogonality with respect to the off-diagonal generators is guaranteed. For the 10 H and 126 H , however, we have in both cases two field components getting a VEV: the neutral components of H u and H d , and those of Σ u and Σ d , respectively. H u has the same quantum numbers as Σ d , meaning identical weights. A similar relation holds for H d and Σ d . Since the orthogonality condition can be checked in terms of weights, it suffices to consider the Σ components. From table 2 one gets their quantum numbers under 1 R , B − L. The fact that the states are colour singlets implies w 3 = w 4 = 0. The table gives B − L = 0, so that equation (B.9) implies then w 5 = 0. Charge neutrality, together with (B.10), fixes T 3 = −1 R . In the Dynkin basis, the SU (2) weights are twice the usual ones, as follows from relation (B.6), and the fact that there is a unique simple SU (2) root given by the number 2, in the conventional normalisation. Then the 4 C 2 L 2 R weights of the neutral Σ u and Σ d states in the Dynkin basis are One has λ Σ 0 u − λ Σ 0 d = (0, 0, −2, 2, 2) SO(10) , which is not a root of the Lie Algebra, as neither it nor its opposite are within the list in (B.7). This means then that the orthogonality condition (3.43) is satisfied for all non-diagonal generators in the 126 H representation. Identical results apply for the 10 H .

C Coupling evolution
As usual, we can write the renormalisation group equations for the gauge couplings as where i, j indices refer to different subgroups of the unified gauge group at the energy scale µ and The β-function of a gauge coupling g i associated with the gauge group G i at two-loop order in the MS scheme is given by [98] In the above equation, the irreducible fermion and scalar representations are labelled by a and m, respectively. An irreducible representation of a product of groups can contain several copies of irreducible representations of the individual groups. For a fermion representation ρ a we denote the multiplicity of representations of the group i as n a,i ; similarly, for a scalar representation ρ m we use the notation and n m,i . S i (ρ a ) is a shorthand for the the Dynkin index of the irreducible representation of the group i contained within a given fermion representation ρ a . S i (ρ m ) is the analogue for a scalar representation ρ m . C 2 (G i ) = S i (ad) designates the quadratic Casimir for the gauge fields in the adjoint representation of the gauge group i, while C 2,i (ρ a ), C 2,i (ρ m ) are the quadratic Casimirs of the representation of the group i contained in ρ a and ρ m , respectively. Finally, in equation (C.3) one has κ = 1, 1 2 for Dirac and Weyl fermions, respectively, and η = 1, 2 for real and complex scalar fields. At each scale, one has to take care as to which multiplets have to be included in the running. As described in section 5, we consider for the scalars an extended survival hypothesis, modified so as to allow for a 2HDM limit at low energies, while still having electroweak VEVs for all doublets in the 10 and 126, as needed to achieve realistic fermion masses. According to the extended survival hypothesis, fields contribute to the running only if they acquire a VEV at lower scales. The exceptions are the doublets Σ u , Σ d in the (15, 2, 2) PS component of the 126, which are assumed to have a mass of the order of M BL . A list of the scalar components that get VEVs is given in table 2. The decomposition of the fermions is given in table 1. With the previous assumptions, between M W and the lowest intermediate scale, the beta functions for all models mentioned in this paper are the beta functions of a two-Higgs doublet model, with gauge groups given in the order where we used the GUT normalisation for the hypercharge gauge coupling, g Y = 5/3g , which ensures that the generator T Y enters the Lagrangian in the combination g y 3/5T Y , with 3/5T Y a generator with the appropriate GUT normalisation. For a consistent analysis at the two-loop order, at each symmetry breaking scale one needs to impose matching conditions for the gauge couplings that account for finite one-loop thresholds. For a symmetry breaking scale in which each ultraviolet group G UV i is broken down to a subgroup G IR i , the matching conditions for the gauge couplings g i are of the form [99,100]: where, assuming diagonal mass matrices compatible with the infrared gauge symmetries -that is, with a common mass for each IR multiplet-one has For each value of i in the above equation, the V j designate the G IR i representations of gauge bosons that receive a mass at the corresponding threshold, leading to the breaking of the UV group G UV i . S k phys designate the G IR i representations of heavy scalars that are integrated out at the threshold, omitting the unphysical Goldstone bosons. Finally, F l are the G IR i representations of heavy Dirac fermions that decouple at the threshold. The notation of η and κ is as in equation (C.3). We will apply the former matching conditions at the threshold scale µ corresponding to the masses of the heavy gauge bosons, so that the contributions λ V i can be ignored (up to subleading effects from possible lack of degeneracy of the massive gauge bosons from different groups, if the UV gauge group is not simple).
Next we consider the case in which a U IR (1) group arises by combining two U (1) subgroups in the UV, denoted as U 1 (1) ∈ G UV 1 and U 2 (1) ∈ G UV 2 . The associated U (1) generators T IR , T UV 1 , T UV 2 are all part of the Lie Algebra of the GUT group, and for GUT multiplets in representations ρ of the GUT group with Dynkin index S GUT (ρ), they satisfy The k i encode the normalisation of the U(1) generators when embedded into the GUT group, such that define GUT generators with the usual normalisation. Assuming that G UV 1 and G UV 2 become broken at the threshold to G IR 1 and G IR 2 , respectively -so that part of the symmetry breaking is given by G UV 1 ⊗ G UV 2 → G IR 1 ⊗ G IR 2 ⊗ U (1) IR -the matching of couplings goes as: 1 k IR α IR (µ) = 1 k 1 α UV 1 (µ) sin 2 θ 12 −λ (µ) 12π = 1 k 1 α UV 1 (µ) (C.9) In the above equation, g 1 and g 2 are the couplings of the groups G UV 1 and G UV 2 , respectively, and Q 2 IR represent the U (1) charges under the generator T IR . The coupling g IR arising from the previous matching is in the GUT-compatible normalisation, as ensured by the factors of k i . Before moving on to the different models, we provide in table 7 a summary of the decompositions of the different scalar multiplets under the gauge groups appearing in our breaking chains. The table also lists the scales at which the different representations decouple; decompositions are only provided for gauge groups which emerge at or above the decoupling scale, with the exception of fields decoupling at M PQ , since the latter may or may not break the gauge group. For fermion fields, the reader is referred to table 4.  1, 1, 1) VEV v U in the 210 representation, which, given its nonzero PQ charge (see (4.14)), is taken as complex, as are the scalar representations 126 (complex to start with) and the 10. There are 24 broken generators, and correspondingly 24 Goldstone bosons inside the 210 representation, with the same quantum numbers as the broken generators. These Goldstones reside in the real part of the (6, 2, 2) ⊂ 210. According to the extended survival hypothesis, the scalar multiplets which don't get VEVs at lower scales should be integrated out. These are the multiplets Table 7. Decomposition of the scalar multiplets according to the various subgroups in our breaking chains. We list the scales at which the different representations decouple, and for a given representation we don't provide the decomposition under gauge groups that emerge below its decoupling scale, except for fields decoupling at M PQ (depending on the model, M PQ can lead to the breaking of the gauge group, or not). (C. 12) In the previous equations, we have defined We have omitted threshold corrections depending on the masses of the heavy gauge bosons, as we assumed a choice of µ = M U for which these contributions cancel; we will proceed similarly in the rest of the section. (C.14) The matching of the couplings of the groups 3 C and 2 L follows equations (C.5) and (C.6), which yield 2R ⊂ SU (2) R , with associated k 1 = 3/2, k 2 = 1. On the other hand, the GUT-normalised hypercharge coupling g Y has an associated k Y = 3/5. Then the matching goes as , Things go as in C.1, but with the following differences: in models such as the presently analyzed -and the ones that will follow-the 210 is not charged under PQ and can be taken as real, which reduces the threshold corrections. Also, there are new fields in the 45 (which is charged under PQ and thus complex) which have to be integrated out, as they don't get VEVs at lower scales. The scalar multiplets to be integrated out at the M U threshold are then (see table 7 (C.20) The matching conditions(C.5) and (C.6) give now (using the same notation as before) This gives threshold corrections (with notation that should be clear from the above cases) (λ P Q 4C , λ P Q 2L , λ P Q 1R ) = (0, 0, 2) + (0, 0, 1) log P Q M 1,1,1 M 1,1,−1 + (6, 0, 0) log P Q M 10,1,0 + (6, 0, 20) log P Q M 10,1,−1 . (C.23) This case is similar to that in section C.1, with the following differences. First, the Goldstones from the breaking of SU (2) R are now shared between the 45 (whose Goldstones were integrated out at the previous thresholds) and the 126. Additionally, now one must exclude from the loops the heavy gauge bosons that were decoupled at the scale M PQ . All the 45 fields were integrated out at the latter scale, so that the fields that acquire a mass at the scale M BL are (see table 7 Note the difference in Goldstone mode counting with respect to (C.14), and the absence of the descendants of the (10, 1, 0) 4 C 2 L 1 R , (10, 1, −1) 4 C 2 L 1 R . The values of λ BL 3C , λ BL 2L are, as follows follows from equations (C.5) and (C.6), At scales below M PQ and above M BL , the beta coefficients are given by (C.10). At the lowest scales above M Z , we have two Higgs doublet running given by (C.4) as usual. The matching will only differ from that of Model 1 due to the effects of the fermions, as summarised next.   As in (C.33).

D Higher dimensional PQ-violating operators
In models where the Peccei-Quinn symmetry is a low-energy remnant of a discrete global symmetry -which can protect the axion sector from gravitational corrections [85,86,103]one can derive constraints on charges of the scalar fields under such discrete symmetries, see e.g. [69]. (For other works using discrete symmetries to protect the axion's interactions in models with extended gauge groups, see for example [104][105][106], and the recent [107,108]). For the derivation we need to know how the higher-order Peccei-Quinn violating operators that are allowed by the discrete symmetry enter in the axion effective potential. In order to keep in line with the observations of the electric dipole moment of the neutron, one has to ensure that the contributions of these higher order operators are small enough. In the models described in [69], the VEV that breaks the accidental Peccei-Quinn symmetry is the VEV of the additional scalar σ whose phase eventually becomes the axion. The dominant contribution to the axion potential then comes from the PQ violating operator σ N M N −4 p . In Models 2.1 and 3.1 of the present paper, it is not quite so obvious which operator