Probing Minimal Grand Unification through Gravitational Waves, Proton Decay,
and Fermion Masses

Shaikh Saad222E-mail<EMAIL_ADDRESS>
Department of Physics, University of Basel, Klingelbergstrasse 82, CH-4056
Basel, Switzerland

###### Abstract

Motivated by the direct discovery of gravitational waves (GWs) from black
holes and neutron stars, there is a growing interest in investigating GWs from
other sources. Among them, GWs from cosmic strings are particularly
fascinating since they naturally appear in a large class of grand unified
theories (GUTs). Remarkably, a series of pulsar-timing arrays (PTAs) might
have already observed GWs in the nHz regime, hinting towards forming a cosmic
string network in the early universe, which could originate from phase
transition associated with the seesaw scale emerging from GUT. In this work,
we show that if these observations from PTAs are confirmed, GWs from cosmic
strings, when combined with fermion masses, gauge coupling unification, and
proton decay constraints, the parameter space of the minimal SO(10) GUT
becomes exceedingly restrictive. The proposed minimal model is highly
predictive and will be fully tested in a number of upcoming gravitational wave
observatories.

###### Contents

  1. 1 Introduction
  2. 2 The Yukawa sector
  3. 3 Fit to the fermion masses and mixings
  4. 4 Symmetry breaking chain and topological defects
  5. 5 Gravitational waves from cosmic strings
    1. 5.1 Observational signals from strings
    2. 5.2 Fit to fermion mass with fixed intermediate scale
  6. 6 Proton decay and Gauge coupling unification
    1. 6.1 Proton decay
    2. 6.2 Unification
  7. 7 Conclusions
  8. 1 Best fit parameters without restricting the seesaw scale
  9. 2 Best fit parameters for fixed seesaw scale
  10. 3 Decomposition of representations

##  1 Introduction

The basic idea of grand unified theories (GUTs) [1, 2, 3, 4, 5, 6, 7] is
embedding the Standard Model (SM) gauge group in a simple non-Abelian group.
Consequently, the three gauge couplings associated with the SM gauge group
$SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$ that have different values at the
low energies are unified into a single gauge coupling at a very high energy
scale, namely, the GUT scale, $M_{\textrm{GUT}}$. In this setup, all the SM
fermions of each generation can also be unified within a single representation
under the chosen GUT group.

The minimal simple group that can serve the aforementioned purposes is the
$SO(10)$ group [5, 6]. Interestingly, the $16$-dimensional spinorial
representation of $SO(10)$ group, in addition to all the SM fermions, contains
a right-handed neutrino $\nu_{R}$, which via seesaw mechanism can provide the
SM neutrinos a tiny non-zero masses. To evade the current proton decay
constraints, the GUT breaking scale has to be quite large, roughly,
$M_{\textrm{GUT}}\gtrsim 10^{16}$ GeV, whereas explaining the light neutrino
oscillation data requires the right-handed neutrinos to have masses much
smaller than the GUT scale that fixes the intermediate breaking scale to be of
order $v_{R}\sim 10^{12}-10^{14}$ GeV.

Remarkably, this intermediate stage symmetry breaking that generates the
right-handed neutrino masses by breaking the $B-L$ ($U(1)_{B-L}$) by two units
leads to the formation of a cosmic string network in the early Universe [8, 9,
10]. Among all the topological defects formed during the subsequent breaking
of the GUT symmetry to the SM group, cosmic strings are particularly
interesting because, unlike monopoles and domain walls that are ruled out, the
existence of cosmic strings has compelling observational consequences. These
cosmic strings intercommunicate and eventually emit gravitational radiation
forming a stochastic gravitational wave (GW) background that ongoing
experiments and future GW interferometers can test.

Inspired by the direct discovery of GWs from black holes and neutron stars by
the LIGO and Virgo Collaborations [11, 12, 13, 14, 15, 16, 17, 18] at
frequencies $f\gtrsim 10$ Hz, there is widespread interest in experiments
exploring GW signals in other frequency ranges. Among others, Pulsar-timing
arrays (PTAs) are sensitive to GW signals with frequencies in the nHz range
and ideal for detecting GWs originating from cosmic strings. Intriguingly,
pulsar-timing arrays: NANOGrav [19], PPTA [20], EPTA [21], and IPTA [22] might
have already seen evidence of a primordial gravitational wave background that
could originate from the breaking of the seesaw scale emerging from $SO(10)$
GUT. Moreover, determining the intermediate symmetry breaking scale from
neutrino oscillation data as well as from gravitational wave signal strongly
constrains the GUT scale, which subsequently has a profound effect on
determining the proton lifetime.

This work considers a minimal grand unified theory based on $SO(10)$ gauge
symmetry. The presence of no other global or discrete symmetries is assumed.
Within this framework, the most minimal Yukawa sector consists of a real
$10_{H}$, a real $120_{H}$, and a complex $126_{H}$ dimensional
representations [23]. The Yukawa sector of this theory has only a limited
number of parameters to reproduce the charged fermion and neutrino masses and
mixing angles. Owing to the minimality criterion, we restrict ourselves to
small representations and employ a $45_{H}$ and a $54_{H}$ dimensional Higgs
representations to break the GUT symmetry. By acquiring vacuum expectation
values (VEVs), the SM singlet components of $45_{H}+54_{H}$ fields break the
GUT symmetry to an intermediate symmetry that contains $U(1)_{B-L}$ symmetry.
Subsequent breaking by $126_{H}$ multiplet gives super-heavy masses to the
right-handed neutrinos and is responsible for the formation of the cosmic
string network. Our detailed study reveals that only a small range for the
seesaw scale is permitted when constraints from GWs detectors, SM fermion mass
spectrum, gauge coupling unification, and proton decay are imposed, making the
proposed model highly predictive and fully testable.

This paper is organized in the following way: in Secs. 2 and 3 we discuss the
details of the Yukawa sector and perform numerical fits. In Sec. 4, symmetry
breaking chain and the associated topological defects are elaborated.
Gravitational waves from cosmic string networks is discussed in Sec. 5.
Furthermore, unification and proton decay constraints are examined in Sec. 6.
Finally, we conclude in Sec. 7.

##  2 The Yukawa sector

In a renormalizable $SO(10)$ theory with no new fermions beyond the three
families of chiral 16s, the Higgs multiplets that can contribute to the
fermion masses can be identified from the following fermion bilinears:

$16\times 16=10_{s}+120_{a}+126_{s}.$ (1)

The 10 and the 120 are real representations in $SO(10)$, whereas the 126 is
complex. In the above equation, the subscripts $s$ and $a$ stand for symmetric
and antisymmetric components in the family space. With this, the most general
renormalizable Yukawa sector takes the following form:

$\mathcal{L}_{yuk}=16_{F}(Y_{10}10_{H}+Y_{120}120_{H}+Y_{126}\overline{126}_{H})16_{F}.$
(2)

It turns out that the Yukawa sector of $SO(10)$ GUTs are very predictive,
which have been extensively analyzed in the literature [24, 25, 26, 27, 28,
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 23, 40, 41, 42, 43, 44, 45].

The up-quark, down-quark, charged leptons, Dirac neutrino, and Majorana
neutrino mass matrices derived from Eq. (2) can now be written down [23]:

$\displaystyle M_{U}=\underbrace{v_{10}Y_{10}}_{\equiv
D}+\underbrace{v^{u}_{126}Y_{126}}_{\equiv
S}+\underbrace{(v^{(1)}_{120}+v^{(15)}_{120})Y_{120}}_{\equiv A},$ (3)
$\displaystyle
M_{D}=v^{\ast}_{10}Y_{10}+v^{d}_{126}Y_{126}+(v^{(1)\ast}_{120}+v^{(15)\ast}_{120})Y_{120},$
(4) $\displaystyle
M_{E}=v^{\ast}_{10}Y_{10}-3v^{d}_{126}Y_{126}+(v^{(1)\ast}_{120}-3v^{(15)\ast}_{120})Y_{120},$
(5) $\displaystyle
M_{\nu_{D}}=v_{10}Y_{10}-3v^{u}_{126}Y_{126}+(v^{(1)}_{120}-3v^{(15)}_{120})Y_{120},$
(6) $\displaystyle M_{\nu_{R}}=v_{R}Y_{126}.$ (7)

One crucial feature of the above set of mass matrices is that the same factor
of $|v_{10}|$ appears both in the up (up-type quark and Dirac neutrino) and
the down (down-type quark and charged lepton) sectors. Since $10_{H}$ is a
real field in this model, the Higgs doublet residing inside $10_{H}$ is self-
conjugate. As a result, the bi-doublet $\Phi=(1,2,2)\subset 10_{H}$ under the
Pati-Salam group satisfies the following condition,

$\displaystyle\Phi^{*}=\tau_{2}\Phi\tau_{2}$ (8)
$\displaystyle\mathrm{i.e.,}\;$
$\displaystyle\Phi=\begin{pmatrix}\phi_{0}&\phi^{+}\\\
\phi^{-}&\phi_{0}^{*}\end{pmatrix},$ (9)

leading to identical vacuum expectation values both in the up and the down
sectors, $v_{u}=v^{*}_{d}=v_{10}$.

By defining

$\displaystyle r_{1}$ $\displaystyle=\frac{v^{d}_{126}\textbullet}{\textbullet
v^{u}_{126}},\;\;r_{2}=\frac{\textbullet
v^{(1)\ast}_{120}-3v^{(15)\ast}_{120}}{\textbullet
v^{(1)}_{120}+v^{(15)}_{120}},\;\;e^{i\phi}=\frac{\textbullet
v^{(1)\ast}_{120}+v^{(15)\ast}_{120}}{\textbullet
v^{(1)}_{120}+v^{(15)}_{120}},\;\;c_{R}=\frac{\textbullet v_{R}}{\textbullet
v^{u}_{126}},$ (10)

we re-write the above mass matrices as,

$\displaystyle M_{U}=D+S+A\equiv vY_{U},$ (11) $\displaystyle
M_{D}=D+r_{1}S+e^{i\phi}A\equiv vY_{D},$ (12) $\displaystyle
M_{E}=D-3r_{1}S+r_{2}A\equiv vY_{E},$ (13) $\displaystyle
M_{\nu_{D}}=D-3S+r^{\ast}_{2}e^{i\phi}A\equiv vY_{\nu_{D}},$ (14)
$\displaystyle M_{\nu_{R}}=c_{R}S.$ (15)

The mass matrices are written in a basis $f^{c}_{i}M_{ij}f_{j}$, and we have
chosen a phase convention where $v_{10}$ is made real by an $SU(2)_{L}$
rotation. Without loss of generality, we choose to work in a basis where the
matrix $S$ is diagonal. As usual, $v=174.104$ GeV. When charged fermion Yukawa
coupling matrices are bi-diagonalized, we denote the corresponding diagonal
couplings as $Y^{\textrm{diag}}_{U}=\left(y_{u},y_{c},y_{t}\right)$,
$Y^{\textrm{diag}}_{D}=\left(y_{d},y_{s},y_{b}\right)$, and
$Y^{\textrm{diag}}_{E}=\left(y_{e},y_{\mu},y_{\tau}\right)$, for the up-type,
down-type quarks, and charged lepton sectors, respectively. The light neutrino
mass matrix, obtained from the seesaw formula [46, 47, 48, 49, 50], assuming
the typical type-I dominance, is given by

$\displaystyle M_{N}=-M^{T}_{\nu_{D}}M^{-1}_{\nu_{R}}M_{\nu_{D}}.$ (16)

Due to minimal number of parameters, in this work, our focus is on type-I
dominance. A full exploration of the type-I+II case is left for a future
exploration.

It is crucial to identify the intermediate symmetry breaking scale for our
analysis. First, note that the mass of the heaviest right-handed neutrino is
given by $M_{3}=v_{R}(Y_{126})_{33}$, see Eq. (7). Moreover, this model has
the specific property that the top quark mass is solely determined by the
matrix $S$ that provides $m_{t}=v^{u}_{126}(Y_{126})_{33}$, see Eq. (3). This
can be explained as follows. First, note that at the GUT scale
$m_{t}/m_{b}\sim 75$ (see Table I). Moreover, the dominant contribution to the
bottom-quark mass in Eq. (14) comes from $\left(D_{33}+r_{1}S_{33}\right)$,
which demands $\left(D_{33}+r_{1}S_{33}\right)/S_{33}\sim 10^{-2}$ to
correctly reproduce the above-mentioned ratio between the top and the bottom
quark masses. Hence, $D_{33},r_{1}S_{33}\ll S_{33}\sim m_{t}$, as anticipated.
As will be shown later, the GUT scale value of the top quark mass is
$m_{t}\approx 85$ GeV. Consequently, the lowest value of $(Y_{126})_{33}$
corresponds to the maximum value allowed for $v^{u}_{126}\lesssim 174$ GeV,
which provides $(Y_{126})_{33}^{\mathrm{low}}=0.488$. On the other hand, to
ensure perturbativity of the Yukawa couplings, we allow the highest value of
$(Y_{126})_{33}^{\mathrm{high}}=2$. Since the mass of the heaviest right-
handed neutrino is given by $M_{3}=v_{R}(Y_{126})_{33}$, the lowest and the
highest symmetry breaking scales producing the cosmic string network are
determined to be

$\displaystyle v_{R}=\begin{cases}v_{R}^{\textrm{min}}=0.5\times M_{3}.\\\
v_{R}^{\textrm{max}}=2.05\times M_{3}.\end{cases}$ (17)

As will be discussed, a fit to fermion mass spectrum prefers a large $M_{3}$
value, hence is highly constrained from GW detectors. (A choice of
$(Y_{126})_{33}^{\mathrm{high}}=\sqrt{4\pi}$ would lead to
$v_{R}^{\textrm{min}}=0.28\times M_{3}$.)

##  3 Fit to the fermion masses and mixings

In this section, we provide numerical fits to fermion masses and mixing
utilizing the mass matrices derived in the previous section. We substantially
improve the fitting procedure compared to Ref. [23]. Our fitting procedure is
summarized as follows:

  1. 1.

The fermion mass matrices Eqs. (11)-(16) contain only a limited number of
parameters: 15 magnitudes and 12 phases. All these parameters are randomly
varied at the GUT scale, which for fitting purposes, we fix at
$M_{GUT}=2\times 10^{16}$ GeV.

  2. 2.

A significant improvement over the previous study is that we take into account
the full SM+Type-I seesaw renormalization group equation (RGE) running from
$M_{GUT}$ to $M_{Z}$ scale. For this task, we utilize the package REAP [51],
which integrates out heavy right-handed neutrinos successively at their
respective mass thresholds. This is very crucial due to the extremely
hierarchical mass structure in the model under consideration. Since the
intermediate scale is quite close to the GUT scale, to the precision we are
working, it is good enough to consider only SM+Type-I RGEs. Incorporating
full-fledged RGEs from the intermediate to the GUT scale, which is expected to
provide a small correction, is, however, beyond the scope of this work.

Observables |  Values at $M_{Z}$ scale |  Values at $M_{GUT}$ scale
---|---|---
($\Delta m^{2}_{ij}$ in $eV^{2}$) | Input | Best Fit: NO | Best Fit: IO | NO | IO
$y_{u}/10^{-6}$ | 6.65$\pm$2.25 | 6.65 | 6.71 | 2.85 | 2.87
$y_{c}/10^{-3}$ | 3.60$\pm$0.11 | 3.60 | 3.60 | 1.54 | 1.54
$y_{t}$ | 0.986$\pm$0.0086 | 0.986 | 0.986 | 0.48 | 0.48
$y_{d}/10^{-5}$ | 1.645$\pm$0.165 | 1.645 | 1.675 | 0.73 | 0.74
$y_{s}/10^{-4}$ | 3.125$\pm$0.165 | 3.125 | 3.146 | 1.38 | 1.39
$y_{b}/10^{-2}$ | 1.639$\pm$0.015 | 1.639 | 1.639 | 0.637 | 0.637
$y_{e}/10^{-6}$ | 2.7947$\pm$0.02794 | 2.7947 | 2.7899 | 2.8873 | 2.8817
$y_{\mu}/10^{-4}$ | 5.8998$\pm$0.05899 | 5.8998 | 5.9021 | 5.924 | 5.894
$y_{\tau}/10^{-2}$ | 1.0029$\pm$0.01002 | 1.0029 | 1.0012 | 0.985 | 0.989
$\theta_{12}^{\textrm{CKM}}/10^{-2}$ |  $22.735\pm$0.072 | 22.735 | 22.739 | 22.73 | 22.74
$\theta_{23}^{\textrm{CKM}}/10^{-2}$ | 4.208$\pm$0.064 | 4.208 | 4.204 | 4.79 | 4.79
$\theta_{13}^{\textrm{CKM}}/10^{-3}$ | 3.64$\pm$0.13 | 3.64 | 3.64 | 4.15 | 4.15
$\delta^{\textrm{CKM}}$ | 1.208$\pm$0.054 | 1.208 | 1.204 | 1.207 | 1.204
$\Delta m^{2}_{21}/10^{-5}$ | 7.425$\pm$0.205 | 7.425 | 7.433 | 9.714 | 950.84
$\Delta m^{2}_{31}/10^{-3}$ (NO) | 2.515$\pm$0.028 | 2.515 | - | 12.909 | -
$\Delta m^{2}_{32}/10^{-3}$ (IO) | -2.498$\pm$0.028 | - | -2.497 | - | -12.515
$\sin^{2}\theta_{12}$ | 0.3045$\pm$0.0125 | 0.3045 | 0.3053 | 0.308 | 0.177
$\sin^{2}\theta_{23}$ (NO)∗ | 0.5705$\pm$0.0205 | 0.5726 | - | 0.484 | -
$\sin^{2}\theta_{23}$ (IO)∗ | 0.576$\pm$0.019 | - | 0.5819 | - | 0.542
$\sin^{2}\theta_{13}$ (NO) | 0.02223$\pm$0.00065 | 0.02223 | - | 0.007 | -
$\sin^{2}\theta_{13}$ (IO) | 0.02239$\pm$0.00063 | - | 0.02238 | - | 0.0223
$\chi^{2}$ | - | $3\times 10^{-8}$ | 2.77† | - | -

Table I: Best fit values of the observables without restricting the
$U(1)_{B-L}$ breaking scale. ∗Note that experimental measurements of
$\theta_{23}$ have two local minimum, one corresponds to
$\theta_{23}<45^{\circ}$ and the other $\theta_{23}>45^{\circ}$, see Ref.
[52]. The Yukawa couplings, $y_{f}$, at the GUT scale are defined below Eq.
(15). †It is important to point out that for the inverted neutrino mass
ordering, a total $\chi^{2}\approx 2.7$ is entirely coming from
$\left(\Delta\chi^{2}\right)_{\theta_{23}}\approx 2.7$, for details, see Ref.
[52].

  3. 3.

We perform a minimization of a $\chi^{2}$-function that includes: 18
observables: 6 quark masses, 3 quark mixing angles, 1 CKM phase, 3 charged
lepton masses, 2 neutrino mass squared differences, and 3 mixing angles in the
neutrino sector. Since the Dirac CP phase in the neutrino sector has not been
measured yet, we do not include it in the $\chi^{2}$-function. All these
observables are fitted at the $M_{Z}=91.8176$ GeV scale after running them
down from the GUT scale. Low scale measured values of the observables in the
charged and neutral fermion sectors are taken from Refs. [53] and [52, 54],
respectively.

  4. 4.

Finally, for a chosen benchmark scenario (with $M_{3}\sim 4\times 10^{13}$
GeV), a Markov chain Monte Carlo (MCMC) analysis is performed to explore the
parameter space.

To demonstrate the viability of this theory, first, we perform a fit to
fermion masses and mixings without restricting the intermediate $U(1)_{B-L}$
breaking scale. The resulting best fit values are presented in Table I. It is
interesting to point out that the minimal $SO(10)$ Yukawa sector admits normal
as well as inverted mass ordering for the neutrinos. Allowing a solution to
the inverted neutrino mass ordering is a very distinct feature of the current
model, which was not realized in the literature when this model was first
introduced in Ref. [23]. A collection of additional predicted quantities in
the neutrino section are given in Table II. Moreover, best fit parameters (at
the GUT scale) are given in Eqs. (62)-(68) and in Eqs. (69)-(75) for normal
ordering (NO) and inverted ordering (IO) cases, respectively.

Quantity | Best fit prediction
---|---
| NO | IO
$\left(m_{1},m_{2},m_{3}\right)$ |  $\left(0.00014,0.0086,0.0501\right)$ eV |  $\left(0.04922,0.04997,0.00038\right)$ eV
$\left(\sum_{i}m_{i},m_{\beta},m_{\beta\beta}\right)$ |  $\left(0.0589,0.0088,0.0016\right)$ eV |  $\left(0.099,0.041,0.033\right)$ eV
$\left(\delta,\varphi_{1},\varphi_{2}\right)$ | $\left(326.4,109.0,94.7\right)^{\circ}$ | $\left(209.5,164.4,343.4\right)^{\circ}$
$\left(M_{1},M_{2},M_{3}\right)$ |  $\left(2.13\times 10^{5},6.46\times 10^{11},2.28\times 10^{14}\right)$ GeV |  $\left(1.31\times 10^{4},6.42\times 10^{11},2.37\times 10^{14}\right)$ GeV

Table II:  Best fit predictions of several quantities in the neutrino sector.
We have defined
$U_{\mathrm{PMNS}}=V_{\mathrm{PMNS}}\;\mathrm{diag}\\{e^{-\varphi_{1}/2},e^{-\varphi_{2}/2},1\\}$,
where $V_{\mathrm{PMNS}}$ is the CKM like matrix that contains the Dirac
phase. Here $m_{\beta}=\sqrt{\sum_{i}|U_{ei}|^{2}m_{i}^{2}}$ is the effective
mass parameter for beta decay and $m_{\beta\beta}=|\sum_{i}U_{ei}^{2}m_{i}|$
is the effective mass parameter for neutrinoless double beta decay.

Figure 1:  Plots depicting the importance of the RGE running on the tau lepton
mass and neutrino mass-squared as well as mixing parameters. Scales, where the
heavy right-handed neutrinos $\nu_{R_{i}}$ decouple from the theory, are
presented via different shades of gray.

As can be seen from Table II, this theory predicts right-handed neutrinos that
show a very strong hierarchy, which is a generic feature of this model. To
have an understanding of this behavior, let us consider the best fit
parameters given in Eqs. (62)-(68) (this behavior, however, is realized for
any parameter set that provides a realistic fermion fit). One finds that
$(D_{22},\\{D_{33},r_{1}S_{33}\\})/v\sim(10^{-4},10^{-3})$ ($r_{1}S_{33}$ is
smaller in order than $D_{22}$). These are the correct orders to reproduce the
Yukawas in the down sector, namely, $y_{s,\mu}\sim 10^{-4}$ and
$y_{b,\tau}\sim 10^{-2}$ for the second and the third generations,
respectively. However, these are too small to provide correct masses to the
up-type quarks. Consequently, the charm- and top-quark masses are entirely
coming from the matrix $S$. This can be seen by looking at the relevant
entries, i.e., $(S_{22},S_{33})/v\sim(10^{-3},0.5)$, as needed to replicate
$y_{c}\sim 10^{-3}$ and $y_{t}\sim 0.5$. On the other hand, $D_{11}/v\sim
10^{-6}$ (in association with other relevant entries in $D$ and $A$ matrices
not listed here) reproduces the first generation fermion masses.

Moreover, naively speaking, neutrino mass from the type-I seesaw, with order
one coupling demands the seesaw scale to be $M_{3}\sim 10^{14}$ GeV. Since
right-handed neutrino mass matrix is fully determined by the matrix $S$, which
is diagonal, $S_{22}/S_{33}\sim y_{c}/y_{t}\sim 10^{-3}$ uniquely determines
the second lightest right-handed neutrino mass $M_{2}\sim 10^{-3}M_{3}\sim
10^{11}$ GeV.

With the entries of order
$(D_{11},\\{D_{12,21},A_{12}\\},\\{D_{13,3,1},A_{13}\\})/v\sim(10^{-6},10^{-5},10^{-4})$,
the structure of the Dirac Yukawa coupling of the neutrino is also
hierarchical. As a result, from the type-I seesaw formula Eq. (16), elements
in $M_{N}$, associated with $M^{-1}_{2,3}$ provide too small (when compared to
the 2-3 sector of $M_{N}$) contributions to be compatible with neutrino
oscillation data. A proper fit to neutrino data requires non-significant
entries in the first row and first column of $M_{N}$ demanding $M_{1}\ll
M_{2,3}$; meaning, all the dominant terms in the first row and first column
result from terms that are suppressed only by $M_{1}^{-1}$ in the seesaw
formula, which is, however, not true in general, in the 2-3 block. For the
fits provided, we find $M_{1}\sim 10^{5}\ll M_{2,3}$ GeV to ensure a
consistent fit to the neutrino sector for the case of normal mass ordering. On
the other hand, for the inverted mass ordering, one must make sure that
entries in the first row and column are more dominating compared to the 2-3
sector. This is precisely obtained by ensuring the $M_{1}$ mass to be slightly
smaller than the normal ordering solution. As can be seen from our fits (see
Table II), $M^{\mathrm{IO}}_{1}\sim 10^{4}$ GeV as required to reproduce
inverted neutrino mass ordering, which is one order smaller than the solution
we obtain for normal ordering, $M^{\mathrm{NO}}_{1}\sim 10^{5}$ GeV. These
features in the fermion mass fit as explained above are very distinct from
other $SO(10)$ models proposed in the literature.

Before concluding this section, we emphasize on the impact of the RGE running
on the neutrino parameters. This is clearly demonstrated in Fig. 1 for the two
mass-squared differences and three mixing angles in the neutrino sector. This
figure illustrates that a naive fit to the fermion masses and mixing by
ignoring the RGEs can lead to significant errors. Moreover, the top plots in
Fig. 1 show that the running effects on the Tau Yukawa coupling are also very
crucial and cannot be neglected.

##  4 Symmetry breaking chain and topological defects

Since GUT unifies all the known forces present in the SM, it comprises of a
larger symmetry group $G$ that contains the SM gauge group ($SU(3)_{C}\times
SU(2)_{L}\times U(1)_{Y}\equiv G_{321}\subset G$). A series of spontaneous
summery breaking (SSB) must lead to the SM at low energies:

$\displaystyle G\to H\to...\to G_{321}\to G_{31},$ (18)

where $G_{31}\equiv SU(3)_{C}\times U(1)_{em}$, $H\subset G$, and $...$
represents the possibility of multi-step SSB. Subsequent phase transitions in
the early Universe, as shown in Eq. (18), have direct observational
consequences; since various topological defects are produced during
cosmological phase transitions via Kibble mechanism [55] that depend on the
topology of the vacuum manifold.

Types of topological defects for $G\to H$ are determined by the $n^{th}$
homotopy group $\pi_{n}(\mathcal{M})$ of the vacuum manifold
$\mathcal{M}=G/H$, and defects are formed for $\pi_{k}(\mathcal{M})\neq 0$
[56]. For example, $\pi_{0}(\mathcal{M})\neq 0$, $\pi_{1}(\mathcal{M})\neq 0$,
and $\pi_{2}(\mathcal{M})\neq 0$ lead to domain walls, strings, and monopoles,
respectively. Symmetry breaking patterns that lead to either monopoles or
domain walls (or both) are ruled out. If stable topological defects of these
types are formed, they would dominate the energy density of the Universe and
overclose it, hence must be eliminated to create our observable Universe.

To solve a number of cosmological puzzles (flatness, homogeneity, and isotropy
observed in the CMB data, including the monopole problem), the idea of
inflation was put forward [57, 58, 59, 60]; for a review, see Ref. [61]. If
these unwanted topological defects are generated before the inflationary era,
they will be diluted away by the Universe’s rapid expansion. Production of
monopoles is a generic feature of GUTs due to $U(1)\subset G_{321}$; on the
contrary, domain walls may or may not appear. To be consistent with
observation, in this work, we assume inflation to take place after the last
stage of the SSB chain that involves any unwanted defect. Inflation can be
successfully implemented with the addition of a GUT singlet, and the detail is
not relevant to our discussion here.

Cosmic strings, however, are not ruled out by measurements and may lead to
stimulating cosmological signatures. As aforementioned, fascinatingly, pulsar-
timing arrays NANOGrav [19], PPTA [20], EPTA [21], and IPTA [22] might have
already seen evidence of a primordial gravitational wave background that could
originate from the last stage spontaneous symmetry breaking in $SO(10)$ GUT.
In such a scenario, determining the intermediate symmetry breaking scale from
low energy fermion masses and mixings data as well as gravitational wave
signal strongly constrains the GUT scale, which subsequently has a profound
effect on determining the proton lifetime. To accommodate cosmic strings
within a minimal (i.e., requiring the least number of Higgses and lowest
dimensional representations) framework, we consider the following symmetry
breaking chain:

$\displaystyle SO(10)$
$\displaystyle\xrightarrow[54_{H}]{M_{X}}SU(4)_{C}\times SU(2)_{L}\times
SU(2)_{R}\times D\equiv G_{422D}$ (19)
$\displaystyle\xrightarrow[45_{H}]{M_{I}}SU(3)_{C}\times SU(2)_{L}\times
SU(2)_{R}\times U(1)_{B-L}\equiv G_{3221}$ (20)
$\displaystyle\xrightarrow[126_{H}]{M_{II}}SU(3)_{C}\times SU(2)_{L}\times
U(1)_{Y}\equiv G_{321}$ (21)
$\displaystyle\xrightarrow[10_{H}+126_{H}]{M_{EW}}SU(3)_{C}\times
U(1)_{em}\equiv G_{31}$ (22)

such that $M_{X}>M_{I}>M_{II}>M_{EW}$ (here $M_{X}\equiv M_{\mathrm{GUT}}$,
$M_{I}\equiv M_{\mathrm{PS}}$, $M_{II}\equiv M_{\mathrm{LR}}$).

At the first step, the GUT symmetry is spontaneously broken by the VEV of
$54_{H}$ into Pati-Salam symmetry with the discrete $D$-parity [62]. In the
next symmetry breaking scale, the VEV of $45_{H}$ breaks Pati-Salam symmetry
into left-right symmetry. At this stage, the $D$-parity is also broken since
$45_{H}$ contains a parity-odd singlet under $G_{3221}$. The remaining
symmetry is then subsequently broken to the Standard Model gauge group via the
VEV of $126_{H}$ Higgs that contains a singlet under $G_{321}$. To be precise,
the breaking by $126_{H}$ leaves a remnant $Z_{2}$ symmetry (it is the $Z_{2}$
subgroup of $Z_{4}$ , the center of SO(10) [63], which is not broken by tensor
representations, regardless of the breaking chain). If this $Z_{2}$ were
absent, embedded strings would appear, which are in general unstable under
small perturbations; for a list of breaking chains leading to embedded
strings, see, for example, Ref. [64].

Since stable monopoles are formed during the Pati-Salam symmetry breaking, we
assume inflation to take place right after $M_{I}$ breaking scale that gets
rid of all unwanted topological defects. The last symmetry breaking scale
(prior to the appearance of the SM gauge group) is carried out by $126_{H}$
Higgs, and topologically stable strings with mass scale related to
$M_{II}=v_{R}$ are formed during this phase transition. Finally, the SM gauge
symmetry is broken by the VEV of the Higgs doublet, which comes from a
combination of the bi-doublets in $10_{H}$ and $126_{H}$ dimensional Higgses.
For related works on GWs arising from cosmic strings within $SO(10)$
framework, see, for example, Refs. [65, 66, 67, 68, 69, 70, 71, 72, 73].

##  5 Gravitational waves from cosmic strings

###  5.1 Observational signals from strings

We assume the ideal Nambu-Goto string approximation, where the dominant
radiation emission is in the form of GWs [74] (for different modeling
approaches, see, for example, Ref. [75]). Macroscopic properties of strings
are characterized by their energy per unit length $\mu$, known as the string
tension parameter, which is given by [76],

$\displaystyle\mu\sim 2\pi v^{2}_{R}B(2\lambda/e^{2}),$ (23) $\displaystyle
B(x)\sim\begin{cases}1,\;\;\;x=1.\\\ 1.04x^{0.195},\;\;\;10^{-2}<x\ll
1.\end{cases}$ (24)

With the above, we obtain $G\mu=4.22\times 10^{-38}v^{2}_{R}$, where
$G=m^{-2}_{\textrm{pl}}$ denotes the Newton’s constant
($m_{\textrm{pl}}=1.22\times 10^{19}$ GeV). This quantity, $G\mu$,
parametrizes the gravitational interactions of the string. While the existing
loops oscillate and emit radiation dominantly to GWs, the expansion of the
Universe causes the long strings to intersect to form new loops. The power
radiated in GWs is given by

$\displaystyle P_{GW}=\Gamma G\mu^{2},$ (25)

where $\Gamma=50$ [77, 78, 79, 80, 81] is a dimensionless quantity. For Nambu-
Goto strings, only large loops give a dominant contribution to GWs, which has
the initial length of $\ell_{i}=\alpha t_{i}$, where $\alpha$ characterizes
the loop size at the time of formation, and numerical simulation typically
corresponds to $\alpha=0.1$ [80, 75]. Due to energy loss, their lengths
decrease with time $\ell=\alpha t_{i}-\Gamma G\mu(t-t_{i})$, until loops
disappear completely. The frequency associated to the GW from the total energy
loss from a loop is $f_{k}=2k/\ell$, where $k=1,2,3,...$ is the mode number.
Then the relic GW background from a cosmic string network has the following
form [81, 82]:

$\displaystyle\Omega_{GW}(f)=\sum_{k}\Omega_{GW}^{(k)}(f),$ (26)

where the sum over all modes are performed, and the explicit expression of
$\Omega_{GW}^{(k)}(f)$ is given by

$\displaystyle\Omega_{GW}^{(k)}(f)=\frac{1}{\rho_{c}}\frac{2k}{f}\frac{\mathcal{F}_{\alpha}\Gamma^{(k)}G\mu^{2}}{\alpha(\alpha+\Gamma
G\mu)}\int^{t_{0}}_{t_{F}}d\widetilde{t}\frac{c_{eff}(t^{(k)}_{i})}{{t^{(k)}_{i}}^{4}}\left(\frac{a(\widetilde{t})}{a(t_{0})}\right)^{5}\left(\frac{a(t^{(k)}_{i})}{a(\widetilde{t})}\right)^{3}\Theta\left(t^{(k)}_{i}-t_{F}\right),$
(27)

here $\Gamma^{(k)}=\Gamma k^{-4/3}/3.6$ and the integral is taken over the GW
emission time $\widetilde{t}$, $t_{F}$ is the time of the formation of cosmic
string network, and the formation time of loops contributing with $k$-th mode
number is defined as

$\displaystyle t^{(k)}_{i}(\widetilde{t},f)=\frac{1}{\alpha+\Gamma
G\mu}\left(\frac{2k}{f}\frac{a(\widetilde{t})}{a(t_{0})}+\Gamma
G\mu\widetilde{t}\right).$ (28)

Furthermore, $\rho_{c}=3H^{2}_{0}/(8\pi G)$ is the critical energy density of
the Universe, $\mathcal{F}_{\alpha}=0.1$ [80], and $c_{eff}=5.7$ (0.5) for
radiation (matter) dominance are determined numerically.

Figure 2:  Gravitational wave signals from stable cosmic string networks as a
function of frequency observed today together with the corresponding
sensitivities of the indicated observatories. See text for details.

Gravitational wave signals from stable cosmic string networks as a function of
frequency observed today together with the sensitivities of various
observatories are presented in Fig. 2 for three chosen values of $G\mu$. CMB
measurements directly provide constraints on string tension $G\mu<1.1\times
10^{-7}$ for Nambu-Goto strings [83, 84], which corresponds to $M_{II}\lesssim
2\times 10^{15}$ GeV, hence, GUT scale appearance of strings are ruled out (if
no fine-turning is done, then the current experimental bounds on proton
lifetime requires the GUT scale to be $M_{\mathrm{GUT}}\gtrsim 5\times
10^{15}$ GeV). Recently, LIGO updated their O3 data [85] that rules out
$G\mu\gtrsim 10^{-8}$ in the $f\sim\mathcal{O}(10)$ Hz regime, providing
$M_{II}\lesssim 5\times 10^{14}$ GeV. Moreover, non-observation of SGWs in the
$f\sim$nHz region by EPTA [86] and 11 year data of NANOGrav [87] rule out
string formation up to almost $10^{14}$ GeV.

On the other hand, 12.5 years of NANOGrav’s data indicates SGWs for string
tension in the range $G\mu\sim[4\times 10^{-11}-10^{-10}]$ at $68\%$
confidence level (CL), while $95\%$ CL corresponds to the range
$G\mu\sim[2\times 10^{-11}-3\times 10^{-10}]$ [88]. As seen from Fig. 2, since
GW spectrum of cosmic strings is very flat over a wide frequency range, the
current hints of SGW backgrounds at PTAs for sure can be confirmed by the
measurement of various upcoming experiments. Frequencies in the window
$10^{-5}-10^{3}$ Hz will be probed by ground-based interferometers: Einstein
Telescope (ET) [89] and Cosmic Explorer (CE) [90]; space-based laser
interferometers Laser Interferometer Space Antenna (LISA) [91], Taiji [92],
TianQin [93], BBO [94], and DECIGO [95]); atomic interferometers MAGIS [96],
AEDGE [97], AION [98]. Whereas pulsar timing arrays such as EPTA and NANOGrav
are currently covering the nHz regime. Future pulsar timing array SKA [99]
will also cover a similar frequency regime which is sensitive to much smaller
GW energy densities. Furthermore, large-scale surveys of stars such as Gaia
[100] and the proposed upgrade, THEIA [101] will be covering $1-100$ nHz
regime, where the latter is sensitive to GW energy densities as low as
$10^{-16}$. These experiments are sensitive to GWs arising from cosmic strings
with breaking scales $v_{R}\gtrsim 10^{9}$ GeV. It is remarkable that these GW
experiments will probe GUT scale physics and can rule out various symmetry
breaking chains; such large energy scales, however, are impossible to probe in
the colliders.

Figure 3:  Best fit values of the parameters $(A,\gamma)$ for the various PTAs
with a fixed frequency, $f=5.6$ nHz. Good fits (i.e., $\sum\chi^{2}<2$) are
obtained for $a\sim[1,6]$. See text for details.

The observation of GWs at the PTAs NANOGrav [19], PPTA [20], EPTA [21], and
IPTA [22] hint towards gravitational wave signals with the energy density
given by [72],

$\displaystyle\left(\Omega_{\rm GW}(f)h^{2}\right)_{\rm PTA}\approx 2.02\cdot
10^{-10}\left(\frac{A}{10^{-15}}\right)^{2}\times\left(\frac{f}{f_{\rm
yr}}\right)^{5-\gamma},$ (29)

where $f_{\rm yr}=1yr^{-1}$. $A$ is the characteristic strain amplitude and
$\gamma$ is the spectral index. In fact, NANOGrav [19] finds strong evidence
of a stochastic process modeled as a power-law, $f^{-2/3}$, with common
amplitude and spectral slope across pulsars. This power-law spectrum of the
characteristic GW strain corresponds to $\gamma=13/3$, which is motivated by
searches for mergers of supermassive black hole binaries. Within this fiducial
model, the amplitude parameter $A$ that relates the correlation between
pulsars has a median $1.92\times 10^{-15}$ at $95\%$ confidence level. Without
fixing the value of $\gamma$, other pulsar timing arrays also find the
amplitude of similar values in the same frequency range. Fitted values of
these pulsar timing array experiments are summarized in Table III.

Experiment | $A\times 10^{15}$ | $\gamma$
---|---|---
NANOGrav | $1.92_{-0.55}^{+0.75}$ | $4.333$
PPTA | $2.82_{-1.16}^{+0.73}$ | $4.11_{-0.41}^{+0.52}$
EPTA | $5.13_{-2.73}^{+4.20}$ | $3.78_{-0.59}^{+0.69}$
IPTA | $3.8_{-2.5}^{+6.3}$ | $4.0_{-0.9}^{+0.9}$
Table III: Experimentally allowed values of the characteristic strain
amplitude and the spectral index at $95\%$ CL.

Intriguingly, these data from various pulsar timing array experiments strongly
favor stochastic gravitational-wave background in the nHz range with energy
density

$\displaystyle\left(\Omega_{\rm GW}h^{2}\right)_{\rm PTA}\sim 10^{-10},$ (30)

which is in the ballpark, as expected from $SO(10)$ GUTs.

Figure 4:  A combined fit (cyan curve) to NANOGrav, PPTA, EPTA, and IPTA data
in the $(A,\gamma)$ plane. The diamond and the star shaped points are the two
chosen benchmark scenarios of the proposed model. See text for details.

By using the experimental data shown in Table III, we present the best fit
values of the parameters $(A,\gamma)$ in Fig. 3 for the separate PTAs by
fixing $\Omega_{\rm GW}h^{2}$ in the range $8\times 10^{-11}-10^{-9}$. For
this analysis, we have fixed the frequency at $f=5.6$ nHz, which is within the
frequency range $f\in(2.5,12)$ nHz considered by NANOGrav [19] (see also Ref.
[88]). For our numerical study, we allow a $5\%$ fluctuations of
$\gamma_{\mathrm{NANOGrav}}=13/3$.

To illustrate that the proposed model can simultaneously be consistent with
these experimental measurements, we perform a combined fit to these data. In
doing so, first, we symmetrize the experimental errors on the two observables,
i.e., $A^{(k)}\pm\sigma^{(k)}_{A}$ and $\gamma^{(k)}\pm\sigma^{(k)}_{\gamma}$;
$k=$NANOGrav, IPTA, PPTA, and EPTA. Since these experimental values of the
same physical quantity are uncorrelated, we get the best estimate via [102]

$\displaystyle A^{\rm
comb.}=\frac{\sum_{k}\frac{A^{(k)}}{\left(\sigma^{(k)}_{A}\right)^{2}}}{\sum_{k}\frac{1}{\left(\sigma^{(k)}_{A}\right)^{2}}},\;\;\;\sigma^{\rm
comb.}_{A}=\frac{1}{\sqrt{\sum_{k}\frac{1}{\left(\sigma^{(k)}_{A}\right)^{2}}}},$
(31) $\displaystyle\gamma^{\rm
comb.}=\frac{\sum_{k}\frac{\gamma^{(k)}}{\left(\sigma^{(k)}_{\gamma}\right)^{2}}}{\sum_{k}\frac{1}{\left(\sigma^{(k)}_{\gamma}\right)^{2}}},\;\;\;\sigma^{\rm
comb.}_{\gamma}=\frac{1}{\sqrt{\sum_{k}\frac{1}{\left(\sigma^{(k)}_{\gamma}\right)^{2}}}}.$
(32)

The result of our numerical fit is presented in Fig. 4. This figure shows that
an excellent fit to the data from PTAs can be obtained for
$\Omega_{\mathrm{GW}}h^{2}\sim 10^{-10}$. We find that for $f=5.6$ nHz, in the
range of $\Omega_{\mathrm{GW}}h^{2}\in(2,6)\times 10^{-10}$, all these PTAs
observations can be satisfied within $2\sigma$ CL. This range roughly
corresponds to $G\mu\in(4.9,6.9)\times 10^{-11}$. In other words, the seesaw
scale from the observation of GWs is highly restricted within
$v_{R}\in(3.4,4.1)\times 10^{13}$ GeV.

###  5.2 Fit to fermion mass with fixed intermediate scale

As discussed in the previous sub-section, data from PTAs prefers $v_{R}\sim
4\times 10^{13}$ GeV. On the contrary, as we have seen in Sec. 2 (see Table
II), this theory, for both the NO and IO scenarios, prefers $M_{3}\sim 2\times
10^{14}$ GeV, which is in direct conflict (assuming perturbative Yukawa
couplings) with observations of pulsar-timing arrays. To lower the seesaw
scale from the best fit value, we perform an extensive numerical analysis and
find that for $M_{3}$ as low as $M_{3}\sim 7\times 10^{12}$ GeV,
$\chi^{2}\sim\mathcal{O}(10)$, which arises since the up-quark mass cannot be
fitted to its observed value. This corresponds to more than $3\sigma$
discrepancy in the up-quark mass measurement, which is unacceptable. The
$\Delta\chi^{2}$ contribution to the up-quark mass reduces to about $2\sigma$
for $M_{3}=2\times 10^{13}$ GeV, which is in the acceptable range. Starting
from this value, an excellent fit to all observables can be found up to about
$M_{3}\sim 10^{15}$ GeV. Above this scale, neutrino masses cannot be correctly
incorporated. Hence, the low energy observables are consistent with
experimental measurements provided that $M_{3}=[2\times 10^{13}-10^{15}]$ GeV.

To be consistent with both low energy fermion masses and mixings data as well
as results of pulsar-timing arrays, as a benchmark, we take $M_{3}=4\times
10^{13}$ in what follows, which corresponds to $v_{R}\in(2,8)\times 10^{13}$
GeV depending on the value of the chosen Yukawa coupling, see Eq. (17). Two
benchmark values, $v_{R}=2\times 10^{13}$ and $v_{R}=4\times 10^{13}$ GeV
corresponding to $G\mu=1.688\times 10^{-11}$ (diamond shaped point) and
$G\mu=6.754\times 10^{-11}$ (star shaped point) are highlighted in Fig. 4.

Other than fixing the seesaw scale, the fitting procedure performed in this
section is identical to the prescriptions outlined in Sec. 2. This analysis is
performed by assuming normal ordering; however, repeating it for inverted
ordering is expected to result in similar conclusions. The best fit result of
the observables is presented in Table IV. This solution corresponds to a total
$\chi^{2}=0.52$, which receives the largest contribution from the up-quark
mass corresponding to $\Delta\chi^{2}\sim 0.34$. The best fit parameters are
presented in Eqs. (76)-(82). Proton decay branching fraction predicted from
this best fit is summarized in Table V.

Observables | Best fit ($M_{Z}$ scale) | Observables | Best fit ($M_{Z}$ scale)
---|---|---|---
$y_{u}/10^{-6}$ | 5.33 | $\theta_{12}^{\textrm{CKM}}/10^{-2}$ | $22.737$
$y_{c}/10^{-3}$ | 3.60 | $\theta_{23}^{\textrm{CKM}}/10^{-2}$ | 4.208
$y_{t}$ | 0.986 | $\theta_{13}^{\textrm{CKM}}/10^{-3}$ | 3.63
$y_{d}/10^{-5}$ | 1.698 | $\delta^{\textrm{CKM}}$ | 1.207
$y_{s}/10^{-4}$ | 3.126 | $\Delta m^{2}_{21}/10^{-5}(eV^{2})$ | 7.403
$y_{b}/10^{-2}$ | 1.639 | $\Delta m^{2}_{31}/10^{-3}(eV^{2})$ | 2.517
$y_{e}/10^{-6}$ | 2.7942 | $\sin^{2}\theta_{12}$ | 0.3028
$y_{\mu}/10^{-4}$ | 5.8961 | $\sin^{2}\theta_{23}$ | 0.5723
$y_{\tau}/10^{-2}$ | 1.0028 | $\sin^{2}\theta_{13}$ | 0.02233

Table IV: Best fit values of the observables by restricting the intermediate
symmetry breaking scale. In particular, we fix $M_{3}=4\times 10^{13}$ GeV
that results in $\chi^{2}=0.52$, hence representing an excellent agreement
with measurements. See text for details.

$p$ decay modes | BR[$\%$] predictions
---|---
$\left(p\rightarrow e^{+}\pi^{0},p\rightarrow e^{+}K^{0},p\rightarrow e^{+}\eta\right)$ | $\left(32.4,0.52,0.03\right)$
$\left(p\rightarrow\mu^{+}\pi^{0},p\rightarrow\mu^{+}K^{0},p\rightarrow\mu^{+}\eta\right)$ | $\left(0.42,4.3,0.002\right)$
$\left(p\rightarrow\overline{\nu}\pi^{+},p\rightarrow\overline{\nu}K^{+}\right)$ | $\left(61.5,0.81\right)$

Table V: Proton decay branching fraction predictions of the best fit with the
fixed seesaw scale, $M_{3}=4\times 10^{13}$ GeV.  Figure 5:  Right-handed
neutrino masses in the benchmark scenario resulting from MCMC analysis. The
vertical solid (dashed) line corresponds to $1\sigma$ ($2\sigma$) highest
posterior Density (HPD) for the masses $M_{i}$. Figure 6:  Neutrinoless double
beta decay in the benchmark scenario resulting from MCMC analysis.
Cosmological exclusion bound (on the sum of the neutrino masses) is taken from
Ref. [103].

For this benchmark scenario, after finding the best fit, we perform a MCMC
analysis using an adaptive Metropolis-Hastings algorithm, the results of which
are presented in Figs. 5, 6, 7, and 8 for selected physical quantities. Fig. 5
shows the Histogram plot of the right-handed neutrino masses $M_{i}$ with
their associated $1\sigma$ and $2\sigma$ highest posterior Density (HPD)
intervals. In Fig. 6, we present the neutrinoless double beta decay parameter
along with sensitives from a number of experiments. This includes current
bound from KamLAND-Zen [104] and future sensitives of next-generation
experiments such as GERDA Phase II [105] and nEXO [106]. Fig. 7 demonstrates
the interrelationship between the two most dominating proton decay modes in
this theory. Finally, Fig. 8 depicts the correlation between the 2-3 mixing
angle and the CP violating Dirac phase in the neutrino sector, where the
latter is yet to be measured in the experiments. As mentioned above, in our
analysis, we do not include this phase in the $\chi^{2}$-function.

Figure 7:  Proton decay branching ratios of the two most dominating channels
in the benchmark scenario resulting from MCMC analysis.  Figure 8:
Correlation between the Dirac CP phase and $\theta_{23}$ in the neutrino
sector in the benchmark scenario resulting from MCMC analysis.

##  6 Proton decay and Gauge coupling unification

Here, we discuss the constraining arising from proton decay and gauge coupling
unification in the proposed model.

###  6.1 Proton decay

Nucleon decay is the smoking gun signal of grand unification (for a recent
review, see Ref. [107]). Proton decay in nonsupersymmetric grand unified
models is dominated by the exchange of the superheavy gauge bosons,
interactions of which violate baryon number conservation. Since these gauge
bosons have masses at the GUT scale, observation of proton decay would probe
the grand unification scale. In our model, the GUT scale is intrinsically
connected to the intermediate symmetry breaking scales. As we will show, there
is a non-trivial correlation between the proton lifetime and the observation
of gravitational waves within the proposed model.

The superheavy gauge bosons lead to proton decays of the form: $p\to
M+\bar{l}$, where the meson $M$ can be $\pi^{0}$, $\pi^{+}$, $K^{0}$,
$K^{+},\eta$ and the lepton can be $l=e,\mu,\nu$ [108]. The partial decay
width for such nucleon decay process is given by the following formula:

$\displaystyle\Gamma(p\to
M+\bar{l})=\frac{m_{p}}{32\pi}\left[1-\left(\frac{m_{M}}{m_{p}}\right)^{2}\right]^{2}A_{L}^{2}\bigg{|}\sum_{n}A_{Sn}W_{n}F_{0}^{n}(p\to
M)\bigg{|}^{2}.$ (33)

In the above equation, $m_{p}$, and $m_{M}$ are the masses of the proton and
Mesons, respectively. Furthermore, $W_{n}$ are the Wilson coefficients of the
operators that give rise to a particular decay channel $p\to M+\bar{l}$;
$A_{Sn}$’s are the short-range enhancement factors (to be discussed shortly),
and $F_{0}^{n}=\langle M\rvert{(qq^{\prime})}_{n}q^{\prime\prime}_{L}\lvert
p\rangle$ is the relevant form factor. Here, $q$, $q^{\prime}$ and
$q^{\prime\prime}$ are the light quarks, namely, $u,d,s$.

The most important proton decay channel within our framework corresponds to
$p\to e^{+}\pi^{0}$, which is experimentally the most constraining one. The
current experimental bound from Super-Kamiokande provides $\tau_{p}(p\to
e^{+}\pi^{0})>2.4\times 10^{34}$ yrs. [109]. The expression of the lifetime of
this process is then given by [110, 111]:

$\displaystyle\tau_{p}=$
$\displaystyle\Bigg{[}\frac{m_{p}}{32\pi}\left(1-\frac{m_{\pi^{0}}^{2}}{m_{p}^{2}}\right)^{2}A_{L}^{2}\frac{g_{\mathrm{GUT}}^{4}}{4M_{X}^{4}}\times$
$\displaystyle\bigg{\\{}A_{SR}^{2}\left|\left[V_{1}^{11}V_{3}^{11}+\left(V_{4}V^{\dagger}_{UD}\right)^{11}\left(V_{1}V_{UD}V_{4}^{\dagger}V_{3}\right)^{11}\right]\langle\pi^{0}\rvert(ud)_{R}u_{L}\lvert
p\rangle\right|^{2}$
$\displaystyle+A_{SL}^{2}\left|\left[V_{1}^{11}V_{2}^{11}+\left(V_{1}V_{UD}\right)^{11}\left(V_{2}V^{\dagger}_{UD}\right)^{11}\right]\
\langle\pi^{0}\rvert(ud)_{L}u_{L}\lvert
p\rangle\right|^{2}\bigg{\\}}\Bigg{]}^{-1},$ (34)

where, $M_{X}$ is the mass of the superheavy gauge boson, and $A_{SL}$ and
$A_{SR}$ are the short-range enhancement factors associated with the left-
handed $\mathcal{O}^{d=6}_{L}\left(e^{C},d\right)$ and right-handed
$\mathcal{O}^{d=6}_{R}\left(d^{C},e\right)$ operators, respectively. Moreover,
in our calculation, we use the following values of the matrix elements [112]:

$\displaystyle\langle\pi^{0}\rvert(ud)_{R}u_{L}\lvert p\rangle=-0.131,\ \
\,\,\langle\pi^{0}\rvert(ud)_{L}u_{L}\lvert p\rangle=0.134.$

The rest of the matrix elements required to compute proton decay to other
channels can also be found in Ref. [112]. The mixing matrices $V_{1},V_{2}$
etc are defined as follows: $V_{1}=U_{R}^{T}U_{L}$, $V_{2}=E_{R}^{T}D_{L}$,
$V_{3}=D_{R}^{T}E_{L}$, $V_{4}=D_{R}^{T}D_{L}$, $V_{UD}=U_{L}^{\dagger}D_{L}$,
$V_{EN}=E_{L}^{\dagger}N_{L}$ and $U_{EN}=E^{T}_{R}N_{L}$, where $U,D,E$
define the diagonalizing matrices given by

$\displaystyle U^{\dagger}_{R}\ M_{U}\ U_{L}$ $\displaystyle=$ $\displaystyle
M_{U}^{diag}$ (35) $\displaystyle D^{\dagger}_{R}\ M_{D}\ D_{L}$
$\displaystyle=$ $\displaystyle M_{D}^{diag}$ (36) $\displaystyle
E^{\dagger}_{R}\ M_{E}\ E_{L}$ $\displaystyle=$ $\displaystyle M_{E}^{diag}$
(37) $\displaystyle N^{\dagger}_{R}\ M_{N}\ N_{L}$ $\displaystyle=$
$\displaystyle M_{N}^{diag}.$ (38)

with $N_{R}=N^{*}_{L}$.

Finally, the long distant enhancement factor $A_{L}$ [113] roughly takes the
value $A_{L}\approx 1.2$, and the short range renormalization factors $A_{S}$
are written as follows [114, 115, 116, 117]:

$\displaystyle A_{S}=\prod_{j}^{M_{Z}\leq M_{j}\leq
M_{X}}\prod_{i}\left[\frac{\alpha_{i}\left(M_{j+1}\right)}{\alpha_{i}\left(M_{j}\right)}\right]^{\frac{\gamma_{i}}{b_{i}}},$
(39)

where, as usual, $\alpha_{i}=g_{i}^{2}/4\pi$ and $\gamma_{i}$’s are the
relevant anomalous dimensions, values of which are shown in Table VI. Whereas
the $\beta$-coefficients $b_{i}$ are computed in the next subsection.

Gauge group | Anomalous dimensions
---|---
$\mathcal{O}^{d=6}_{L}\left(e^{C},d\right)$ | $\mathcal{O}^{d=6}_{R}\left(e,d^{C}\right)$
$G_{321}$ | $\\{2,\frac{9}{4},\frac{23}{20}\\}$ | $\\{2,\frac{9}{4},\frac{11}{20}\\}$
$G_{3221}$ | $\\{2,\frac{9}{4},\frac{9}{4},\frac{1}{4}\\}$ | $\\{2,\frac{9}{4},\frac{9}{4},\frac{1}{4}\\}$
$G_{422}$ | $\\{\frac{15}{4},\frac{9}{4},\frac{9}{4}\\}$ | $\\{\frac{15}{4},\frac{9}{4},\frac{9}{4}\\}$
Table VI: Relevant anomalous dimensions for the considered breaking chains.

###  6.2 Unification

The two-loop renormalization group equations (RGE) for the gauge couplings can
be written as,

$\displaystyle\dfrac{d\alpha^{-1}_{i}(\hat{\mu})}{d\ln\hat{\mu}}=-\dfrac{b_{i}}{2\pi}-\sum\limits_{j}\dfrac{b_{ij}}{8\pi^{2}\alpha^{-1}_{j}(\hat{\mu})},$
(40)

where $i,j$ indices refer to different daughter groups of the parent group at
the energy scale $\hat{\mu}$ and we have defined,

$\displaystyle\alpha^{-1}_{i}=\dfrac{4\pi}{g_{i}^{2}}.$ (41)

Using the generalized formula for the beta-function,
$\beta=\hat{\mu}dg/d\hat{\mu}$, one finds [118],

$\displaystyle b_{i}$
$\displaystyle=-\dfrac{11}{3}C_{2}(G_{i})+\dfrac{4}{3}\kappa
S_{2}(F_{i})+\dfrac{1}{6}\eta S_{2}(S_{i}),$ (42) $\displaystyle b_{ij}$
$\displaystyle=-\dfrac{34}{3}\left[C_{2}(G_{i})\right]^{2}\delta_{ij}+\kappa\left[4C_{2}(F_{j})+\dfrac{20}{3}\delta_{ij}C_{2}(G_{i})\right]S_{2}(F_{i})$
$\displaystyle+\eta\left[2C_{2}(S_{j})+\dfrac{1}{3}\delta_{ij}C_{2}(G_{i})\right]S_{2}(S_{i}).$
(43)

Here $S_{2}$ and $C_{2}$ are the Dynkin indices of the representations with
the appropriate multiplicity factors and the quadratic Casimir of a given
representation. For Dirac and Weyl fermions $\kappa=1$ and $\frac{1}{2}$,
respectively, and $\eta=1,2$ for real and complex scalar fields. $G$, $F$ and
$S$ stand for gauge multiplets, fermions, and scalars.

The decomposition of the $SO(10)$ representations under the relevant sub-
groups are presented in Tables 7, 8, 9, 10, and 11. In these tables,
multiplets highlighted in red, blue, and green contribute to the RGE running
of the gauge coupling in the energy scale in between $SO(10)$ and $G_{422}$,
$G_{422}$ and $G_{3221}$, and $G_{3221}$ and $G_{321}$, respectively. On the
other hand, multiplets highlighted in gray represent the Goldstone bosons.
Utilizing these information from Tables 7-11 and using Eqs. (42)-(43), in the
following, we compute the relevant beta coefficients required for gauge
coupling unification.

The beta coefficients of the SM sector are well-known, which are:

$\displaystyle\bigg{\\{}b^{\textrm{SM}}_{1Y},b^{\textrm{SM}}_{2L},b^{\textrm{SM}}_{3C}\bigg{\\}}=\left\\{\frac{41}{10},-\frac{19}{6},-7\right\\},$
(44) $\displaystyle
b^{\textrm{SM}}_{ij}=\left(\begin{array}[]{ccc}\frac{199}{50}&\frac{27}{10}&\frac{44}{5}\\\
\frac{9}{10}&\frac{35}{6}&12\\\ \frac{11}{10}&\frac{9}{2}&-26\\\
\end{array}\right).$ (48)

For our model, the beta coefficients relevant to the energy range between
left-right symmetry and Pati-Salam symmetry are:

$\displaystyle\bigg{\\{}b^{\textrm{LR}}_{2L},b^{\textrm{LR}}_{2R},b^{\textrm{LR}}_{3C},b^{\textrm{LR}}_{BL}\bigg{\\}}=\left\\{-\frac{17}{6},-\frac{13}{6},-7,\frac{11}{2}\right\\},$
(49) $\displaystyle
b^{\textrm{LR}}_{ij}=\left(\begin{array}[]{cccc}\frac{61}{6}&\frac{9}{2}&12&\frac{3}{2}\\\
\frac{9}{2}&\frac{173}{6}&12&\frac{27}{2}\\\
\frac{9}{2}&\frac{9}{2}&-26&\frac{1}{2}\\\
\frac{9}{2}&\frac{81}{2}&4&\frac{61}{2}\\\ \end{array}\right),$ (54)

and in the energy range between Pati-Salam and GUT symmetries are:

$\displaystyle\bigg{\\{}b^{\textrm{PS}}_{2L},b^{\textrm{PS}}_{2R},b^{\textrm{PS}}_{4C}\bigg{\\}}=\left\\{\frac{17}{2},\frac{17}{2},\frac{4}{3}\right\\},$
(55) $\displaystyle
b^{\textrm{PS}}_{ij}=\left(\begin{array}[]{ccc}\frac{515}{2}&\frac{93}{2}&\frac{1245}{2}\\\
\frac{93}{2}&\frac{515}{2}&\frac{1245}{2}\\\
\frac{249}{2}&\frac{249}{2}&\frac{3775}{6}\\\ \end{array}\right).$ (59)

Figure 9: Plot without incorporating threshold correction, i.e., $R=1$.
Current Super-K bound requires $M_{II}<4\times 10^{11}$ GeV, and future
Hyper-K bound would push this scale to below $M_{II}<10^{11}$ GeV.

At a symmetry breaking scale, $\hat{\mu}$, transitioning from the parent group
$G_{p}$ to the daughter gauge group $G_{d}$ happens. This requires a proper
matching condition between the gauge couplings $\alpha_{d}^{-1}$ of the
daughter gauge group and the gauge couplings $\alpha_{p}^{-1}$ of the parent
group, which reads

$\displaystyle\alpha_{d}^{-1}(\hat{\mu})-\frac{C_{2}(G_{d})}{12\pi}=\left(\alpha_{p}^{-1}(\hat{\mu})-\frac{C_{2}(G_{p})}{12\pi}\right)-\frac{\Lambda_{d}(\hat{\mu})}{12\pi},$
(60)

where one-loop threshold corrections after integrating out the superheavy
fields are given by [119, 120, 121, 122],

$\displaystyle\Lambda_{d}(\hat{\mu})=-21\;{\mathrm{Tr}}(t_{dV}^{2}\ln\frac{m_{V}}{\hat{\mu}})+2\;\eta\;{\mathrm{Tr}}(t_{dS}^{2}\ln\frac{m_{S}}{\hat{\mu}})+8\;\kappa\;{\mathrm{Tr}}(t_{dF}^{2}\ln\frac{m_{F}}{\hat{\mu}}).$
(61)

Here, $t_{dV}$, $t_{dS}$, and $t_{dF}$ are the generators for the
representations under $G_{d}$ of the superheavy vector, scalar, and fermion
fields, respectively, $m_{V}$, $m_{S}$, and $m_{F}$ are their respective
masses.

First, we consider a scenario without threshold corrections, meaning we set
$R=m_{S}/\hat{\mu}=1$; the result is presented in Fig. 9. As can be seen from
this figure, the current Super-K [109] bound demands the seesaw scale to be
smaller than $M_{II}<4\times 10^{11}$ GeV. On the other hand, future Hyper-K
[123] bound would push this scale to below $M_{II}<10^{11}$ GeV, leaving a
small window for which gauge coupling unification can be obtained. However,
the seesaw scale of this magnitude is inconsistent with the findings of the
PTAs and the fermion mass fit described above. This shows the need for
considering the threshold corrections.

When threshold corrections are turned on, we follow the common approach in
which the relevant superheavy gauge bosons at each symmetry breaking scale
have the masses identified with the corresponding breaking scale. On the other
hand, the scalars are assumed to be non-degenerate; hence, only scalars
contribute to the threshold correction. Masses of the scalars are varied
within the range $[1/R-R]$, with $R=2-5$. In this procedure, all the
superheavy degrees of freedom have the masses around the breaking scale within
the span $[1/R-R]$. We perform a $\chi^{2}$-analysis with a fixed seesaw scale
(corresponding to $M_{\mathrm{LR}}=2\times 10^{13}$ GeV and $4\times 10^{13}$
GeV), where the GUT scale is maximized by varying the associated threshold
corrections within a given range. The result is depicted in Fig. 10. As can be
seen from this figure, allowing a small threshold correction is enough for the
compatibility with gravitational wave observations, fermion masses and
mixings, as well as unification and proton decay constraints.

Figure 10:  Plot with threshold corrections for $R=[2-5]$. As can be seen from
this plot, a threshold correction with $R=2$ is enough to achieve high scale
gauge coupling unification which is also consistent with the current lower
bound on proton decay. Hence no additional fine-turning is required in this
model.

##  7 Conclusions

Grand unified theories are attractive candidates for physics beyond the
Standard Model. In this work, we have considered a minimal realistic grand
unified theory based on $SO(10)$ gauge group. The Yukawa sector of this theory
consists of a real $10_{H}$, a real $120_{H}$, and a complex $126_{H}$
dimensional representations, which is the most minimal option without
introducing any additional symmetries. The symmetry breaking sector that
employs a real $54_{H}$ and a real $45_{H}$ dimensional representations is
also taken to be the most economical that can give rise to stable cosmic
strings when the last step of the symmetry breaking takes place. Intriguingly,
a series of pulsar-timing arrays might have already observed GWs in the nHz
regime, hinting towards forming a cosmic string network in the very early
universe, which we identified with the seesaw scale generating tiny masses for
the active neutrinos. To be consistent with these observations, we find that
the seesaw scale is restricted to be within a narrow range
$v_{R}\in(3.4,4.1)\times 10^{13}$ GeV. More data, however, is required to
confirm if these signals are from a stochastic gravitational wave background.

Since within this framework, the charged fermion masses are inherently
interconnected with the neutrino masses, low energy data of fermion masses and
mixing angles allow this seesaw scale in the range $v_{R}\in(10^{13},10^{15})$
GeV. Contrastingly, if no threshold correction from the scalar spectrum is
considered, constraints from gauge coupling unification and proton decay
demand a much smaller seesaw scale $v_{R}\in(10^{10},10^{11})$ GeV. This
latter range for the seesaw scale is in direct conflict with observations of
pulsar-timing arrays as well as from low energy data. Our detailed study shows
that if the scalars are not assumed to be degenerate with the corresponding
symmetry breaking scales, then a small amount of threshold correction is
sufficient for compatibility with gravitational wave observations, fermion
mass fit, gauge coupling unification, and proton decay bounds. The proposed
minimal model is highly predictive and will be fully tested in a series of
experiments in the near future.

## Acknowledgments

S.S would like to thank Valerie Domcke for useful discussion.

##  Appendix 1 Best fit parameters without restricting the seesaw scale

The best fit parameters for the NO solution is as follows:

$\displaystyle r_{1}=4.99994\times 10^{-3}\,+6.64062\times 10^{-5}i,$ (62)
$\displaystyle r_{2}=-1.33796-0.257534i,$ (63) $\displaystyle\phi=-0.0504299,$
(64) $\displaystyle c_{R}=2.72334\times 10^{12},$ (65)

$S=\left(\begin{array}[]{ccc}7.83123\times 10^{-8}&0&0\\\ 0&0.237256&0\\\
0&0&83.9454\\\ \end{array}\right)\;\rm{GeV},$ (66)

$D={\mbox{$\displaystyle\left(\begin{array}[]{ccc}0.00051669\,+0.000140812i&0.00470802\,+0.00142105i&0.0293596\,+0.00308667i\\\
0.00470802\,+0.00142105i&0.0299984\,-0.0299305i&0.457075\,-0.132563i\\\
0.0293596\,+0.00308667i&0.457075\,-0.132563i&0.0118382\,-0.456846i\\\
\end{array}\right)$}}\;\rm{GeV},$ (67)

$A={\mbox{$\displaystyle\left(\begin{array}[]{ccc}0&0.00363066\,+0.00172004i&0.0220297\,+0.00302919i\\\
-0.00363066-0.00172004i&0&0.42873\,-0.0919318i\\\
-0.0220297-0.00302919i&-0.42873+0.0919318i&0\\\
\end{array}\right)$}}\;\rm{GeV}.$ (68)

For the IO case, the best fit parameters are given below,

$\displaystyle r_{1}=5.26833\times 10^{-3}\,-3.51051\times 10^{-4}i,$ (69)
$\displaystyle r_{2}=-1.42403-0.253046i,$ (70) $\displaystyle\phi=0.0569906,$
(71) $\displaystyle c_{R}=2.82679\times 10^{12},$ (72)

$S=\left(\begin{array}[]{ccc}4.62836\times 10^{-9}&0&0\\\ 0&0.227095&0\\\
0&0&83.8975\\\ \end{array}\right)\;\rm{GeV},$ (73)

$D={\mbox{$\displaystyle\left(\begin{array}[]{ccc}0.000544454\,+0.0000187733i&0.00466318\,+6.02921\times
10^{-6}i&0.0232066\,-0.0017959i\\\ 0.00466318\,+6.02921\times
10^{-6}i&0.0417626\,-0.0141477i&0.497413\,-0.00289039i\\\
0.0232066\,-0.0017959i&0.497413\,-0.00289039i&0.0700564\,-0.185603i\\\
\end{array}\right)$}}\;\rm{GeV},$ (74)

$A={\mbox{$\displaystyle\left(\begin{array}[]{ccc}0&0.00322916\,+0.000429876i&0.0162086\,-0.000403031i\\\
-0.00322916-0.000429876i&0&0.460026\,-0.0329099i\\\
-0.0162086+0.000403031i&-0.460026+0.0329099i&0\\\
\end{array}\right)$}}\;\rm{GeV}.$ (75)

##  Appendix 2 Best fit parameters for fixed seesaw scale

As described in the main text, in achieving this fit, we fixed $M_{3}=4\times
10^{13}$ GeV. The best fit parameters in this scenario are given by,

$\displaystyle r_{1}=4.79222\times 10^{-3}\,+5.85329\times 10^{-5}i,$ (76)
$\displaystyle r_{2}=-0.88177-1.51667i,$ (77) $\displaystyle\phi=0.682943,$
(78) $\displaystyle c_{R}=4.7\times 10^{11},$ (79)

$S=\left(\begin{array}[]{ccc}9.97395\times 10^{-7}&0&0\\\ 0&0.28853&0\\\
0&0&85.6185\\\ \end{array}\right)\;\rm{GeV},$ (80)

$D={\mbox{$\displaystyle\left(\begin{array}[]{ccc}-0.000178778+0.000261172i&0.00621503\,+0.000737064i&0.0124298\,-0.000291788i\\\
0.00621503\,+0.000737064i&-0.0135252+0.0153914i&0.337599\,+0.0444584i\\\
0.0124298\,-0.000291788i&0.337599\,+0.0444584i&0.0394304\,-0.832908i\\\
\end{array}\right)$}}\;\rm{GeV},$ (81)

$A={\mbox{$\displaystyle\left(\begin{array}[]{ccc}0&0.00297024\,+0.00195566i&0.0115613\,-0.0108631i\\\
-0.00297024-0.00195566i&0&0.235703\,-0.157565i\\\
-0.0115613+0.0108631i&-0.235703+0.157565i&0\\\
\end{array}\right)$}}\;\rm{GeV}.$ (82)

##  Appendix 3 Decomposition of representations

The decomposition of the $SO(10)$ representations under the relevant sub-
groups are presented in Tables 7, 8, 9, 10, and 11.

$SU(4)_{C}\times SU(2)_{L}\times SU(2)_{R}$ | $SU(3)_{C}\times SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$ | $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$
---|---|---
$\mathbf{(6,1,1)}$ | $\mathbf{(3,1,1,-2/3)}$ | $\mathbf{(3,1,-1/3)}$
| $\mathbf{(\overline{3},1,1,2/3)}$ | $\mathbf{(\overline{3},1,1/3)}$
$\mathbf{(1,2,2)}$ |  $\mathbf{(1,2,2,0)}$ | $\mathbf{(1,2,1/2)}$
|  | $\mathbf{(1,2,-1/2)}$

Table 7:  Decomposition of $\mathbf{10}$-Higgs, which is a real field.
Multiplets highlighted in red, blue, and green contribute to the RGE running
of the gauge coupling in the energy scale in between $SO(10)$ and $G_{422}$,
$G_{422}$ and $G_{3221}$, and $G_{3221}$ and $G_{321}$, respectively.
Multiplets highlighted in gray represent the Goldstone bosons.

$SU(4)_{C}\times SU(2)_{L}\times SU(2)_{R}$ | $SU(3)_{C}\times SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$ | $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$
---|---|---
$\mathbf{(1,1,3)}$ | $\mathbf{(1,1,3,0)}$ | $\mathbf{(1,1,1)}$
|  | $\mathbf{(1,1,0)}$
|  | $\mathbf{(1,1,-1)}$
$\mathbf{(1,3,1)}$ | $\mathbf{(1,3,1,0)}$ | $\mathbf{(1,3,0)}$
$\mathbf{(6,2,2)}$ | $\mathbf{(3,2,2,-2/3)}$ | $\mathbf{(3,2,1/6)}$
|  | $\mathbf{(3,2,-5/6)}$
| $\mathbf{(\overline{3},2,2,2/3)}$ | $\mathbf{(\overline{3},2,5/6)}$
|  | $\mathbf{(\overline{3},2,-1/6)}$
$\mathbf{(15,1,1)}$ | $\mathbf{(1,1,1,0)}$ | $\mathbf{(1,1,0)}$
|  $\mathbf{(3,1,1,4/3)}$ | $\mathbf{(3,1,2/3)}$
|  $\mathbf{(\overline{3},1,1,-4/3)}$ | $\mathbf{(\overline{3},1,-2/3)}$
| $\mathbf{(8,1,1,0)}$ | $\mathbf{(8,1,0)}$

Table 8:  Decomposition of $\mathbf{45}$-Higgs, which is a real field.
Multiplets highlighted in red, blue, and green contribute to the RGE running
of the gauge coupling in the energy scale in between $SO(10)$ and $G_{422}$,
$G_{422}$ and $G_{3221}$, and $G_{3221}$ and $G_{321}$, respectively.
Multiplets highlighted in gray represent the Goldstone bosons.

$SU(4)_{C}\times SU(2)_{L}\times SU(2)_{R}$ | $SU(3)_{C}\times SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$ | $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$
---|---|---
$\mathbf{(1,1,1)}$ | $\mathbf{(1,1,1,0)}$ | $\mathbf{(1,1,0)}$
$\mathbf{(1,3,3)}$ | $\mathbf{(1,3,3,0)}$ | $\mathbf{(1,3,1)}$
|  | $\mathbf{(1,3,0)}$
|  | $\mathbf{(1,3,-1)}$
$\mathbf{(20^{\prime},1,1)}$ | $\mathbf{(\overline{6},1,1,4/3)}$ | $\mathbf{(\overline{6},1,2/3)}$
| $\mathbf{(6,1,1,-4/3)}$ | $\mathbf{(6,1,-2/3)}$
| $\mathbf{(8,1,1,0)}$ | $\mathbf{(8,1,0)}$
$\mathbf{(6,2,2)}$ | $\mathbf{(3,2,2,-2/3)}$ | $\mathbf{(3,2,1/6)}$
|  | $\mathbf{(3,2,-5/6)}$
| $\mathbf{(\overline{3},2,2,2/3)}$ | $\mathbf{(\overline{3},2,5/6)}$
|  | $\mathbf{(3,2,-1/6)}$

Table 9:  Decomposition of $\mathbf{54}$-Higgs, which is a real field.
Multiplets highlighted in red, blue, and green contribute to the RGE running
of the gauge coupling in the energy scale in between $SO(10)$ and $G_{422}$,
$G_{422}$ and $G_{3221}$, and $G_{3221}$ and $G_{321}$, respectively.
Multiplets highlighted in gray represent the Goldstone bosons.

$SU(4)_{C}\times SU(2)_{L}\times SU(2)_{R}$ | $SU(3)_{C}\times SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$ | $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$
---|---|---
$\mathbf{(1,2,2)}$ | $\mathbf{(1,2,2,0)}$ | $\mathbf{(1,2,1/2)}$
|  | $\mathbf{(1,2,-1/2)}$
$\mathbf{(10,1,1)}$ | $\mathbf{(1,1,1,-2)}$ | $\mathbf{(1,1,-1)}$
| $\mathbf{(3,1,1,-2/3)}$ | $\mathbf{(3,1,-1/3)}$
| $\mathbf{(6,1,1,2/3)}$ | $\mathbf{(6,1,1/3)}$
$\mathbf{(\overline{1}0,1,1)}$ | $\mathbf{(1,1,1,2)}$ | $\mathbf{(1,1,1)}$
| $\mathbf{(\overline{3},1,1,2/3)}$ | $\mathbf{(\overline{3},1,1/3)}$
| $\mathbf{(\overline{6},1,1,-2/3)}$ | $\mathbf{(\overline{6},1,-1/3)}$
$\mathbf{(6,3,1)}$ | $\mathbf{(3,3,1,-2/3)}$ | $\mathbf{(3,3,-1/3)}$
| $\mathbf{(\overline{3},3,1,2/3)}$ | $\mathbf{(\overline{3},3,1/3)}$
$\mathbf{(6,1,3)}$ | $\mathbf{(3,1,3,-2/3)}$ | $\mathbf{(3,1,2/3)}$
|  | $\mathbf{(3,1,-1/3)}$
|  | $\mathbf{(3,1,-4/3)}$
| $\mathbf{(\overline{3},1,3,2/3)}$ | $\mathbf{(\overline{3},1,4/3)}$
|  | $\mathbf{(\overline{3},1,1/3)}$
|  | $\mathbf{(\overline{3},1,-2/3)}$
$\mathbf{(15,2,2)}$ | $\mathbf{(1,2,2,0)}$ | $\mathbf{(1,2,1/2)}$
|  | $\mathbf{(1,2,-1/2)}$
| $\mathbf{(3,2,2,4/3)}$ | $\mathbf{(3,2,7/6)}$
|  | $\mathbf{(\overline{3},2,1/6)}$
| $\mathbf{(\overline{3},2,2,-4/3)}$ | $\mathbf{(\overline{3},2,-7/6)}$
|  | $\mathbf{(3,2,-1/6)}$
| $\mathbf{(8,2,2,0)}$ | $\mathbf{(8,2,1/2)}$
|  | $\mathbf{(8,2,-1/2)}$

Table 10:  Decomposition of $\mathbf{120}$-Higgs, which is a real field. For
the simplicity of the analysis, we assume the entire $120_{H}$ lives near to
the GUT scale and does not contribute to the RGEs. Putting the bi-doublet at
the intermediate scale, however, can be straightforwardly incorporated.

$SU(4)_{C}\times SU(2)_{L}\times SU(2)_{R}$ | $SU(3)_{C}\times SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$ | $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}$
---|---|---
$\mathbf{(6,1,1)}$ | $\mathbf{(3,1,1,-2/3)}$ | $\mathbf{(3,1,-1/3)}$
| $\mathbf{(\overline{3},1,1,2/3)}$ | $\mathbf{(\overline{3},1,1/3)}$
$\mathbf{(\overline{10},3,1)}$ | $\mathbf{(1,3,1,2)}$ | $\mathbf{(1,3,1)}$
| $\mathbf{(\overline{3},3,1,2/3)}$ | $\mathbf{(\overline{3},3,1/3)}$
| $\mathbf{(\overline{6},3,1,-2/3)}$ | $\mathbf{(\overline{6},3,-1/3)}$
$\mathbf{(10,1,3)}$ |  $\mathbf{(1,1,3,-2)}$ |  $\mathbf{(1,1,0)}$
|  |  $\mathbf{(1,1,-1)}$
|  | $\mathbf{(1,1,-2)}$
| $\mathbf{(3,1,3,-2/3)}$ | $\mathbf{(3,1,2/3)}$
|  | $\mathbf{(3,1,-1/3)}$
|  | $\mathbf{(3,1,-4/3)}$
| $\mathbf{(6,1,3,2/3)}$ | $\mathbf{(6,1,4/3)}$
|  | $\mathbf{(6,1,1/3)}$
|  | $\mathbf{(6,1,-2/3)}$
$\mathbf{(15,2,2)}$ |  $\mathbf{(1,2,2,0)}$ |  $\mathbf{(1,2,1/2)}$
|  |  $\mathbf{(1,2,-1/2)}$
| $\mathbf{(3,2,2,4/3)}$ | $\mathbf{(3,2,7/6)}$
|  | $\mathbf{(\overline{3},2,1/6)}$
| $\mathbf{(\overline{3},2,2,-4/3)}$ | $\mathbf{(\overline{3},2,-7/6)}$
|  | $\mathbf{(3,2,-1/6)}$
| $\mathbf{(8,2,2,0)}$ | $\mathbf{(8,2,1/2)}$
|  | $\mathbf{(8,2,-1/2)}$

Table 11:  Decomposition of $\mathbf{\overline{126}}$-Higgs, which is a
complex field. Multiplets highlighted in red, blue, and green contribute to
the RGE running of the gauge coupling in the energy scale in between $SO(10)$
and $G_{422}$, $G_{422}$ and $G_{3221}$, and $G_{3221}$ and $G_{321}$,
respectively. Multiplets highlighted in gray represent the Goldstone bosons.

## References

  * [1] J. C. Pati and A. Salam, _Is Baryon Number Conserved?_ , _Phys. Rev. Lett._ 31 (1973) 661.
  * [2] J. C. Pati and A. Salam, _Lepton Number as the Fourth Color_ , _Phys. Rev. D_ 10 (1974) 275.
  * [3] H. Georgi and S. Glashow, _Unity of All Elementary Particle Forces_ , _Phys. Rev. Lett._ 32 (1974) 438.
  * [4] H. Georgi, H. R. Quinn and S. Weinberg, _Hierarchy of Interactions in Unified Gauge Theories_ , _Phys. Rev. Lett._ 33 (1974) 451.
  * [5] H. Georgi, _The State of the Art—Gauge Theories_ , _AIP Conf. Proc._ 23 (1975) 575.
  * [6] H. Fritzsch and P. Minkowski, _Unified Interactions of Leptons and Hadrons_ , _Annals Phys._ 93 (1975) 193.
  * [7] P. Langacker, _Grand Unified Theories and Proton Decay_ , _Phys. Rept._ 72 (1981) 185.
  * [8] A. Vilenkin, _Cosmic Strings and Domain Walls_ , _Phys. Rept._ 121 (1985) 263.
  * [9] R. R. Caldwell and B. Allen, _Cosmological constraints on cosmic string gravitational radiation_ , _Phys. Rev._ D45 (1992) 3447.
  * [10] M. B. Hindmarsh and T. W. B. Kibble, _Cosmic strings_ , _Rept. Prog. Phys._ 58 (1995) 477 [hep-ph/9411342].
  * [11] LIGO Scientific, Virgo collaboration, _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116 (2016) 061102 [1602.03837].
  * [12] LIGO Scientific, Virgo collaboration, _GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence_ , _Phys. Rev. Lett._ 116 (2016) 241103 [1606.04855].
  * [13] LIGO Scientific, VIRGO collaboration, _GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2_ , _Phys. Rev. Lett._ 118 (2017) 221101 [1706.01812].
  * [14] LIGO Scientific, Virgo collaboration, _GW170608: Observation of a 19-solar-mass Binary Black Hole Coalescence_ , _Astrophys. J. Lett._ 851 (2017) L35 [1711.05578].
  * [15] LIGO Scientific, Virgo collaboration, _GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence_ , _Phys. Rev. Lett._ 119 (2017) 141101 [1709.09660].
  * [16] LIGO Scientific, Virgo collaboration, _GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral_ , _Phys. Rev. Lett._ 119 (2017) 161101 [1710.05832].
  * [17] LIGO Scientific, Virgo collaboration, _GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object_ , _Astrophys. J. Lett._ 896 (2020) L44 [2006.12611].
  * [18] LIGO Scientific, Virgo collaboration, _GW190521: A Binary Black Hole Merger with a Total Mass of $150M_{\odot}$_, _Phys. Rev. Lett._ 125 (2020) 101102 [2009.01075].
  * [19] NANOGrav collaboration, _The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background_ , _Astrophys. J. Lett._ 905 (2020) L34 [2009.04496].
  * [20] B. Goncharov et al., _On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 917 (2021) L19 [2107.12112].
  * [21] S. Chen et al., _Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: inferences in the stochastic gravitational-wave background search_ , _Mon. Not. Roy. Astron. Soc._ 508 (2021) 4970 [2110.13184].
  * [22] J. Antoniadis et al., _The International Pulsar Timing Array second data release: Search for an isotropic gravitational wave background_ , _Mon. Not. Roy. Astron. Soc._ 510 (2022) 4873 [2201.03980].
  * [23] K. Babu, B. Bajc and S. Saad, _Yukawa Sector of Minimal SO(10) Unification_ , _JHEP_ 02 (2017) 136 [1612.04329].
  * [24] K. Babu and R. Mohapatra, _Predictive neutrino spectrum in minimal SO(10) grand unification_ , _Phys. Rev. Lett._ 70 (1993) 2845 [hep-ph/9209215].
  * [25] B. Bajc, G. Senjanovic and F. Vissani, _How neutrino and charged fermion masses are connected within minimal supersymmetric SO(10)_ , _PoS_ HEP2001 (2001) 198 [hep-ph/0110310].
  * [26] B. Bajc, G. Senjanovic and F. Vissani, _b - tau unification and large atmospheric mixing: A Case for noncanonical seesaw_ , _Phys. Rev. Lett._ 90 (2003) 051802 [hep-ph/0210207].
  * [27] T. Fukuyama and N. Okada, _Neutrino oscillation data versus minimal supersymmetric SO(10) model_ , _JHEP_ 11 (2002) 011 [hep-ph/0205066].
  * [28] H. S. Goh, R. N. Mohapatra and S.-P. Ng, _Minimal SUSY SO(10), b tau unification and large neutrino mixings_ , _Phys. Lett._ B570 (2003) 215 [hep-ph/0303055].
  * [29] H. S. Goh, R. N. Mohapatra and S.-P. Ng, _Minimal SUSY SO(10) model and predictions for neutrino mixings and leptonic CP violation_ , _Phys. Rev._ D68 (2003) 115008 [hep-ph/0308197].
  * [30] S. Bertolini, M. Frigerio and M. Malinsky, _Fermion masses in SUSY SO(10) with type II seesaw: A Non-minimal predictive scenario_ , _Phys. Rev._ D70 (2004) 095002 [hep-ph/0406117].
  * [31] S. Bertolini and M. Malinsky, _On CP violation in minimal renormalizable SUSY SO(10) and beyond_ , _Phys. Rev. D_ 72 (2005) 055021 [hep-ph/0504241].
  * [32] K. S. Babu and C. Macesanu, _Neutrino masses and mixings in a minimal SO(10) model_ , _Phys. Rev._ D72 (2005) 115003 [hep-ph/0505200].
  * [33] S. Bertolini, T. Schwetz and M. Malinsky, _Fermion masses and mixings in SO(10) models and the neutrino challenge to SUSY GUTs_ , _Phys. Rev._ D73 (2006) 115012 [hep-ph/0605006].
  * [34] B. Bajc, I. Dorsner and M. Nemevsek, _Minimal SO(10) splits supersymmetry_ , _JHEP_ 11 (2008) 007 [0809.1069].
  * [35] A. S. Joshipura and K. M. Patel, _Fermion Masses in SO(10) Models_ , _Phys. Rev._ D83 (2011) 095002 [1102.5148].
  * [36] G. Altarelli and D. Meloni, _A non supersymmetric SO(10) grand unified model for all the physics below $M_{GUT}$_, _JHEP_ 08 (2013) 021 [1305.1001].
  * [37] A. Dueck and W. Rodejohann, _Fits to SO(10) Grand Unified Models_ , _JHEP_ 09 (2013) 024 [1306.4468].
  * [38] T. Fukuyama, K. Ichikawa and Y. Mimura, _Revisiting fermion mass and mixing fits in the minimal SUSY $SO(10)$ GUT_, _Phys. Rev._ D94 (2016) 075018 [1508.07078].
  * [39] K. Babu, B. Bajc and S. Saad, _New Class of SO(10) Models for Flavor_ , _Phys. Rev. D_ 94 (2016) 015030 [1605.05116].
  * [40] S. Saad, _SO(10) models for flavor with vector-like fermions_ , _AIP Conf. Proc._ 1900 (2017) 020010.
  * [41] K. S. Babu, B. Bajc and S. Saad, _Resurrecting Minimal Yukawa Sector of SUSY SO(10)_ , _JHEP_ 10 (2018) 135 [1805.10631].
  * [42] K. Babu, T. Fukuyama, S. Khan and S. Saad, _Peccei-Quinn Symmetry and Nucleon Decay in Renormalizable SUSY $SO$(10)_, _JHEP_ 06 (2019) 045 [1812.11695].
  * [43] T. Ohlsson and M. Pernow, _Fits to Non-Supersymmetric SO(10) Models with Type I and II Seesaw Mechanisms Using Renormalization Group Evolution_ , _JHEP_ 06 (2019) 085 [1903.08241].
  * [44] K. S. Babu and S. Saad, _Flavor Hierarchies from Clockwork in SO(10) GUT_ , _Phys. Rev. D_ 103 (2021) 015009 [2007.16085].
  * [45] V. S. Mummidi and K. M. Patel, _Leptogenesis and fermion mass fit in a renormalizable SO(10) model_ , _JHEP_ 12 (2021) 042 [2109.04050].
  * [46] P. Minkowski, _$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?_, _Phys. Lett. B_ 67 (1977) 421.
  * [47] T. Yanagida, _Horizontal gauge symmetry and masses of neutrinos_ , _Conf. Proc. C_ 7902131 (1979) 95.
  * [48] S. Glashow, _The Future of Elementary Particle Physics_ , _NATO Sci. Ser. B_ 61 (1980) 687.
  * [49] M. Gell-Mann, P. Ramond and R. Slansky, _Complex Spinors and Unified Theories_ , _Conf. Proc. C_ 790927 (1979) 315 [1306.4669].
  * [50] R. N. Mohapatra and G. Senjanovic, _Neutrino Mass and Spontaneous Parity Nonconservation_ , _Phys. Rev. Lett._ 44 (1980) 912.
  * [51] S. Antusch, J. Kersten, M. Lindner, M. Ratz and M. A. Schmidt, _Running neutrino mass parameters in see-saw scenarios_ , _JHEP_ 03 (2005) 024 [hep-ph/0501272].
  * [52] _Nufit webpage, available online:http://www.nu-fit.org (october 2021 data)_, .
  * [53] S. Antusch and V. Maurer, _Running quark and lepton parameters at various scales_ , _JHEP_ 11 (2013) 115 [1306.6879].
  * [54] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, _The fate of hints: updated global analysis of three-flavor neutrino oscillations_ , _JHEP_ 09 (2020) 178 [2007.14792].
  * [55] T. W. B. Kibble, _Topology of Cosmic Domains and Strings_ , _J. Phys. A_ 9 (1976) 1387.
  * [56] A. Vilenkin and E. P. S. Shellard, _Cosmic Strings and Other Topological Defects_. Cambridge University Press, 7, 2000.
  * [57] A. H. Guth, _The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems_ , _Phys. Rev. D_ 23 (1981) 347.
  * [58] A. Albrecht and P. J. Steinhardt, _Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking_ , _Phys. Rev. Lett._ 48 (1982) 1220.
  * [59] A. D. Linde, _A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems_ , _Phys. Lett. B_ 108 (1982) 389.
  * [60] A. D. Linde, _Chaotic Inflation_ , _Phys. Lett. B_ 129 (1983) 177.
  * [61] E. W. Kolb and M. S. Turner, _The Early Universe_ , vol. 69. 1990.
  * [62] D. Chang, R. Mohapatra and M. Parida, _Decoupling Parity and SU(2)-R Breaking Scales: A New Approach to Left-Right Symmetric Models_ , _Phys. Rev. Lett._ 52 (1984) 1072.
  * [63] T. W. B. Kibble, G. Lazarides and Q. Shafi, _Strings in SO(10)_ , _Phys. Lett. B_ 113 (1982) 237.
  * [64] R. Jeannerot, J. Rocher and M. Sakellariadou, _How generic is cosmic string formation in SUSY GUTs_ , _Phys. Rev. D_ 68 (2003) 103514 [hep-ph/0308134].
  * [65] W. Buchmuller, V. Domcke, H. Murayama and K. Schmitz, _Probing the scale of grand unification with gravitational waves_ , _Phys. Lett. B_ 809 (2020) 135764 [1912.03695].
  * [66] S. F. King, S. Pascoli, J. Turner and Y.-L. Zhou, _Gravitational Waves and Proton Decay: Complementary Windows into Grand Unified Theories_ , _Phys. Rev. Lett._ 126 (2021) 021802 [2005.13549].
  * [67] S. F. King, S. Pascoli, J. Turner and Y.-L. Zhou, _Confronting SO(10) GUTs with proton decay and gravitational waves_ , _JHEP_ 10 (2021) 225 [2106.15634].
  * [68] M. A. Masoud, M. U. Rehman and Q. Shafi, _Sneutrino tribrid inflation, metastable cosmic strings and gravitational waves_ , _JCAP_ 11 (2021) 022 [2107.09689].
  * [69] G. Lazarides, R. Maji and Q. Shafi, _Cosmic strings, inflation, and gravity waves_ , _Phys. Rev. D_ 104 (2021) 095004 [2104.02016].
  * [70] G. Lazarides, R. Maji and Q. Shafi, _Gravitational waves from quasi-stable strings_ , _JCAP_ 08 (2022) 042 [2203.11204].
  * [71] R. Maji and Q. Shafi, _Monopoles, Strings and Gravitational Waves in Non-minimal Inflation_ , 2208.08137.
  * [72] B. Fu, S. F. King, L. Marsili, S. Pascoli, J. Turner and Y.-L. Zhou, _A predictive and testable unified theory of fermion masses, mixing and leptogenesis_ , _JHEP_ 11 (2022) 072 [2209.00021].
  * [73] G. Lazarides, R. Maji, R. Roshan and Q. Shafi, _A Predictive $SO(10)$ Model_, 2210.03710.
  * [74] T. Vachaspati and A. Vilenkin, _Gravitational Radiation from Cosmic Strings_ , _Phys. Rev._ D31 (1985) 3052.
  * [75] P. Auclair et al., _Probing the gravitational wave background from cosmic strings with LISA_ , _JCAP_ 04 (2020) 034 [1909.00819].
  * [76] M. Hindmarsh, _Signals of Inflationary Models with Cosmic Strings_ , _Prog. Theor. Phys. Suppl._ 190 (2011) 197 [1106.0391].
  * [77] A. Vilenkin, _Gravitational radiation from cosmic strings_ , _Phys. Lett._ 107B (1981) 47.
  * [78] N. Turok, _Grand Unified Strings and Galaxy Formation_ , _Nucl. Phys. B_ 242 (1984) 520.
  * [79] J. M. Quashnock and D. N. Spergel, _Gravitational Selfinteractions of Cosmic Strings_ , _Phys. Rev. D_ 42 (1990) 2505.
  * [80] J. J. Blanco-Pillado, K. D. Olum and B. Shlaer, _The number of cosmic string loops_ , _Phys. Rev._ D89 (2014) 023512 [1309.6637].
  * [81] J. J. Blanco-Pillado and K. D. Olum, _Stochastic gravitational wave background from smoothed cosmic string loops_ , _Phys. Rev._ D96 (2017) 104046 [1709.02693].
  * [82] Y. Cui, M. Lewicki, D. E. Morrissey and J. D. Wells, _Probing the pre-BBN universe with gravitational waves from cosmic strings_ , _JHEP_ 01 (2019) 081 [1808.08968].
  * [83] Planck collaboration, _Planck 2015 results. XIII. Cosmological parameters_ , _Astron. Astrophys._ 594 (2016) A13 [1502.01589].
  * [84] T. Charnock, A. Avgoustidis, E. J. Copeland and A. Moss, _CMB constraints on cosmic strings and superstrings_ , _Phys. Rev. D_ 93 (2016) 123503 [1603.01275].
  * [85] LIGO Scientific, Virgo, KAGRA collaboration, _Constraints on Cosmic Strings Using Data from the Third Advanced LIGO–Virgo Observing Run_ , _Phys. Rev. Lett._ 126 (2021) 241102 [2101.12248].
  * [86] L. Lentati et al., _European Pulsar Timing Array Limits On An Isotropic Stochastic Gravitational-Wave Background_ , _Mon. Not. Roy. Astron. Soc._ 453 (2015) 2576 [1504.03692].
  * [87] NANOGRAV collaboration, _The NANOGrav 11-year Data Set: Pulsar-timing Constraints On The Stochastic Gravitational-wave Background_ , _Astrophys. J._ 859 (2018) 47 [1801.02617].
  * [88] J. Ellis and M. Lewicki, _Cosmic String Interpretation of NANOGrav Pulsar Timing Data_ , _Phys. Rev. Lett._ 126 (2021) 041304 [2009.06555].
  * [89] B. Sathyaprakash et al., _Scientific Objectives of Einstein Telescope_ , _Class. Quant. Grav._ 29 (2012) 124013 [1206.0331].
  * [90] LIGO Scientific collaboration, _Exploring the Sensitivity of Next Generation Gravitational Wave Detectors_ , _Class. Quant. Grav._ 34 (2017) 044001 [1607.08697].
  * [91] LISA collaboration, _Laser Interferometer Space Antenna_ , 1702.00786.
  * [92] W.-H. Ruan, Z.-K. Guo, R.-G. Cai and Y.-Z. Zhang, _Taiji program: Gravitational-wave sources_ , _Int. J. Mod. Phys. A_ 35 (2020) 2050075 [1807.09495].
  * [93] TianQin collaboration, _TianQin: a space-borne gravitational wave detector_ , _Class. Quant. Grav._ 33 (2016) 035010 [1512.02076].
  * [94] V. Corbin and N. J. Cornish, _Detecting the cosmic gravitational wave background with the big bang observer_ , _Class. Quant. Grav._ 23 (2006) 2435 [gr-qc/0512039].
  * [95] N. Seto, S. Kawamura and T. Nakamura, _Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space_ , _Phys. Rev. Lett._ 87 (2001) 221103 [astro-ph/0108011].
  * [96] MAGIS collaboration, _Mid-band gravitational wave detection with precision atomic sensors_ , 1711.02225.
  * [97] AEDGE collaboration, _AEDGE: Atomic Experiment for Dark Matter and Gravity Exploration in Space_ , _EPJ Quant. Technol._ 7 (2020) 6 [1908.00802].
  * [98] L. Badurina et al., _AION: An Atom Interferometer Observatory and Network_ , _JCAP_ 05 (2020) 011 [1911.11755].
  * [99] G. Janssen et al., _Gravitational wave astronomy with the SKA_ , _PoS_ AASKA14 (2015) 037 [1501.00127].
  * [100] Gaia collaboration, _Gaia Data Release 2: Summary of the contents and survey properties_ , _Astron. Astrophys._ 616 (2018) A1 [1804.09365].
  * [101] Theia collaboration, _Theia: Faint objects in motion or the new astrometry frontier_ , 1707.01348.
  * [102] C. W. Fabjan and H. Schopper, eds., _Particle Physics Reference Library: Volume 2: Detectors for Particles and Radiation_. Springer Nature, Cham, 2020, 10.1007/978-3-030-35318-6.
  * [103] Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * [104] KamLAND-Zen collaboration, _Search for Majorana Neutrinos near the Inverted Mass Hierarchy Region with KamLAND-Zen_ , _Phys. Rev. Lett._ 117 (2016) 082503 [1605.02889].
  * [105] GERDA collaboration, _Modeling of GERDA Phase II data_ , _JHEP_ 03 (2020) 139 [1909.02522].
  * [106] nEXO collaboration, _nEXO: neutrinoless double beta decay search beyond 10 28 year half-life sensitivity_, _J. Phys. G_ 49 (2022) 015104 [2106.16243].
  * [107] P. S. B. Dev et al., _Searches for Baryon Number Violation in Neutrino Experiments: A White Paper_ , 2203.08771.
  * [108] M. Machacek, _The Decay Modes of the Proton_ , _Nucl. Phys. B_ 159 (1979) 37.
  * [109] Super-Kamiokande collaboration, _Search for proton decay via $p\to e^{+}\pi^{0}$ and $p\to\mu^{+}\pi^{0}$ with an enlarged fiducial volume in Super-Kamiokande I-IV_, _Phys. Rev. D_ 102 (2020) 112011 [2010.16098].
  * [110] P. Fileviez Perez, _Fermion mixings versus d = 6 proton decay_ , _Phys. Lett._ B595 (2004) 476 [hep-ph/0403286].
  * [111] P. Nath and P. Fileviez Perez, _Proton stability in grand unified theories, in strings and in branes_ , _Phys. Rept._ 441 (2007) 191 [hep-ph/0601023].
  * [112] Y. Aoki, T. Izubuchi, E. Shintani and A. Soni, _Improved lattice computation of proton decay matrix elements_ , _Phys. Rev. D_ 96 (2017) 014506 [1705.01338].
  * [113] T. Nihei and J. Arafune, _The Two loop long range effect on the proton decay effective Lagrangian_ , _Prog. Theor. Phys._ 93 (1995) 665 [hep-ph/9412325].
  * [114] A. J. Buras, J. R. Ellis, M. K. Gaillard and D. V. Nanopoulos, _Aspects of the Grand Unification of Strong, Weak and Electromagnetic Interactions_ , _Nucl. Phys. B_ 135 (1978) 66.
  * [115] J. T. Goldman and D. A. Ross, _How Accurately Can We Estimate the Proton Lifetime in an SU(5) Grand Unified Model?_ , _Nucl. Phys. B_ 171 (1980) 273.
  * [116] W. E. Caswell, J. Milutinovic and G. Senjanovic, _Predictions of Left-right Symmetric Grand Unified Theories_ , _Phys. Rev. D_ 26 (1982) 161.
  * [117] L. E. Ibanez and C. Munoz, _Enhancement Factors for Supersymmetric Proton Decay in the Wess-Zumino Gauge_ , _Nucl. Phys. B_ 245 (1984) 425.
  * [118] M. E. Machacek and M. T. Vaughn, _Two Loop Renormalization Group Equations in a General Quantum Field Theory. 1. Wave Function Renormalization_ , _Nucl. Phys. B_ 222 (1983) 83.
  * [119] S. Weinberg, _Effective Gauge Theories_ , _Phys. Lett. B_ 91 (1980) 51.
  * [120] L. J. Hall, _Grand Unification of Effective Gauge Theories_ , _Nucl. Phys. B_ 178 (1981) 75.
  * [121] S. Bertolini, L. Di Luzio and M. Malinsky, _Intermediate mass scales in the non-supersymmetric SO(10) grand unification: A Reappraisal_ , _Phys. Rev. D_ 80 (2009) 015013 [0903.4049].
  * [122] S. Bertolini, L. Di Luzio and M. Malinsky, _Light color octet scalars in the minimal SO(10) grand unification_ , _Phys. Rev. D_ 87 (2013) 085020 [1302.3401].
  * [123] Hyper-Kamiokande collaboration, _Hyper-Kamiokande Design Report_ , 1805.04163.