A Predictive and Testable Unified Theory of Fermion Masses, Mixing and Leptogenesis

We consider a minimal non-supersymmetric $SO(10)$ Grand Unified Theory (GUT) model that can reproduce the observed fermionic masses and mixing parameters of the Standard Model. We calculate the scales of spontaneous symmetry breaking from the GUT to the Standard Model gauge group using two-loop renormalisation group equations. This procedure determines the proton decay rate and the scale of $U(1)_{B-L}$ breaking, which generates cosmic strings and the right-handed neutrino mass scales. Consequently, the regions of parameter space where thermal leptogenesis is viable are identified and correlated with the fermion masses and mixing, the neutrinoless double beta decay rate, the proton decay rate, and the gravitational wave signal resulting from the network of cosmic strings. We demonstrate that this framework, which can explain the Standard Model fermion masses and mixing and the observed baryon asymmetry, will be highly constrained by the next generation of gravitational wave detectors and neutrino oscillation experiments which will also constrain the proton lifetime.


Introduction
Grand Unified Theories (GUTs) have long been an attractive framework for unifying the non-gravitational interactions. The minimal option, which can predict neutrino masses and mixing, uses the gauge group SO (10). Several well-studied symmetries can be embedded in SO (10), including SU (5) [1], flipped SU (5) × U (1) [2][3][4][5] and the Pati-Salam model SU (4) c × SU (2) L × SU (2) R [6]. Thanks to this rich structure, there are many possible symmetry-breaking chains from SO(10) down to the Standard Model (SM) gauge group, G SM , most of them via the Pati-Salam symmetry [7]. An appealing feature of an intermediate Pati-Salam symmetry in non-supersymmetric GUTs is that gauge unification can be achieved, and there is an intermediate U (1) B−L subgroup which is spontaneously broken, generating right-handed neutrino masses. In addition to inducing light neutrino masses via the seesaw mechanism, the CP-violating and out-of-equilibrium decays of the right-handed neutrinos can produce the observed matter-antimatter asymmetry via thermal leptogenesis [8]. Moreover, the U (1) B−L symmetry breaking can also generate cosmic strings in the early Universe, which can intercommute and emit gravitational radiation forming a stochastic gravitational wave (GW) background that future GW interferometers can test.
The connection between GUTs and gravitational waves has been studied in [9] where the simple breaking pattern SO(10) → G SM × U (1) B−L → G SM was shown to be consistent with inflation, leptogenesis, and dark matter, while the U (1) B−L symmetry breaking generates cosmic strings. The connection between high-scale thermal leptogenesis and GWs was also pointed in [10] where it was assumed that the U (1) B−L breaking scale is the same as the seesaw and leptogenesis scales. In Ref. [11], we highlighted the complementarity between proton decay and gravitational wave signals from cosmic strings as a powerful method of probing GUTs. Subsequently, in Ref. [12], we studied all possible non-supersymmetric SO(10) symmetry-breaking chains. We performed a comprehensive renormalisation group (RG) analysis to find the correlations between the proton decay rate and the GW signal. We also identified which chains survived the current non-observation of both proton decay and GWs and could be tested by future neutrino and GW experiments.
In this paper, we go beyond these works by providing a detailed study on a specific SO (10) breaking chain that provides unification and predicts a proton decay width via the channel p → π 0 e + , consistent with the experimental bound of the Super-Kamiokande (Super-K) [13] and can be fully tested by future proton decay searches of Hyper-K [12]. Further, this breaking chain generates cosmic strings at the lowest intermediate scale, M 1 ∼ 10 13 GeV. A GW background generated by such a string network is just around the corner and may be even already hinted at by recent observations in PTA experiments, including NANOGrav [14], PPTA [15], EPTA [16] and IPTA [17]. We determine the minimal necessary particle content to induce the pattern of breaking and perform an RG analysis and numerical fit of our model to SM data to postdict the fermion masses and mixing, including the mass scales of the right-handed neutrinos. As this procedure determines the scales of symmetry breaking of our model and the masses of the right-handed neutrinos, the matter-antimatter asymmetry associated with thermal leptogenesis is predicted. We then show that successful leptogenesis can occur in the regions of the model parameter space consistent with SM fermion masses and mixing and can be correlated with a GW signal and proton decay. Compared with [10], such an approach allows us to go beyond generic considerations and instead to quantitatively account for the hierarchy between the leptogenesis and see-saw scales, as well as with the U (1) B−L breaking scale, thanks to the constraints imposed by reproducing the low energy data. The latter scale is of particular interest since pulsar timing arrays such as PPTA [18] and NANOGrav [19] are sensitive to the predicted GW signals while future large-scale neutrino experiment, Hyper-Kamiokande (Hyper-K) [20], will be able to probe the predicted proton decay rate of this model. The correlation between these two observables will be a crucial test of our GUT model, and such methodology can be applied to other GUT models, presenting a new avenue to try to unveil the physics at very high scales. This paper is organised as follows: in Section 2, we discuss the GUT symmetry breaking pattern and the particle content of our model, including fermionic and Higgs representations of the GUT and our RG running procedure. In Section 3, we discuss how we relate our model to the quark lepton data, our fitting procedure and the ensuing results. In Section 4, we discuss the basics of non-resonant thermal leptogenesis and how we determine the baryon asymmetry produced from the successful points in the model parameter space and in Section 5, we demonstrate that the regions of the model parameter space that yield successful leptogenesis and fermionic masses and mixing will be associated with a GW signal.
Finally, we summarise and discuss in Section 6. As a case study, we consider a benchmark point (referred to as BP1 throughout) and discuss how it satisfies all these experimental constraints in each section.

The framework
We focus on a breaking chain (classified as chain III4 of type (c) in Ref. [12]) that is of particular interest as it is currently allowed and predicts a proton decay rate testable by Hyper-K. We discuss the breaking chain's matter content and gauge unification in this section.

Symmetry breaking of SO(10)
We study the following breaking chain with three intermediate symmetries (G 3 , G 2 , and G 1 ): The boldface number beside the arrow indicates the Higgs representation of SO (10), triggering the symmetry breaking. In this work, we follow the same convention as Ref. [12] where the GUT symmetry breaking scale is denoted as M X and the mass scale of the subsequent breaking of the group G I (for I = 1, 2, 3) is denoted as M I . 1 All particles, except the gauge fields of the model, are listed in Table 1. We note that Z C 2 refers to the parity symmetry between left and right conjugation (L ↔ R c , where c indicates charge conjugation) and that . The correlations between U (1) charges are given by Y = To achieve each step of symmetry breaking, i.e., SO(10) → G 3 → G 2 → G 1 , we include three Higgs multiplets, 54, 210, and 45 of SO(10), respectively. These Higgs fields spontaneously break the GUT symmetry as follows: • 54 contains a parity-even singlet (1, 1, 1) of G 3 ≡ SU (4) c × SU (2) L × SU (2) R where each entry in the bracket (r 1 , r 2 , · · · ) refers to the field representation transforming in the group G ≡ H 1 × H 2 × · · · . Once (1, 1, 1) gains a non-trivial vacuum expectation value (VEV) at scale M X , SO(10) is spontaneously broken to G 3 .
• In G 3 , 210 can be decomposed to (15, 1, 1) 1 of G 3 , which is further decomposed into a parity-even and trivial singlet (1, 1, where the last entry in the bracket is the field charge in the U (1) symmetry and the subscript is used to distinguish from another field with the same representation discussed below. The VEV of this singlet breaks G 3 to G 2 at scale M 3 .
• The breaking of G 2 to G 1 is realised via a 45 of SO(10), which decomposed into (15, 1, 1) 2 of G 3 and a further (1, 1, 1, 0) 2 of G 2 and G 1 . This singlet is parity-odd, and its VEV induces the breaking of G 2 → G 1 at scale M 2 .
We summarise the decomposition of Higgses, which triggers the breaking of SO(10) and intermediate symmetries, in Table 2.

Matter Field Decomposition and Fermion Masses
In order to assess if our model can predict the measured fermionic masses and mixing, it is important to understand the matter content of the breaking chain. Fermions are arranged as a 16 of SO(10) and follow the decomposition given in Table 3 where L (R) denote the left-handed (right-handed) fermions of G 3 which contains the SM left-handed (righthanded) fermions where Q L(R) and L(R) are the quark and leptonic SU(2) L(R) doublets, respectively, and u R , d R , e R , and ν R are the quark and lepton SU(2) L singlets, respectively.
Three Higgs multiplets, 10, 126 and 120, are required to generate the Standard Model fermion masses. Compared to Ref. [12], where we considered a minimal survival hypothesis [21], we include one additional Higgs (120) which is required to generate all fermion mass spectra, mixing angles, and CP-violating phases in the quark and lepton sectors. Here, 126 is the same Higgs used in the breaking G 1 → G SM . For this breaking chain, we list decompositions of Higgs, which are responsible for mass generation in Table 2.
and we assume that the heavy one gains a mass ∼ M X and thus decouples at scales below M X . The same assumption applies to the other two (1, 2, 2)'s. Using these assumptions, we have two bi-doublets (1, 2, 2, 0) at the scale M 3 > Q > M 2 and M 2 > Q > M 1 , where G 2 and G 1 are preserved, respectively. We retain them as the physically relevant degrees of freedom in this range of scales, following the logic of Ref. [12]. Four electroweak doublets remain at energies below M 1 but above the electroweak scale. Naively, one can assume all massive states are sufficiently heavy that they decouple at scale M 1 , except for the lightest electroweak doublet, which is the SM Higgs and should be massless before electroweak symmetry breaking. Without loss of generality, we can write these Higgses as superpositions of mass eigenstates,ĥ i = j V ij h j , with h SM ≡ĥ 1 , where V is a unitary matrix and the heavy doublets that decouple at M X have also been taken into account. With this treatment, all physical degrees of freedom present at the relevant scale are the same as those of chain III4 in Ref. [12]. For the second Higgs multiplet, 126, we retain another decomposed representation (10, 3, 1) of G 1 , which contains a SU (2) L triplet, (1, 1, 3, −1), of G 2 and G 1 which contains the singlet S ∼ (1, 1, 0) of G SM that is important not only in its role in symmetry breaking, but also in the generation of neutrino masses. (10, 3, 1) is retained due to the requirement of left-right parity symmetry, Z C 2 , and it is decomposed to a (1, 3, 1, 1) of G 2 . After G 2 breaking, i.e., the breaking of the left-right parity symmetry, we assume that this particle decouples.
In the non-SUSY case, two further couplings 16 · 16 · 10 * and 16 · 16 · 120 * are allowed by the gauge symmetry; however, we forbid them by imposing an additional Peccei-Quinn U (1) symmetry [22] as described in [23][24][25]. After the final symmetry is broken to G SM , the above Yukawa terms generate the following SM fermion mass terms in the left-right convention: Rotating the Higgs fields to their mass basis, we derive Yukawa couplings to the SM Higgs as A Majorana mass term for the right-handed neutrinos is generated from the second term of Eq. (2.2): once φ S acquires a VEV, v S , which controls the scale of the masses: After the right-handed neutrinos decouple and electroweak symmetry is broken, the light neutrinos acquire their mass via the Type-I seesaw mechanism [26][27][28][29]: where the SM Higgs VEV is v SM = 175 GeV. We emphasise that the electroweak singlet, φ S , is essential for the symmetry breaking G 1 → G SM and thus, its VEV determines the scale of M 1 and the right-handed neutrino masses. As required by perturbativity, Y 126 O(1), the mass of the heaviest right-handed neutrino, M N 3 , should be not heavier than the lowest intermediate scale, M 1 . On the other hand, neutrino oscillation experiments have given relatively precise measurements of light neutrino masses and mixing angles. These data restrict the right-handed neutrino mass spectrum via the seesaw formula, and a realistic GUT model should survive all such constraints.

Gauge Unification
Given an arbitrary gauge symmetry G, which can be expressed as a product of simple Lie groups, G = H 1 × · · · × H n , the two-loop renormalisation group running equation for group H i , for i = 1, 2, · · · , is given by where α i = g 2 i /(4π) and the β function is determined by the particle content of the theory: Here, i ∈ [1, · · · , n] for H n , g i is the gauge coefficient of H i , and b i and b ij refer to the normalised coefficients of one-and two-loop contributions, respectively. In the following, we neglect the Yukawa contribution to the RG running equations as it gives a subdominant contribution. Given two scales Q 0 and Q, if the conditions Q 0 < Q and b j α j (Q 0 ) log(Q/Q 0 ) < 1 are both satisfied then an analytical solution for these equations can be obtained [30]: In the case that both H i and H j are non-abelian groups, the coefficients b i and b ij are where the ψ and φ indices sum over the fermions and complex scalar multiplets, respectively, and ψ i and φ i are their representations in the group H i , respectively. C 2 (R i ) (for R i = ψ i , φ i ) denotes the quadratic Casimir of the representation R i in group H i and C 2 (H i ) is the quadratic Casimir of the adjoint presentation of the group H i . In particular, for SU (N ), C 2 (SU (N )) = N and the quadratic Casimir of the fundamental irrep N of SU (N ) is given by C 2 (N) = (N 2 − 1)/2N ; for SO(10), C 2 (SO(10)) = 8, and the quadratic Casimir of the fundamental irrep 10 of SO(10) is given by C 2 (10) = 9/2. The spinor representation of SO(10), 16, has C 2 Explicit values of b i and b ij depend on the degree of freedoms introduced by the gauge, matter and Higgs fields. The gauge fields are directly determined by the gauge symmetry in the breaking chain. In regards to the matter fields, we assume they are the minimal extension which includes all the SM fermions, i.e., minimally a 16 of SO(10) as in Table 3. The most significant uncertainty contributing to RG running comes from the Higgs sector as one has to account for all the Higgses used to generate fermion masses and the GUT and intermediate symmetry breaking. Given the decomposition of Higgs fields in Table 2 and the discussion in Section 2.2, the Higgs fields included in each step of the RG running are: • For G 1 → G SM , we include only the SM Higgs. Although we arrive at a series of electroweak doublets after field decomposition, we assume that all degrees of freedom except the SM Higgs are sufficiently heavy that they are integrated out by this breaking step and thus have a negligible effect on the RG running.
• For G 2 → G 1 , we include three Higgses in the running, two (1, 2, 2, 0)'s and one (1, 1, 3, −1) of G 1 . The former includes the SM Higgs, and the latter includes the gauge singlet φ S of G SM which is used to achieve the breaking of G 1 → G SM and right-handed neutrino masses.
in the RG running. Two further Higgses are included compared to the above item as (1, 3, 1, 1) is required for the matter parity symmetry Z C 2 and (1, • For SO(10) → G 3 , we include (1, 2, 2), (15, 2, 2), (10, 1, 3), (10, 3, 1) and two (15, 1, 1)'s in the RG running. The former two are required to obtain the two (1, 2, 2, 0)'s above. (10, 1, 3) and (10, 3, 1) are required for (1, 1, 3, −1) and (1, 3, 1, 1). One (15, 1, 1), decomposed from 45 is for (1, 1, 1, 0) 2 , and the other, decomposed from 210, includes the singlet (1, 1, 1, 0) 1 to achieve the breaking G 3 → G 2 . By including the above particle content in the RG running, we obtain the coefficients b i and b ij at the two-loop level, which we list in Table 5 and are the same as in the chain III4 of Ref. [12]. Although we include one more Higgs multiplet 120, the contribution of induced new particles can be ignored, as explained in the previous subsection, by assuming heavy mass eigenstates heavier than the breaking scale, M X [31]. In order to keep the treatment of the RG running economical, the scalar multiplets which are unnecessary for the breaking chain are assumed to be as massive as the SO(10) breaking scale M X . Therefore, these scalars will not affect the RG running or provide threshold corrections.
During the symmetry breaking at an intermediate scale (M 3 , M 2 or M 1 ), gauge couplings of the larger symmetry and those of the residual symmetry after spontaneous symmetry breaking (SSB) must satisfy matching conditions. Here we list one-loop matching conditions that appear in the GUT breaking chains. For a simple Lie group H i+1 broken to subgroup H i at the scale Q = M I , the one-loop matching condition is given by [32] Table 5: Coefficients b i and b ij of gauge coupling β functions appearing in the specified breaking chain.
For G 1 → G SM , we encounter the breaking, SU (2) R × U (1) X → U (1) Y , which has the matching condition [33]: Applying the matching conditions of the above two equations, all gauge couplings of the subgroups unify into a single gauge coupling, α X ≡ g 2 X /4π, of SO(10) at the GUT scale, M X . This condition restricts both the GUT and intermediate scales for each breaking chain. We denote the mass of the heavy gauge boson masses associated with SO(10) breaking as M X and M 3 , M 2 and M 1 are associated to the breaking of G 3 , G 2 and G 1 , respectively. Correlations among M 1 , M 2 , M 3 and M X are determined numerically using the following procedure for the breaking chain SO(10) 1. Begin the evaluation from the scale M Z with the SM gauge couplings α 3 = 0.1184, α 2 = 0.033819 and α 1 = 0.010168 [34]. Evolve these couplings using the RGE of the SM to scale M 1 , where G 1 is recovered. Apply the matching conditions for the SM 3. Repeating this same procedure, to evolve all couplings to the GUT scale, M X , to unify to a single value α X with the matching condition at M X fully accounted for.
The above RG running procedure involves four scales M 1 , M 2 , M 3 and M X . Gauge unification requires that three SM gauge couplings meet each other at the GUT scale, up to matching conditions, and enforces two constraints; thus, there are only two free scales. The remaining scales and the gauge coupling, α X , are then determined via gauge unification. General restrictions on the parameter space of scales i.e., M 2 , M 3 and M X varying with M 1 , are shown in the left plot of Fig. 1. Given the gauge unification scale, M X , and its gauge coupling at that scale, the proton lifetime via the decaying process p → π 0 e + is predicted. Due to the scale correlations imposed by gauge unification, the correlation between the proton lifetime, τ p , and M X can be transformed into a correlation between τ p and any intermediate scale. Following the formulation of Ref. [12], we derive the allowed parameter space of τ p versus intermediate scales. By varying the lowest intermediate scale M 1 , we obtain general regions of τ p . The bound of τ p versus the lowest scale M 1 is shown in the right panel of Fig. 1. The Super-K experiment set a lower bound on the proton lifetime, τ (p → e + π 0 ) > 2.4 × 10 34 years at 90 % confidence level [13]. In the future, Hyper-Kamiokande (Hyper-K) is expected to improve the measurement of proton lifetime by almost one order of magnitude [20]. If proton decay is not observed, the entire parameter space of this breaking chain will be excluded.  This benchmark point will be considered throughout this paper. Its associated proton decay rate, τ (p → e + π 0 ) ∼ 5.1 × 10 34 years, is consistent with the current Super-K bound and will be tested by Hyper-K. We note that BP1 is consistent with SM fermion masses and mixing to a high statistical significance, and this requires M 1 ∼ 10 13 GeV. Such a high value for M 1 leads a compressed hierarchy between M 1 , M 2 and M 3 and this comes from the constraint of gauge unification (from the left panel of Fig. 1 this region is in the right corner of the blue triangle.)

Fermion masses and mixing
As all the SM fermions are embedded in the same SO(10) multiplet (16), their masses are correlated with each other. Therefore, it is a non-trivial task to find regions of the GUT model parameter space that predict the SM fermion masses and mixing consistent with the precisely measured (particularly in the quark sector) experimental data. This section presents the correlations of masses and mixing between quarks and leptons and predicts heavy neutrino masses using the model we discussed in the previous section. We parametrise the up, down, neutrino, charged lepton Yukawa couplings and right-handed neutrino mass matrix, respectively, as follows [35]: 2) and V ji denotes the mixing between the mass and interaction basis of the electroweak Higgs doublets. The light neutrino mass matrix, M ν , is obtained by

Parametrisation using Hermitian Yukawa matrices
The most general form of Yukawa couplings and neutrino mass matrix includes many free parameters. A considerable reduction in the number of parameters can be achieved by considering only the Hermitian case for all fermion Yukawa couplings matrices Y u , Y d , Y ν and Y e (and M R should be real as a consequence of the Majorana nature for right-handed neutrinos). Such a reduction can result from spontaneous CP violation [36,37] which assumes that there exists a CP symmetry above the GUT scale, leading to real-valued Y 10 , Y 126 and Y 120 , and the CP is broken by some complex VEVs of Higgs multiplets during GUT or intermediate symmetry breaking. For the particular chain we applied in the last section, one can consider, for example, the parity-odd singlet of , gains a purely imaginary VEV. Then, via couplings such as 45 · 10 · 120 (and 45 · 126 · 120) which generate purely imaginary offdiagonal mass terms between h u,d 10 and h u,d 120 (and those between h u,d 126 and h u,d 120 ) and further purely imaginary mixing entries V 13 , V 14 (and V 17 , V 18 ) are obtained. As a result, h, f and h , as well as all parameters on the right-hand side of Eq. (3.1), are real. Since h is antisymmetric, we arrive at Hermitian Dirac Yukawa coupling matrices Y u , Y d , Y ν and Y e . This texture has widely been applied in the literature, see e.g., Refs. [23,31,38]. The resulting fermion mass matrices conserve parity symmetry L ↔ R [31] and following from the assumption that there is no CP violation in the Higgs sector, apart from that of 120, r 1 , r 2 , r 3 , c e , and c ν are all real parameters resulting in a real symmetric right-handed neutrino mass matrix, M ν R . The CP symmetry in the Yukawa coupling is spontaneously broken after the Higgses gain VEVs.
For simplicity, we assume that r 3 = 0, which implies that the imaginary part of Y u vanishes. It is convenient to write the up-type Yukawa in the diagonal basis which can be achieved via a real-orthogonal transformation on the fermion flavours without changing the Hermitian property of Y d , Y e , and Y ν . In the above, η u,c,t = ±1 refer to signs that cannot be determined by the real-orthogonal transformation. While η t = +1 can be fixed by making an overall sign rotation for all Yukawa matrices, the remaining signs, η u and η c , cannot be fixed and are randomly varied throughout our analysis. In the basis of the diagonal up-quark mass matrix, Y d is given by where s ij = sin θ q ij , c ij = cos θ q ij and P a = diag{e ia 1 , e ia 2 , 1}. The matrices h, f and h are then expressed in terms of Y u and Y d The light neutrino mass matrix can be expressed as Using this parametrisation, all six quark masses and four CKM mixing parameters are treated as inputs, and we are then left with seven parameters (a 1 , a 2 , r 1 , r 2 , c e , c ν , and m 0 ) to fit eight observables, including three Yukawa couplings y e , y µ , y τ , two neutrino masssquared differences ∆m 2 21 , ∆m 2 31 and three mixing angles θ 12 , θ 13 , θ 23 , where the leptonic CP-violating phase, δ, will be treated as a prediction 2 .

Procedure of numerical analysis
This section describes how we identify regions of our model parameter space consistent with fermion masses and mixing while evading the existing proton decay limit. In our numerical analysis, we use the following experimental data: • We fix the Yukawa couplings (y) of charged fermions and CKM mixing angles (θ) at their best-fit (bf) values [23,39,40] y bf u = 2.54 × 10 −6 , y bf c = 1.37 × 10 −3 , y bf t = 0.43 , y bf d = 6.56 × 10 −6 , y bf s = 1.24 × 10 −4 , y bf b = 5.7 × 10 −3 , y bf e = 2.70 × 10 −6 , y bf µ = 5.71 × 10 −4 , y bf τ = 9.7 × 10 −3 ,  3 However, as we will later see, fixing them at the best fit values is sufficient to reproduce all mixing data. Thus, in this current discussion, we will ignore them for simplicity.
• In the neutrino sector, we use the best-fit values from NuFIT 5.1 [41] and include the 1σ uncertainty. Those data with and without Super-K atmospheric data are, respectively, given by The atmospheric mixing angle, θ 23 , is restricted to first octant (0 < θ 23 < 45 • ) and the second (45 • < θ 23 < 90 • ), respectively, in the two cases. In both cases, normal ordering (i.e., m 1 < m 2 < m 3 ) of neutrino masses is assumed. Inverted ordering (i.e., m 3 < m 1 < m 2 ) will not be discussed as a preliminary scan indicates that our model does not favour the inverted ordering. We do not consider the small flavour-dependent RG running effect due to the suppression of charged lepton Yukawa coupling.
The statistical analysis is performed in the following way: • As quark masses and mixing parameters are fixed at their best-fit values, Y u is fully determined except for the signs of η u and η d (note that η t = +1 is fixed by an overall sign rotation). Y d depends on two free model parameters, a 1 and a 2 , and signs (η d , η s , η b ).
• Based on Eq. (3.7), Y e depends on the two phases a 1 , a 2 and three ratios r 1 , r 2 , c e up to the above sign differences. Note that Y e must satisfy three equations simultaneously: and as the right hand side is fixed, r 1 , r 2 and c e are fully determined by the phases a 1 , a 2 and the signs η q (for q = u, c, d, s, b). We scan the phase parameters in the range a 1 , a 2 ∈ [0, 2π] and vary the signs η q = ±1 randomly and solve for r 1 , r 2 and c e . Then, we substitute these values into Eq. (3.7) and determine the unitary matrix V e used in the diagonalisation V † e Y e Y † e V e = diag{y 2 e , y 2 µ , y 2 τ }.
• In Eq. (3.8), the neutrino mass matrix, M ν , is determined by two further parameters c ν and m 0 . The former determines the flavour structure and the latter the absolute mass scale, and by scanning these parameters, we determine M ν . The diagonalisation V † ν M ν V * ν = diag{m 1 , m 2 , m 3 } provides the neutrino mass eigenvalues and unitary matrix V ν .
• The PMNS matrix is given by U PMNS = V † e V ν , and the three leptonic mixing angles are derived via In summary, once the charged fermion masses and quark mixing parameters are fixed, we are left with only four free model parameters a 1 , a 2 , c ν , m 0 and signs η q : We scan the model parameter space, P m , to fit five observables: O n ∈ {θ 12 , θ 13 , θ 23 , ∆m 2 21 , ∆m 2 31 } . In this way, we efficiently reduce the dimensionality of the parameter space from 17 to 5 dimensions. Following the above simplified treatment, we scan two phases a 1 , a 2 in the range [0, π]. The coefficient |c ν | is logarithmically scanned in the range [10 −3 , 10 3 ], and we randomly assign its ± sign. m 0 (meV) is solved by minimising the χ 2 function, which is used as a measure of how well our model fits the data, being defined as (3.17) Given the predefined theory model parameter space, P m , and scanning in the relevant ranges of these parameters, we determine which regions fit the experimental data by setting an upper bound of χ 2 value. This procedure of the scan is divided into two steps: We first perform a preliminary scan by setting the upper bound of χ 2 < 100 and then perform a subsequent scan to find the points with χ 2 < 10. The results of the first scan which uses the neutrino oscillation data of Eq.
where the mass states of ν R , from the lightest to heaviest, are denoted as N 1 , N 2 and N 3 . We impose an upper bound by requiring M N 3 M 1 , and this is approximately equivalent to requiring that the largest eigenvalue of Y 126 1 such that the perturbativity is respected.
Since the maximal value of M 1 allowed by proton decay measurements is given by 4.4 × 10 13 GeV [12], viable points in the model parameter space require that Naively, by assuming the magnitude of the Dirac Yukawa coupling Y ν ∼ O(1), we know from the seesaw formula that the RHN mass scale is around 10 15 GeV. Thus, one can expect that the condition of Eq. (3.19) rules out most points. This is confirmed by the bottom-left panel Fig. 3 where most of the points predict the heaviest neutrino mass, M N 3 , to be heavier than 4.4 × 10 13 GeV. Therefore, these points are not consistent with the requirement of gauge unification. We then perform a second more dense scan around the former points by requiring χ 2 < 10 and gauge unification, e.g., the bound of the heaviest right-handed neutrino mass satisfying Eq.
where θ C is the Cabibbo angle. Taking these relations into account, we derivẽ and from Eq.
which parametrises small deviations from 1 due to the suppressions of y µ /y τ and y s /y τ . We keep this deviation due to the fact that y τ /y b ≈ 1. We rewrite it in the form  which provides an almost linear correlation between r 1 and r 2 , which is clearly shown in Fig. 5 for both |ỹ τ −ỹ b | = y τ − y b and y τ + y b . To satisfy the experimental constraint that y τ /y b 1, most of the points with low χ 2 values satisfy the following hierarchical relation: namely that r 2 + 3/r 2 − 1 ≈ 1 and |4r 1 /r 2 − 1| 1. Using Eqs. (3.20) and (3.24), one can further derive approximate expressions for y µ and y e : y µ sin 2 θ q 13 sin 2 θ q 23 (cos a 1 cos a 2 + c 2 e sin a 1 sin a 2 ) 2 + c 2 e sin 2 (a 1 − a 2 ) , where a 1 = a 1 − δ q . However, the relation shown in Eq. (3.24) is not always satisfied as a few points with r 2 ∼ 0.13 are also allowed by χ 2 < 10 in our scan and this leads to (r 2 + 3)/(r 2 − 1) ≈ −3.6.

Inputs
which neutrinoless double-β decay experiments can determine if neutrinos are Majorana as predicted by our model. m 1 is predicted to be less than 2 meV and m ββ ∼ (2, 20) meV. The above numerical analysis applies to only the normal ordering of light neutrino masses. We completed a preliminary scan for the inverted ordering and found that very few points had χ 2 < 100 and did not proceed with a more dense scan. Although we have not proved that the GUT model disagrees with the inverted mass ordering, our model shows a preference for normal mass ordering given our scan.

The benchmark study
In this subsection, we study BP1 focusing on the flavour sector of quarks and leptons. Inputs and predictions of fermion Yukawas and mixing parameters are shown in Table 6 from which we obtain From this, the Yukawa and neutrino mass matrices are obtained: Diagonalisation of Y e and M ν gives rise to the lepton masses and mixing, and the above benchmark provides χ 2 = 0.33. We show predictions of charged lepton Yukawa couplings, mixing angles and the Dirac CP-violating phase, as well as the lightest neutrino mass m 1 and m ββ in Table 6. Applying the type-I seesaw formula, the right-handed neutrino mass matrix is Moreover, the three associated eigenvalues are given in Table 6. Finally, we note that the heaviest right-handed neutrino mass is given by 1.6 × 10 13 GeV. This is lower than the lowest intermediate scale M 1 and thus consistent with proton decay measurement.

Leptogenesis
As SO(10) GUTs predict very massive right-handed neutrinos, leptogenesis is a natural consequence of such a framework. SO(10) leptogenesis was initially studied in [42,43]. Later, the importance of flavour effects was emphasised [44] and in [45] the same authors investigated the correlations between the viable SO(10) leptogenesis parameter space and low-scale observables. Leptogenesis within a specific SO(10) model, which fit low energy SM fermionic data [46], was investigated in [47]. Fitting leptogenesis together with all the fermion mass observables by solving density matrix equations was carried out in [48]. In this section, we calculate the baryon asymmetry produced via thermal leptogenesis for the points of the parameter space scan, which are consistent with the quark and lepton experimental data with χ 2 < 10 and later emphasise the connection between SO(10) leptogenesis and observables such as proton decay and GWs.
To calculate the associated baryon-to-photon ratio, we solve the Boltzmann equations which determine the time evolution of the lepton asymmetry that manifests from the CP-violating and out-of-equilibrium decays of the right-handed neutrinos and associated washout processes. In the simplest formulation, these kinetic equations are in the one-flavoured regime, in which only a single flavour of charged lepton is accounted for. This regime is only realised at sufficiently high temperatures (T 10 12 GeV) when the rates of processes mediated by the charged lepton Yukawa couplings are out of thermal equilibrium. Therefore, a single charged lepton flavour state is a coherent superposition of the three flavour eigenstates. However, if leptogenesis occurs at lower temperatures (10 9 T 10 12 GeV), scattering induced by the tau Yukawa couplings can cause the single charged lepton flavour to decohere, and the dynamics of leptogenesis must be described in terms of two flavour eigenstates. In such a regime, a density matrix formalism [49][50][51][52][53] allows for a more general description than the one-flavoured semi-classical Boltzmann equations. We begin by rotating the Yukawa coupling of the leptonic and Higgs doublet to the right-handed neutrinos mass basis. For example, the benchmark case in Eq. (3.30) is rotated to (4.1) The CP-asymmetry matrix, describing the decay asymmetry generated by N i is denoted by αβ , and may be written as [54]: where x i ≡ M 2 N i /M 2 N 1 and Greek and Roman indices denote charged lepton flavour and right-handed neutrino generation indices, respectively, and N B−L is a density matrix parametrising the lepton asymmetries in flavour space, and P (i)0 denotes the projection matrices which describe how a given flavour of lepton is washed out: Finally, the density matrix equations which track the time evolution of the lepton asymmetry generated by the decays and washout of the three right-handed neutrinos are given by: where z = M N 1 /T is the evolution parameter, H is the Hubble expansion rate, and the decay and washout terms are given by where z i ≡ √ x i z, K 1 and K 2 are modified Bessel functions of the second kind with the decay asymmetry K i given by respectively. The thermal widths of the charged leptons, Λ τ , Λ µ , are given by the imaginary part of the self-energy correction to the charged lepton propagators in the plasma, and this mediates flavour correlations (for further details of these equations and their application we refer the reader to Ref. [56]). For each point in our parameter scan, we have solved Eq. (4.5) which provides the baryon-to-photon ratio using the publicly available tool ULYSSES [57] and the associated "3DME" code which accounts for the decays and washout of all three right-handed neutrinos. We show our results in Fig. 6 where all the points have χ 2 < 10 and satisfy the proton decay bounds. We observe that for both octants, many low χ 2 points achieve thermal leptogenesis successfully, and their baryon-to-photon ratio is η B ≈ 10 −10 . Interestingly, the predicted baryon-to-photon ratio shows little dependence on δ but has a very constrained prediction for the effective Majorana mass, 4 m ββ (meV) 10. This indicates that the predicted Majorana phases are highly constrained within our model.
For all the points with χ 2 < 10, we found leptogenesis is always in the strong washout regime K 1 1 (for the benchmark point K 1 ≈ 130) since the Yukawa couplings (Eq. (3.30)) are not very small. As we are in the strong washout regime, neglecting finite density effects is a good approximation [58]. For the above benchmark point, the lightest two masses of right-handed neutrinos are M N 1 = 4.23 · 10 11 GeV and M N 2 = 5.32 · 10 11 GeV and the baryon-to-photon ratio is ∼ 6.11 × 10 −10 ). Next generation 0νββ experiments such as . The colour of the points denotes the ratio of the predicted baryon-to-photon ratio to the experimentally observed best-fit value as measured using CMB data η CMB B = 6.15 × 10 −10 [55]. Consistency with gauge unification is considered. In the leftmost plots, the dashed line labels the sensitivity of the next generation experiments on 0νββ decay.
Legend-1000 [59], nEXO [60], NEXT-HD [61], DARWIN [62], SNO+II [63] and CUPID-Mo [64] will allow to test further test the results of our scan. The most sensitive experiments will probe m ββ down to 8 − 10meV which, as we can see from Fig. 6, covers a significant fraction of the points of our scan.

Testability of SO(10) GUT and leptogenesis using gravitational waves
The creation of topological defects in GUT symmetry breaking is ubiquitous [7] and in our breaking chain produces unwanted defects in the breaking of SO(10) → G 2 (monopoles and domain walls), but in the final breaking, G 1 → G SM , cosmic strings are formed. We assume that inflation occurs after the unwanted defects are formed but before string formation. 4 Under this assumption, the string network can intersect to form loops which oscillate and emit energy gravitational radiation that can constitute a stochastic gravitational wave background. Importantly, this background can, in principle, be observed by currently running and future GW experiments.
We assume the Nambu-Goto string approximation, where the string is infinitely thin with no couplings to particles [68], and the amplitude of the relic GW density parameter is: where ρ c is the critical energy density of the Universe and ρ GW depends on a single parameter, Gµ where G = M −2 pl is Newton's constant and µ is the string tension. For strings generated from the gauge symmetry Gµ is approximately given by [69]: and hence we can relate the string tension parameter to the lowest intermediate scale of GUT symmetry breaking. The spectrum of Ω GW (f ) is calculated numerically by solving the formulation in [70] and [11,12]. A well-known feature is at the high-frequency band, Ω GW (f ) has a Gµ-dependent plateau [71]. In our calculations, we derive a semi-analytical correlation where numerical values α 2R ∼ α 1X ∼ 1/40 around the scale M 1 have been accounted for. A GW background with (Ω GW h 2 ) plateau ∼ 10 −10 is naturally predicted from cosmic strings from the breaking scale M 1 ∼ 10 13 GeV. 5 There is an important relationship between the string loop length (l) at the time of formation (t i ), which is usually taken as a linear relation: l = αt i . While α has a distribution of values, it peaks at around 0.1 for both matter and radiation-dominated eras [72].
Gµ, or equivalently M 1 , should be restricted to a specific range as its value must be consistent with proton decay measurement and fermion masses and mixing data. For the breaking chain we consider, the upper bound is M upper 1 4.4 × 10 13 GeV from the proton decay. We account for a perturbativity ansatz for Majorana Yukawa coupling, i.e., The prediction of Ω GW h 2 depends on the loop size parameter α. While simulations [71,72] show this parameter peaks around α ≈ 1. Deviation could be considered as there could be some uncertainties in the peak value and distribution around it. A smaller value of α reduces the GW intensity as Ω GW ≈ √ α. In the plot, we consider the case with the value α = 0.01 (light red) as a comparison to the main result derived at α = 0.1 (red). To obtain the bound explicitly, one has to consider the RG behaviour of Yukawa couplings and include the influence of Yukawa couplings on gauge unification. This procedure is complicated, and we relegate this for future work. In this work as an alternative, we consider a simplified perturbative ansatz by requiring the heaviest right-handed neutrino mass lighter than the U (1)B−L scale, i.e., MN 3 < M1.  Figure 8: Comparison between SGWB signals produced by cosmic strings with Gµ from 10 −11 to 10 −9 the possible 1σ and 2σ regions hinted by EPTA, PPTA, IPTA and NANOGrav. The SGWB signal has been fitted for three frequencies: 2.4 nHz, 5.4 nHz and 12 nHz. The simulation shows that the signal is compatible with all experiments at 2σ. The orange star indicates the prediction of (γ, A) in BP1.
Data from pulsar timing arrays (PTA) such as EPTA [83] and NANOGrav (11-year data set) [19] have already probed the nHz regime and provided upper limits through the nonobservation of GWs. Future PTA such as SKA [84] will cover even more parameter space, and large-scale surveys of stars such as Gaia [85] and the proposed upgrade, THEIA [86], can be powerful and complementary probes of gravitational waves in the same frequency regime.
The PTA experiment NANOGrav released its 12.5-year data set [14] which announced the detection of a common-spectrum process that has a characteristic strain of the form where f yr = 1/year, and A is an amplitude parameter referring to the correlation between pulsars. If this signal is confirmed as a cosmic GW background, the characteristic strain can be transformed into the GW energy density: (IPTA) , respectively. These data hint at a GW energy density Ω GW h 2 ∼ 10 −10 in the nHz band. The hinted ranges of A and γ are highly correlated in all these experiments. We include 1σ and 2σ regions of (γ, A) restricted by these measurements in Fig. 8. In the same figure, we also present three dotted curves which show how the Gµ, associated with cosmic strings, overlays on the A-γ regions of various experiments. Following the treatment used in [88] the curves are obtained as follows: for a given value of Gµ, the spectrum of Ω GW (f )h 2 is obtained from our numerical calculation, A and γ are reversely calculated via Eq. (5.5) by fixing the frequency f = f * . In particular, γ is solved via γ = 5 − ∂ log Ω GW (f ) ∂ log f f * . We set f * = 2.4, 5.6 and 12 nHz, respectively, vary Gµ from 10 −11 to 10 −9 , and obtain three curves respectively, shown in Fig. 8 where the second frequency is the one considered in [88] and the first and third refer to the first and the fifth lowest frequencies measured in NANOGrav12.5. The overlap between these curves and the PTA hinted regions indicate the consistency between cosmic-string-sourced SGWB signal and the experimental searches. For f * = 5.6 nHz, we found that the signal is compatible with NANOGrav12.5, PPTA and IPTA in 1σ and EPTA in 2σ ranges. A and γ in 2σ ranges are restricted in respectively. Here, upper bounds of A and γ in NANOGrav and IPTA are not considered as they predict Ω GW h 2 > 10 −9 which are too large for our interests. We apply these ranges to Eq. (5.5), and in the left panel of Fig. 9 we show a zoom-in of the upper bound and BP1 GW signal intersecting the various PTA experimental sensitivities. We observe that our benchmark point can explain the possible signal observed. In the right panel of Fig. 9 we show our GUT model's predicted proton lifetime as a function of M 1 .
In the left panel of Fig. 10, we show the range of M 1 values predicted by our model, where the yellow shows the region consistent with the NANOGrav 12.5-year data set. As there is a degree of freedom which relates M 1 and M N 1 that we cannot constrain, there is a range of lightest right-handed neutrino masses that can accommodate this GW signal. In the right panel, we show how our leptogenesis prediction lies in the M N 1 − M N 3 plane.
In Fig. 10, we present the plot of leptogenesis vs RHN neutrino masses M N 1 and M N 3 and include the M 1 upper bound required by NANOGrav for the first (left panel) and second (right panel) octant. The heaviest RHN neutrino mass M N 3 should not be larger than M 1 and thus not larger than the upper bound of M 1 . Importantly, the RHN masses which give viable leptogenesis are in the region of the PS testable by NANOGrav12.5.
Finally, we briefly discuss the GW signal predicted by the benchmark point. Given the lowest and the second lowest intermediate scales M 1 = 2 × 10 13 GeV and M 2 = 5 × 10 13 GeV and the consistency with gauge unification, Gµ is predicted to be 1.1 × 10 −10 which corresponds to Ω GW (f )h 2 ∼ 4.2 × 10 −10 at the peak around 6.3 × 10 −7 Hz. This is a critical value to be tested in the PTA experiments. In the high frequency band, Ω GW (f )h 2 ∼ 1.6 × 10 −10 is predicted. The future space-based, atomic and undergroundbased GW observatories will give an excellent test of this value.

Discussion and conclusion
Grand Unified Theories offer an attractive ultraviolet completion of the Standard Model and can explain the masses and mixings of the known particles. In our former works, we proposed gravitational wave and proton decay measurements as complementary windows to identify possible breaking chains of GUTs [11] and provided a systematic analysis for all SO(10) chains via the Pati-Salam-type breaking. In addition, we focused on the connection between the proton lifetime and stochastic GW background energy density [12].
In this paper, we perform a detailed analysis of one such SO(10) GUT breaking chain involving a realistic and minimal SO(10) GUT model which breaks to the Standard Model gauge group via three intermediate gauge symmetries. Four Higgs multiplets are required to induce the desired pattern of breaking. An additional two Higgs multiplets are introduced for our model to predict the SM fermion masses and mixing to a high statistical significance. From the constraint of gauge unification and proton decay, and including the relevant particles in our two-loop renormalisation group equations, we determine the scale of spontaneous symmetry breaking from our GUT model to the Standard Model. The final intermediate symmetry breaking induces a network of cosmic string, which results in the emission of GWs, assuming that inflation ends earlier than string formation. Our numerical procedure uses the quark masses and mixing as constraints on our GUT model parameter space (four-dimensional). We performed a grid scan of this space and found that neutrino masses and mixing (five observables) can be accommodated to a high statistical significance. Once the viable regions of the model parameter space were found and the intermediate scales determined, we have a fixed prediction for leptogenesis, gravitational waves and proton decay.
Our scan result shows that the right-handed neutrino spectrum can be strongly restricted as the model has to fit all flavour data, including fermion masses, CKM mixing and PMNS mixing. The heaviest right-handed neutrino mass is predicted to be in a narrow region above 10 13 GeV. Given the chosen breaking chain, this region is still consistent with the upper bound of the lowest intermediated scale restricted by proton decay and GW measurements. It is worth emphasising that recently released data from a series of Pulsar Timing Arrays measurements hints at this region. However, more data is required to confirm if these signals are from a stochastic GW background. In the future, both proton decay measurement experiment Hyper-K and GW observatories can exclude or confirm this region. Furthermore, the latter has the potential to provide powerful constraints to many GUT models.
Throughout, we present a benchmark point, BP1, which satisfies all experimental constraints, including proton decay, quark masses and mixing, lepton masses and mixing and can address the observed matter-antimatter asymmetry in the Universe. By fixing intermediate scales at 2 × 10 13 and 5 × 10 13 GeV, respectively, the GUT breaking scale is solved to be around 5.68 × 10 15 GeV, and the proton lifetime is predicted to be around 5.1 × 10 34 years, which survives the upper bound restricted by Super-K and will be tested in Hyper-K in the future. This benchmark predicts three right-handed neutrino masses (M N 1 , M N 2 , M N 3 ) = (4.23 · 10 11 , 5.32 · 10 11 , 1.66 · 10 13 ) GeV. These heavy neutrinos, together with their CP-violating Yukawa couplings with active neutrinos, can explain the baryon-antibaryon asymmetry of the Universe quantitatively via leptogenesis. The mass of the heaviest right-handed neutrino is lower than the lowest intermediate scale M 1 , thus consistent with the Super-K constraint on proton decay.
In summary, we are entering an exciting era where the culmination of data from proton decay, gravitational wave and neutrinoless double beta decay experiments will allow us to test Grand Unification and its connection with the matter-antimatter asymmetry at unprecedented levels. To illustrate this, we have analysed in full detail a specific example of a realistic and minimal SO(10) model broken to the SM via three intermediate scales, and shown that it will be tested very soon by such experiments. Given the upcoming synergy in experimental data it is of interest to extend our methodology to further GUT models.