The Interplay Between GUT and Flavour Symmetries in a Pati-Salam x S4 Model

Both Grand Unified symmetries and discrete flavour symmetries are appealing ways to describe apparent structures in the gauge and flavour sectors of the Standard Model. Both symmetries put constraints on the high energy behaviour of the theory. This can give rise to unexpected interplay when building models that possess both symmetries. We investigate on the possibility to combine a Pati-Salam model with the discrete flavour symmetry $S_4$ that gives rise to quark-lepton complementarity. Under appropriate assumptions at the GUT scale, the model reproduces fermion masses and mixings both in the quark and in the lepton sectors. We show that in particular the Higgs sector and the running Yukawa couplings are strongly affected by the combined constraints of the Grand Unified and family symmetries. This in turn reduces the phenomenologically viable parameter space, with high energy mass scales confined to a small region and some parameters in the neutrino sector slightly unnatural. In the allowed regions, we can reproduce the quark masses and the CKM matrix. In the lepton sector, we reproduce the charged lepton masses, including bottom-tau unification and the Georgi-Jarlskog relation as well as the two known angles of the PMNS matrix. The neutrino mass spectrum can present a normal or an inverse hierarchy, and only allowing the neutrino parameters to spread into a range of values between $\lambda^{-2}$ and $\lambda^2$, with $\lambda\simeq0.2$. Finally, our model suggests that the reactor mixing angle is close to its current experimental bound.


Introduction
Gauge coupling unification suggests that the Standard Model (SM) gauge group is generated when a unified larger gauge group is broken at a very high energy scale compared to the electroweak (EW) one. In this picture SM fermions are accommodated in representations of the unified gauge group G and an appropriate scalar Higgs sector is introduced both to trigger the spontaneous breaking of G down to SU (3) C × U (1) em and to reproduce the fermion mass matrices. Typically Grand Unified Theories (GUTs) are quite constrained and non-trivial relations are obtained among the SM fermion mass matrices. On the other hand it is a common belief that fermion mass and mixing matrix patterns can be explained by appealing to a flavour symmetry G f and in the last years particular attention has been devoted to the study of discrete flavour symmetries, thanks to their simplicity in recovering realistic lepton mixing patterns.
The present neutrino oscillation data [1][2][3] are summarised in table 1, where we display the results of two independent global fits. The pattern of the mixings is characterised by two large angles and a small one: θ l 23 is compatible with a maximal value, but the accuracy admits relatively large deviations; θ l 12 is large, but about 5σ far from the maximal value; θ l 13 has only an upper bound. According to the type of the experiments which measured them, the mixing angle θ l 23 is called atmospheric, θ l 12 solar and θ l 13 reactor. We underline that there is a tension among the two global fits presented in table 1 on the central value of the reactor angle: in [1] we can read a suggestion for a positive value of sin 2 θ l 13 0.016 ± 0.010 [1.6σ], while in [2] a best fit value consistent with zero within less than 1σ is found. Therefore we need a direct measurement by the future experiments like DOUBLE CHOOZ [4], Daya Bay [5], MINOS [6], RENO [7],T2K [8] and NOvA [9].
The closeness of the leptonic atmospheric angle θ l 23 to the maximal value gives relevant indications on the flavour symmetry: it is well known [10,11] that a maximal θ l 23 is not achievable with an exact realistic symmetry. This forces to study models based on the breaking of the flavour symmetry and a promising choice is the kind of realizations based on non-Abelian discrete groups which reproduce the lepton mixing matrix of the so-called the data; the perturbations are naturally constrained to get the "weak" complementarity relation, θ l 12 + O(λ) ∼ π/4, and sin θ l 13 ∼ λ in most of the parameter space. The model only deals with the lepton sector and an extension in the quarks sector is lacking. In this paper we aim at revisiting the model in [23] in order to include a realistic description of quarks.
Our starting point lies in the complementarity relations: For the third mixing angles we know that in the quark sector it is very small, θ q 13 = O(λ 3 ), while in the lepton sector, as already mentioned, it has only an upper bound, sin θ l 13 λ.
We will see that in the model described below, we predict it to be O(λ).
In non-GUT contexts, no compelling model leading to the exact complementarity has been produced so far and indeed in [23] a weaker version of eq. (1) has been used in which the quark contributions are substituted by similar terms originating from the charged lepton sector. On the other hand, exact complementarity is possible in cases where the flavour symmetry group is combined with a GUT group, as in the Pati-Salam 1 context [24], where the charged lepton and the down-quark mass matrices are similar, Other popular GUTs are SU (5) and SO (10), which however are less appealing when trying to recover the QL complementarity relations. Indeed, in the minimal SU (5) [25] one has M e ∼ M T d and as result, a correction of order λ to the solar angle does not correspond to the Cabibbo angle of the CKM matrix.
On the other hand a reason to prefer Pati-Salam over SO (10) is related to the type-I and type-II See-Saw mechanisms. In these two GUT contexts, we expect both left-handed (LH) and right-handed (RH) neutrino Majorana mass matrices to be present. As a result, the effective LH neutrino mass matrix will get the contributions through the type-I as well as the type-II See-Saw mechanisms. In general, and this happens also in our proposal, this interplay introduces two mass scales and a highly non-trivial flavour structure for the effective neutrino mass matrix, which difficultly leads to a realistic description of the PMNS matrix. For this reason, a hierarchy between the two contributions is usually assumed. As we will see we can reproduce quark-lepton complementarity in our model, if the type-II See-Saw is dominant. This possibility has already been investigated in several flavour GUT models, for example in [26] in the context of the SO(10) GUT. However, in [27] it has been agued that the type-II dominance in the context of minimal SO(10) models is highly disfavoured. Even so, when restricting to particular supersymmetric parameter space, the type-II dominance could be possible both in the minimal and non-minimal SO(10) approaches, [28]. In the Pati-Salam context, there is much more freedom and the eventual dominance of one of the two contributions could be realized. In this paper we study the gauge Higgs potential and we verify that a type-II dominance can be justified, even if it puts strong constraints on the model building realization.
After a detailed analysis of the fermion phenomenology, we move to the study on the Higgs gauge sector, which is responsible for the gauge symmetry breaking steps to finally get SU (3) C × U (1) em . Thus we consider the renormalization group equations (RGEs) both for the gauge couplings and for the fermion masses and mixings from the cutoff of the theory down to the low-energy scale. Here we anticipate that the RGEs analysis is crucial: the model turns out to be viable only in a small region of the parameter space. On one side the gauge coupling RGEs analysis is deeply modified by the non-minimal Higgs field content required by the presence of the flavour symmetry G f and puts strong constraints on the scales of the model. On the other side the GUT flavour mass matrix structures interfere with each other, through the Yukawa RGEs, with non-negligible consequences for the neutrino phenomenology.
The paper is organized as follows. In section 2 we describe the flavour structure of the fermion mass matrices in order to recover realistic lepton and quark mixing matrices and fermion mass hierarchies. In section 3 we enter in the details of the model building construction, specifying the transformations of all the fields of the model under the gauge and flavour groups and discussing the flavon vacuum misalignment necessary to the flavour symmetry breaking chain. Afterwards we deal in section 4 with the study of the gauge Higgs potential and we analyze the constraints coming from the flavour symmetry on the scalar field content. In section 5 we perform the analysis of the renormalization group running of fermion masses and mixings and of the gauge coupling constants from the GUT scale down to the low-energy scale. In section 6 we report the phenomenological analysis after the running evolution. Finally in section 7 we conclude. In the appendices we report technical details.

Outline of the model
In this section we present a general discussion on the flavour structure of fermion masses and mixings resulting form the analysis of the complementarity relations. Eq. (1) suggests that the angles in the CKM and PMNS matrices may have a common origin. Looking at their definitions where V u , V d , U e , U ν diagonalize M u M † u , M d M † d , M e M † e and m ν respectively, we see that this common origin can be motivated in a GUT context where some relations among the mass matrices are present. In PS models, we have the following expression which links down-quarks and charged leptons U e ∼ V d ; we will see in a while how this enters in the model construction.
In [23] the two large lepton angles in the PMNS matrix arise only from the neutrino sector: in the first order approximation, both the solar and the atmospheric angles are maximal and after that some corrections from the charged lepton sectors lower the value of the solar angle in order to accommodate the data. In our realization one maximal angle originates from the charged leptons and the other from the neutrinos: afterwards some corrections are necessary to make the solar angle agree with the measurements and we require that these corrections come from the charged leptons. This is an other possible choice with respect to [23] and in the next section we will show that it is easily realized in our context.
In the approximation of small λ, the lepton mixing matrix can schematically be written as where R ij (θ) stands for a rotation in the (ij) plane of an angle θ. 2 As a result, after a suitable commutation of matrices, U e can naively be written as The CKM matrix is given in first approximation as Because of eq. (4), V d has the same structure as U e in eq. (6) and therefore we can obtain the CKM matrix as We see that the angles of the (23) and (13) rotations in V † u should be the opposite of those in V d , while the angles in the (12) sector should be different. We have schematically indicated this via the α coefficient. Analogous to eq. (6), we write V u as Moving to the explicit form of the mass matrices, the generic effective Majorana neutrino mass matrix m ν which is diagonalized by U ν as in eq. (5), U ν = R 12 In the charged lepton sector, we are looking for the most general form of the mass matrix, whose square M e M † e is diagonalized by the action of U e of eq. (6), Inverting eq. (11) we find in the limit m e → 0 where the dots stand for suppressed contributions. This can be obtained if M e is given by It is interesting to note that, moving to the basis of diagonal charged leptons and considering only the leading order terms, the neutrino mass matrix results to be of the BM type, already presented in the introduction. As already stated, the PS gauge structure implies the relation M e ∼ M d and therefore the down-quark matrix has a similar structure as in eq. (13). Looking at eq. (9) we find that also M u should have a similar structure, but here there is a slight difference due to the α coefficient.
In order to reproduce the result in eq. (8), i.e. exact cancellations in the (23) and (13) sectors but not in the (12) sector, we need to study the origin of the different rotations in (9). It is easy to see that in the diagonalization of the mass matrices in eq. (12), the (12) rotation originates from the second families, while the (13) rotation comes from the third families. We conclude that the up-quark mass matrix has the same form as (13), in which the third column is proportional to that of the down-quarks, while the second column is not.
As already stated in the introduction, we expect contributions on the effective neutrino mass matrix coming from both the type-I and the type-II See-Saw mechanisms. For the moment we just assume that the type-II is the only responsible of eq. (10) and we verify in section 4 that it is indeed dominating with respect to the type-I terms.
In the next section we enter in the details of the model building, explaining the origin of the mass matrices displayed above.

The flavour model building
In this section we define our framework. The model is based on the PS gauge symmetry, SU (4) C × SU (2) L × SU (2) R , present at high energy where a supersymmetric context is assumed. The complete flavour group G f is the same as in [23] given by the product of the following different terms: The group S 4 is the permutation group of four distinct objects, isomorphic to the group O which is the symmetry group of a regular octahedron. It has 24 distinct elements filled in five conjugate classes and therefore it has five irreducible representations, two one-dimensional denoted as 1 1 and 1 2 , one two-dimensional labelled as 2 and two threedimensional written as 3 1 and 3 2 . Here we recall the multiplication rules, while the Clebsch-Gordan coefficients, the explicit structures of the generators, the list of the elements and 3.1 The matter, Higgs and flavon content of the model In the PS context, the five matter multiplets of each family of the SM plus a RH neutrino and their superpartners are unified in only two supermultiplets: a LH and a RH ones as follows, The three copies of the LH supermultiplet are combined in the three-dimensional representation 3 1 of S 4 , while the three families of the RH supermultiplet are in 1 2 , 1 2 and 1 1 respectively. The fact that we can put different representations within one family in different representations of the family symmetry group is essential here. Note that this would not be possible in (minimal) SO (10), where all Standard Model particles are in one sixteen dimensional representation. The first two families are also charged under U (1) F N by two units. This suppresses their masses with respect to the third family ones. Further suppression of the first family with respect to the second is due to their different Z 4 charges. All the properties of the matter fields are summarized in table 2.
Matter Table 2: Transformation properties of the matter fields. Notice that the PS assignments should be read in agreement with SU (4) C × SU (2) L × SU (2) R .
Our model contains five flavon fields: two S 4 triplets (ϕ and ϕ ) that, because of their Z 4 charge, deal at LO only with the Dirac Yukawa couplings of quarks and leptons, and two fields, one singlet (σ) and one triplet (χ), that, by Z 4 , deal at LO only with the Majorana masses of neutrinos. The fifth flavon is the Froggatt-Nielsen messenger, which we indicate with θ. Their properties are shown in table 3. Under the continuous R-symmetry, the matter superfields transform as U (1) R = 1, while all the flavons are neutral.
Fermion masses and mixings arise from the spontaneous breaking of the flavour symmetry by means of the flavons which develop VEVs according to the following configuration: at LO we have In this section we simply assume this VEV alignment and we will prove it to be a natural solution of the minimization of the scalar potential in section 3.4. Furthermore we assume that the FN messenger and the other flavons have VEVs of the same order of magnitude: it results partly from the minimization procedure and partly from the constraints coming form the comparison with the measured mass hierarchies, as it will be clearer in the following. Without loss of generality it is useful to keep the notation compact and write where V EV refers to the vacuum expectation value of any flavon of the model and λ 0.2 is close to the Cabibbo angle. Λ is the flavon cut-off that we assume to be the largest scale in the model: it corresponds to the scale of the flavon dynamics. The insertion of flavons in the mass terms, leads to non-renormalizable operators, as we will see in the following section. This means that we can write the superpotential as an expansion in powers of flavon/Λ and we can stop the expansion after the first orders.
The Higgs fields of our model relevant to build the fermion mass matrices transform under the gauge group and under the Z 4 factor of the flavour symmetry: in table 4 we summarize their transformation rules. The first three fields, φ, φ and ρ, deal at LO only with the Dirac Yukawas. Due to the Z 4 charges, φ and φ are responsible of the third family and the charm quark masses, while ρ is responsible for the strange and µ masses. The field ρ ∼ (15,2,2) being in the adjoint of SU (4) C may develop VEV along the SU (4) C direction diag (−3, 1, 1, 1). This implies that the leptons which get mass via this field are a factor 3 heavier than the corresponding quarks and therefore this field is very useful to describe the second family, at least in the down sector, reproducing the well known Georgi-Jarlskog relation [30], m µ ≈ 3m s , at the high energy scale. As we will see in the next sections, in order to recover the up-quark mass hierarchies the ρ projection along the light doublet Higgses, v u ρ and v d ρ , has to point only in the down direction: the requirement v u ρ = 0 can be realized only if the Higgs field content contains two identical copy of the Higgs field (1,2,2) and this justify the presence of φ and φ .
Finally, as we will see in detail in the following sections, the field ∆ R is necessary to conclude the PS symmetry breaking pattern and to recover the SM gauge group through its spontaneous symmetry breaking VEV. At the same time, when ∆ R develops a VEV, it gives a Majorana mass to the right-handed neutrinos thus contributing to the effective neutrino mass matrix through the usual type-I See-Saw mechanism. A second source for the neutrino mass matrix comes from ∆ L , in terms of the type-II See-Saw mechanism.

PS
(1, 2, 2) (15, 2, 2) (10, 3, 1) (10, 1, 3) In our scheme the neutrino mass matrix is dominated by type-II See-Saw. As already stated in the introduction, the flavour structure of the effective neutrino mass matrix in eq. (10) arises from an interplay between the two See-Saw sources. In the PS context m D ∼ M u and this suggests a hierarchical structure for the type-I contribution, which does not agree with the flavour structure in eq. (10). We will show that, given m D ∼ M u , the required flavour structure for the Majorana RH neutrino mass necessary to recover eq. (10) is not allowed in our model. This suggests to find a construction in which the type-II contributions are dominating over the type-I and we will show that such a feature puts strong constraints on the model building.

Fermion mass matrices at leading order
In this section we study the matter superpotential, W Y , and the resulting mass matrices for all the fermions. The Yukawa superpotential can be written as a sum of two pieces, containing respectively terms with Dirac and Majorana structures. Furthermore, making a power expansion in terms of flavon/Λ, we distinguish between leading and subleading couplings. We refer as "leading order" operators to those terms which provide M e,d,u and m ν as in eqs. (10,13) in the limit λ → 0. The subleading orders corresponds to operators of higher dimensions. In this way it is easier to identify each entry of the mass matrices with the corresponding operator in the superpotential and furthermore it underlines the relevance of the subleading contributions.

Dirac mass terms
We first study the Dirac matter superpotential at LO which reads Here we use a compact notation to avoid the proliferation of coefficients: the term X (1,2) i represents a list of products defined as where each term represents all the different S 4 contractions which can be constructed with those flavons; furthermore we indicate with y 1 (φ + φ ) the combination y 1 φ and similarly for y 3, (i) and y 5 .
When the flavour symmetry is broken the model describes a non-minimal PS model in which the Yukawa couplings present a well defined structure. Then, when the PS gauge symmetry is broken to the SM gauge group φ, φ and the colour singlet component of ρ mix in four light Higgses, two up-type and two down-type, h u,d and h u,d . Thus at the EW scale φ, φ and ρ have non-vanishing projections to the light Higgs components that acquire a VEV. We will indicate these components as v u,d φ , v u,d φ and v u,d ρ . As already said, we need to impose that the ρ field has no projection along the up direction: v u ρ = 0. This can be realized because in terms of the light Higgs up-VEVs, v u where the matrix U is introduced in the appendix B. The constraint v u 2 = −U 13 /U 23 v u 1 can be imposed thanks to the freedom we have in the superpotential and in the soft potential. 3 Note that we could relax the condition v u ρ = 0 allowing a mild hierarchy between v u ρ and v d ρ , for example of order λ 2 , without affecting the final mass hierarchies, but in the following, for simplicity, we work under the assumption that v u ρ = 0. The final Dirac fermion mass matrices we get are given by where we used the compact notation i v u/d φ and absorbed all the non-relevant CG coefficients. Note that y 3 is the sum of the different y 3, (i) and that, by construction, y 1,2,3 can be considered complex coefficients with modulus of order 1. Note also that the different numerical coefficients between charged leptons and downquarks originate from to the presence of ρ instead of φ (φ ) in the superpotential. The operators which should give contributions to the first families (those proportional to y 4 and y 5 ) are vanishing, thanks to the special flavon VEV alignment. As a final comment, we are neglecting at this level of approximation the contributions to M LO e,d coming from the operators proportional to y 3 : these terms, which preserve the anti-alignment of the second and third entries of the second columns, are λ 2 suppressed with respect to the LO ones proportional to y 2 .
In order to identify the mass matrices in eqs. (25)-(26)- (27) with those in eqs. (10,13), we need to define the (complex) fermion masses as follows: Now we can comment on the masses as well as on the mixing matrices which can be driven at this approximation level. Notice that the top-quark Yukawa does not come from a renormalizable coupling and it presents the same suppression as the other third family fermion masses: as a result the dominance of the top-quark mass can be justified by the hierarchy between v u φ and v d φ . The mass matrices in eqs. (25)-(26)- (27) are diagonalized by a maximal rotation in the sector (23), i.e. U e = V d = V u = R 23 (π/4), while the fermion mass hierarchies are given by Furthermore, at the cutoff, we recover some relations among the masses of different fermions: the b − τ unification and the Georgi-Jarlskog relation [30] Finally we should comment of the relative value of the top and of the bottom masses: Notice that usally this ratio is proportional to tan β, the ratio between the up-and down-VEVs of the light Higgses, but this is not the case: indeed the light Higgses are combinations of φ, φ and ρ and therefore we may define tan β as

Majorana mass terms
We now discuss the part of the superpotential which contains the Majorana couplings. At LO 4 it is given by where we used the compact notation This superpotential is responsible for giving the following Majorana LH and RH neutrino mass matrices: where k 0,1,2 and z 1 are coefficients of order 1 and v L , v R are the VEVs of ∆ L,R respectively. In particular k 0,1,2 are defined in terms ofk 0 ,k 1,(i) andk 2,(i) : 4 Regarding the terms which contribute to M R , we consider at the LO only the first non vanishing term in the superpotential.
While M L corresponds to the type-II contribution to the effective neutrino mass matrix, M R provides a type-I term. Even at this approximation level, we note the tension between the two See-Saw contributions: m (type-II) ν ≡ M L presents a maximal rotation in the (12) sector, while it is easy to verify that m (type-I) ν ≡ m D M −1 R m T D shows a democratic structure on the (23) sector which corresponds to a maximal rotation in this sector. To recover the mass matrix in eq. (10) it is necessary that the type-II contribution dominates over the type-I terms: in this case it is sufficient to identify a with k 0 , b with k 1 λ and c with k 0 + k 2 λ 2 .
The neutrino masses can be written as with θ k i the argument of the complex number k i . We need cos(θ k 0 − θ k 1 ) > 0 in order to have |m 1 | smaller than |m 2 |. We see that in most of parameter space the spectrum is quasi degenerate, as the term with |k 0 | 2 that appears in all three masses dominates over the other terms that are λ or λ 2 suppressed.
We also see that in most of the parameter space |m 3 | is the central eigenvalue: indeed L by a term proportional to λ, while |m 3 | stays closer to this central value as it is only shifted by a term proportional to λ 2 . Having |m 3 | as the central eigenvalue is obviously in contradiction with experimental data. We conclude that our model is only viable when the term 2|k 0 | |k 1 | cos(θ k 0 − θ k 1 )λ is suppressed. This is clearly possible if either |k 0 | or |k 1 | or cos(θ k 0 −θ k 1 ) (or a combination of them) is small. In particular, the latter condition means that k 0 and k 1 are almost perpendicular in the complex plane. We now investigate on these different scenarios, calculating the solar and atmospheric mass squared differences: taking as definition of ∆m 2 atm the mass squared difference between the heaviest and the lightest neutrinos, we have different results for normal (NO) and inverse (IO) mass ordering as given by On the right-hand side of the last equation, the first term is suppressed by λ, while the second one by λ 2 , and therefore we could conclude that it is the dominant contribution: it is however exactly the term which must be suppressed in order to avoid |m 3 | as the central eigenvalue. Furthermore, this term is equal to half of the solar mass squared difference that is about 30 times as small as the atmospheric splitting (see table 1). As a result, to recover a value for ∆m 2 atm close to the measured one, we need that the second term on the right-hand side of eq. (39) is the dominant one. We can estimate the ratio between the two terms, calculating r, the ratio of the two mass squared differences: The natural range of this quantity r −1 would be something like [0.3 − 0.7] (central value 0.5 and corrections of order λ). However, measurements give r = 0.032 +0.006 −0.005 ≈ λ 2 , or in other words 1/r ∼ 30 ∼ λ −2 . We conclude that The most natural explanation may be assuming that cos (θ k 0 − θ k 1 ) is very small. In that case the absolute values of all parameters can still be of order one, which was part of the naturalness requirement of the model. Neutrinos present a quasi degenerate (QD) spectrum and both normal and inverse ordering are possible. The typical scale of neutrino masses v L is given in this case by Possible alternative solutions of eq. (41) that give a non-QD spectrum can be obtained only by admitting the parameters belong to a larger, but less natural, range, λ 2 − λ −2 . In this case when k 0 is of order λ −2 and cos (θ k 0 − θ k 1 ) ∼ λ while k 1,2 still of order one a inverse hierarchical spectrum can be obtained. Another possibility to get an inverse hierarchical (IH) spectrum is having k 0 and cos (θ k 0 − θ k 1 ) of order one, k 1 ∼ λ −1 and k 2 ∼ λ 2 . Finally the normal hierarchical (NH) spectrum can be obtained only in the case in which k 2 ∼ λ −2 , cos (θ k 0 − θ k 1 ) and k 1 of order 1 and k 0 ∼ λ 2 .

Mixing matrices at Leading Order
Apart from the charged fermion masses of the second and third families, mixing angles and masses at leading order do not fit the experimental data: the charged fermion first families are massless, the CKM matrix is the unity matrix and the PMNS matrix can account for two maximal rotations in the (12) and (23) sectors: At this level of approximation, the PMNS matrix corresponds to the BM pattern discussed in the introduction. As already stated, this mixing scheme does not fit the data, due to the large value of the bimaximal solar angle which is out of about 5σ. Only considering the NLO contributions, which we will study in the next section, the model agrees with the measurements.

Fermion mass matrices at higher orders
The NLO contributions in the mass matrices originate from two sources: the first are the higher order terms in the superpotential, while the others come from the insertion of the NLO flavon VEVs in the operators in eqs. (22,33). In section 3.4 we will show how the flavons develop VEVs and how they are corrected. Here we anticipate the results, reporting the flavon VEVs in a form which is useful for the discussion in this section: Some comments are noteworthy: the subleading corrections are suppressed with respect the LO terms as δv/v ∼ λ and δ 2 v/v ∼ λ 2 for each flavon; NLO corrections to the second and third entries of ϕ and ϕ are present, but they present the same structure of the LO terms and can be re-absorbed; similarly, the NLO corrections to the third entry of χ and to σ are present, but they can be re-absorbed into the LO terms; the other entries of χ do not receive any corrections at NLO. If we consider the NNLO approximation level, i.e. corrections of relative order λ 2 with respect the LO terms, we see that the second and the third entries of ϕ ( ϕ ) are not (anti-)aligned anymore and that the second entry of χ is filled in. It is interesting that the first entry of χ is still vanishing at this level. We will see the relevance of this structure in a while.

Dirac mass terms
The Dirac matter superpotential at NLO is given by where, as for eq. (22), we adopt a compact notation: Note that not all of these terms are non-vanishing when the flavons develop VEV: in particular X (9) i for any i do not give a contribution when the LO VEVs are considered: only when the corrections to the VEVs are introduced, they contribute to the mass matrices. For this reason also the terms with X (10) i must be taken into account, even if they are suppressed by an additional Λ. The other terms which are vanishing at this order of approximation are X and m N LO where we used the definitions In the previous definitions we can see that eachỹ is the sum of two pieces: the first refers to the terms in eq. (45) when the LO flavon VEVs are considered; the second comes from the terms in eq. (22) where the NLO flavon VEVs are introduced. The only exception is F 3 which refers to the term proportional to X (9) i in eq. (45) and that give contribution only when the NLO flavon VEVs are considered. These two parts are of the same order of magnitude, since δv/v ∼ λ and δ 2 v/v ∼ λ 2 . Note that F i are distinct linear combinations of the arguments in the squared brackets. The expressions in eqs. (47,48,49) are valid at NLO level and note that the (anti-)alignment between the second and third entries of the (second) third families are still preserved. When considering higher order terms, this feature is lost and the (1, 1) entry of each mass matrix is filled in.
The values for the charged fermion masses given in eq. (28) are modified only by substituting the coefficients y i with their tilde-versions: At this approximation level, the first family masses are not yet well-described, because they are too small.

Majorana mass terms
Moving to the Majorana part of the matter superpotential at NLO, we get where as usual we used the compact notation with A few comments are in place. Note that all the terms proportional to k 3 can be reabsorbed by a redefinition ofk 0,1,2 and that the only new structure which corrects M L comes from the term F L F L ∆ L χ when we consider the subleading corrections of χ : as a result the entries (1,3) and (3,1) of M L are filled in by terms proportional to λ 3 . Regarding the contributions to the Majorana mass matrix for the RH neutrinos, it is important to say that the terms proportional to z 3,(i) with i = 1, . . . , 3 and to z 4,(i) with i = 1, . . . , 4 are vanishing, due to the particular flavon VEV alignment of the model. As a consequence all the NLO contributions to M R are of the order of λ 6 , apart that one to the entry (2,2) which is of the order of λ 4 . Finally, we note that each entry of M R is independent from all the others, being F c i singlets of the flavour symmetry, and therefore all the z i are free parameters with modulus of order 1. We listed only the dominant contributions, but the higher order terms would correspond to subleading corrections, which we can safely neglect in the following.
As a result of this analysis the Majorana masses for LH and RH neutrinos are given by where with the notation k i we account for all the redefinitions done on the parameters. As already stated when discussing the LO mass matrices, the contributions to the effective light neutrino mass matrix come from the type-I and type-II See-Saw mechanisms. The resulting NLO m (type-I) ν is given by and it is diagonalized by a maximal rotation in the (23) sector and not in the (12) sector as demanded by eq. (5). As a result, we need that the type-I See-Saw contribution is at least O(λ 2 ) the type-II See-Saw one, which we can express using eq. (51) as We remember here that v L , v R are the VEVs of ∆ L and ∆ R respectively. In particular v L is the VEV developed by the SM (1, 3, 1) triplet component of ∆ L and it is induced once the EW symmetry is broken. As we will see in detail in the next sections the physical SM triplet T ∼ (1, 3, 1) arises by the mixing between the SM (1, 3, 1) components of ∆ L and of two additional fields, Σ and Σ transforming under the PS gauge symmetry as (1,3,3). The VEV T of the SU (2) L triplet T is related to its mass M T through the following expression where α ij are numerical coefficients arising by the details of the scalar potential and v u 1 = h u , v u 2 = h u are the VEVs of the two up-type light Higgs doublets needed for the realization of our model (see appendix B for details). Since v L is the projection of T along ∆ L neglecting fine tuned cases it holds that v L ∼ T .
From eq. (56) we need therefore that M T and v R satisfy At the same time neutrino mass data imply that Combining the constraints of eqs. (59)-(60) we see that in the most natural scenario, Nevertheless if we allow the numerical factors α ij laying in the range 0.1 − 10, then M T , and consequently v R , can be reduced even of two and one orders of magnitude respectively In the following discussion of lepton mixing and in the phenomenological analysis, we will assume that indeed type-II See-Saw is dominating and we will neglect the type-I contributions. In the section devoted to the study of the scalar potential we will justify and find out the region of the parameters space where indeed type-II See-Saw dominates over type-I.

Mixing angles at the Next-to-Leading Order
Looking at eqs. (47)-(49)-(54) we see that the fermion mass matrices are of the required form as in eqs. (10)- (13). The resulting mixing matrices are modified with respect to the LO approximation and interesting new features follow. On the quark sector, the CKM matrix receives deviations from the unity and at NLO the angle θ q 12 is not vanishing anymore: Looking at this result, the meaning of the parameter λ is clear: it is defined as the ratio of the flavon VEVs over the cut-off of the theory, but it also determines the order of magnitude of the Cabibbo angle. This justify our initial assumption of λ = 0.2. We note that this result is in concordance with the discussion below equation (13). If the second columns of the up-and down-quark matrices are not proportional to each other, we can generate a non-vanishing Cabibbo angle θ q 12 , while the two other angles in the CKM matrix are still vanishing.
In the lepton sector, the PMNS matrix corresponds to the BM scheme which receives the corrections as illustrated in eq. (5) and as a result the solar and reactor angles are given by As it is easy to see, the reactor angle and the deviation from the maximal value of the solar angle are of order λ and therefore the model is now in agreement with the experimental data 5 . This corresponds to a weaker version of eq. (1). At this approximation level, the atmospheric angle remains maximal. As already stated, the forthcoming neutrino appearance experiments will provide a good test for the model, indeed they will be able to verify if θ l 13 is close to its present upper bound.

Higher order effects
At the next-to-next-to-leading order (NNLO) and even higher orders, many new terms appear in the superpotential. However, only few of them lead to new terms in the mass matrices, while the rest can be absorbed in redefinitions of the parameters as in eq. (50). For this reason we do not report the full list of NNLO contributions, but we just comment on the physical consequences. Three effects are worth mentioning: -as expected, the masses of the first families are strongly suppressed. We find that the down-quarks and the electron masses are suppressed by a factor of λ 6 and the up-quark mass by a factor of λ 8 . This leads to the following mass hierarchies: As it is easy to see, the electron and the up-quark masses perfectly fit the experimental values and there is only a small tension in the down-quark sector, where we expected a smaller suppression. However, the large number of parameters, in particular from φ, φ as well as ρ couplings, which enter in the definition of the first families allows to correctly fit the down mass, without affecting other observables.
-the (anti-)alignment of the (23) and (33) elements of the Dirac mass matrices in eqs.
(47-49) gets broken as the new terms appear. The new elements are λ 2 suppressed with respect to the older terms. As a result, the matrix that diagonalizes M i M † i , with i = e, u, d, has no longer an exact maximal mixing in the (23) sector. In the lepton sector, this translates to a λ 2 deviation from maximality in the atmospheric angle of the PMNS matrix. In the quark sector, the angle θ q 23 becomes of order λ 2 . It is interesting to note that θ q 13 remains vanishing at this order. It only appears when even stronger suppressed terms are taken into account and is of order λ 3 , in accordance with the Wolfenstein parametrization [31].
-the third columns of the mass matrices in eqs. (47) Therefore, also at NLO, eqs. (30) are fulfilled. At the NNLO level, terms proportional to v d ρ appear in the third columns of charged lepton and down-type quark matrices and terms proportional to v d φ in the second columns. The new terms are λ 2 ≈ 5% suppressed with respect to the old entries. We thus expect deviations from the relations |m τ | = |m b | and |m µ | = 3|m s | at the 5% level.

Flavon scalar potential
In this section we comment on the vacuum alignment mechanism which explains the flavon VEVs as in eqs. (17)-(18)- (19). It turns out that our desired alignment is exactly the one presented in [23], once we transform all the fields in our basis. Notice indeed that it is possible to identify each flavon of table 3 with the flavons in [23], by simply comparing the transformation properties under the full flavour group: Using the following unitary matrix to move from the basis in [23] to our basis, we find (up to irrelevant phases) the following flavon VEV alignment 6 . We will comment in a while about the presence of equivalent solutions.
which correspond to eqs. (17)- (18). In [23], a set of driving superfields has been introduced to guarantee the correct flavon vacuum alignment: these new fields are gauge singlets and transform only under the flavour group, but differently from the flavons, they do not develop vacuum expectation values 7 . Under the continuous R-symmetry all the driving fields transform as U (1) R = 2 and therefore, constructing the superpotential, they appear only linearly. For the same purpose, we introduce in our model a set of driving fields which recall those in [23]. In table 5 we show the driving fields and their transformation properties under S 4 × Z 4 .
It is easy to determine the correspondence between our set of driving fields and those of [23]: We now construct the driving superpotential w d , which contains only flavons and driving fields and in particular neither matter fields nor Higgses, and look for the conditions that minimize the scalar potential, where Φ i denote collectively all the scalar fields of the theory, m 2 i are soft masses and dots stand for D-terms for the fields charged under gauge group and possible additional soft breaking terms. Since m 2 i are expected to be much smaller that the mass scales involved in w d , it is reasonable to minimize V in the supersymmetric limit and to account for soft breaking effects subsequently.
Since our flavon and driving field content exactly corresponds to the one in [23], we already know that the VEV alignment in eqs. (17)-(18) represents an isolate minimum of the scalar potential. We only need to identify the relations which link the VEVs v ϕ , v ϕ , v χ and v σ among each other in our model. To this purpose we write the driving superpotential, which reads as where the first line deals with only the fields ϕ and ϕ and the other two with χ and σ. As also explained in [23], this leads to the alignment in eqs. (17,18) where the VEVs satisfy to The solution in eqs. (17,18,73) is not unique, but it is possible to introduce a set of soft supersymmetric breaking parameters, which selects this solution as the lowest minimum of the scalar potential.
It is interesting to note the presence of an other source of uncertainty in our solution, which minimizes V . Given the symmetry of w d and the field configurations of eqs. (17,18,73), by acting on them with elements of the flavour symmetry group S 4 × Z 4 , we can generate other minima of the scalar potential. These alternative solutions however are physically equivalent to those of the original set and it is not restrictive to analyze the model by choosing as local minimum that one in eqs. (17,18,73).
For the FN field θ to acquire a VEV, we assume that the symmetry U (1) F N is gauged such that θ gets its VEV through a D-term. The corresponding potential is of the form: where g F N is the gauge coupling constant of U (1) F N and M 2 F I denotes the contribution of the Fayet-Iliopoulos (FI) term. Dots in eq. (74) represent e.g. terms involving the fields F c 1 and F c 2 which are charged under U (1) F N . These terms are however not relevant to calculate the VEV of the FN field and we omit them in the present discussion. V D,F N leads in the supersymmetric limit to It is relevant to underline that the VEVs in eqs. (73, 75), depend on Mass parameters: all these mass scales naturally have the same order of magnitude and as a result V EV s/Λ f ∼ λ. The only exceptions are the VEVs of ϕ and ϕ , which depend on a flat direction. In the model, we simply assume that their VEVs have values of the same order of all the other flavon VEVs.

Higher order contributions
In this section we briefly comment on the corrections which enter in the flavon VEVs, once the higher order contributions are taken into account. We leave all the details to the appendix C.
In the superpotential w d , the flavons which contribute to the Dirac mass terms, ϕ and ϕ , and those which contribute to the Majorana mass terms, χ and σ, at LO belong to two separated sectors, indeed any mixing term is prevented due to the Z 4 symmetry. This situation is not preserved at NLO, since the fields χ and σ are neutral under the Z 4 symmetry and therefore we can add each of them to all the terms in w d . This leads to modifications to the LO VEV alignment of eq. (17) and it turns out that the first entries of ϕ and ϕ are filled in, while the second and third entries are corrected by terms which can be however absorbed into the LO ones, without spoiling the alignment. Also the VEVs in eq. (18) receive some corrections: the first and second entries of χ still vanish and the NLO contributions to the third entry can again be absorbed into the LO term. This discussion justifies the results showed in eq. (44).
Corrections from the NNLO contributions are without particular alignments: in particular, the second and third entries of ϕ and ϕ are no longer related, and also the second entry of χ gets a non-zero value. It is interesting to note the the first entry of χ remains zero.

Higgs scalar potential
In this section we present the study of the Higgs potential in our model. It is an interesting example of how the introduction of flavour symmetries and the assumptions done to get the correct mass matrices have non-negligible consequences on the Higgs sector. As a result, the study of the Higgs scalar potential and of the gauge and the Yukawa coupling runnings in a general non-flavour PS context does not strictly hold. Notice that even if the following analysis refers to our particular choice of fields and symmetries, our conclusions can be taken as a general hint for a very large class of models that combine a discrete flavour symmetry with a grand unified scenario: indeed our model building strategy shares common features with other constructions. In particular the Higgs fields usually transform under the flavour symmetry G f and this has direct consequences on the implementation of the grand unified symmetry breaking. Moreover type-II See-Saw dominance and particular patterns of vanishing projections of the heavy Higgs fields on the light Higgs doublets are frequently required to get the correct fermion mass matrices.
In table 6 we list all the Higgs fields which are necessary to reproduce the correct mass matrices and to implement the desired PS symmetry breaking pattern: in the first part of the table we report the Higgses already introduced in table 4, while the second part contains all the additional fields.

PS
(1, 2, 2) (15, 2, 2) (10, 3, 1) (10, 1, 3)  Some of these Higgs fields are already present in the minimal version of PS models [34]: typically a (15,1,1) multiplet -as the fields A and B in table 6 -is used to break SU (4) C down to SU (3) C × U (1) B−L and to induce the VEVs of the couple of (10, 1, 3) ⊕ (10, 1, 3)corresponding to the fields ∆ R ⊕ ∆ R in table 6. The latter VEVs breaks SU (2) R × U (1) B−L into the SM hypercharge U (1) Y concluding the symmetry breaking chain from the PS gauge group to the SM SU (3) C × SU (2) L × U (1) Y . The field that triggers the EW symmetry breaking is usually a bidoublet (1,2,2) -as φ or φ in table 6. The fields (10, 3, 1)⊕(10, 3, 1) -the fields ∆ L ⊕ ∆ L in table 6 -do not develop VEVs at tree level in the usual minimal PS model, but only when next to leading order terms are taken into account; these are typically suppressed by the Planck scale, as already stated in [34]. For this reason in the minimal PS the type-II See-Saw contributions to the effective neutrino masses are almost negligible.
We identify three main reasons for which the existent studies of the symmetry breaking patterns in the PS context [24,34] have to be modified and this automatically justifies the presence of the new fields in table 6: -the assumption v u ρ = 0 necessary to distinguish the up-quark sector from the others and to recover the up-quark mass hierarchies can be realized only if we include two identical copies of bidoublet (1,2,2), φ and φ (see details in appendix B), and we then impose that four SU (2) L doublets (2 up-type and 2 down-type) remain light; -since the fields ρ, ∆ R and ∆ R transform non-trivially under the flavour symmetry Z 4 , it is necessary to introduce two copies of (15,1,1) multiplets, A and B, with opposite Z 4 charges, 1 and -1 respectively: A is responsible of inducing the breaking of SU (2) R through its coupling with ∆ R and ∆ R ; B allows the coupling of the bidoublets φ, φ with the bidoublet ρ. In this way all of these three fields have a non-vanishing projection on the light Higgs SU (2) L doublets; -the component of ∆ L which corresponds to the usual SM triplet (1,3,1) can develop a VEV once the EW symmetry is broken only in the presence of a trilinear coupling with the SU (2) L Higgs doublets. This coupling cannot be originated with only the Higgs fields ∆ L ⊕ ∆ L and the field content given in table 6. For this reason we need an additional field which mediate this coupling and the simplest choice would be a bitriplet Σ ∼ (1, 3, 3) that can couple with the fields φ, φ or ρ and at the same time can mixes with ∆ L , when ∆ R develops VEV at the SU (2) R breaking scale. However, once more, the presence of the Z 4 symmetry obliges the introduction of two distinct (1,3,3) Higgs fields, Σ and Σ , with opposite Z 4 charges, 1 and −1 respectively. In this way Σ can couple to the bidoublets φ, φ or ρ, while Σ can mix with ∆ L . Finally, we need a new ingredient that mixes Σ with Σ : a PS singlet ξ charged −1 under Z 4 can do the job.
The scalar part of the superpotential is then given by 8 The vacuum configuration at the GUT breaking scale is given by 9 8 Since all the Higgs fields are neutral under the continuous U (1) R , the scalar superpotential explicitly breaks it, while preserving the usual R-parity. The terms in eq. (76) could be generated from a U (1) Rconserving superpotential in which the breaking is mediated by additional fields, which are U (1) R = 2 and develop non-vanishing VEVs. For instance the mass term M φ φφ could originate from a trilinear term X φφ, when X = M φ . Similarly the trilinear coupling λ ξ ΣΣ ξ could originate by the non-renormalizable term X ΣΣ ξ/Λ , when X /Λ = λ ξ and Λ is the energy scale of the dynamics of the field X. In our model we simply assume the existence of the terms in eq. (76) in the superpotential and allow for an explicit breaking of the U (1) R symmetry in this sector. 9 We redefine ∆ R = v R as ∆ R = M R in order to adapt to the usual notation.
The VEVs A and B break the SU (4) C to SU (3) C × U (1) B−L , while ∆ R and ∆ R break SU (4) C × SU (2) R into SU (3) C × U (1) Y . Therefore given our field content the colour breaking scale, M C = Max(M C 1 , M C 2 ), is never smaller than the SU (2) R breaking scale, M R . We may expect that M C 1 ∼ M C 2 ≡ M C but in principle they could be different. Finally V ξ is expected to be close to the flavour breaking scale, due to the ξ gauge singlet nature. The PS breaking pattern so far sketched is therefore summarized by The F-derivative system obtained by the superpotential in eq. (76) reads as By solving the previous equations, we can express the mass parameters that enter in the superpotential in term of the adimensional parameters λ i and the physical breaking scales. All the details regarding the mass spectrum are reported in the appendix B, but some comments are in place. As in the minimal supersymmetric PS [34] when the singlet component of A develops a VEV, there is an accidental SU (3) symmetry involving ∆ R and ∆ R . When the singlet components of these fields acquire a VEV the accidental symmetry is broken to SU (2) giving rise to 5 Goldstone Bosons (GBs). At the same time SU (2) R × U (1) B−L is broken down to U (1) Y , eating 3 of the 5 GBs. Therefore 2 of them, namely δ ++ and δ ++ are left massless, down to the SUSY soft breaking scale ∼ 1 TeV. This is a well known prediction of SUSY PS theories, which can be tested at LHC [35]. On the other hand, contrary to the minimal case described in [34], due to the mixing between the Higgs fields A and B, no colour octet is lighter than M R . In order to assure type-II dominance and to get the correct PS symmetry breaking pattern, it is necessary that M T ≤ M R ≤ M C ≤ V ξ . We will see in a while the constraints on the parameters which can be derived from this special mass ordering. For the moment, we just assume this scheme and we describe the mass spectrum in terms of the SM gauge group. Starting from the heaviest energy scale With the Higgs field content given in tables 4 and 6 and the scalar spectrum so far sketched we can plug the gauge coupling runnings and see if the conditions given in eq. (61) or those in eq. (62) can be satisfied. Furthermore, in the study of the gauge coupling runnings (see appendix D for details) from the M GU T to the EW scale we have to impose the following constraints: -recover the EW values for α 3 , α 2 and α 1 , related to the gauge couplings of the SM gauge group ; -define the GUT scale as the scale at which the largest α i = 1. In this way we are sure to be in a perturbative regime up to the GUT scale and thus we are allowed to adopt the one-loop renormalization group equations (RGEs).
Even using all these constraints we are left with two more freedoms, the value of the SU (4) C and SU (2) R breaking scales, i.e. M C and M R respectively. We adopt two distinct approaches, that we will indicate as the more constraining and the less constraining ones. In the first case we define M C the scale at which the largest α i is equal or smaller than 1/4π. In this way all the gauge coupling at M C are smaller than 1. In the second case we allow the largest α i to correspond to a gauge coupling in the range 1 < g i < 3. Then M R should satisfies eq. (61) or eq. (62), but its exact value is not fixed yet.
A few general comments are needed. The non-minimal PS field content affects the gauge coupling runnings in a non-negligible way. In particular the presence of the charged singlets δ ++ and δ ++ down to M SU SY deeply modifies the U (1) Y and SU (2) R gauge coupling evolution. It turns out that the largest α i above the M R scale is always α R . Therefore the two approaches we described can be formulated as follows: -More constraining approach ⇐⇒ α R ≤ 1/4π, at M C -Less constraining approach ⇐⇒ 1 < g R < 3, at M C .

More constraining approach
In this case there are no solutions neither for the ranges of values of M T and M R given in eq. (61) nor for those in eq. (62). Indeed we find that M R ≤ 10 12 GeV, as can be seen in fig. 1. In other words if we adopt this constraining approach to fix the value of M C , the type-I and type-II See-Saw scales require Yukawa parameters which are at least 2 order of magnitude far from their natural values to reproduce the correct neutrino mass scale. We cannot be satisfied with such a solution because in this case the type-II See-Saw dominance is obtained by increasing the Yukawa couplings in the right neutrino sector and at the same time reducing the coupling of the left-handed neutrinos with the scalar triplet. Even if this may be considered a solution, our challenge was to provide a justification of type-II See-Saw dominance through the analysis of the Higgs scalar potential and not by tuning the Yukawa parameters. Moreover we introduced a FN Abelian symmetry to explain the small (≤ 10 −2 ) Yukawa parameters necessary in the charged fermion sector to reproduce the correct mass hierarchies. The presence of Yukawa parameters of this order in the purely left-handed neutrino sector makes the introduction of the FN symmetry questionable. Notice that in the discussion of the neutrino spectrum we have enlarged the parameter range up to values not smaller than λ 2 .

Less constraining approach
In this case there are solutions only for the second range of values of M T and M R given in eq. (62). As can be seen in fig. 2, M R can now reach the value of 10 13 GeV. However the three scales M R , M C and M GU T are compressed in a narrow region around 10 13 GeV and therefore our model is described by an extended MSSM model almost up to the GUT scale 10 . Nevertheless, its PS origin is reflected in the non-trivial relations between the Yukawa couplings. In conclusion, by admitting the Yukawa parameters span in a range 0.1-10 and switching M R , M C and M GU T very close to each other, we have found a narrow 10 Above this energy scale an other gauge structure could be active. We do not take in consideration an high-energy completion of the model, but it is reasonable that larger gauge groups or particular constructions could be present at these energies: for example an SO(10) inspired approach in which fermions do not belong to a unique representation. region of the parameter space where our model could still give a realistic description of fermion masses and mixings and in which type-II See-Saw dominance is not imposed by hand.
To finally consider our model viable, we should study the stability of the flavour structure of the mass matrices under the RGEs from the GUT scale down to the EW one. The study of the full set of the RGEs of the model presented is beyond the purpose of this paper. For this reason we will neither run the parameters of the scalar superpotential nor include and run the parameters of the soft SUSY breaking potential. Under these approximations the EW vacuum expectation values do not change from the GUT scale down to the EW one. However this does not affect our conclusions for what concerns the stability of the mass matrix structures since the EW VEV shifts due to the running factorize out and leave the Yukawa flavour structure unchanged. The study of the Yukawa matrix running will be done in the next section.

Yukawa coupling running
In the previous section we analyzed the constraints on the scalar Higgs sector coming from the requirement of type-II See-Saw dominance and from the presence of the flavour symmetry under which the Higgs fields non-trivially transform. We found that the model is viable only in a small region of the parameter space for which M R , M C and M GU T ∼ 10 13 GeV are very close to each other. At the same time M T lies at only one order of magnitude below M R . For these reasons to study the stability of the flavour structure of the fermion mass matrices at low scale we can consider only the running from M T onwards, thus neglecting the running from higher energies. Furthermore we also neglect the running from M SU SY to the EW scale, which would introduce only minor corrections. We work under the assumption that type-II See-Saw is dominating over type-I and moreover that the effects from the type-I terms under the RGEs are negligible. Therefore in the studying of the Yukawa coupling running we do not take into account the Weinberg operator originating by integrating out the right-handed neutrinos. The error we introduce in this way is less than λ 2 and we will see in a while that these contributions do not modify our results. Furthermore, we study the stability under the Renormalization Group (RG) running, in the approximation corresponding to the NLO, i.e. considering the mass matrices introduced in eqs. (47)-(49).

Yukawa matrices at M T
Since we start the renormalization group running at M T , at this scale we integrate out the SU (2) L scalar triplet T obtaining an effective Weinberg operator responsible of the type-II See-Saw contribution. We recall here the origin of this effective operator. The Majorana parts of the matter superpotential given in eqs. (33), (52) contain terms with the coupling while the scalar part of the superpotential in eq. (76) contains the terms that ensures the mixing between the (1,3,1) ((1,3,-1)) components of ∆ L (∆ L ), Σ and Σ , whose lighter combination is identified with T (T ), and provides the coupling of T (T ) with the light doublets h d and h d (h u and h u ). The effective Weinberg operator at M T is given by where L i is the SU (2) L lepton doublets, h u 1 = h u , h u 2 = h u and α ij are coefficients arising by the scalar potential and Y L is given by Notice that we are neglecting the λ 3 terms in Y L , because they would be irrelevant for the following analysis.
For what concerns the charged fermion Yukawa, at M T the Dirac part of the superpotential is written as The U matrix defines the light SU (2) L Higgses in term of the PS Higgs field components, as explicitly written in appendix B, andỹ i has to be read as U 11ỹ (2) i . With respect to the mass matrices given in eqs. (47-49) we have reabsorbed a power of λ.

Analytical approximations
In the appendix E, we report all the RGEs for the Yukawa matrices, while here we discuss the results. The RGEs present the general compact expressions where the index f runs over {e, u, d}, the parameter t is defined as When we fix µ 0 ∼ M T and µ ∼ M SU SY these formulas can be approximated by For what concerns the quark sector we find the following approximated expressions for the masses of the last two families and the Cabibbo angle where γ = m 2 t /(v u 1 + βv u 2 ) 2 and the masses and the Cabibbo angle on the right of the previous expressions are intended at the M T scale. Note that the demand that m 2 b is still positive at M SU SY gives an upper bound on γ of 0.7. The charged lepton masses are very similar to the down-quark masses and indeed we have We now consider the neutrino sector (see [36] for a general approach at RGEs with or without flavour symmetries) and the modification due to the RG running. We remember that our model at the GUT scale naturally predicts the QD spectrum, with both normal and inverse ordering, while choosing a less natural range of the parameters space we may have both NH and IH spectrum. To analyze the effect of the RGEs on the neutrino mass matrix M ν , it is worth rotate M ν of a maximal rotation in the (12) sector. Then, at the high scale, M ν is diagonal and reads as where without loss of generality k 0 can be taken real, by a redefinition on the phases. After the running at M SU SY eq. (91) gets a correction ∆M ν that is given by with w ± = (±γ/ √ 2 + √ 2λ), for the QD and IH case when k 2 ∼ O (1). For the NH case characterized by k 2 ∼ λ −2 , ∆M ν assumes the following form with k ± = (k 1 λ ± k 2 λ 2 ). We can now consider the three different cases QD, NH and IH, which the model accounts for at the GUT scale.
The correction given by the running induces a rotation in the (23) sector characterized by being ∆t ∼ 2λ 2 . This is a large contribution, which deviates the atmospheric angle from the initial maximal value, spoiling the agreement with the experimental data at 3σ. A possible way-out would be if this large correction is erased by a corresponding large correction in the charged lepton mass matrix. However this is not the case, because the maximal θ e 23 in the charged lepton mixing matrix is stable under the RG running. As a result, the QD case is not viable.
The corrections both for the atmospheric and the reactor angles are of order ∆t γ/2 √ 2 ∼ 3λ 3 and can be safely neglected. Analogously, also the mass splittings receive deviations which can be neglected. On the other hand, the charged lepton mixing matrix is stable under the RG running. As a result the three mixing angles at M SU SY can be well approximated with their initial values at M T .
All the corrections to the neutrino mixing are of order γ∆t / √ 2. While the solar mass splitting receives negligible contributions, the atmospheric mass splitting is corrected as follows Combining the neutrino mixing with the charged lepton one we get the following mixing angles for the lepton sector: where the corrections coming from the RG running are proportional to ∆t . In particular the corrections to θ l 12 (M SU SY ) and θ l 13 (M SU SY ) come from the charged lepton sector, while that one to θ l 23 (M SU SY ) arises only by the neutrino sector. This distinction strictly holds only at this order of approximation, while considering contributions of the order of λ 3 all the angles receive corrections from both the sectors.
With respect to the previous case, the solar mass splitting gets a non-negligible correction, that can be written as For what concerns the atmospheric mass splitting and the lepton mixing angles we recover the same shifts as in the previous case.

Neutrino Phenomenological Analysis
In the previous section, it was concluded that the quasi-degenerate spectrum is unstable under the RG running and becomes phenomenologically inviable due to too large corrections to the atmospheric mixing angle. Only the normal and inverse hierarchies were shown to be stable and phenomenologically viable. This is different from [23], where only the normal hierarchy and the quasi-degenerate spectrum with normal ordering were found.
In this section, we will discuss neutrino phenomenology in more detail, focussing on the value of the reactor mixing angle θ l 13 and the possibility of neutrinoless double beta decay. The analytical expression of the reactor mixing angle is given by eq. (64) at NLO and it is corrected by the RG running as in eq. (96) for the two IH scenarios studied in the previous section. We see that θ l 13 is typically of order λ ≈ 0.2, so sin 2 θ l 13 ≈ 0.04, which is rather large, but still allowed at the 3σ level as shown in table 1. Here, we complete the study on θ l 13 , performing a numerical analysis. As can be seen from eqs. (47)-(54), the neutrino and the charged lepton mass matrices at NLO are a function of many parameters. However the GUT nature of the model allows us to fix most parameters that occur at LO because they enter in the low energy expressions for quark and charged lepton masses, as can been seen comparing eq. (28) with eq. (89). Note that all dimensionless parameters, i.e. theỹ i , can be fixed of order 1.
The other free parameters can be fixed as random numbers of order 1, except for the cases where we have argued in the previous section that they should have slightly larger or smaller values. Because of the use of random numbers, the predictions of our model are no longer single valued. The two vertical lines are the 3σ bounds for sin 2 θ l 12 according to [1]. The upper horizontal line is the 3σ upper bound for sin 2 θ l 13 and the middle line is the best fit value. The lower line is the 95 % exclusion confidence level after 3 years of Daya Bay data taking [5].
We plot the reactor angle versus the solar angle in figure 3. At NLO, eq.(65), θ l 12 is driven away from the maximal value π/4 by a term proportional to λ (note that we take only the corrections which decrease the value of the solar angle, neglecting those which increase it). We see that this deviation is not for all values of the parameters large enough to bring it in the observed region, although this happens for a significant number of them. As explained above, larger values of sin 2 θ l 13 are favoured and almost all points are in the sensitive region for experiments.
To study neutrinoless double beta decay, we consider the effective 0νββ parameter m ee , defined as In figure 4 we plot m ee against the lightest neutrino mass, which is m 1 and m 3 in the NH and IH case respectively. The future experiments are expected to reach good sensitivities: 90 meV [38] (GERDA), 20 meV [39] (Majorana), 50 meV [40] (SuperNEMO), 15 meV [41] (CUORE) and 24 meV [42] (EXO). As a result, looking at figure 4, the whole IH band will be tested in the next future and with it the two cases of our model which allow for the IH spectrum.

Conclusion
In this paper we have addressed several aspects of the interplay between a GUT based model and a discrete flavour symmetry. The paper should indeed be considered as the combination of two distinct parts: in the first one we mainly discussed the building of the model from the flavour point of view, while in the second one we faced the problem to justify the assumptions made in the first part and to achieve the correct gauge symmetry breaking chain. More in detail, the symmetry group of our model is P S × G f , where P S stands for the GUT Pati-Salam gauge group SU (4) C × SU (2) L × SU (2) R and G f for the flavour group Within this GUT context one has the relationship between the down-quark and charged leptons mass matrices, M d ∼ M e , which can easily be used to revise the old idea of quark-lepton complementarity. In the model this is obtained by the use of the non-Abelian discrete flavour symmetry S 4 properly broken through the VEVs of a set of flavon fields, which transform as triplets under S 4 . The additional Abelian symmetries, which enters in G f , play different roles: Z 4 keeps quarks separated from leptons and neutrinos from charged leptons and prevents dangerous couplings in the superpotential of the model; U (1) F N helps to justify the charged fermion mass hierarchies; U (1) R is a common ingredient of supersymmetric flavour models. It contains the discrete R-parity and is useful to build a suitable flavon superpotential that allows the correct S 4 breaking pattern.
Already at the leading order, the model shows nice features: we are able to reproduce the mass hierarchy between the third and the second charged fermion families, the bottomtau unification, the Georgi-Jarlskog [30] relation |m µ | = 3|m s | and, under the assumption of type-II See-Saw dominance at the GUT scale, a realistic neutrino spectrum. However at this level of approximation, both the CKM and the PMNS mixing matrices are not correct: the quark mixing matrix coincides with the identity matrix, while the lepton one is given by the BM pattern. It is worth to recall here that the BM mixing corresponds to maximal solar and atmospheric angles and to a vanishing reactor angle: only the solar angle is not in agreement with the data as it deviates from the experimental central value by a quantity close to the Cabibbo angle, λ ∼ 0.2.
At next-to-leading order, the wrong predictions for the fermion mixing angles are corrected: in the CKM matrix, the mixing angle θ q 12 receives contributions of the order of λ, fitting the value of the Cabibbo angle; analogously, in the PMNS matrix, the solar angle is corrected by the same amount and we find the nice result that θ l 12 ∼ π/4 − O(λ). At the same time, also the reactor angle receives significant contributions and indeed at this level of approximation it results θ l 13 ∼ O(λ): this is an interesting feature of our model, because this value is close to its present upper bound and it will be tested in the forthcoming neutrino experiments [4][5][6][7][8][9].
Once we consider the higher order terms, we find the other two CKM angles of the correct order of magnitude, θ q 23 ∼ O(λ 2 ) and θ q 13 ∼ O(λ 3 ), and small corrections are introduced in the PMNS angles: in particular the atmospheric angle becomes θ l 23 ∼ π/4 + O(λ 2 ), justifying a small deviation from the maximality. For what concerns the masses, all the fermions are massive and the mass hierarchies fit the experimental observations.
On the other hand, the neutrino spectrum could be either quasi degenerate or normal or inverse hierarchical. Only the first case corresponds to a completely natural choice of the parameters, which, in the absence of an explanation coming form a higher energy theory, should be of order 1: in order to allow the NH and the IH, the parameters should span in a larger range of values, namely λ −2 ÷ λ 2 .
In the second part of the paper we have studied the Higgs scalar potential and the running of both the gauge couplings and the Yukawa mass matrices under the RGEs. With this analysis we looked for the constraints which arise to justify the Higgs field VEV pattern used in the flavour section and the assumption of the type-II See-Saw dominance. The presence of the flavour group G f modifies the Higgs field content necessary to implement the classical breaking pattern of the PS gauge group and as consequence not all the results obtained by studying minimal versions of PS are recovered. In particular we need the presence of two PS multiplets (15,1,1), A and B, responsible to break the unified colour symmetry SU (4) C to SU (3) C × U (1) B−L . Two copies of the (1,2,2) multiplet, φ and φ , and one (15,2,2) field, ρ, are necessary to implement the condition v u ρ = 0 . Lastly, we need the new fields Σ, Σ ∼ (1,3,3) and ξ ∼ (1, 1, 1) to have a type-II See-Saw contribution at tree level. Gauge coupling runnings are affected by the large field content and in particular we found that the requirement of having type-II See-Saw dominance constrains the model in a small region of the parameter space, in which all the heavy mass scales are sandwiched between 10 11 GeV and 10 13 GeV. At the same time Yukawa mass matrix running shows that while the CKM Cabibbo angle is stable under the RGEs evolution, the PMNS mixing angles are stable only if neutrinos present a NH or an IH spectrum, ruling out the QD case. As already stated, the QD spectrum would be the most natural and probable case at the GUT scale, but the Yukawa RGEs analysis further reduces the allowed region of the parameter space.
In section 6, we performed a brief phenomenological analysis on neutrino observables considering all the constraints which come from the flavour and the Higgs sectors. We have first considered the value of the reactor angle in terms of the deviations of the solar angle from the maximal value and the numerical analysis confirmed the analytical results: in our model, θ l 13 naturally acquires a value not far from its present upper bound. After this, we studied the neutrinoless double beta decay effective mass, m ee , and we have seen that the next future experiments are expected to reach sufficiently good sensitivities to test our model in the IH regime.
As a final comment, we underline that there is a strong tension in combining a GUT model with a (discrete non-Abelian) flavour symmetry and therefore a parallel study is not only interesting but also recommended to provide a viable model. In the first part of the paper we produced a flavour-GUT model that is quite realistic and viable. Only the study done in the second part, regarding the Higgs sector reveals that the model is restricted into a small region of the parameter space, reducing the freedom in the choice of the parameter values. Even if our results are model dependent, our construction shares many features with other models present in literature, where often a detailed discussion of the Higgs sector is missing. In our opinion, this is a serious drawback and we would suggest to consider the interplay between GUTs and flavour symmetries in this kind of models as well.

Acknowledgments
We thank Guido Altarelli and Ferruccio Feruglio for useful comments. The work of RdAT and FB is part of the research program of the Foundation for Fundamental Research of Matter (FOM). The work of FB has also been partially supported by the National Organization for Scientific Research (NWO). The work of LM has been partly supported by the European Commission under contracts MRTN-CT-2006-035505 and by the European Program "Unification in the LHC Era", contract PITN-GA-2009-237920 (UNILHC).
A The symmetric group S 4 In this appendix we report the character table and the Clebsch-Gordan coefficients of the S 4 discrete group in our basis. The generators, S and T , obey to the following rules and are of the following form for the five different representations: Using these generators we calculate the Clebsch Gordan coefficient for all the Kronecker products. In the following we use α i to indicate the elements of the first representation of the product and β i to indicate those of the second representation. We start with all the multiplication rules which include the 1-dimensional representations: The multiplication rules with the 2-dimensional representation are the following: The multiplication rules with the 3-dimensional representations are the following:

B Higgs scalar spectrum
We now proceed to presenting the scalar mass matrices for all the fields introduced in sec. 4 according to the SM∼ SU (3) C × SU (2) L × U (1) Y representations, indicating their origin with respect to the PS and the colour-broken Pati-Salam (CbPS) phase. We recall that in the colour broken phase the symmetry group is given by SU (3) For completeness for each field we also indicate the corresponding T 3R value, with T 3R the diagonal generator of SU (2) R . We use the same notation as in [37], in which we write the Dirac scalar mass matrices as they could be read directly from the superpotential of eq. (76) at the scale M R . We label the mass matrices with S, D and T according to the singlet, doublet or triplet representations of the SU (2) L gauge symmetry.
M S1 has a vanishing eigenvalue that is eaten by the corresponding gauge boson. Moreover it can be checked that one of the singlets, which we call ξ 0 , has a mass ∼ M R while all the others masses appear as combination of M C1 , M C2 and V ξ .   These states correspond to the last two massless GBs eaten by the respective gauge bosons. Indeed with these last two massless states the total amount of GBs eaten is 9 = 1+2×3+2, as it should be for the breaking SU (4) When M 2 D1 is diagonalized according to we get three up-type (down-type) Higgs doublets h u , h u , H u (h d , h d , H d ).
Therefore for the up (down) projections we have where v u,d In general a light doublet, i.e. massless at the M C scale, gets a mass term at M SU SY and its VEV, at the EW scale, is of order of the EW scale ∼ v W . On the contrary for a heavy doublet of mass M , its induced VEV at the EW scale is ∼ v 2 W /M , that for M ∼ M C is completely negligible with respect to v W . Consider now the condition ρ u = 0 that we imposed to get the correct fermion mass matrices. In the standard case we would have only one up-and one down-type light Higgs doublets, being all the other doublets heavy. Assume now that h u,d in eq. (106) are the upand down-type light doublets at M C . Then we have v u,d From eq. (107) we see that in this case it would be impossible to make vanishing the ρ projection along the up-direction (that implies U 13 = 0) still maintaining a non-vanishing ρ d . 11 For this reason the condition ρ u = 0 implies a non-standard scenario and the presence of two light doublets of up-type at the M C scale, namely h u and h u . However the symmetric nature of M D1 ensures that as consequence we are also left with two down-type light Higgs doublets, h d and h d . Nevertheless this does not imply that ρ d vanishes because we have since v u,d i depend on the soft terms. In conclusion we have to impose two constraints on the free parameters that enter in M D1 corresponding to require that M D1 has two vanishing eigenvalues. Notice that the condition of having M D 1 of rank 1 is fine-tuned but it is not more fine-tuned than imposing M D 1 of rank 2, which is universally accepted whenever the MSSM has to be recovered. In our case the fermion mass matrix structures impose a slightly different condition -M D 1 of rank 1 -but both the requirements are satisfied by fine-tuning the parameters that enter in the mass matrix.

CbPS
T 3R C 1 φ (1, 2, 2) (1, 2, 2, 0) 1/2   The superpotential w d , linear in the driving fields D R , ϕ R , χ R and σ R , is modified into: where δw d contains the NLO contributions, suppressed by one power of 1/Λ with respect to w d . The corrective term δw d is given by the most general quartic, S 4 × Z 4 -invariant polynomial linear in the driving fields, and can be obtained by inserting an additional flavon field in all the LO terms. The Z 4 -charges prevent any addition of the flavons ϕ and ϕ at NLO, while a factor of σ or χ can be added to all the LO terms. The full expression of δw d is the following: where x i , w i , s i and v i are coefficients and I σ R i , I χ R i , I D R i , I ϕ R i represent a basis of independent quartic invariants: I χ R 1 = (χ R χ)(χχ) I χ R 4 = (χ R (χχ) 3 1 ) σ I χ R 2 = ((χ R χ) 2 (χχ) 2 ) I χ R 5 = (χ R χ)σσ I χ R 3 = ((χ R χ) 3 1 (χχ) 3 1 ) (114) In the previous terms we indicate with (. . .) the singlet 1 1 , with (. . .) the singlet 1 2 and with (. . .) R (R = 2, 3 1 , 3 2 ) the representation R.
The NLO flavon VEVs are obtained by imposing the vanishing of the first derivative of w d + δw d with respect to the driving fields σ R , χ R , D R and ϕ R . We look for a solution that perturbs eqs. (17) and (18) to first order in the 1/Λ expansion: for all components of the flavons Φ = (σ, χ, ϕ, ϕ ), we denote the shifted VEV's by where Φ LO are given by eqs. (17) and (18). It is straightforward to verify the following results. In the Majorana mass sector the shifts δσ, δχ turn out to be proportional to the LO VEV's Φ LO and can be absorbed in a redefinition of the parameters v χ and v σ . Instead, in the Dirac mass sector, the shifts δϕ, δϕ have a non-trivial structure, so that the LO texture is modified: where v ϕ and v ϕ satisfy a relation similar to that in eq. (73) and the shifts δv ϕ and δv ϕ are suppressed by a factor λ with respect to the LO entries v ϕ and v ϕ , respectively.

D Gauge coupling running
In this appendix, we provide the coefficients of the β-functions for the gauge coupling running in the different regimes. The complete matter fields run from the GUT scale down to the M SU SY scale, where the SUSY partners decouple. For what concerns the scalar fields in section 4, we have already outlined the scalar spectrum according to the different scale at which the fields decouple. As a result the computation of the β-functions is straightforward. Calling µ the generic scale, we have -for M C < µ < M GU T : all matter is in the left and right handed multiplets (4,2,1) and (4, 1 , 2) as mentioned in table 2. In the Higgs sector, we have all the fields mentioned in table 6. This leads to the coefficients β SU (4) C = 54 , β SU (2) L = 69 , β SU (2) R = 69 .
Due to the large matter content these coefficients are very large and the β-functions are very steep. As consequence the theory is in the Pati-Salam regime only for a very small range of energies, as can indeed be seen in figure 2. Almost after passing the scale M C , the SU (2) R coupling constant enters the non-perturbative regime.
-M R < µ < M C : the theory undergoes to the SU (3) C × SU (2) L × SU (2) R × U (1) B−L symmetry. Considering the matter (left and right handed doublets of quarks and leptons characterized by different U (1) B−L charges) and the scalar fields, we find the coefficients -M T < µ < M R : in this regime we have all the usual MSSM matter particles, four light higgs doublets (two up-type and two down-type), a couple of SM triplets (1, 3, 1) ⊕ (1, 3, −1) and two extra charged singlets (δ ++ andδ ++ ). The coefficients of the β-functions are The hypercharge that appears in the last term is related to SU (2) R and the B − L charges in the previous regime by -M SU SY < µ < M T : in this regime we have all the usual MSSM matter particles, four light higgs doublets (two up-type and two down-type) and two extra charged singlets (δ ++ andδ ++ ). The β-function coefficients are This should be compared with the (-3, 1, 33/5) coefficients of the ordinary MSSM.
-v W < µ < M SU SY : we have the particle content of the standard model, with the exception that there are four Higgs doublets. We have therefore the following βfunction coefficients This should be compared with the (-7, -19/6, 41/10) coefficients of the ordinary SM.

E Yukawa running
As we already reported in section 5, because of the closeness of the intermediate scales between M GU T and M T we can consider only the running between M T and M SU SY to provide analytical approximations for the evolution of fermion masses and mixing under the RGEs effect. At M T the scalar SU (2) L triplet has been already integrated out giving rise to the effective Weinberg operator given by with h u 1 = h u , h u 2 = h u and α ij coefficients arising by the scalar potential and Y L given in eq. (83).
For what concerns the charged fermion Yukawa, at M T the Dirac part of the superpotential is written as with the Yukawa mass matrices given in eq. (85).
The Yukawa matrices RGEs are therefore given by