Quantum nature of the minimal potentially realistic $\mathrm{SO}(10)$ Higgs model

We study several aspects of the quantum structure of the minimal potentially realistic renormalizable $\mathrm{SO}(10)$ Higgs model in which the $\mathbf{45}\oplus \mathbf{126}$ scalars spontaneously break the symmetry down to the Standard Model group $\mathrm{SU}(3)_{c}\times \mathrm{SU}(2)_{L}\times \mathrm{U}(1)_{Y}$. With a complete information about the one-loop corrections to the masses of all scalars in the theory and the one-loop beta functions governing the running of all dimensionless scalar self-couplings, the domains of the parameter space where the model can be treated perturbatively are established, along with improved bounds from the requirements of the SM vacuum stability and gauge coupling unification. We demonstrate that the model is fully consistent and potentially realistic only in very narrow regions of the parameter space corresponding to the breaking chains with well pronounced $\mathrm{SU}(4)_{C}\times \mathrm{SU}(2)_L\times \mathrm{U}(1)_R$ and $\mathrm{SU}(3)_{c}\times \mathrm{SU}(2)_L\times \mathrm{SU}(2)_R\times \mathrm{U}(1)_{B-L}$ intermediate symmetries, with a clear preference for the former case. Barring accidental fine-tunings in the scalar sector, this makes it possible to provide a very sharp prediction for the position of the unification scale and the value of the associated gauge coupling, with clear implications for the phenomenology of grand unified models based on this structure.


I. INTRODUCTION
With the upcoming generation of large-volume experiments aiming to test the potential instability of baryonic matter (DUNE [1,2], Hyper-K [3,4]), one can expect at least an order-of-magnitude improvement of their sensitivity in most of the relevant nucleon decay channels (p → π 0 e + , p → π +ν , p → K +ν , etc.) with respect to the current limits.
Unfortunately, on the theory side, these efforts are notoriously difficult to meet with good enough estimates that would, at least in principle, make it possible to distinguish among different scenarios. To this end, even the most popular models of baryon number (B) violation based on the idea of grand unification, the so called Grand Unified Theories (GUT) [5], often come short when better than several-orders-of-magnitude predictions are demanded. This has to do with various types of obstacles plaguing an accurate determination of some of the critical inputs to such calculations, namely, (i) the proximity of the unification scale to the Planck scale which, in general, enhances the effects of higher-dimensional operators inducing out-of-control shifts in, e.g., the GUT-scale matching conditions [6][7][8], (ii) the need to go beyond the leading-order approximation in the high-scale mediator masses in order to deal with the associated theoretical uncertainties in the proton lifetime estimates, (iii) the need to model the flavour structure of the relevant baryon and lepton number violating (BLNV) currents and, finally, (iv) the lack of accurate information about the B = 0 hadronic matrix elements.
While the last two issues may be alleviated to some degree by, e.g., focusing on specific observables with less sensitivity to flavour uncertainties such as branching ratios and/or neutrino production channels (case iii) and, perhaps, investing more resources to accurate lattice QCD modelling (case iv), the first two are difficult in principle. As for point (ii), higher-order calculations in the GUT context are, by definition, complicated by the typically large number of degrees of freedom in the loops, raising questions about the stability of the results obtained at any given order of the perturbative expansion. Concerning (i), there is hardly anything one can do about this issue in general.
Nevertheless, there are very particular model scenarios in which both (i) and (ii) can be addressed in a relatively satisfactory manner. Among these, a prominent role is played by the minimal renormalizable nonsupersymmetric SO (10) GUT with the adjoint 45 triggering the first stage of symmetry breaking (followed by a second stage where the rank of the gauge group is reduced to that of the SM) [9,10]. Remarkably, the structure of the scalar sector of this theory is such that the most troublesome Planck-scale-associated effective operators are entirely absent [11] and, at the same time, the underlying Higgs model is simple enough to admit a comprehensive numerical analysis. Let us note that both these features are not only vital for any sensible physics scrutiny, but they also enable one to overcome [12] a peculiar pathology that the classical version of the model has been known to suffer from for decades [13][14][15][16], namely, the instability of its SM-like vacua. For these reasons, the minimal SO(10) GUT has attracted a lot of attention in recent years with a number of interesting works touching upon its specific aspects, often within the bigger phenomenological picture; see, e.g., [17][18][19][20][21][22].
To this end, detailed studies of the minimal renormalizable SO(10) Higgs model(s) play a central role as precursors to essentially all other activities. To date, these have focused predominantly on the leading quantum corrections to the masses of the SU(2) L -triplet (1, 3, 0) and SU(3) c -octet (8, 1, 0) pseudo-Goldstone bosons (PGBs) [12], which were identified long ago as the main culprits behind the tree-level vacuum instability issues [12,[23][24][25]. It has recently been noted [26] that a third potentially problematic singlet pseudo-Goldstone mode worth detailed scrutiny pops up along the potentially realistic symmetry-breaking chains with the U(1) B−L -breaking (seesaw) scale parametrically smaller than the GUT scale M GUT ∼ 10 16 GeV. This further complicates matters, since the field in question is a member of a rich SM-singlet family of four scalars, and thus the analysis requires a thorough inspection of a 4 × 4 mass matrix along with the associated quantum corrections.
In this paper, we aim to provide the ultimate synthesis of these (and several new) aspects into a decisive and self-contained analysis of the one-loop quantum structure of the minimal potentially realistic renormalizable SO(10) Higgs model. Besides complementing the previous studies by a refined account of the pseudo-Goldstone sector, including issues related to the previously unnoticed instability also plaguing one of the SM-singlet scalars (which, in certain limits, behaves as a third pseudo-Goldstone boson in the spectrum), we calculate the leading quantum corrections to the masses of all other fields in the scalar sector along with the one-loop beta functions of all the dimensionless scalar-potential couplings. This not only makes it possible to verify the convexity of the local extrema supporting the potentially realistic SM-like vacua, but at the same time opens the door to another important aspect ignored to a large degree so far, namely, that of the perturbative stability of all the results.
Remarkably enough, such perturbativity requirements turn out to be extremely powerful in eliminating large patches of the formerly allowed parameter space. As we shall demonstrate, the model entertains a certain level of perturbative stability only in very specific limits corresponding to the breaking chains with well-pronounced SU(4) C × SU(2) L × U(1) R and SU(3) c × SU(2) L × SU(2) R × U(1) B−L intermediate-level symmetries, with a clear preference for the former. This has to do with an interesting interplay between the three SM-compatible vacuum expectation values (VEVs) available in the scalar sector, which ubiquitously appear in the form of the universal structure This structure may give large massive contributions when a hierarchy between the GUT scale (represented by the larger of ω R and ω BL ) and the seesaw scale (denoted by σ) is assumed. One possible way to retain perturbativity would be to suppress the structure's dimensionless prefactors, as was often done in previous accounts [24,25]. However, due to its ubiquitous appearance also in higher-order loop corrections with different dimensionless prefactors, we would be hard pressed to suppress all the relevant coefficients simultaneously, not least due to the presence of the gauge coupling g, whose value is dictated by unification constraints, and as such cannot be suppressed. We thus conclude that the structure of Eq. (1) itself needs to be kept under control. The possibility of a small (ω BL + ω R ) turns out to be unviable for phenomenological reasons [unification through an intermediate flipped SU (5) is unattainable], so the smaller of the two ω VEVs must therefore be hierarchically smaller than the seesaw scale; i.e., it is merely an induced VEV. This implies one of the two mentioned intermediate symmetries must be realized.
The work is organized as follows: In Sec. II, we recapitulate the salient features of the model of interest, specify its field content and scalar potential, as well as recognize the possible breaking patterns. In Sec. III, we discuss at a conceptual level various theoretical constraints bounding the allowed parameter space -nontachyonicity of the scalar spectrum, perturbativity, and one-loop gauge-coupling unification. Preliminary analysis of the parameter space based on analytical considerations, the results of our numerical scans and the accompanying discussion are presented in Sec. IV. In Sec. V, we summarize our main conclusions and provide an outlook. All technical details related to the one-loop spectrum computation (including the resulting masses in both symmetry-breaking scenarios of interest), decomposition of the relevant SO (10) representations under intermediate-scale effective symmetry groups, detailed treatment of the parameterspace constraints, and running of gauge and scalar couplings (including their one-loop beta functions) are relegated to a set of Appendices.

II. THE 45 ⊕ 126 HIGGS MODEL
The minimal potentially realistic renormalizable Higgs model of our interest features a scalar sector transforming as 45 ⊕ 126 of the SO (10). In what follows, we write the 45 as a set of real components φ ij , while the 126 is parametrized in terms of complex components Σ ijklm , with latin indices running from 1 to 10. Both tensors are completely antisymmetric, and Σ is a self-dual tensor; cf. [26] for more details. The decompositions of these multiplets into their irreducible components with respect to several subgroups of SO(10) relevant for our analysis are given in Table VII in Appendix D 1. The complexconjugate representation of Σ is denoted by Σ * . The gauge fields (including those of the SM, as well as the extra components with leptoquark/diquark characteristics relevant for proton decay) are accommodated in the 45-dimensional adjoint representation.
It is perhaps worth noting that a fully realistic symmetry-breaking pattern supporting the observed SM fermion spectrum at the renormalizable level requires at least one more scalar multiplet [27], typically the 10 of SO (10). The electroweak (EW) VEVs carried by this representation, however, do not impact the high-scale symmetry breaking. Moreover, the one-loop effective-mass contributions coming from the 10 are subdominant due to the small dimensionality of the representation. Hence, we mostly ignore such an extra vector representation in the current analysis, since its absence typically makes little difference in the highscale spectrum and the associated gauge unification constraints. Any possible implications are discussed later as the need arises. Note that the Higgs-doublet mass eigenstate of the SM in an extended scenario must live partly in the 10 and partly in the 126, which must be consistent with the doublet extended mass matrix, whose new columns and rows contain new scalar-potential parameters introduced by the extension.
A. The classical-level setup

The Lagrangian
Conforming to the notation of [26], the most general form of the Lagrangian in the unbroken phase can be written as L = L kin − V 0 , where the kinetic part is defined as for We use the definition A µ := A a µ T a , where T a denotes the SO(10) generators in the representation 10. The fundamental (latin) indices refer to the real basis of the SO(10) vector 10, and they are always written in the lower position. Summation over repeated indices is implicitly assumed.

Field content
For later convenience, we gather in Table I a list of all scalar fields in our 45 ⊕ 126 Higgs model in terms of their SM symmetry representations. Alongside the representation type, we provide the information on whether each representation is real or complex (R/C), its multiplicity # in the model, as well as the SO(10) origins of each instance.
Note that for each complex representation, we could have equivalently chosen its complex conjugate as the canonical label; our choices are purely conventional in this regard. In this paper, we shall denote mass eigenstates by the SM representation labels and add a numbered index when the state has a multiplicity greater than 1. The value of the index increases with the mass eigenvalue. This labelling scheme will be convenient in our numerical analysis, since the masses can be computed explicitly and ranked for each parameter point.

Symmetry breaking and VEVs
The scalar spectrum contains three SM singlets: two real singlets residing in φ and one complex in Σ, for a total of 4 real SM-singlet degrees of freedom; cf. Table I. Their vacuum expectation values are parametrized as (1,1,1,0) (1, 1, 3, 0) φ ≡ √ 2 ω R , I: Field content of the 45 ⊕ 126 Higgs model. Each SM representation R has its reality / complexity, multiplicity #, and SO(10) origins indicated. A dagger ( †) indicates the presence of a massless would-be Goldstone mode.
For unambiguous identification of these states, we referred to their SU(3) c × SU(2) L × SU(2) R × U(1) B−L transformation properties. The values ω BL and ω R are real because φ is a real representation, while σ is, in general, complex. Since the overall phase of Σ can be redefined without loss of generality, σ can be taken real and positive. Although this freedom is utilized in our scans of the parameter space, we retain a notation consistent with complex σ in all analytical expressions.
Because of phenomenological requirements (gauge-coupling unification and a need for a seesaw scale), the GUT symmetry is assumed to be broken spontaneously in two stages. At the unification scale M GUT , the VEVs in φ (ω BL and ω R ) break SO(10) down to one of its subgroups of rank five. The subsequent breaking to the SM, preferably well below M GUT , is then accomplished by the rank-reducing VEV of Σ (σ) which is identified with the seesaw scale. Breaking patterns associated with various VEV directions are summarized in Table II.   TABLE II: Residual gauge symmetries (in self-explanatory notation) for various VEV configurations. The 5 1 Z refers to an intermediate flipped-SU(5) stage [28,29], while in the last column, the SU(5) symmetry remains unbroken due to the SU(5)-singlet nature of σ.
The classical vacuum structure and pseudo-Goldstone modes The three mass parameters {µ, ν, τ } are connected to the three VEVs {ω R , ω BL , σ} by the vacuum stationarity conditions 1 , which take the following form at tree level: Notice the presence of the VEV structure of Eq. (1) in both ν 2 and τ . For later convenience, we define χ as the dimensionless universal ratio of VEVs present in that structure: Given the relations of Eqs. (12)- (14), the VEVs and dimensionless scalar couplings can be taken as independent input parameters that fully determine the (tree-level) scalar and gauge spectra. The key observation made in [14,15] was that the tree-level masses of scalars transforming as (1, 3, 0) and (8,1,0) under the SM group take the simple form and can thus be simultaneously made non-tachyonic if and only if i.e., in the vicinity of the intermediate flipped-SU(5) configuration; cf. Table II. This, however, triggers the usual issues with gauge unification and/or baryon number violation -either one respects the proton lifetime limits and breaks the residual flipped-SU(5) immediately by lifting σ to the vicinity of M GUT (which corresponds to the problematic one-stage symmetry-breaking pattern), or one postpones the flipped-SU(5) symmetry breaking and faces light (3, 2, + 1 6 ) gauge leptoquarks in the spectrum with all the implications for matter instability. In either case, the phenomenological constraints are practically impossible to meet. The model has thus been discarded as non-viable and it took almost 30 years to bring it back from oblivion by invoking radiative effects [12]: for small a 2 , these may remedy the tree-level tachyonicity of Eqs. (16) and/or (17) along the potentially viable breaking chains well outside the (near)flipped-SU(5) region of Eq. (18).
Remarkably enough, only recently, another tachyonic instability was revealed in the tree-level mass matrix of the SM singlets [26] assuming the seesaw-compatible regime |σ| max[|ω BL |, |ω R |]. In the σ → 0 limit 2 , one of the masses of the physical SM-singlet scalars takes the form 1 Note that for special configurations of VEVs the number of non-trivial relations can be reduced. For instance, only two independent conditions exist in the ωBL = ωR case. Thus, one of the {µ, ν, τ } parameters remains unspecified [meaning that stable points of the scalar potential with SU(5) symmetry exist for any possible value of this parameter]. 2 This limit needs to be taken carefully due to the appearance of the χ-structure in the trilinear scalar coupling τ of Eq. (14).
One assumes a2 or the smaller of the two ω scales to be taken to zero alongside σ in such a way that τ is kept fixed and sub-Planckian, i.e., under perturbative control; cf. Sec. II B.
We refer to this state as the pseudo-Goldstone boson singlet. A companion SM singlet to the PGB singlet has the same mass expression as that of Eq. (19), except for changing the sign in front of the square root, while the remaining two SM singlets are true would-be Goldstone boson (WGB) modes 3 when σ → 0. In the well-motivated |a 2 | |a 0 | limit 4 , the PGB-singlet mass can be expanded as In the same |a 2 | |a 0 | limit, the companion singlet has the mass 8a 0 3ω 2 BL + 2ω 2 R + O(a 2 ) controlled by the a 0 parameter, providing further justification of the limit a posteriori. Incidentally, the opposite |a 2 | |a 0 | regime leads to the masses of the two physical singlets [based on Eq. (19)] equal to the PGB triplet and octet masses of Eqs. (16) and (17), respectively.
The expression of Eq. (20) is positive only for ω BL −ω R assuming a 2 > 0. Hence, one reveals again that the tachyonic instabilities are absent from the tree-level mass spectrum only in the vicinity of the phenomenologically problematic flipped-SU(5)-breaking direction.

B. The quantum-level situation
As argued in Sec. II A, potentially viable scenarios require dealing with three rather than the two previously identified instabilities in the scalar spectrum. Since the SM-singlet PGB mass is buried within a 4 × 4 mass matrix, a far more elaborate account of radiative corrections to the scalar spectrum of the model is required than that available in the existing literature [26]. Hence, for the purposes of this study, we have developed a numerical code that calculates the one-loop quantum corrections to all scalars of the model, i.e., including the modes that should not suffer from any issues inherent to the pseudo-Goldstone nature of the three culprits of Eqs. (16), (17), and (19). Even though the quantum effects should not significantly affect the heavy non-tachyonic part of the tree-level spectrum, this additional information enables a comprehensive analysis of perturbativity, a feature that is seldom addressed in the existing GUT literature.
The first perturbativity issue to be addressed is the potentially large terms with the universal ratio χ defined in Eq. (15). Remarkably, this structure pops up not only in the tree-level vacuum conditions of Eqs. (12)- (14) (and by extension in the tree-level spectrum), but independently also at the level of quantum corrections.
At tree level, its tendency to diverge in various limits can be compensated by taking the accompanying a 2 parameter appropriately small (which in addition helps to keep the tachyonic instabilities of PGB states under control). All symmetry-breaking patterns identified in Table II can therefore be, at least in principle, consistently attained. 3 More precisely, one is a true WGB and one is a B − L breaking Higgs field whose mass is proportional to |σ| 2 to all orders in the perturbative expansion; see [30]. 4 The tree-level triplet and octet PGB masses are proportional to a2, so taking a2 small enables loop corrections to overwhelm them and cure their tree-level tachyonic instability.
The situation changes dramatically at the loop level where the same χ-structure appears in the (polynomial part of the) one-loop stationarity conditions [26] but with parameters other than a 2 present in the prefactors. It is difficult to keep all these contributions simultaneously under control merely by suppressing the value of relevant scalar parameters due to the presence of the (relatively large) gauge coupling g among them. The only way to retain control over such χ-terms is by (i) sticking to small fine-tuned patches of the parameter space where the scalar couplings just cancel the effects of g and/or (ii) keeping the universal ratio χ itself under control by pushing the VEVs into several "prophylactic" corners of the parameter space.

Landau poles in scalar couplings
To this end, case (i) is generally difficult to achieve because cancellations of the gauge coupling effects necessarily invoke relatively large scalar couplings. This, unfortunately, brings in another aspect of the overall perturbativity issue, namely, the potential proximity of the scalar-sector Landau pole(s) to M GUT . In order to address this, we have derived one-loop beta functions for all scalar couplings at play and used them to look for and inspect the regions where the scalar couplings are stable enough to support such a regime. Remarkably, these constraints turn out to be extremely powerful in excluding large patches of the parameter space that were formerly thought to be viable; cf. Sec. III C.

Perturbative VEV configurations
Consequently, one can expect that at the quantum level the omnipresent factor χ will have to be dealt with along the lines of option (ii) above. Hence, in what follows, we shall require which confines the viable VEV configurations to four distinct classes corresponding to four different breaking patterns in Table II: 1. |σ| ≈ max[|ω BL |, |ω R |] corresponding to approximate single-stage spontaneous symmetry breaking As already mentioned, the first two options are strongly phenomenologically disfavoured either by gauge unification constraints or by proton longevity. Hence, we shall focus predominantly on the latter two scenarios and present the results obtained in the corresponding ω BL → 0 and ω R → 0 limits. From them, it will become evident that the ω R → 0 case is disfavoured in several aspects. Therefore, a clear quantum-level preference for the symmetry-breaking chains passing through a well-pronounced intermediate SU(4) C × SU(2) L × U(1) R symmetry stage can be identified.

III. ANALYSIS OF CONSTRAINTS
In this section, we provide a systematic account of the different types of constraints that a consistent GUT theory amenable to a perturbative expansion must satisfy. We start with those for which rigorous criteria can be implemented more easily (non-tachyonicity of the scalar spectrum and gauge-coupling unification) and follow with those requiring a subjective choice of the used criterion (perturbativity).
The constraints of this section can be considered for any given parameter point of the theory. The goal ultimately is to identify the parameter-space region(s) which pass all the viability criteria. This analysis is carried out later in Sec. IV. To facilitate the readability of the paper and streamline the main text to reach our numerical results quicker, we discuss in this section the constraints used later only at a conceptual level. The interested reader is kindly referred to Appendices A and B for technical details of the implementation of the criteria presented in this section.
A. Non-tachyonicity of the scalar spectrum A consistent broken-phase perturbative expansion is developed around the true vacuum, i.e., around a (possibly local) minimum of the scalar potential at which all physical scalar masses are non-negative. Hence, a parameter point is not considered viable if some of the masses are found to be tachyonic.
To this end, we provide in Appendix A a detailed description of the procedure used to calculate the one-loop scalar spectrum of the model in any given parameter point. Two conceptual considerations are important to note here: 1. The masses of the fields associated with the U(1) B−L -breaking scale are proportional to the |σ| VEV; i.e., they naturally live at the intermediate (seesaw) scale rather than in the vicinity of the GUT scale. 5 Because of symmetry reasons, not only their tree-level masses but even the corresponding loop corrections are |σ|-proportional [30]. Since, conceptually, these should be computed in the effective field theory at the intermediate scale -either in the SU(4) C × SU(2) L × U(1) R model (for ω BL → 0) or in the SU(3) c × SU(2) L × SU(2) R × U(1) B−L theory (for ω R → 0), where both of these setups contain a significantly smaller number of degrees of freedom than the full SO(10) Higgs model -the seesaw-scale fields are expected to receive rather small loop corrections for any values of couplings in the perturbative regime. We thus simplify the analysis by taking only the tree-level expressions for the |σ|-proportional part of the spectrum.
2. In the minimal realistic scenario, the scalar SM multiplets (3, 1, + 1 3 ) and (1, 2, + 1 2 ) are eventually mixed with their counterparts from additional scalar 10's of SO (10). Although one can neglect its impact in almost all aspects of our analysis, 6 the triplet and doublet mass matrices in the model without the 10 should be treated only as subparts of larger structures. Nevertheless, it is still possible to formulate a necessary condition for the non-tachyonicity of these states even with just partial information of the complete mass matrices in this sector.
According to Sylvester's criterion [31], a Hermitian matrix M 2 is positive definite (it has positive eigenvalues) if and only if all its leading principal minors are positive, i.e., if the determinants of all upper-left submatrices of M 2 (its upper-left 1×1, 2×2, 3×3, . . . blocks up to M 2 itself) are positive. 7 Since the 10 introduces no new SM-singlet VEV, the mass matrices of (3, 1, + 1 3 ) and (1, 2, + 1 2 ) used here represent upper-left blocks of the corresponding full mass matrices in the extended models. The positivity of all their leading principal minors (and thus of the eigenvalues presented in Table IX) forms a necessary condition for the non-tachyonicity of the (3, 1, + 1 3 ) and (1, 2, + 1 2 ) SM multiplets within.

B. Gauge unification constraints
Consistency of a GUT model demands that the SM gauge couplings unify at some high scale, starting with their measured EW-scale values. Since our SO(10) Higgs model is envisioned to be employed as a subsector of such a realistic GUT model, we include gauge unification as a viability criterion for our points.
We perform the unification test top down, i.e., starting with the value of the unified gauge coupling g. We use the renormalization group equation (RGE) to compute the EW-scale values of the SM gauge couplings and accept only points whose coupling values match the experimental ones.
Technically, the computation is done using standard techniques; cf. Appendix B 1 b. Although a reasonable account of the relevant BLNV phenomenology (such as proton lifetime calculations) in potentially realistic models requires a detailed two-loop gauge running analysis (such as [25]), a one-loop approximation is sufficient for the purposes of this Higgs-model study.
Despite considering a 2-stage breaking SO(10) We therefore need to consider how best to mimic the spectrum in a realistic model, e.g., an extension with an extra 10 scalar representation. The only non-trivial issue arises for SM representation types introduced by the extension: the doublets (1, 2, + 1 2 ) and triplets (3, 1, + 1 3 ). Recall that in the current Higgs-model setting we have access only to incomplete doublet and triplet mass matrices of the realistic theory. The relevant masses (as inputs to the RGE analysis) thus cannot be fully determined, yet we can still make use of the (positive) eigenvalues M 2 S (1, 2, + 1 2 ) 1,2 and M 2 S (3, 1, + 1 3 ) 1,2,3 as computed in the Higgs model; cf. Table IX in the Appendix. The best approximation to the realistic case involves the following 2 considerations: 1. As one of the doublets plays the role of the light SM Higgs doublet, it should be removed from the heavy RGE-contributing scalar spectrum. However, there is no point in imposing the corresponding fine-tuning on either of the M 2 S (1, 2, + 1 2 ) 1,2 eigenvalues of the incomplete doublet mass matrix. What we instead do is to model their effect in the full setting by taking into account only one copy of a doublet (not two) and assign it a mass corresponding to the geometric mean of the two M 2 S (1, 2, + 1 2 ) 1,2 . As explained above, the other doublet is taken at the EW scale.
2. There exist some ambiguities related to the possible admixture of additional doublet and triplet fields from the extra 10's into the physical mass eigenstates in models with a fully realistic Yukawa sector. If these additional multiplets came exactly degenerate in mass at around M GUT [i.e., as complete SO(10) 7 An analogue of Sylvester's criterion for positive semi-definite matrices requires all principal minors of M 2 to be non-negative.
The (1, 2, + 1 multiplets], they would inflict no change at all to the position of the unification scale and only a very small (practically irrelevant) shift to the value of g. In the realistic case, the new doublets and triplets from the 10's mix with the old ones; hence, they are not exactly degenerate, and even a shift in the doublet and triplet Higgs-model eigenvalues is induced. However, the net effect on M GUT and g is still expected to be small for at least two reasons. First, the beta-function contributions of these new scalar states (both of them in the vector representations of their associate gauge factors) are minute, and thus the corresponding changes to the renormalization group (RG) running are generically subleading. Second, as all the heavy doublets and triplets are clustered around the GUT scale, the interval of scales between which the running is non-trivial (corresponding to the mass differences between the heavy doublets and triplets) is very short. Hence, in most cases the associated uncertainties in the RG evolution are negligible, and we shall not consider the effects of the extra 10's here. To summarize, we use the computed spectrum of the triplets from the Higgs model, while the treatment of doublets was described in the previous point.

C. Perturbativity aspects
Since all calculations in the model rely on perturbative methods, some type of perturbativity test needs to be performed to check for their self-consistency. The loss of order-by-order robustness in a perturbative calculation can manifest itself in different ways, so we consider a number of different perturbativity constraints. Needless to say, their definitions are typically subject to some arbitrariness, so we shall often test them at different levels of strictness (producing different datasets).
We conceptually discuss the considered perturbativity criteria one at a time in the numbered subsections below. The technical details of their implementation are found in Appendix B 1 c.

The global-mass-perturbativity test
The first obvious restriction that we impose concerns the relative size of the one-loop shifts to the treelevel scalar masses. As simple as it sounds, it is not necessarily trivial in practice for at least two reasons: • There are accidentally light pseudo-Goldstone modes in the tree-level scalar spectrum for which a large one-loop shift is not only admissible but, in most cases, even mandatory; cf. Sec. II A. Thus, we should exclude the relative shifts to these fields' masses from the assessment. In practice, the one-loop scalar-mass correction largest in magnitude is compared to the average of the heavy tree-level masses. 8 • The relatively simple effective potential methods that we use for the computation of the leading quantum corrections to the scalar masses do not, in fact, provide the fully physical one-loop masses but rather their counterparts calculated in one of the unphysical schemes such as MS. These, however, suffer from several drawbacks such as sensitivity to potentially large IR or UV logs and residual renormalization-scale dependence. As for the former, we work with the regularized one-loop effectivemass spectrum [see Appendix A, Eq. (A64)] which, in the current situation, is perhaps the closest attainable approximation to the actual physical spectrum. However, even in such a case there is a residual renormalization-scale dependence that should be kept under control.
Considering the above, we define the quantity ∆, which represents an overall measure of mass shifts: This quantity effectively compares the largest one-loop correction in the heavy fields' masses to their average; see Appendix B 1 c for further details. A necessary condition for perturbativity can be imposed by only accepting parameter points with ∆ below a chosen threshold.

Renormalization-scale dependence and stability under the RG running
Since the earlier "global-mass-perturbativity" test is not entirely renormalization-scale independent, we need to ensure that the computed one-loop scalar masses are kept under control under a change of renormalization scale. Note that this issue can be rather severe in the busy environment of grand unified models with typically many degrees of freedom "flying around" the loops. Technically, such pathologies exhibit themselves as Landau-pole instabilities in the RG flows which, at the given level of perturbative expansion, can be studied in terms of the corresponding beta functions. To this end, the complete system of the one-loop beta functions for dimensionless scalar couplings has been derived (see Appendix C) and used as a basis for the study of RGE stability of the scalar-mass spectrum.
In this context, we label the initial renormalization scale by µ R , the upper and lower scales where the RGE system blows up by µ R+ and µ R− , respectively, and define the useful perturbativity measures The quantity t + (t − ) tells us how many orders of magnitude above (below) the initial scale µ R the theory can be run in its full form (i.e., with no degrees of freedom integrated out), whilet as the geometric mean tells us the average amount of allowed running up or down. Further details are given in Appendix B 1 c. Note that checking RG stability above the Planck scale is physically not required, and a switch to an effective theory should be performed below the scale where most of the spectrum lies. Nevertheless, persistence of perturbativity under RGE demonstrates numerical robustness of the calculation. RG stability can be checked by imposing a minimum threshold of t ± ort for viable points.

Vacuum position stability
Another aspect of perturbativity, though perhaps even more arbitrary than the two discussed so far, concerns the stability of the location of the broken-phase-theory vacuum in the VEV space. On one hand, it deals with quantities which do not have a clear physical interpretation unlike masses or couplings 9 but, on the other hand, it is still quite intuitive and can be seen as a one-point complement to the two-and four-point Green's functions' constraints above. Since the vacuum position is used in the computation of one-loop scalar masses, a big shift in vacuum typically causes also a large numerical shift in the masses. Technically, the requirement that the position of the one-loop vacuum in the VEV space should not be "too far" from the tree-level one is also one of the easiest conditions to test in practice (the vacuum position is 9 Moreover, the criterion even depends on the rescaling of unphysical parameters {µ, ν, τ }, which have a "natural" normalization chosen so that −µ 2 and −ν 2 are exactly the masses of the scalar fields in the unbroken phase. determined by one-loop stationarity condition expressions that represent only a subset of those that enter the mass corrections) and, as such, it can be used as a fast first perturbativity check. The reader is referred to Appendix B 1 c for details of implementation.

Iterative pseudo-Goldstone masses
The last perturbativity constraint arises purely from the technical aspects of the one-loop mass calculations (see Appendix A), in particular, whether the regularized effective mass of Eq. (A64) in Appendix A 3 a is a good approximation of the physical mass. Essentially, the main concern arises from diagrams with pseudo-Goldstone bosons in both the outer legs and the loop, which lead to one-loop mass corrections of PGBs proportional to logs of masses of those same PGBs. The regularized-effective-mass approach breaks down for points overly sensitive to these contributions; i.e., our method is unreliable at those points, and hence, we do not accept them as valid. We check for stability by iterative computation, initially feeding the unreliable tree-level PGB masses into the logs: The first and final converged iteration should not be "too far separated" from each other.

IV. RESULTS
Having established the Higgs model in Sec. II and presented the vital considerations required for its analysis in Sec. III, we now turn to the results.
In Sec. IV A, we first discuss the results of a simplified non-tachyonicity analysis based on analytic considerations, which help to build the initial intuitive picture. This analysis includes one-loop contributions to PGB masses in a simplified a 2 → 0 regime, while allowing for |γ 2 | = 0, which was inaccessible in previous works [26]. This already introduces an important novel result.
From Sec. IV B onward, we proceed with a discussion of the full numerical analysis and its results, implementing all viability criteria from Sec. III. In Secs. IV B, IV C, IV D, and IV E, respectively, the used datasets of points, the viable parts of parameter space, the predictions for the masses, and the analysis of gauge-coupling unification are described.

A. Analytical aspects of non-tachyonicity
Remarkably enough, in the σ M GUT regime of the two relevant scenarios (ω BL → 0 and ω R → 0), one can get good insight into the full numerical results of subsequent sections by means of a semi-analytic account of the non-tachyonicity criterion.
We further assume a 2 1 in order to suppress the potentially large tachyonic contributions to the treelevel masses of the PGBs (8, 1, 0), (1, 3, 0), and/or (1, 1, 0). These then become dominated by the one-loop effects given in Table X. Note that in the full numerical scans of Secs. IV B-IV E, this assumption on a 2 is discarded, but the obtained results nevertheless still strongly prefer small values of a 2 . The a 2 → 0 limit in FIG. 1: Regions in the β 4 -β 4 plane where PGB masses (solid shapes) and non-PGB masses (areas enclosed by thin contours) are non-tachyonic in the a 2 → 0, σ → 0, and ω BL → 0 approximation for different values of |γ 2 | (indicated by different colors). In the left panel, the solid purple region depicts the β 4 -β 4 domain where PGB masses are non-tachyonic for γ 2 = 0, and one can notice no overlap with the corresponding area supporting the non-tachyonic non-PGB spectrum (stretching down and right from the thin purple line). Retreating colored contours labelled by the corresponding |γ 2 | values display areas supporting non-tachyonic non-PGBs. In the right panel, the color code denotes the minimal |γ 2 | required for all PGB masses to be non-tachyonic for each β 4 and β 4 . The black contour in both panels encloses the area in which a |γ 2 | value exists so that both PGBs and non-PGBs are non-tachyonic. For comparison, we present the results of the full numerical scans described in Secs. IV B-IV E (discrete color-coded points in the left panel). The white dots along the β 4 = 1 4 β 4 line in the solid purple region are numerical artifacts related to additional zero-mass scalars inflicting large logs.
the semi-analytic approximation is thus justified a posteriori. Furthermore, the gauge coupling at the GUT scale is fixed to g = 0.5, which is consistent with the numeric results of Sec. IV E.
With all this in hand, there are only 4 relevant parameters driving the shape of the entire scalar spectrum, namely, a 0 , β 4 , β 4 (all real), and γ 2 (complex). Among these, a 0 appears solely in the mass of the heaviest singlet [the SO(10)-breaking Higgs field], which can thus always be made positive by a suitable choice of a 0 . The rest of the scalar spectrum is then in this limit determined by only three parameters: β 4 , β 4 , and |γ 2 |. We consider two γ 2 regimes separately.
In this case, everything depends predominantly on β 4 and β 4 , and complete analytical results are available not only for the tree-level masses of non-PGBs (Table IX), but also for the one-loop masses of PGBs (Table X).
• The (tree-level) non-PGB masses are all non-tachyonic in the β 4 < 0 and β 4 < 1 2 β 4 area depicted in the left panels of Figs. 1 and 2 -the allowed region stretches down and right from the thin purple contour.

FIG. 2:
The same panels of the β 4 -β 4 plane as in Fig. 1, but for the case ω R → 0.
• The PGB masses (which for a 2 → 0 do not obtain any tree-level contribution) are all non-tachyonic in the solid purple region in the first and third quadrants in the same plots. In both scenarios (ω BL → 0 and ω R → 0), the viable regions are typically bounded from above by the non-tachyonicity of the PGB singlet and from below by the PGB triplet. Note also that the tips of the purple triangular shapes do not extend all the way to the origin of the β 4 -β 4 plane. The reason is that gauge loop contributions to the triplet and octet PGB masses cannot be made simultaneously non-negative, and for small β 4 and β 4 cannot be overcome in the a 2 → 0 limit.
The main lesson to be learned here is that the two listed regions do not overlap at all, and there is thus no way to make the entire scalar spectrum non-tachyonic in the γ 2 → 0 limit.

The γ 2 = 0 regime
For |γ 2 | = 0, analytic formulae for the tree-level non-PGB masses retain their relatively simple form, in which the complex phase of γ 2 plays no role. The one-loop PGB masses, however, have to be calculated numerically. The situation then changes as follows: • The non-tachyonicity regions for the (tree-level) non-PGB masses are given by the inequalities They are derived by applying the a 2 → 0, σ → 0 limit to the masses listed in Table IX. The above conditions introduce the boundaries depicted by a set of |γ 2 |-labelled colored contours in the left panels of Figs. 1 and 2 (as before, the viable regions stretch down and right of these contours). Interestingly, the non-tachyonic region recedes toward the lower-right corner of the β 4 -β 4 plane with increasing |γ 2 |.
• The rather complicated non-tachyonic regions for numerically calculated PGB masses are displayed in the right panels of Figs. 1 and 2 for various |γ 2 | values. It is perhaps worth noting that they retain their (β 4 , β 4 ) → (−β 4 , −β 4 ) symmetry because the relevant radiative corrections are still quadratic in both β couplings (i.e., there are no β 4 γ 2 or β 4 γ 2 mixed terms). In the right panels of Figs. 1 and 2, we plot for a given β 4 and β 4 the value of minimal |γ 2 | for which a non-tachyonic PGB spectrum is attainable. It is particularly interesting that for |γ 2 | 0.2, one can find such points even in the 4th quadrant of the β 4 -β 4 plane into which also the non-tachyonic region for non-PGBs retreats.
• The last observation provides a clear hint where to look for a fully non-tachyonic scalar spectrum.
The black contours depict the overlap of the regions corresponding to the non-tachyonic non-PGB spectrum (the receding polygons in the left panels of Figs. 1 and 2) with these newly emerging β 4 -β 4 areas supporting non-tachyonic PGB masses (the colored shapes in the right panels). Note that in doing so, we need to look for the overlap of the corresponding viable regions for each value of |γ 2 | separately; only then can these be superimposed and projected onto the β 4 -β 4 plane. Remarkably, in the ω BL → 0 case, a fully consistent region exists for 0.19 |γ 2 | 0.47 within a relatively wide β 4 -β 4 range, while for ω R → 0, a valid region is obtained for 0.14 |γ 2 | 0.29 in only a very narrow sliver in the β 4 -β 4 plane corresponding to small β 4 and relatively large β 4 . This indicates that the ω R → 0 scenario is far more restrictive, and it is correctly anticipated that this remains so even in the full-fledged numerical scans performed later.
To demonstrate the relevance of the simplified picture we have just outlined, we add into the left panels of Figs. 1 and 2 the results of the full numerical scans of Secs. IV B-IV E (where the entire spectrum has been treated numerically at one loop). One can see that the viable points are essentially where they are expected to be based on the black contours (i.e., in the fourth quadrant of the β 4 -β 4 plane with a clear affinity toward the larger β 4 and the smaller β 4 values). The slight discrepancy between the results of the simplified semi-analytic account given here and the data from scans can be attributed to a non-zero a 2 value admitted in the latter case. It is interesting that for ω R → 0 a non-zero a 2 is actually enforced (cf. Sec. IV C 4), yet the overlap of the results of the two methods is almost perfect.

B. Data from numerical scans
We now turn to the full numerical analysis and its results. We explore the space of parameters defined by the dimensionless couplings which are all assumed to be within the O(1) domain, and the dimensionful VEVs whose values are restricted by the perturbativity constraint of Eq. (21) and unification.
We evaluate the suitability of a parameter point by its viability with respect to non-tachyonicity, gaugecoupling unification, and perturbativity, as discussed in Sec. III (the technical procedure is described in all detail in Appendix A 2). The suitability criteria are numerically implemented as a penalization function, which gives zero when all criteria are satisfied. Furthermore, the penalization function rises monotonically with the quantitative size of the violation of any suitability criterion. We use a stochastic version 10 of the differential evolution algorithm to find and explore viable regions of the parameter space.
Global mass perturbativity Global mass perturbativity Since the threshold values in the perturbativity criteria are to some degree arbitrary, we performed a number of numerical scans with varying degrees of strictness. We consider two main perturbativity measures: 1. The persistence of perturbativity at different RG scales, referred to as RG perturbativity, is encoded in the quantityt; cf. Eq. (24). Intuitively, it tells us how many orders of magnitude (in powers of 10) a point can be run either up or down via RGEs before at least one of the couplings blows up. A similar measure is also t + [cf. Eq. (23)], which considers only RG running upward in scale.
2. The ratio of the largest one-loop correction to the average of the heavy masses is denoted by ∆; cf. Eq. (22). This measures global mass (GM) perturbativity.
With these definitions, a biggert (or t + ) and smaller ∆ imply a better perturbativity of the point. We impose t + > 0.5 and ∆ < 1 in all our datasets, which are conveniently listed in Table III. The main datasets B + and R + do not have any additional constraints, while B 1,2,3 and R 1 have a stricter RG-perturbativity criterion imposed in the form of an acceptance threshold fort. There are no datasets R 2 and R 3 because no points witht > 2 were found in the ω R → 0 case. Note that all datasets consist only of viable points, i.e., those passing all criteria from Sec. III.
For some datasets, we used an additional penalization of how well a perturbativity criterion is satisfied, so as to push the parameter scan to be biased with respect to this quantity; i.e., new points are accepted only when they are at least as good as the old ones with respect to that criterion. In such cases, we refer to the scans as biased. The biased datasets searching for the best values oft and ∆ are labelled as RG and GM, respectively; cf. Table III. For each dataset in that table, we denote its label, the VEV regime explored (either ω BL → 0 or ω R → 0), the RG range in terms of t + ort, the bias criterion used for optimization (if any), and the number of points in the dataset.
All numerical results are based on the datasets from Table III and are presented in the form of figures. A list of figures, alongside the used datasets for each figure and a brief description, are gathered in Table IV. For readability, we separate the results into three sections: viable regions for input parameters are identified in Sec. IV C, results for the observables (masses) are collected in Sec. IV D, and sample patterns of the unification of gauge couplings for selected points are presented in Sec. IV E.  figure, while commas separate datasets whose information is presented separately in the same figure.
Comparison of scales (dimensionful parameters) 9 B + , B + ; R + , R + Effect of non-tachyonicity of doublets and triplets 10 Points in Table V Gauge-coupling unification

C. Viable regions of the parameter space
In this subsection, we present the viable regions of the parameter space for both ω BL → 0 and ω R → 0 scenarios. As we are limited to 2-dimensional projections, the information contained in the plots can never be complete. In what follows, we thus provide two complementary perspectives: planar correlation plots for chosen pairs of parameters in Sec. IV C 1 and likelihood σ-ranges for each individual parameter in Sec. IV C 2.

Correlation plots for different pairs of scalar parameters
Altogether, there are 11 real dimensionless scalar parameters of interest: We hence choose 6 correlation pairs (with β 4 -one of the two main parameters of interest [cf. Sec. IV A] -included twice) in a way that best demonstrates the salient features of our results. Since γ 2 and η 2 are complex, they also carry phases δ γ 2 and δ η 2 . However, we omit these phases from the plots, as it turns out that the distributions of both are practically uniform on the entire [0, 2π) interval. Interesting patterns correlating the two phases appear only by employing stricter RG-or GM-perturbativity constraints in the ω R → 0 limit, in which case the parameter space strongly prefers a 2δ γ 2 = δ η 2 relation that prevents both phases from changing under one-loop RG running; see Eqs. quantifies the relative size of loop corrections to masses (Figs. 5 and 6; lower ∆ is better). The plots are produced by merging the main datasets denoted by +, which consist of unbiasedly sampled viable points, with the RG or GM biased datasets; cf. Tables III and IV. When points are overlapping, those considered better with respect to the relevant perturbativity measure are drawn in front. This allows for identification of hot spot regions, where the best points (those colored toward red) were found.
Note that a 0 does not appear in any tree-level mass apart from the (1, 1, 0) 4 and (1, 1, 0) 2 , with the latter being the U(1) B−L -breaking SM-singlet-Higgs field. The mass of this field is |σ|-proportional and only its tree-level value is relevant; see Sec. III A. It is effectively governed by the λ 0 parameter: For small a 2 that is needed to change the tachyonic character of PGBs by loop corrections, nontachyonicity requires λ 0 (α+β 4 ) 2 4a 0 in both scenarios; cf. Table IX. Then, a 0 > 0 implies λ 0 > 0.
• |a 2 | 1: As expected, a 2 is small since it controls the PGB tree-level masses (note the different scaling of the associated axes in the relevant panels). While a 2 can be of either sign in the ω BL → 0 case and can even vanish, it turns out to be strictly negative in the ω R → 0 case. We explicitly confirmed this by an unsuccessful dedicated search for viable points in the a 2 > 0 region of the ω R → 0 case. The main obstruction turns out to be the non-tachyonicity of the doublets and triplets; cf. Sec. IV C 4. Incidentally, the a 2 ∈ (−0.05, −0.01) in the ω R → 0 case implies that the triplet PGB is always non-tachyonic because M 2 S (1, 3, 0) ≈ −2a 2 ω 2 BL at tree level. Interestingly, for ω BL → 0 the points with larger RG-perturbativity ranges prefer the a 2 ≈ 0 region (cf. Fig. 3), while global mass perturbativity prefers pushing a 2 toward 0.05 (cf. Fig. 5), generating a slight tension if the scans are biased simultaneously toward both these criteria.
• 0.1 |γ 2 | 0.4: As discussed in Sec. IV A, a compact range for |γ 2 | with a lower bound of around 0.1 is expected for a non-tachyonic scalar spectrum. While RG perturbativity strongly prefers smaller values of |γ 2 | near this bound (see lower-left panels in Figs. 3 and 4), the GM-perturbativity criterion is optimized in the higher |γ 2 | region.
• λ 4 ∼ −λ 2 : This pair of quantities exhibits the strongest visible linear correlation among all parameter combinations. Its appearance is mostly due to the shape of the intermediate-scale (|σ|-proportional) scalar masses.
• General remarks on scalar parameters' domains: Except for λ 2 , λ 4 , and β 4 , the allowed ranges of the scalar parameters are typically much smaller than the standard [−1, 1] domain. 11 On the other hand, the region where all scalar couplings almost vanish is not viable. The main reason is the need to compensate for the large gauge coupling contributions in their beta functions (cf. Appendix C) that would otherwise lead to a rapid breakdown of their RG perturbativity. 12 Moreover, the tachyonicity issues when β 4 and β 4 simultaneously vanish in the |γ 2 | → 0 regime have been discussed in Sec. IV A.
Nevertheless, smaller-coupling regions are still preferred from the point of view of RG perturbativity, as seen from highert values on the color scale in Figs. 3 and 4. Global mass perturbativity in Figs. 5 and 6, on the other hand, prefers some parameters (e.g., |γ 2 | or β 4 ) to be on the larger side of their allowed ranges, indicating a complicated interplay between the tree-level and one-loop contributions to scalar masses. This makes the numerical analysis presented here not only technically necessary but also highly non-trivial.
Note that the lack of red points with hight in Fig. 4 (compared with Fig. 3) makes the ω R → 0 case significantly less favourable than the ω BL → 0 scenario from the perturbativity point of view. Remarkably enough, this is indeed consistent with the results of the highly simplified semi-analytic account of Sec. IV A.

Ranges for individual scalar parameters
An alternative way of presenting the viable regions for the scalar parameters at hand is to show the individual range each of them can take. We present these results in Fig. 7 for both the ω BL → 0 (left panel) and ω R → 0 (right panel) cases. g ω BL → 0  Table III.
Let us note that the datasets B + and R + used therein (see Table III) essentially correspond to uniform sampling of points from the viable subregion of the parameter space (due to the stochastic nature of the differential evolution sampler). Projecting such a dataset to one parameter thus represents an approximation of a marginal probability distribution in the Bayesian interpretation, effectively providing information about the volume of viable parameter space associated with a particular parameter attaining values close to a certain point. Borrowing tools from Bayesian statistics, we thus present the ranges of each parameter in terms of their highest density intervals (HDI): The vertical extent of the bars of decreasing opacity and same horizontal position represents the 1σ, 2σ, and 3σ HDIs.
Furthermore, the plots include the information obtained from multiple datasets (cf . Table IV), which is encoded by different colors. We make use of our main datasets with t + > 0.5 (in blue), as well as those where the viability criterion with stricter threshold values for the RG-perturbativity measure was imposed: t > 1,t > 2, andt > 3 colored, respectively, by light blue, green, and orange.
Note that the best points in the ω R → 0 case havet ≈ 1.86, so there are not > 2 datasets R 2,3 ; i.e., no points can be run up and down by 2 orders of magnitude on average in the renormalization scale without blowing up. As expected, increasing the strictness with respect to the RG-perturbativity measure shrinks the allowed parameter ranges, as can be consistently seen in the narrowing of the vertical bars in Fig. 7 for more constrained datasets. The complex quantum-level interplay between different parameters generates severe constraints even for couplings of seemingly little impact on the observables of our main interest if highest-level RG perturbativity is required. For instance, |η 2 | is pushed to 0 fort > 3 (in both scenarios), despite appearing only in one-loop corrections to the heavy fields' masses and the scalar-sector beta functions.
The final observation concerns the fact that the allowed ranges of certain parameters within a stricter dataset may be in an unlikely region from the point of view of less strict datasets; i.e., the HDIs of a strict dataset may not overlap with even the 3σ HDI of a less strict one. This implies that enhancing RG perturbativity sometimes requires a push toward a very particular corner of the allowed parameter space. Note that this was already indicated by the positions of the hot spots appearing at the very edges of the clouds of viable points for some parameters in Figs. 3 and 4. A prominent example of this effect is β 4 in the ω BL → 0 case.

The VEVs and the renormalization scale
We now turn our attention to dimensionful input parameters. The tree-level potential in Eq. (6) contains three dimensionful parameters µ, ν, and τ , which we compute via one-loop stationarity conditions from the three VEVs ω BL , ω R , and σ of Eq. (11). Together with the renormalization scale µ R , one has four dimensionful parameters in total.
a. The VEVs: The complex VEV σ can be made positive and real by a phase redefinition of the 126 tensor of SO(10), while the bigger of the two real VEVs ω BL and ω R can be made positive by a sign redefinition of the (real) adjoint 45. Since we are interested only in the ω BL → 0 or ω R → 0 regimes (see Sec. II B), the bigger of the ω's sets the GUT scale, and σ plays the role of the intermediate U(1) B−L -breaking (seesaw) scale. The smaller of the ω's must then be small enough to keep the universal VEV ratio χ defined in Eq. (15) under perturbative control, i.e., ensure that |χ| 1. The subdominant ω thus plays the role of an induced VEV. Since it is far smaller than the other two VEVs and it is not associated with any distinct physical scale either, we shall not pay much attention to it in what follows.
The allowed ranges for the relevant VEVs, i.e., max(|ω BL |, |ω R |) and σ corresponding to the viable points in both the ω BL → 0 (left panel) and ω R → 0 (right panel) limits, are given in Fig. 8. As before, the data corresponding to 1-, 2-, and 3σ HDIs are represented by decreasing opacity, while the colors code different levels of strictness imposed on the RG-perturbativity side; cf. Sec. IV C 2.
FIG. 8: 1-, 2-, and 3σ HDIs for the dominant ω VEV, the σ VEV, and the relevant renormalization scale µ R in the cases ω BL → 0 (left) and ω R → 0 (right). Colors encode different strictness of the RGperturbativity measure defined in Appendix B 1 c; see the legend and Table III.
From the perturbativity and tachyonicity perspective, the absolute sizes of max[|ω BL |, |ω R |] and σ play no role, as nothing changes if these parameters were freely rescaled by a common factor. Thus, the main constraint here comes from the gauge-coupling unification (cf. Sec. III) in which the max[|ω BL |, |ω R |] plays the role of the GUT scale, while σ sets the seesaw scale. To this end, one can expect that the freedom of choosing these two scales together with the value of the unified gauge coupling should, in principle, always admit good fits to the low-energy data given in Appendix B 2.
The results in Fig. 8 show that the two scenarios of interest are rather different from this perspective.
requires the GUT scale to be almost as high as the Planck scale and a very low (yet more constrained) seesaw scale. Consequently, the GUT-to-seesaw-scale hierarchy ratio is rather large. In the opposite case [i.e., for ω BL → 0 with SU(4) C × SU(2) L × U(1) R as the intermediate symmetry], this hierarchy is generally milder, and the GUT scale of ∼ 10 15 GeV is rather close to the lower bound implied by proton lifetime limits [33][34][35][36][37]. These results agree very well with previous estimates [9,10,[38][39][40] based on the minimal survival hypothesis 13 [41,42]. A more detailed account of the unification constraints is provided in Sec. IV E.
b. The renormalization scale µ R : For each point, the quantum-level scalar spectrum computation is performed at a specific renormalization scale µ R which, in order to tame potentially large logs, we choose to be the square root of the average of all heavy scalar tree-level masses-squared (weighed by the numbers of the corresponding real degrees of freedom; for technical details of the procedure, 14 see Appendix B 2). All couplings then depend on the selection of such a µ R for any particular point.
Different parameter-space points can be directly compared only when taken at the same µ R which, in principle, requires RG evolution from their specific renormalization scale(s) to the universal one (using, among other things, the beta functions given in Appendix C). This procedure would be further complicated if some of the points began diverging before they reached the common µ R or ceased satisfying some of the other viability criteria, some of which are not RG invariant. 13 It is the presence of lighter gauge bosons in the RGE that crucially contributes to gauge-coupling unification. The scalars that are accidentally light then mostly just shift the seesaw and GUT scales; see, e.g., [24,25]. 14 Note that µR is subject to iterative changes throughout the procedure because the overall scale of all the heavy spectrum must be re-adjusted to attain gauge unification, and as such, it cannot be anticipated in advance.  Table III. As it turns out, this is more of an academic interest rather than a real hurdle to our analysis because the range of µ R 's corresponding to fully viable points does not exceed half-an-order of magnitude in either of the two limits; see Fig. 8. Thus, different points can be compared right away as they are calculated at nearly identical scales. Moreover, the choice of our RG-perturbativity requirements ensures that the viable parameter points could, if desired, all be run safely to a common scale without blowing up.
Numerically, the resulting ranges of µ R are close to those of the largest VEV for both the ω BL → 0 and ω R → 0 limits. Finally, let us discuss the effect of imposing the non-tachyonicity condition on the SM multiplets (1, 2, + 1 2 ) and (3, 1, + 1 3 ), which in realistic settings should mix with extra components in order to allow for a phenomenologically viable Yukawa sector. In the minimal version, such extra degrees of freedom come from an additional 10 in the scalar sector. Consequently, the doublet and triplet mass matrices we have been working with in the 45 ⊕ 126 context are incomplete. Nevertheless, as described in Sec. III A, even in such a situation the non-tachyonicity conditions can be applied using Sylvester's criterion, and the datasets that we have been working with so far (e.g., B + and R + ; cf. Table III) were all derived with these constraints in play.
It is very interesting though to see what happens if these constraints are not taken into account. 15 For this purpose, special datasets denoted by B + and R + have been produced. Technically, these satisfy the same requirements as B + and R + , but without imposing non-tachyonicity on the doublet and triplet mass matrices. The effect of this change is best seen in the a 2 -|γ 2 | correlation plot in Fig. 9, where viable regions with and without Sylvester's criterion are compared.
One can see that the impact of Sylvester's criterion is much bigger in the ω R → 0 regime (the right panel in Fig. 9) where it leads to a significant reduction of the viable parameter space. In particular, positive a 2 is no longer available in this case (in fact, a 2 −0.01). At the same time, the lower bound on |γ 2 | (expected 15 The tachyonicity of the (1, 2, + 1 2 ) and (3, 1, + 1 3 ) multiplets was not discussed in detail in previous attempts [24,25].
on the analytical grounds in Sec. IV A) is pushed even higher, thus excluding the interesting low-|γ 2 | regions corresponding to the most favourable values oft of the RG-perturbativity measure. Note that this is not the case for ω BL → 0 (the left panel in Fig. 9) where the doublet and triplet nontachyonicity criterion does not affect the lower limit on |γ 2 | at all. This can be understood analytically by noticing that in such a limit the critical doublet and triplet fields become members of larger representations 16 of the intermediate SU(4) C ×SU(2) L ×U(1) R symmetry (cf. Appendix D and Table VIIb). Their masses must thus be identical (up to subdominant corrections from σ) to the companion fields whose non-tachyonicity is always checked.
Hence, one can conclude that the non-tachyonicity constraints imposed on the triplet and doublet scalars play a very important role in determining the shape of the viable parameter space, and they are at the core of the aforementioned preference of the ω BL → 0 scenario with respect to the ω R → 0 one.

D. Results for the mass spectrum
Next, let us turn our attention to the bosonic (i.e,. scalar and vector) spectrum of the model. As we have already seen in Sec. IV C, the criteria of non-tachyonicity, gauge-coupling unification, and perturbativity (cf. Sec. III) shrink the viable parameter space to rather small patches, and the resulting mass ranges typically turn out to be quite narrow as well. The results are given in a series of Figs. 10-12, which correspond to three distinct classes of fields with respect to their characteristic mass scales: 1. The masses of the PGB scalars (see Sec. II A 4) and those of the associated fields 17 in the two limits of our interest (ω BL → 0 and ω R → 0) are covered in Fig. 10. This class of fields is especially prone to tachyonic instabilities and thus the main motivation behind the one-loop analysis carried out in this study.
2. The spectrum of the heavy GUT-scale fields (both scalars and vectors), i.e., those associated with the first stage 18 of the unified symmetry breaking, are shown in Fig. 11.
3. The masses of the intermediate σ-scale fields associated with the U(1) B−L (i.e., second stage) symmetry breaking are displayed in Fig. 12.
In all these figures, M S or M G indicate the scalar or vector (gauge) boson nature of the multiplet whose SM transformation properties are given in the adjacent bracket (in the case of degeneracy, the multiplicity subscript follows an ascending mass order). For completeness, the masses of the (1, 2, + 1 2 ) and (3, 1, + 1 3 ) scalars are also included here, despite the fact that these may be subject to further changes in Yukawarealistic scenarios with an additional 10 in the scalar sector (cf. Secs. III A and III B). Moreover, one doublet mass (the SM Higgs doublet) must be fine-tuned to the EW scale. We simulate this effect in our spectrum by replacing ad hoc the two doublets of the model with the SM Higgs doublet and a heavy companion, whose mass M S is computed as the geometric mean of the two eigenvalues of the (1, 2, + 1 2 ) mass matrix. There are several points and observations of the results worth making here.
1. In both cases of interest, i.e., in the ω BL → 0 and ω R → 0 limits, the shapes of the bosonic spectra confirm the expectations based on the structure of the associated symmetry-breaking patterns: 16 For instance, the doublet becomes a member of the (15, 2, + 1 2 ) multiplet along with two other propagating fields transforming as (3, 2, + 7 6 ) and (8, 2, + 1 2 ). 17 Associated fields are those that in the two limits of interest belong to the same multiplet as one of the "genuine" PGBs discussed in Sec. II A 4. 18 Let us use this simplified terminology here despite the fact that we envision the breaking to occur in a single step, albeit with a hierarchy of VEVs, rather than a true multi-stage breaking due to a dynamical mechanism.
FIG. 10: The allowed 1-, 2-, and 3σ HDI ranges (corresponding to decreasing opacity) for the PGB masses in the ω BL → 0 (left) and ω R → 0 cases (right). Colors (see, e.g., Fig. 7 for the legend) represent different datasets in Table III  2. Interestingly, the GUT-scale bosonic spectrum of the ω BL → 0 scenario is significantly lighter than that of the ω R → 0 case (see Fig. 11), while the opposite holds true for the σ-associated masses in Fig. 12. This is in accordance with the VEV hierarchy given in Fig. 8. The gap between the GUT and the seesaw scale is thus much more pronounced in the latter case, amounting to about 10 orders of magnitude, than in the ω BL → 0 setting, where it is just about 4 orders of magnitude. Note that this behaviour is in accordance with the previous estimates based on the minimal survival hypothesis; cf. [40]. From a model-building perspective, the ω BL → 0 scenario is therefore again far more attractive, as one does not need to resort to large fine-tunings to attain potentially realistic flavour patterns (including realistic neutrino masses). Moreover, the proximity of ω BL to the Planck  As before (cf. Fig. 7), the color code represents different datasets of Table III corresponding to different levels of strictness of the RG-perturbativity measuret defined in Appendix B 1 c. Left arrows point to masses which fall outside the displayed ranges (see Fig. 12), and PGB labels the pseudo-Goldstone field whose mass is plotted in Fig. 10. The dashed lines denote the sample points used in Fig. 13.
scale in the ω R → 0 case raises issues with theoretical uncertainties due to enhanced contributions from D > 4 operators.
3. Concerning the relative positions and widths of the ranges corresponding to different datasets, one can see several effects in Figs. 10-12:  • Enhancing RG perturbativity typically lowers the allowed ranges of the scalar masses, especially those of the PGBs in Fig. 10. This happens mostly due to the preference for smaller scalar couplings as indicated, e.g., in Figs. 3 and 4, where perturbativity generally improves toward the origin.
• The mass ranges for the heaviest fields are relatively narrow; cf. Fig. 11. For the gauge fields, this is due to gauge unification constraining the values of the GUT-scale VEV and SO(10) gauge coupling; see Secs. IV C 3 and IV E. As for the heavy scalars, the effect can be attributed to the structure of their mass formulae, which are often dominated by a coupling that is significantly constrained by the perturbativity criteria of Sec. IV C.
• For the fields whose mass origin is less definite (such as the PGBs, for which the tree-level mass contributions often compete with the loop effects), the main effect of increasing the RGperturbativity strictness often corresponds to a shift rather than a compression of their mass ranges (the orange or green bars are just as wide as the blue or light blue ones).  Fig. 11. This, however, is beyond the scope of the current study and will be elaborated on elsewhere. Nonetheless, it is reassuring that in the obtained spectra all potentially harmful states [including the (3, 2, − 5 6 ) G and (3, 2, + 1 6 ) G vector leptoquarks] have masses well above 10 14 GeV, and thus, they do not trivially violate any direct phenomenological bounds. Finally, let us present a couple of examples of how the mass patterns described in previous sections satisfy the gauge unification criterion. For this purpose, we select two representative points from the B RG and R RG datasets of Table III that correspond to the ω BL → 0 and ω R → 0 limits. The relevant parameter-space points are specified in Table V, and the one-loop gauge running (and unification) patterns are depicted in Fig. 13.
Several remarks are perhaps worth making here: • There is a clear qualitative difference between the ω BL → 0 and ω R → 0 scenarios in the positions of the two characteristic scales (namely, the M GUT and seesaw scale) and in the clustering of the relevant states around these. This is in accord with the discussion in Secs. IV C 3 and IV D. It should also be pointed out that besides perturbativity, unification represents another important argument in favor of considering only the symmetry-breaking chains along the two special "maximally hierarchical" directions. Assuming only a single (non-SM) light threshold admitted in the bulk that can aid the (one-loop) unification, there are then just two viable possibilities 19 : either having a (3, 1, + 2 3 ) gauge boson at ≈ 10 12.5 GeV with couplings unifying at M GUT ≈ 10 15 GeV, or a (1, 1, +1) gauge boson of mass ≈ 10 10 GeV and unification achieved at M GUT ≈ 10 17 GeV. This agrees reasonably well with the results in Fig. 12. If at the same time we require that the proton-decay-mediating (3, 2, + 1 6 ) vector leptoquark remains heavy, that implies a very strong preference for either the |ω BL |, |σ| |ω R |-or |ω R |, |σ| |ω BL |-breaking pattern; cf. Table VIII for gauge boson masses. The produced scales ω BL , ω R , and σ are in very good agreement with the results of [40].
• Given the relatively shallow angle under which the three gauge couplings eventually unify, 20 one can expect that the two-loop effects (including contributions from the Yukawa couplings that we ignore here 21 ) may cause significant shifts in both GUT and seesaw scales. The figures given in this study should thus be understood as a mere first approximation to the fully physical picture. 19 Note that the contribution of vector states was crucial for unification even in scenarios with either the (6, 3, + 1 3 ) scalar in the desert [25] or the exceptionally light scalar (8, 2, + 1 2 ) [24]. 20 Interestingly, the two non-Abelian couplings actually intersect twice in the ωR → 0 scenario. 21 A rough estimate of the size of the two-loop Yukawa contributions to the relevant beta functions can be found, e.g., in [40].

V. CONCLUSIONS AND OUTLOOK
The minimal renormalizable SO(10) Higgs model with 45 ⊕ 126 has attracted much attention in the past decade [24][25][26]. It is well known that at tree level, possible SM-like vacua suffer from tachyonic instabilities either in the color-octet or weak-triplet directions. We refer to these states as PGBs, since their tree-level masses are proportional to only a single scalar-potential coupling (a 2 ). The tree-level tachyonic instabilities thus force us to consider the model at the quantum level, significantly complicating the analysis.
While the previous state of the art was the derivation of analytic PGB singlet, octet, and triplet one-loop mass formulas in the a 2 → 0, γ 2 → 0, σ → 0 regime [26,43], in this work we have developed a numerical procedure for the computation of one-loop masses of all scalar fields associated with the GUT scale. Another important result is the analytic formulae for the one-loop beta functions of all dimensionless parameters in the scalar potential.
The new computational tools have allowed us to perform a comprehensive analysis of the Higgs model taking into account the following considerations:

Non-tachyonicity:
This criterion is a rigorous requirement for the consistency of the theory. Most prone to develop a tachyonic instability are the PGB states: the octet, the triplet, and the SM singlet.

Perturbativity:
In order for the perturbative calculation in a given parameter point to be valid, loop corrections to masses, as well as the coupling values under RG running, need to be under perturbative control.
The developed numerical tools have allowed us to consider both issues by constructing appropriate perturbativity measures.

Gauge-coupling unification:
This last criterion is phenomenological and puts requirements on the mass spectrum of the theory. We consider one-loop unification only.
The results of our analysis are as follows: • We have argued in Sec. II B that perturbativity requires the structure (15) to be kept under control, concluding that the model is perturbative only in a regime where the VEV of the 45 is dominantly aligned along the SU(4) C × SU(2) L × U(1) R or SU(3) c × SU(2) L × SU(2) R × U(1) B−L direction. These two scenarios are referred to as ω BL → 0 and ω R → 0, respectively, indicating which of the two VEVs in the adjoint gets merely an induced value. The same two scenarios are preferred from the point of view of unification and proton lifetime; cf. Sec. IV E.
• Although we have found viable parameter points in both scenarios, there seems to be a preference for the ω BL → 0 case. This holds true both from the perturbativity point of view, since the more stable points (especially with respect to RG perturbativity) were found in that scenario, as well as phenomenologically due to the better-suited seesaw and GUT scales at 10 11 and 10 15 GeV compared to 10 8 and 10 18 GeV in the ω R → 0 case.
• The viable part of parameter space does not admit all parameter values to vanish. In particular, a non-tachyonic spectrum requires that γ 2 is not close to zero in either scenario (cf. Sec. IV A), contrary to assumptions in previous studies [24][25][26]. This shows that an implementation of the one-loop PGB mass calculation crucially requires this parameter to be present.
This work represents significant progress in the efforts to elucidate the viability and consequences of the SO(10) Higgs model containing 45 ⊕ 126. Having identified the viable parameter regions, the obvious next step would be to implement the full model with a realistic Yukawa sector by adding a scalar representation 10 and upgrade its unification analysis to two-loop order. The additional states would also allow us to treat properly the EW symmetry breaking. Because of the small number of added fields and numerous new parameters in the scalar potential, the Higgs-model results presented here are expected to be a good approximation for the fully realistic case. However, the addition of fermions allows for the implementation of further phenomenological constraints, e.g., a proton decay rate prediction and computation of neutrino masses (with both type-I and type-II seesaw contributions present in general).

Appendix A: Technical details of the one-loop mass computation
Many of the results throughout this paper require the computation of quantum corrections to either the vacuum of the theory or the masses of the particles. This appendix section provides the reader a detailed description and commentary of the procedure we use for their calculations.

The general procedure
The one-loop contribution to the effective potentialà la Coleman and Weinberg in the zero-momentum scheme has the following compact form (cf. [44]): where the expression V was defined for later convenience via for a matrixÂ and numeric coefficients c 1 and c 2 with the hat on A used to denote its field dependence. The bold font is adopted for all matrix quantities. We use Φ as a generic label for the vector of all scalar fields in the theory, µ R is the renormalization scale, and M 2 S (Φ) and M 2 G (Φ) are the tree-level field-dependent mass matrices of scalars and gauge bosons, respectively. The explicit expressions for their entries are given by where the indices i and j run over all the fields in the scalar sector,T a Φ represents the action of the a-th SO(10) generator T a on the (reducible) scalar representation Φ, and (a ↔ b) symbolizes the symmetric part of the expression with respect to the indices a and b, i.e., [X ab ] (a↔b) := 1 2 (X ab + X ba ). Before any further discussion of our procedure, some important technical considerations for the explicit computation of these expressions are given below.
• We reiterate that the matrices M 2 S,G (Φ) are field dependent, which means that they have not been evaluated in vacuum; i.e., no expectation values of the fields have been inserted.
• The usual way to write the scalar-mass-square matrix M 2 S (Φ) is in the basis of all real scalar degrees of freedom, so we need to consider the real and imaginary components of complex fields separately. In the above expression, we have instead written the scalar-mass matrix in a more convenient holomorphic and anti-holomorphic basis. This implies that Φ i in the derivative first runs over all holomorphic fields and then over all anti-holomorphic fields. Conversely, Φ * j are conjugates of all the fields in Φ j ; thus, they first run over the anti-holomorphic fields and then over holomorphic ones. The used expression is valid also for the special case of real scalar fields: They need to be counted only once, and since Φ j = Φ * j , factors of 1/2 for real mass matrices are correctly taken into account. In our particular model, the fields consist of Φ = (φ, Σ, Σ * ) and Φ * = (φ, Σ * , Σ). (A5) The number of (real) degrees of freedom over which the indices i and j run is 45 + 2 × 126 = 297.
• If Φ is a reducible representation, as is the case in this paper, the expression in Eq. (A4) involves a sum over irreducible representations: In contrast to the scalar-mass case, the complex degrees of freedom are taken into account only once, and a factor 1/2 is inserted for the real representation φ. The symmetrization in Eq. (A4) also assumes that the basis of the generators is real, i.e., that all matrices T a are Hermitian matrices. 22 Lastly, we use the standard GUT normalization of generators, in which the Dynkin index of the representation 10 of SO(10) equals 1. We do not write further technical details on the tensor methods or conventions in this appendix, but invite the interested reader to check the appendices of [26], to which we adhere in this paper. Also, the appendices of [45] elaborate on different bases one can use for the representation 10 of SO(10).
• The field-dependent mass-square matrices are manifestly Hermitian, as is obvious from the expression in Eqs. (A3) and (A4). Furthermore, the Hermitian expression for M 2 G (Φ) is symmetric with respect to a and b, thus resulting in a real matrix. The situation in our SO(10) Higgs model for any value of Φ is thus the following: M 2 S (Φ) is a 297 × 297 Hermitian matrix and M 2 G (Φ) is a 45 × 45 real symmetric matrix.
Our initial procedure for loop corrections follows [26]: Expanding in powers of , the one-loop expansion of the potential and the fields' VEVs read V = V 0 + V 1 and v = v 0 + v 1 , respectively, where we used units with = 1. The stationarity condition at tree level is and the one-loop stationarity condition ∂ x V | v = 0 is simplified into where we ignored the O( 2 ) terms, which formally contribute to two-loop order in the -expansion. We used an abbreviated notation where ∂ x denotes ∂/∂Φ x , and v represents the insertion of the VEV v for the fields, i.e., taking Φ = v in the result. Furthermore, the one-loop corrected scalar-mass matrix is determined by the second derivative. Expanding up to O( 1 ) yields An alternative route of applying vacuum conditions is to fix the VEVs as input parameters and solve the conditions for Lagrangian parameters. In our Higgs model, we choose that stationarity conditions solve for the parameters µ 2 , ν 2 , and τ , while we take the VEVs ω BL , ω R , and σ as inputs. This same approach was implicitly used in Sec. II A 4, as well as in previous work on the model [26]. In this context, the -expansion is applied to the {µ 2 , ν 2 , τ } parameters instead of the VEVs.
3. Insert the tree-level and one-loop vacuum into Eq. (A9) to obtain the one-loop mass matrix (entry by entry).
The challenging parts of this procedure are the last two steps. Assuming one has the full potential V 0 implemented in terms of all fields Φ, it is straightforward to obtain the tree-level field-dependent mass matrices from Eqs. (A3) and (A4) and insert them into the Coleman-Weinberg expression for V 1 in Eq. (A1). The difficult task is evaluating the derivatives ∂ x V 1 and ∂ x ∂ * y V 1 due to the matrix logarithm in the Coleman-Weinberg expression. Therefore, we ultimately seek a method to efficiently evaluate the derivatives ∂ x V and ∂ x ∂ * y V. One possible route to dealing with the matrix logarithm in V is to expand it into a power series of matrices around the identity and then apply derivatives to the series. This leads to an infinite series of nested commutators for the mass matrix. This approach was used in [26] and can provide partial analytic insights, e.g., when the series of commutators terminates. However, it is not suitable for an efficient numeric calculation in the general regime.
Notice that the evaluation of the expression V is greatly simplified assuming that the Φ-dependent matrix A can be smoothly diagonalized viaÂ where (Λ) ij =λ i δ ij (no sum over i) is the field-dependent diagonal matrix, andR is the field-dependent transition matrix. We again remind the reader that throughout this appendix a hat on top of any quantity indicates its field dependence. The expression V and its derivatives then yield Note that we write ∂ y instead of ∂ * y for simplicity -the reader should perform this trivial replacement in all formulas involving the index y as well, so as to obtain ∂ x ∂ * y V instead of ∂ x ∂ y V. Luckily, the procedure in Eqs. (A8) and (A9) requires merely V, ∂ x V and ∂ x ∂ * y V evaluated in the treelevel vacuum Φ = v 0 , so we only need to know the eigenvalues and its derivatives in vacuum; i.e., we require just , and λ i,xy =λ i,xy (v 0 ), assuming that the diagonalization of Eq. (A10) is smooth in a neighborhood of Φ = v 0 . Notice that we use a notation where all non-hatted symbols denote quantities evaluated in vacuum, i.e., numeric quantities. We specify a numeric algorithm for calculating λ i , λ i,x , and λ i,xy in Appendix A 2.
Once the numeric values for λ i , λ i,x , and λ i,xy are available, their use in vacuum-evaluated Eqs. (A11)-(A13) still involves certain subtleties. Since the field-dependent matrices M 2 S,G (Φ) are Hermitian, the eigenvalues λ i are always real. This is expected, as they correspond to mass squares of particles. However, the eigenvalues come as arguments into logs, so questions arise on how to properly deal with cases when they vanish or are negative. There is an issue even when they are positive and very small compared to µ 2 R , since log contributions are then "unphysically" enhanced. We address how to deal with these log cases in Appendix A 3.
Having established the calculational procedure, we now reflect on how to optimize the calculation. Notice that the expression V defined in Eq. (A2) and its derivatives are invariant under a basis change ofÂ, and ifÂ is block diagonal, they can be evaluated on each diagonal block ofÂ separately. The matrixÂ stands for either M 2 S (Φ) or M 2 G (Φ). Vacuum can already be inserted for all fields except for the Φ x and Φ * y states due to derivatives, so we need a list of matrices M 2 S,G (Φ x , Φ * y ) for all relevant pairs of indices x, y compatible with gauge symmetry. We permute a basis of each M 2 S,G (Φ x , Φ * y ) in such a way that it is block diagonal. Then, we traverse the entire list of blocks for all x, y and collect only those which are formally different (treating field labels Φ x and Φ * y as dummy variables). Hence, the evaluation of V on each formal block needs to be computed only once, provided we keep track of their multiplicities for each fixed x and y.
To recapitulate, our calculation of one-loop masses for a particular parameter point thus involves the computation of ∂ x V(v 0 ) and ∂ x ∂ * y V(v 0 ) for a few thousand formally different blocks of various sizes. If the diagonalization of Eq. (A10) for that block is smooth at Φ = v 0 , we compute λ i , λ i,x , and λ i,xy with the algorithm described in Appendix A 2. They are then inserted into Eqs. (A11)-(A13) evaluated at v 0 , with proper logarithm treatment described in Appendix A 3.
Only a few small blocks were found not to be smoothly diagonalizable at v 0 , all of them located in the gauge mass matrix M 2 G (Φ x , Φ * y ). Fortunately, it was possible to compute the expression V for these blocks analytically. 23 2. Details of the numerical procedure a. Preliminaries and notation As discussed in Appendix A 1, we are interested in finding a procedure where the starting point is a Φdependent matrixÂ, and the desired result are the eigenvalues and their first-and second-order derivatives evaluated in the tree-level vacuum v 0 : λ i , λ i,x , and λ i,xy . Equivalently, expanding the diagonal fielddependent matrixΛ in Eq. (A10) into a Taylor series around vacuum viâ we seek the diagonal-matrix coefficients Λ, Λ ,x , and Λ ,xy . Note that Φ x goes over all real degrees of freedom when the indices x, y, etc., are summed over.
Obtaining the eigenvalues λ i is not difficult, since the tree-level vacuum can simply be inserted intô A and then the numeric matrixÂ(v 0 ) can easily be diagonalized. Finding the eigenvalue derivatives, however, is complicated in the most general case by possible degeneracies. In particular, suppose we have an eigenspace of Λ that corresponds to a particular eigenvalue. While any basis of this eigenspace can be used for diagonalizing Λ, derivatives such as Λ ,x may impose a particular preferred basis in this eigenspace -in other words, the transition matrixR is field dependent and higher-order coefficients in its Φ-expansion matter. The proper (Φ-dependent) eigenbasis therefore gets progressively resolved only when higher derivatives of degenerate eigenvalues are taken into account. For a more detailed discussion, we refer the interested reader to the literature [46][47][48]. We mimic the procedure from [48], in particular to derive the algorithm for the numeric evaluation of the expansion coefficients in Eq. (A14) up to 2nd order.
Before describing the algorithm and providing a quick derivation, we set up our compact but convenient notation and definitions. As part of the algorithm, we shall obtain a sequence of progressive basis transformations P n . We denote the cumulative transformation up to the n-th level by R n : For any matrix X from the set we define its form in a different basis by We combine the basis and derivative labels in the subscript via, e.g., A n,x , and refer to such a matrix as being in the n-basis. The matrix P n transitions from the n-basis to the (n − 1)-basis. It turns out each matrix P n is obtained in the algorithm by diagonalizing some particular numeric matrix; cf. Eqs. (A25) and (A27)-(A29). We write this diagonalization generically as where book-keeping subscripts under A and D correspond also to the underlying basis of the matrix, but the transformations among them analogous to Eq. (A18) do not apply. The diagonal matrices D n are the matrix coefficients of interest in the expansion ofΛ. Explicitly, for fixed x and y we have Crucially, at each step n there is an associated structure that A n−1 possesses.
• There exists a partition of the (n − 1)-basis already from the previous step. It is referred to as the (n − 1)-partition, and the associated block structure is composed of (n − 1)-blocks. Edge case: The initial 0-basis has the trivial partition corresponding to a single block.
• The matrix A n−1 turns out to be (n − 1)-block diagonal. Consequently, the transition matrix P n that diagonalizes it can also be taken as (n − 1)-block diagonal.
• The transition matrix P n defines the n-basis. Since P n is (n − 1)-block diagonal, the (n − 1)-partition also applies to the n-basis. We order the n-basis so that within each (n − 1)-block the eigenvectors belonging to the same eigenspace of D n are grouped together. This subdivision of the n-basis defines the n-partition, which is clearly a refinement of the (n − 1)-partition of the same basis. The transition matrix P n is arbitrary only up to basis changes within n-blocks, within which P m for m > n can make changes to the basis in later steps.
• Incidentally, knowledge of P n suffices to obtain the full transition matrix R n = R n−1 P n transforming the n-basis to the initial 0-basis, where for n = 1 we have R 1 = P 1 . Intuitively, R n simultaneously diagonalizes A k for all k = 0, . . . , (n − 1) if the A matrices are rewritten in the 0-basis. Note that R n is in general not block diagonal.
The above underlying structure therefore yields a progressive series of basis partitions, where the n-partition acts on the n-basis. The n-partition is obtained in the diagonalization procedure of A n−1 and is subordinate to the m-partition for any m < n. Consequently, P n is m-block diagonal for m < n. The inductively defined structure is self-perpetuating from step to step as long as A n is indeed n-block diagonal, which is true by construction.
For any matrix X n in the n-basis, we denote the block-diagonal and block-off-diagonal parts of the n-partition in the matrix by [n] and (n) in the superscript, respectively, where capital indices I and J are n-block labels. For m > n we further define the notation which represents taking the n-block (off-)diagonal part of X n and further transforming it from the n-basis into the m-basis. Equivalently, we can obtain X

[n]
m or X (n) m by first fully rotating X to the m-basis, and taking the (off)-diagonal part of the n-partition only in the end. This works due to the hierarchy in block structure; i.e., the n-partition can be applied to any m-basis for m > n.

b. The algorithm
Starting with a field-dependent matrixÂ ≡ A(Φ), we want to obtain λ i , λ i,x , λ i,y , and λ i,xy for fixed x and y. Following the notation introduced in Appendix A 2 a, perform the following steps: 1. Given the matrixÂ, compute (in the tree-level vacuum) the numeric matrices A 0 , A 0,x , A 0,y , A 0,xy .
6. Obtain Λ ,xy by diagonalizing the matrix on the left side of the expression The dot (·) denotes entry-wise multiplication. This gives us the sought-after λ i,xy , as well as the incidental transition matrix P 4 and the corresponding 4-partition.

c. Derivation of the algorithm
To derive the algorithm we have specified in Appendix A 2 b, we start by considering the field-dependent diagonalization in Eq. (A10), which we rearrange intô The matrixÂ is provided, while bothR andΛ are not known explicitly, but their full form is not necessary to compute the desired Λ, Λ ,x , Λ ,y , Λ ,xy for fixed x and y. We remind the reader that all hatted quantities imply field dependence, while the non-hatted quantities are evaluated in the tree-level vacuum and are numeric for given input parameters. We assume all notation defined so far in Appendix A 2 a, in particular, Eqs. (A14)-(A24). For later convenience, we further amend the definitions witĥ where X denotes any combination of derivative variables X ∈ {x, y, xy}. It turns out we need N = 5 progressive transformations in total, with all field dependence then pushed toP 5 , while the prior P transformations are used in the process of determining Λ, Λ ,x , Λ ,y , Λ xy . The above definitions implŷ for any n = 1, 2, 3, 4, along with the trivial edge caseΓ 1 =R. Intuitively, R n is the rotation from the initial into the n-basis, while the remaining part of the fullR diagonalization is represented byΓ n+1 . The useful quantitiesŜ n;X represent the X-derivative ofR unrotated byR from the left, with the result taken in the n-basis, while the unhatted S n;X is the same quantity evaluated in the tree-level vacuum. Lastly, we define To perform algebraic manipulations in the subsequent derivation, the following set of commutators will likely prove useful to the reader: The commutators' identities with the P n matrices follow directly from the block structure of P n and the hierarchies among n-partitions described in Appendix A 2 a. The P m is n-block diagonal for m > n, while the matrices D n identified in Eq. (A20) are diagonal with the same eigenvalues in a given n-block, leading to vanishing P-commutators in Eqs. (A36)-(A39) in particular.
As for commutators with Γ n , the argument needs to be a bit more subtle, since they are used in deriving the A n matrix at a step where P m for m > n has not yet been defined. Consider the lowest non-trivial Γ, which is Γ 2 . Since A 0 = A 0 , we have both R =R(v 0 ) and R 1 = P 1 diagonalizing it according to Eq. (A30) evaluated in the vacuum and the definition of P 1 in Eq. (A19), respectively. Thus, with the two rotation matrices related via R = R 1 Γ 2 due to Eq. (A35). The matrix Γ 2 thus represents the arbitrariness in the transition matrix, which must commute with the matrix being diagonalized (cf. [48]): [Λ, Γ 2 ] = 0. The matrices Γ m for m > 2 are related to Γ 2 recursively via Γ n = P n Γ n+1 , so their commutators with Λ can be derived from the Γ 2 -and P-commutators. Analogous arguments can be used to derive Γcommutators with derivatives of Λ, since Γ n represents the arbitrariness in the simultaneous diagonalization of the matrices A m for all m between 0 and n − 2.
We now perform a sequence of steps.
2. Take the derivative ∂ x of Eq. (A30): Multiply from the left withR −1 , rearrange and insert suitable definitions forŜ to obtain We now formally multiply Eq. (A42) withΓ 3 from the left andΓ −1 3 from the right, thus undoing the field-dependent rotationR up to the as-yet-unknown R 2 = R 1 P 2 . We then evaluate that expression in vacuum and multiply it with P 2 from the left to obtain We consider this matrix equation in terms of 1-blocks, using labels K 1 and L 1 for 1-block rows and 1-block columns on each side, respectively. Since P 2 is 1-block diagonal, we can write (P 2 ) K 1 L 1 = (P 2 ) K 1 K 1 δ K 1 L 1 . Similarly, (Λ) K 1 L 1 = λ K 1 δ K 1 L 1 ; i.e., each diagonal 1-block in Λ is proportional to the identity matrix, while the off-diagonal blocks are zero.
We first examine the 1-block-diagonal part of Eq. (A43). For K 1 = L 1 , the right-hand side vanishes, since Λ acts as an identity block multiplied with the same proportionality factor in both terms, i.e., λ K 1 = λ L 1 . Hence, Since P 2 and Λ ,x are already block diagonal with respect to the 1-partition, the equation can be written in full-matrix form as Diagonalizing A [1] 1,x yields P 2 (and thus R 2 ) and the first derivatives Λ ,x . We now return to the 1-block off-diagonal parts in Eq. (A43), i.e., K 1 = L 1 . Since P 2 , Λ, and Λ ,x are all 1-block diagonal, the second term on the left-hand side vanishes, and we get Solving for (S 2;x ) K 1 L 1 gives where the dot (·) denotes entry-wise multiplication and the auxiliary matrix Ω is defined in Eq. (A26). This can be written as a full-matrix equation via Note that the knowledge of the 1-block off-diagonal parts of S 2;x is sufficient for all further steps in our derivation of the quantities of interest.
3. Perform the following sequence of operations: Take the derivative ∂ y of Eq. (A30), multiply witĥ R −1 from the left, then formally multiply withΓ 4 from the left andΓ −1 4 from the right, and finally evaluate in vacuum. These steps give the analogue of Eq. (A43): Similar to the case of ∂ x , we consider this equation in terms of 2-blocks. For diagonal blocks, we get yielding the full-matrix equation indicating that diagonalization of A [2] 2,y produces the transition matrix P 3 and the eigenvalue derivatives Λ ,y . The block off-diagonal part of Eq. (A49) analogous to the ∂ x case gives in terms of the coarser 1-blocks with K 1 = L 1 , resulting in In full-matrix notation, this can be written as 4. We now turn to the 2nd derivative and take ∂ x ∂ y of Eq. (A30): A ,xyR +Â ,xR,y +Â ,yR,x +ÂR ,xy =R ,xyΛ +R ,xΛ,y +R ,yΛ,x +RΛ ,xy .
Multiplying this equation from the left withΓ 5 , from the right withΓ −1 5 , evaluating it in vacuum, making use of the commutators in Eq. (A39) and expressing Λ ,xy yields Since Λ ,xy is diagonal (due toΛ being diagonal by definition), knowledge of diagonal 3-blocks in Eq. (A57) suffices: To derive the first line, we used the fact that the diagonal matrices Λ, Λ ,x , and Λ y respect the 3-partition; i.e., each of their diagonal blocks is proportional to an identity matrix block, and thus they commute with any matrix on the 3-block-diagonal part. The second line was derived by taking Eqs. (A43) and (A49) into account, properly transformed to the 3-basis.
Eq. (A58) gives the diagonal 3-blocks of Λ ,xy in terms of an expression that is sandwiched between P −1 4 and P 4 . Since P 4 is 3-block diagonal, and (Λ) K 3 L 3 = λ K 3 δ K 3 L 3 , we can write the result more explicitly as with an implied sum over L 3 on the right-hand side. For K 3 and L 3 belonging to the same 1-block, i.e., K 3 , L 3 ∈ K 1 , the factor λ L 3 − λ K 3 becomes zero, since the eigenvalues in 1-blocks are identical. Hence, only the 1-block off-diagonal entries of S 3;x and S 3;y need to be considered. These are available from Eqs. (A47) and (A54) transformed to the 3-basis if necessary (which is subordinate to the 1partition of Ω). This finally gives The λ L 3 − λ K 3 factor removes one entry-wise Ω multiplication while preserving only the 1-block offdiagonal parts. The entire equation for diagonal 3-blocks of Λ ,xy can then be extended to a full-matrix equation by using the block-diagonal extraction operator [3] on both sides. By taking into account that Λ ,xy and P 4 are 3-block diagonal already, we get 3,y − (Ω · A 3,y ) A We have thus derived the final result of Eq. (A29). Note that this expression is manifestly symmetric with respect to x and y.

Detailed handling of logs
The procedure of Appendix A 1 computes one-loop corrections only to the mass parameters in the scalar potential, so the resulting quantity is the one-loop effective mass rather than the physical mass (also referred to as the pole mass). The latter requires also knowledge of momentum-dependent contributions coming from self-energy diagrams, which we do not consider. It is possible, however, to improve upon the effective-mass calculation by carefully considering its logarithmic contributions. This improved version is referred to as the regularized effective mass.
One-loop effective scalar masses in Eq. (A9) inevitably contain contributions of the form where m l is the tree-level mass of a scalar or gauge field in the loop part of a Feynman diagram. The following situations require further attention: 1. The particles in the loop are much lighter than the rest of the spectrum, e.g., when their masses are near the intermediate scale.
In such a case, the regime m 2 l µ 2 R causes logarithmic contributions to be unphysically large. This effect is cancelled in physical masses by self-energy contributions, which we do not have access to.
2. For WGBs, m 2 l = 0, so the logs contain an IR divergence unless their prefactor expressions vanish as well. 24 3. At least one of the PGBs has a tree-level tachyonic instability, i.e., m 2 l < 0, which causes the argument of the log to be negative.
We prepare the ground for solving these issues in Appendices A 3 a and A 3 b and then present a comprehensive list of how to practically treat all the cases in Appendix A 3 c.

a. The regularized effective mass
We define a quantity called the regularized one-loop effective scalar mass M 2 S,reg : We take the expression for the one-loop effective scalar mass M 2 S but solve the problem of IR divergences by emulating the shift from the effective to the physical mass.
In particular, we take inspiration from the Abelian Higgs model [49], which suggests using the replacement where m S is the tree-level mass of the particle on the outer legs (i.e., the particle whose one-loop mass we are computing) and c := We refer to this procedure as taming the logs.
b. Physical vs. regularized effective mass The regularized one-loop effective mass is the best possible approximation to the actual physical mass without knowledge of self-energy contributions. A crucial observation for inferring the non-tachyonicity of the physical spectrum is that in the perturbative regime, the non-tachyonicity of the regularized effective mass implies that the physical mass is also non-tachyonic. We schematically argue this point below.
Let us consider the continuous Fourier-transformed two-point one-particle irreducible (OPI) Green's function Γ (2) MS (p 2 ) evaluated in the MS renormalization scheme, where p 2 is square of the outer leg particle's four-momentum. In what follows, we suppress the matrix structure in mass expressions for simplicity. In the effective potential approach, we obtain the effective mass of a field, which is simply the zero-momentum value of the Green's function: Subsequently, schematically holds for the regularized effective mass. The physical mass though is defined as a solution to the equation In addition, the behaviour of the two-point Green's function is constrained in such a way that dZΓ (2) MS dp 2 where Z is the field-strength renormalization. We require that there exists only one solution to Eq. (A67), and we stay in the perturbative regime, where higher-loop corrections are subdominant with respect to the tree-level values. The latter tree-level case corresponds to For valid parameter points, which pass the perturbativity and non-tachyonicity check, the situation for all mass eigenstates thus corresponds to the red-line scenario.
c. Practical aspects of the one-loop scalar-mass calculations All potential issues related to the presence of vastly different scales in the one-loop mass calculation (due to lighter fields in the loop), as well as methods of their resolution are identified in Table VI. This resolves cases 1 and 2 in Appendix A 3 proper. If the light mass is way lighter than m 2 S (denoting the tree-level mass the loop correction is calculated to) potentially large logs may emerge.

A:
The tree-level mass will be used for all practical purposes. B: Large logs will get replaced (as for WGBs) via substitution (A64). PGBs (pseudo-Goldstone bosons), such as those discussed in Sec. II A 4 Light fields in the loops may produce large logs as they are not properly tamed by using substitution (A64) with tree-level PGB mass inserted in.
The same as above for the light fields associated with intermediate-symmetry breaking.
A: An iterative approach is employed, a with the PGB tree-level masses in (A64) substituted by one-loop PGB masses calculated in the previous iteration; cf. Appendix B 1 c. B: Via Eq. (A64).
a Some potentially large logs vanish in the two limits of interest due to the Ti and Oi prefactors of Ref. [26], Appendix B.
The arguments of the logs are always taken in absolute values, so that the issues with a negative log argument are avoided (case 3 of Appendix A 3 proper). Since we check the non-tachyonicity of the spectrum for every point at the one-loop level, and our perturbativity requirements demand the loop contributions to GUT-scale particles to be smaller than the tree-level values, the absolute value is relevant only for PGB fields in the loop. Furthermore, since the PGB contribution is tamed in the computation of masses for heavy fields, we only need to worry about this problem when computing corrections to PGBs coming from PGBs in the loop. In such a case, the tree-level PGB approximation is considered bad anyway, and we perform a further iterative procedure with an associated perturbativity check; cf. item 4 in Appendix B 1 c.

Appendix B: Details of the numerical analysis
In this appendix, we address the technical details of the viability criteria formulated in Sec. III. In particular, we outline the implementation of various parameters describing the quality of fits presented in Sec. IV B along with the penalties associated with the potential violation of the tachyonicity and perturbativity criteria in the step-by-step numerical procedure.

The viability criteria of Sec. III
a. Non-tachyonicity of the scalar spectrum Non-tachyonicity of the mass spectrum is an essential consistency criterion. To this end, we use the regularized one-loop effective mass defined in Appendix A 3 for all the scalars, except for those fields which are parametrically lighter due to their association with the intermediate-symmetry scale σ. Since none of these states is a pseudo-Goldstone boson, we can afford using the tree-level formulae for their masses instead; see Table VI.
The contribution to the penalty of a parameter-space point from the possible scalar spectrum tachyonicity is defined as where A 1 ∈ R + is a weight parameter, j runs over all the scalar fields, and s j stands for the expected scale of M 2 j , i.e., s j = ω 2 max (ω max := max [|ω BL |, |ω R |]) for the heavy fields associated with the GUT scale, s j = ω 2 In this study, the SM gauge couplings are expected to be unified at the one-loop level; cf. Sec. III B. The evolution of the usual α-factors with the EW-scale boundary conditions is then driven by the equations The dimensionless running parameter t is defined as with µ R denoting the running renormalization scale, and correspond to the beta functions of the "reduced" one-loop SM gauge couplings. The three sums above run over the gauge, fermion, and scalar fields, respectively, T (R i ) are the Dynkin indices of representations R i with respect to the group factor G i ∈ {SU(3) c , SU(2) L , U(1) Y }, and D(R i ) denote the relevant multiplicities of R i . We take κ F = 1 2 or 1 for Weyl or Dirac fermions, respectively; likewise η S = 1 2 or 1 for real or complex scalar fields.
The solution of Eq. (B3) is trivial: where we used the abbreviation α −1 The gauge-coupling-unification constraint thus reads α −1 i (t GU T ) = α −1 (t GU T ), i = 1, 2, 3, where t GU T is any scale at or above the heaviest threshold non-trivially influencing the gauge coupling evolution. Thus, relation (B6) can be rewritten as where the last subleading term takes care of the fact that the gauge boson masses are calculated with the SO(10) gauge coupling evaluated at the renormalization scale µ R , which is optimally chosen for each point separately (to be specified in Appendix B 2). In addition to all details of the bosonic spectrum, a good unification pattern is characterized by 3 main parameters, namely, t GU T , t σ = 1 2π log |σ| M Z , and α −1 . These three quantities correspond to three basic degrees of freedom which can be manipulated by: (i) scaling all dimensionful quantities {ω BL , ω R , σ, µ R } by a common factor, (ii) shifting the intermediate-scale VEV σ, and (iii) modifying the initial condition for the SO(10) gauge coupling. Parametrizing the associated shifts by Eqs. (B7) become the following set of three independent conditions for these parameters 26 : These are rather easy to solve for most of the initial choices of t GU T , t σ , and α −1 , which together with the resulting shifts yield the desired gauge unification pattern.

c. Perturbativity constraints
As far as the perturbativity constraints of Sec. III C are concerned, we implement four simple tests that quantify the level of our satisfaction 27 with the overall perturbativity. The associated considerations in the order of Sec. III C are as follows: 1. Global-mass-perturbativity test The penalty p 2 associated with the Global-mass-perturbativity test is defined as where A 2 ∈ R + is a weight factor, H(x) := x θ(x), θ(x) is the Heaviside step function, ∆ is defined in a repeat of Eq. (22) by with M 2 ij,one-loop and M 2 ij,tree denoting one-loop and tree-level scalar-mass-matrix elements, respectively. The symbol M 2 heavy denotes the average (over real degrees of freedom) of the heavy tree-level scalar masses. In physical terms, a penalty is awarded if the maximal one-loop correction of the mass square is larger than the average tree-level mass square.

Perturbativity of the RG evolution
For a given individual point in the parameter space, all calculations are done at a certain "optimal" renormalization scale µ R . For measurable quantities (such as pole masses), the specific choice of µ R should not matter. However, due to our use of the effective potential approach and proxy quantities such as regularized effective masses, the µ R -dependence is technically not eliminated [50,51].
The residual µ R -dependence needs to be kept under control to ensure that the perturbativity constraints will not get out of hand once the renormalization scale is changed. For that purpose, we perform one-loop RG running of dimensionless parameters and check for the stability of the vacuum position under loop corrections (cf. next point) at a renormalization scale other than µ R as well (specifics are given later in Appendix B 2).
For further convenience, we introduce some auxiliary quantities, repeating here the definitions of Eqs. (23) and (24) for completeness: where µ R is the starting "optimal" renormalization scale specific to every point in the parameter space, and µ R+ (µ R− ) denotes the renormalization scale for which the couplings hit a Landau pole when running upward (downward). The geometric average encodes, roughly speaking, how many orders of magnitude the scalar couplings can run up and down before encountering a Landau-pole-type singularity. Finally we define an optional penalty p 5 in order to test levels of robustness with respect to the RG running: where A 5 ∈ R + is a weight factor, H(x) is defined below Eq. (B10), andt thr is an acceptance threshold fort.
3. Stability of the vacuum position We start by defining a vector of dimensionful scalar parameters w = {µ 2 , ν 2 , τ 2 }. It is connected to the VEVs via one-loop stationarity conditions, which can be written in the form where w 0 denotes the tree-level part from Eqs. (12)- (14), and the function f ( w) represents one-loop corrections. Eq. (B15) is solved iteratively via taking w (0) = w 0 . The iterative procedure is stopped, and w is deemed to not have converged if where we set ζ = 0.3, and the result of the iterative procedure is labelled as w iter . The maximal number of iterations is chosen to be 30.
The penalization p 3 that assesses vacuum stability is the following: if λ 2 12 > 10 : else if w iter does not converge : else w iter converges: where A (i) 3 ∈ R + are weight factors, H(x) is defined below Eq. (B10), and the vector of scalar couplings is λ := {a 2 , a 0 , λ 0 , λ 2 , λ 4 , λ 4 , α, β 4 , β 4 , γ 2 , η 2 }. (B19) The penalization in Eq. (B18) is constructed so that each parameter point falls under exactly one of three conditions. The first condition and the associated A 3 penalize points for which the vector of couplings λ is sufficiently small, but w iter does not converge under the criterion in Eq. (B17); i.e., the iterative procedure is stopped prematurely. Analogously, the third condition and weight A 3 penalize points if a converged w iter differs too much from the tree-level w 0 and the initial approximation of the one-loop result w (1) . Note that the expressions are such that the penalization is smaller if a later condition applies, effectively ranking the criteria in descending order of importance.

Iterative pseudo-Goldstone masses
It turns out the computation of regularized one-loop effective masses for PGBs involves certain subtleties. As laid out in Table VI, the presence of particles with masses much below the GUT scale in the loop is handled by the replacement in Eq. (A64) in order to tame unphysical large-log contributions [52]. When computing mass corrections to PGBs, Eq. (A64) implies replacing log arguments involving the mass of lighter (e.g., intermediate-scale) particles with the PGB mass, leading to an equation of the form where C 1 and C 2 are numeric coefficients independent of the expression M 2 P GB . The tree-level masses of the PGBs are accidentally small and not very close to the chosen renormalization scale. This implies that the log contribution is additionally enhanced, and the C 1 term, usually dominated by the tree-level mass, is reduced. The relative importance of the C 2 term is thus increased, and the solution of Eq. (B20) can be far away from C 1 . This is problematic conceptually, since the C 2 term itself is merely an approximation arising from taming the logs. Note that the computation of heavy masses does not suffer from the same problem, despite an analogous equation, since their tree-level masses are large and the C 1 term there dominates.
In practice, we solve Eq. (B20) iteratively via starting with the value M 2 P GB,(0) as the tree-level mass and requiring that all iterations after i ≥ 1 are non-tachyonic. Note the absolute value in the argument of the log placed there to deal with the possible initial tachyonic instability. We perform a maximum of 30 iterations. We deem that convergence has been achieved and stop the process if the relative size of the shift between M 2 P GB,(i+1) and M 2 P GB,(i) of two consecutive steps is smaller than 10 −3 , a size corresponding to corrections at two-loop level. We denote the result of the iterative process as M 2 P GB,iter , while M 2 P GB, (1) is what we refer to in Appendix A 3 as the regularized one-loop effective mass.
Because of the issues discussed above, a point is considered valid not only if convergence in Eq. (B21) is achieved, but also that the C 2 -proportional log contribution is sufficiently small, ensuring the logtaming approximation works as intended and the computed regularized effective mass is close to the physical mass. Hence, we demand that the relative difference between the regularized one-loop effective PGB mass M 2 P GB, (1) and the iterative solution M 2 P GB,iter is less than 28 10%. We define the penalty p 4 associated with these considerations, applied now to multiple PGBs, by if M 2 P GB,iter does not converge : 4 arctan H(∆M 2 − 10 −1 ) +H(δM 2 iter − 10 −3 ) , else M 2 P GB,iter converges : with i max = 30.
2. Implementation of the numerical analysis: Step-by-step procedure Eventually, the criteria discussed above are used to find viable points, i.e., those which pass all the consistency checks and hence have zero penalization. Conceptually, the numerical procedure for assessing a single parameter point consists of two steps: (a) First, a point in the parameter space is selected. The VEVs and unified gauge coupling g are adjusted so that the resulting one-loop scalar spectrum is consistent with gauge unification requirements.
(b) Second, the penalization of the adjusted point is calculated. It is used for determining its viability.
The procedure of choosing new candidate points in the parameter space for assessment, i.e., a scan of the parameter space, follows the differential evolution algorithm [32] as the minimization algorithm of choice. We use a stochastic implementation with a mutation factor F randomly sampled from the interval [0.5, 2]. The scanning process first discovers viable points by minimizing the penalization function until it reaches zero, and then it explores the zero penalty region. In this way, datasets of viable points from Table III are  produced. Returning to part (a), the detailed steps are the following: 1. The input parameters for a point consist of where λ is defined in Eq. (B19) and g is the unified gauge coupling. These parameters are deemed to be the running values at the yet to be determined "optimal" renormalization scale µ R . For later convenience, we label the SO(10)-breaking VEV as ω max (ω R or ω BL in the regime ω BL → 0 or ω R → 0, respectively), and the subdominant induced VEV as ω min (ω BL or ω R in the ω BL → 0 or ω R → 0 limit, respectively).
2. We initially fix g = 0.5 and ω max = 10 16 GeV. The λ, σ, and ω min inputs are provided by the selection of the parameter point, for which we demand that they conform to Eq. (21).
3. With all the initial input parameters now set, the tree-level gauge and scalar spectra are calculated. 4. The "optimal" renormalization-scale square µ 2 R is determined as the arithmetic mean (over real degrees of freedom) of the heavy tree-level scalar-mass squares.
5. With µ R in hand, the one-loop scalar-mass spectrum for heavy and PGB fields is computed using the method described in Appendix A, replacing their tree-level values. The corrections to gauge fields and the intermediate-scale scalars are not calculated. We refer to this collection of our best available results for all gauge and scalar masses as the initial spectrum.
6. The initial spectrum is used for thresholds in the one-loop gauge evolution analysis. 29 The complete SM-gauge-coupling unification with the EW-scale boundary conditions [53] α −1 EM (M Z ) = 127.952 ± 0.009, is achieved as described in Appendix B 1 b. Technically, this is accomplished by three simultaneous actions: (i) rescaling of all dimensionful quantities {ω BL , ω R , σ, µ R } by a common factor, (ii) further adjustment 30 of σ and ω min , and (iii) optimizing the value of g. This procedure provides new values for g, µ R , and the three VEVs -defining the adjusted parameter point.

7.
For the adjusted parameter point, the tree-level spectrum for gauge and scalar fields is computed. Furthermore, the heavy and PGB masses are improved with the one-loop correction. This collection of the best available mass values is referred to as the updated spectrum.
In part (b), the remaining consistency criteria of Appendices B 1 a and B 1 c (tachyonicity and perturbativity) are imposed on the adjusted parameter point and the relevant penalizations are calculated. The detailed steps for part (b) are as follows: 1. The stability of the vacuum position is investigated -the penalization p 3 is computed. For a viable point, the one-loop corrections should not shift the position of the tree-level vacuum by more than 100%.
2. All the dimensionless parameters are RG run half-an-order of magnitude upward using one-loop beta functions (cf. Appendix C). Note that this running distance in scale turns out to be consistent with the spread of µ R values for viable points (see Sec. IV), implying that all viable points could be adjusted to a common scale for comparison if desired.
3. RG-run couplings are used to recheck the stability of the vacuum position -the penalization p 3 is acquired. All viable points satisfy t + > 0.5, and their shift of the tree-level vacuum position due to the one-loop corrections is still, even after running, not more than 100%.
4. The tachyonicity of the updated spectrum is inspected -the penalization p 1 is obtained.

5.
Global mass perturbativity is examined -the penalization p 2 is computed. For viability, the maximal one-loop scalar-mass correction is required to be smaller than the average of heavy tree-level scalar masses, i.e., ∆ ≤ 1. 29 For tachyonic masses, we take their absolute values. Such points are later penalized due to tachyonicity anyway. Note that the tachyonic property is generally preserved even after adjustment described in this step. 30 The VEVs σ and ωmin are adjusted in such a way that the ratio χ of Eq. (15) stays constant.
6. The iterative pseudo-Goldstone masses are investigated -the penalization p 4 is determined. For a viable point, the relative difference between the initial and final values in the iterative procedure for the PGB one-loop mass is constrained to less than 10%, i.e., ∆M 2 < 0.1.

(Optional)
Stricter RG-perturbativity criterion is imposed -the penalization p 5 is added. All viable points satisfyt >t thr , wheret thr is an acceptance threshold of a given search (see Table III in Sec. IV B).
8. The overall penalization of a point is the following: The relative significance of a particular criterion in the minimization algorithm, such as tachyonicity or perturbativity, can be tuned by changing the value of the A i weight factors present in penalizations p i . We chose the configuration A i = 1, which leads to p i ∈ [0, 1] for all i = 1, 2, 3, 4, 5. If the point has zero penalization, it satisfies all the imposed criteria and belongs to the set of viable points obtained in the parameter-space scan. Note that the shape of the viable part of the parameter space is independent of the values A i .

Appendix C: Beta functions for scalar parameters
The one-loop beta functions for scalar-potential parameters can be extracted from the field-dependent effective potential (A1) following the procedure outlined in [44].

General principles
Suppose that at tree level a scalar coupling λ can be written as a derivative of the scalar potential via for a suitable choice of scalar fields Φ i . The complete one-loop effective potential (A1) has to satisfy the Callan-Symanzik renormalization group equation where µ R is the running renormalization scale, and the j index runs over all scalar fields. Consequently, Equations (C1), (C3), and (A1) then imply where M 4 S (Φ) and M 4 G (Φ) are defined in Eqs. (A3) and (A4), β λ denotes the desired one-loop quartic scalar beta function satisfying and is the field-strength-dependent part.

Field-strength-dependent part
The β λ,F S part of Eq. (C4) is closely connected to the Φ-field anomalous dimension defined as where is the field-strength-renormalization factor of Φ in the MS renormalization scheme, and Σ M S Φ (p 2 ) is the corresponding self-energy. The momentum-squared-proportional part of Σ M S Φ (p 2 ) originates from the diagram which yields where C 2 (R) is the quadratic Casimir of the irreducible representation R that Φ belongs to. Hence, Combining Eqs. (C6) and (C10), the β λ,F S term in Eq. (C4) becomes  (10) into multiplets of SU(4) C × SU(2) L × U(1) R (the ω BL → 0, σ → 0 limit) in the left columns and SU(3) c × SU(2) L × SU(2) R × U(1) B−L (the ω R → 0, σ → 0 limit) in the right columns and their subsequent breaking into Standard Model representations after engaging non-zero σ; cf. Table II. The color scheme indicates whether a representation is real (blue), complex (black) or a complex conjugate (red) of some other multiplet within the same SO(10) representation. Since the 126 is a complex representation, all the multiplets that belong to it are also complex, while the 45 is real and thus consists of real multiplets and complex-conjugate pairs. The multiple copies of the same SM multiplets in the 126, e.g., the (1, 2, + 1 2 ) and (1, 2, − 1 2 ) weak doublets, therefore represent independent degrees of freedom and not just the complex-conjugate counterparts of each other.
(a) 45 of SO(10) Because of the relation in Eq. (15), the effect of the smallest VEV -either ω BL or ω R -on tachyonicity of the scalar spectrum should not be immediately dismissed. The distribution of χ presented in Fig. 16 shows a slight inclination for negative values in the ω BL → 0 case, while favouring positive values for ω R → 0. In other words, the first scenario prefers ω BL and ω R of different sign, while in the other limit, both VEVs are more often of the same sign. The eigenvalues of tree-level mass matrices for scalars are presented in both limits in Table IX. We show only their leading-order terms, which are proportional either to the SO(10)-breaking VEV (|ω R | or |ω BL |) or the intermediate-symmetry-breaking scale |σ|. When the dominant contributions to masses are M 2 GUT -proportional, the corresponding states belong to approximate SU(4 VIII: Tree-level spectrum of massive gauge bosons produced in SO(10) → SU(3) c × SU(2) L × U(1) Y breaking is listed for the most general case and for both perturbative limits of interest. In the last column, the changes of the corresponding one-loop β coefficients ∆a i for gauge-coupling running are gathered; cf. Eq. (B5). Note that massive gauge contributions implicitly include also those of their longitudinal components treated as adding the WGB scalar counterparts at the same scale.
• The masses of generic (heavy) states can be accurately estimated by their dominant M 2 GUT -proportional tree-level contributions. These are ω 2 R -or ω 2 BL -proportional in the ω BL → 0 and ω R → 0 limits, respectively.
• A separate category from the generic states are the pseudo-Goldstone bosons, whose tree-level masses are accidentally small, since their leading contributions are proportional to the suppressed value of a 2 . Therefore, their one-loop corrections can be of comparable size or even dominate. As in the generic case, the subleading |σ| 2 -proportional terms can still be safely neglected, and they are in fact not present at all at tree level for the triplet and octet. The states in this category consist of (8, 1, 0), (1, 3, 0), and (1, 1, 0) 3 , as well as (3, 1, − 2 3 ) 2 for ω BL → 0 and (1, 1, +1) 2 for ω R → 0.
As an interesting aside, note that the actual mass eigenvalues are not necessary to check for either nontachyonicity of the spectrum or gauge-coupling unification -in both cases, the quantities of interest can be directly extracted from mass matrices by computing their characteristic polynomial or from knowledge of their principal minors. Positive (semi-)definiteness can be examined by applying Sylvester's criterion on the matrices; cf. Sec. III A. This, for example, leads to conditions in Eqs. (25) and (26) for the a 2 → 0 limit. Alternatively, alternating signs of coefficients in the characteristic polynomial also imply positive definiteness of the spectrum. For positive semi-definite matrices with n zero-eigenvalues, the n lowest-order coefficients in the characteristic polynomial vanish. Unification constraints at one loop, on the other hand, require only knowledge of the product of non-zero eigenvalues n R k=1 m 2 R k for every SM representation R with multiplicity n R . This number can be extracted as the absolute value of the lowest degree coefficient in the characteristic polynomial that does not vanish. The reason that this product suffices is that all n R eigenvalues of multiplets in the representation R contribute with the same β coefficient ∆a to Eq. (B7).

One-loop pseudo-Goldstone masses
An accurate assessment of the PGB masses in the small a 2 regime requires also the computation of their one-loop corrections. In addition to the recurring octet, triplet, and singlet pseudo-Goldstone fields, the two limits ω BL → 0 and ω R → 0 provide another PGB in the form of the (3, 1, − 2 3 ) 2 or (1, 1, +1) 2 state, respectively.
An interesting observation is that the weak-triplet and SM-singlet PGBs receive the same one-loop gauge contributions to their masses in the ω R → 0 limit. Note that the two mentioned states belong to different representations of the intermediate left-right symmetry, namely, to (1, 3, 1, 0) and (1, 1, 3, 0), respectively. However, these intermediate-stage representations are connected by D-parity (L ↔ R exchange), which corresponds to a transformation ω BL → −ω BL at the VEV level. In the ω R → 0 case, the dominant oneloop gauge contributions to masses must be proportional to g 4 ω 2 BL and are hence the same for pairs of states that exchange under D-parity. The one-loop contributions from scalars, on the other hand, are not parity invariant due to the implicit presence of massive scalar-potential parameters µ, ν, and τ obtained by solving stationarity conditions, where a change of sign in ω BL does have an effect. The dominant contributions to tree-level scalar masses in both perturbative scenarios and the associated contributions ∆a i to one-loop β coefficients for gauge-coupling running.
Mass ω BL → 0 ω R → 0 (∆a 3 , ∆a 2 , ∆a 1 )  The one-loop contributions to masses of PGB states in the a 2 → 0, γ 2 → 0 limit 31 are presented in Table X. Although this limit is apparently unphysical since the numerical scans show that |γ 2 | 0.1 for valid points (cf. Fig. 9), it is well suited for obtaining analytical expressions. For both the ω BL → 0 and 31 Note that in this limit, the tree-level PGB masses exactly vanish (see Table IX), so the one-loop corrections actually represent their entire masses at the one-loop level. ω R → 0 scenario, we show only the dominant, SO(10)-breaking terms proportional to the largest VEV, while neglecting those with σ. The qualitative features of PGBs discussed in the beginning of this section are confirmed by the explicit formulae in Table X, at least in the given limit. For more general formulae with an arbitrary ω BL and ω R hierarchy, but still in the a 2 → 0, γ 2 → 0, σ → 0 limit, the interested reader is referred to [26,43]. The dominant one-loop contributions to masses of PGB scalars in the a 2 → 0, γ 2 → 0 limit for both perturbative scenarios and containing polynomial as well as logarithmic terms with gauge bosons and scalar fields running in the loop.