Confronting SUSY SO(10) with updated Lattice and Neutrino Data

We present an updated fit of supersymmetric SO(10) models to quark and lepton masses and mixing parameters. Including latest results from lattice QCD determinations of quark masses and neutrino oscillation data, we show that fits neglecting supersymmetric threshold corrections are strongly disfavoured in our setup. Only when we include these corrections we find good fit points. We present $\chi^2$-profiles for the threshold parameters, which show that in our setup the thresholds related to the third generation of fermions exhibit two rather narrow minima.


Introduction
Since the advent of Maxwell's theory of electromagnetism it is a common dream in physics to unify all interactions into one single theory. This dream is particularly persistent in particle physics where Grand Unified Theories (GUTs) -which unify three of the fundamental forces of nature -have been a major guiding principle in the last couple of decades. To be precise in this paper we will focus on GUTs based on the SO(10) gauge group first proposed in the 1970s [1,2]. This choice contains all the essential features interpreted as hints towards unification, like charge quantization, anomaly cancellation or smallness of neutrino masses.
SO (10) GUTs work especially well in the context of low energy supersymmetry (SUSY) [3] which helps significantly to unify the gauge couplings, see, e.g., [4]. In fact, supersymmetry is part of the dream of unification as it unifies the concepts of bosons and fermions or the concepts of matter and forces. It also provides an elegant mechanism to stabilize mass hierarchies and provides a dark matter candidate on top.
In this paper though we will not discuss all features and aspects of SUSY SO(10) models in general. We focus on the Yukawa sector, which is interesting because this sector is constrained by low energy observables, namely the observed fermion masses and mixing parameters. There has been some tremendous progress in the last couple of years in particular for light quark masses from lattice computations and neutrino masses and mixing from oscillation experiments. We will use this updated information to provide an updated fit to a minimal SUSY SO(10) Yukawa sector.
There is a long history of fitting SO (10) to fermion masses, see, e.g., . Some works focussed on particular SO(10) models, while others focussed on a general SO(10) Yukawa sector only. These efforts have always been a difficult task due to the large number of parameters and observables. For that reason many groups had to resort to some simplifications or estimates. For instance, the recent fits by Dueck and Rodejohann [47] neglect in the SUSY case threshold corrections which are known to give sizeable, important corrections, e.g., Refs. [3,7,11,48,49]. One of the main improvements of our work compared to previous studies is to include these threshold corrections in terms of three parameters and to provide χ 2 -profiles for them for the first time, to our knowledge.
Our paper is organised as follows: In Section 2 we describe in detail how we match the SO(10) Yukawa couplings at the high scale to the fermion masses and mixing parameters at low energies including the fit procedure. In Section 3 we discuss our results. First we present the global minima of our fits and the pulls to identify which observables are potentially driving some tensions of the data with our setup. Then we present the χ 2 -profiles for the SUSY threshold corrections before we summarise and conclude in Section 4. We provide additional information in two appendices. Appendix A gives additional numerical values for Standard Model (SM) observables, and Appendix B the values of the GUT scale fit parameters at the global and local minima.

Fitting GUT Scale Yukawa Couplings to Fermion Masses
Besides the unification of the gauge sector, SO(10) GUTs also predict the unification of the Yukawa sector by arranging all matter fields in the spinorial 16 representation. Of course, these constraints hold only at the GUT scale and effects of GUT breaking and of the renormalisation group equations (RGEs) will alter the relations among the Yukawa couplings.
Below the GUT scale, important physics happens around the seesaw scale where the heavy right-handed neutrinos are integrated out one after another. Moreover, SUSY breaking effects have to be taken into account in supersymmetric GUTs. Especially for large values of tan β, threshold corrections from SUSY particles have sizeable effects [3]. In our analysis, we choose the matching scale between the SM and its minimal supersymmetric extension (MSSM) to be M SUSY = 1 TeV. At this scale we define as well the χ 2 -function to fit the GUT parameters to the SM Yukawa couplings. Therefore we first need to evolve the fermion masses and mixing parameters to M SUSY using the SM β-functions. We describe the procedure in greater detail in the following.

Yukawa Couplings in SO(10) GUTs
In this work we restrict ourselves to renormalisable SO (10) GUTs. This has a stronger predictive power than SU(5) or non-renormalisable GUTs. Furthermore we assume the standard embedding of the matter fields into the spinorial 16 of SO (10). This restricts the Higgs representations relevant for fermion masses to a 10, 120 and/or 126 of SO (10).
In SUSY SO(10), the most general renormalisable superpotential describing Yukawa interactions is then given by are symmetric while Y ij 120 is antisymmetric in the flavour indices i and j. Without loss of generality we apply an unphysical flavour rotation to choose Y 10 real and diagonal.
After breaking SO(10) at the high scale M GUT = 1.353 × 10 16 GeV, given by gauge coupling unification 1 , the Yukawa matrices are matched to the MSSM Yukawa matrices Y x , x = u, d, e, ν, the Majorana mass matrix for the right handed neutrinos M and the Wilson coefficient of the Weinberg operator κ, see for example Ref. [47] The parameters s, t u , t d , t ν , t e correspond to the mixing of GUT scale Higgs doublets into the MSSM Higgs doublets H u , H d whereas r, r R , r L correspond to the vacuum expectation values (vevs) of the GUT scale Higgs fields and can be chosen real without loss of generality. A minimal realistic choice omits the 120 representation (i.e., the parameters t u , t d , t ν , t e ) and r L . Previous fits to fermion observables have shown that the 120 does not significantly improve the fit to fermion masses in the SM, e.g., Ref. [47], so we will not include it here (Y 120 = 0). Furthermore, in Ref. [50] it was argued that while a generic fit to fermion masses prefers a mix of type-I seesaw (parametrised by r R ) and type-II seesaw (parametrised by r L ), a fit which includes the GUT potential favours scenarios where type-I is dominant over type-II. Hence, we will assume r L = 0. Thus we end up with a set of 19 GUT scale parameters: Y 10 (three parameters), Y 126 (twelve parameters), r R (one parameter), r (one parameter) and s (two parameters).
Before we discuss our numerical implementation and the results we want to give a few comments on Yukawa couplings from higher-dimensional operators. In this study we explicitly assume that no higher-dimensional operators correct the Yukawa coupling relations beyond the ones present at renormalisable level. Especially, in the context of supergravity one might expect higher-dimensional operators to be present in the superpotential suppressed by powers of M GUT /M Planck = O(0.01). Such higher-dimensional operators could alter the Yukawa relations if they have a non-trivial gauge structure.
It is an open question if quantum gravity introduces GUT non-singlets at the Planck scale which couple to ordinary matter. The presence of non-singlets would imply that quantum gravity knows something about flavour. If it would be flavour-blind all Yukawa couplings should be at least of order M GUT /M Planck = O(0.01). However, this is not the case unless there are some large cancellations at work. Hence, we interpret the smallness of the electron mass as an indication that these higher-dimensional operators are negligibly small and we do not consider them here. Nevertheless, in a supersymmetric GUT theory of flavour they can be the dominant operators and very sensible, for a recent review, see, e.g., Ref. [51].

RGE Running from the GUT Scale to the SUSY Scale
The values of the Yukawa couplings at different energy scales are related by the renormalisation group equations via with the beta function being defined as β Y (µ) = dY /d log µ. At one-and two-loop level they are well-known for both the SM and the MSSM with additional right-handed neutrinos, see, e.g., Ref. [52]. We also include the running of the Weinberg operator at the one-loop level. For every point in the GUT parameter space the RGEs have to be solved, i.e. numerically integrated. For a correct treatment of the seesaw mechanism, we consecutively integrate out the right-handed neutrinos as described in Ref. [53] and implemented in the Mathematica package REAP presented therein. REAP provides a useful tool for cross checking our calculations. We implement the numerical procedure in C++ using the template libraries odeint [54] and eigen3 [55]. This speeds up the calculation of the RGE running by up to a factor of 100 compared to the Mathematica version of REAP on our test desktop machine (Intel i5-4590 CPU, 3.30 GHz).

Fermion masses at the SUSY scale
Our low scale input for the quark masses and α s is given in Table 1. We perform the QCD RGE evolution to M Z = 91.1876 GeV [57] and the matching from 4 to 5 as well as from 5 to 6 active flavours with the Mathematica package RunDec [58] to four-loop accuracy. Our results are shown in Table 2. Errors from variation of the bottom and top quark scales are          negligible compared to the experimental uncertainties. The Yukawa couplings at M Z can be obtained from y = m/v, with v = 174.104 GeV. As additional input at M Z , we use the values for the gauge couplings g 1 and g 2 as well as the values for the lepton Yukawa couplings from Ref. [59]. Regarding the input for the CKM matrix, we use the ICHEP2016 update from CKMfitter [60]. Corresponding results from UTfit are in good agreement [61]. In order to perform the RGE evolution from M Z to M SUSY , we use the two-loop RGEs for the SM Yukawas and gauge couplings [62]. For solving the RGEs at two-loop also the Higgs quartic coupling is needed which we get from the PDG average of the Higgs mass measurement m Higgs = 125.09 ± 0.24 GeV [57]. For the off-diagonal elements of the Yukawa matrices we use the standard parametrisation of the PDG. Further details can be found in Appendix A.
We treat the uncertainties of the experimental data as Gaussian and symmetrise them arithmetically where necessary. For the propagation of uncertainties we sample normally distributed random numbers at M Z . Each sample point is then evolved to the respective energy scale. There, we fit a normal distribution to the set of sample points. Thus, also nonlinear effects, especially the influence of the top quark Yukawa coupling, and interdependencies in the RGE running are taken into account. This has the effect that the relative errors are larger at higher energy scales due to the error of the top quark mass measurement.
For our global fit we set M SUSY = 1 TeV. However, for reference, in Table 3 we additionally list also the SM Yukawa and gauge couplings in the MS scheme at M Z , 3 TeV and 10 TeV for reference. The same results converted to the DR scheme can be found in Table 4.

SM vs. MSSM Yukawa couplings
At the SUSY scale (which we set to 1 TeV) two important things happen at the same time.
First of all, loop calculations in the SM are usually done in the MS scheme which is unsuitable for supersymmetry. For supersymmetry the DR scheme is preferred. We will use two-loop RGEs so we have to perform the matching at the one-loop level. The relevant formulae for gauge couplings and Yukawa couplings are [63]  with C(G) = N and C(r) = (N 2 − 1)/2N for SU (N ). In the MSSM the matching conditions for the Yukawa matrices then read The second thing is that the Yukawa couplings in the SM and the MSSM are not the same, not even on tree-level since the MSSM is a two Higgs doublet model. This introduces a dependence on tan β, the ratio of the two Higgs vevs. Beyond that also finite one-loop corrections have to be taken into account in the matching. In particular there are certain pieces which are enhanced by tan β [3] and thus can easily make changes of O(10%) or even larger on the Yukawa couplings.
The approach we will take here is well documented, for instance, in Refs. [64][65][66], see also Ref. [59], so that we will not go into much detail here. Note that we include the tan β enhanced parts only. This allows a simple parametrisation of the SUSY threshold corrections in terms of three parameters, b , q and l , only using some rather mild, plausible assumptions on the SUSY breaking parameters. To be more precise we assume that the squark and slepton mass matrices are very close to being proportional to the unit matrix and that the trilinear couplings are hierarchical and dominated by the 3-3 element.
For the up-type quarks we can use the tree-level matching relation for SM and MSSM Yukawa couplings in the DR scheme, since their threshold corrections are proportional to cot β 1. For the Yukawa couplings of the charged leptons and down-type quarks there are tan β enhanced threshold corrections in the matching formulas where we have used the conventions as in Ref. [66] but we defined here b ≡ A + q . Note that Y e MSSM and Y d MSSM are not diagonal in our approach and we determine the corresponding low energy mixing parameters from these matrices after the matching.
In the following, if we do not specify SM or MSSM we refer to MSSM DR quantities.

Fitting Procedure
For every point in the GUT parameter space that is scanned over we extract the fermion observables O theo where σ i are the standard deviations of the experimental observables. The first term in Eq. (2.17) includes the masses of quarks (6 observables) and charged leptons (3 observables) as well as the CKM parameters (4 observables) and the mass squared differences of the light neutrinos (2 observables). For these we use a Gaussian error. The The other Gaussian observables are listed in Table 4 in the column '1 TeV'. For reference we also show the corresponding results in the MS scheme in Table 3.
The Gaussian treatment of the errors would be yet inadequate for the PMNS parameters. Therefore, for these we use the corresponding χ 2 -profiles from NuFit 3.0 [67], which are represented in the second term in Eq. (2.17). Our implementation also includes a check for inverted or normal mass ordering to ensure the right distribution is used. In particular we also include the second (higher) χ 2 minimum of θ l 23 . In that way, also the CP phase δ l CP can be included although it is not yet measured directly. Note that neutrino observables are treated here as tree-level observables and consequently by definition do not take part in the running under the renormalisation group. Considering this, we employ the results of NuFit 3.0 directly at M SUSY = 1 TeV. Altogether, we have 19 observables and 22 parameters in the fit, including the three threshold parameters b , q and l . As these stem from 1-loop corrections, we allow for them the following ranges: see, for example, Refs. [59,64,66,68]. The corrections to quarks are generally larger since they receive SUSY QCD corrections. Furthermore, b receives another correction from the potentially large stop trilinear SUSY breaking coupling. One might worry here that the corrections become non-perturbative, since our scan allows, e.g., | b tan β| > 1. But this is not the case since there are no higher order corrections O(( b tan β) n ) with n ≥ 2, see, e.g., [69]. In order to perform the fits we link our C++ code to the Sbplx/Subplex [70,71] and ISRES algorithms [70,72,73] implemented in the NLopt 2.4.2 library [70]. The extensive numerical fits are performed on the TTP computing cluster. As these are highly time consuming we do not include two minor online-updates of the NuFit results which appeared in the meantime.

Results of the SUSY SO(10) Fit to Flavour Data
In the following, we present our results for the χ 2 fit of SUSY SO(10) models to the quark and lepton flavour data given at 1 TeV as described in Section 2.5. First, we present the global minima, which show that threshold corrections are essential to get a good fit in our setup. Afterwards, we present the χ 2 -profiles for the SUSY threshold parameters.

Comparison of the Global Minima
We search for the global minima for three different values of tan β = 10, 38 and 50 as well as both with and without SUSY threshold corrections so that we have six global minima in total. The results of our global minimisation can be found in Table 5. The table shows the minimal χ 2 from our scan, the pulls of the SM observables as well as the predictions for the light neutrino masses and the needed value of the threshold corrections in case they are included. We define the pull for an observable O i with a Gaussian error by For the PMNS parameters we use the χ 2 -distribution and the best-fit values of NuFit O best−fit i to define Additionally, these pulls are also visualised in the bar charts of Fig. 1. The corresponding detailed fit parameters at the GUT scale are given in Appendix B.
Without including SUSY threshold corrections we find a poor description of the data by our SO(10) setup, with best fit points having χ 2 > 70. It turns out that this is especially triggered by the high precision of quark mass determinations on the lattice, in particular, for the down-type quark Yukawa coupling the tension is for all three values of tan β more than 6σ. There is another major tension in the atmospheric mixing angle where the modulus of the pull is always larger than 2 and even larger than 5 for tan β = 10.
Note that the fit without threshold corrections has as many parameters as observables, and including threshold corrections we have three parameters more than observables so that one might expect to see a perfect fit. However, the dependence of the observables on the parameters is highly non-linear, so that the expectation of a perfect fit in this case is misleading. Indeed, the observables are singular values of complex matrices with a strong hierarchical structure and are coupled via non-linear differential equations. This also leads to the fact that not every GUT scale fit parameter enters the same low-scale observables. For instance, the 1-1 elements of the Yukawa matrices will practically not affect the masses of the third generation fermions, whereas the 3-3 elements affect all other couplings and together with the 2-3 elements determine the mixing between second and third generation.
Having more parameters than observables the statistical interpretation of χ 2 is not straight forward. Nevertheless, we can still use the χ 2 information in order to compare different fit scenarios. First of all, our results clearly prefer a large tan β which is a well-known result. Our minimal SO(10) setup is hardly viable without threshold corrections and for tan β = 10 even including that corrections, see Table 5. Therefore, we conclude that sizeable, i.e., percent order SUSY threshold corrections and large values of tan β are needed in order to be in accordance with data. Note that for tan β = 10 at the global minimum we have q = 0.05 and l = −0.03 which are both at the edges of the respective allowed ranges for these parameters. That means the fit would prefer larger values of q and l or larger tan β.
Comparing the pulls of the different observables, it is remarkable that besides y d we find large discrepancies in the neutrino observables. This makes a thorough treatment of the right-handed neutrinos necessary: We treat them in an effective field theory picture where each of them is integrated out at its mass scale, as described in detail in Ref. [53].
Within our setup we always find the neutrinos to have a normal mass hierarchy with the heaviest neutrino having a mass of about 50 meV. The sum of light neutrino masses is below the cosmological bound m ν 0.23 eV/c 2 [74] for all global minima. The same holds obviously for the effective β-decay mass where the current upper bound is m ν,β < 2.05 eV [75]. Our best-fit results are also far below the sensitivity of KATRIN [76] which is expected to improve this bound by an order of magnitude.

Likelihood Profiles for Threshold Corrections
Given the importance of the SUSY threshold corrections for the goodness of fit of SO(10) models to fermion masses, see Section 3.1, it is an interesting question how sensitive the data is to the magnitude of the threshold corrections. This has in particular non-trivial implications for potential concrete SUSY spectra.
In order to answer this question, we determine the χ 2profiles of b , q and l , which are shown in Fig. 2. On the left-hand side we show the results over the full range that we allow, see Eqs. (2.20)-(2.22), while on the right-hand side we zoom into the interesting regions around the minima for tan β = 38, 50. Around the minima we increased the number of points to exclude possible numerical artefacts and to obtain a better resolution of the profile. Keeping in mind that the threshold corrections get multiplied by tan β, see Eqs. (2.15) and (2.16), we first like to discuss the case of tan β = 10. In this case, the χ 2 -profiles have in general less features than for tan β = 38, 50. Still, there is a distinct minimum in the b -profile and a clear preference for smaller values of l and larger values of q . Actually, the global minimum lies on the border of the allowed ranges for q and l as we mentioned before.
Interestingly, for tan β = 38 and tan β = 50 we see two competing narrow minima in the χ 2 -profile of b . In fact, our fit fixes the viable values for b to narrow intervals which can be an interesting input for SUSY model building and phenomenology. Furthermore, we obtain that b < 0 whereas q > 0. This means that the fit favours to decrease the bottom Yukawa coupling and to increase the Yukawa coupling for the strange quark and down quark.
As a curiosity, we note that in our scans the deeper of the two minima for tan β = 38 is the one with smaller b while for tan β = 50 it is the other way around. In particular for tan β = 50 the minima are almost degenerate with χ 2 = 3.7 and χ 2 = 4.7. It might be an interesting observation for SUSY phenomenology, that the two minima competing correspond to either b tan β < −1 or −1 < b tan β < 0, respectively.
Altogether, it is interesting to note that the largest tensions, cf. Figs. 1 and 3, are in the down quark mass, up quark mass, θ q 13 and the leptonic CP phase δ l CP . Especially, for the light quark masses and the leptonic CP phase we expect a significant reduction in the uncertainties in the next decade which might help to disfavour one of the two minima or even both minima. The importance of the leptonic CP phase for SO(10) fits was also noticed already, e.g., in Ref. [45]. This again shows that, although we have more parameters than observables, we can still make qualitative statements.

Summary and Conclusions
SUSY SO(10) is one of the best motivated extensions of the SM and due to its aesthetics it served as a guiding principle in particle physics for quite a while. Yet, as no signs of SUSY or proton decay have been seen in experimental searches, it still remains uncertain if SUSY GUTs are just a dream or if they have indeed something to do with nature.
Fermion masses and mixing can give a partial answer to this question. In this work we have fitted a minimal SO(10) Yukawa sector to recent data for fermion masses and mixing. Regarding the available data, compared to most previous studies there is a quite significant reduction of uncertainties for light quark masses and neutrino mixing parameters. These changes are so drastic that minimal SUSY SO(10) without the inclusion of sizeable SUSY threshold corrections is highly disfavoured nowadays.
Although it was known for quite a while already that these corrections are generically large and important [3] they were usually neglected in fits. Even in our rather minimal setup the Yukawa sector has already 19 parameters which we fit to 19 observables. SUSY thresholds plausibly add at least another three parameters which we have included here increasing the computational needs quite drastically.
In fact, we did not only include them but also provide for them χ 2 -profiles for the first time to our knowledge. These turn out to be particularly interesting for the third generation parameter b . We find two rather narrow minima for tan β = 38 and 50 each, all of them with an acceptable χ 2 < 5. These precise predictions for b have consequences for the SUSY spectrum and phenomenology which nevertheless goes beyond the scope of the current work. Unfortunately, we could not identify here a general set of observables which allows us to distinguish between the two solutions for both tan β values. In principle, once the SUSY spectrum is known we can directly calculate the -parameters which would provide a definite answer.
Until then a more precise determination of the light quark masses and neutrino mixing observables will help to challenge some of the minima which we have found here. The good thing is that we expect a significant progress for these observables in the next years before a detailed knowledge of the SUSY spectrum could realistically be available.
In particular, we would like to highlight here the role of the leptonic Dirac CP violating phase δ l CP which creates a considerable contribution to the χ 2 for most of the found minima. Hence, with the measurement of low energy leptonic CP violation we can test SUSY SO (10).
In summary, we see that with the emergence of precision neutrino data and improved lattice calculations SUSY SO(10) can be challenged. So that in the near future we might find out, if we have been foolishly dreaming or following a path towards a deeper understanding of nature.

A SM Observables at Various Scales
The quark masses (except for the top quark) and the value of α (5) s are taken from Ref. [56], as described in Section 2.3. For the mass of the top quark we use the PDG average m t,pole = 174.2 ± 1.4 GeV [57]. At two-loop accuracy also the quartic Higgs coupling contributes to the RGE running of the Yukawas in the SM This convention is in accordance with the convention used in the two-loop RGEs, i.e. the quartic term in the SM Higgs potential is given by λ(φ † φ) 2 /4. We calculate λ from the Higgs mass measurement m Higgs = 125.09 ± 0.24 GeV [57]. The values for the charged lepton Yukawa couplings and the electroweak gauge couplings have already been given in Ref. [59]  For the CKM matrix we use the updated results of the CKMfitter group as presented at ICHEP 2016 [60]. We translate from the Wolfenstein parametrisation to the parametrisation of the PDG and obtain In Table 3 all these values are given at 1 TeV, 3 TeV and 10 TeV, respectively. The procedure of RGE running is described in Sec. 2.3. For use within SUSY models the values have also been converted to the DR scheme, cf. Table 4.

B GUT Scale Parameters for the Global Minima
For reference, we give here the coordinates of our global minima to a high numerical precision, which allows the reproduction of the results given in Table 5. We determine the GUT scale by minimising the sum of the squared differences of the gauge couplings. The central values of the Yukawa couplings are included in the running to the GUT scale. At M GUT = 1.35276 × 10 16 we find g 1 = 0.706132, g 2 = 0.708008, g 3 = 0.705465.
In all fits, these values remain fixed. The parametrisation of the MSSM Yukawa matrices is given in Eqs. (2.2)-(2.7). If SUSY threshold corrections are included, their value at M SUSY is also given below. The parametrisation of the SO(10) Yukawa couplings is

B.1 Global Minima without SUSY Threshold Corrections
We list here our fit results for the coordinates of the global minima without SUSY threshold corrections for different values of tan β. For reference, we use a high numerical precision.

B.2 Global Minima with SUSY Threshold Corrections
Here we give our fit results for the coordinates of the global minima including SUSY threshold corrections, again with a high numerical precision, see above.

B.3 Local Minima with SUSY Threshold Corrections
Here we give our fit results for the coordinates of the local minima of b including SUSY threshold corrections, again with a high numerical precision, see above.