Twin Modular $S_4$ with $SU(5)$ GUT

We discuss the $SU(5)$ grand unified extension of flavour models with multiple modular symmetries. The proposed model involves two modular $S_4$ groups, one acting in the charged fermion sector, associated with a modulus field value $\tau_T$ with residual $Z_3^T$ symmetry, and one acting in the right-handed neutrino sector, associated with another modulus field value $\tau_{SU}$ with residual $Z_2^{SU}$ symmetry. Quark and lepton mass hierarchies are naturally generated with the help of weightons, which are SM singlet fields, where their non-zero modular weights play the role of Froggatt-Nielsen charges. The model predicts TM$_1$ lepton mixing, and neutrinoless double beta decay at rates close to the sensitivity of current and future experiments, for both normal and inverted orderings, with suppressed corrections from charged lepton mixing due to the triangular form of its Yukawa matrix.


Introduction
The Standard Model (SM) of particle physics, though highly successful, does not explain why there are three families of quarks and leptons, and does not account for neutrino mass and mixing. The SM also offers no insight into the pattern of quark and lepton (including neutrino) mass and mixing. Moreover, the three gauge groups of the SM are unrelated, with anomaly cancellation and charge quantisation being unexplained, a seemingly accidental consequence of having complete families of quarks and leptons.
Promising attempts to go beyond the SM in order to address these shortcomings are often based on family symmetry on the one hand, and grand unified theories (GUTs) [1] on the other, where the combination of these two symmetries [2] provides a powerful and constrained framework. Neutrino mass and mixing can be included by the addition of gauge singlet states, the right-handed neutrinos, which together with the type I seesaw mechanism [3][4][5][6][7][8][9], provides an elegant understanding of the extreme smallness of neutrino mass compared to charged fermion masses. However understanding the observed approximate tri-bimaximal (TB) nature of neutrino mixing, in which the lepton mixing matrix has a very symmetric structure with U e3 = 0 and equal non-zero elements in the second and third column, suggests an underlying discrete non-Abelian family symmetry [2,10]. For example, the group S 4 , with generators S, T, U ‡ , is capable of enforcing partial symmetries in the neutrino and charged lepton sector capable of enforcing TB mixing, with T associated with a Z T 3 subgroup preserved in the charged lepton sector, and S, U associated with Z S 2 × Z U 2 controlling the neutrino sector. Following the discovery at reactor experiments that U e3 takes a non-zero value [11], exact tribimaximal mixing became excluded, while the larger mixing measured initially by atmospheric and solar neutrino experiments, and later by long baseline oscillation experiments, continue to take the tribimaximal form. This observation suggests possible relaxed forms of mixing known as trimaximal (TM) mixing, in which the first or second column of the lepton mixing matrix maintains the tri-bimaximal form, while the third column is unconstrained, allowing a non-zero element U e3 . There is a phenomenological preference for the first case, called TM 1 mixing [12,13], in which follows from S 4 where T associated with a Z T 3 subgroup preserved in the charged lepton sector, as usual, while the product SU associated with Z SU 2 controlling the neutrino sector. Despite the motivation from neutrino physics for introducing a discrete non-Abelian family symmetry such as S 4 , this raises the question of the origin of such symmetry, and also the nature of its spontaneous breaking, with particular subgroups preserved in different sectors of the theory, all of which seems to require a large number of flavon fields whose vaccuum alignment requires further driving fields, plus ad hoc shaping symmetries [12,13]. The S 4 , for example, could originate from a continuous non-Abelian symmetry [14][15][16][17][18][19][20][21], or have an extra-dimensional origin [22][23][24][25][26][27][28][29][30][31][32][33], being an accidental symmetry of orbifolding [30,[34][35][36], or as a subgroup of the so-called modular symmetry [37,38]. Indeed it has been suggested that a finite subgroup of modular symmetry may be well suited to describe neutrino mass and mixing [39][40][41], where modular forms play the role of the vacuum alignments of flavon models, but without requiring flavons or driving fields [42].
In general the complex modulus τ can take any complex value in the upper half complex plane, but it can also be restricted to take particular values, associated with preserved residual symmetries of the finite modular group, which is referred to as the stabiliser in the framework of finite modular symmetry [51,60,73]. However it is clear from the previous discussion that realistic neutrino models require more than one stabiliser to generate interesting models. This implies that more than one modular symmetry must be considered requiring the framework to be extended to multiple modular symmetries [74]. In a recent paper [75] we considered an example with twin modular S 4 symmetries, leading to TM 1 mixing. The idea was that one modular S 4 controls the neutrino sector, associated with a modulus field value τ SU , while another S 4 acts in the charged lepton sector, associated with a modulus field value τ T . Apart from the predictions of TM 1 mixing, the model led to a novel neutrino mass sum rule with stringent lower bounds on neutrino masses close to current limits from neutrinoless double beta decay experiments and cosmology.
Such models, while interesting, hitherto do not address the question of quark mass and mixing, nor the possibility of extending the framework of modular symmetry to GUTs. Indeed there has been progress of both questions in the literature, with modular models of both quark and lepton mass and mixing having been proposed [53,65,82]. In one of these approaches [65] the quark and lepton mass hierarchies were explained in the framework of modular symmetry without introducing any additional symmetry such a the Froggatt-Nielsen (FN) U (1) [83] which is commonly used to account for charged fermion mass hierarchies. It was pointed out that the role of the FN flavon may be played by a singlet under both the gauge symmetry and the finite modular symmetry which carries a non-zero modular weight, called the weighton [65], leading to natural quark and lepton mass hierarchies. The question of quark and lepton mass and mixing is necessarily simultaneously addressed in GUTs, and various SU (5) models have been proposed [45,49,[84][85][86] which include finite modular symmetries which act as family symmetries. Modular symmetry in the context of SU (5) GUTs was first studied in an (Γ 3 A 4 )×SU (5) model in [49], then (Γ 2 S 3 )×SU (5) [45,84], and (Γ 4 S 4 )×SU (5) [85]. Most recently a comprehensive study of (Γ 3 A 4 ) × SU (5) has been peformed [86]. All these models consider a single modular group.
In this paper, we shall propose an SU (5) GUT with Trimaximal TM 1 neutrino mixing arising from twin modular S 4 groups. The resulting model is based on a similar strategy that we developed for our previous lepton model [75], and indeed the resulting neutrino mass and Yukawa matrices are identical to those considered previously. However, the quark and charged lepton sectors are quite different, with both necessarily having a non-diagonal structure in order to account for quark mixing, unlike the previous model in which the charged lepton Yukawa matrix was diagonal. Also, the hierarchy of the charged lepton masses which was previously obtained by the tuning of Yukawa couplings, is now accounted for by the weighton fields, associated with the twin S 4 symmetries, along with the quark mass hierarchies and quark mixing parameters which of course were not considered at all before, since the previous model was purely a model of leptons. Thus we shall construct a flavour model based on SU (5) with two modular symmetries S F 4 × S N 4 , broken to a single S 4 by a bi-triplet scalar, leading to the effective theory invariant under one single S 4 but involving two modulus fields. The modulus fields gain different VEVs, leading to the breaking of S 4 to a residual Z T 3 in the charged fermions sector, and a residual Z SU 2 symmetry in the neutrino sector. We analyse the resulting model numerically and show that, in the quark sector, the correct CKM matrix can be achieved, while in the lepton sector, there are small deviations from the TM 1 mixing. These corrections are small, but may be tested in the future high-precision neutrino oscillation experiments.
The layout of the remainder of the paper is as follows. In section 2, we briefly review multiple S modular symmetries as the direct origin of fermion masses and mixing. In section 3, we propose the SU (5) GUT model based on S F 4 × S N 4 with two moduli fields. We perform a numerical scan in section 4 to see how good the model match with experimental data. Section 5 concludes the paper. The Appendix A includes some background information about the modular group and modular forms of level 4, while Appendix B gives some details about the vacuum alignment required in our model.

Multiple S 4 modular symmetries
We review the approach of using multiple finite modular symmetries to explain flavour mixing, with Γ 4 S 4 as the representative in the following discussion.
We give a brief introduction of S 4 as a finite modular group here. For more details of Γ 4 and its connection with the infinite modular group, see Appendix A. In the modular S 4 symmetry, each element γ acting on the complex modulus τ (Im(τ ) > 0) as a linear fractional transformation: where a, b, c, d are integers mod 4 and satisfy ad − bc = 1. It is convenient to represent each element of the modular S 4 group by a 2 × 2 matrix, namely, The finite modular group has two generators, S τ and T τ , which satisfy S 2 τ = (S τ T τ ) 3 = T 4 τ = 1. § They acting on the modulus τ take the following forms respectively. They could be represented by 2 × 2 matrices as We consider flavour model construction in the modular invariance approach where the modular symmetry is regarded as the direct origin of flavour mixing. Each field φ i , which transforms as an irreducible represenation (irrep) I i of S 4 in the flavour space, has a modular weight 2k i . Under the modular transformation, (5) § Relaxing the last condition T 4 τ = 1 leads to the infinite modular group Γ. Replacing it with T N τ = 1 with N = 2, 3, 5 leads to finite modular groups Γ2 S3, Γ3 A4 and Γ5 A5, respectively. For N > 5, additional restrictions are required to make ΓN finite. For example, at N = 7, one such restriction is (Sτ Tτ ) 4 = 1, leading to the finite Γ7 Σ(168) [81].
The superpotential in an N = 1 supersymmetry can be generically written in powers of the fields as where the subscript 1 is the singlet contraction of the coupling and fields. Here, the coupling Y I Y does not have to be a trivial coefficient, but any modular form with the representation contributed to the singlet contraction and the modular weight k Y must satisfy k Y + k 1 + k 2 + · · · + k n = 0. With these requirements, the superpotential is invariant under any tranformation of S 4 . Besides, k Y is required to be an even number, as an intrinsic property of modular forms. It transforms as under the action of γ. Given suitable arrangements for representations and weights for particles in the flavour space, a modular-invariant flavour model can be constructed by including a set of modular forms {Y I Y } which keep the superpotential invariant under the modular transformation. After the modulus field τ gains a vacuum expectation value (VEV), a specific flavour texture of Yukawa couplings is achieved.
We extend the discussion to theories including more modular symmetries. Without loss of generality and for convenience of the latter discussion, we include two, labeled as S F 4 and S N 4 , respectively. Let us assume S F 4 and S N 4 to be independent of each other, and the moduli fields are denoted as τ F and τ N , respectively. Any modular transformation where the field φ i and coupling Y (I Y,F ,I Y,N ) transform as respectively. Here, we have arranged φ i and Y (I Y,F ,I Y,N ) as matrices such that γ F acts on them vertically and γ N acts on them horizontally. The modular weights 2k Y,F and 2k Y,N are even numbers.
In the perspective of model building, as we have shown in [75], including two modular symmetries 1) allows the modular symmetries to be broken into different residual symmetries in the charged lepton sector and neutrino sector, respectively, and 2) leads to new lepton flavour mixing patten. In the rest of the paper, we extend the multiple modular invariance approach into the quark sector in the framework of grand unification.

An SU (5) model with twin S 4 modular symmetries
Below we construct a flavour model with two S 4 modular symmetries. The particle contents and their representation and modular weight assignments are given in Table 2. In the gauge space, the SM matter fields belong to5-and 10-plets and the latter are represented as The right-handed neutrino field N is introduced to generate neutrino masses. F and N are arranged as triplets in S F 4 and S N 4 , respectively. T 1 , T 2 and T 3 are arranged as singlets of S F 4 and S N 4 but with different modular weights. The SM Higgs is embedded into heavy Higgs multiplets 5, 5 and 45 of SO (10). They are arranged as trivial singlets in the flavour space. We introduce a bi-triplet scalar Φ, which is supposed to break the twin modular S F 4 and S N 4 symmetries into a single S 4 symmetry [75]. We further introduce two scalars φ 1 and φ 2 . These particles are singlets in both the gauge and flavour space but take non-trivial modular weights, and therefore named as weightons by authors in [65]. In the present work, they will be responsible for hierarchies of quark masses and charged lepton masses.
Superpotential terms to generate down quark and charged lepton masses are w ⊃ w5 + w4 5 with and H4 5 , respectively. Here, Y ··· , · · · are triplet (either 3 or 3 ) modular forms of modular weights 2, 4, 6, · · · , respectively. We fix the VEV of τ F at τ F = τ T ≡ ω, which preserve a residual Z T 3 symmetry. Then, we obtain where y e1,5 , y e2,5 , y e3,5 , y µ,5 and y τ,5 are overall factors which are free parameters in the model. Similarly, (τ T ) up to overall factors. These textures lead to Yukawa coupling matrices for downtype quarks and charged leptons in the LR convention as y dd = y e1,5 s + y e1,45 c, y ds = y e2,5 s + y e2,45 c, y db = y e3,5 s + y e3,45 c, y ss = y µ,5 s + y µ,45 c, y bb = y τ,5 s + y τ,45 c, y ee = y e1,5 s − 3y e1,45 c, y µe = y e2,5 s − 3y e2,45 c, y τ e = y e3,5 s − 3y e3,45 c, , and * represents the complex conjugation of matrices as we present the result in the left-right notation. The above triangular form of the Yukawa matrices in the LR convention ensures that charged lepton corrections to PMNS mixing are very suppressed, while allowing the down quark matrix to contribute dominantly to Cabibbo mixing at order 1 . Superpotential terms to generate up-type quark masses are given by where y uu , y cc , y tt , y uc , y ut , and y ut are dimensionless coefficients, and y where y (4) uu (τ T ) = 0 was used. In the neutrino sector, the superpotential terms to generate neutrino masses are nothing new but a GUT extension of our former work [75]. They are given by where the subscript of a bracket, e.g., (N N ) 3 , represents the irreducible representation contraction of the fields in the bracket. Φ is a bi-triplet of S F 4 × S N 4 . Its VEV is not responsible for special Yukawa textures for leptons, but used to break two modular S 4 's to a single modular S 4 symmetry [74], The VEV of Φ takes the form The technique of how to achieve this VEV structure was first developed in [88] and has been applied to achieve the breaking of multiple modular symmetries in [74]. For details of how to derive it in this model, we refer the reader to Appendix B of this paper and [74]. After the breaking S F 4 × S N 4 → S 4 , the resulted Yukawa coupling between left-handed and right-handed neutrinos is given by On the other hand, the Majorana mass matrix for right-handed neutrinos is straightforwardly obtained from the singlet, doublet and triplet modular forms By assuming the stabiliser in the neutrino sector at τ N = τ SU ≡ − 1 2 + i 2 . This stabiliser is invariant under the action of SU . Namely, we are left with a residual Z SU where a, b and c are complex parameters. And the Majorana mass matrix for right-handed neutrinos are written in the form of four matrix patterns, i.e., This is explicitly the same mass matrix as that obtained in [75]. Without considering the correction from the charged lepton sector, a mass texture as in Eq. (24), together with a Dirac neutrino matrix in Eq. (21), leads to the trimaximal TM 1 lepton mixing which preserves the first column on the tri-bimaximal mixing matrix, Three equivalent relations are predicted from the TM 1 mixing, More precisely, we can write U TM 1 explicitly in the form which includes three free parameters θ R , α 1 , and α 2 . The α 3 is not an independent parameter but have a complicated correlation with the rest parameter as seen in [75]. The mixing angles and Dirac-type CP-violating phase are determined to be [13] sin θ 13 = sin θ R √ 3 , From these equation, we see that the CP-violating phase δ and the angle θ 23 are correlated via the phase difference α 1 − α 2 . In the case of θ R = 0 and maximal atmospheric mixing θ 23 = 45 • , the phase α 1 − α 2 must equal ±90 • , leading to maximal CP violation δ = 90 • or 270 • . We know that that the oscillation data shows a small deviation from the maximal atmospheric mixing. Therefore, we expect δ has a small deviation from its maximal CP-violating value. This will be checked numerically in the next section. The mass matrix in Eq. (24) further constrains coefficients for the third and fourth mass matrix patterns on the right hand side with the ratio − 2/3. This is a feature different from classical flavour models without modular symmetry, e.g., in [13], where all parameters in front of different patterns are arbitrary. We have obtained a neutrino mass sum rule for the light neutrino mass eigenvalues Furthermore, we predicted m ee , appearing as the effective mass in neutrino-less double beta decays, correlated with other mass parameters as which is beyond those reported in [89].
In the present work, the GUT extension induces small corrections from the non-diagonal charged lepton Yukawa coupling Y e . We will check how large is the deviation from the TM 1 mixing numerically in the next section.

Numerical study
We perform a χ 2 analysis and discuss numerical prediction of the model.
The charged fermion Yukawa coupling matrices can be diagonalised via For the eigenvalues of the Yukawa couplings, we use the following numerical results, for the normal ordering (NO, i.e., m 1 < m 2 < m 3 ) of neutrino masses and ∆m 2 21 = (7.42 ± 0.21) × 10 −5 eV 2 , ∆m 2 3l = −(2.497 ± 0.028) × 10 −3 eV 2 , for the inverted ordering (IO, i.e., m 2 < m 3 < m 1 ), where ∆m 2 3l = ∆m 2 31 (for NO) and ∆m 2 32 (for IO). To check the phenomenological viability of the model, we perform a simple χ 2 analysis. In the quark sector, the χ 2 function is defined via , Para = {y uu , y uc , y ut , y cc , y ct , y tt , y dd , y ds , y db , y ss , y bb , 1 , 2 } , Here, in the Para set, 1 and 2 are scanned in the range (0, 0.3) and all coefficients y αβ are scanned with absolute values in the range (0, 1) and phases in (0, 2π). Due to the large parameter space, the model can fit the experimental data very well. We have checked that in all points with small χ 2 , a hierarchy 2276 is the Cabibbo angle. A scan plot for 1 and 2 with χ 2 < 1 (blue) and χ 2 < 10 (light blue), respectively, is shown in Fig. 1. Below is a sample of parameters and predictions of observables with χ 2 = 0.001, We discuss the phenomenological prediction in the lepton sector. We have recovered the same neutrino mass matrix as in [75]. Such a texture leads to the TM 1 mixing if the charged lepton mass matrix is diagonal. However, as we see in Eq. (14), some small off-diagonal entries of Y l are allowed by the modular invariance. Therefore, deviations from the TM 1 mixing are expected. Off-diagonal entries of the unitary matrix V e , which is used to diagonalise Y e Y † e (c.f. Eq. (31)), are estimated to be (V e ) 12 ∼ 1 2 , (V e ) 23   entries can maximally reach (V e ) 12 ∼ 4 × 10 −3 , and (V e ) 23,13 ∼ 10 −5 . Therefore, a deviation less than one percent may be induced by (V e ) 12 . We perform a χ 2 analysis to compare the prediction with and without the correction from charged lepton sector. The χ 2 function is defined following Eq. (36), but the set of Obs and Para are given differently, Para = {y ee , y µµ , y τ τ , y µe , y τ e , θ R , α 1 , α 2 , ∆m 2 21 , ∆m 2 3l , 1 , 2 } , Obs = {ỹ e ,ỹ µ ,ỹ τ , ∆m 2 21 , ∆m 2 3l , θ 12 , θ 23 , θ 13 } .
Note that the CP-violating phase δ has been measured to be 195 •+51 • −25 • at the best-fit value ±1σ errors for NO and 286 •+27 • −32 • at the best-fit value ±1σ errors for IO, respectively. We have not included the it in the Obs set as its current global fit is far away from a Gaussian distribution. Instead, we treat δ as a prediction of the model. The lepton sector shares the same weightons and thus have the same parameters 1 and 2 in the scan. All the other free parameters are independent of those in the quark sector. In the charged lepton mass matrix, we can rotate phases of all coefficients {y ee , y µµ , y τ τ , y µe , y τ e } away without loss of generality and scan them in the range (−1, 1). In the neutrino sector, as the mass matrix is explicitly the same as that genrated in [75], we follow the same parametrisation therein and left with five free parameters, one angle θ R , two phases α 1 , α 2 , and two mass square differences ∆m 2 21 and ∆m 2 3l . We scan θ R in [0, π/2), α 1,2 in [0, 2π), and two mass square differences in their 3σ ranges. Points for χ 2 < 10 are abandoned. Note that the two mass square differences appear in both the Para and Obs sets. Points with ∆m 2 21 and ∆m 2 3l in their 3σ ranges but leading to χ 2 > 10 are discarded as required. Prediction of mixing angles θ 13 and θ 12 with χ 2 < 10 are shown in Fig. 2. Predictions without considering corrections from the charged lepton are listed in the figure as a comparison. A deviation of 0.02 • from the TM 1 mixing could be predicted. Therefore, the TM 1 mixing is still a very good approximation. The CP-violating phase δ is listed as a prediction of the model, as shown in Fig. 3. The model in both mass orderings support small deviation from the maximal CP violation, 270 • < δ < 305 • for χ 2 < 10. Below is a sample of parameters and predictions of observables with χ 2 = 4.515 for the NO of neutrino mass ordering, where 1 and 2 take same values as in Eq. (37). The CP-violating phase δ in this sample is predicted to be δ = 285.55 • . We also check the prediction of the effective neutrino mass parameter m ee which is under measurements of the neutrinoless double beta decay experiments. The scanned results are show in Fig. 2 for both the NO and IO of neutrino mass ordering. In particular, the sample point with inputs in Eq. , CUORE [96] and GERDA [97], and cosmological constraints from PLANCK 2018 (disfavoured region 0.12 eV < m i < 0.60 eV and very disfavoured region m i > 0.60 eV) [98] are shown for comparison.
CUORE [96] and GERDA [97] at 90% CL, respectively, are shown in the figure by comparison. We also include the cosmological constraints on the neutrino masses from Planck 2018 [98]. The two vertical grey bands refer to the disfavoured region 0.12 eV < m i < 0.60 eV and very disfavoured region m i > 0.60 eV, respectively. Due to neutrino mass sum rule in Eq. (29), the lightest neutrino mass is constrained to be m lightest 0.03 eV, leading to most of the parameter space disfavoured by the cosmological constraints.

Conclusion
We have constructed an SU (5) GUT with twin S F 4 × S N 4 modular symmetries, accompanied by two moduli fields. This is a grand unified extension of our previous work in [75], allowing quark mass and mixing to be included, while preserving the good predictions in the lepton sector. The two modular symmetries S F 4 × S N 4 , one acting on charged fermions and one on right-handed neutrinos, are broken to a single S 4 symmetry by a bi-triplet scalar, leading to the effective theory invariant under one single S 4 but involving two modulus fields. The two modulus fields gain different VEVs, leading to the breaking of S 4 to a residual modular symmetry Z 3 in the charged fermion sectors, and a residual modular symmetry Z 2 in the neutrino sector, ensuring the leading order TM 1 lepton mixing.
In the neutrino sector, the model has the same flavour structure as that in [75]. In the charged fermion sector, there are two main differences from our former work: 1) charged fermion mass hierarchies are explained due to the coupling to the two weightons; 2) small mixings are induced in the charged fermion mass matrices. The triangular form of the charged lepton and down quark Yukawa matrices plays a special role in this model, ensuring suppressed charged lepton corrections to the PMNS matrix, while allowing the down quark Yukawa matrix to dominantly contribute to Cabibbo mixing, where we have shown that a good fit may be achieved to the CKM matrix. The model predicts TM 1 lepton mixing to very good approximation, and neutrinoless double beta decay at rates close to the sensitivity of current and future experiments, for both normal and inverted orderings. A Modular group and modular forms of level 4 The modular group Γ is a set of modular transformations acting on the upper complex plane. Given a modulus viable τ with Im(τ ) > 0, each element γ of Γ appears as a linear fractional transformation with a, b, c, and d are any integers satisfying ad − bc = 1. As γ can be represented by a 2 × 2 matrix in Eq. (1), the group Γ is expressed as The modular group is an infinite group. It has two generators, S τ and T τ , satisfying S 2 τ = (S τ T τ ) 3 = 1. These generators act on the modulus τ as, respectively, which, represented as 2 × 2 matrices, are given by With the requirement a, d = 1 (mod 4) and b, c = 0 (mod 4), a subset of Γ, which is also infinite, is obtained, Γ 4 is the quotient group Γ 4 = Γ/Γ(4). It is equivalently obtained by imposing the identity T 4 τ = 1. As a subgroup of Γ, its elements can also be represented as 2 × 2 matrices, but the representation matrices are not unique. As the quotient group Γ/Γ(4), each element γ of Γ 4 satisfy the equality where k a , k b , k c and k d are integers and satisfy 4k a k d + ak d + dk a = 4k b k c + bk c + ck b and η = ±1. The finite modular group Γ 4 is isomorphic to S 4 . The latter is the permutation group of 4 objects, see e.g. [99]. It has 5 irreducible representations, 1, 1 , 2, 3 and 3 . In the studies of flavour symmetries, a common set of generators of S 4 are S, T and U which satisfy S 2 = T 3 = U 2 = (ST ) 3 = (SU ) 2 = (T U ) 2 = 1. We work in a widely used basis with representation matrices of S, T and U listed in Table 3. Table 3: Representation matrices for the S 4 generators T , S and U used in the paper, where ω = e 2πi/3 .
As S 4 Γ 4 , S, T and U can be represented by S τ and T τ , and vice versa, Any element of S 4 can be represented by S τ and T τ . For example, SU , which is crucial to achieve the special mass structure in the neutrino sector, is given by SU = S τ T τ S τ T −1 τ S τ . Given representation matrices of S τ and T τ in Eq. (45), it is straightforward to obtain 2 × 2 matrices of S, T , U and SU as We emphasise that the representation matrix for each element is not unique. An alternative representation matrix is obtained by the equality in Eq. (47). A main difference of the modular invariance approach from the classical flavour symmetry approach is that the Yukawa couplings are formed as a consequence of modular forms, instead of combinations of a series of flavon VEVs. Modular forms are holomorphic functions of τ under modular transformations.
Modular forms of level 4 are classified by modular weights. The latter must be even and positive integrals, which we label as 2k. There are 4k + 1 linearly independent modular forms of level 4 and weight 2k. They are all decomposed to irreducible representations of S 4 . For 2k = 2, there are 5 modular forms. They are decomposed into a doublet 2 and a triplet 3 of S 4 , A non-linear algebra satisfied among three of the modular forms, is satisfied [72]. This constraint is essential to cover the modular space of Γ 4 . Modular forms of higher weights are contracted from Y 1 , · · · , Y 5 . For 2k = 4, there are 9 linearly independent modular forms, forming a singlet 1, a doublet 2 and two triplets 3 and 3 of S 4 , Y (4) For 2k = 6, the number is increased to 13. They are decomposed to More modular forms with higher modular weights are listed in [72]. In particular, singlet modular forms at weights 8, 10 and 12 are respectively given by A modular form at a stabiliser takes an interesting weight-dependent direction which satisfies [74] ρ I (γ)Y I (τ γ ) = (cτ γ + d) −2k Y I (τ γ ) .
Namely, a modular form at a stabiliser τ γ is an eigenvector of the representation matrix ρ I (γ) with respective eigenvalue (cτ γ + d) −2k . In the special case (cτ γ + d) −2k = 1, leading to ρ I (γ)Y I (τ γ ) = Y I (τ γ ), the residual modular symmetry is reduced to the residual flavour symmetry. Otherwise, the residual modular symmetry is different from the latter. We discuss which directions triplet modular forms Y (2k) 3 ( γ ) (τ ) may take at τ = τ T and τ SU . The eigenvalue (cτ γ +d) −2k at these stabilisers is given by (−τ T −1) −2k = ω 2k at γ T and (2τ SU +1) −2k = (−1) k at γ SU , respectively. Given triplet (3, 3 ) representation matrices for S, T and U in Table 3, it is straightforward to obtain These results are directly obtained from the symmetry argument, without knowing explicit expressions of modular forms. However, there are some exceptions of modular forms whose directions cannot be directly obtained from the above argument, e.g., Y 3 (τ SU ), which correspond to eigenvectors of degenerate eigenvalues. Their directions have been calculated in [75] with the help of the identity in Eq. (51), For the singlet modular forms, some of them may vanish at the stabilisers, e.g., Y 1 (τ T ) = Y 1 (τ T ) = Y (10) 1 (τ T ) = 0. They do not contribute to fermion masses.

B Vacuum alignments
We discuss how to achieve the required VEVs of the bi-triplet scalar and weightons using the general driving-field approach in the supersymmetry.
The vacuum alignment for the bi-triplet scalar Φ can be realised by introducing two driving fields χ F N ∼ (3, 3) and χ N ∼ (1, 3) of S F 4 × S N 4 with zero modular weights. The general superpotential terms for driving fields are given by where M is a mass-dimensional coefficient. Minimisation of the superpotential gives rise to (ΦΦ) (3,3) As proven in [74], these identities guarantee the VEV of Φ as in Eq. (20) up to an unphysical basis transformation, thus leading to the breaking S F 4 × S N 4 → S 4 . The non-vanishing VEVs of weightons can be realised by introducing two driving fields χ 1 and χ 2 , both of which are trivial singlets of S F 4 × S N 4 with zero modular weights. Superpotential terms are where M 1 and M 2 are two mass-dimensionful parameters, y where we have taken τ F = τ T and τ N = τ SU and used Y