A natural $ S_4 \times SO(10) $ model of flavour

We propose a natural $ S_4 \times SO(10) $ supersymmetric grand unified theory of flavour with an auxiliary $\mathbb{Z}_4^2 \times \mathbb{Z}_4^R$ symmetry, based on small Higgs representations (nothing larger than an adjoint) and hence a type-I seesaw mechanism. The Yukawa structure of all fermions is determined by the hierarchical vacuum expectation values of three $ S_4 $ triplet flavons, with CSD3 vacuum alignments, where up-type quarks and neutrinos couple to one Higgs $\mathbf{10}$, and the down-type quarks and charged leptons couple to a second Higgs $\mathbf{10}$. The Yukawa matrices are obtained from sums of low-rank matrices, where each matrix in the sum naturally accounts for the mass of a particular family, as in sequential dominance in the neutrino sector, which predicts a normal neutrino mass hierarchy. The model accurately fits all available quark and lepton data, with predictions for the leptonic $CP$ phase in 95$\%$ credible intervals given by $ 281^\circ<\delta^\ell<308^\circ $ and $ 225^\circ<\delta^\ell<253^\circ $. The model reduces to the MSSM, with the two Higgs doublets emerging from the two Higgs $\mathbf{10}$s without mixing, and we demonstrate how a $\mu$ term of $\mathcal{O}$(TeV) can be realised, as well as doublet-triplet splitting, with Planck scale operators controlled by symmetry, leading to acceptable proton decay.


Introduction
The flavour problem, the origin of the three families of quarks and leptons with their observed pattern of masses and mixing, remains one of the longest-standing and deepest mysteries left unanswered by the Standard Model (SM) [1]. Among its features are very hierarchical charged fermion masses, with the up-type quark mass hierarchy m u m c m t being stronger than for the down-type quark masses m d m s m b which resemble more closely the charged lepton masses m e m µ m τ . The lightest charged fermion is the electron, with m e ∼ 0.5 MeV, while the heaviest (third family) masses satisfy m τ ∼ m b m t at high energies. Quark mixing is small and hierarchical: V ub ∼ λ 3 , V cb ∼ λ 2 , V us ∼ λ, where λ ≈ 0.225 [1].
The discovery of neutrino mass and mixing [2], elucidated over the past twenty years [3], makes the flavour problem more acute [4], but also provides new features, namely small neutrino masses, and large lepton mixing resembling tri-bimaximal (TB) mixing, but with non-zero reactor angle: [5]. The origin, nature and ordering of the neutrino masses remain open questions, but cosmology suggests that all neutrino masses must be below about 100 meV [6], making them by far the lightest (known) fermions in nature. There are hints of CP violation in the lepton sector [7] but so far it has not been established, unlike in the quark sector where the CP phase is well measured [1].
One of the most attractive and popular possibilities for generating small neutrino masses is the type-I seesaw mechanism involving three right-handed (RH) neutrinos [8]. A natural way to obtain large lepton mixing and normal neutrino hierarchy within type-I seesaw is to assume the sequential dominance (SD) of RH neutrinos [9]. SD postulates three RH neutrinos, where one of them, usually the heaviest one, is almost decoupled from the seesaw mechanism, and is responsible for the lightest physical neutrino mass m 1 . Of the remaining two, one gives the dominant seesaw contribution and is mainly responsible for the (heaviest) atmospheric neutrino mass m 3 and mixing, while the other gives a subdominant contribution, responsible for the (second-heaviest) solar neutrino mass m 2 and mixing. SD therefore predicts m 1 m 2 m 3 ∼ 50 meV.The magnitude of atmospheric and solar mixing is determined by ratios of Yukawa couplings, which can easily be large, while the reactor mixing is typically U e3 O(m 2 /m 3 ) ≈ 0.17. This successful prediction was made over a decade before the reactor angle was measured [10].
To obtain precise predictions for mixing one can go further and impose constraints on the Yukawa couplings, as in constrained sequential dominance (CSD) [11][12][13][14][15]. A particularly successful scheme is known as CSD3 [12,[14][15][16][17] where the neutrino Yukawa matrix is controlled by particular vacuum expectation values (VEVs) of three triplet flavon fields, φ i , as discussed later. The flavon vacuum alignments are fixed by a superpotential which we do not specify here, but may be enforced by an S 4 symmetry, as discussed in [16]. After implementing the seesaw mechanism, the above flavons yield a light effective left-handed (LH) Majorana neutrino mass matrix, where Y ij ∼ φ i φ j T , up to S 4 Clebsch-Gordan (CG) factors. Each of the matrices Y ii is quadratic in φ i and therefore has rank 1. The SD condition implies that µ 2 > µ 1 and hence maximal atmospheric mixing is controlled by Y 22 , solar mixing is controlled by Y 11 , while Y 33 plays no important role in neutrino physics due to the smallness of µ 3 , which implies that m 1 is similarly small. In the limit µ 3 = m 1 = 0, this effectively reduces to a two RH neutrino model, with only two real parameters µ 1 and µ 2 . The resulting model is known as Littlest Seesaw (LS) [15][16][17]. 5 It is remarkable that, with µ 1,2 fitted to the neutrino masses, the entire PMNS matrix is then uniquely determined with no free parameters, giving predictions for mixing angles and the CP phase in agreement with current data.
In order to extend the above ideas into a unified model of flavour, the CSD3 model with two RH neutrinos may be embedded into an A 4 × SU (5) Grand Unified Theory (GUT) [18]. However, from a theoretical point of view, the SO(10) GUT is preferred, since it predicts three RH neutrinos and makes neutrino mass inevitable. This motivated the ∆ 27 × SO(10) model [19], where the mass matrices for all quarks and leptons had the same universal form as the neutrino mass matrix in Eq. 1, but with different coefficients multiplying each rank-1 matrix. While successful and rather complete, there were several unsatisfactory features of the ∆ 27 × SO(10) model. It was quite a complicated model involving a large number of fields. Furthemore, ∆ 27 cannot enforce the desired CSD3 vacuum alignments by symmetry alone, unlike S 4 , which therefore seems to be preferred by CSD3. Finally, a universal mass matrix structure, which seems quite appealing at first sight, led to problems in the quark sector, which were fixed by adding an extra non-universal term in the up-type quark sector, together with some degree of fine-tuning between matrix coefficients in order to obtain the correct quark masses and mixing angles.
In the present paper we propose a natural S 4 × SO(10) grand unified theory of flavour in which the CSD3 model of neutrinos is embedded. Our guiding principles are firstly simplicity, involving the fewest number of low-dimensional fields, secondly naturalness, and thirdly completeness, in particular addressing the doublet-triplet splitting problem. What does natural mean? For us it means that we have a qualitative explanation of charged fermion mass and mixing hierarchies, as for neutrino mass and mixing, with all dimensionless parameters O(1), and in particular that the Yukawa matrices are obtained from sums of low-rank matrices, as in Eq. 1, where each matrix in the sum naturally accounts for the mass of a particular family, analogous to SD in the neutrino sector. This qualitative picture of "universal sequential dominance" is underpinned by a detailed quantitative fit of the fermion spectrum.
In order to achieve this, we shall introduce two Higgs 10s, H u 10 and H d 10 , which will give rise, at low energy, to the minimal supersymmetric standard model (MSSM) Higgs doublets, h u and h d , respectively, with no appreciable Higgs mixing effects. Neutrinos and up-type quarks, which couple to H u 10 , have Yukawa matrices with a universal structure as in Eq. 1, dictated by CSD3. The charged leptons and down-type quarks, which couple to H d 10 , have Yukawa matrices with a different universal structure where Y 11 is replaced by Y 12 ∼ φ 1 φ 2 T . Then quark mixing originates primarily in the down-type quark sector, with the down and strange quark masses successfully realised by having a zero entry in the (1,1) element of the down-type quark Yukawa matrix Y d , as in the GST approach [20], with a milder hierarchy among down-type quarks as compared to up-type quarks.
The model accurately fits all available quark and lepton data, and predicts the leptonic CP phase in 95% credible intervals given by 281 • < δ < 308 • and 225 • < δ < 253 • , with a significant deviation from maximal CP violation. Since quark mixing dominantly originates from Y d , analytical estimates for the quark mixing angles are given, revealing some tension in the predicted observables, which can be ameliorated by assuming rather large SUSY threshold corrections. A hierarchy in the flavon VEVs fixes the scales of all but one parameter, with all dimensionless couplings in the renormalisable theory naturally O(1). The model reduces to the MSSM, and we demonstrate how a µ term of O(TeV) can be realised, as well as doublet-triplet splitting, with Planck scale proton decay operators suppressed. In order to achieve the above we also require an auxiliary Z 2 4 and Z R 4 symmetry and a spectrum of messenger fields.
We would like to emphasise that the model presented here is very different from earlier models based on S 4 × SO(10) [21][22][23][24] (see also [25][26][27]). 6 Firstly, the full symmetry is different, since we invoke an extra Z 2 4 ×Z R 4 symmetry, while earlier works use a Z n [22][23][24]. Furthermore, we only allow small Higgs representations 10 (fundamental), 16 (spinor) and 45 (adjoint) and do not allow the large Higgs representations such as the 126 and 120 which are used in the other approaches. As a consequence our neutrino masses follow from a type-I seesaw mechanism, rather than a type-II seesaw employed in other papers. In further contrast, we do not allow Higgs mixing: the MSSM Higgs doublets h u and h d emerge directly from H u 10 and H d 10 , respectively, whereas in [21][22][23][24] they arise as unconconstrained linear combinations of doublets contained in 10-and 126-dimensional Higgs fields. In addition we consider doublet-triplet splitting. These features are largely absent from earlier works.
Another important difference is that we have used the CSD3 vacuum alignments in [16], whereas the vacuum aligments used in most previous works were geared towards TB mixing, and do not naturally provide a large reactor angle. For instance, in [22] they predict U e3 ∼ 0.05, which is now excluded. Indeed our model is motivated by the success of CSD3 and LS in the neutrino sector, although the LS predictions will be slightly modified due to charged lepton mixing corrections which are necessarily present due to the relation between charged leptons and down-type quarks. We should also mention that the first attempts to unify the three families as ψ ∼ (3,16) with Yukawa matrices resulting from flavons with particular vacuum alignments were in [31].
The remainder of the paper is laid out as follows: in Section 2 we describe the symmetries and field content of the model, the diagrams which yield the fermion mass matrices, and discuss proton decay. In Section 3 we define the mass and Yukawa matrices, derive analytical estimates for quark mixing parameters, and perform a detailed numerical fit to data. Section 4 concludes. In appendices we describe the mechanism giving doublettriplet splitting and a µ term, and the correct treatment of the third family Yukawa couplings, as well as a detailed derivation of the mass matrices.

The model
We present a model with quarks and leptons unified in a single ψ in the (3 , 16) representation of S 4 × SO (10), and with H u,d 10 in (1,10) and φ i in (3 , 1) representations. The idea is that the up-type quark Yukawa matrix Y u and neutrino Yukawa matrix Y ν arise from effective terms like where the group contraction in each bracket is into an S 4 singlet. These nonrenormalisable operators will have denominator scales of order M GUT , determined by the VEVs of additional Higgs adjoint 45s, leading to various CG factors. This leads to the Yukawa matrices Y u and Y ν being a sum of rank-1 matrices as in Eq. 1, with independent coefficients multiplying each rank-1 matrix, where Y ij ∼ φ i φ j T , up to S 4 CG factors. In the present model we assume the CSD3 flavon vacuum alignments [16], with VEVs driven to scales with the hierarchy 7 so that each rank-1 matrix in the sum contributes dominantly to a particular family, giving a rather natural understanding of the hierarchical Yukawa couplings y u ∼ v 2 1 /M 2 GUT , y c ∼ v 2 2 /M 2 GUT , y t ∼ v 2 3 /M 2 GUT , and similarly for the neutrino Yukawa couplings. Since the expansion breaks down for the third family, in the complete model we shall find a renormalisable explanation of the third-family Yukawa couplings. The RH neutrino Majorana mass matrix will also have the same universal form, leading to the seesaw mass matrix as in Eq. 1.
The down-type quark Yukawa matrix Y d and charged lepton Yukawa matrix Y e arise from terms like introducing a mixed term involving φ 1 and φ 2 , leading to a new rank-2 Yukawa structure Y 12 ∼ φ 1 φ 2 T . This leads to the Yukawa matrices Y d and Y e having Y 11 replaced by Y 12 , which has two consequences: it enforces a zero in the (1,1) element of Y d , giving the GST relation for the Cabibbo angle, i.e. θ q 12 ≈ y d /y s , and also leads to a milder hierarchy in the down and charged lepton sectors. Both features are welcome.
We need additional symmetries and fields to ensure the above structures, provide renormalisable third family Yukawa couplings, give the desired CG relations to distinguish down-type quarks from charged leptons, achieve doublet-triplet splitting, and obtain the MSSM Higgs doublets h u and h d from H u 10 and H d 10 , respectively. Two Z 4 shaping symmetries help to forbid any mixed flavon Yukawa terms. We also assume a discrete R symmetry Z R 4 , under which the superpotential has total charge two, and which is broken at the GUT scale by the H B−L 45 VEV to Z R 2 , the usual R (or matter) parity in the MSSM, ensuring a stable lightest supersymmetric particle (LSP). It also controls the µ term and helps ensure that only two light Higgs doublets (and no Higgs triplets) are present in the effective MSSM. Z R 4 is the smallest R symmetry that can achieve the above, and is specially motivated within SO(10) [32]. We shall also assume a spontaneously broken canonical CP symmetry at the high scale.
The full superfield content of the model is given in Table 1. It contains the following: a "matter" superfield ψ containing all known SM fermions, three triplet flavons φ which acquire CSD3 vacuum alignments, two Higgs 10s containing one each of the electroweakscale Higgs SU (2) doublets, 8 a spinor H 16 which breaks SO(10) (and gives masses to the RH neutrinos), as well as several Higgs adjoints. The χ superfields are messengers that are integrated out below the GUT scale, and are given GUT-scale masses by the VEV of Matter, Higgs and flavon superfields.

Field
Representation At the GUT scale, the renormalisable Yukawa superpotential is given by where we sum over indices a = 1, 2, 3 and b = 2, 3, and have suppressed O(1) coefficients λ that multiply each term. Furthermore, there are several crucial terms that appear suppressed by one Planck mass M P . These are where a = 1, 2, 3. The first term couples H 16 to fermions via the messengers. The second is allowed by the symmetries and will be shown to contribute at the order of the smallest GUT-scale terms to the fermion Yukawa matrices, and thus cannot be ignored. gains a VEV in the direction that preserves B − L, generating GUT-scale masses for Higgs triplets via the Dimopoulos-Wilczek (DW) mechanism [33]. Our implementation of the DW mechanism is described in Appendix A.
The VEVs of φ 1 and φ 2 are assumed to acquire VEVs well below the GUT scale, i.e. φ 1,2 M GUT , while φ 3 ∼ M GUT , which is therefore also the scale at which the flavour symmetry is broken, along with CP . We note that no residual CP symmetry remains at low scales, but CP does play a role in fixing phases in the mass matrices. As φ 3 is near the messenger scale, the process of integrating out messengers χ 3 , χ 3 is not trivial. The correct procedure and the consequences of having a flavon VEV near M GUT are discussed in detail in Appendix B, where we verify also that the third family Yukawa couplings are renormalisable at the electroweak scale.
The diagrams giving the mass and Yukawa matrices are drawn in Figs. 1-3. 9 The three diagrams in Fig. 1 correspond to the ultraviolet completion of the three terms in Eq. 2, while those in Fig. 2 are the completion of the terms in Eq. 5. The diagrams ensure correct S 4 group theory contractions and introduce CG coefficients due to the H X,Y,Z

45
VEVs. These diagrams are analogous to how the seesaw mechanism replaces the Weinberg operator for neutrino mass. Of course neutrino mass itself in this model is more subtle, since both the Dirac and RH Majorana masses arise from these diagrams.
Each diagram leads to a 3 × 3 matrix, whose internal structure is dictated by the vacuum alignment of the relevant flavon VEVs in Eq. 3. The Yukawa and mass matrices are consequently given as a sum over these matrices. A prominent feature is a texture zero in the (1,1) element of Y d and Y e , which realises the GST relation for the Cabibbo angle. The exact matrices that we fit to data are given in Section 3, with a full derivation in Appendix C. Further Planck-scale operators suppressed by one power of the Planck mass M P are forbidden by the symmetries. However we expect additional effective operators arising in the model, suppressed by at least two powers of the Planck mass M 2 P . These include terms involving all possible contractions of S 4 multiplets ψ and φ i , which are forbidden at the renormalisable level, but allowed by the symmetries. The largest of these terms Figure 2: Diagrams coupling ψ to H d 10 . When flavons acquire VEVs, these give the down-type quark and charged lepton Yukawa matrices. Figure 3: Diagrams coupling ψ to H 16 . One copy of the right diagram may be drawn for each of a = 1, 2, 3, although for a = 3, its contribution is negligible compared to the left diagram. When flavons acquire VEVs, these give the RH neutrino mass matrix.
can be O(M 2 GUT /M 2 P ) ∼ 10 −6 . We will assume these contributions are negligible, but note that such corrections may pollute the texture zero in Y d .
An adjoint of SO (10) can acquire a VEV aligned in the direction of any of the four U (1) subgroup generators that commute with the SM, or a combination thereof. 10 The VEVs of H X,Y,Z 45 may be written as linear combinations of these alignments. Fermions couple to these VEVs with strengths that depend on their associated U (1) charges, which are different for quarks and leptons.
Up-type quarks and Dirac neutrinos couple to H Z 45 (see Fig. 1). As H Z 45 is arbitrary, there is no hard prediction for the ratio between quark and neutrino Yukawa couplings within a family. However, as all flavons φ a couple to this VEV in the same way, flavour unification demands that the same ratio hold for all families. Therefore, once Y u is determined, Y ν is also fixed, such that Y ν ∝ Y u , to good approximation, up to an overall CG factor, with small deviations for the third family. While the elements of Y ν are not accessible at low energy scales (e.g. in neutrino oscillations), neutrino Yukawa couplings may be probed by considering leptogenesis. Since the lightest right-handed neutrino N 1 is too light, we anticipate that N 2 thermal leptogenesis will be required, where the secondto-lightest right-handed neutrino in our model has a mass of O(10 10 ) GeV, which is in the preferred range.
Meanwhile, the down-type quarks and charged leptons couple to two adjoints H X 45 and H Y 45 (see Fig. 2). Unlike the up sector, where matter always couples to the same SO(10) VEV, each diagram like Fig. 2 involving a different flavon will couple to a different linear combination of VEVs. This introduces CG factors non-trivially into Y d and Y e . As such, there is no fixed relationship between down-type quark and charged lepton Yukawa couplings, neither within a family, nor across families. They are nevertheless expected to be of the same order.
One of the characteristic features of GUTs is the prediction of proton decay. It has not been observed and the proton lifetime is constrained to be τ p > 10 34 years [1]. Proton decay can be mediated by the extra gauge bosons and by the triplets accompanying the Higgs doublets. In SUSY SO(10) GUTs the main source for proton decay comes from the triplet Higgsinos. The decay width is dependent on SUSY breaking and the specific coupling texture of the triplets. In general the constraints are barely met when the triplets are at the GUT scale [35] (as shown in Appendix A, this is our case).
The existence of additional fields in the model may allow proton decay from effective terms of the type Such terms must obey the constraint g X < 3 × 10 9 GeV. In our model, the largest contribution of this type comes from the term The constraint on X is easily met, so proton decay from such terms is highly suppressed.
3 Mass matrices, estimates and fits

Mass matrices
We present here the Yukawa and mass matrices, which will be used in the numerical fits below. The full derivation of these matrices, taking into account S 4 products, CG factors and third-family mixing, is given in Appendix C. We begin by defining numerical matrices We note that all matrices derive from triplet products like (ψφ i )(ψφ j ), with S 4 singlet contractions in each bracket, except Y P which derives from the Planck-suppressed operator ψψφ 3 H d 10 .
The up, down, charged lepton and Dirac neutrino Yukawa matrices (Y u , Y d , Y e and Y ν respectively) and RH neutrino mass matrix M R arising from Figs. 1-3, assuming that the MSSM Higgs doublets h u and h d arise from H u 10 and H d 10 , respectively, as shown in Appendix A, may then be expressed, as shown in Appendix C, as Y e = y e 12 e i η 2 Y 12 + y e 2 e iαe Y 22 + y e 3 e iβe Y 33 + y P e iγ Y P .
The flavon VEVs v a are complex, with the fixed phase relation given (up to a sign) by the superpotential that fixes the alignments. The remaining phase η is determined by the fit.
The light neutrino mass matrix is obtained by the seesaw mechanism. Both Y ν and M R have the same structure, namely both are sums over the same rank-1 matrices Y 11 , Y 22 and Y 33 . By a proof given in [36], the light neutrino matrix m ν will also have this structure, i.e.
where the parameters µ i are given in terms of the parameters y ν i and M R i simply by As shown in the introduction, the flavons yield a light neutrino mass matrix m ν , where the normal hierarchy m 1 m 2 m 3 then corresponds to µ 3 µ 1 µ 2 . Achieving this hierarchy after seesaw implies that the RH neutrino masses are very hierarchical, as we will see below. 11

Analytic estimates
The mass matrices involve the following real free parameters: y u i , y d i , y e i , µ i , and y P (a total of 13). Recalling that η is fixed by flavon vacuum alignment, we have the following further free parameters: η , α d,e , β d,e , and γ (a total of 6). The scales of the real parameters are mostly fixed by the scales of the flavon VEVs, v 1,2,3 . We set the flavon VEV scales to some appropriate values, where we set M GUT 10 16 GeV. The terms giving M R 1,2 and y P in Eqs. 13 and 14-15, respectively, derive from terms suppressed by one Planck mass M P . As they arise 11 While the model does not mathematically forbid an inverted hierarchy, we have checked that the corresponding predictions for neutrino masses and mixing angles would always give a bad fit to data. It would also require parameter choices that strongly violate the naturalness principle employed here. from unspecified dynamics, the scale of these parameters is not very well defined. For definiteness, we set M P 10 19 GeV and again assume the associated coefficients are close to one. We assume that M ρ ∼ M GUT .
We may estimate the parameters of the matrices defined in Eqs. 11-15 as follows: set all O(1) coefficients to exactly one, and ignore CG factors by setting all adjoint Higgs VEVs to M GUT 10 16 GeV. Then the Yukawa couplings are estimated to be The RH neutrino mass parameters are estimated to be This very strong hierarchy implies negligible RH neutrino mixing, such that the mass eigenvalues closely correspond to the above values. As each parameter contains several O(1) coefficients λ and CG factors, the above numbers only represent order of magnitude estimates.
As we will see in the numerical fit below, the above estimates are in good agreement with the values that produce a good fit to data, with a single exception: the parameter M R 1 , which is primarily responsible for the lightest RH neutrino mass, should be a factor O(0.01) times the estimate above in order to give the correct light neutrino mass spectrum. This can be understood by inserting the above estimates for y ν 1 and M R 1 into the expression for µ 1 in Eq. 18, which suggests µ 1 ∼ 0.01 meV, whereas we will see the fit prefers a value of O(1) meV. The necessary factor can be achieved by assuming one or more coefficients deviates from unity.
One may also obtain approximate expressions for the quark mixing angles in terms of quark Yukawa couplings as follows. The very strong hierarchy in the three real parameters of Y u is correlated with that in the physical Yukawa eigenvalues of up, charm and top quarks. We therefore expect negligible contributions from the up sector to quark mixing. This implies that not only do the four real parameters in the down sector, y d i and y P , fix the down-type Yukawa eigenvalues, they also must reproduce the observed CKM mixing angles.
Let us consider Y d , keeping only the leading terms in each element. For simplicity, we ignore free phases. As noted above, y d 12 ∼ y P < y d 2 y d 3 . We also define y 2 = y d 2 + 2y d 12 + 2y P . Then In the small angle approximation, the mixing angles can be estimated by The down-type Yukawa eigenvalues are given by y d ≈ (y d 12 ) 2 /y 2 , y s ≈ y 2 , y b ≈ y d 3 . Solving for y d 12 , y 2 and y d 3 , we have, to good approximation, y d 12 ≈ √ y d y s , y 2 ≈ y s , y d 3 ≈ y b . Reintroducing these into our estimates for mixing angles, we get Note that the first equality is exactly the GST relation [20], which is in good agreement with data. In fact, the GST relation, which predicts θ q 12 0.224 for the central values of y d and y s , is in mild tension with experimental data, which gives θ q 12 0.227. Possible modifications to the GST result have been proposed [37], e.g. adding a correction like y u /y c , which can be realised by a texture zero also in Y u . Alternatively, one may exploit the statistical uncertainties on each of the down and strange quark masses. A small deviation from their central values can predict a slightly different θ q 12 .
On the other hand, the mixing angles θ q 13 and θ q 23 are less precisely estimated, as the parameter y P can be as large as y d 12 , and the final result will depend on the relative phase between y d 12 and y P . Note however that both mixing angles depend in the same way on y d 12 − y P . Generally, the approximations in Eq. 24 predict some tension between θ q 13 and θ q 23 , which are too large and too small, respectively. This tension cannot be resolved simply by tuning y P .

Numerical fit
Our model determines the Yukawa couplings and mixing parameters at the GUT scale, which is also the highest flavour-breaking scale. The values from experiments must therefore be run up to the GUT scale. Moreover, when matching the SM to the MSSM at the scale M SUSY , supersymmetric radiative threshold corrections have to be included. Such an analysis has been performed in [38]. The parametrisation of these corrections is summarised in Appendix D. Most parameters do not significantly affect the fit, so are simply set to reasonable values. Specifically, we set M SUSY = 1 TeV, tan β = 5 andη q =η = 0. We also find that a good fit can be achieved for a rather large valueη b = −0.8. The choices of SUSY parameters tan β andη b are here empirically determined to give a good fit of the model to data. It is clear from the fit that large (negative)η b is required, affecting primarily the bottom quark Yukawa coupling y b . In order to keep y b perturbative, we must assume reasonably small tan β. In the region of 5 < tan β < 10 or so, the fit is rather insensitive to the exact choice. Neutrino data is taken from the NuFit global fit [7].
To find the best fit of the model to data, we minimise a χ 2 function, defined in the standard way: for a given set of input parameters x, we calculate the n observables P n (x). These are then compared to the observed values P obs n , which have associated statistical errors σ n . 12 Then For our model, the input parameters are x = {y u i , y d i , y e i , y P , µ i , η , α d,e , β d,e , γ}, and the observables are given by P n ∈ {θ q ij , δ q , y u,c,t , y d,s,b , θ ij , y e,µ,τ , ∆m 2 ij }. Note that as the lepton CP phase δ is not yet well measured, we do not include it in the fit, rather we prefer to leave it as a pure prediction. Furthermore, only the neutrino mass-squared differences are measured in oscillation experiments (as opposed to the masses themselves), while our model predicts the masses outright, including the lightest neutrino mass m 1 .  Table 2: Model predictions in the lepton sector for tan β = 5, M SUSY = 1 TeV andη b = −0.8. The observables are at the GUT scale. The lepton contribution to the total χ 2 is 0.03. δ as well as the neutrino masses m i are pure predictions of our model. The model interval is a Bayesian 95% credible interval. The bound on m i is taken from [6].
We present the best fit (minimum χ 2 ) of the model to physical observables (Yukawa couplings and neutrino mass and mixing parameters) in Tables 2 and 3, which also include the central values and 1σ ranges from data. Fig. 4 shows the associated pulls, and Table 4 shows the corresponding input parameter values. The fit gives χ 2 ≈ 3.4. 13 A second minimum with χ 2 ≈ 4 was also found, leading primarily to a different prediction for δ , as discussed below, although we shall not present the full fit parameters for this case. 12 In order for a minimum χ 2 to correspond to the maximum likelihood, the statistical uncertainties should be symmetric (Gaussian). This is essentially satisfied for all parameters except θ 23 , where current experimental data cannot conclusively resolve the octant, i.e whether it is larger or smaller than 45 • . Currently, the data favours θ 23 < 45 • , with a central value 41.6 • [7]. We will assume this is the true value. 13 The best fit predicts a strong neutrino hierarchy, with m 1 < 1 meV. It is possible to achieve a milder hierarchy, although the numerical fit gives χ 2 20 in such cases, predicting neutrino masses of approximately 5, 10 and 51 meV. Additionally it predicts δ ≈ +25 • , currently disfavoured by experiment.   We see from Tables 2, 3 and Fig. 4 that both quark and lepton sectors are fitted to within 1σ of the values predicted by global fits to experiment. The biggest pulls are in downtype quark Yukawa couplings y d,s,b and θ q 23 . As shown in Section 3.2, θ q 23 is approximately given by the ratio y s /y b , which is typically too small. Furthermore, attempts to increase θ q 23 , e.g. by tuning y P , tends to increase θ q 13 , which is then too large. This tension can be ameliorated by assuming large threshold corrections, i.e. by settingη b = −0.8, although some tension remains among the above parameters, which deviate by about 1σ. Tables 2 and 3 also include a Bayesian 95% credible interval for each observable, 14 which was found by performing a Markov Chain Monte Carlo (MCMC) analysis. The interval for a given parameter corresponds to the region of highest posterior (probability) density (hpd), marginalised over the other parameters, and may be interpreted as follows: given the data, there is a 95% probability that the true model value of that observable resides in the stated interval. For many observables, these probability distributions are essentially Gaussian, centred around the best fit value. This is not always the case: the distributions for θ 12 and θ 23 are asymmetric, consisting of two partially overlapping peaks. Moreover, the hpd region for δ consists of two completely distinct intervals, which contain the best fit values 300.9 • (as seen in Table 2) and 233.9 • (corresponding to a second best fit point with χ 2 ≈ 4). Their associated 95% credible intervals are given by 280.7 < δ < 308.3  Table 4: Best fit input parameter values. The model has 13 real parameters: y u i , y d i , y e i , µ i and y P . While η is fixed by flavon alignment to −2π/3, there are six additional free phases: η , α d,e , β d,e and γ. The total χ 2 is 3.4. and 225.1 < δ < 253.2, respectively. We note that neither region includes maximal CP violation δ = 270 • , which is close to the prediction from CSD3 with diagonal charged leptons. In short, charged-lepton corrections induce a deviation from maximal CP phase, which can either be positive or negative, depending on the phases of Y e .
One may be tempted to calculate a reduced chi-squared χ 2 red , i.e. the χ 2 per degree of freedom (d.o.f.), where the number of d.o.f. is naively given by the number of observables minus the number of input parameters. In the conventional picture, a good fit has χ 2 red 1. However, as discussed in [39], this interpretation is only valid for linear models, which our model is not. Indeed, when evaluating χ 2 we fit 19 inputs to 18 observables, which in a linear model would suggest a perfect fit is always possible; this is certainly not the case. While χ 2 is a valid tool for comparing models to each other, since it is not possible to establish an exact number of d.o.f., we cannot reliably define χ 2 red .

Conclusion
The flavour puzzle in the SM is the source of a majority of the SM free parameters, characterised by different mixing behaviours for quarks and leptons, and very hierarchical masses. The most minimal solution to the problem of neutrino masses remains the seesaw mechanism with heavy RH neutrinos, which arise automatically in SO (10), with naturally large masses. This motivates SO(10) above other popular gauge groups, such as SU (5), where RH neutrinos are added by hand. All three families of SM fermions in the 16 of SO(10) are here also unified in a single triplet of S 4 . This very elegant picture presents model-building challenges, many of which we have tackled in this paper.
We have constructed a rather simple, natural and complete SO(10) model of flavour with a discrete S 4 × Z 2 4 × Z R 4 symmetry, where all Yukawa matrices derive from the VEVs of triplet flavons, in the CSD3 alignment. It is simple in the sense that the field content is reasonably minimal, with small Higgs representations of SO(10) consisting of two 10s which contain the MSSM doublets, a Higgs spinor pair 16 and 16 responsible for Majorana masses and four adjoint Higgs 45s, which provide necessary Clebsch-Gordan factors that distinguish charged leptons and down-type quarks. It is natural in the sense that Yukawa and mass matrices consist of sums of low-rank matrices, each of which contributes dominantly to a particular family, i.e. "universal sequential dominance". It is complete in the sense that we address the µ-problem, Higgs mixing and doublettriplet splitting, and provide an ultraviolet renormalisable model, with Planck-suppressed operators controlled by symmetry. However, we do not discuss the origin of the hierarchy of flavon VEVs, nor do we repeat the discussion of flavon vacuum alignment, which can be found in [16]. We believe this model represents a signficant step forward in the quest for a complete and correct description of fermions within SUSY GUTs. For instance, we have demonstrated the correct procedure for treating the third family couplings and how to generate an electroweak-scale renormalisable third-family Yukawa coupling. We also emphasise that the principle of universal sequential dominance is a simple and effective way to understand fermion hierarchies. Although the origin of such family hierarchies has not been fully resolved, as the scales of flavon VEVs v i are assumed rather than proven, the problem has been ameliorated, since the hierarchy is given by the squares of these VEVs.
The model successfully reproduces the observed fermion masses and mixing, even in the quark sector, where the CKM parameters are measured to very high precision. Analytical estimates are underpinned by a detailed numerical analysis, demonstrating the viability of the model. Moreover, there is no tuning of O(1) parameters necessary to explain the mass hierarchies of charged fermions, accounting also for the milder hierarchy in down-type quarks compared to up-type quarks. The model simultaneously realises large lepton mixing and small quark mixing, as well as the GST relation for the Cabibbo angle, θ q 12 ≈ y d /y s via a texture zero in the down-type Yukawa matrix Y d . In the lepton sector an excellent fit to data is found, predicting a normal neutrino hierarchy and lightest neutrino mass m 1 0.5 meV. The CP phase δ was not fitted, but left as a pure prediction. Two distinct regions are preferred, with corresponding best fit values δ ≈ 301 • and 234 • . We emphasise that the model predicts significant deviation from both zero and maximal CP violation.

A Doublet-triplet and doublet-doublet splitting
As is the case for every broken GUT, the Higgs sector of our model contains more fields than the usual MSSM. The H u,d 10 multiplets contain colour triplets that mediate proton decay. Since we have two 10s, there is an additional pair of doublets that, if light, could spoil gauge coupling unification. For these reasons, those extra fields need to be heavy, while ensuring the MSSM doublets are massless. This splitting can be achieved in our model.
The splitting mechanism involves superfields given in Table 1. The singlet ξ obtains a VEV slightly above the GUT scale and ensures the correct structure to the masses. The H 16 generates a mass for the H 16 and also gets a VEV in the RH neutrino (ν c ) direction.
is the only R-charged field that gets a VEV, breaking Z R 4 to the usual R parity. This splitting mechanism needs three extra messenger pairs which are listed in table 5.

Field
Representation where we assume that the VEV ξ M GUT , so that we may integrate out the messenger fields and obtain effective superpotential where we have suppressed dimensionless couplings, and the final term involves all combinations of adjoints allowed by the symmetries, i.e. either (H Z 45 ) 4 or any combination of powers of H X 45 and H Y 45 totalling four. The three terms suppressed by ξ are allowed by the integration of three messenger pairs.
We assume that the superfields H 16,16  gets a VEV aligned in the B − L direction, which splits doublet and triplet Higgs masses through the Dimopoulos-Wilczek (DW) mechanism [33]. This can be understood by considering the decomposition of the H u,d 10 into the Pati-Salam group. The triplets behave as a sextuplet of SU (4) while the doublets are singlets. Since U (1) B−L ⊂ SU (4), the triplets get a mass from the first term of Eq. 27 while the doublets don't. In the last term, all the SO(10) adjoints can be contracted to a singlet, so they affect doublets and triplets equally.
To demonstrate the mechanism, we construct the doublet and triplet mass matrices. We define some dimensionless scale parameters y = M GUT /M P , z = M GUT / ξ . We label the up-type doublets inside a given Higgs representation H by 2 u (H), and down-type doublets by 2 d (H). We define triplets 3 u (H) and 3 d (H) analogously. H can be either H u 10 , H d 10 or H 16,16 . The mass matrices M D and M T are given by The triplets mass matrix M T has three eigenvalues of O(M GUT ). The doublets mass matrix has two eigenvalues at O(M GUT ) and one at O(y 4 M GUT ), which we identify with the µ term. Since y ≈ 10 −3 we have µ ∼ 1 TeV, which is the desired order. Furthermore, the light eigenvectors of M D define the MSSM doublets h u,d as where the contribution of O(y) is negligible, so that the MSSM doublets are located as required by the Yukawa structure of the model.

B Renormalisability of the third family
In this section we show that naive integration over messenger fields is not possible for the third family, due to the large VEV of φ 3 . We reiterate that there is an assumed hierarchy of flavon VEVs, such that v 1 v 2 v 3 ∼ M GUT , implying it is not possible to formally integrate out the messengers χ 3 which couple to the flavon φ 3 .
To explore this further, let us single out the terms in W Y involving these fields and H u 10 (the same method applies to terms coupling to H d 10 ). Suppressing O(1) couplings, the relevant terms are W After fields acquire VEVs (with φ 3 = v 3 (0, 0, 1)), we have These two terms are of comparable order.
Naively, ψ 3 may be interpreted as the set of third-family particles. The problem with this picture is that it has a large coupling to χ 3 , which induces a mass for ψ 3 via the second term in Eq. 31. This clearly does not correspond to the physical third-family states (top quark and third Dirac neutrino), which are massless above the electroweak scale. To obtain the physical (massless) states, which we label t, we rotate into a physical basis (ψ 3 , χ 3 ) → (t, χ), such that t does not couple to χ 3 . This basis change is given by Physically, it may be interpreted as follows: inside the original superpotential W Y lie the terms which generate renormalisable mass terms for the top quark and the third Dirac neutrino at the electroweak scale.

C Complete derivation of mass matrices
In this section we derive the precise forms of the Yukawa and Majorana mass matrices, and cast them in terms of a minimum number of free parameters, taking into account the vacuum alignments of the adjoint Higgs superfields as well as the induced rotation in the third family couplings due to a large φ 3 .
The renormalisable superpotential in Eq.6, including the dominant Planck-suppressed terms in Eq.7, and writing explicitly all O(1) couplings, becomes, The flavon vacuum alignments, obtained from [16], which preserve the generator product SU , are (as in Eq. 3) with the hierarchy v 1 v 2 v 3 . The singlet product which occurs in ψφ a above, i.e. 3 × 3 → 1, is given by To account for this nontrivial product as well as the field redefinition ψ 2 → −ψ 2 (this overall sign is unphysical), we define the vectors In the new variables, the alignments become, Next, we introduce notation to specify the relevant components of a VEV H k 45 , corresponding to unique CG factors. The index k labels the adjoint, i.e. k = X, Y, Z, or B − L. ψ couples to the adjoint VEVs via χ messengers. After the GUT is broken and ψ is decomposed into multiplets of the SM gauge group, the part of an adjoint VEV which couples to a given multiplet f is then denoted where f = Q, u c , d c , L, e c , or ν c . The H 16 gets a VEV in the direction which preserves SU (5), which we call the (singlet) ν c direction. Its VEV only affects the RH neutrino mass matrix and is simply denoted v 16 .
We extract the Yukawa matrices from diagrams in Figs. 1-3. Taking into account nontrivial S 4 products (as above), we have where v 3 = | φ 3 |, v u Y ν ij is the Dirac neutrino mass matrix and M R ij is the RH Majorana matrix, respectively. The last term in Eq. 34 is a singlet coming from three S 4 triplets and gives rise to the last terms in Eqs. 43 and 44, respectively, where Y P is the numerical matrix defined in Eq. 10.
Finally, we take into account the effect of mixing between the state ψ 3 and messenger χ 3 . This mixing provides additional contributions to the fermion mass matrices in the form of coefficients multiplying the third rows and columns. The size of each coefficient depends on the CG factors and the ratio(s) of v 3 to adjoint Higgs VEVs v k 45 , for k = X, Y, χ. In the limit where v 3 v k 45 , all these coefficients are 1, corresponding to a negligible amount of χ 3 being mixed into the physical state. This is exactly what occurs for the other two families: the massless states are aligned almost exactly with the states ψ 1,2 . Generally, any significant deviation would require a tuning among CG factors and O(1) parameters λ. We do not expect these factors to have a large effect on mixing, hence we set them all to one for simplicity.

D SUSY threshold corrections
An analysis of the running of MSSM Yukawa parameters up to the GUT scale has been performed in [38], where they propose a useful parametrisation of tan β-enhanced corrections to the charged fermion Yukawa couplings and quark mixing angles.
The matching conditions at the SUSY scale M SUSY are parametrised in terms of four parametersη q,b, andβ, as  to a good approximation only onη b and tanβ. In the limit where threshold corrections to y τ are negligible, β reduces to the usual β. We will assume just such a scenario. We will also setη q =η = 0 for simplicity, as these are found not to affect the quality of the fit. The associated shift in down-type and charged lepton Yukawa couplings can be largely subsumed into available free parameters y d,e 12 y d,e 2 . Similarly, we fix M SUSY = 1 TeV. Slightly larger values are allowed, up to O(10) TeV, but the effect on the fit is minor. Conversely, the remaining SUSY parameter,η b , will be important and prefers a large (negative) value. A very good fit is found forη b = −0.8, for which results are presented in this paper. Meanwhile, the neutrino masses and mixing angles are expected to be largely insensitive to group running.