Peccei-Quinn symmetry and nucleon decay in renormalizable SUSY SO(10)

We suggest simple ways of implementing Peccei-Quinn (PQ) symmetry to solve the strong CP problem in renormalizable SUSY SO(10) models with a minimal Yukawa sector. Realistic fermion mass generation requires that a second pair of Higgs doublets survive down to the PQ scale. We show how unification of gauge couplings can be achieved in this context. Higgsino mediated proton decay rate is strongly suppressed by a factor of (MPQ/MGUT)2, which enables all SUSY particles to have masses of order TeV. With TeV scale SUSY spectrum, p→v¯K+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\to \overline{v}{K}^{+} $$\end{document} decay rate is expected to be in the observable range. Lepton flavor violating processes μ → eγ decay and μ − e conversion in nuclei, induced by the Dirac neutrino Yukawa couplings, are found to be within reach of forthcoming experiments.


Introduction
Grand unified theories (GUTs) [1][2][3] are some of the best motivated extensions of the Standard Model, and have been extensively studied in the literature. Not only do they unify the various forces of nature, they also unify quarks with leptons and particles with antiparticles. Unification of gauge couplings has been known to work very well in the context of TeV scale supersymmetry (SUSY), motivated independently from the Higgs mass hierarchy perspective. SUSY GUTs based on SO(10) gauge symmetry [4,5] would unify all members of a family of quarks and leptons, including the right-handed neutrino, into a single 16-dimensional multiplet. Owing to such a grouping, SUSY SO(10) models are capable of explaining various features of the observed fermion mass spectrum. In particular, renormalizable SUSY SO(10) models which utilize a single 10 and a single 126 of Higgs fields to generate fermion masses provide an excellent fit to all of quark and lepton masses and mixings, including neutrino oscillation data, with a relatively small number of parameters [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. The value of the reactor neutrino mixing angle was predicted in these models before it was measured by the DayaBay collaboration [22], which turned out to be consistent with the prediction. These models may be referred to as SUSY SO(10) models with a minimal Yukawa sector.
GUTs do not shed much insight to the strong CP problem -why the QCD parameter θ takes a value less than 10 −10 . Perhaps the most compelling explanation of the strong CP problem is in terms of the Peccei-Quinn (PQ) symmetry, which is a global U(1) symmetry broken spontaneously by a Higgs field, and also explicitly by the QCD anomaly [23]. The -1 -

JHEP06(2019)045
breaking of the U(1) PQ leads to a near massless scalar, the axion [24][25][26][27][28][29], which may constitute a fraction or the entire dark matter in the universe. The spontaneous PQ symmetry breaking scale should be of order 10 11 -10 12 GeV, to be consistent with direct experimental limits as well as indirect limits from astrophysics and cosmology [30]. It would be of great interest to combine PQ symmetry with SUSY SO(10) models with the minimal Yukawa sector, which has not been done to date. 1 We undertake this task in this paper.
SUSY SU(5) GUT has been extended to include the U(1) PQ symmetry [33]. To achieve this goal in SUSY SO(10) with a minimal Yukawa sector, we identify two key ingredients: (i) an additional 10 of Higgs field, and (ii) a singlet sector that breaks the PQ symmetry in the SUSY limit. Without the additional 10 Higgs field the color triplet partners of the Higgs bosons would survive down to the PQ scale, mediating relatively rapid proton decay (assuming TeV scale masses for SUSY particles). While the SO(10) multiplets 126 and 126 can be utilized to break the PQ symmetry, these fields with their PQ charges being opposite (so that they can have a mass term) would leave a linear combination of U(1) X and U(1) PQ unbroken, where U(1) X is part of SO(10) gauge symmetry. This surviving global U(1) will only be broken spontaneously at the electroweak scale, leading to a weak scale axion model [24,25], which is excluded by direct experiments such as K L → πa searches [30].
One feature that results in the PQ extension of SUSY SO (10) is that an extra pair of Higgs doublets (H u , H d ) survives down to the PQ scale. If their masses were at the GUT scale, the masses of either the up-type quarks or the down-type quarks (and charged leptons) would be suppressed by a factor (M PQ /M GUT ). This would lead to unacceptably small values of b-quark and τ -lepton masses. With (H u , H d ) having masses of order the PQ scale, m b and m τ would not be suppressed and can be fit to their observed values. With an intermediate scale mass for (H u , H d ), the unification of gauge couplings that works well in the MSSM would however be spoiled. We show that lowering the masses of certain colored multiplets slightly from their GUT scale values (by a factor of 10 or so) can compensate for the effects of the extra doublet pair at the PQ scale. These correlated features are also present in the SUSY SU(5) × U(1) PQ models [33].
One important consequence of the PQ embedding of SUSY SO(10) models is that Higgsino-mediated d = 5 proton decay operators [34,35] become suppressed compared to the corresponding non-PQ models. This is a great bonus, as it has been shown that these d = 5 proton decay operators lead to rather fast proton decay in SUSY SO (10) with the minimal Yukawa sector, assuming that all the SUSY particles have masses of order TeV [20]. This problem prompted the suggestion of a mini-split SUSY spectrum in ref. [20] with the gauginos having masses of order TeV and scalars having masses of order 100 TeV. The decay rate of the proton would be suppressed by a factor of (M PQ /M GUT ) 2 compared to the results of ref. [20] in the SUSY SO (10) with PQ symmetry that we present here. This suppression would enable all SUSY particles to have masses of order TeV. This statement would be quantified later on in this paper. Within this setup we find that the decay rate for p → νK + is within reach of ongoing and proposed experiments. With TeV scalars, we also find that lepton flavor violating (LFV) decays µ → eγ and µ−e conversion in nuclei lie

JHEP06(2019)045
in the range that may be observed in forthcoming experiments. Such flavor violations have their origin in the neutrino Dirac Yukawa couplings which are active between M GUT and the B − L symmetry breaking scale v R ∼ 10 12 GeV. Renormalization group flow of SUSY parameters in the momentum range v R ≤ µ ≤ M GUT where ν R 's are active would transfer LFV information to the sleptons, which have masses of order TeV. This in turn would lead to LFV processes such as µ → eγ. The Dirac neutrino Yukawa couplings and the B − L breaking scale are fixed in these models owing to the minimality of the Yukawa sector, leading to crisp predictions for LFV, which depend only on the SUSY particle masses. In this framework, the dark matter candidate can be a mixture of a neutralino-like (wino-like or Higgisno-like) WIMP and a SUSY DFSZ axion [36,37]. An alternative possibly is to have axino as the dark matter [38].

SUSY SO(10) with U(1) PQ
In this section we present a viable model that combines a global U(1) PQ symmetry with SO(10) gauge symmetry in the supersymmetric context. The PQ symmetry solves the strong CP problem; it also enables us to realize TeV scale super-particles consistent with proton decay limits. The model we present is an extension of the renormalizable SUSY SO(10) which preserves the minimality of the Yukawa sector. Fermion families, which belong to the 16-dimensional representations, have Yukawa couplings in these models with a single 10 and a single 126 of Higgs superfiels. The Yukawa superpotential of these models is given by: Here Ψ i stand for the three families of fermions in the 16, H is the 10-plet of Higgs and ∆ is the 126-plet of Higgs. Y 10 and Y 126 are complex symmetric Yukawa coupling matrices, of which one can be chosen diagonal without loss of generality. It was shown in ref. [6] that this Yukawa sector can generate fermion masses consistently, since the ∆ field acquires a large vacuum expectation value (VEV) v R along its SM singlet direction, thus breaking B − L gauge symmetry and supplying Majorana masses to the right-handed neutrinos, as well as electroweak scale VEVs along its SU(2) L doublet directions, contributing to the charged fermion and Dirac neutrino masses. Subsequent analyses have shown that this Yukawa sector can explain small quark mixings and large leptonic mixings simultaneously [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. The model has 12 real parameters and 7 phases to fit 18 measured quantities including the neutrino oscillation parameters, leading to certain predictions. In particular, the model prediction for the reactor neutrino angle was borne out by experiments [22]. To complete the symmetry breaking, a 210 and a 126 Higgs fields need to be employed [39][40][41]. Such a model, with {210 + 126 + 126 + 10} of Higgs fields, can separately be consistent with SO(10) symmetry breaking and fermion mass generation. However, when these two requirements are combined, the model does not fare well [13,[42][43][44] for the following reason. A fit to the neutrino oscillation parameters sets the overall right-handed neutrino mass scale, and the (B − L) symmetry breaking scale v R to be (10 12 -10 13 ) GeV. The breaking of SO(10) symmetry requires, on the other hand, v R ∼ (10 15 -10 16

JHEP06(2019)045
is chosen to be in the phenomenologically viable range of (10 12 -10 13 ) GeV, certain colored Higgs multiplets would acquire masses of order v 2 R /M GUT ∼ (10 8 -10 9 ) GeV. This would spoil perturbative unification of gauge couplings, making the model inconsistent. This problem can be resolved, while maintaining the minimality of the Yukawa sector of eq. (2.1), by the introduction of a 54 Higgs field, which cannot couple to the 16 fermions [20]. In this case, SO(10) symmetry can break down to SU(3) c × SU(2) L × U(1) Y × U(1) B−L at the GUT scale once the 54 and 210 Higgs fields acquire VEVs. The (B − L) symmetry is broken subsequently at a lower scale v R 10 13 GeV when the singlet components of 126 + 126 acquire VEVs. No exotica survives below the GUT scale, except for the singlets from 126 + 126 to facilitate (B − L) symmetry breaking. Consequently, gauge coupling unification is maintained as in the MSSM in this model. 2 The model we present is a PQ symmetric extension of the aforementioned SUSY SO(10) model of ref. [20] with the minimal Yukawa sector of eq. (2.1). The model thus contains 210 + 126 + 126 + 10 + 54 of Higgs fields. To implement the PQ symmetry, we also use an additional 10 and a set of SO(10) singlet fields for breaking the PQ symmetry in the SUSY limit. These fields and their PQ charges are listed in table 1. The U(1) PQ charges are chosen such that the Yukawa superpotential of eq. (2.1) is allowed, and SO(10) symmetry breaking can be achieved consistently. Note that the fields that acquire GUT scale VEVs, viz., 210 and 54, carry no PQ charge.
Without the PQ symmetry, SUSY SO(10) models with the minimal Yukawa sector would require a mini-split SUSY spectrum [20] to suppress the decay rate for p → νK + arising via dimension-5 operators mediated by the Higssinos. In the PQ symmetric model, these baryon number violating operators are induced only after the PQ symmetry is spontaneously broken, leading to a suppression factor of (M PQ /M GUT ) 2 in the decay rate. Thus, neither the mini-split SUSY spectrum of ref. [20] with the minimal Yukawa sector, nor the cancellation mechanism adopted in the color-triplet Higgs Yukawa couplings with an extended Yukawa sector that also includes couplings to a 120 of Higgs boson [45,46], is necessary.
To see the viability of the model and to understand the importance of introducing 10 H -multiplet, we first write down the superpotential involving the Higgs fields consistent with the PQ charge assignment of In the next subsection we will discuss the breaking of PQ symmetry in the SUSY limit that involves SO(10) singlet fields S 1,2,3 . There we will see that all these singlet fields acquire VEVs which will be taken to be of order (10 11 -10 12 ) GeV. Among these fields, owing to our charge assignment, only S 3 has coupling with the SO(10) non-singlet fields.  With the superpotantial of eq. (2.2) the mass spectrum of the SM non-singlet fields can be found readily from the results given in ref. [47]. Certain couplings are absent in our case however, due to the PQ symmetry. One must set m 3 = λ 4 = λ 11 = λ 12 = λ 13 = 0 in the results of ref. [47]. Due to the presence of an additional 10 H -multiplet, which contains SU(2) L doublet and SU(3) c triplet fields, the weak-doublet and color-triplet mass matrices get altered compared to the results of ref. [47]. We present these matrices here, which play important roles in the fermion mass fit and dimension-5 baryon number violation. The explicit form of the doublet and triplet mass matrices would also help understand the need for the 10 H -multiplet in the theory. From the superpotential eq. (2.2) we find the SU(2) L doublet mass matrix to be

JHEP06(2019)045
Here the Φ 1 = (1, 1, 1) , Φ 2 = (1, 1, 15) and Φ 3 = (1, 3, 15) are the VEVs of Φ(210 H ) and the 54 H VEV is E = (1, 1, 1) under the Pati-Salam group SU(2) L × SU(2) R × SU(4) C decomposition. And the color-triplet mass matrix is found to be It is evident from these matrices that if the 10 H field is not present, then the last row and the last column in both the mass matrices would be absent, which would result in one massless state in each sector. Electroweak symmetry breaking would generate mass to JHEP06(2019)045 these states, which is however phenomenologically unacceptable since the light color-triplet would mediate rapid proton decay. One could avoid this by assigning S 3 a PQ charge such that the superpotential coupling S 3 10 2 H is allowed. In this case, the (1,1) entry will be nonzero in the doublet and triplet mass matrices. The color triplet would then acquire a mass of order the PQ breaking scale, of order 10 11 GeV, still leading to rapid proton decay. To give large mass to the color triplet fields and avoid rapid proton decay, the simplest choice is to extend the Higgs sector by the addition of a 10 H with a PQ charge shown in table 1. Note that 10 H has no Yukawa coupling with the fermions owing to the PQ charge, so the Yukawa superpotential remains minimal as in eq. (2.1).
Another important distinction between this PQ symmetric version and the conventional SUSY SO(10) GUTs without PQ symmetry is that a second pair of Higgs doublets would have mass at the PQ scale. Note that in the PQ symmetric limit S 3 → 0, the doublet mass matrix given in eq. (2.4) would take a block-diagonal form: Similarly, if the matrix M D II is fine-tuned to generate a weak scale doublet pair, H u of MSSM would be composed of (H u , ∆ u ) fields, with neither component having Yukawa couplings with the up-type quarks. It would be necessary to do two fine-tunings, one in each sector of M D I and M D II , to generate both up-type and down-type quark masses, which would however result in two pairs of (H u , H d ) at the weak scale, ruining gauge coupling unification. To resolve this issue, mini fine-tunings are required to bring down two pairs of Higgs doublet to the PQ-scale -one from M D I and one from M D II -followed by a large fine-tuning to make one combination of these two pairs at the weak scale. The weak scale doublets so obtained will have components from (H u , ∆ u ) as well as (H d , ∆ d ), leading to realistic masses for both up-type quarks, down-type quarks and charged leptons. Note that the mixing of these two pairs of Higgs doublets at the PQ scale is of order one, induced by the λ 5 S 3 term, which is of the same order as the masses of the two pairs. This mixing cannot be much smaller than unity, or else some of the fermion masses will be too small.
To identify the MSSM Higgs doublets, let us denote the light Higgs doublet pair emerging from M D I as H u The remaining Higgs doublets will all have GUT scale masses, which we shall denote asH. Then the states H u and H d , contained in the 10 H , can be written as where (α, β) and α i are elements of two unitary matrices obeying |α| 2 + |β| 2 = |α 1 | 2 + fields, the superpotential relevant for the light MSSM Higgs doublets can be written as: (2.9) Here µ 11 and µ 22 are the effective masses of the doublet pairs from M D I and M D II respectively, taken to be at the PQ scale. Furthermore, the Yukawa superpotential with the doublets at the PQ scale is given by: Here the matrix elements µ ij in eq. (2.9) and the effective Yukawa couplings y i in eq. (2.10) can be computed straightforwardly. From the Yukawa couplings in eq. (2.10) it is clear that H u 1 and H d 2 must be part of light MSSM states to generate realistic fermion mass. This can be achieved by performing a large fine-tuning by setting the determinant of the mass matrix in eq. (2.9) equal to zero, which reuires µ 11 µ 22 = 0. If µ 22 = 0 is chosen, the MSSM Higgs doublet pair will be identified as: The normalization factors are defined as: The expressions for p i and q i can be found in a straightforward way from the left and the right eigenvectors corresponding to the zero eigenvalue of eq. (2.4). We compute these numerically in our estimate of proton decay rate. Two important parameters, denoted as r and s, which appear in the fermion mass fits are determined in terms of these superpotential parameters: (2.14) The best fit values of these parameters from fermion mass spectrum are given in eq. (3.8).

Symmetry breaking chain
In our set-up the symmetry breaking proceeds in three steps (Chain A): The first step of symmetry breaking occurs at M GUT ≈ 2×10 16 GeV, which leaves the rank of SO(10) group intact. At the next step, the VEVs of 126 H + 126 H multiplets break the U(1) B−L symmetry spontaneously. And in the final step, Higgs doublets from the 10, 10 etc break the electroweak symmetry. We assume this symmetry breaking chain in order to achieve a realistic fit for the fermion masses and mixings, while preserving the successful gauge coupling unification realized in the MSSM. Note that the U(1) B−L breaking scale is determined to be v R ∼ 10 12 GeV M GUT from a fit to neutrino oscillation data with the minimal Yukawa couplings of eq. (2.1). Hence, decoupling of the B − L scale from the GUT scale is demanded by phenomenology. We shall discuss below other possible breaking chains and show that these scenarios are either not viable, or would require large threshold corrections to maintain gauge coupling unification. In this discussion we focus only on the gauge symmetry; the breaking of the global U(1) PQ symmetry will be addressed in the next subsection.
• Chain B. Let us consider the case with 54 H = M GUT 210 H , in which case the symmetry breaking scheme would be SO(10) This is the Pati-Salam gauge symmetry along with a discrete Dparity [48,49], which interchanges SU(2) L and SU (2)  Furthermore, the generation of realistic fermion masses and mixings also requires (2, 2, 1) ⊂ 10 H and (2, 2, 15) ⊂ 126 H to reside at the intermediate scale. Note that, due to the D-parity, the parity partner of (1, 3, 15), viz., (3,1,15) multiplet, would also remain light at the scale M I in this scenario. If M I ≈ v R , then due to the presence of too many particles at the scale of ∼ 10 12 GeV, perturbative unification of the gauge couplings would be spoiled. A scenario with of v R M X < M GUT could be viable; however, in this case large threshold corrections would be needed to maintain gauge coupling unification. symmetry breaking chain is certainly not viable, since it would lead to rapid proton decay induced by the SU(5) gauge bosons, resulting in a lifetime that is in conflict with experimental limits.
• Chain D. If only the (1, 1, 15) ⊂ 210 H field acquires a VEV at the GUT scale, then symmetry breaking proceeds via the left-right symmetric vacuum: left unbroken by the 54 H field, but is broken by the 126 H + 126 H fields. This scenario would lead to pseudo-Goldstone Higgs multiplets carrying SM quantum numbers with masses of order ∼ v 2 R /M GUT ∼ 10 10 GeV. Perurbative unification of gauge couplings would be spoiled in this scenario.
From the above analysis it is clear that our choice of the symmetry breaking chain given in eq. (2.15) (Chain A) is the only one that generates a realistic fermion mass spectrum, while being consistent with gauge coupling unification without relying on excessive threshold effects.

PQ symmetry breaking via singlet fields
Now we focus on the PQ symmetry which also needs to be broken around this intermediate scale, M PQ ≈ v R , although these two symmetries are a priori unrelated. Note that the SM singlet fields σ + σ from 126 H + 126 H carrying two units of B − L charge that break the gauged U(1) B−L cannot simultaneously break the global U(1) PQ symmetry. Additional SM singlet fields carrying non-zero PQ charge are necessary for this purpose. The two sets of fields jointly break U(1) B−L × U(1) PQ symmetry completely. If only the SM singlets from 126 H + 126 H were to acquire VEVs, one linear combination of B − L and P Q symmetries would remain unbroken. To see the unbroken global symmetry we first list the SO(10) → SU(5) × U(1) X decomposition of all the fields:  Table 2. U(1) P Q charges of the color-triplets with P Q = (5P Q − X)/4. Note that except ∆ c and ∆ c multiplets, there is an associated isospin-doublet partner originating from the same Higgs field and carrying identical U(1) P Q charge as that of the color-triplet.
From these, we find that the 54 H , 126 H , 126 H and 210 H VEVs leave a global symmetry, U(1) P Q unbroken, given by P Q = (5P Q − X)/4. It is because all the fields that get VEVs, including σ + σ carry zero charge under P Q . The charges of the isospin doublets and the color triplets under this P Q symmetry are shown in table 2. Since the original PQ symmetry and the U(1) X symmetry commute with SU(5), the P Q charges are the same for all members of a given SU (5) multiplet. The existence of the unbroken global P Q symmetry shows that SO(10) singlet fields are necessary for consistent symmetry breaking. Without these singlets, the P Q symmetry would be broken by Higgs doublets, leading to weak scale axion excluded by experiments.
To break the PQ symmetry consistently, we introduce three singlet fields S 1,2,3 that carry non-trivial charges under U(1) PQ . Their PQ charges are listed in table 1. From various experimental bounds, the PQ breaking scale f a is restricted to be within the range 10 10 GeV f a 10 12 GeV. The superpotantial consisting of the SO(10) singlet fields is given by: (2.21) In the SUSY preserving limit all the F -terms must vanish. They are given by This sets S 1 = 0, and one combination of S 2 and S 3 is fixed. The undetermined VEV leads to a flat direction in S 2,3 . This flat direction is lifted once soft SUSY breaking terms are included. The full potential, including soft SUSY breaking terms of the singlet sector is given by: Straightforward calculation shows that including soft SUSY breaking, all VEVs of the singlet fields are fixed. We find (treating SUSY breaking terms perturbatively)

JHEP06(2019)045
Here we have defined This shows that the singlet sector can consistently break U(1) PQ symmetry. All three singlets acquire masses of order the PQ scale. We shall use the VEV of S 3 as an independent parameter for our fermion mass fit and proton decay calculations, assuming that S 3 ∼ 10 11 GeV. The VEVs of 126 + 126 multiplets along with the singlets carrying non-zero PQ charge together break the symmetry in such a way that a discrete Z 2 symmetry remains unbroken. This is the well known Z 2 subgroup that is identified as the R-parity which stabilize the dark matter in this setup. This R-parity can be identified with the well known matter parity P M = (−1) 3(B−L) under which all the fermions are odd and all the scalars are even.

PQ symmetry and the domain wall problem
Here we briefly discuss the cosmological consequences of the PQ symmetry breaking. It is well known that spontaneous breaking of the PQ symmetry leads to N DW distinct degenerate vacua in axion models due to a residual discrete Z N symmetry. This Z N symmetry, which is a subgroup of U(1) PQ , is left unbroken by non-perturbative QCD effects. In QCD, the vacuum has non-trivial structure and for N f number of flavors has a global SU(N f ) × SU(N f ) × U(1) V × U(1) A symmetry. Of this global symmetry, the axial U(1) A is anomalous and is broken by the QCD instanton effects [50] down to a discrete symmetry Z N , with N = 2N g (N g is the number of quark generations). Topological objects called domain walls [51] are formed as a consequence of this non-trivial structure of the QCD vacuum [52]. The associated domain wall number, N DW , can be computed from the [SU(3) C ] 2 × U(1) PQ anomaly coefficient [30].
For clarity of discussion, let us assume that the B − L symmetry breaking scale v R is slightly above the PQ symmetry breaking scale. Once the singlet components of 126 H and 126 H acquire VEVs, the SO(10)×U(1) P Q symmetry breaks down to MSSM ×U(1) P Q . As discussed in section 2.2, the matter fields at this scale consist of the SM fermions, two-pairs of Higgs doublets, and three MSSM singlet superfields. The charges of the fermion fields under U(1) P Q are listed in table 3, and those of the singlet and doublet scalars are listed in table 4. The superpotential involving these fields is given by Note that the Lagrangian, which includes SUSY breaking terms, would contain bilinear and trilinear scalar couplings analogous to the superpotential terms given above. Table 3. U(1) P Q charges of the fermions with P Q = (5P Q − X)/4.
We now proceed to compute the domain wall number following ref. [32]. The symmetry left unbroken by the QCD instanton effect is Z N . However, N is not necessarily the number of domain walls, since the transformation φ i → φ i + 2πv i is trivial when the fields are expressed in the exponential parameterization. To determine the domain wall number, we write the fields appearing in eq. (2.30) as: We first identify the axion field which should be orthogonal to the massive fields, as well as the Goldstone field eaten up by the Z boson. Eq. (2.30) and its soft SUSY breaking counterparts contain the following relevant terms: The Goldstone boson eaten up by the Z boson is the combination The axion field, which should be orthogonal to these fields is then found to be where we have defined x as Note that 0 ≤ x ≤ 1. The axion decay constant is identified as

JHEP06(2019)045
It is convenient to define the axion field as: 37) in which case the charges q i are found to be: x and from eq. (2.10) they obey the relations: This gives for the integer N in the unbroken Z N for our model to be Then the domain wall number is given by [32]: After a straightforward algebra we find: (2.41) To get integer value for N DW , the n i are strongly restricted, as the monomials in the numerator and denominator should be proportional. We find n 1 = −2n 3 , n 2 = 2n 3 , n d 1 = −n u 1 , n d 2 = −n u 2 . (2.42) Putting these relations back in the original expression gives: This restricts n 3 = n u 1 − n u 2 , and substituting this back gives: So the number of domain walls we have in our theory is N g = 3.

JHEP06(2019)045
Note that, as a result of the PQ phase transition around the temperature of order T PQ ≈ f PQ one-dimensional topological defects called axionic strings are formed [51]. These axion strings radiate axions in between temperatures T PQ and T QCD and at T QCD each string becomes the boundry of N DW domain walls and form stable topological defects that are known as the string-domain wall networks [56][57][58]. It should be pointed out that for N DW = 1, the string-wall network is unstable, decaying into axions soon after their fomration [59,60]. In this case there is no cosmological problem. On the other hand, for N DW ≥ 2 these networks are quite stable. Such stable topological objects are disastrous and would create a serious problem because energy density in the axion domain walls would dominate the universe soon after their formation and overclose the universe. We assume that the axionic string-wall network problem is solved by having inflation with reheating temperature less than the PQ phase transition temperature [61]. When the PQ symmetry is broken before or during the inflation, each domain with a specific value of PQ phase θ a expands to a size larger than the size of the present observable universe resulting in no topological defects. The value of the axion field becomes homogeneous throughout the whole observable universe as a result of the exponential cosmic expansion during the inflationary time. However, further complexity arises since the axion field obtains fluctuations during inflation. When the axions acquire mass at the QCD scale these fluctuations lead to isocurvature density perturbations [62][63][64][65][66][67][68]. Observation of the cosmic microwave background (CMB) puts stringent constraints on the amplitude of these isocurvature perturbations. Solutions of these cosmological problems depend on the details of the inflationary model. By performing lattice simulations, it has been shown that the domain wall problem as well as the isocurvature perturbation problem can be solved [69] simultaneously in general inflationary models that allows intermediate PQ symmetry breaking scale. These numerical computations also show that chaotic inflation can solve both these problems together, although this requires the initial misalignment angle θ a to be somewhat small O(10 −2 ) and axion cannot be the main component of the dark matter in the universe. Similar lattice simulations have been carried out recently within the supersymmetric framework assuming chaotic inflationary scenario in [70] where it is shown that an intermediate scale PQ symmetry breaking, such as the one needed in our model, is completely consistent and devoid of all cosmological problems discussed above.

Suppression of proton decay rate
The Peccei-Quinn symmetry present in the model leads to a strong suppression of proton decay rate induced by the color triplet Higgsinos. We now turn to this issue. In SUSY GUTs the d = 5 baryon number violating operators induced by color triplet Higgsinos are more dominant compared to the d = 6 operators mediated by the gauge bosons [34,35]. As a result, the most dominant modes of proton decay are typically p → νK + and p → νπ + . The Feynman diagram responsible for d = 5 proton decay operators mediated by the color-triplets in our PQ symmetric SUSY SO(10) is shown in figure 1. Note that the d = 5 superpotential operator 16 i 16 j 16 k 16 l is not invariant under PQ symmetry. The net charge −4 of this operator is compensated by charge +4 of S 3 field that breaks PQ symmetry. The singlet sector superpotential allows for S 3 ∼ 10 11 GeV, consistent with constraints As a result, compared to the conventional SUSY SO(10)-GUT, proton decay rate in our theory is suppressed by a factor of (M PQ /M GUT ) 2 ∼ 10 11 /10 16 2 = 10 −10 . Although there is no quantitative prediction, this suppression would still allow p → νK + to be within reach of ongoing and proposed experiments, with a lifetime of order 10 34 yrs., assuming that all SUSY particles that enter into the dressing of the d = 5 operator have masses of order TeV. 4 The lack of a definite prediction for the Higgsino mediated proton decay lifetime is understandable, as it depends on the parameters appearing on the superpotential. Furthermore, the lifetime depends strongly on the SUSY particle masses, with the rate more suppressed with increasing sfermion masses. The case of TeV scale SUSY particle masses is of most interest, as it can provide a solution to the gauge hierarchy problem. We have therefore focused here on the case of TeV scale SUSY spectrum.
Contrary to the d = 5 Higgsino mediated proton decay, the d = 6 gauge mediated proton decay rate can be predicted in terms of the GUT scale gauge boson masses. The dominant decay mode resulting from the d = 6 baryon number violating operators is p → e + π 0 . If the threshold corrections arising from the GUT scale are not significant, then the prediction for proton lifetime is similar, but somewhat shorter (due to new contributions from the X and Y gauge bosons), compared to SUSY SU(5) prediction. However, the mass of the (X, Y ) and (X , Y ) gauge bosons depend on the GUT scale threshold effects, which can cause considerable uncertainty in the prediction. The d = 6 proton lifetime can JHEP06(2019)045 be in the range 10 34-36 years, if the GUT scale threshold corrections allow the gauge boson masses to deviate from the GUT scale of M GUT ∼ 2 × 10 16 GeV by a factor of 3.
We can in fact calculate the Higgsino mediated proton decay rate exactly in our model. Dressing of the d = 5 LLLL operators by the Wino gives the most dominant contribution to proton decay, so we shall focus on this operator. To compute the decay rate, we define the relevant amplitude functions A ν ijkρ as [20,[72][73][74][75][76] as: where the amplitude functionsÂ ν 1bcρ and the Yukawa couplings H (i) , F (i) are defined in ref. [20]. The parameters x, y, z are defined as [20,77]: Here, H and F are Yukawa coupling matrices obtained within the model from a fit to fermion masses and mixings. The best fit values of these matrices from fermion data are given in eqs. (3.10)-(3.9) in the next section. With these as input one can compute the proton decay rate following the procedure explained in ref. [20].

Results
In this section we present our results. First we discuss the fit to the fermion masses and mixings, then we present the details of the proton decay rate calculation, then we discuss the gauge coupling unification within our framework, and finally address lepton flavor violation predicted by the model.

Fit to fermion masses and mixings
The Yukawa couplings of the model are given in eq. (2.1). This is identical to the case of SUSY SO(10) without PQ symmetry. So we closely follow the definitions and the parametrizations of ref. [20] in the context of non-PQ SO(10) model. Following the same notation as ref. [20], at the GUT scale, one has the SO(10) relations:

JHEP06(2019)045
Here Y U,D,E are the MSSM Yukawa couplings of up-quarks, down-quarks and charged leptons. We have defined Here v u and v d are the VEVs of the MSSM fields H MSSM u and H MSSM d , and we define as usual tan β = v u /v d . Although, the light neutrino masses receive contributions from type-I seesaw as well as type-II seesaw, the magnitude of the type-II contribution comes out to be small, because the weak triplets in the model have masses of order the GUT scale. Thus we work in the type-I seesaw dominance scenario [78][79][80][81][82].
The minimal Yukawa sector of SUSY SO(10)-GUT has 12 real parameters and 7 phases to fit 18 observables. To fit the fermion masses and mixings we perform a χ 2 -analysis at the GUT scale, which we take to be M GUT = 2 × 10 16 GeV. We minimize the function χ 2 = i P 2 i to achieve the best fit, where the pull is defined as P i = (O i th − E i exp )/σ i . Here σ i represents experimental 1σ uncertainty and O i th and E i exp represent theoretical prediction and experimental central value of observable i. To get the GUT scale values of the observables in the charged fermion sector we take the central values at the M Z scale from table-1 of ref. [83]. For neutrino observables, the low energy values are taken from ref. [84]. With these inputs, we do the RGE running of the Yukawa couplings [85,86], the CKM parameters [87] and the effective couplings of the neutrino d = 5 operator, κ [88][89][90] within the SM up to the SUSY-scale, which we choose to be 1 TeV. Above 1 TeV the full MSSM is restored, so we use the relevant MSSM RGEs [91][92][93] and evolve them up to the GUT scale. For the charged fermion observables, we fit to these evolved values at the GUT scale. Since the requirement from the neutrino data is to have right-handed neutrinos at some lower scale than the GUT scale in the type-I seesaw scenario, in the fit procedure we include the threshold corrections due to the right-handed neutrinos from the v R scale to the GUT scale [20] (see also ref. [94]). So above the SUSY scale, the effective couplings of the neutrino d = 5 operator running is performed up to the intermediate scale instead of the GUT scale.
For our numerical analysis we fix tan β = 10. Our best fit observables are presented in table 5. We see from this table that the model gives an excellent fit to the data, with a total χ 2 = 6.3. Most of the contributions to χ 2 arises from the d-quark mass, which is 2.3 σ below the central value. All other observables are within 1 σ, providing an excellent overall fit.
In table 6 we list the predictions of the model in the fermion sector for quantities that are currently unknown. The central value of the CP violating Dirac phase relevant for neutrino oscillations is found to be δ P M N S = 327 0 . However, it should be noted that large deviations from the best fit are possible, owing to existence of nearby local minima with acceptable χ 2 . The possible variation of this phase as a function of χ 2 was presented in figure 1 of ref. [20]. Similar results should also hold in our current setup, although we have not investigated this in detail.
In the neutrino sector, the ordering of masses is normal, with the effective mass for neutrinoless double beta decay found to be m ββ 5 meV. While this is well below the  Table 5. Best fit values of the observables corresponding to the type-I dominance seesaw scenarios for tan β = 10. This best fit corresponds to total χ 2 = 6.3. For the associated 1 σ uncertainties of the observables at the GUT scale, we keep the same percentage uncertainty with respect to the central value of each quantity as that at the M Z scale. For the charged lepton Yukawa couplings at the GUT scale, a relative uncertainty of 0.1% is assumed in order to take into account the theoretical uncertainties arising for example from threshold effects.  Table 6. Predictions corresponding to the best fit values presented in table 5 for type-I dominance seesaw scenario. m i are the light neutrino masses, M i are the right handed neutrino masses, α 21,31 are the Majorana phases following the PDG parametrization, m cos = i m i , m β = i |U ei | 2 m i is the effective mass parameter for beta-decay and m ββ = | i U 2 ei m i | is the effective mass parameter for neutrinoless double beta decay. current limit from Kamland-Zen experiment of (61-165) meV [95], future improvements can potentially probe this prediction.

Quantity
The Yukawa coupling parameters corresponding to the best fit are given by:

Proton decay calculation
For the proton decay rate calculation, we proceed in the following way: -We first solve the stationary conditions corresponding to the superpotential of eq. (2.2). These conditions, which can be found in refs. [47,96], are solved for the mass parameters E, Φ 1 , m 1 , m 2 and m 5 . Note that m 5 does not play any role for proton decay calculation. We follow the symmetry breaking pattern shown in eq. (2.15).
-We set λ 1 = 1 and Φ 3 = M GUT . We choose M GUT = 2 × 10 16 GeV. Furthermore, we fix v R = v R = S 3 = 10 12 GeV. Then the dimensionless free parameters of the theory are: λ 2 , λ 3 , λ 4 , λ 5 , λ 10 and λ 13 which we take to be real. There is one free parameter with dimension of mass Φ 2 , which is taken to be complex. We scan over this parameter set {λ 2 , λ 3 , λ 4 , λ 5 , λ 10 , λ 13 , Φ 2 } such that the dimensionless couplings are restricted to be |λ -The numerical values of r and s parameters defined in eq. (2.14) are fixed by the best fit of the fermion spectrum given in eq. (3.8). While computing the proton decay rate, we also impose the constraints that the superpotential parameters correctly reproduce r and s values given in eq. (3.8).
-With each parameter set, satisfying all of the above constraints, we compute the proton decay amplitude function defined in eq. (2.47) that contains the Yukawa couplings, for which we use the best fit values given in eqs. (3.9)-(3.10) by going to the physical basis of fermions (for details see ref. [20]).
-By using the amplitude function eq. (2.47) calculated as above, we estimate the minimum value of the sfermion mass, m S that satisfies all the experimental proton decay bounds: τ p (p → K + ν) ≥ 5.9 × 10 33 years and τ p (p → π + ν) ≥ 3.9 × 10 32 years [97]. For this computation we fix the Wino mass to be m W = 1 TeV. We have taken the relevant nuclear matrix elements from ref. [98].
In figure 2, we present our result on the allowed range for SUSY scalar masses arising from proton lifetime limits obtained by following the procedure outlined above. This figure shows the probability distribution as a function of the minimum required sfermion mass, m S , to satisfy the experimental lower bounds on proton lifetime when m S is varied in the range 3 TeV ≤ m S ≤ 100 TeV. We see that, owing to the suppression in the rate from PQ symmetry, SUSY scalar masses of 3 TeV is fully consistent. Even lower value, e.g., 1 TeV is found to be consistent with proton decay. In the analysis performed above, for simplicity,  Note that, for the above sample point, λ 5 S 3 ∼ 10 12 GeV. Below, we present a second consistent sample point for which λ 5 S 3 ∼ 10 11 GeV. This second sample point that is consistent with the proton decay bounds for m S = 3 TeV and reproduces r and s exactly is given by the following parameter set:

JHEP06(2019)045
Multiplet, φ Running coefficient  Table 7. Three example scenarios of gauge coupling unification which invoke small threshold correction from a single color-sextet multiplet and its conjugate and the corresponding d = 6 proton decay lifetime. Concerning the quoted values of the proton decay lifetime, we stress the fact that these are just estimates, since we have not included the threshold corrections from the other fields except the color-sextet. Each of these scenarios can be considered viable once all the threshold corrections are taken into account.
We see from the analysis above that proton lifetime for decay into νK + is close to the current experimental limit. While it is difficult to make this statement more precise owing to uncertainties in the superpotential parameters, we expect proton decay to be within reach of ongoing and forthcoming deep underground experiments.

Gauge coupling unification
We now comment on the gauge coupling unification in our model. Note that the automatic gauge coupling unification of the MSSM is not realized in our scenario. This is because of the appearance of an additional SU(2) L doublet Higgs fields at the PQ breaking scale. However, since the doublet contributions to the running of gauge couplings are not very significant, small threshold corrections near the GUT scale can restore the coupling unification. With the MSSM spectrum all the way to the PQ symmetry scale, and an extra pair of doublets at the PQ scale, unification can be realized by lowering the mass of a single GUT-scale multiplet carrying color. It is sufficient to lower the mass by an order of magnitude below the GUT scale. To demonstrate gauge coupling unification within our set-up, we stick to this simple scenario where only one pair of fields has a mass around 10 15 GeV. We choose this multiplets to be a color-sextet, isospin-singlet field and its conjugate. In our model three such fields are present: (6, 1, 1/3), (6, 1, −2/3) and (6, 1, 4/3), along with their conjugate fields. The possibilities of coupling unification with these color-sextets are summarized in table 7 and the corresponding unification plots are presented in figure 3. To obtain this results, we first run the 1-loop MSSM gauge coupling constants from 1 TeV to the PQ scale that we fix to be 10 12 GeV, where a pair of Higgs doublet lies; then from the PQ scale to the color-sextet mass scale including the RGE effects of the extra doublet pair; and from there upto the unification scale including the effects of the color sextet fields. The color-sextet mass is appropriately chosen to achieve unification. At the TeV scale the gauge couplings are taken to be g 1 = 0.46774, g 2 = 0.63990 and g 3 = 1.0597 [83]. Among these three different options, lowering the mass of the (6, 1, −2/3)-multiplet is the most favorable. The mass of this multiplet is given by (m 5 − 2 3/5λ 8 E), which can be tuned to a value below the GUT scale without making other multiplets light. Also, in this mass formula, the dimensionless parameter λ 8 is free and does not play any role in the proton decay   calculation. (The mass parameter m 5 also does not enter in proton decay calculation and is fixed by one of the stationary conditions). The corresponding unification scale is found to be ∼ 10 16 GeV and the unified gauge coupling is α −1 GUT = 24.46, with these, the expected d = 6 proton decay p → e + π 0 lifetime can be estimated [99] to be τ p ∼ 9 × 10 34 yrs, which is above the current experimental lower limit set by Super-Kamiokande, τ p > 1.6 × 10 34 yrs [100]. In table 7, along with the unification constraints, we also estimate the lifetime for the other two cases. For the case of (6, 1, 4/3), the lifetime for the p → e + π 0 mode came out to be smaller than the present experimental lower limit. However, these are just estimates, since we have not included the threshold corrections from the other fields except the color-sextet. Each of these scenarios can be considered viable once all the threshold corrections are taken into account.

Lepton flavor violation
In this section we discuss lepton flavor violating (LFV) rare decays predicted by our model. Assuming flavor universality of the SUSY breaking mass parameters at the GUT scale, the main source of LFV in SUSY SO (10) is the presence of the Dirac neutrino coupling Y ν ν c H u of eq. (3.4), which is a combination of Yukawa coupling of the 10 H and 126 H multiplets. These ν c fields have masses of order v R ∼ 10 12 GeV, and the renormalization group running of the SUSY mass parameters in the momentum range v R ≤ µ ≤ M GUT would impart the flavor structure of the neutrino Dirac Yukawa coupling to the sfermions. This RGE-induced off-diagonal entries of the left-handed slepton mass matrix in the leading log approximation can be expressed as: where, Y ν is the Dirac neutrino Yukawa coupling in a basis where the charged lepton and right-handed neutrino mass matrices are diagonal. Assuming scalar mass universality and gaugino mass unification as is usually done in constrained MSSSM, the number of parameters in the SUSY breaking sector is reduced to the set: {m 0 , m 1/2 , A 0 , sgn(µ), tan β}, where, m 0 is the common scalar mass, m 1/2 is the common gaugino mass, A 0 the universal tri-linear coupling and µ the Higgs mass term. For LFV calculation, we restrict ourselves to the cMSSM scenario that will also fix the µ parameter from electroweak symmetry breaking condition. We choose sgn(µ) > 0 (although the LFV results are essentially the same for negative µ), and set tan β = 10 corresponding to the fermion fit. The fit to fermion spectrum fully determines the Dirac coupling matrix Y ν and the right-handed neutrino masses M R k in our framework. As a result we can now compute the rates for LFV processes as functions of only the SUSY breaking parameters. With our choice, this set is reduced to {m 0 , m 1/2 , A 0 }. As is well recognized, cMSSM with TeV scale scalar masses cannot reproduce a Higgs mass of 125 GeV, unless A 0 takes rather large values of order 10 TeV. Here we are interested in LFV processes within MSSM that can correctly produce the Higgs mass as well. We shall thus allow for some new contributions, for example from stop quarks which have non--23 -JHEP06(2019)045 universal masses, for the Higgs mass. For our numerical analysis, we will fix few different values of A 0 of order TeV.
For the calculation of LFV, first we write down the relevant MSSM sparticle masses at the M Z scale. The gaugino, sleptons and squark masses are given by [101]: m 2 Q 1,2 = 6.79m 2 1/2 +m 2 0 , m 2 U 1,2 = 6.37m 2 1/2 +m 2 0 , m 2 D 1,2 = 6.32m 2 1/2 +m 2 0 , (3.13) m 2 L 1,2,3 = 0.52m 2 1/2 +m 2 0 , m 2 E 1,2,3 = 0.15m 2 1/2 +m 2 0 . (3.14) Furthermore, the µ term is fixed from the symmetry breaking constraint: Now, the branching ratio of i → j γ is given by: The amplitude functions A L,R can be found in refs. [102][103][104]. The amplitudes are evaluated in the mass insertion approximation and we include loop contributions from both the neutralinos and charginos. For our model, the relevant amplitude is shows the same variation for different values of m 1/2 and A 0 fixed at 3.5 TeV. From these figures we see that part of the parameter space is excluded by the current experimental limit, and that future improvement can probe SUSY spectrum as large as 10 TeV. We note that the branching ratio for the decay τ → µγ is suppressed in our framework. When we choose parameters to satisfy the current limit on µ → eγ branching ratio, we obtain Br(τ → µγ) = 8×10 −11 , which is well below the experimental limit. The rare decay processes i → j γ provides the most stringent constraints on the parameter space of this class of models. Other LFV processes such as i → 3 j and µ−e transition in nuclei get the main contribution from the i → j γ penguin, with the photon attached to a lepton pair or a quark pair. Such diagrams are enhanced by a tan β factor compared to box diagrams, and will dominate. As a result, simple relations for their branching ratios can be derived [103]:  Br(µ − e in Ti) ≈ 6 × 10 −3 Br(µ → eγ) and Br( i → 3 j ) ≈ 7 × 10 −3 Br( i → j γ). Future experiments can probe µ − e transition in nuclei that will be sensitive to the level of O(10 −18 ). These models can be experimentally probed by future experiments in µ − e conversion in nuclei, µ → 3e decay as well as µ → eγ transition.

Discussion and conclusion
In this paper we have presented a U(1) PQ embedding of SUSY SO(10) models. We have adopted renormalizable SO(10) with a minimal Yukawa sector consisting of two matrices in flavor space. Such a Yukawa structure has been known to fit all of fermion masses and mixings, including neutrino oscillation data, in terms of a relatively small number of parameters. Without the PQ symmetry, these models would require a mini-split SUSY -25 -

JHEP06(2019)045
spectrum with gauginos having masses of order TeV and sfermions having masses in the 100 TeV range. Successful implementation of the PQ symmetry enables us to lower the SUSY spectrum, which now allows all SUSY particles to have masses of order TeV. Proton decay rates are suppressed in this scenario by a factor of (M PQ /M GUT ) 2 ∼ 10 −10 compared to the same model without the PQ symmetry. One feature of the PQ embedding is the appearance of an extra pair of Higgs doublets at the PQ scale. We have shown that unification of gauge couplings can be maintained by small threshold effects arising from color sextet scalars near the GUT scale. Although suppressed relative to SUSY SO(10) models without the PQ symmetry, proton lifetime is still within observable range. Both the Higgsino mediated d = 5 decay mode p → νK + and the gauge boson mediated d = 6 mode p → e + π 0 are within reach of ongoing and forthcoming experiments, with the rate for e + π 0 mode enhanced due to the lowering of the GUT scale in the PQ model compared to other GUTs.
The proposed model also predicts observable lepton flavor violation, notably in the decay µ → eγ and in µ − e conversion in nuclei. These processes originate from the Dirac neutrino Yukawa couplings which are effective in the momentum range 10 12 GeV ≤ µ ≤ 10 16 GeV, which transfer the LFV information to the sleptons and in turn to the leptons. The flavor structure of these LFV operators is completely determined in the model in terms of fermion masses and mixings.