Testing Realistic SO (10) SUSY GUTs with Proton Decay and Gravitational Waves

.


I. INTRODUCTION
SO (10) Grand Unified Theories (GUTs) are widely studied, ultraviolet complete frameworks that unify three of the fundamental forces and have unique features [1][2][3]: • Highly correlated fermion masses and mixing as quarks and leptons are arranged in a single representation in the GUT gauge space.• Inclusion of a U (1) B−L gauge symmetry, whose spontaneous breaking gives raise to cosmic strings [4] and right-handed neutrino masses, that can generate light Majorana masses for neutrinos and a baryon asymmetry.• Generically, baryon and lepton number violating operators that induce proton decay and allow to set stringent limits on the GUT breaking scale.These features of SO (10) GUTs allow to constrain the models with present and future data from a variety of complementary approaches.Next-generation large-scale neutrino experiments, including JUNO [5], DUNE [6], and Hyper-Kamiokande [7], are expected to measure the majority of neutrino oscillation parameters at percent level precision.An additional goal of these experiments will be to constrain, or possibly even measure, the pro-ton lifetime at an unprecedented level.Both measurements will probe the theory parameter space of SO (10) GUTs.Moreover, as discussed above, a generic feature of many GUTs is the breaking of a U (1) gauge symmetry which can generate a network of cosmic strings [8].If the string network is not completely diluted by inflation, it may be a source of stochastic gravitational wave background (SGWB) with a broad spectrum of frequency from nanohertz to kilohertz.Due to this broad spectrum, a variety of currently running and upcoming gravitational wave (GW) experiments, including those from Pulsar Timing Arrays (PTAs), space-and ground-based laser interferometers, as well as atomic interferometers, will be able to constrain such GUTs [9][10][11][12][13].Very recently, the NANOGrav collaboration with 15 years of data activity (NANOGrav15) [14][15][16][17][18] reported the evidence for quadrupolar correlations [19], indicating a GW origin of the signal.Similar evidences have been independently reported by EPTA [20][21][22][23][24][25], PPTA [26][27][28], and CPTA [29].Based on the Bayesian analysis of NANOGrav15, it was shown that the Nambu-Goto (NG) strings do not provide a good fit to the signal [18] and sets a stringent bound on the symmetry breaking scale for stable strings.The possibility that metastable cosmic strings [30][31][32] or superstrings [33] may be the source of the observed SGWB remains open.
We have shown that this rapid progress in neutrino and gravitational wave measurements provides complementary probes of SO (10) GUTs [34,35], offering a unique opportunity to indirectly test very high energy scales.In Ref. [36], we carried out a detailed study of SO (10), showing that all Standard Model (SM) fermion masses and mixing parameters and the baryon asymmetry can be matched to their observed values in this model.In that work, we constructed a model with a U (1) B−L symmetry breaking scale around 10 13 GeV, that is not excluded by PTAs but cannot explain the newer indications of SGWB that may originate from metastable strings.
In this Letter, we continue our roadmap on the testability of SO (10) GUTs by extending the analysis to a supersymmetric (SUSY) version.As we will show, SUSY GUT can have marked differences with respect to the non-SUSY models considered already.It offers a natural framework for metastable strings consistent with the NANOGrav signal, as it favours a U (1) B−L breaking very close to the GUT scale.In addition to being currently one of the favoured explanation to the SGWB signal recently observed by the PTA observatories, thanks to its rather broad frequency spectrum the resulting SWGB will be within the sensitivity of future GW experiments at higher frequencies.Moreover, one of the key predictions of SUSY GUTs is kaonic proton decay (p → K ν), in addition to the other non-SUSY decay channels.In the coming years, JUNO will provide the best sensitivity to this channel and place a key constraint on SUSY GUTs.
We perform a comprehensive analysis of this model by determining each scale of symmetry breaking by solving the renormalisation group equations (RGEs) and fitting our model to SM fermion masses and mixing data.We concretely demonstrate that the predictions of our model can be tested by next-generation cosmic microwave background (CMB) observations, neutrinoless double beta decay, oscillation measurements, GW experiments, and searches for proton decay.Moreover, our model can predict the observed baryon asymmetry and accommodate a viable dark matter candidate.

II. MODEL FRAMEWORK
The model we present and confront with flavour, proton decay, and gravitational wave data is a SUSY SO(10) GUT, that is spontaneously broken to the SM as follows: At scale M GUT , the GUT symmetry is broken, and dimension-six operators which mediate proton decay are induced.This GUT symmetry breaking also leads to the production of monopoles which we assume are removed through a period of rapid inflation.The next breaking step occurs at M B−L , where the U (1) B−L gauge symmetry is spontaneously broken.This leads to the production of a network of cosmic strings, that can decay gravitationally, as well as to the generation of right-handed neutrino (RHN) masses.These RHNs can decay to produce a matter-antimatter asymmetry via thermal leptogenesis [37].In the final step, SUSY is broken at a scale M SUSY , defined to be equal to the common masses of the squarks and sleptons, and dimension-five operators which mediate proton decay are induced via wino and higgsino exchange, whose masses, M W , are allowed to be below M SUSY .This is typical of widely studied split SUSY scenarios [38,39].The boldface numbers to the left of the arrows of Eq. ( 1) denote the representations of Higgs superfields of SO (10) required for the symmetry breaking.
In Ref. [35], we analysed all symmetry-breaking patterns of a non-SUSY SO (10) GUT to the SM via the Pati-Salam path by solving the two-loop RGEs and assuming gauge coupling unification at M GUT .We found that the majority of breaking patterns constrain M B−L ≲ 10 13 GeV with no symmetry breaking patterns attaining M B−L ≳ 10 14 GeV.This leads to the formation of cosmic strings with a relatively small string tension which is more difficult to be tested by GW interferometers.One notable motivation to extend from non-SUSY GUTs to SUSY GUTs is that, in the SUSY version, M B−L can naturally reach 10 14 − 10 15 GeV, which leads to the formation of very heavy metastable strings that can accommodate the GW signal detected by the PTAs.To determine the B − L breaking scale in the SUSY version, we follow the same procedure as our previous work Ref. [35], by solving the two-loop RGEs, including threshold effects from gauginos and higgsino masses, which may be a few orders of magnitude lower than M SUSY .The β coefficients at each intermediate scale are listed in Appendix A. From our RGE analysis, we find that the lower the scale of SUSY breaking, the higher the U (1) B−L symmetry breaking scale as shown in Fig. 1.As the wino mass increases, the B − L scale is suppressed while the GUT scale remains roughly the same, enlarging the hierarchy between the two scales, that leads to a stable network of strings.We use the RGE solutions as input for determining the proton decay rate and gravitational wave signal.In the following, we discuss the testability of this model.

A. Fermion masses and mixing angles
At the GUT scale, the Yukawa superpotential is given by: where the asterisk denotes complex conjugation, 16 is the SO(10) matter multiplet and 10, 126 and 120 are the For an observable GW signal, the inflationary scale must lie between the red and blue lines.
Higgs superfields.See Appendix B for more details.In flavour space, the Yukawa matrices are 3×3 with Y 10 and Y 126 complex, symmetric and Y 120 real, antisymmetric.We treat the quark masses and CKM mixing parameters [40,41] as inputs and hence this model has seven free parameters which we vary to predict eight observables in the lepton sector.We perform the procedure using MultiNest [42] to minimise the χ 2 statistical measure as detailed in Appendix C..

B. Leptogenesis and 0νββ decay
For the points in the model parameter space scan which fits the flavour data at high statistical significance, there is a prediction for the Yukawa matrix, Y ν , which couples the RHNs to the leptonic and Higgs doublets and the mass spectrum for the RHNs where the latter is shown in the upper panel of Fig. 2. The heaviest right-handed neutrino mass is constrained at around M N3 ∼ 10 13 GeV and we observed that there is a mild mass hierarchy predicted.We can calculate the baryon asymmetry generated from the decays of these RHNs using ULYSSES [45,46] to solve the density matrix equations.All points of the scan which fit the flavour data allow for viable leptogenesis with a baryon asymmetry of the same order or larger than the observed value [47].We found that the model parameter space highly favours normally ordered neutrino masses, with the lightest neutrino mass in the range 5 ≲ m ν1 (meV) ≲ 15.In the lower panel of Fig. 2, we show the predictions of our scan for the effective Majorana mass (m ββ ) as a function of the sum of neutrino masses ( i m νi ).Both observables are testable by the next generation of ν0ββ [48][49][50][51][52] and CMB experiments [43,44,53] with their sensitivities shown in grey and blue, respectively.We note that all the points shown fit the flavour data well (χ 2 ≲ 10), FIG. 2. All coloured points show regions of the parameter space that fit the neutrino sector with χ 2 < 10.The yellow star shows a benchmark point in the scan with ηB ∼ 6.2 × 10 −10 .All the points in the scan predict −10 ≲ log 10 (ηB) ≲ −8.In the upper plot, we see that the heaviest right-handed neutrino mass is predicted to be MN 3 ∼ 10 13 GeV, with the lightest right-handed neutrino mass MN 1 an order of magnitude less, and the mass hierarchy is mild, with the lightest right-handed neutrino mass MN 1 an order of magnitude less.In the lower plot, we show the effective neutrino mass (m ββ ) predicted from our model as a function of the sum of neutrino masses.The grey regions shows the reach of the next generation ν0ββ experiments and the blue region shows the sensitivity for EUCLID [43,44].
and our benchmark point (indicated by the associated yellow star in Fig. 2 see Appendix C 2 for the Yukawa matrices) is consistent with the flavour data (χ 2 ≲ 3) and provides a prediction of the baryon-to-photon ratio η B ∼ 6.2 × 10 −10 .

C. Proton Decay
The proton decay bound on the GUT scale can be translated to the restriction on the SUSY breaking scale and wino mass from the solutions of the RGEs (see Fig. 1) and the assumption of gauge coupling unification.Increasing m W suppresses M GUT , and therefore increases the proton decay rate to a level which can be tested in Hyper-K, as later summarised in Fig. 5.Such low mass sparticles may be within reach at the FCC [54], offering another avenue to test the model.As the GUT we study is supersymmetric, additional contributions to proton decay from the colour-triplet Higgs superfields can mediate the baryon-antilepton transition.The consequence is that some decay channels, such as p → K + ν, are enhanced.Although this channel is less well experimentally constrained, τ K ν ≳ 5.9 × 10 33 years in Super-K [55], the SUSY GUT provides τ K ν ∝ M 2 GUT M 2 SUSY , and thus this channel can lead to stronger constraints on the SUSY GUT.Since we consider split SUSY spectrum, there is an additional enhancement to the partial lifetime by a factor M 2 SUSY /m 2 W in the case m W ≪ M SUSY .Moreover, as this channel originates from the Yukawa superpotential in Eq. ( 2), the partial lifetime is determined by the Yukawa coupling matrices, which are almost entirely fixed, up to overall order-one factors, by our fit to the fermion flavour data.Further details of the proton lifetime calculations are provided in Appendix D. We will discuss how the pionic and kaonic decay channels can constrain the GUT model parameter space and the non-trivial interplay with the GW predictions in the discussion section.

D. Dark Matter
If R-parity is conserved after SUSY breaking, the lightest SUSY particle (LSP) would be stable and thus can be a dark matter candidate if it is a neutralino [56,57].Due to the observed relic abundance, an upper limit on the mass of LSP can be obtained to avoid over-abundance.Generally, the lightest neutralino can be the mixture of wino, bino and higgsino.The maximal dark matter mass can be as high as 10 TeV with resonant heavy Higgs annihilation or enhancement in annihilation rate through next-to-LSP [58].On the other hand, a pure bino LSP or bino-dominated LSP commonly leads to an over-abundance unless (co)annihilation processes further reduce the relic density.The upper limit of pure wino dark matter mass is around 3 TeV [59], and for pure higgsino the limit is around 1 TeV [38] which can account for the correct dark matter relic density.

E. Gravitational Waves
When a U (1) symmetry is broken at a certain scale M B−L , a cosmic strings network is generated.The long strings in the network can intersect to form loops that oscillate and emit energy via gravitational radiation.Such radiation is not coherent and can be seen as a stochastic background of gravitational waves.Importantly, this background can, in principle, be observed by currently running and future GW experiments.
We begin with the Nambu-Goto string approximation, where the string is infinitely thin with no couplings to particles [60], 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 approximately given by [61]: In the case that M B−L is not far away from M GUT , deviation of α 2R (M B−L ) and α 1X (M B−L ) from α GUT (M GUT ) due to RG running is small.So we can approximate Moreover, hence we can relate the string tension parameter to the intermediate scale As the model we consider has U (1) breaking energy scale close to the GUT symmetry breaking (which generates monopoles), the cosmic strings network can decay via monopole-antimonopole nucleation as studied in Refs.[11,62,63].After the stage of string formation, the cosmic strings network and the consequent loops start to decay producing monopoles-antimonopoles pairs.The exponential suppression of the loop number density is characterised by the decay width per unit length Γ d , for strings and the consequent loops, where m = M V α f m the monopole mass, and f m is an undetermined factor depending on the model's detail usually assumed of order one [64].In practice, it translates into a cutoff in the low-frequency spectrum of the gravitational waves background which avoids the region tested by Pulsar Time Arrays and provides the appropriate tilt in the spectrum observed by the PTA experiments.These signals can be tested by high-energy gravitational waves experiments such as the Einstein Telescope and, at high energies, by LIGO-Virgo-KAGRA.
Although the stable cosmic strings are disfavoured, the source of the signal detected by the PTAs [14,18,20,28,29] may be a network of metastable cosmic strings, which occurs when the string network decays to monopole-antimonopole pairs [65][66][67][68][69].The hierarchy between the GUT scale and the string formation scale can be parametrised by where an order-one coefficient is ignored.The smaller κ, the closer the GUT and string scales and the more efficient the annihilation of the string network.As we have shown, this SUSY SO(10) GUT prefers a small hierarchy between M GUT and M B−L , which naturally leads to a prediction of metastable cosmic strings.Metastable strings have been suggested as a possible explanation of the PTA observations which favour √ κ ≈ 8.We can use the PTA observations as one of the strongest constraints for a stable network of cosmic strings requiring that Gµ < 2 × 10 −10 [16] when √ κ ≫ 9, which corresponds to M B−L ≲ 6 × 10 13 GeV due to Eq. (5).
By assuming gauge unification, we can directly connect the intermediate scale of the GUT symmetry breaking with the SUSY breaking scale, which we assume to be of the same order of magnitude as the sfermions masses.Therefore, we can constrain the SUSY breaking scale using gravitational waves.It is worth noticing that the link between the SUSY breaking scale and M B−L depends on the mass of the gauginos, which, as we can see in Fig. 1, heavily affects the running of the gauge coupling.

F. Results and Discussion
Here, we consider the various GUT observables and their interplay and focus on three benchmark points, representative of the three key behaviours of the model: BP1 This has a high SUSY-breaking scale, (M SUSY ∼ 10 9 GeV), as well as wino masses, leading to a B−L breaking scale which is lower than BP2 and BP3.The model exhibits characteristics very similar to the non-SUSY SO(10) GUT.
BP2 The SUSY and wino mass scales are still quite high but lower than BP1, and a prediction for proton decay via the SUSY channel p → Kν can be achieved close to current bounds.Therefore, this case could be differentiated from non-SUSY SO (10).
BP3 This is the most characteristic case for this model.Thanks to the low wino mass and relatively low SUSY mass scale, the B − L breaking scale is very close to M GUT leading to the possible generation of metastable strings and a viable explanation of NANOGrav15.A sizable prediction for the SUSY proton decay channel also emerges.
In Fig. 3, we study in detail the GW predictions in the three cases, confronting them with NANOGrav15 and the LIGO-Virgo-KAGRA (LVK) bounds [70] on the spectrum at the high-frequency band.For BP1, the green curve shows the GW spectrum from stable strings, consistent with all constraints, including the upper bound set by PTA.We note that this benchmark does not provide an explanation of the PTA observation.For BP2, the red curve shows a GW spectrum from stable strings with a larger Gµ value, that is inconsistent with the PTA observations, showing the importance of GW observations in constraining the model (this assumes a sufficiently high inflationary scale).Finally, for BP3, the blue curve shows the GW spectrum from metastable strings diluted by the inflation in the high-frequency band, which fits NANOGrav15 very well and is testable with future GW observations.It is important to note that since the generation of a metastable string network requires M B−L ∼ M GUT , without any other assumptions, this model tends to favour a signal that would be excluded by the current observations from LVK [70].Therefore, to provide an explanation of the PTA observations, we require the cosmic string network to be partially diluted by inflation [31,71], thereby suppressing the signal in the higher frequency regime.Since M B−L , M GUT , and inflation are at approximately the same scale in such a scenario, it is plausible to have the cosmic string network generated towards the end of inflation and, therefore, are slightly inflated away with a typical 1/f suppression in the region in the high-frequency regime, as shown in Fig. 3, for a typical time t F ∼ 10 −10 s.Comparison between SGWB produced by a metastable network of cosmic strings with Gµ from 10 −7 to 10 −5 and three different values of κ, and the 1σ, 2σ, 3σ regions of the signal detected by NANOGRAV [14], EPTA [20], and PPTA [28].
Assuming a power-law spectrum of the characteristic GW strain, the GW energy density spectrum can be parameterised by two parameters: the amplitude parameter, A, and power parameter, γ.Indeed, the characteristic strain of the signal detected by Pulsar Time Arrays can be parameterised as and the GW energy density spectrum , in the nHz region, is given by [72] where We perform the fit in the interval 2 − 59 nHz following the procedure of Refs.[36,72], and the results are shown in Fig. 4, where the value of A and γ favoured by observation is compared with the prediction of the metastable cosmic string network.The reference frequency of the PTA results is chosen to be 1yr −1 .We found that the signal produced by a metastable network of cosmic strings is compatible with the EPTA and PPTA in the 1σ range and with NANOGrav in the 2σ region.From Fig. 4 we note that the predicted signal is very sensitive to the value of κ.Further, the values of Gµ consistent with the observations from the Pulsar Time Array observatories get higher as κ becomes smaller.For the values of κ we have shown, κ = 7.85, 7.9, 7.95, the values of Gµ for which the predicted signal is within the 2σ region, are between 10 −5 and 10 −7 , depending on κ, and they are consistent to the scenario in which M B−L ≃ M GUT .
FIG. 5. Constraints and sensitivities on proton decay in the MSUSY −MWino plane, (where MSUSY is defined as the squark and slepton mass scale).The orange lines show the constraint of Super-K (solid) and sensitivity of Hyper-K (dashed) on the p → π 0 e + channel proton decay [7].Super-K excludes the parameter space to the left of the purple solid line due to strong p → K + ν channel decay.To the left of the purple dashed line, the model can predict a proton lifetime that is within the sensitivity of JUNO [5].The narrow cyan region shows where the NANOGrav result can be explained at 95% confidence [16].
The benchmark points can be located also in Fig. 5, which presents the testability of this model in the M SUSY − m W plane.The hatched region on the top-left indicates (m W > M SUSY ) is disfavoured in most SUSY scenarios.Gauge unification excludes the grid region on the bottom-left and the purple and orange solid lines are excluded by the current bounds on τ K + ν and τ π 0 e + from SUSY and non-SUSY contributions, respectively.The orange dashed line indicates the sensitivity of Hyper-K to π 0 e + channel decay, while the purple dashed line shows the potential target region of JUNO on K + ν channel decay.A summary of the predictions of the model for different benchmark points is given in Table I.BP1 has a phenomenology similar to non-SUSY SO(10) models [36].BP2 predicts proton decay rates that Hyper-K and JUNO can probe, although it has a string tension too high to be compatible with PTA observations.The region of the parameter space that predicts metastable strings that can explain the PTA observations, such as BP3, is in the region Gµ ≳ 10 −6 with √ κ ∼ 8 and can predict proton decay rates in the kaon channel at reach at JUNO, but interestingly, not at Hyper-K.Moreover, we notice that this region requires wino masses that are at reach at present or future colliders Refs.[54,[73][74][75], with lower masses the lower the proton decay rate.Interestingly, most of the parameter space that predicts proton decay rates too small to be observed by Hyper-K or JUNO is inconsistent with PTA observations.In the non-supersymmetric SO(10) unification, it is also possible to obtain gauge unification, predict fermion masses and mixing and the dark matter [76].However, to achieve successful gauge unification and leptogenesis in the non-SUSY SO(10) framework with a simple mass spectrum, the intermediate scale, that generates righthanded neutrinos masses and cosmic strings, is at least an order of magnitude smaller than the GUT scale [35,36].In that case, intermediate U (1) symmetry breaking can only lead to stable cosmic strings and thus the resulting gravitational wave cannot provide an explanation of the spectrum observed by PTAs.On the contrary, the intermediate U (1) symmetry-breaking scale can be naturally close to the GUT scale in SUSY SO(10) GUTs, which leads to the prediction of a metastable cosmic string network that can explain the PTA observations.Moreover, the proton decay predicted by SUSY GUTs is in general faster than that predicted by non-SUSY GUTs due to the additional kaon channel.As a result, the SUSY GUTs is more testable than the non-SUSY ones.As shown in Fig. 5, the current bound on kaon channel proton decay rules out the possibility of low scale SUSY, which explains the absence of superpartners at LHC. Apart from the phenomenological perspectives, the hierarchy problem also motivates SUSY theoretically, further making SUSY SO (10) an attractive candidate for unification.

III. SUMMARY
We have presented a SUSY SO (10) GUT which can successfully predict fermion masses and mixing angles, leptogenesis, dark matter and can be tested via proton decay and GW signatures at next-generation experiments.The B − L breaking scale correlates with the SUSY breaking scale (the squark and slepton mass scale) via the gauge unification.The natural proximity of the GUT breaking scale and the B − L breaking scale leads to metastable cosmic strings decaying to monopole-antimonopole pairs, and we find that the GW signal from metastable strings can be consistent with the NANOGrav 15-year data.Considering a split-SUSY scenario we found that proton decay measurements and PTA observations cover complementary regions of the parameter space.An eventual observation of proton decay from both the pion and kaon channels is not consistent with the current PTA observations.Exploiting this complementarity, we can, therefore, test the majority of the parameter space.Regarding the interpretation of the observed GW signal as generated by a metastable cosmic string network, we found that it is consistent with our model if the signal is partially inflated away, and that it is possible to fully test this possibility with the nextgeneration GW observatories and JUNO and/or collider searches.
A necessary condition for a realistic GUT model is that all gauge couplings unified to a single gauge coupling at a certain scale given by g GUT and M GUT , respectively, up to matching conditions.We begin with the two-loop RGEs.Given an interval of energy scale Q ∈ (Q 0 , Q 1 ), where the gauge symmetry is G = H 1 ×H 2 ×• • •×H n and no particles decouple in this period, the coupling α i = g 2 i /(4π) for each gauge symmetry ) is described by the following differential equation The β function, up to the two-loop level, is given by where 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 R.G. running equations as it gives a subdominant contribution.
If the conditions b j α j (Q 0 ) log(Q/Q 0 ) < 1 is satisfied, an analytical solution for these equations can be obtained [77]: In non-SUSY and SUSY, the coefficients b i and b ij are respectively given by non-SUSY: where ϕ, ψ and Φ represent any complex scalar, chiral fermion and chiral superfield evolving in the scale between Q 0 and Q 1 , respectively, C 2 (H i ) is the quadratic Casimir invariant of the adjoint presentation in the group is quadratic Casimir invariant of the representation of the field ϕ (superfield Φ) in the group is the Dynkin index of the field ϕ (superfield Φ) in group H i , and the summation goes over all fields (superfields).In our model, values of coefficients b i and b ij in each interval of the energy scale from M GUT down to M Z are shown in Table II.At the intermediate scale, where a larger symmetry is broken to its sub-symmetry, gauge couplings between them 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 in the MS scheme is given by [78] H i+1 → H i : C 2 (H i ) .

(A6)
Above the SUSY scale, the couplings must run in the D.R. scheme to preserve the supersymmetry [79].The relation of couplings in the MS scheme and DR scheme is described by Thus the one-loop matching condition in the DR scheme is simply For G int → G MSSM , we encounter the breaking, SU (2) R × U (1) X → U (1) Y , where U (1) X is identical to U (1) B−L with the charge normalised as X = 3/8(B − L).The matching condition in DR scheme is given by Applying the matching conditions, all gauge couplings of the subgroups unify into a single gauge coupling of SO( 10) at the GUT scale, M GUT .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 GUT while M B−L and M SUSY are associated to the breaking of G int and G MSSM , respectively.
SO (10) broken at   II.Coefficients bi and bij of gauge coupling β functions appearing in the specified breaking chain.In this table, we identify the scale of gauge symmetry as the corresponding heavy gauge boson mass from the EFT point of view.The SUSY breaking scale MMSSM is regarded as the unified mass of all sfermions.In the split SUSY, gauge bosons and Higgs superpartners are allowed to be different from and usually a few orders of magnitude lower than the SUSY-breaking scale.We also consider the possibility of having a further gap between the wino and gluino masses and the bino and higgsino masses.We include their threshold effect in the RG running.
. respectively, as follows: where and V ij and U ij denotes the mixing between the mass and interaction basis of Higgs doublets in the up and down SO(10) 10 126 120 IV.Decomposition of Higgses responsible for the fermion mass generation.126 is the same Higgs as shown in Table III and it is responsible for both the breaking Gint → GSM and right-handed neutrino mass generation.(1, 1, 0)S is the same singlet given in Table III.SO (10) 16 sectors, respectively, before the SUSY breaking [82,83], and v S is the VEV of singlet component of 126 that gives mass to the RHNs.The light neutrino mass matrix, M ν , is obtained by where Once the scale v S is determined, U 12 can be solved as 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 [85,86], 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.As a result, h, f and h ′ , as well as all parameters on the righthand side of Eq. (C1), are real.Since h ′ is antisymmetric, we arrive at Hermitian Dirac Yukawa coupling matrices Y u , Y d , Y ν and Y e .The Higgs mixing elements V 13 , V 14 and U 13 , U 14 are also purely imaginary with the relations This texture has been widely applied in the literature, e.g., Refs.[82,87,88].The resulting fermion mass matrices conserve parity symmetry L ↔ R [86] and following from the assumption that there is no C.P. 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.

Procedure to fit the quark and lepton flavour data
For simplicity, we assume that r 3 = 0, which implies that the imaginary part of Y u vanishes.As a result, the relation between V 13 and V 14 in (C5) is no longer valid.Instead, there is a simpler relation that reads 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 the real-orthogonal transformation cannot determine.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 ia1 , e ia2 , 1}.The matrices h, f and h ′ are then expressed in terms of where Y ν , Y e are 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 mass-squared 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.
By fitting the fermion mass and mixing, the matrices h, f , h ′ and parameters a 1 , a 2 , r 1 , r 2 , c e , c ν , and m 0 can be fully determined.To perform the parameter scan, and find viable regions of the model parameter space that postdict the quark and predict the leptonic data, we run all the SM Yukawa couplings to the GUT scale (using two-loop RGEs, appropriate matching between scales and threshold corrections) using REAP [89,90] and SARAH [91].We then scan the free parameters of the GUT model as described above and assess how well they fit the leptonic data using the statistical measure below: where P m ∈ {a 1 , a 2 , r 1 , r 2 , c e , c ν , m 0 , η q } and O n ∈ {m e , m µ , m τ , θ 12 , θ 13 , θ 23 , ∆m 2 21 , ∆m 2 31 }.Then Yukawa couplings for SO(10) 16 multiplets can be expressed as this will be relevant for the subsequent discussion on proton decay, tightly linked with the scan of fermion masses and mixing since it is mediated by the Higgs colour triplet and the operators computed from the same superpotential.Apart from the matrices and parameters that can be fixed by fermion mass and mixing, there are three parameters: two Higgs mixing elements V 11 and V 13 and tan β.Apart from the equations above, there are also while U 12 can be solved using (C4).Each element in the mixing matrices has to satisfy the unitarity.We show a subset of our predictions from the scan in Fig. 6.In the left (right) panel, we show the predictions for the muon (δ phase) Yukawa versus the electron Yukawa (θ 23 ).We note that the predictions are consistent with the experimental values at the 3σ level.While all values of δ can be accommodated, the model prefers the atmospheric mixing angle to be in the lower octant.Moreover, the model strongly prefers normally ordered light neutrino masses.

A benchmark study
We considered a benchmark point of our scan (yellow star), achieving successful leptogenesis, giving η B ∼ 6.16 × 10 −10 .Inputs and predictions of fermion Yukawas and mixing parameters are shown in Table VI.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 = 8.22.

Appendix D: Proton decay
Unified theories may contain gauge or scalar bosons that mediate baryon (B) number and lepton number (L) violating processes and hence can cause the decay of nucleons.The dominant mediator is the heavy colour gauge boson which leads to proton decay where the most constrained channel is p → π 0 e + with the lifetime τ π 0 e + ≳ 2.4 × 10 34 years set by Super-K [92].Due to the direct correlation τ π 0 e + ∝ M 4  GUT , this channel usually provides the most constraining limit on the GUT scale, M GUT ≳ 3 × 10 15 GeV, almost regardless of the breaking chains of SO(10) [35].
The following section discusses proton decay induced by GUT and SUSY symmetry breaking, respectively.

a. Pion channel
The computation of the proton lifetime in this channel is identical in both the SUSY and the non-SUSY case, and it is carried by the heavy SO(10) gauge bosons.The relevant dimension-six operators are written as where colour and flavour indices have been suppressed and Λ ≃ √ 2M GUT /g GUT denotes the UV completion scale.The proton lifetime is [35]: (D2) The enhancement factors denoted as A L , A S.L. , and A S.R. , correspond to the influence of long and short-range effects on proton decay.The relevant hadronic matrix element for this particular decay mode is ⟨π 0 |(ud) L,R u L |p⟩, which has been determined through a QCD lattice simulation [93].The long-range effect incorporates the renormalisation enhancement from the proton decay process (at scales ∼ 1 GeV) to the electroweak scale (defined as the mass of the Z boson at the scale of M Z ).This enhancement factor, calculated at the two-loop level, for is A L = 1.247 in the studies [94].In the SUSY case, we used A L = 0.4 [95,96].The short-range factors are determined through the renormalisation group equations, which span from the scale M Z to M GUT : where γ i are the anomalous dimensions and b i the oneloop β coefficients.

Kaon channel
For the computation of the proton decay in the kaon channel, we followed the treatment of Ref. [97,98].The baryon number violating interaction is mediated by the Higgs colour triplet with mass M T ≃ M GUT .The terms in the superpotential, which we consider, are the same as the Yukawa sector in Eq. ( 2).Let us consider the effective superpotential from which we can infer all the 5dimensional operators contributing to the kaon channel, where the matter should be understood ) αγ (Y 120 ) βδ , (D4) where x i and y i are the matrices' elements that diagonalise the colour triplet mass matrix and can be treated as free parameters that we vary randomly in the interval [0, 1].In Eq. (D4) all the Yukawa couplings are rotated to be in the mass basis.The coefficients C L αβγδ and C R αβγδ induces uncertainties for the partial lifetime.Fitting of fermion masses and mixing helps to obtain h, g and h ′ , overall factors between them and Y 10 , Y 126 and Y 120 , in particular, V 11 between h and Y 10 , are undermined.By requiring the perturbativity of the theory, we scan and obtain the maximal and minimal contributions of these coefficients to the proton decay.In Fig. ( 4), the exclusion region for mass scales set by Super-K and the future sensitivity of JUNO is obtained by considering the maximal and minimal contribution of coefficients, respectively.In order to be able to predict proton decay, the squarks or the sleptons in the 5-dimensional operators of Eq. (D4) need to be dressed with gaugino or higgsino vertices.We note that from the model's symmetry, namely that a bino or a gluino dressing gives zero contribution due to the Fierz identity, only the contributions from the wino and the higgsino dressing significantly promote decay in this channel.In contrast, the dressing from gluino and binos is negligible.The sub-operators we consider are [97,98] where we followed the notation of [98] and the superscripts label the Feynman diagram, which contributes to the kaon channel, and the subscript label the gaugino, which dresses the Feynman diagram.Considering the above sub-operators, we note that they are proportional to the mixing matrices and the Yukawa couplings at low energies.Therefore the proton lifetime is tightly linked to the predictions of the scan, and therefore is possible to test the parameter space of the Yukawa sector of the theory with proton decay experiments.In particular, the free parameters on which the overall lifetime depends are r 1 , r 2 , a 1 , a 2 , c e , c ν which are fixed from the scan, x i , y i , and the Higgs mixing parameters U 12 , V 11 and U 13 .We found that by randomly varying x i and y i , we only achieve a few percent difference in the lifetime while the influence of the Higgs mixing parameters is much more important.The parameter U 12 can be fixed since it depends on the known parameters m 0 , r 1 , and M B−L .The other two parameters are constrained by relations of Eq. (C13), and we are considering the highest allowed value.We can compute the proton decay width considering two different operators proportional to the sum of the coefficients C W and C h , which are called respectively O W and Oh.These two operators are expressed as: where M T ∼ M GUT is the mass of the Higgs triplet and I(a, b) is the loop contribute; we have I(a, b) ≃ a/b 2 when a ≪ b.Finally, the proton decay width can be expressed in terms of these two operators as From Eq. D6 and Eq.D7, we can also see that the lifetime is highly dependent on the value of the ratio between the gaugino mass and the SUSY breaking scale and from the SUSY breaking scale itself.

FIG. 1 .
FIG. 1.The GUT scale, MGUT (red), and U (1)B−L breaking scale, MB−L (blue), as a function of the SUSY breaking scale MSUSY, defined as the common squark and slepton mass scale.The solid and dashed lines show the effect on the RGEs solutions for various wino masses, assumed to be below MSUSY.For an observable GW signal, the inflationary scale must lie between the red and blue lines.

5 Ω GW h 2 BP1:FIG. 3 .
FIG.3.GW spectra of the benchmark points.The model naturally accommodates a signal generated by a partially inflated metastable network (BP3, blue), which supports NANOGrav15.BP1 (BP2) shown in green (red) predicts stable strings which is consistent (inconsistent) with the PTA observations.
and higgsino decoupling happening at Q = M B

c 12 c 13 s 12 c 13 s
C7) where again η d,s,b = ±1 represent the signs of eigenvalues, and V CKM is the CKM matrix parametrised in the following form 13 e −iδq −s 12 c 23 − c 12 s 13 s 23 e iδq c 12 c 23 − s 12 s 13 s 23 e iδq c 13 s 23 s 12 s 23 − c 12 s 13 c 23 e iδq −c 12 s 23 − s 12 s 13 c 23 e iδq c 13 c 23

FIG. 6 .
FIG. 6.All coloured points show regions of the parameter space which fit the lepton sector with χ 2 < 10.The yellow star shows a benchmark point in the scan with ηB ∼ 6.16 × 10 −10 .All the points in the scan predict −11 ≲ log 10 (ηB) ≲ −8.

TABLE I
. Complementary predictions of BP in the nextgeneration proton decay measurements and NANOGrav15.

TABLE III .
Decomposition of the Higgses which induce spontaneous symmetry breaking at each step of the breaking chain.

TABLE V .
Decomposition of the matter multiplet 16 in each step of the breaking chain.

TABLE VI .
10 11 GeV 1.53 • 10 12 GeV 4.67 • 10 13 GeV Inputs and predictions of neutrino masses and the benchmark point mixing parameters fully satisfy all experimental data.Charged fermion masses and CKM mixing are all fixed at experimental best-fit values.Neutrino masses with normal ordering are predicted.