Data-driven analysis of a SUSY GUT of flavour

We present a data-driven analysis of a concrete Supersymmetric (SUSY) Grand Unified Theory (GUT) of flavour, based on SU(5) × S4, which predicts charged fermion and neutrino mass and mixing, and where the mass matrices of both the Standard Model and the Supersymmetric particles are controlled by a common symmetry at the GUT scale. This framework also predicts non-vanishing non-minimal flavour violating effects, motivating a sophisticated data-driven parameter analysis to uncover the signatures and viability of the model. This computer-intensive Markov-Chain-Monte-Carlo (MCMC) based analysis includes a large range of flavour as well as dark matter and SUSY observables, predicts distributions for a range of physical quantities which may be used to test the model. The predictions include maximally mixed sfermions, μ → eγ close to its experimental limit and successful bino-like dark matter with nearby winos (making direct detection unlikely), implying good prospects for discovering winos and gluinos at forthcoming collider runs. The results also demonstrate that the Georgi-Jarlskog mechanism does not provide a good description of the splitting of down type quark masses and charged leptons, while neutrinoless double beta decay is predicted at observable rates.


Introduction
The Minimal Supersymmetric Standard Model (MSSM) remains an appealing extension of the Standard Model of particle physics, as it provides solutions to the most prominent shortcomings of the latter.In addition to solving the hierarchy problem related to the mass of the Higgs boson [1,2], the model includes a viable candidate for the observed Cold Dark Matter (CDM) in the Universe, namely the lightest of the four neutralinos.Furthermore, the masses of the Standard Model neutrinos can be generated through the Seesaw mechanism [3][4][5][6][7] by including heavy right-handed neutrinos which can easily be implemented in the MSSM.
While collider searches for new physics have remained unsuccessful so far, additional information can be obtained from examining precision observables involving flavour transitions.For example, in the hadronic sector, the branching ratios of rare decays such as b → sγ are sensitive to new physics contributions, especially if they involve non-minimal flavour violation (NMFV), i.e. sources of flavour violation beyond the Cabbibo Kobayashi Maskawa (CKM) matrix [8,9].The same holds for the lepton sector [10,11], where NMFV contributions can induce the branching ratios like µ → eγ and µ → 3e, or µ − e conversion rates in nuclei.Flavour precision observables therefore provide an interesting handle towards the new physics spectrum, and in particular towards the underlying flavour structure.
A particularly interesting feature of the MSSM and related models is that they can be successfully embedded into Grand Unified Theories (GUTs).Such a framework allows a unification of the gauge couplings at a scale of about 10 16 GeV with better precision than the Standard Model alone [33].In the same spirit, the soft breaking parameters related to squarks and sleptons stem from a common origin.In the most simple realizations this allows for the reduction of the number of parameters of the model.Starting from the imposed values at the GUT scale, the phenomenological aspects are obtained through renormalization group running to the TeV scale, where the physical masses and related observables are computed.
The two aforementioned aspects can be addressed by considering SUSY-GUTs including flavour symmetries, such as SU (5) × A 4 [34] or SU (5) × S 4 [35,36], to cite only two examples.In such a situation, the flavour structure of the theory is defined at the GUT scale by the imposed symmetry.Renormalization group running then translates the GUT-scale structure into the observable mass spectrum at the TeV scale.The TeV-scale phenomenology thus inherits a footprint of the imposed flavour structure at the GUT scale.
In a previous study [37], some of the authors have explored the phenomenology of NMFV within a SU (5) × A 4 implementation of the MSSM suggested first in Ref. [34].Based on the variation of the NMFV parameters around a MFV reference scenario taken from Ref. [34], it has been shown that these parameters need to be varied simultaneously in order to cover all phenomenological aspects, in particular since cancellations between different contributions may occur.Moreover, it has been demonstrated that a model may feature a reasonable amount of flavour violation while satisfying the stringent constraints of rare decays such as b → sγ or µ → eγ while sufficient lepton flavour violation may also address baryon asymmetry in the universe [38].It is therefore interesting to pursue the study of GUT implementations of flavour violation, e.g., via flavour symmetries, in the context of low-energy and precision constraints as well as TeV-scale phenomenology.
In the present paper, we shall focus on the case of SU (5) unification combined with an S 4 flavour symmetry, as suggested by one of the authors in Refs.[35,36].In a similar way as in Ref. [37], we will explore in detail the TeV-scale aspects of this model, including observables related to flavour violation and dark matter phenomenology.More precisely, in this study, we aim at a complete exploration of the associated parameters, i.e. including a variation of all relevant parameters at the GUT scale.For the sake of an efficient exploration, we make use of the Markov Chain Monte Carlo technique [39][40][41].
The paper is organised as follows.In Section 2, we review the assumed model.Section 3 is dedicated to the discussion of the Markov Chain Monte Carlo algorithm that we employ to efficiently explore the model parameter space.Results are then presented in Section 4. Our conclusions are given in Section 5.

Fields and symmetries
The model developed in Refs.[35,36,42,43] is based on the grand unifying group SU (5) combined with an S 4 family symmetry, and supplemented by a U (1) symmetry.The lefthanded quarks and leptons are unified into the representations 5, 10 and 1 of SU (5) according to, where the superscript c stands for CP -conjugated fields (which would be right-handed without the c operation), and α = 1, 2, 3 is a family index.The three families are controlled by a family symmetry S 4 , with F and N each forming a triplet and the first two families of T forming a doublet, while the third family T 3 (containing the top quark) is a singlet, as summarised in Table 1.The choice of the third family T 3 being a singlet, permits a renormalisable top quark Yukawa coupling to the singlet Higgs discussed below.The S 4 singlet Higgs fields H 5 , H 5 and H 45 , each contain a doublet SU (2) L × U (1) Y representation that eventually form the standard up (H u ) and down (H d ) Higgses of the MSSM (where the H d emerges as a linear combination of doublets from the H 5 and H 45 ) [44]. 1 The VEVs of the two neutral Higgs fields are where GeV.Just below the SU (5) breaking scale to the usual SM gauge group, the flavour symmetry is broken by the VEVs of some new fields: the flavons, Φ f ρ , which are labelled by the corresponding S 4 representation ρ as well as the fermion sector f to which they couple at leading order (LO).Two flavons, Φ u 2 and Φu 2 , generate the LO up-type quark mass matrix.
energy Higgs doublet H d must necessarily break U (1).Although the discussion of any details of the SU (5) GUT symmetry breaking (which, e.g., could even have an extra dimensional origin) are beyond the scope of our paper, we remark that a mixing of H5 and H 45 could be induced by introducing the pair H ± 24 with U (1) charges ±1 in addition to the standard SU (5) breaking Higgs H 0 24 . Field Table 1.Field content of the model and associated charges and representations.
Three flavon multiplets, Φ d 3 , Φd 3 and Φ d 2 , are responsible for the down-type quark and charged lepton mass matrices.Finally, the right-handed neutrino mass matrix is generated from the flavon multiplets Φ ν 3 , Φ ν 2 and Φ ν 1 as well as the flavon η which is responsible for breaking the tri-bimaximal pattern of the neutrino mass matrix to a trimaximal one at subleading order.An additional U (1) symmetry must be introduced in order to control the coupling of the flavon fields to the matter fields in a way which avoids significant perturbations of the flavour structure by higher-dimensional operators.

Flavon alignments
The vacuum alignment of the flavon fields is achieved by coupling them to a set of so-called driving fields and requiring the F -terms of the latter to vanish.These driving fields, whose transformation properties under the family symmetry are discussed in Refs.[35,36,42,43], are SM gauge singlets and carry a charge of +2 under a continuous R-symmetry.The flavons and the GUT Higgs fields are uncharged under this U (1) R , whereas the supermultiplets containing the SM fermions (or right-handed neutrinos) have charge +1.As the superpotential must have a U (1) R charge of +2, the driving fields can only appear linearly and cannot have any direct interactions with the SM fermions or the right-handed neutrinos.
Using the driving fields, the flavour superpotential may be constructed, resulting in the following vacuum alignments (for details see Refs.[42,43]), where λ = 0.22 is approximately equal to the Wolfenstein parameter [45] and the φ's are dimensionless order one parameters.Imposing CP -symmetry of the underlying theory [46], all coupling constants can be taken real, so that CP is broken spontaneously by generally complex values for the φs.The remaining phases are independent and there is no residual CP -symmetry.M denotes a generic messenger scale which is common to all the nonrenormalisable effective operators and assumed to be around the scale of grand unification.

Yukawa matrices
Because of the non-trivial structure of the Kähler potential, non-canonical kinetic terms are generated.For a proper analysis of the flavour structure, one needs to perform a canonical normalisation (CN) operation, swapping the misalignment of the kinetic terms to the superpotential.Therefore, in the model proposed in Refs.[35,36,42,43], contributions to the flavour texture from both the superpotential and the Kähler potential are taken into account.In this subsection, we shall begin by ignoring such corrections, and also only consider the leading order Yukawa operators, in order to clearly illustrate the origin of the flavour structure in the model.However, all such corrections are taken into account in the phenomenological treatment of the Yukawa matrices in the following subsection.We remark that the model is highly predictive, as the parameters entering the flavour structure are expected to be of O(1) but the overall flavour texture is provided as a function of the expansion parameter λ = 0.22.

Up-type quarks
The Yukawa matrix of the up-type quarks can be constructed by considering all the possible combinations of a product of flavons with T T H 5 for the upper-left 2 × 2 block, with T T 3 H 5 for the (i3) elements, and with T 3 T 3 H 5 for the (33) element.The most important operators which generate a contribution to the Yukawa matrix of order up to and including λ 8 are where the parameters y t and y u i are real order one coefficients.Inserting the flavon VEVs and expanding the S 4 contractions of Eq. (2.6), with T T and Φ u 2

Down-type quarks and charged leptons
The Yukawa matrices of the down-type quarks and the charged leptons can be deduced from the leading superpotential operators where the y d i are real order one coefficients.For the operators proportional to y d 2 and y d 3 , specific S 4 contractions indicated by (• • • ) 1 and (• • • ) 3 have been chosen (justified by messenger arguments) such that the Gatto-Sartori-Tonin (GST) [47] and Georgi-Jarlskog (GJ) [48] relations are satisfied.Separating the contributions of H 5 and H 45 , the S 4 contractions give rise to The parameters in these expressions are related to the flavon VEVs as defined in Eqs. (2.3)-(2.5)via (2.11) The Yukawa matrices of the down-type quarks and the charged leptons are linear combinations of the two structures in Eq. (2.10).Following the construction proposed by Georgi and Jarlskog, we have CKM mixing is dominated by the diagonalisation of the down-type quark Yukawa matrix.Note that in some models it is possible to go beyond the simple case m b = m τ at the GUT scale, by including larger Higgs representations [49,50].

Neutrinos
The neutrino masses originate from a standard Supersymmetric Type I Seesaw mechanism, where the heavy right-handed fields, N, are turning the tiny observed neutrino effective Yukawa couplings into natural parameters.The Lagrangian for the neutrino sector is therefore given by where Y ν is the Dirac Yukawa coupling and M R is the right-handed Majorana mass matrix.The Dirac coupling of the right-handed neutrinos N to the left-handed SM neutrinos is dominated by the superpotential term where y D is a real order one parameter.The corresponding Yukawa matrix is determined as The mass matrix of the right-handed neutrinos is obtained from the superpotential terms where w i denote real order one coefficients.This results in a right-handed Majorana neutrino mass matrix M R of the form with (2.17) The first matrix of Eq. (2.16) arises from terms involving only Φ ν 1,2,3 .As their VEVs respect the tri-bimaximal (TB) Klein symmetry Z S 2 ×Z U 2 ⊂ S 4 , this part is of TB form.The second matrix of Eq. (2.16), proportional to D, is due to the operator w 4 1 M N N Φ d 2 η which breaks the Z U 2 at a relative order of λ, while preserving the Z S 2 .The resulting trimaximal TM 2 structure can accommodate the sizable value of the reactor neutrino mixing angle θ l 13 .It is instructive to show the effective light neutrino mass matrix which arises via the type I seesaw mechanism, and has the form with a ν , b ν , c ν and d ν being functions of the real parameters A, B, C and D. The deviation from tri-bimaximal neutrino mixing is controlled by d ν ∝ D. Due to the three independent input parameters (w 1 ∝ A, w 2 ∝ B, w 3 ∝ C), any neutrino mass spectrum can be accommodated in this model.However, in our numerical analysis, we shall restrict our scans to the case of a normal neutrino mass ordering, which is preferred by the latest global fits of neutrino oscillation data.Note that this expression, while generally providing a good estimation of the effective neutrino mass matrix, is only an approximation valid at order D × λ and therefore, it does not strictly hold when considering a potential O(10) D-parameter.This is why this expression is more for illustrative purposes while we performed a rigorous treatment of the full seesaw mechanism in the numerical analysis.

Phenomenological Yukawa couplings at the GUT scale
The true model predictions at the high scale differ from those shown previously, since they also involve other higher order corrections to the Yukawa terms, and one must also include the effects of canonical normalisation (CN) leading to the matrices in Ref. [36].For simplicity, while keeping the phenomenology indistinguishable from the constructed model, we allow for minor approximations, and here we summarise the form of the Yukawa matrices that we actually assume at the GUT scale.
Concerning the up-type quark Yukawa matrix, we shall continue to take it to be diagonal as the off-diagonal entries are much more λ suppressed than the diagonal ones.We may also absorb the phases into a redefinition of the fields.
Since the CKM matrix is controlled by the down-type quark Yukawa matrix, we shall include some of the higher-order terms and some of the effects of CN, in order to obtain a perfect fit to quark data.Therefore there are some corrections to the GST relations.
The charged lepton Yukawa matrix is, like in any standard SU (5) model, closely related to the down-quark Yukawa matrix as per Y Y T d , together with a modified via the GJ mechanism through the incorporation of 45 and 5 Higgs representations in order to generate a more reasonable relation between m s and m µ .
The explicit Yukawa matrices we will use for the charged fermionic sector are therefore provided by the following expressions: where δ = 2θ d 2 + 3θ d 3 , and θ d ρ (ρ = 2, 3) corresponds to the phase of a ρ-representation flavon in the original model.Note that our analysis, including the soft masses discussed below, relies only on the two phases θ d 2 and θ d 3 .These Yukawa matrices are obtained using Eq.(2.10).However, one needs to perform a canonical normalisation of the kinetic terms.This procedure has been taken care of in Ref. [36].Therefore, our phenomenological Yukawa couplings involve more parameters than the ones exposed in Eq. (2.10).The up-type quark Yukawa has been approximated with respect to Ref. [36].
In the neutrino sector, the effects of CN are negligible, and we therefore take these matrices to have the same form as given previously, The Dirac neutrino coupling, neglecting the O(λ 4 ) terms, compared to the original paper [36], that is also of the form given in the previous subsection, (2.21)

SUSY breaking terms
We now consider the SUSY breaking sector of the low energy scale MSSM generated after integrating out the heavy degrees of freedom.In the context of the standard phenomeno-logical R-parity conserving MSSM, the soft Lagrangian is parametrised as where φ denotes the generic superpartner of a generic SM particle φ.Assuming that the SUSY breaking is controlled by some hidden sector mediated by a superfield X, the soft parameters described in (2.22) are generated when X develops a VEV in its F-term at the SUSY breaking scale.Furthermore, we consider that the SUSY breaking mechanism is independent of the flavour breaking one.Note that we assume non universal gaugino masses in our analysis.While this is the case in simplistic SU (5) models, many extensions can account for non universal masses, see e.g.Ref. [51].
The new flavour structure arising from the SUSY breaking sector is also controlled by the flavour symmetry, in a similar fashion as the SM texture is.Extracting the results from Ref. [36], we first summarise the predictions for the soft trilinear terms A f , (2.23) The trilinear soft couplings exhibit the same structure as the Yukawa terms, except that the O(1) parameters are now different.Note that in our numerical analysis we neglect the phases appearing in the trilinear matrices.We have verified that such phases have a negligible affect on the CP -conserving constraints.We further use the same approximations as the ones considered for the Yukawa couplings.
Similarly, we summarise the results on the soft scalar mass matrices, where N is the right handed sneutrino term, Because of the unification, all the sfermion soft matrices are linked to the soft matrix of the SU (5) representation they belong to.
3 Data-driven model exploration

Algorithm
The full analysis of the parameter space relies on a Markov-Chain Monte Carlo technique [41], and more specifically the Metropolis-Hastings algorithm [39,40].This technique allows one to perform a sophisticated data-driven exploration of an high-dimension parameter space.The idea behind the algorithm is to estimate the likelihood L of a given set of parameter values θ with respect to the set of observables O.For simplicity and the rest of the analysis we assume that the observables are not correlated, i.e.
where σ i is the uncertainty associated to the observable O i .Successively, random values of the parameters, picked around the previous ones, are evaluated at each iteration.In our implementation, the new proposed parameter value θ n+1 is obtained through a Gaussian jump, where G (a, b) is a Gaussian distribution centered around a with width b, κ is a parameter that needs to be tuned empirically for the algorithm and θ max i and θ min i stands for the extrema values of the θ i considered range. If ), the point is accepted and the chain continues from this point.Otherwise, the new point is accepted with probability This avoids falling into local minima, and thus allows for a better parameter space exploration.In practice, we randomly choose a number µ ∈ [0, 1] such that the test succeeds Fail: restart at n Success: if µ < p. Otherwise, the point is rejected, and we reevaluate the step n + 1 for another proposal set of parameters deduced from step n.Within this framework, the algorithm can move across larger regions while still converging to highest likelihood regions.
In high-dimensional parameter space, the quality of the exploration relies more on the total of chain numbers than the length of the chain themselves.Indeed, different starting points (chosen randomly) can lead to different likelihood maximums.A summary of the algorithm is given in Figure 1.

Constraints, tools and setup
We now develop on the numerical tools and constraints employed in our analysis.As the model is defined at the GUT scale, we first perform the evolution of the renormalisation group equations (RGEs) to the low scale, to derive low energy observables.For this purpose we employed the SARAH v4.14.1 Mathematica package [52][53][54][55][56] in order to generate a type I Seesaw GUT MSSM model based on SPheno v4.0.4 [57,58].Right handed Majorana neutrinos, which typically live near the GUT scale are therefore consistently integrated out at their mass scale.Furthermore, the Flavour kit [59] available within SARAH / SPheno computes a wide range of flavour observables, simplifying our framework as both one-loop masses, two-loop Higgs mass [60][61][62][63] and flavour observables are evaluated within a single executable.
However, modifications of this model have been realised: In the usual SPheno instances, SM fermion masses are enforced to match the experimental data by several runs up and down between the GUT scale and the low scale, rendering our model predictions impossible to estimate.To overcome this, we have removed the SM fermions from this iterative convergence process while keeping the massive gauge bosons.To consistently implement such restrictions, attention must be payed concerning several features.An extended discussion regarding our modified SPheno version can be found in Appendix A.
We have also decided to include dark matter constraints in our analysis.Restricting ourselves to neutralino dark matter, we imposed a step dark matter candidate likelihood (1 if the LSP is the lightest neutralino, 0 otherwise).In order to derive relic density and direct detection constraints, we have used micrOMEGAs v5.2, which accepts the spectrum files generated by SPheno through the SUSY Les Houches Accord (SLHA) [64,65] interface.
It should be noted that SARAH generated micrOMEGAs models are in general limited to real Lagrangian parameters, that is all couplings need to be real.This caused problems in our relic density calculations due to the presence of phases in multiple sectors.To overcome this, we maintained a full calculation including phases within SPheno but recast the model with real valued couplings (taking the modulus) for the relic density calculation.As, in general, the CP -violating contribution to relevant (co-)annihilation channels are limited, and our CP -violating parameters are also numerically rather small, this is a valid approximation.We have verified for a few cases, that the effect of the phases had little impact on the amplitudes squared of the relevant processes.
Linking these tools together, we are able to investigate a wide range of constraints.The list of the Standard Model parameters to be fitted is given in Table 3 while the flavour and dark matter constraints are listed in Table 4.We give the list of input parameters and their respective scanning range in Table 2.
In all tables, the upper bounds constraints are given at the 90% confidence level and in order to help the process of chain convergence and initialization, we postulate a smoothing step function for the upper limits constraints where we chose, somehow arbitrarily, a common value of On the other hand, we associate a Gaussian likelihood function for all experimentally measured observables σ i being the uncertainty given in Tables 3 and 4.
Regarding the dark matter direct detection constraints, we extracted and extrapolated the curves from Refs.[73,74] as shown in Figure 2, while for the relic density we have used the results from Ref. [69] adding a 10% uncertainty due to micrOMEGAs precision in combination with underlying cosmological assumptions.For the other constraints we are using Table 2. GUT scale input parameters for the model and their scanning range.For all parameters, the step size for a Markov chain iteration is given as 0.5% of the total range length of the allowed values.This step size was procured by trial and error in order to balance scan efficiency and a search of the parameter space.In addition, we set a fixed value for the following parameters: sign(µ) = −1; and λ = 0.22 and M GUT = 2 the current experimental uncertainties associated to the values given in the different tables while adding in quadrature a theoretical constraints on the different standard model masses.
The theory uncertainty on the Higgs mass is fixed at 2 GeV [75] while we assume a common 1% uncertainty on the different fermion masses because of RGE fixed order precision and changes from the DR to the on-shell renormalisation scheme.If no experimental constraints is present for a given value in the tables it is understood that theoretical uncertainties are by far dominant with respect to the experimental ones.Finally, the different quark masses are extracted at different scales, i.e.Q = 2 GeV for q = (u, d, c, s) and Q = m b for q = b.
Having implemented and executed the above, 197 chains were recovered.As the parameter space was so vast and such a large number of very precise constraints were used, the efficiency of these scans was very low, requiring weeks of computer time to complete.Therefore, the scans were allowed 2000 Markov chain steps each.After this process we collect the data and applied a likelihood cutoff such that only points with relatively high likelihood are left in the final data set.This was in order to prevent the distributions presented here-in from misleading the reader into thinking some parts of parameter space BR(τ − → µ + e − e − ) < 1.5 Ωh 2 0.12 ± 0.012 th.[69][70][71][72] Direct detection cf. Figure 2 [73,74] Table 4. Leptonic flavour and dark matter constraints.These upper limits numbers are given at the 90% confidence level.For the dark matter relic density we assume 10% theoretical uncertainties because of cosmological assumptions.were viable when, in fact, they produce excessively low likelihoods.The likelihood cutoff applied was 10 −150 .Although this is tiny, much of the poor likelihood comes from poor convergence of the fermion masses (see discussion in Section 4.1).In general, the remaining constraints converged very well to the observed values.As an example of two constraints, see Figure 3, which shows how the Higgs boson mass and relic density are centered around the expected value.

Results
The subsequent results are based on the Markov Chain Monte Carlo (MCMC) study following the methods elucidated previously.Having already presented two illustrative plots showing the constraints used to guide the MCMC, we now present the resultant spectra and phenomena.We begin with a discussion of the fermion masses, mixing, and a general discussion of the model's success in recreating the Standard Model observables.We then look at the supersymmetric (SUSY) spectrum, the dark matter sector, and further phenomenological results.Finally, we give a discussion of the effects on collider physics and experimental physics more generally.

Fermion masses and mixing
In this subsection we present the results of our scan for fermion masses and mixing parameters, which are put in as constraints as shown in Table 3.The results for the fermion masses are shown in Figure 4, while those for the mixing parameters are shown in Figures 6 and 7.These results follow from the charged fermion Yukawa matrices at the GUT scale shown in Eqs.(2.19), together with the neutrino Dirac Yukawa matrix in Eq. (2.21) and the heavy right-handed neutrino mass matrix in Eq. (2.20).Note that the (3, 3) entries of the charged lepton and down type quark Yukawa matrices are equal at the GUT scale (yielding approximate bottom-tau unification m b = m s ), while the (2, 2) entries of these matrices differ by the Georgi-Jarlskog (GJ) factor of 3 (yielding an approximate strange to muon mass ratio m s = m µ /3 at the GUT scale).
The results for fermion masses in Figure 4 show that the above GJ relations do not lead to phenomenologically viable charged lepton and down type quark masses at low energy, in particular m s , m b and m τ are not well fitted.This problem has also been noted by other authors, and possible solutions have been proposed based on various alternative choices of GUT scale Higgs leading to different phenomenologically successful mass ratios at the GUT scale [49,50].Since the purpose of this paper is to perform a comprehensive phenomenological analysis on an existing benchmark model, we shall not consider such alternative solutions here, but simply note that such solutions exist and could be readily applied.
Note that the absolute values of the neutrino masses in Figure 4 are genuine predictions of the model, since only the experimentally measured mass squared differences in Table 3 were put in as constraints.In particular, the lightest neutrino mass distribution is peaked around a few times 10 −3 eV.This leads to an interesting prediction for neutrinoless doublebeta becay.In Figure 5 we give the model prediction of the neutrinoless double-beta decay parameter m ee against the mass of the lightest neutrino m ν1 .m ee is given by where m νi are the light neutrino masses and U is the PMNS matrix including the Majorana phases.We use the same convention as in Ref. [66].Future projections for CUORE rule The first generations are very well fitted.However, due to the link between the down type quarks and the leptons at the GUT scale, the second and third generation masses are slightly off.Furthermore, the top mass is also slightly poorly aligned due to the Higgs mass constraint.Note that few points exhibit neutrino masses above the visible part of the histograms.
out approximately 52% of the data set.Indeed, the model tends to favour relatively high values of m ee as compared to the theoretically allowed region.
The results for the CKM and PMNS mixing parameters in Figures 6 and 7 show a good fit to the constraints.The PMNS mixing parameters sin θ PMNS  [76], while the coloured horizontal lines show the limits on |m ee | from KamLAND [77], EXO-200 [78], CUORE [79], and GERDA [80].We also indicate future prospects for CUORE [81].
the phases from δ CKM and δ P M N S .Therefore, the high scale phases, who determine the Majorana phases, must be correlated.However, the striking nature of the correlation is surprising.It seems that, roughly speaking, the Majorana phases must sum to a multiple of π.Of course, these phases are important with regards to CP -violating processes but this should suggest that their individual values should be multiples of π, not their sum.For the time being we leave this puzzle as a comment to be understood in greater detail.

SUSY spectrum
Figure 10 shows the distribution of the lightest charged sfermion masses.We recall from table 1 that in the sfermion sector the SU (5) 5-plets, the first two generation 10-plets and the third generation 10-plet transform as a triplet, doublet and singlet under S 4 , respectively.The nature of the lightest sfermions depens on the mass hierarchy of these states at the GUT scale.In case that S 4 doublet is the lightest one, then the lightest sfermion consists of a strong admixture of the first two generation sfermions.More precisely, the lightest slepton, d-type squark and u-type squark will be an admixture of ẽR and μR , dL and sL , and ũR and cR , respectively.In all other cases these states will be sfermions of the third generation.On average the corresponding mass splitting is larger in these cases compared to the first one.This can be traced back to the large top Yukawa coupling as well as to a possible difference between the masses of the 10-plet and the 5-plet at the GUT scale.The neutralino, chargino, and gluino masses are displayed in Figure 11.Although the gluino mass is not particularly constrained, the two lightest neutralinos and the lightest chargino have a very constrained spectrum.As, a priori, the relic density is too high, the model requires a specific mechanism to reduce the dark matter relic density to phenomenological values.As much of the rest of the spectrum is large, the neutralinos and charginos supply an alternative mechanism via co-annihilation.In order to allow for such contributions, the lightest gauginos must be comparable in mass as will be seen in the dark matter dedicated section below.
The last panel of Figure 11 demonstrates the "decoupling limit" for two Higgs doublet models.In this limit, the three additional Higgs states have very large masses and are approximately degenerate.
Finally, Figure 9 depicts the mass distribution of the heavy neutrinos.Overall, the masses are very large and approximately of the same order of magnitude (10 14 GeV).

Dark matter
We now come to the discussion of dark matter aspects of the model under consideration.As we have seen in Figure 3, the dark matter relic density given by the latest Planck results is well accommodated for in the parameter regions surviving the numerous imposed   constraints.The corresponding parameter configurations feature essentially bino-like dark matter, which can be understood from Figure 12, where we depict the relevant bino and wino content of the lightest neutralino.
Looking at the gaugino masses (Figure 11), it can be seen that the second-lightest neutralino as well as the lighter chargino lie very close to the lightest neutralino.In other words, the bino and wino mass parameters are, at the SUSY scale, almost equal, the bino lying just below the wino mass.This feature is driven by co-annihilations needed to achieve the required relic density.Note that this corresponds to a situation where the GUT-scale values of the bino and wino mass differ roughly by a factor of two.
It is interesting to note that, although the bino-like lightest neutralino χ0 1 is the dark matter candidate, the (co-)annihilation cross-section is dominated by the (co-)annihilaton of the wino-like states χ0 2 and χ± .This is explained, on the one hand, by the very small mass difference between the wino-like states, and, on the other hand, by the enhanced annihilation cross-section for the latter as compared to the bino-like state.Typical final states of these (co-)annihilation channels are quark-antiquark pairs and gauge boson pairs (including Z 0 , W ± , and γ).Neutrino final states are subdominant.The presence of co-annihilation can also be understood through Figure 13 showing the correlation of the three lightest gaugino masses and the dark matter relic density.Let us finally note that scenarios with wino-like dark matter would give rise to insufficient relic density to align with the experimental evidence, as the wino (co-)annihilation cross-section is numerically more important as the one for the bino.
Coming to the direct dark matter detection, we can see from Figure 14 that this constraint is also well satisfied in the model under consideration, both for the spin-dependent and the spin-independent case.It is important to note that all points shown in Figure 14 lie also below the projected limits of the XENONnT experiment [82].The fact that all points are found below this limit can be traced to the fact that we have applied a cut on the global likelihood value as explained in Section 3.This procedure discards the points which are too close to the current XENON1T limit, since they typically feature a somewhat lower likelihood value.This means that parameter configurations with reasonably high global likelihood values may not be challenged by direct detection experiments in a near future.
In summary, the relic density constraint implies a relatively small mass difference between the lightest neutralino and the next-to-lightest states, leading to final states with soft pions and leptons which are difficult to detect.The current bounds depend on the nature of the NLSP go up to masses of about 240 GeV [83][84][85][86].This cuts slightly into the allowed parameter space.

Collider related aspects
As already mentioned in Section 4.2 the flavour structure of the lightest sfermions falls into two extreme case: For each sfermion type, we observe either a first and second generation fully mixed state or a strict third generation state.This feature is illustrated in Figure 15 where we show the distribution of the sum of the square of the mixing matrix entries.While at first glance it seems that this particular prediction of the model would be quite interesting from a collider perspective, enabling potential flavour mixed search channels, the model also predicts rather high masses for sfermion states.The lightest squark masses are peaked around 4 TeV and stand well beyond any potential collider sensitivity reach (see for e.g.Ref. [87] where the limits are around 1.8 TeV squarks using simplifying assumptions that do not hold in our case).The lightest slepton on the other hand can be as light as 1 TeV.However, the production cross section drops significantly with respect to the QCD dominated squark ones.Furthermore, our model naturally predicts right-handed slepton states to be the lightest ones which further decreases the production cross section [88].As an example, recent searches for flavour conserving channels with nearly helicity degenerated slepton states are excluding masses of the order 600 -700 GeV [89] if m χ0 1 < 150 GeV (and limits are even weaker if this is not the case).In addition, current exclusion on τ flavoured sleptons masses are of order 390 GeV [90].
It turns out that our model predictions regarding the sfermionic sector are not very promising for potential collider searches.However this flavour mixing particularity leads to some interesting features from indirect searches perspectives.It is well known that slepton mixing can generate significant contributions to flavour violating decay constraints; in particular, if the mass of the slepton is rather light.In our case, BR(µ → eγ) illustrates very well this feature as being one of the most stringent test for lepton flavour violation.Figure 16 shows the distribution of this branching ratio as a function of the mass of the lightest slepton.While the points are within the current experimental limits, it appears that future prospects in the current MEG II experiment [91] could rule out a vast majority of our light slepton points, implying the interesting conclusion that this particular constraint and other flavour violating constraints would have more discriminatory power than classic direct slepton searches and propose typical smoking guns for our framework.We also note that, while we can have light staus in the model, τ → eγ limits are far from constraining BR( e ) our model with branching ratios below 10 −13 while current limits are around 10 −8 .
Despite the sfermionic states being unreachable in near future collider searches, this is not necessarily the case for the other SUSY particles in our model.In particular, the model can predict rather light electroweakinos with sufficient mass gap for collider considerations.As a comparison, recent searches from ATLAS [92] put a lower bound on m χ0 1 of 270 GeV when the mass gap with m χ± 1 is of order 50 GeV.While this very light masses are not present in our current framework, we can hope for more stringent limits from future LHC runs.
Similar conclusion can be derived regarding gluino searches.While the predicted masses lie in the 2 -6 TeV range, ATLAS and CMS limits using simplified models [87,93] reaches 1.8 TeV exclusion.Again, we can argue that future LHC runs might restrict our light spectrum parameter distributions.
To illustrate the collider phenomenology we present three benchmark points where relevant information is given in Table 5.The detailed information is added as supplementary material to the arXiv submission of this paper.In particular, we list the masses of the electroweakinos for an illustrative point with low masses for the light gauginos and a relatively small mass gap.We also give the dominant decay channels and decay widths.All other particles are too heavy to be detected in the upcoming LHC run.We give a benchmark for a point which features, besides the light chargino and the neutralinos of BP1, a gluino with mass 2 TeV which should be in reach of the upcoming LHC run.Again, we give the masses, decay widths, and decay channels for the particles under examination.Regarding the third benchmark point, while having similar features as BP2 in view of collider physics, it further highlights the fully mixed nature of the lightest sleptons and the potential relevance for µ → eγ.Its branching of 4 • 10 −13 is close to the current experimental bound.However, even if this rare decay is discovered in an upcoming experiment, this example shows that it will be rather challenging to detect lepton flavour at a high energy collider as the sleptons are quite heavy.χ1 q * q χ± 1 q * q χ2 q * q 13% 50% 37% Table 5. Selected benchmark points (BP) for phenomenology.The first BP exhibits rather light electroweakinos and might represent a challenge for future colliders.The second BP has a rather light gluino, which is already close to the current experimental limits from the LHC collaborations.Finally, the last BP is an example of maximally mixed lightest sleptons.All masses and widths are given in GeV.

Conclusion
In this paper we have presented a detailed phenomenological analysis of a concrete Supersymmetric (SUSY) Grand Unified Theory (GUT) of flavour, based on SU (5) × S 4 .The model predicts charged fermion and neutrino mass and mixing, and where the mass matrices of both the Standard Model and the Supersymmetric particles are controlled by a common symmetry at the GUT scale, with only two input phases.The considered framework predicts small but non-vanishing non-minimal flavour violating effects, motivating a sophisticated data-driven parameter analysis to uncover the signatures and viability of the model.The computer-intensive Markov-Chain-Monte-Carlo (MCMC) based analysis performed here, the first of its kind to include a large range of flavour as well as dark matter and SUSY observables, predicts distributions for a range of physical quantities which may be used to test the model.The predictions include maximally mixed sfermions, µ → eγ close to its experimental limit and successful bino-like dark matter with nearby winos (making direct detection unlikely), implying good prospects for discovering winos and gluinos at forthcoming collider runs.The results also demonstrate that the Georgi-Jarlskog mechanism does not provide a good description of the splitting of down type quark masses and charged leptons.However neutrinoless double beta decay, which depends on a curious pattern of Majorana phases resulting from the two input phases, is predicted at observable rates.
The analysis here may be repeated for any given SUSY GUT of flavour, leading to corresponding predictions for fermion masses and mixing as well as SUSY masses and flavour violating physical observables at colliders and high precision experiments.The results here exemplify the synergy between the theory of quark and lepton (including neutrino) mass and mixing, dark matter and the SUSY particle spectrum and flavour violation, that is possible within such frameworks.It is only by systematically confronting the detailed predictions of concrete examples of SUSY GUTs of flavour with experiment that the underlying unified theory of quark and lepton flavour beyond the Standard Model may eventually be discovered.
• at M GUT : Y , Y ν , Y d , Y u , arg(µ) (in practice the sign of µ), as well as the soft SUSY breking terms in a non-universal form: scalar mass squares m 2 f ( f = . . .), trilinear couplings A f ( f = . . . ) and non-universal gaugino mass parameters M 1 , M 2 , M 3 .All phases can be non-zero in principle.
• at M SUSY : tan β The calculation is done in an iterative way: 1.The gauge couplings are evolved from the electroweak scale using the SUSY RGEs at the one-loop level to M GUT = 2 • 10 16 GeV.We do not require unification at this scale.These couplings at this stage serve only as a starting point for our iteration.
2. All parameters are evolved from M GUT to M SUSY = (M Q ) 33 (M U ) 33 using RGEs at the two-loop level.The right-handed neutrinos are decoupled at their respective mass scale during this evaluation and the contributions to the Weinberg operator are calculated at these scales as well.The running of this operator is taken into account as well.
3. The SUSY spectrum is calculated at the scale M SUSY at the one-loop level and the heavy Higgs masses at two-loop level taking into account the contributions from third generation sfermions and fermions, see Ref. [75] for a summary.In addition the matching to the SM-parameters is performed as described in Ref. [94].
4. The SM-parameters are evolved from M SUSY to m Z using the two-loop SM RGEs.At m Z , the masses of the SM fermions are calculated and the mass of the Higgs boson at the two-loop level.For these calculations we however take G F , m Z , α(m Z ) and α s (m Z ) to calculate the gauge couplings and the VEV.
5. The gauge and Yukawa couplings as well as the quartic Higgs couplings are then evolved to M SUSY using the two-loop RGEs.At this scale the SUSY threshold corrections to the gauge and Yukawa couplings are taken into account.The resulting couplings are evolved to M GUT = 2 • 10 16 GeV.Then steps 2. to 4. are repeated until a relative precision of all masses at the level of 10 −5 is reached.
6. Once the precision goal for the spectrum has been achieved, the flavour observables are calculated.Also here we have modified the procedure slightly: we have a quite heavy spectrum leading potentially large logs of the form ln(M SUSY /m t ).For this reason the calculation is done in two steps: (i) Calculate the SUSY contributions to the Wilson coefficients at Q = M SUSY .(ii) Calculate the SM contributions to the Wilson coefficients at Q = m t .(iii) Add both contributions to calculate the relevant observables.

Figure 1 .
Figure 1.Illustration of the MCMC algorithm utilisation.

Figure 2 .Figure 3 .
Figure2.Dark matter direct detection limits in plane of dark matter mass and spin-(in)dependent nucleon scattering cross section.The dots correspond to data extracted from Refs.[73,74], while the solid line is the extrapolation we performed.

Figure 4 .
Figure 4.The fermion masses are displayed where the red (light red) region indicating the 1σ (2σ) limits.The first generations are very well fitted.However, due to the link between the down type quarks and the leptons at the GUT scale, the second and third generation masses are slightly off.Furthermore, the top mass is also slightly poorly aligned due to the Higgs mass constraint.Note that few points exhibit neutrino masses above the visible part of the histograms.

12 and sin θ PMNS 13 alsoFigure 5 .
Figure 5.The dark blue points are the model prediction of the neutrinoless double beta decay parameter |m ee | vs the mass of the lightest neutrino m ν1 .The light blue shaded region shows the allowed range in this plane for a normal hierarchy as predicted by the model.The vertical grey shaded bands to the right show the current Planck disfavoured region[76], while the coloured horizontal lines show the limits on |m ee | from KamLAND[77], EXO-200[78], CUORE[79], and GERDA[80].We also indicate future prospects for CUORE[81].

Figure 6 .
Figure 6.The CKM parameters are displayed where the red (light red) region indicating the 1σ (2σ) limits.All parameters in the CKM are fitted very well with an approximately Gaussian distribution around the expected value.

Figure 7 .Figure 8 .Figure 9 .
Figure 7.The PMNS parameters are displayed where the red (light red) region indicating the 1σ (2σ) limits.Unlike the CKM parameters, the PMNS shows more variation due to the less stringent experimental constraints.

2 Figure 10 .
Figure 10.Distribution of masses for the lightest up-type, down-type sfermions and the lightest slepton.The slepton is the lightest of these particles.

Figure 11 . 2 Figure 12 . 2 Figure 13 .
Figure 11.Distribution of the masses of the gauginos and higgses.The mass spectrum of the lightest chargino and two lightest neutralinos is compressed to provide co-annihilation mechanism for dark matter.

)Figure 14 .
Figure 14.The direct detection limits for spin dependant and spin independent cross-sections are shown with the experimental limit plotted.The solid line indicates the XENON1T limits[74], while the dashed line in the first panel indicates the expected limit for the XENONnT[82] experiment.We have precluded a representation of the neutron direct detection calculation as all data points are far away from the exclusion limit, much like the proton calculation of the same.

Figure 15 .
Figure 15.The mass of the lightest sparticle for up and down type squarks and sleptons is plotted against the third generation content for the given particle on the left panels.When the third generation content is close to 0, the sfermion features a second and first generation maximally mixed state.For low mass down squark and slepton states, an admixture of first and second generation is favoured.The right panels illustrate the actual proportion of points that belong to these two extreme flavour cases.The sfermion mixing matrices R f ( f = ũ, d, ẽ) are defined according to the SLHA 2 standard[65].

Figure 16 .
Figure 16.BR(µ → eγ) as a function of m ẽ1 .Current experimental limit is represented by the red solid line while the dashed one corresponds to the future prospects of MEG II (6 • 10 −14 ) [91].

Table 3 .
νSM parameters, masses and EWSB constraints for our model exploration.All masses are given in GeV and are pole masses, except for the bottom and light quarks: the bottom (light quarks) one is the MS mass given at the scale Q = m b (µ = 2 GeV).Theoretical uncertainties of 1% are assumed for the different masses and are added in quadrature with the experimental ones.Note that the charged lepton and Higgs boson mass experimental uncertainties are negligible with respect to the theoretical ones and are therefore omitted.