Modular Invariant Quark and Lepton Models in Double Covering of $S_4$ Modular Group

We perform a comprehensive analysis of the homogeneous finite modular group $\Gamma'_4\equiv S'_4$ which is the double covering of $S_4$ group. The weight 1 modular forms of level 4 are constructed in terms of Dedekind eta function, and they transform as a triplet $\mathbf{\hat{3}'}$ of $S'_4$. The integral weight modular forms until weight 6 are built from the tensor products of weight 1 modular forms. We perform a systematical classification of $S'_4$ modular models for lepton masses and mixing with/without generalized CP, where the left-handed leptons are assigned to triplet of $S'_4$ and right-handed charged leptons transform as singlets under $S'_4$, and we consider both scenarios where the neutrino masses arise from Weinberg operator or type I seesaw mechanism. The phenomenological implications of the minimal models for lepton masses, mixing angles, CP violation phases and neutrinoless double decay are discussed. The $S'_4$ modular symmetry is extended to quark sector, we present several predictive models which use nine or ten free parameters including real and imaginary parts of $\tau$ to describe quark masses and CKM mixing matrix. We give a quark-lepton unified model which can explain the flavor structure of quarks and leptons simultaneously for a common value of $\tau$.


Introduction
The origin of fermion masses and the mixing matrices is one of the greatest challenges for modern particle physics. Neutrino oscillation provides new clues for the understanding of the flavor problem. It is known that neutrino mixing angles show a pattern which is completely different than that of quark mixing: all quarks mixing angles are small, while for the lepton sector two mixing angles θ 12 , θ 23 are large, the third one θ 13 is small and it is comparable to the size of the quark Cabibbo mixing angle [1]. The evidence of CP violation in neutrino oscillation is reported recently [2]. Given the successful use of symmetries in various fields of physics, it was conceived that the flavor structure of quarks and leptons is dictated by certain flavor symmetry, and different kinds of flavor symmetry groups (abelian, non-abelian, continuous, discrete, global, local, linearly or non-linearly realized) have been considered so far. In particular, it turns that the discrete non-abelian flavor symmetry is quite suitable to reproduce the large lepton mixing angles, a huge number of models have been constructed, see [3] for recent review. If discrete flavor symmetry is combined with generalized CP symmetry [4,5], one can predict leptonic CP violation phase. It is notable that a unified description of the observed structure of the quark and lepton mixing can be achieved if the flavor and CP symmetries are broken down to Z 2 × CP in neutrino, charged lepton, up quark and down quark sectors, and the minimal flavor group is the dihedral group D 14 [6][7][8][9].
In any realistic model based on discrete flavor symmetry, the flavor symmetry is spontaneously broken by the vacuum expectation values (VEVs) of a set of scalar fields called flavons which are standard model singlet albeit transforming non-trivially under the flavor symmetry group. The VEVs of flavon are typically aligned along certain directions in flavor space, and the vacuum alignment determines the flavor structure of quarks and leptons. One has to intelligently design the flavon energy density to achieve the required vacuum alignment as the global minimum of the scalar potential. In most models, discrete flavor symmetry is accompanied by additional symmetries, either discrete like Z N or continuous like U (1), to ensure the needed vacuum alignment and to reproduce the observed mass hierarchies. Hence the flavor symmetry breaking sector introduces many independent parameters, makes the flavor model rather complicated. Moreover, high dimensional operators compatible with symmetry the model can lead to corrections to leading order results such the predictability of the moded is spoiled in some sense.
Recently modular invariance as flavor symmetry has suggested to understand the neutrino masses and lepton flavor mixing [10]. Modular symmetry naturally appears in torus and orbifold compactifications of string theory. In this approach, flavon fields are not absolute requirement, the flavor symmetry symmetry can be uniquely broken by the VEV of the modulus τ . Hence the vacuum alignment problem is simplified considerably although a moduli stabilization mechanism is needed. In modular invariant models, the Yukawa couplings transform nontrivially under the modular symmetry and they are just modular forms which are holomorphic functions of τ . In the limit of exact supersymmetry, the superpotenial is completely fixed by modular symmetry. Furthermore, modular invariant models can be quite predictive, typical minimal modular models describe the neutrino masses, mixing angles and CP violating phases in terms of five free real parameters including the real and imaginary parts of τ .
The finite modular group Γ N = SL(2, Z)/Γ(N ) arising from the quotient of the SL(2, Z) modular group by congruence subgroups Γ(N ) have been utilized for the flavor symmetry of quarks and leptons. Some models for lepton masses and flavor mixing have been constructed at level N = 2 [11][12][13][14], level N = 3 [10-12, , level N = 4 [29,[36][37][38][39][40][41][42], level N = 5 [41,43,44] and level N = 7 [45]. The quark masses and mixing parameters can also be addressed by using modular symmetry [18,21,33,35], and the fermion mass hierarchies can naturally arise as a result of a weighton which is a standard model singlet field with non-zero modular weight [34]. Modular symmetry has been discussed in the context of SU (5) grand unification theory [13,17]. It is notable that the dynamics of modular symmetry could be tested at present and future neutrino oscillation experiments [46]. The modular symmetry has been extended to consistently include generalized CP symmetry under which the complex modulus τ transforms as τ → −τ * [47][48][49][50][51]. The interplay between flavor symmetry, CP symmetry and modular invariance was recently analyzed in string theory [48,52]. Extension to the direct product of multiple modular symmetry has been proposed [38,40]. We have generalized the modular invariance approach to include the odd weight modular forms which can be organized into irreducible representations of the homogeneous finite modular group Γ N [23]. Γ N is generally the double covering of the inhomogeneous finite modular group Γ N . Texture zeros of fermion mass matrices can be naturally obtained from Γ N , the masses and mixing of quarks and leptons can be addressed in Γ 3 ∼ = T [44]. The integral weight modular forms provide new opportunity for modular symmetry model building, and Γ 3 ∼ = T has been studied in [23,44]. In the present work, we shall consider the next homogeneous finite modular group Γ 4 ≡ S 4 .
We intend to perform a systematical analysis of lepton models based on S 4 modular symmetry with/without generalized CP. For normal ordering neutrino masses, we find that fifteen viable models which can describe the neutrino masses, mixing angles and CP violation phases in terms of five real parameters |g 2 /g 1 |, arg (g 2 /g 1 ), Re(τ ), Im(τ ), and g 2 1 v 2 u /Λ. After imposing the generalized CP symmetry, the phase arg (g 2 /g 1 ) is constrained to be 0 or π, six out of the fifteen models can produce a good fit to the data. The neutrino mass spectrum tends to be quasi-degenerate in previous models based on inhomogeneous finite modular group Γ N , nevertheless the neutrino masses are much lighter in these S 4 models. Moreover, we extend the S 4 modular symmetry to the quark sector. The rich structure of the modular forms at level 4 allows many possibilities to accommodate the experimental on quark masses and CKM matrix. After scanning the assignments for quark and lepton fields, we find a model which can describe the flavor structure of quarks and leptons simultaneously for a common value of τ .
The remainder of the paper is organized as follows. In section 2, we briefly review the basic aspects of modular symmetry, we show that the inhomogeneous finite modular group Γ N is isomorphic to the quotient of the homogeneous finite modular group Γ N over the center {I, −I}, i.e. Γ N ∼ = Γ N /{I, −I}. The integral weight modular forms at level 4 are constructed up to weight 6 in section 3, and they are arranged into different irreducible representations of S 4 . The generalized CP symmetry compatible with S 4 modular symmetry is discussed in section 4. We find that the generalized CP symmetry requires all the coupling constants real in our working basis. In section 5, we perform a systematical classification of S 4 modular models for lepton masses and mixing, where the left-handed leptons are assigned to triplet of S 4 and right-handed charged leptons transform as singlets under S 4 , and the neutrino masses are described are described by the Weinberg operator or through the type I seesaw mechanism. The S 4 modular symmetry is utilized to address the flavor problem of quark masses hierarchies and CKM mixing matrix, and several models with small number of free parameters are presented in section 6. We give a quark-lepton unification model in section 7, which can explain the masses and mixing patterns of quark and lepton for a common value of τ . Section 8 concludes the paper. Appendix A gives the necessary group theory of S 4 as well as the Clebsch-Gordan (CG) coefficients.

Modular symmetry and modular forms
The modular group Γ is isomorphic to the projective special linear group P SL(2, Z) of 2 × 2 matrices with integer coefficients and unit determinant, where the pairs of matrices A and −A are identified. Hence P SL(2, Z) is the quotient of the two-dimensional special linear group Γ ≡ SL(2, Z) over the integers by its center {I, −I}, i.e., Γ = P SL(2, Z) ∼ = SL(2, Z)/{I, −I}, where I is two dimensional unit matrix. The modular group acts on the upper-half complex plane H = {τ ∈ C|Im(τ ) > 0} by fractional linear transformations, which implies Hence every fractional linear transformation corresponds to a modular group element a b c d , and a b c d and − a b c d represent the same fractional linear transformation. The modular group Γ has infinity group elements which can be obtained as a combination of the two fundamental transformations with the corresponding matrices We check immediately that in Γ we have the relations and also (T S) 3 = I which is equivalent to (ST ) 3 = I if S 2 = I. The corresponding relations in Γ are S 2 = −I, (ST 3 ) = I so that S 4 = (ST ) 3 = I. The Γ orbit of every τ ∈ H has a unique representative in the standard fundamental domain D, which is bounded by the vertical lines Re(τ ) = − 1 2 , Re(τ ) = 1 2 and the circle |τ | = 1 in the upper half plane H. The transformations S and T can map any point in H into the fundamental domain D, and no two points inside D differ by a linear fraction transformation. The transformation T pairs the two vertical lines Re(τ ) = ± 1 2 , and the transformation S maps the arc of |τ | = 1 from i to e πi/3 into the arc from i to e 2πi/3 . The principal congruence subgroup of level N is defined as which are normal subgroups of Γ and Γ respectively. Obviously we have T N ∈ Γ(N ), Γ = Γ(1), Γ = Γ(1), Γ = Γ/{±I}, Γ(2) = Γ(2)/{±I}. For N > 2, we have −1 = 1 (mod N ) and thus I / ∈ Γ(N ), consequently Γ(N ) = Γ(N ). The finite modular group is the quotient of modular group over its principal congruence subgroup [23,53], inhomogeneous finite modular group : Γ N ≡ Γ/Γ(N ) , homogeneous finite modular group : We see Γ 2 ∼ = Γ 2 , Γ N is the double covering of Γ N for N > 2 and Γ N is isomorphic to the quotient of Γ N over its center {I, −I}, i.e., Γ N ∼ = Γ N /{I, −I} [23]. Hence Γ N has double the number of group elements as Γ N with |Γ N | = 2|Γ N |. The homogeneous finite modular group Γ N can be obtained from Γ N by including another generator R which is related to −I ∈ SL(2, Z) and commutes with all elements of the SL(2, Z) group. For N ≤ 5, the multiplication rules of the finite modular groups are 1 [23], It is remarkable that Γ N and Γ N for N ≤ 5 is isomorphic to permutation groups and their double covering, i.e., For N > 5, additional relations besides those in Eq. (10) are needed to render the groups Γ N and Γ N finite [45,53].
From the identities of the eta function in Eq. (16), we know that e 1,2,3 (τ ) transform under the actions of S and T as follow, As shown in Eq. (13), it is always possible to choose a set of basis in M k (Γ(4)) such that the basis vectors can be arranged into several modular multiplets which transform in irreducible representations of Γ 4 ≡ S 4 . Thus for the weight 1 modular forms of level 4, solving the condition of Eq. (13), we find the original basis e 1,2,3 (τ ) can be arranged into triplet modular form Y transforming as a triplet3 of S 4 , The weight 5 modular forms of level 4 decompose as2 ⊕3 ⊕3 ⊕3 under S 4 , and they are given by where Y 3 ,I and Y (5) 3 ,II denote two weight 5 modular forms transforming as triplet3 of S 4 , and they can also taken to be any two linearly independent combinations of Y (5) 3 ,I and Y (5) 3 ,II . Finally there are 13 independent weight 6 modular forms of level 5, and they can be arranged into the following irreducible representations of S 4 , We note that the results of even weight modular forms obtained here coincide with those of our previous work [29] up to some overall constants. Specifically the following relations are fulfilled: 4 Generalized CP consistent with S 4 modular symmetry In order to consistently implement CP symmetry in the context of modular symmetry, the complex modulus τ should transform under the action CP as [47][48][49][50][51] up to modular transformations. A generic chiral superfield Φ(x) assigned to an irreducible representation r of the finite modular group Γ N transforms under the action of Γ N as where −k is the modular weight of Φ. We impose CP symmetry on the modular invariant theory. A generalized CP transformation acts on the chiral superfield Φ(x) as where Px = (t, − x), and a bar denotes the hermitian conjugate superfield, X r is not necessarily diagonal and it in general acts in a non-trivial way on the flavor space. As has been shown in [47], constraints on the choice of X r arise from the requirement that the subsequent application of the CP transformation, the modular symmetry and the inverse CP transformation should be represented by another element of the modular symmetry group, i.e. X r ρ * r (γ)X −1 r = ρ r (u(γ)) , where u(γ) is an outer automorphism of the modular group, Eq. (36) is the so-called consistency condition which CP and modular symmetries have to obey in order to give a consistent definition of generalized CP transformations in setting with modular symmetry. It is notable that the consistency condition Eq. (36) should be satisfied for all irreducible representations of the finite modular group Γ N . We see that the CP transformation X r maps the modular group element γ onto another element u(γ) and the group structure of the modular symmetry is preserved, i.e. u(γ 1 γ 2 ) = u(γ 1 )u(γ 2 ). Hence it is sufficient to impose Eq. (36) on the generators S and T , where the identities u(S) = S −1 and u(T ) = T −1 are used. The consistency condition in Eq. (36) determines the CP transformation X r up to an overall phase for a given irreducible representation r. As regards the double covering group S 4 with the basis given in table 7, solving the consistency conditions of Eq. (38), we find that the generalized CP transformation X r coincides with the representation matrix of S, which is a combination of the modular symmetry transformation S and the canonical CP transformation. Furthermore, we have checked that the modular forms Y Since the theory is required to be modular invariant, the CP transformation in Eq. (39) is essentially equivalent to the canonical CP transformation. As shown in Appendix A, all the CG coefficients in our working basis are real, therefore the generalized CP symmetry would constrain all the couplings in the Lagrangian to be real.

Lepton models based on S 4 modular symmetry
We work in the framework of the modular invariant supersymmetric theory [10,58,59]. In the setting of N = 1 global supersymmetry, the action can be generally written as where K(Φ I ,Φ I , τ,τ ) is Kähler potential, it is real gauge invariant function of the chiral superfields Φ and their hermitian conjugatesΦ. W (Φ I , τ ) refers to the superpotential, and it is a holomorphic gauge invariant function of the chiral superfields Φ. The whole action S should be modular invariant. The transformation properties of Φ I are specified by its modular weight −k I and the representation r I under Γ N , Following Ref. [10], we take the Kähler potential to be the minimal form, where h is a positive constant. After the modulus τ gets a vacuum expectation, this Kähler potential gives the kinetic terms for the scalar components of the supermultiplet Φ I and the modulus field τ . The Kähler potential is strongly constrained in some models based on string theory [60][61][62], and the above minimal Kähler potential as the leading order contribution could possibly be achieved. The superpotential W can be expanded into power series of Modular invariance requires the function Y I 1 ...In (τ ) should be a modular form of weight k Y of level N and in the representation r Y of Γ N : where k Y and r Y should satisfy the conditions In the present work, we shall study the modular symmetry group of level N = 4, and a comprehensive analysis of lepton models with S 4 modular symmetry is performed in the following. Without mentioning explicitly, we shall consider the cases involving weight 1, weight 2 and weight 3 modular forms. The cases with higher weight modular forms can be discussed in the same fashion.

Charged lepton sector
The left-handed lepton doublet fields are assigned to transforms as triplet 3, 3 ,3,3 of S 4 , and we assume the right-handed charged leptons transforms as singlet 1, 1 ,1 or1 of Γ 4 . Thus the most general superpotential for the charged lepton masses can be written as : The modular forms should transform as three dimensional irreducible representations under S 4 , and their explicit forms depend on the weight and representation assignments for L and E c 1,2,3 . In order to charged lepton mass matrix with rank less than three otherwise at least one charged lepton would be massless, and f E 3 (Y ) must be different modular multiplets. For illustration, we consider modular forms of weight less than four, consequently 3 . It is remarkable that the CG coefficients for the contraction triplet⊗ triplet→singlet are all the same in our basis. As a consequence, there are only four different structures of charged lepton mass matrix if the weights of the relevant modular forms are less than four.
In this case, there are four different representation assignments which give rise to the same charged lepton mass matrix: The superpotential for the charged lepton masses are given by, The condition of modular weight cancellation requires There are also four representation assignments for the lepton fields, The superpotential for the charged lepton masses takes the following form, Modular invariance imposes the following constraints on modular weights, Similar to previous cases, the lepton fields can be assigned to The superpotential for the charged lepton masses is of the form, with the modular weights Likewise, we have four different representation assignments which gives the same superpotential W e as well as the the same charged lepton mass matrix, The superpotential for the charged lepton masses reads as, The modular weights k L and k E 1 ,E 2 ,E 3 satisfy the constraints It is straightforward to read out the predicted charged lepton mass matrix for each case discussed above, and results are summarized in table 2. We can exchange the assignments for the right-handed charged lepton fields E c 1,2,3 , accordingly the rows of the charged lepton mass matrix would be permutated. However, the hermitian combination M † e M e are left invariant such that the predictions for charged lepton mass and the unitary rotation U e are unchanged, where U e diagonalize the charged lepton mass matrix via U † e M † e M e U e = diag(m 2 e , m 2 µ , m 2 τ ).

Neutrino sector
In neutrino sector, we assume that neutrinos are Majorana particles, and we consider two scenarios that the neutrino masses are described by the effective Weinberg operator or arise from the type-I seesaw mechanism. The left-handed lepton doublets would be assigned to transform as triplet under S 4 . Guided by the principle of minimality and simplicity, we shall consider modular multiplets with weight less than four similar to the charged lepton sector. For the cases involving higher weight modular forms, more modular invariant operators accompanied by free coupling constants would be allowed, and the predictive power of the models would be reduced.

Weinberg operator
to form a S 4 singlet. At the lowest order, the weight 2 modular multiplets Y enter into the Weinberg operator, and the superpotential for neutrino masses are as follows.
The modular weight k L should be equal to 1, i.e. k L = 1. From the CG coefficients of 3 ⊗ 3 → 3 and 3 ⊗ 3 → 3, we know that the contraction (LL) 3 is a antisymmetric combination of lepton fields L, while Lorentz invariance requires that the Majorana mass term ((LL) 3 Y 3 ) 1 H u H u should be symmetric with respect to L. As a result, the term proportional to g 2 is vanishing, and the corresponding neutrino mass matrix M ν read as where with the weight k L = 1. The light neutrino mass matrix M ν is of the form

Type-I seesaw mechanism
Three generations of right-handed neutrinos are introduced in the present work and they are assumed to transforms as triplet under the S 4 . Then the most general superpotential in neutrino sector can be written as where f N (Y ) and f M (Y ) are modular multiplets. Similar to the case of Weinberg operator, from the Kronecker products of two triplets, we know that f M (Y ) can be τ −independent constant 2 or weight 2 modular form.
In this case, the right-handed neutrinos can transform as 3 or 3 under S 4 (i.e., ρ N c = 3 or 3 ) and their modular weight should be vanishing with k N c = 0. The heavy neutrino mass term is which leads to the following heavy neutrino mass matrix, If the right-handed neutrinos are assigned to transform as unhatted triplet ρ N c = 3 or 3 with k N c = 1, we have 3 ] 1 is vanishing because the contractions for both 3 ⊗ 3 → 3 and 3 ⊗ 3 → 3 are antisymmetric combinations. The corresponding heavy Majorana mass matrix M N can be easily read out as On the other hand, we can also assign the right-handed neutrinos to hatted triplets ρ N c =3 or3 with k N c = 1. Then the superpotential W N is There are no non-trivial modular forms of weight zero.
which gives rise to Now we proceed to discuss the neutrino Yukawa interaction term g (N c LH u f D (Y )) 1 . The modular form f D (Y ) is fixed by the assignments for L and N c , it can 1, Y up to weight 3. We shall report the predictions for Dirac neutrino mass matrix for each possible cases.
In this case, left-handed lepton doublet L and right-handed neutrinos N contract to a singlet, hence their assignments can be (ρ N c , ρ L ) = (3, 3) or (3 , 3 ) or (3,3 ) or (3 ,3). The Dirac neutrino mass term is with k N + k L = 0. Consequently the Dirac neutrino mass matrix read as There are eight possible assignments for ρ L and ρ N c , and they can divided into two categories. In the case of (ρ N c , ρ L ) = (3,3), (3 ,3 ) or (3, 3) , (3 , 3 ), we have with the modular weights k N c + k L = 1. We can read out the Dirac neutrino mass matrix is For the second type of assignments (ρ N c , ρ L ) = (3,3 ), (3 ,3), which leads to The modular weights of L and N c should compensate that of f D (Y ), they satisfy the condition k N c + k L = 2. For the assignments (ρ N c , ρ L ) = (3, 3) , (3 , 3 ) , (3,3 ) or (3 ,3), we have Accordingly the Dirac neutrino mass matrix is of the following form We can also assign N c and L to the S 4 triplets (ρ N c , ρ L ) = (3,3), (3 ,3 ) , (3,3 ) or (3 , 3), and thus The Dirac neutrino mass matrix read as The weight cancellation requires k L and k N c fulfill the condition k N c +k L = 3. Invariance of the neutrino Yukawa coupling under S 4 entails N c and L should contract to1,3 and 3 . Therefore N c and L can be assigned to (ρ N c , ρ L ) = (3,3), (3 ,3 ), (3,3) or (3 , 3 ), then the superpotential W D is of the form which gives rise to Cases Eq. (70)   Summary of neutrino models with less than four free parameters excluding τ , W 1,2 and S i (i = 1, . . . , 16) denote the models in which neutrino masses arise from Weinberg operator and type-I seesaw mechanism respectively. We can also assign N c and L to transform as (ρ N c , ρ L ) = (3,3 ), (3 ,3), (3 , 3) or (3, 3 ), then we have The Dirac neutrino mass matrix is determined to be For all the above type-I seesaw models, the effective light neutrino mass matrix is given by the seesaw formula, We are interested in the models with less free parameters, and we list the possible neutrino models in table 3 for which the resulting light neutrino mass matrices contain less than four free parameters excluding the modulus τ .

Numerical results
In short, the charged lepton can take four possible forms shown in table 2 if only modular forms of weight less than 4 are considered, and there are eighteen neutrino models with parameters less than 4, as summarized in table 3. Combining charged lepton sector with neutrino sector, we obtain totally 4 × 18 = 72 lepton models which are denoted as C i -W 1 , C i -W 2 and C i -S j with the indices i = 1, . . . , 4, j = 1, . . . , 16. We see that for the four cases C 1,2,3,4 the charged lepton mass matrix M e depends on three parameters α, β and γ which can be made real by redefining the phases of the right-handed charged leptons E c 1,2,3 . The three parameters α, β and γ are in one-to-one correspondence with the charged lepton masses. The electron, muon, tau masses can be reproduced by adjusting the parameters α, β and γ. We confront each model with the neutrino oscillation data and charged lepton masses, we perform a conventional χ 2 analysis to optimize the model parameters and determine how well each model can be compatible with the observations. The overall mass scale αv d in the charged lepton mass matrix and g 2 v 2 u /Λ in the neutrino mass matrix can be fixed by requiring the electron mass and the mass splitting ∆m 2 21 are reproduced. Since the overall factor of the mass matrix doesn't affect the predictions for mass ratios, mixing angles and CP violating phases, we construct the χ 2 function using the lepton mixing angles θ 12 , θ 13 , θ 23 and the mass ratios m e /m µ , m µ /m τ , ∆m 2 21 /∆m 2 31 . The neutrino oscillation parameters are taken from the latest global fit results of NuFIT v4.1 including the atmospheric neutrino data from Super-Kamiokande [63]. Since the current data somewhat prefer normal ordering (NO) over inverted neutrino (IO) mass ordering, we shall focus on NO neutrino masses in the numerical analysis. The best fit values and 1σ ranges of the three lepton mixing angles CP violating phase δ CP and the neutrino mass squared differences are as follows The ratios of charged lepton masses taken from [64], m e /m µ = 0.0048 ± 0.0002, m µ /m τ = 0.0565 ± 0.0045 The leptonic Dirac CP phase δ l CP is not measured precisely at present and the indication of a preferred value of δ l CP from global data analyses is rather weak, we don't include the information of δ l CP in the χ 2 function. The absolute values of all coupling constants are scanned in the region [0, 10 4 ] and the phases are freely varied in the range [0, 2π], and the modulus τ is restricted in the right-hand part of the fundamental domain D with 0 ≤ Re(τ ) ≤ 0.5, the reason of why not scanning the complete fundamental domain is explained below. We numerically minimize the χ 2 function by using the minimization algorithms incorporated in the package MINUIT developed by CERN to determine the optimum values of the input parameters. We find that 15 models can give very good fit to the data for certain values of input parameters. We display the best fit values of the input parameter for which the χ 2 function reach a global minimum χ 2 min in table 4, and we also give the predictions for lepton mixing parameters and neutrino masses at the best fitting point in table 4. For all the 15 phenomenologically viable models, the light neutrino mass matrix M ν depends on a single complex parameter g 2 /g 1 and the complex modulus τ besides the overall scale g 2 1 v 2 u /Λ. Hence the three lepton mixing angles, Dirac and Majorana CP phases and three light neutrino masses are completely determined by five real parameters |g 2 /g 1 |, arg (g 2 /g 1 ), Re(τ ), Im(τ ), and g 2 1 v 2 u /Λ. The number of free parameters is four less than the that of observables, therefore these models are quite predictive. Note that the models C i -S 9 and C i -S 10 contain more free parameters, consequently we don't not show the numerical results of these model here. As can be seen from table 4, the three lepton mixing angles θ 12 , θ 13 , θ 23 and the neutrino mass squared difference ∆m 2 21 , ∆m 2 31 fall in the 1σ experimental range for these 15 viable models except the models C 1 -S 3 and C 4 -S 3 where sin 2 θ 23 = 0.4949 is outside the 1σ region but still in the 3σ region.
Furthermore, we notice that the modular forms have the property Y (k) r (−τ * ) = ρ r (S)[Y (k) r (τ )] * as shown in Eq. (40). Therefore if we make the replacment τ → −τ * , g 1,2 → g * 1,2 and perform the S transformation on both lepton and right-handed neutrino fields, the charged lepton and neutrino mass matrices would become their complex conjugate. Therefore under such transformation, lepton masses and mixing angles are unchanged while the signs of all CP violating phases are flipped. As a consequence, the complex modulus τ is limited in the right-hand part of the fundamental domain D with 0 ≤ Re(τ ) ≤ 0.5 when we scan over the parameter space in numerical minimisation. The predictions of the mixing parameters in the left-hand part of D with −0.5 ≤ Re(τ ) ≤ 0 can be easily obtained by reversing the overall signs of the Dirac and Majorana CP phases. Hence all the numerical results given in table 4 should understand to come in pair with opposite CP violating phases.
It is known that the neutrino mass spectrum tends to be nearly degenerate in modular invariant models based on inhomogeneous finite modular group. As can be seen from table 4, a remarkable feature of these modular S 4 models is that the neutrino masses are hierarchical except the models C 2 -W 2 and C 3 -W 2 . From the predictions for neutrino masses, mixing angles and CP violating phases in table 4, we can pin down the effective neutrino mass |m ee | relevant to neutrinoless double beta decay. We displayed the the lightest neutrino mass and |m ee | of each viable model in figure 1, where the experiemental bound of KamLAND-Zen [65] and the expected sensitivities of future experiments [66][67][68][69][70][71] are indicated by the horizontal lines. For the models C 1 -S 6 and C 3 -S 6 , the effective Majorana mass are |m ee | 2.864 × 10 −6 meV and |m ee | = 8.620 × 10 −7 meV respectively with the lightest neutrino mass m 1 6.953 meV. Hence the corresponding points are not visible in the figure. As indicated in the figure, we expect that some of these predictions will be tested in next generation neutrinoless double beta decay experiments.
If we require the theory to be invariant under both S 4 modular symmetry and the generalized CP symmetry, all the couplings would be restricted to be real in our working basis, as shown in section 4. Thus the number of free parameters in a model would be reduced further. For the 15 viable models listed in table 4, the generalized CP symmetry enforces both coupling constants g 1 and g 2 to be real. As a consequence, the light neutrino mass matrix would be described by four real parameters g 2 /g 1 , Re(τ ), Im(τ ) and the overall scale g 2 1 v 2 u /Λ. In particular, the complex modulus τ and modular forms would be sources all CP violation phases. We find only six out of the 15 models are compatible with data, the numerical results are shown in table 5. In order to show the predictive ability of the model after adding generalized CP symmetry, for comparison, we take the model C 1 -S 5 as an example and give detailed numerical scanning results , the correlation graph with or without generalized CP symmetry are shown below. From these figure 3, figure 4, figure 5, figure 6 and figure 7, it can be seen that for the model C 1 -S 5 , the points compatible with experiment data are all located near the boundary of the fundamental domain D, and there are five independent and disconnected regions in the parameter space. The correlation between some observables in each region is very strong, and the patterns of different regions show very different behaviors. After gCP symmetry is added, only two independent and disconnected regions are left in the allowed parameter space of model C 1 -S 5 . And they all shrink in the τ plane, one of them is very small, almost a point. These shows the strong restriction of gCP on the parameter space.

Quark models based on S 4 modular symmetry
In this section, we will exploit the S 4 modular symmetry to understand the quark mass hierarchies and the observed pattern of hierarchial quark mixing angles and CP violating Models Best fit values of the input parameters for NO χ 2 min Re τ Im τ β/α γ/α |g 2 /g 1 | arg (g 2 /g 1 )/π αv d /MeV  phase encoded in the CKM matrix. We aim to construct viable quark mass models with a minimal amount of free parameters. The quark fields can be assigned to triplet of S 4 , the direct product of a doublet and a singlet, or the direct sum of three singlets. Similar to what we have done in the charged lepton sector, we can classify the structures of the quark mass matrix for each assignment. We have constructed the tens of thousands possible quark models by using Wolfram Mathematica, and then find the best fit points by using the package MINUIT.
The modular symmetry is broken by the vacuum expectation value of the modulus τ at high energy scale. We assume that modular invariance breaking scale is around the GUT scale 2 × 10 16 GeV. From the up type and down type quark mass matrices M u and M d we can calculte the quark masses, mixing angles and CP violation phase in terms of the input  [63]. The vertical grey exclusion band denotes the bound on the lightest neutrino mass coming from the cosmological data Σ i m i < 0.120 eV at 95% confidence level obtained by the Planck collaboration [72]. The values of |m ee | in the models C 1 -S 6 and C 3 -S 6 are too tiny to be visible. parameters of the model. In order to find the point in parameter space which optimizes the agreement between predictions and data, we generalize the numerical analysis strategy to the quark sector, and we search the minimum of the χ 2 contributions from quark mass ratios and CKM parameters. For the calculation of χ 2 min , we use the values of Yukawa couplings and the CKM parameters calculated at the GUT scale from a minimal SUSY breaking scenario, with SUSY breaking scale M SUSY = 1 TeV and tan β = 7.5,η b = 0.09375 [73], whereη b denotes the contribution from SUSY threshold corrections which mainly affects the bottom quark Yukawa coupling. After examining many possible constructions, we succeeded in finding some models which can accommodate the experimental data of quark masses and CKM matrix. In the following, we present eight benchmark models with small number of free parameters. The transformation properties of the quark fields under S 4 and their modular weights are summarized in table 6.
1 ○ Model I with gCP: 9 free real parameters including Re(τ ) and Im(τ ) In this model, left-handed quarks Q and right-handed up quarks u c , c c , t c are assigned to a direct sum of doublet and singlet 2 ⊕ 1 of S 4 , the right-handed down quarks d c , s c , b c are assigned to2 ⊕ 1 of S 4 . For convenience, we use the subscript "D" to denote the doublet assignment, i.e.
The modular weights of these quark superfields are set to Models with gCP  Table 5: The best fit values of the input parameters after imposing generalized CP. We give the predictions for lepton mixing parameters and neutrino masses at the best fit points. Notice in the CP dual point τ → −τ * , the signs of Dirac and Majorana CP phases are reversed while the predictions for lepton mixing angles and neutrino masses are unchanged.
Thus the modular invariant superpotentials for quark masses read as follows,  The quark mass ratios and mixing parameters are determined to be 2 ○ Model II with gCP: 9 free real parameters including Re(τ ) and Im(τ ) In this model, the representation assignments for quark fields are the same as Model I except changing the right-handed quarks t c , b c to 1 , 1 of S 4 . The modular weights of these quark superfields are set to Thus the modular invariant superpotentials for quark masses read as follows, From the CG coefficients of S 4 group in Appendix A, we find the up and down quark mass matrices are given by  We see that this model make use of five real parameters α u1,d , β t,d , γ d and two complex parameter α u2 , γ u to describe quark masses and CKM matrix. If we impose gCP symmetry on this model, α u1 , γ u are either positive or negative. A good agreement between data and predictions is obtained for the following values of input parameters The quark mass ratios and mixing parameters are determined to be 3 ○ Model III with gCP: 9 free real parameters including Re(τ ) and Im(τ ) In this model, the Irrep assignments are same as Model I except that the assignments of the right-handed up quarks and down quarks are interchanged. The modular weights of these quark superfields are set to Thus the modular invariant superpotentials for quark masses read as follows, From the CG coefficients of S 4 group in Appendix A, we find the up and down quark mass matrices are given by We see that this model make use of five real parameters α u,d1 , β u,d , γ u and two complex parameter α d2 , γ d to describe quark masses and CKM matrix. If we impose gCP symmetry on this model, α d2 , γ d are either positive or negative. A good agreement between data and predictions is obtained for the following values of input parameters The quark mass ratios and mixing parameters are determined to be 4 ○ Model IV with gCP: 9 free real parameters including Re(τ ) and Im(τ ) In this model, we assume the left-handed quarks Q transform as triplet 3 of S 4 , the righthanded up quark fields u c , c c , t c transform as a direct sum of doublet and singlet 2 ⊕1 ,, the the right-handed down quark fields d c , s c , b c transform as2 ⊕ 1 of S 4 . The modular weights of these quark superfields are set to Thus the modular invariant superpotentials for quark masses read as follows, From the CG coefficients of S 4 group in Appendix A, we find the up and down quark mass matrices are given by The quark mass ratios and mixing parameters are determined to be

5
○ Model V with gCP: 9 free real parameters including Re(τ ) and Im(τ ) In this model, we assume the left-handed quarks Q transform as triplet 3 of S 4 , the righthanded up quark fields u c , c c , t c transform as a direct sum of doublet and singlet2 ⊕ 1 ,, the the right-handed down quark fields d c , s c , b c transform as 2 ⊕1 of S 4 . The modular weights of these quark superfields are set to Thus the modular invariant superpotentials for quark masses read as follows, From the CG coefficients of S 4 group in Appendix A, we find the up and down quark mass matrices are given by fields d c , s c , b c transform as1,1 ,1 of Γ 4 . The modular weights of these quark superfields are set to Thus the modular invariant superpotentials for quark masses read as follows, The phases of α u , β u , γ u , α d , β d and γ d1 can be absorbed into the quark fields while the phase of γ d2 can not be removed by field redefinition. From the CG coefficients of S 4 group in Appendix A, we find the up and down quark mass matrices are given by We see that this model make uses of six real parameters α u,d , β u,d , γ u,d1 and one complex parameter γ d2 to describe quark masses and CKM matrix. A good agreement between data and predictions is obtained for the following values of input parameters The quark mass ratios and mixing parameters are determined to be which are compatible with the experimental data in Eq. (88) except that θ q 23 is somewhat larger. Notice that the top and bottom quark masses can be reproduced by adjusting the parameters α u and α d .

7
○ Model VII with gCP: 10 free real parameters including Re(τ ) and Im(τ ) The left-handed quarks Q are assigned to triplet 3 of S 4 , u c , c c and t c transform as 1, 1 and 1 respectively under S 4 , down type quarks d c , s c , b c transform as1, 1, 1 respectively. Note that u c and c c are distinguished by their different modular weight, and similarly for s c and b c . We choose of the weights of quark fields to fulfill The superpotentials of the quark sector are given by, 3,II ) 1 H d ,(119) which lead to the quark mass matrices, symmetry can also help to understand the quark mass hierarchies and CKM mixing matrix. In this section, we shall investigate the quark-lepton unified models which can explain the experimental data of quarks and leptons simultaneously for certain common value of modulus τ . Such kind of models at level N = 3 has been studied [18,21,33,35] in the literature. It is generally not an easy task to construct a quark-lepton unified model. After trying many possibilities, we find out a realistic quark-lepton unified model with small number of free parameters, the quark sector is described by Model VIII of section 6, we extend the lepton models by introducing higher weight modular forms, and the generalized CP symmetry in combination with S 4 modular symmetry is imposed to match the quark sector.
We assign the three generations of left-handed lepton doublets L and of right-handed neutrinos N c to two triplets 3 of S 4 , while the right-handed charged leptons E c 1 , E c 2 and E c 3 transform as1,1 and 1 respectively. We choose the modular weights of lepton fields as k N c = 2, k L = −2, k e c = 3, k µ c = 5 and k τ c = 6. Then the modular invariant superpotential of the lepton sector is given by where all couplings α e , β e , γ e , g 1 and Λ 1,2,3 are real because of the generalized CP invariance. We can read out the lepton mass matrices as follows, We see that the solar mixing angle θ l 12 is within the 2σ experimental region, and all other observables fall in the 1σ ranges. The sum of neutrino masses is determined to be 80.293 meV, this is compatible with the latest bound i m i < 120 meV at 95% confidence level from Planck [72].

Conclusion
The homogeneous finite modular group Γ N provides new opportunity for understanding the flavor structure of quarks and leptons based on modular invariance. Γ 2 is identical to Γ 2 ∼ = S 3 , Γ N is the double covering of the inhomogeneous finite modular group Γ N for N > 2, and Γ N is isomorphic to the quotient of Γ N over its center {I, −I}, i.e., Γ N ∼ = Γ N /{I, −I}. It is notable that texture zeros of fermion mass matrices can be naturally obtained from Γ N , and Γ 3 ∼ = T has been studied in [23,44]. In the present work, we have considered the modular group Γ 4 ≡ S 4 in the setup of modular invariance approach.
The weight 1 modular forms of level 4 are constructed in terms of the Dedekind eta function, and they can be arranged into a triplet Y (1) ) T which transforms as3 of S 4 . The higher weight modular forms up to weight 6 are built from the tensor products of Y (1) 3 (τ ), and they are homogeneous polynomials of Y 1,2,3 . The odd weight modular forms can be decomposed into the hatted representations1,1 ,2,3 and3 of S 4 while the even weight modular forms can be organized into the other representations 1, 1 , 2, 3 and 3 in common with S 4 . The results are summarized in table 1. Solving the consistency condition, we find the generalized CP transformation corresponding to τ → −τ * is X r = ρ r (S) which is a combination of the modular symmetry transformation S and the canonical CP transformation. All couplings in the Lagrangian would be real if the generalized CP symmetry is imposed.
We perform a systematical analysis of S 4 modular models for lepton masses and mixing with/without generalized CP. We assume that the left-handed leptons transform as triplet of S 4 , and the right-handed charged leptons are assigned to singlets under S 4 , and we consider both the case where neutrino masses are described by the Weinberg operator and the case where neutrino masses arise from the type I seesaw mechanism. The charged lepton mass matrix can only take four possible forms in table 2 and the forms of the neutrino mass  matrices are summarized in table 3 if the weights of the involved modular forms are less than 4. The charged lepton masses m e , m µ and m τ are in a one-to-one correspondence with the parameters α, β and γ which can be taken real without loosing generality. We look for phenomenologically viable models with a minimal amount of free parameters. We find fifteen predictive lepton models which can describe the neutrino masses, mixing angles and CP violation phases in terms of five real parameters |g 2 /g 1 |, arg (g 2 /g 1 ), Re(τ ), Im(τ ) and the overall scale g 2 1 v 2 u /Λ. If generalized CP symmetry is imposed, the phase arg (g 2 /g 1 ) is restricted to be 0 or π. Thus only four real input parameters g 2 /g 1 , Re(τ ), Im(τ ) and g 2 1 v 2 u /Λ are left, and we find six out of the fifteen models can fit the charged lepton masses and neutrino oscillation data very well, as shown in table 5. A remarkable feature of these models is that the light neutrino masses can be very tiny, while the neutrino masses are typically quasi-degenerate in previous models based on Γ N modular group.
We have extended the S 4 modular symmetry to quark sector, different possible assignments (triplet, the direct sum of a doublet and a singlet, or the direct sum of three singlets) of the quark fields under S 4 are considered. Because of the rich structure of the S 4 modular group, we find many models can accommodate the observed patterns of quark masses and CKM mixing matrix. For illustration, we select eight benchmark models in which all the best fit values of observables fall in the 1σ experimental ranges. It is notable that the hierarchical quark masses, quark mixing angles and CP violation phase can even be described very well by models with only nine real parameters including real and imaginary parts of the modulus τ . Note that models at level 3 achieve this with ten [33] or more free parameters [18,21,35].
Finally we present a quark-lepton unified model which can explain the masses and mixing of quarks and leptons simultaneously for a common value of the complex modulus τ . The value of τ is mainly fixed by the precisely measured quark masses and mixing, then the entire neutrino sector including the three neutrino masses as well as the lepton mixing matrix only depends on three real parameters Λ 2 /Λ 1 , Λ 3 /Λ 1 and g 2 v 2 u /Λ 1 . We summarize that S 4 modular symmetry is a promising framework to understand the flavor structure of quarks and leptons.
Note added: During the final preparations of this work, a paper [74] dealing with the same topic appeared on the arXiv. We use different representation basis of S 4 , the Clebsch-Gordan coefficients in our basis is simpler, and this basis is convenient to classify the S 4 modular models. Modular forms of level N = 4 up to weight 6 are constructed in this work, and higher weight modular forms until weight 10 are given in [74]. The authors of [74] present one Weinberg operator model and one type I seesaw model, and the right-handed charged leptons E c are assigned to a triplet of S 4 in [74]. We perform a systematical classification of modular S 4 symmetry models for leptons with/without generalized CP symmetry, E c are assumed to transform as singlets under S 4 in this work. We also apply the S 4 modular symmetry to the quark sector, and construct a quark-lepton unified model. Our work significantly extends the model construction of [74].
where kC n denotes a conjugacy class with k elements of order n. Note that one half of these conjugacy classes can be written as the product of the other half with R. There are four one-dimensional irreducible representations 1, 1 ,1 and1 , two two-dimensional irreducible representations 2 and2, and four three-dimensional irreducible representations 3, 3 ,3 and 3 . We have summarized the explicit matrix representations in table 7. In the representations 1, 1 , 2, 3 and 3 , the generator R = 1 is identity matrix, the representation matrices of S and T coincide with those of S 4 [29], consequently S 4 can not be distinguished from S 4 in these representations since they are represented by the same set of matrices. In the representations1,1 ,2,3 and3 , the generator R = −1. The character table of S 4 can be obtained directly as shown in table 8. Moreover, the Kronecker products between all irreducible representations are given as follows:  Corresponding to the above direct product rule, we give the Clebsch-Gordan (CG) coefficients of S 4 one by one in our working basis. All CG coefficients can be expressed in the form of α ⊗ β, we use α i (β i ) to denote the component of the left (right) basis vector α(β).
For direct products involving singlet1 ( )( ) , their CG coefficients are as followŝ where we have introduced the notation p to distinguish between different products, it makes the results more compact. The CG coefficients for the direct product involving doublet2