Tracking Down the Route to the SM with Inflation and Gravitational Waves

We explore supersymmetric SO(10) models predicting observable proton decay and various topological defects which produce different shapes and strengths of gravitational wave backgrounds depending on the scales of intermediate symmetry breaking and inflation as well. We compare these to their non-supersymmetric counterparts. By identifying the scales at which gravitational wave signals appear, we would be able to track down a particular breaking chain and discern if it has a supersymmetric origin or not. It would also be useful to observe gravitational waves from more than one source among all possible topological defects and first order phase transitions for a realistic breaking chain. For these purposes, we work out specific examples in which the grand unification and relevant intermediate scales are calculable explicitly. It turns out that examples with gravitational waves from different sources are quite difficult to obtain, and the predicted gravitational wave profiles from domain walls and first order phase transitions obtained in some examples will require detectors in the kHz to MHz region.


Introduction
The first ever direct detection of a gravitational wave (GW) signal in 2015 [1] opened the era of gravitational wave astronomy and revived interests in cosmological GW signals. The study of GW signals from cosmological origin, either topological defects or phase order first transitions (FOPT), is a long standing subject, but only recently the current and future experiments have the potential to detect these kind of signals. During the next years, KAGRA [2] and LIGO-India will be integrated in the LIGO [3,4,5,6] and Virgo [7] interferometer network (LVC). Thanks to the improved sensitivity, a huge number of GW events will be detected. The detection could be so frequent that the associated signals would be too entangled to be identify individually. These events will then appear as a Stochastic Gravitational Wave Background (SGWB), a signal continuous in time, loud in a very broad frequency band, and coming from the whole sky dome. In fact, these are also the characteristics of GW signals created at the earliest moments of the universe. Such signals, called SGWB of Cosmological Origin (SGWBoCO), exist but we do not know their strength and the frequencies at which we could detect them. It is then crucial to design SGWB searches suitable for scenarios where the SGWBoCO and SGWB of astrophysical origin (SGWBoAO) are both present at some level. This kind of searches will need to be implemented by all current and projected experiments, such as LISA 1 TianQin [9], Taiji [10], LISA [11], Einstein Telescope [12,13], Cosmic Explorer [14], BBO [15] and DECIGO [16], as well as Pulsar Timing Array (PTA) experiments SKA [17], EPTA [18], PPTA [19], IPTA [20]. In fact, in September 2020 the NANOGrav Collaboration announced their PTA analysis [21] showing that we are likely on the verge of a SGWB detection. In particular, it has hinted at the detection of a cosmic string in the range Gµ ∈ [2 × 10 −11 , 3 × 10 −10 ] [22]. A firm SGWB observation may require only some extra years of measurements, and it is encouraging to have the recent positive development of the PPTA analysis [23] and the collaborative effort of IPTA [24]. If confirmed, the NANOGrav SGWB signal in the range Gµ ∈ [2 × 10 −11 , 3 × 10 −10 ] would be compatible with the SGWBoCO produced by (Nambu Goto) cosmic strings coming from a phase transition that occurred between 10 −22 and 10 −18 seconds after the Big Bang.
On the theoretical side, model building within the context of Grand Unified Theories (GUT) needs to provide a compelling science case for scenarios where the cosmological signal can be clearly identified. In this respect, there has been a large number of recent studies reviving the study of GW coming from topological defects 2 . These include non-supersymmetric GUT theories in connection to proton decay and unification [31,32,33] and different interconnections to B − L [34], neutrino masses and leptogenesis [35], as well as in connection to addressing the NANOGrav hint for cosmic strings [36,37] and combined effects with FOPT and GW [38].
In this paper, we address the question if different breaking chains of the supersymmetric SO(10) GUT predicting a certain combination of topological defects and first order phase transition (FOPT) can be identified through the detection of SGWB. This task of course can be done taking into account low energy phenomenological constraints (proton decay, masses, flavor effects, etc.) and assumptions about how inflation should take place.
The classification of breaking routes of SO(10) down to Standard Model (SM), or Minimal Supersymmetric SM (MSSM), producing topological defects, started long time ago [27]. As it is well-known, non-supersymmetric models can yield gauge coupling unification with different intermediate scales [39] and therefore it is appealing to study the GW panorama in this context. Making a catalogue of surviving breaking chains coming from SO(10) models as in [32, 33] 3 is considerably more challenging in supersymmetry due to the fact that not only the number of constraints and parameters is considerably bigger but also to the assumptions for the way supersymmetry is broken and the way it is UV completed. For these reasons we present in this work only a few typical examples, rather than a comprehensive catalogue.
Supersymmetric models can be compared with their non-supersymmetric counterparts because they differ in their breaking scales and proton decay channels and this can be contrasted in the future with both GW and proton decay experiments. We provide a study for this comparison by considering examples based on the SO(10) breaking routes containing the SU(5) and SU(3) C × SU(2) L × SU(2) R × U(1) B−L sub-groups. These models are representatives of different subgroups and low energy phenomenology and therefore can shed light on discerning routes down to the SM. We envision plots for frequency vs. density showing all possible combination of topological defects and FOPT for a particular breaking chain. Unfortunately, this task is quite model-dependent and will require future developments of simulations, in particular, for hybrid topological defects since the relative tensions of the topological defects is crucial to the evolution of the networks. Nevertheless, it is worth to establish possible breaking routes where hybrid topological defects can leave a distinct GW imprint. Recently, these observations have been also revived in [42] with a similar motivation to ours.
The paper is organized as follows. In §2 we put into context the appearance of topological defects according to the different breaking routes down to the SM. In §3 we mention how do we obtain the scale of breaking and how do we compute the proton decay ratios. In §4 we analyze our model examples and compare to the non-supersymmetric counterparts. In §5 we conclude. For the sake of completeness, we provide the information that we use to generate the GW signals in the appendices.

GUT breaking and topological defects
In this section, we make a brief summary of topological defects appearing in a GUT breaking chain which we use and refer to. Then, we introduce typical SO(10) GUT breaking patterns for which we study various phenomenological features and GW signals from relevant topological defects and possible phase transitions.

Theory of topological defects in breaking chains
A schematic breaking chain of a GUT model down to the SM can be depicted as where the letters p i , q i and r i represent topological defects. The conditions for the formation, evolution and stability of topological defects has been thoroughly studied 4 . The k−homotopy group, π k (G/H), of the vacuum manifold M = G/H determines the appearance of topological defects since π k (G/H) classifies the distinct topological spaces of G/H that appear after the breaking of G → H. They correspond to k = 3 for textures, k = 2 for monopoles, k = 1 for cosmic strings (CS) and k = 0 for domain walls (DW), although formally π 0 is not a group: it represents merely the number of connected components of the manifold. If the homotopy groups are trivial, that is, π k (G/H) = I, then there are no associated defects. The general features of topological defects from the breaking sequence of Eq. (1) can be summarized as follows.
corrections to achieve gauge coupling unification and acceptable proton decay rates, and also to consider or not multiple parameters controlling topological defects [40,41].
1. The topological defects will be stable up to the H n group if the k-homotopy group of the manifold G GUT /H n is not trivial, that is if π k (G GUT /H n ) = I; and unstable up to the H n group if the k-homotopy group of the manifold G GUT /H n is trivial, that is if π k (G GUT /H n ) = I. The same applies for the manifold G GUT / (SU (3) C × U Y (1)) which will determine stable or unstable topological defects all the way down to the EW scale.
2. There can appear metastable defects that decay quantum mechanically with a decay rate ∝ e −A/ , where A is the "tunneling action" of the particular defect [43]. For the strings decaying via a pair of monopole/anti-monopole the decay rate, which is specifically the tunneling rate per unit string length, ℓ, is given by where µ is the CS tension and m is the scale at which the monopoles are begin created. This means that the creation of strings should not take place much below the formation of monopoles as only decays for values of k = O(1) render a non negligible value of Γ, since the exponential in Eq. (2) decays rapidly. For values much greater than 1, the decaying cosmic strings become indistinguishable from stable strings. In certain cases, the probability per unit area of a domain wall decaying by the nucleation of strings, goes like [44,43] dP dA where P is the probability of nucleation and A is the unit area, σ and µ are respectively the tension of the domain wall and the cosmic string. Hence, depending on the tunneling process or the decay time, the metastable defects could appear as (almost) stable or disappear quickly and hence leave different GW background imprints.
3. Hybrid topological defects appear when a defect associated with the homotopy π k (H i /H j ) is produced at a breaking stage and at a subsequent breaking stage a defect with a homotopy π k−1 (H j /H ℓ ) is generated. Then, the topological defects associated to π k (H i /H j ) interact with those associated with π k−1 (H j /H ℓ ) [45]. Examples of these phenomena include the well known fact that monopoles could be the seeds of monopole-antimonopole decay of a string via the Schwinger mechanism and also the fact that strings could become the boundaries of the DW produced at a later stage [45,44,46]. The interaction and evolution of these defects depend on the tensions and sizes of the interacting defects and on when inflation and reheating take place (see for example [41] for a comprehensive review).
4. If monopoles and cosmic strings are produced at the same breaking step or without any sizable gap, they are unstable [47], but could leave observable GW spectra [48,40].
5. If domain walls appear together with cosmic stings after inflation, they become unstable due to the attachment of strings and thus leave GW signals without causing a cosmological problem. Furthermore. a discrete symmetry associated with a domain wall could be lifted by Planck mass suppressed terms [49] and thus the domain wall decays to leave an observable GW signature while its cosmological evolution does not conflict with the standard cosmology [50].
Due to the recent development in GW detection, there have been extensive studies of the signatures from the GUT breaking. Distinct GW signals from cosmic strings produced in association with/without monopoles, or around/away from the end of inflation could be detected to provide an important clue of both inflation [51,52] and GUT theories. A typical U(1) sub-group of SO (10) GUT is U(1) B−L whose breaking at different scales would be tracked by distinguishing the resulting SGWB [53,54]. On the other hand, the formation of domain walls and the GW production could be a consequence of the spontaneous breaking of discrete R symmetries [55,56]. There have been some arguments against the possibility of finding R symmetries [57] in 4D GUT theories. Lastly, the breaking of supersymmetry itself it can lead to FOPT [58] that can generate GW [59].
This makes the study of GW in the context of GUT a fascinating arena that can probe some of the fundamental aspects and predictions of GUTs. As mentioned before, the study of topological defect formation, evolution and production of GW signals has a long history, however more elaborated simulations and treatment of signals are taking place recently given the opportunity for current and future experiments. 5

Typical examples of SO(10) breaking
Here we present possible breaking chains having more than one source of GW signals and that we will analyze in §4. As it is well-known, uncertainties in determining the parameters involved in the processes of GW production are yet to be greatly reduced. Nevertheless, we can still make qualitative analyses revealing important features. Let us now depict the distinct breaking chains for which we determine the intermediate and unification scales and estimate proton decay ratios: The numbers (2,1,0) in parentheses denote respectively again topological defects (monopole, string, domain wall), and (pt) means phase transition. The subscripts of the subgroups in Eqs. (4)-(6) are as follows: C for color, L for Left, R for right, V = 4 I 3 R − 3(B − L) and Z = −I 3 R + 1/2(B − L). Where I 3 R is the third component of the right-handed isospin and B − L the baryon minus the lepton number. The notation D represent the D parity also known as Z C 2 . Note that Z 2 , which is used in supersymmetry as R-parity, is never broken. Here M GUT , M R and M B−L refer respectively to the GUT scale, the scale where SU(2) R is broken and the scale where B − L is broken. The forthcoming results of JUNO [61], DUNE [62] and Hyper-Kamiokande [63] experiments will be able to tell us which models are going to be ruled out and therefore this will be able to exclude the models as sources of GW at a particular scale.
One may consider the breaking route of Eq. (5) or Eq. It is indeed interesting to look for breaking chains where combined sources of gravitational waves could appear. It turns out, however, that this does not occur generically. For supersymmetric models with one intermediate step of breaking it is difficult to separate the intermediate scale too much from the unification scale without incurring proton decay problems. With more than one intermediate scale, an intermediate scale could be lowered, e.g., between 10 10 GeV and 10 13 GeV. Then a further splitting of scales can be achieved and combined effects can take place. This is the case for the example in (6). For non-supersymmetric models, the separation of the GUT scale M GUT and any other intermediate scales, M I , M II , is more feasible basically because the number of particles in the theory, which change the coefficients in the beta functions, are much larger than in their supersymmetric counterparts and therefore one can allow more running of the gauge couplings between the intermediate scales. However, it is precisely this hierarchy, i.e. M GUT ≫ M I (or any other intermediate scale M II ), that erases potential interactions of combined effects. Specifically, the decay probabilities of Eqs. (2) and (3) decrease very rapidly because the arguments in the exponential go like −1× ratios of powers of the masses involved in the breaking chains. For example, for a decaying string at a scale ∼ M I , the factor k in Eq. (2) goes like the ratio of M 2 GUT /M 2 I , assuming that monopoles are created at M GUT . We also note that the example in ( §2.2) allows the signatures both from cosmic strings and from FOPT. However, the signals from domain walls and FOPT turn out to be away from the current sensitivity as it will be shown in §4.

Phenomenological considerations
The minimum set up for studying GUT, and its different breaking chains, calls for computation of unification scales and proton decay depending on the choices of matter content and intermediate scales.
Detailed aspects for this are well studied in the context of non-supersymmetric models. When supersymmetry is considered, one needs to make more assumptions about the spectra and boundary conditions. We will take no-scale supergravity boundary conditions following the treatment in [64].
Gauge coupling unification Except for the cases studied previously in the literature, we compute the beta functions following [65] and express them in the usual form µ dg a /dµ = b (1) where we keep the top Yukawa coupling, y t , since due to its relative size in comparison to the gauge couplings it cannot be neglected.
We assume the "survival hypothesis" implying that only the SM particles contribute to the running below the intermediate symmetry breaking scale and all other particles have masses around either M GUT or the intermediate scales.
Proton decay For non-supersymmetric theories the most sensitive channel is the dimension-6 operator induced decay p → π 0 e + mediated by gauge fields. We estimate the proton decay for this channel using [66] where m p and m π 0 are the proton and the neutron pion, π 0 , masses respectively. The amplitude at the weak scale is computed from   where the Wilson Coefficients C RL ((ud) R u L ) and C LR ((ud) L u R ) corresponds to C 1 and C 2 of [66], respectively. For numerical values we use the inputs given in Tab. 1. The most constraining channel for the supersymmetric theories is p → K +ν , mediated by dimension-5 operators. We use where m K is the mass of the kaon and the decay amplitude is given in terms of the Wilson Coefficients effective at the hadronic scale (2 GeV) It is useful to remember that all of the Wilson Coefficients in Eq. (10) are suppressed by the color-triplet Higgs mass and the supersymmetry breaking scale. The hadronic matrix elements are given in Tab. 1 and the evolution of the Wilson coefficients from M GUT is given in [64].
The forthcoming experiments DUNE [62], JUNO [61] and Hyper-Kamiokande [63] are expected to improve the current sensitivity by an order of magnitude (see Tab. 1).

SU(5) routes
Non-supersymmetric models via SU(5) routes are quite restricted by the proton decay bound in the channel p → π 0 e + since the Wilson coefficients are suppressed only by the masses of heavy gauge fields: C i ∝ 1/M 2 X (for a review see for example [69].) Thus, the discussion below is specific to supersymmetric SO(10). The following routes are candidates for the appearance of combined effects: The first one, Eq. (11), produces cosmic strings at the first step of breaking, and monopoles at the second stage of breaking. As monopoles need to be inflated away, the only GW signal would be due to a possible breaking of R-parity (Z 2 ) or FOPT. The second possibility, Eq. (12), produces monopoles at the first stage of breaking and cosmic strings at the second stage of breaking. However, the proton decay bound requires both breaking scales to be at the same scale close to 10 16 GeV. Thus, inflation needs to take place after that, leaving no imprint from GW, other than the possible Z 2 breaking. Note however that the situation can be different in the flipped models. The breaking of SU(5) F × U(1) V (and thus one U(1) factor) can occur at a much lower scale than 10 16 GeV, if we introduce a pair of vector-like particles [70]. For the third route, Eq. (13), monopoles will be produced in the first two breaking steps, while cosmic strings in the last one.
Here unification also requires all the three scales to be close to 10 16 GeV. However, just as in the case of of the flipped SU(5) pair of vector-like fields could be introduced to split the scales [70].
In supersymmetric models, the dimension-5 operator in the channel p → K +ν is mediated by the exchange of colour-triplet Higgs and suppressed by the scale of supersymmetric particles, M S , so effectively the operators mediating proton decay are proportional to ∝ 1/( M GUT M S ). This means that we can push up the supersymmetric scale to aid overcoming the present bounds [71,64]. To have an idea we can write, [69], where M S is the supersymmetry breaking scale, M HC is the colour-triplet Higgs mass, A R is a hadronic parameter and β is the usual angle determined by the ratio of the vacuum expectation values of the two Higgs doublets. For M S = 10 TeV, we get a suppression of O(10 −2 ) in Eq. (14) making the proton decay just compatible with the current bound τ (p → K + ν) > 6.6 × 10 33 yrs and thus accessible by the upcoming experiments [61,62,63]. However, it is nontrivial and it is quite model dependent to achieve the measured Higgs mass m h = 125 GeV. This puts additional constraints on the model. When U(1) V and U(1) Z are broken, there can appear B − L proton decay modes through, for example, the operator d c d c u c LH u in supersymmetry [72]. However, the resulting proton decay rate depends on the flavour structure of the model and can be suppressed sufficiently. For our analysis, we adopt the results of Fig. 3 of [64] showing the regions where the proton decay rate is compatible with the bound of the channel p → K +ν . These regions correspond generally to the mass range of M 1/2 ≥ 5 TeV and M 0 ≥ 7 TeV. Taking (M 1/2 , M 0 , tan β, µ) = (5.5 TeV, 7.0 TeV, 5, > 0), we have the GUT scale of and a proton decay rate of This would requires to break all the chains of Eqs. (11)(12)(13) directly to G SM × Z 2 at the GUT scale. The intermediate scales can be separated away from the GUT scale by altering the unification scale and introducing an effective intermediate scale with the addition of multiplets which affect a bit the running but do not have a great impact in the proton decay rate. This would allow detectable hybrid topological defects appearing in the SU(5) route. As an example, we take the chain Eq. (13) and consider the following two cases: Depending on when the inflation takes place, we may observe different gravitational wave signals. If the reheat temperature T RH is higher than M B−L , there appear the usual signal of undiluted stable comic strings. In the opposite case T RH < M B−L , all the defects are diluted away, but sizable string regrowth may occur depending on the conditions. They are basically determined by choosing the parameters that satisfy Eq. (7) of [73]. A brief summary about this is provided in the Appendix (D). Basically we assume a number of long strings with the parameter choice ofz = 10 4 which is used in the plot Fig. 1. Another interesting phenomenon is the appearance of observable signals of decaying cosmic strings attached with monopole-antimonopole pairs. The second breaking route of Eq. (17) can realize this possibility as M R is slightly above than M B−L and the inflation can take place in between, M R > T RH > M B−L . We follow [74] in order to get the profile of the decaying cosmic string. All of these features are illustrated in Fig. 1. Our breaking chain also allows the possibility of T RH > M R > M B−L which has been known to produce undetectable gravitational signals. For this, we refer the readers to the recent study in [42] .
The input values of the quark Yukawa, we use We then employ the well known procedure of running and matching the Wilson Coefficients of the relevant Lagrangian at each step. The running of the coefficients C 5L,5R from M GUT to where β Y uk L,R are the beta functions of SU we have the usual running of the Minimal Supersymmetric Standard Model (MSSM) and use the procedure employed in [64]. This means that the coefficients C * 331i 5R (µ), C jj1k 5L (µ) are evolved using Eq. (24) from the GUT scale to the supersymmetry breaking scale M S . At M S , the sfermions are integrated out via the wino-or Higgsino-exchange one-loop diagrams to get the effective Lagrangian L eff.
where the coefficients CH i and CW jk are proportional to C * 331i 5R (M S ) and C jj1k 5L (M S ), respectively. Finally, the coefficients 6 The dimension-5 effective Lagrangian mediating proton decay is given by L eff iaĒjŪkbDlc, and the supersymmetric multiplets Q, D, E, U and L are integrated over the Grassmann variable θ. of Eq. (10) at the hadronic scale, C RL (usdν τ ), C RL (udsν τ ), C LL (usdν k ) and C LL (udsν k ) are proportional to CH 2 , CH 1 , CW jk and CW jk , respectively (as given in Eq. 38 of [64]). In order to make a definite evaluation of the proton decay in the K ν channel we impose no-scale supergravity conditions at M GUT . Then, we look for values of the supersymmetric particles that satisfy the current limits. Taking M 0 = 1.6 × 10 4 GeV, M 1/2 = 2 × 10 4 GeV and tan β = 3, we get τ (p → K ν) = 8.6 × 10 33 yrs,

Monopoles Cosmic Strings
consistent with the values in Tab. 1.
First breaking stage This breaking happens at M GUT given in Eq. (18). Monopoles are produced at this scale but no CS.
Second breaking stage At this scale SU(2) and U(1) B−L are broken and this can happen at the end of inflation. For this case both decaying CS and diluted CS can occur. In the top plot of Fig. 2 we depict the breaking chain which predicts metastable strings after inflation at 10 13 GeV. The scale at which monopoles are produced is ∼ 10 16 GeV, and thus the string decay rate is exponentially suppressed by the huge number of k = v 2 /µ ∼ 10 6 where v = O(10 16 ) and µ ∼ (10 13 GeV) 2 [see Eq. (2)]. Therefore the strings are almost stable and its signature is plotted in the bottom left panel as an orange band corresponding to (2, 4) × 10 13 GeV. It may happen also that CS are produced during inflation and thus lead to the GW profile of the diluted CS. This is shown in the left bottom panel as a dashed purple line for a redshift ofz = 2 × 10 4 that satisfies Eq. (7) of [73].
First Order Phase Transitions Phase transitions have been analyzed in the context of the present breaking chain and in the context of no-scale supergravity [76], and although they are difficult to realize, there is a possibility and therefore we consider this case. Coleman-Weinberg inflation inflation could also be adopted to achieve FOPT but this has been explicitly carried out only in the context of non-supersymmetric theories [37].
In the bottom-left panel of Fig. 2 we plot a profile of a GW from a FOPT itself at 3 × 10 13 GeV. Note that the predicted GW profile is way out of the present experiment sensitivities. It will require detectors capable to access a frequency band in the kHz to MHz region. In fact, some proposals of exploring up to kHz are starting to take shape, in particular interferometers [77] and optically levitated dielectric sensors spanning a wide frequency band from few kHz to ∼ 300 kHz [78,79]. For the MHz region there are also incipient proposals, [80,81,82] for frequencies up to 100 MHz and in the range 1-250 MHZ [83].
The rest of the parameters used to produce the energy density profiles is given in Tab. 2 and the relevant formulas are given in Appendix C.
Parameters used for the FOPT profiles α β/H g * ǫ 0.   Fig. 3 we use an annihilation temperature also equal to the a breaking scale of 3 × 10 14 GeV.
Inflation framework As mentioned before, this breaking pattern has been analyzed in [84,85], and a Starobinsky-like inflation scenario has been constructed successfully with the inflation energy scale V It is also shown in [32] that the proton decay rate in the channel p → π 0 e + is above the current bound, 1.6 × 10 34 yr even without threshold corrections. Future experiments, JUNO [61], DUNE [62] and Hyper-Kamiokande [63], have the potential to exclude it. Note however that threshold corrections or additional singlets added to the theory can push-up the limits [32]. In the right plot of Fig. 2 we show the GW signal from a CS (orange curve in the lower rightcorner) that would occur at 9.5 × 10 9 GeV corresponding to the tension of Gµ = 7.6 × 10 −20 ). In this non-supersymmetric scenario inflation can take place before that and so the CS will behave like stable string, which is in the reach of BBO.
where we have included uncertainties from the electroweak parameters, M Z , α S and m t , quoted in Tab. 1. We remark that M D and M GUT tend to be closer to each other, than M R and M D , mainly because of the value of the coefficient (b  (14)). At this scale, however the dimension-six operators leading to p → π 0 e + become relevant. To calculate proton decay we proceed as in the previous section §4.2 with two matching scales, M D and M R . For the dimension-six operator, we use the procedure of [64] to obtain the proton lifetime 2 × 10 34 yrs which is just above the current bound (Tab. 1). For the dimension-five operator leading to p → K +ν , we obtain the proton lifetime 9 × 10 33 yrs with the sparticle masses at 100 TeV, which falls in the reach of the next generation proton decay experiments. This, in turns, sets an upper bound for the scale of this kind of inflation and hence serves as a guide to look for the scales for which topological defects are not completely erased by inflation and reheating.
It is known that the domain wall network can overclose the universe if it is generated after inflation. This can be overcome in various ways (see Appendix E). One possibility is that the

DomainWalls
Cosmic Strings formation of cosmic strings takes place just before inflation and the domain wall formation takes place after inflation. Then the strings can nucleate on the wall, which make domain walls to decay through a quantum tunneling process [44,43]. This process is controlled by the ratio µ/σ between the string tension µ and the wall tension σ. If the radius R of a circular string satisfies R > µ/σ, the (nucleating a loop of string) rate per unit area is given by Eq. (3).
In the present case where the D parity breaking produces Z 2 walls, we can use the results of the existing simulations. These simulations assume that all damping forces that could be present in the cosmological evolution of the string-wall network are negligible and hence the wall tension goes like σ ∼ M 3 D [41]. We can then use the formula of Eq. (55) to describe the GW signal left by the domain wall. Note that in this case an order of magnitude for the hierarchy between M D and M GUT is enough 8 to produce a large value of µ 3 /σ 2 which will make the DW evolution indistinguishable from the evolution of DW without the presence of CS.
In this model, the scale M D = 3.8 × 10 14 GeV is of the order of the bound for the Hubble parameter of single-field inflation, as mentioned at the beginning of this section. So in this case the domain walls could have been produced at the end of the inflation leaving their GW as given in Eq. (55). In the left plot of Fig. 3 we plot the DW corresponding to this scale with the rest of the parameters as in Tab. 2.
We recall that our motivation is to look for signals of multiple GW sources coming from topological defects or from FOPT from GUTs. So we ask the following question: suppose that a signal corresponding to the shape and frequency of a DW is reconstructed and has a peak frequency of about 10 5 , would this be then indicative of a DW produced by the breaking of a subgroup of a GUT preceded by a breaking producing CS? To put it in other words, we are asking if such a reconstructed signal would be a clear signal of a DW in the presence of CS. To answer the question we recall the following. It is known that DW could be also produced when a symmetry is lifted in the vacuum after the breaking of a discrete parity, such that the effective bias term drives the annihilation of the DW [87,88,89]. An specific way to do this is to consider a potential of the form V = 1/4λ φ 4 + φ 6 /Λ 2 9 , where Λ is a cutoff scale for the dimension 6 operator and φ represents the Higgs breaking of the D symmetry. In this case, even in the absence of CS, such a term could originate from the 16 of SO (10), which under the group SU(3) C × SU(2) L × SU(2) R ×U(1) B−L × D contains a singlet and hence both of the terms φ 4 and φ 6 would be allowed. Then, the tension of the string can be estimated from where φ f represents the value of the Higgs field where the two minima are quasi-degenerated. At the scale 10 13 GeV the typical values of V and φ will be respectively 10 13 4 and 10 13 , giving then as a result a tension σ of the order 10 13 GeV 3 . In this way, DW can collapse and part of their energy will be converted to GW. We can take the annihilation temperature to coincide with the scale of the breaking 10 13 GeV, assuming of course that inflation takes place above that scale. The signal will be indistinguishable from the DW signal that we illustrate in Fig. 3 we present a DW profile, using the parametric formula of Eq. (55) with the parameters appearing in Tab. 2.

Comparison to Non-Supersymmetric Case
which are to be compared with the supersymmetric case Eq. (15). The big differences in the predicted breaking scales will are maintained even after including the full uncertainties. In this nonsupersymmetric case, domain walls would be formed above 10 15 GeV and hence will not be observable. Cosmic strings can be produced after inflation at a scale 2 × 10 12 GeV, corresponding to approximately to a tension Gµ = 3 × 10 −15 . This is plotted in Fig. 3 with a solid orange line. We also checked that the proton lifetime in the p → π 0 e + channel is 1.0 × 10 35 yrs which is above the reach of the next generation of experiments (Tab. 1).

Conclusions
We studied the possibility of tracking down a route from SO(10) down to the SM in supersymmetric SO(10) models using GW signatures from the topological defects involved in the various stages of breaking, which may differ from the non-supersymmetric counterparts. Non-supersymmetric models have been studied widely as they require various intermediate steps of breaking to achieve gauge unification. Although limited, supersymmetric models may also allow some intermediate breaking scales which could improve gauge coupling unification. The shapes and strengths of gravitational wave signals depend on the scales of intermediate breaking, the scale of inflation, as well as whether or not hybrid topological defects appear. We explored such features considering the breaking routes of SU (5) In the supersymmetric SU(5) example, the only GW signatures that we could identify are from cosmic strings appearing at the intermediate scale around 10 14 GeV. These could decay due to monopoles predicted in a previous breaking step, or become diluted depending on the scale of inflation. Therefore, variant GW signals are predicted. The breaking routes of SO (10)  We conclude by commenting on the reasons why it is difficult to find examples where GW signatures from different sources can appear. Considering only topological defects, we mainly need the breaking of the GUT group down to the SM in at least three steps. The separation of these steps needs to be controlled to assure that the topological defects producing observable GW signals occur after inflation and reheating. These features are illustrated in our breaking chain 4.3. In this case, the appearance of cosmic strings and DW is separated by roughly an order of magnitude, and hence DW decay via nucleation of strings on the wall while avoiding the overclosure of the universe and producing GW at the same time. Thus, two separate GW signals arise: one from DW at a higher scale and the other from CS at a lower scale.
On the other hand, FOPT could occur at any scale and thus two GW signals can arise at a scale below inflation as shown in §4.2. Since monopoles are typically produced at the GUT scale, the scale where GW from FOPT and another source could appear has to be different from the GUT scale and so effectively at least two scales are needed. Ministry of Education through the Center for Quantum Space Time (CQUeST) with Grant No. 2020R1A6A1A03047877.

A Confrontation with GW experiments
The GW signals from stochastic backgrounds are by convention expressed in terms of a GW energy density spectrum Ω signal (f ) h 2 as a function of the GW frequency f , while the instantaneous sensitivity of a GW experiment is expressed as a noise spectrum Ω noise (f ) h 2 . Then to evaluate if a predicted signal would be able to be detected one can use one of the following procedures. I. Compute the associated signal-to-noise ratio (SNR) by integrating over the experiments' total observation time accessible to a given frequency [91,92,93]: where n det , can be 1 or 2 respectively for a cross-correlation or measurement. Then if SNR is bigger than a threshold value, it is assumed that the associated GW experiment will be able to detect the predicted GW signal. II. Construct the power-law-integrated sensitivity curve (PLISC), Ω PLISC , [5] based on Ω signal (f ). If the signal and the PLISC intersect in such a way that Ω signal (f ) > Ω PLISC for a given frequency, then is assumed that the experiment will be able to detect the signal. However, since the PLISCs are constructed on spectra based on pure power laws, in realistic situations where the signal is expected to have a structure different to a pure power law, PLISCs it can be used only as a qualitative visual aid. In cases where the shape of the predicted GW signal is fairly model independent the computation of SNR can be computed only once. This fact was exploited in [94], where new sensitivity curves for SGWB predicted by FOPT were proposed. The shape of these curves is fairly model independent (the computation depends on the contributions from sound wave and turbulence and their associated parameters). The idea is that using a fit function for a peak-integrated sensitivity curves (PISC) for a particular experiment it is not necessary to perform a frequency integration on a parameter-point-by point basis, but simply plot a numerical fit for the PISC against the predicted SWGB. Updating only when the spectral shape functions, the noise spectrum and the functional forms of the peak amplitude change. Unfortunately PISC have been obtained only for few experiments in particular LISA, DECIGO and BBO and for SGWB FOPT. We compare our examples in §4 to the profiles of these last experiments based on the information of [94] and for NANOGrav, PPTA and SKA we use smoothed profiles for the sensitivity curves based respectively on [22,21], [17] and [19], while for LIGO whe use the (PLISCs) [5] approach. As a comparison in our plots we plot the LISA PLISCs and PISC approaches. In all of our plots for LISA, BBO and DECIGO we use t obs = 1 year. For SKA, we present observation times corresponding to one, two and five years, these are shown in the plots appearing in Fig. 1-Fig. 3 respectively as the small, medium and big triangular regions delimited by green dashed lines.  126 is (1, 3, 2, 1) and that of 10 is (2, 2, 0, 1) ⊃ (2, 1, ±1/2), this last decomposition being the decomposition under the SM group, For the matter, all fermions and their corresponding sparticles in the 16 are taken into account: (2, 1, −1/3, 3), (2, 1, 1, 1), (1, 2, 1/3, 3) and (1, 2, −1, 1). Then the beta functions coefficients at one and two loops are respectively while for non-supersymmetric, we have For the non-supersymmetric case our results agree with those of [66,32]. The matching conditions at the scale M R are as those in Eq. (37).
The non-supersymmetric beta function coefficients of this example are Again for these cases, for the non-supersymmetric case, we reproduce the results of [66,32].
The matching of the LR = SU ( here where z is a real number that is determined through the constraint of unification at M GUT . All of the discussion above could change if of course we consider other than minimal models and add matter content at any step that can affect the value of the beta function coefficients without changing the scalars that define the different breaking patterns.

C Production of GW from First Order Phase Transitions
Cosmological FOPT originate from a discontinuity in the entropy when there exist a metastable vacuum that eventually decays into the true vacuum of the theory. This event occurs through bubbling, that is, the process where regions of space get first to the true vacuum state and overcome the barrier separating the two vacua of the theory. Bubble dynamics produces GW through two basic mechanisms: (i) bubble collisions, generating sound waves, and (ii) turbulence, both producing energy that releases into the GW. In the case of turbulence, it happens when bubble expansion causes macroscopic motion in the surrounding plasma. The total stochastic GW background measured from FOPT is the sum of the contributions from bubble collisions and turbulent motions. However, not all the SGWB are detectable, weak production proceeds via vacuum tunnelling and thermal fluctuations, while strong production happens when bubbles are merely nucleated via quantum tunnelling. Only this latter case is detectable.
Given the breaking of the group into another, producing hence the vacua of the effective theory, the effective potential relevant for a phase transition is , where the fields φ are all the fields necessary to parameterize the breaking phase, V 0 (φ i ) is the tree level potential, V CW (φ i ) is the Coleman-Weinberg contribution (one-loop) and V T 1 (φ i , T ) is the finite temperature correction.
The parameters carachterising the FOPT are as follows. The "vacuum-to-thermal energy ratio", α, which is roughly the ratio of the false vacuum energy density to the thermal energy density, measuring hence the amount of energy released during a FOPT (only for strong production of SGWB α ∼ O(1)) which is defined as [95] α = ∆ρ ρ R , where ∆ρ is the released latent heat from the phase transition to the energy density of the plasma background. The "change in nucleation rate", β, is a measure of the bubble nucleation rate per unit volume, its inverse is approximately the decay rate of the nucleation from the metastable vacuum to the true vacuum and hence it characterises the duration of the FOPT. The other parameters characterising the SGWB are the "nucleation temperature", T n (or T * ), and the velocity of the bubbles, v b . The parameter β is given by where S 3 is the Euclidean action for the O(3)-invariant bounce solution. In terms of the effective potential, the α parameter can be written as where ∆V eff (T n ) is the potential energy difference between the true and the false vacuum at T n . As mentioned above, the gravitational waves from the FOPT mainly include two sources 10 : the bubble collisions, producing sound waves, and the MHD (Magnetohydrodynamics) turbulence, with the total energy given by [95] The efficiency factors, k SW and k Turb. characterise the fractions of the released vacuum energy that are converted into the energy of scalar-field gradients, for sound waves and turbulence, respectively.
The bubble wall velocity v b is a function of α [96] , although this should be taken just as a lower bound since in phase transitions there exist a larger class of detonation solutions [97]. The energy density of the sound waves that are created in the plasma is given by where the factor τ sw is given by min 1 H * , R * U f is the time scale of the duration of the phase. It could be either 1/H * or R * /Ū f , where H * R * = v b (8π) 1/3 (β/H) −1 according to [98]. The rootmean-square (RMS) fluid velocity is given roughly by [99,8,100]Ū 2 f ≈ 3 4 κν α 1+α , where κ ν is the fraction of the latent heat transferred into the kinetic energy of the plasma [101] 11 . The peak frequency is given by Hz .
The energy density of the MHD turbulence in the plasma is given by where a * (a 0 ) = and the efficiency factor is given by ǫ ≈ 0.1. The present day Hubble parameter can be expressed as Finally, the peak frequency for GW produced by MHD turbulence is given by [102] f turb = 2.7 × 10 −5 β H 1 v b T * 100 GeV g * 100 1 6 Hz .
In our study we have considered just the contributions from bubble dynamics (leading) and turbulence (sub-leading).

D Production of GW from Cosmic Strings
Stable Cosmic Strings GUT based on simple gauge groups lead to the formation of topologically stable monopoles whose density is about 10 18 times greater than the experimental limit [86], dominating thus the energy density of our universe and closing it. While this kind of topological cosmological defects demands then a wash down effect such as inflation, not all monopoles need to undergo completely washout [87]. Cosmic strings can also have enormous energy. In the simplest case, which we consider, the canonical type of string (Nambu-Goto), the energy per unit length, µ, and the string tension are equal. It is expected that for strings produced at a phase transition at T c , µ ∼ T 2 c [103], where µ is the tension of the string and it characterizes the strength of the gravitational interaction of strings. A grand unified string of length equal to the solar diameter would be as massive as the sun, while such a length of string formed at the electroweak scale would weigh only 10 mg. The gravitational effects of the latter are essentially negligible, though such strings may still be of great interest, because of other types of interactions. The Nambu-Goto cosmic strings are characterized only by 11 In the small and large limit of v b , κν can be approximated as κν ≈ α (0.73 + 0.083 √ α + α) −1 , v b ∼ 1 and κν ≈ v the string tension, µ [41]. Using the Kibble mechanism, the string tension can be estimated to be (with a GUT scale of 10 16 GeV) [104,105] µ ≈ 10 −15 G T p 10 12 GeV where G is the Newton's constant. One simple approximation is to assume T p ≈ T n . When this scale is taken to be the GUT scale, roughly T p = T n = 10 16 GeV, we then get the result of CS produced at GUT scale, we get then the familiar result Another way to express the tension of the CS is to write [106,107] where M CS is the scale of the cosmic string and M P the reduced Planck scale, 2.4 × 10 18 GeV. The function B(x) is B(x) = 1.04 x 0.195 for 10 −2 ≪ x ≪ 10 2 , B(x) = 2.4 log(2/x) for x 0.01. After the formation of strings, the string loops loss energy dominantly through emission of gravitational waves. We compute the relic GW energy density spectrum from cosmic string networks from [108,109] where Ω (k) (52) the critical density is given by ρ c = 3H 2 0 /8/π/G, k is the k-mode at a frequency f. The gravitational loop-emission efficiency factor is Γ ≈ 50 [110] with its Fourier modes for cusps [111] given by [112,110] The factor F α characterizes the fraction of the energy released by long strings. We use F α = 0.1, and α = 0.1 in order to take into account the length of the string loops rendering a monochromatic loop distribution. Θ is the Heavy side step function and a(t) is the cosmological factor at a given time t. The loop production efficiency C ef f is obtained after solving Velocity-dependent One-Scale equations (VOS), with C ef f = 5.4(0.39) in radiation (matter) dominate universe [104]. The VOS equations are where k(v) = 2 1+8 v 6 andc ≈ 0.23 describes loop formation. L is the correlation length parameter (such that the energy density of the long strings is given by ρ ∞ = µ/L 2 [41]) and H is the corresponding Hubble parameter. The scaling regime is reached after three or four orders of magnitude of change in the energy scale of the universe, where we have a stable value of C ef f , see e.g. Fig.3 of [104]. The loop formation time of the k mode is a function of the GW emission time (t) and is given by Assuming the small-scale structure of loops is dominated by cusps, the high mode in Eq. (51) is given by Ω GW (f /k). The low and high frequencies of the spectrum of the GW from cosmic strings are dominated by emissions in the matter and radiation dominated eras respectively.
Diluted Cosmic Strings A standard expectation of primordial cosmological inflation is that it dilutes all relics created to unobservable levels. But this does not need to be the case, for example in [73,113] counter-examples were presented. These correspond to networks of cosmic strings diluted by inflation but that can regrow to a level potentially observable today in gravitational waves (GWs). In contrast to undiluted cosmic strings (from a stochastic GW background), the leading signal from a diluted cosmic string network can be distinctive bursts of GW. In [73] the VOS model was used together with a simplified picture of inflation and reheating to estimate the dilution of cosmic strings. The starting point is to consider the same VOS equations as in Eq. (53). Then take the correlation length L F (t F ) at the time t F , the greater of the beginning of inflation or the formation of the network. After t F , the string network parameters reach an attractor solution during inflation given by L(t) = L F e H I (t−t F ) and v(t) = 2 √ 2/π/H I /L(t) [73]. The solution corresponds to having the long string network with with HL ≫ 1 and v ≪ 1. The conditions under which there is enough regrowth are given in Eq. (7) of [73]. These are written in terms of the redshift,z, [73] at which the condition HL → 1 is achieved. The condition can be satisfied for ∆N = 0 , where ∆N represents the number of e-foldings between t F and t I , the time at which the attractor solution is satisfied. ∆N = 0 corresponds to the string forming at the start of inflation and in this case then only the number of long strings accounts for satisfying the condition of Eq. (7) of [73]. We assume this last possibility and use Eq. (19) of [73] to produce the GW profile.

E Production of GW from Domain Walls
Domain walls are sheet-like topological defects, which might be created in the early universe when a discrete symmetry is spontaneously broken. Since discrete symmetries are ubiquitous in high energy physics beyond the Standard Model (SM), many new physics models predict the formation of domain walls in the early universe. By considering their cosmological evolution, it is possible to deduce several constraints on such models even if their energy scales are much higher than those probed in the laboratory experiments.
In the context of GUT theories they can appear via the D-parity symmetry, denoted also as Z C 2 . They cannot appear after inflation has taken place, unless they are broken by an explicit parameter, by gravitational effects, or decay due to the attachment of CS produced in a previous breaking stage. For DW to which CS do not attach, there could be also the possibility that the discrete symmetry is broken but then restored at a higher temperature [114]. Another possibility is to accompany the GUT theory with an extra PQ symmetry that lifts the symmetry Z C 2 and so it effectively only acts like an effective symmetry (see for example [115] and references within it).
For some ranges of parameters there exists an alternative solution to the domain wall problem [116,117] based on the idea of symmetry non-restoration [114], which does not require any explicit breaking of the discrete symmetry. Simply because the discrete symmetry is never restored at high temperature. In this way the domain walls never form. As it is shown in [116], this can rescue some of the models such as those with spontaneously broken CP and Peccei-Quinn symmetries. Nevertheless, this mechanism is incompatible with renormalizable supersymmetric theories [118]. One way to make the bias compatible with quantum breaking [119,120] and with the de Sitter Swampland program [121,122], would be to make a bias time-dependent in such a way that it disappears after the walls disappear [117]. Nevertheless, we think that identification or lack of signals of domain walls in the expected regions would lead to either constrain of rule out the parameter space of bias parameters.
Assuming that a bias parameter, coming either by Planck suppressed terms, [49,50] or by a soft breaking, breaks the symmetry and hence allows D parity to break below the inflation scale leaving a GW signal described by where we have used the parametric form of the frequency below and above the peak and the assumptions of [115] based on the simulations of [123]. A is a parameter that can be extracted from lattice simulations,ǫ GW is a parameter based on numerical simulations for the energy density of the GW and it has a valueǫ GW = 0.7±0.4 [115], g * s(T Ann. ) are the degrees of freedom at the annihilation temperature, T Ann. . Finally σ is the tension of the DW. Note that a DW which develops after inflation and during the radiation domination era and where cosmic strings formed before inflation will have the same parametric behaviour as Eq. (55), that is Ω DW h 2 (f ) = Ω GW, max.
f p f p Peak , where the peak frequency is determined by the Hubble parameter at the decay time, f Peak ∼ a t dec. /a t 0 H t dec. [123].