Flavor Symmetries in the Yukawa Sector of Non-Supersymmetric SO(10): Numerical Fits Using Renormalization Group Running

We consider a class of $\text{SO}(10)$ models with flavor symmetries in the Yukawa sector and investigate their viability by performing numerical fits to the fermion masses and mixing parameters. The fitting procedure involves a top-down approach in which we solve the renormalization group equations from the scale of grand unification down to the electroweak scale. This allows the intermediate scale right-handed neutrinos and scalar triplet, involved in the type I and II seesaw mechanisms, to be integrated out at their corresponding mass scales, leading to a correct renormalization group running. The result is that, of the 14 models considered, only two are able to fit the known data well. Both these two models correspond to $\mathbb{Z}_2$ symmetries. In addition to being able to fit the fermion masses and mixing parameters, they provide predictions for the sum of light neutrino masses and the effective neutrinoless double beta decay mass parameter, which are both within current observational bounds.


I. INTRODUCTION
One way of extending the Standard Model (SM) is to assume that, at some high energy scale, a Grand Unified Theory (GUT) [1] provides a more fundamental description of particle physics.
Particularly, models based on the gauge group SO(10) [2] are successful in accommodating all SM fermions in one representation per generation and providing an origin for several of the features in the SM. It is also a useful class of models in which to explore the outstanding questions of the SM, such as neutrino masses through the type I and II seesaw mechanisms [3][4][5][6][7][8][9][10].
This work is structured as follows. In Sec. II, we briefly review the models and discuss their parametrization. Then, in Sec. III, we describe the fitting procedure, and in Sec. IV, we present the results of the fits. Finally, in Sec. V, we summarize and conclude. A full list of the 14 models is presented in App. A.

A. Models
In this work, we consider a setup in which the SO(10) Yukawa sector contains scalar fields belonging to the 10 H , 126 H , and 120 H representations. As such, the Lagrangian is given by where the fermions are placed in the 16 F representations and Y i are 3 × 3 matrices of Yukawa couplings in flavor space. Due to the structure of the SO(10) coupling, the Yukawa matrices are such that Y 10 and Y 126 are symmetric, while Y 120 is antisymmetric. The different models that we investigate correspond to different textures of zeros of these three matrices. These textures originate in symmetries under phase transformations of the three fermion generations and the scalars in the Yukawa sector. In this framework, the allowed symmetries are Z 2 , Z 3 , Z 4 , U(1), and Z 2 × Z 2 . The complete derivation and discussion of these symmetries are presented in Refs. [18,38]. Here, we reiterate the main points.
The symmetries are determined by the possible rephasings of each of the three generations independently, as well the three scalar multiplets. Any specific Yukawa matrix texture will impose constraints on these phases.  A : where the symbol "×" denotes a general non-zero matrix element.
The more non-zero matrix elements that a model has, the more constrained the symmetry is.
Similarly, many zeroes imply a large degree of freedom in the symmetry transformations. Some of the models can be derived from other models by setting some of the elements to zero. One can, for example, derive model A 1 from A in this way. All such submodels are included only if they give rise to different symmetry groups. For example, setting the (1, 2) element of Y 10 to zero in model A does not affect the symmetry, so that is not considered to be a separate model.
Although this will, in general, affect the RG running of the fermion observables through the contributions of the gauge couplings to the RGEs, we assume that this change is negligible to the ability to fit the models. This assumption is based on the fact that varying the unification scale by an order of magnitude has a small effect on the goodness of fit (see Fig. 3 [57]. The matching conditions are then given by [16,20,21,23] v where Y u , Y d , Y D , and Y are the Yukawa coupling matrices of the up-type quarks, down-type quarks, neutrinos, and charged leptons, respectively. In addition, M R is the mass matrix of the right-handed neutrinos and M L is the contribution from type II seesaw to the light neutrino mass matrix. The vevs v u,d 10,126,120 are the vevs of the SU(2) L doublets contained in the 10 H , 126 H , and 120 H , while v R is the vev of the SM singlet in 126 H that gives mass to the heavy right-handed neutrinos and v L is the vev of the SU(2) L triplet involved in the type II seesaw contribution to the neutrino masses.
For convenience, we reparametrize using H = /v SM and write the Yukawa and neutrino mass matrices as Since the matrices H, F , and G are proportional to the Yukawa matrices Y 10 , Y 126 , and Y 120 , respectively, they have, of course, the same textures.
The parameters other than the Yukawa matrix elements that are to be sampled are then r, s, t u , t D , t , r R , and r L . Of these, we assume r and r R to be real, while the others are complex, giving a total of twelve parameters related to the vevs. Furthermore, we need v d 126 and the mass M ∆ of the SU(2) L triplet for the matching to the effective neutrino mass matrix, thus adding two parameters.
For more details on the effective neutrino mass matrix, see Sec. III.
To count the number of parameters in the matrices H, F , and G, we first note that they, in general, have complex matrix elements and that the (anti)symmetry of them restrict the number of free parameters. One can, without loss of generality, transform one of H or F to a real and diagonal matrix. This same transformation must also be applied to the other two matrices, so it should only be performed if it decreases the total number of parameters in the model.
Since the vevs of the SU(2) L multiplets all contribute to the electroweak gauge boson masses after electroweak symmetry breaking, they need to add up in quadrature to the SM Higgs vev.
The parametrization used leaves some freedom in this constraint, since not all vevs are sampled.
Further imposing a perturbativity bound of 4π on the Yukawa couplings, we have the inequality where max |H| 2 denotes the largest absolute value of the elements of H and similarly for G. The factor of 2 in front of |r L | 2 is due to how the triplet couples to the gauge bosons. This inequality needs to be checked after each fit to make sure that the found parameter values are acceptable.
Furthermore, there are experimental bounds on the parameters related to the SU(2) L triplet: v L 1 GeV and M ∆ 1 TeV [58,59]. Finally, the vev v R should be set to M GUT . However, we allow this parameter to vary by around an order of magnitude from M GUT , which can be motivated by the existence of threshold corrections and higher-dimensional operators [60][61][62][63][64].

III. FITTING PROCEDURE
When fitting the models to data, we use the known values of the fermion masses and mixing parameters at the electroweak scale M Z and perform top-down RG running for each sampled set of parameters. This is in contrast to the procedure used in Ref. [18], where the fits were performed with values that had been evolved to M GUT using bottom-up RG running prior to fitting. Our approach of solving the RGEs for the sampled parameters has the benefit of correctly taking into account the mass thresholds of the heavy right-handed neutrinos and the scalar triplet for type I and II seesaw mechanisms, respectively. That is not possible when solving the RGEs before the fitting procedure, since the masses of these heavy fields are sampled and are thus unknown prior to the sampling process. However, the approach used in this work significantly increases the computational resources required, since the RGEs are solved for every sampled point in the parameter space. The computational time required is increased by a factor of about 1000 compared to the approximation of extrapolating the data to M GUT before fitting.
For each sampled point in the parameter space of the SO(10) Yukawa matrices and vevs, we transform the parameters to the SM Yukawa matrices and solve the RGEs at one-loop order [16,[65][66][67][68][69][70][71][72] down to M Z . When the RGE solver encounters a threshold corresponding to a heavy righthanded neutrino or the scalar triplet (involved in the type I and II seesaw mechanisms, respectively), the corresponding heavy field is integrated out and all relevant parameters are transformed to the effective field theory, before the RGE solver proceeds. A more detailed description of this procedure may be found in Ref. [28], on which our code is based.
The effective neutrino mass matrix κ is formed by starting with a matrix of zeroes. At an energy scale corresponding to the mass M i of a right-handed neutrino, it is updated according to [69,70] where Y (i) D is the the ith row of Y D in a basis in which M R is diagonal. In this basis, the ith row of Y D is removed as well as the ith column and row of M R . At the energy scale of the scalar triplet mass, κ is similarly updated according to [72]  charged lepton masses, and Higgs quartic coupling are adopted from Ref. [73], neutrino masses and leptonic mixing angles from the July 2020 version of NuFIT [74], and CKM parameters from the Summer 2019 version of CKMFitter [75]. An asterisk denotes an error that has been enlarged to 5 % of the corresponding central value.
and Y ∆ is then removed from the parameters in the RG running.
The data that the parameters are fitted to are displayed in Tab. II. All quark and charged lepton masses as well as the Higgs quartic coupling are adopted from Ref. [73]. Neutrino masses and leptonic mixing angles are the normal neutrino mass ordering 1 values from the July 2020 version of NuFIT [74]. Finally, CKM parameters are taken from the Summer 2019 version of CKMFitter [75].
Some observables (e.g. the charged lepton masses) are known to an extremely high precision, which makes the fit practically very difficult. To remedy this, we have artificially enlarged the errors to be at least 5 % of the central values. This applies to all observables except for m u , m d , and m s .
After performing the RG running down to M Z and computing the values of the 19 SM observables, we form the χ 2 function defined as where X i is the prediction from the sampled point andx i and σ i are the central value and error, respectively, as given in Tab. II. The fitting procedure involves minimizing the χ 2 function using the Diver code from the ScannerBit package [76], which implements a differential evolution algorithm and can be parallelized on a cluster utilizing several CPU cores. Although one can never guarantee a global minimum to be been found, we increase our confidence in the result by running the Diver code twice. More than this was not feasible due to the constraint on the number of core-hours available for computations. After convergence, we look for further improvements of the candidate best-fit point using a Nelder-Mead simplex algorithm [77].
For comparison, we also perform a fit to data at M GUT without the procedure for RG running of the sampled parameters discussed above. To find the data at M GUT , we use the code SMDR [78] to perform bottom-up RG running of the data presented in Ref. [73]. needs to be very close to zero at M GUT . Hence, we set it as such and, in order to make sure that it stays positive for the sake of vacuum stability [79], we perturb it a little after the fit without any noticeable effect on the results.

IV. RESULTS
The results of the fitting procedure for all models are shown in Tab    bounds [82]. On the other hand, the predicted values of δ CP in these two models deviate from the preferred values of global fits [74], as also found in Ref. [28]. However, one should keep in mind that these are not strict predictions of the models and that it is still possible that a fit that includes δ CP is able to accommodate the favored value.
The parameter values that produce the fits in Tab. IV for model A are In both of these models, the dominant contribution to the neutrino masses stems from the type I seesaw mechanism. By comparing the maximum singular values of the contributions in Eqs. (16) and (17) evaluated from the parameters at M GUT , we observe that the type I seesaw mechanism dominates by a factor of about 400 in model A and a factor of about 13 in model B. These numbers may, of course, change due to RG running, but the general conclusion should still hold.

V. SUMMARY AND CONCLUSIONS
In this work, we have analyzed a class of non-supersymmetric SO(10) models with flavor symmetries, following the framework described in Ref. [18]. In the numerical fits, we have solved the RGEs for the sampled parameters from M GUT down to M Z in order to fit to the known parameters at M Z . This approach allows one to properly integrate out the heavy right-handed neutrinos and the scalar triplet, thereby taking into account the RG running of the neutrino parameters in an effective field theory framework, as opposed to approximate approaches in which the known parameters are extrapolated up to M GUT .
Out of the 14 models, only two have been found to accommodate the known values of fermion masses and mixing parameters. These are the ones labeled as model A and model B which both possess a Z 2 symmetry, in agreement with previous results [18] that were found with supersymmetric models and bottom-up RG running. These two models are the ones with the largest number of free parameters, and are hence expected to provide the best fits. The difference between Ref. [18] and this work is that we considered non-supersymmetric models, that we used updated data, and most importantly, that we solved the RGEs of each sampled set of parameters from M GUT down to M Z .
We have only considered the symmetries of the Yukawa potential in this work, and not the effect of the symmetries on the scalar potential. However, imposing a Z 2 symmetry on the scalar potential should not present any major difficulties, and will, in general, mean that the 45 representation involved in the symmetry breaking (not analyzed in detail here) will also inherit this symmetry. It remains to construct a complete model with these flavor symmetries that accounts for symmetry breaking and gauge coupling unification. Finally, it would be interesting to investigate whether the disagreement between the predicted value of the CP-violating phase and that from global fits is a true prediction of the models or if it is possible to accommodate it by including it as one of the parameters in the fit.

ACKNOWLEDGMENTS
We would like to thank Luís Lavoura for useful discussions and Mattias Blennow for supply- In this appendix, we list the Yukawa matrix textures corresponding to the 14 models in Tab. I as Eqs. (A1)-(A14). A "×" symbol signifies a general non-zero matrix element. For the symmetry transformations and groups corresponding to each model, see Ref. [18].