Flavor mixing and CP violation from the interplay of $S_4$ modular group and gCP

We have performed a systematical analysis of lepton and quark masses models based on $\Gamma_4\cong S_4$ modular symmetry with gCP symmetry. We have considered both cases that neutrinos are Majorana particles and Dirac particles. All possible nontrivial representation assignments of matter fields are considered, and the most general form of fermion mass matrices are given. The phenomenologically viable models with the lowest number of free parameters together with the results of fit are presented. We find out nine lepton models with seven real free parameters including the real and imaginary parts of modulus for Majorana neutrinos, which can accommodate the lepton masses and neutrino oscillation data. The prediction for leptogenesis is studied in an example lepton model. The observed baryon asymmetry as well as lepton masses and mixing angles can be explained. For Dirac neutrinos, four lepton models with five real free couplings are compatible with experimental data. Ten quark models containing seven couplings are found to be able to accommodate the hierarchical quark masses and mixing angles and CP violation phase. Furthermore, the $S_4$ modular symmetry can provide a unified description of lepton and quark flavor structure, and a benchmark model is presented.


Introduction
Understanding the hierarchical fermion masses and the flavor mixing structure of quarks and leptons from the first principle is a longstanding challenge in particle physics. The measurement of neutrino mixing parameters provides new clews to the above mentioned flavor puzzle. There is still lack of a guiding principle to explain this flavour puzzle, and one of the most extensively studied schemes is the flavor symmetry, which is traditionally based on continuous Lie groups or discrete finite groups which relate the three generations of fermions. See Ref. [1] for the latest review. In the traditional flavor symmetry, the flavor groups are broken along certain direction in the flavor space by the vacuum expectation value (VEV) of some scalar fields called flavons. In order to realize the desired symmetry breaking pattern, certain shaping symmetry and additional dynamics and fields are generally necessary. As a consequence, the resulting model look quite complex. In order to overcome this drawback, the modular invariance from a bottom-up perspective has been proposed [2]. The modular symmetry plays the role of flavor symmetry, the flavons are replaced by the so-called modulus τ and the Yukawa couplings are modular forms which are holomorphic functions of τ . This framework has been extended to invariance under more general discrete groups and the modular forms become more general automorphic forms [3].
In this scheme, modular symmetry is governed by the infinite discrete group Γ = SL(2, Z). The modular invariant models are classified by the level N which is a positive integer, and the matter fields are assumed to transform in irreducible representations of the finite modular group Γ N ≡Γ/Γ(N ) or its double covering group Γ N ≡ Γ/Γ(N ). For a small number of finite modular groups, some fermion masses models based on the modular invariance and their phenomenology have been studied, such as Γ 2 ∼ = S 3 [4][5][6][7], Γ 3 ∼ = A 4 [2,4,5,, Γ 4 ∼ = S 4 [21,[33][34][35][36][37][38][39][40], Γ 5 ∼ = A 5 [38,41,42], Γ 7 ∼ = PSL(2, Z 7 ) [43], Γ 3 ∼ = T [44,45], Γ 4 ∼ = S 4 [46,47] and Γ 5 ∼ = A 5 [48,49]. There have also been attempts to implement modular symmetry in the Grand Unified Theories (GUTs) to address both the lepton and quark flavor problems [6,10,[50][51][52][53][54]. In the modular invariance approach, the crucial elements are the modular forms of level N . For the even modular weights, the modular forms can be arranged into multiplets of the finite modular group Γ N [2]. If the modular weights are general integers, the modular forms can be organized into irreducible multiplets of the double covering group Γ N [44]. In addition to the integer weight modular forms, there are fractional weight modular forms for some particular level N . Accordingly the modular group SL(2, Z) should be extended to its metaplectic covering group and the finite metaplectic group acts as flavor symmetry [55]. The modular forms of weights k/2 with finite metaplectic modular groupΓ 4 ∼ =S 4 [55] and the weight k/5 modular forms with finite metaplectic modular group Γ 5 ∼ = A 5 × Z 5 [49] have been studied in the bottom-up modular invariance approach. It has been shown that metaplectic flavor symmetries can be derived from compactifications on tori with magnetic background fluxes [56,57]. The VEV of τ is usually treated as a free parameter in modular invariant models in order to match the experimental data, it is remarkable that the hierarchical fermion mass matrices may arise due to the proximity of the modulus to the residual symmetry preserved points τ = i, −1/2 + i √ 3/2, i∞ [32,58,59]. Furthermore, the generalized CP symmetry can be consistently imposed in the context of symplectic modular symmetry for single modulus with g = 1 [60][61][62] and for multi-moduli with g ≥ 2 [63]. Notice that the symplectic group coincides with the modular group SL(2, Z) when g = 1. In the symmetric basis for the modular generators S and T with ρ r (S) = ρ T r (S) and ρ r (T ) = ρ T r (T ), the gCP transformation would coincide with the canonical CP. The gCP invariance enforces all coupling constants to be real if the Clebsch-Gordan (CG) coefficients are also real in the symmetric basis. Thus the gCP symmetry could further reduce the number of free parameters of the modular invariant models and leads to a higher predictive power.
The S 4 modular group has five irreducible representations: two singlets 1 and 1 , a doublet 2, and two triplets 3 and 3 . Similar to the A 4 modular symmetry, the three generations of right-handed lepton fields and right-handed charged leptons are usually assumed to transform as triplet and singlet respectively under S 4 modular symmetry in the known S 4 modular invariant models, and the doublet assignment for the lepton fields has not been considered although it provides new features and possibilities unavailable in the A 4 modular symmetry.
Moreover, S 4 modular symmetry has been used to explain the flavor structure of lepton so far, but it is not clear whether the S 4 modular symmetry can help to address the quark flavor problem except few GUT models [51,53,54]. In this paper, we intend to perform a systematical analysis of lepton and quark models based on Γ 4 ∼ = S 4 modular symmetry and gCP, and we aim at the viable models involving the lowest number of free parameters. For lepton models, we find that thirteen viable models can successfully describe the experimental data of lepton masses and mixing parameters in terms of seven real parameters including Reτ and Imτ . In the quark sector, at least seven real couplings are necessary in order to accommodate the measured values of quark masses and CKM mixing matrix, Furthermore, we find that agreement with the experimental data of quarks and lepton sectors can be achieved for a common value of τ , and a benchmark model is presented. In modular invariant models which also fulfill gCP invariance, the VEV of the modulus τ is the unique source of both modular symmetry breaking and CP violation. Thus imposing gCP symmetry would lead to strong correlations between the low energy CP violation phases and CP asymmetry in leptogenesis. We shall discuss the baryon asymmetry generated via unflavored thermal leptogenesis in an example model of lepton.
This paper is organized as follows. In section 2, we briefly review the modular symmetry and modular forms of level 4. In section 3, we give the most general forms of the Yukawa superpotential and the Majorana mass term for different possible assignments of matter fields, the corresponding mass matrices are presented. Moreover, we show that different assignments can lead to the same fermion mass matrices. In section 4, we find out the phenomenologically viable lepton and quark models with the smallest number of free parameters, and the results of fit are presented. The prediction for leptogenesis is studied in a minimal lepton model. Finally, we draw the conclusions in section 5. The finite modular group Γ 4 ∼ = S 4 and the compact expression of CG coefficients are listed in Appendix A. The concrete forms of modular multiplets at weight 4, 6, 8 are presented in Appendix B. We give the general modular invariant superpotentials and mass matrices for two right-handed neutrino case in Appendix C.
2 Modular symmetry and modular forms of level N = 4 In the modular invariant framework with single modulus, the modular symmetry is described by the modular group which is the special linear group SL(2, Z) of degree two over integers: SL(2, Z) is often denoted Γ as well. It have two generators S and T with which satisfy the relations The modular group Γ has an important class of normal subgroups called principal congruence subgroup of level N which is defined as Note that Γ(1) = SL(2, Z) and T N ∈ Γ(N ). We can obtain the finite modular group from the quotient group: whereΓ andΓ(N ) are the projective groupsΓ = Γ/{±1} andΓ(N ) = Γ(N )/{±1}. Notice that Γ(N ) ∼ =Γ(N ) for N > 2. The group Γ N is usually called inhomogeneous finite modular group of level N , and similarly the homogeneous finite modular group Γ N can be defined as Γ N ≡ Γ/Γ(N ) which can also be generated by S and T obeying the multiplication rules S 4 = (ST ) 3 = T N = 1 [44], and additional relations are need to render the group finite for N ≥ 6 [64]. The group Γ N is isomorphic to quotient group of Γ N over its center {1, S 2 = −1}, consequently Γ N is the double covering group of Γ N , and Γ N has twice as many elements as Γ N . Top-down constructions in string theory generally lead to homogeneous finite modular group Γ N rather than the inhomogeneous finite modular group Γ N [61,62]. The modular group SL(2, Z) acts on the upper half plane H = {τ ∈ C | Im(τ ) > 0} by the linear fractional transformation: Obviously γ and −γ give the same action on the modulus τ , thus each linear fractional transformation corresponds to an element of the projective special linear groupΓ. If all the points in the orbit of a modulus τ is identified, we obtain the coset space H/Γ which is the so-called fundamental domain D of SL(2, Z): which is a hyperbolic triangle bounded by the vertical lines Re(τ ) = 1 2 , Re(τ ) = − 1 2 and the circle |τ | = 1. Every point τ ∈ H is equivalent to a point of D via the action of SL(2, Z), and no two distinct points inside D are equivalent under the action of SL(2, Z) and two points of D are in the same orbit only if they lie on the boundary of D.
The modular form of integral weight k and level N is a holomorphic function of τ , and it transforms under Γ(N ) as follows Therefore the weight k "differential form" f (τ )(dτ ) k/2 is invariant under the action of every element of Γ(N ). The modular forms of weight k and level N span a finite dimensional linear space M k (Γ(N )). The product of a modular form of weight k 1 with a modular form of weight k 2 is a modular form of weight k 1 + k 2 . Thus the set M(Γ(N )) = ∞ k=0 M k (Γ(N )) of all modular forms of level N forms a graded ring. Furthermore, it has been proved that the finite dimensional space M 2k (Γ(N )) can be decomposed into irreducible representations of the finite modular groups Γ N [2,44] up to the automorphy factor (cτ + d) 2k . That is to say, it is always possible to choose a basis in M 2k (Γ(N )) so that Y (2k) r = (f 1 (τ ), f 2 (τ ), . . . ) T transform under the full modular group Γ as where ρ r (γ) is the irreducible representation of quotient group Γ N .
In the present work, we are interested in the level N = 4, the linear space of the modular forms of level 4 is well established, and it can be constructed by making use of Dedekind eta function or the theta constants [21,46,55]: where Dedekind eta function η(τ ) is defined by and the theta constants are defined as Thus the dimension of the modular space M 2k (Γ(N )) is equal to 4k + 1. In the working basis given in Appendix A, all the modular forms of weight k and level 4 can be expressed as the homogeneous polynomials of degree 2k in the modular functions ϑ 1 and ϑ 2 which are linear combinations of θ 2 (τ ) and θ 3 (τ ) as follows [55], with ω = e 2πi/3 . In particular, the weight 2 modular multiplets Y can be written as [55] The weight 2 modular forms of level 4 can also be constructed from the derivative of eta function [33,34] or the products of the Dedekind eta function [21], the resulting q−expansions of the modular forms would be identical when going to the same representation basis of the modular generators S and T . The higher weight modular forms can be generated by the tensor product of the weight 2 modular forms, and their specific forms can be found in Appendix B. We summarize the modular multiplets of level 4 up to weight 8 in table 1.
3 Fermion mass models based on S 4 modular symmetry with generalized CP We shall briefly review the modular invariance approach in the following, then recapitulate on the consistency condition which should be fulfilled to consistently combine modular symmetry with generalized CP symmetry (gCP). Furthermore, we perform a systematical classification of quark and lepton mass models based on the modular symmetry Γ 4 ∼ = S 4 and gCP. Table 1: Summary of the even weight modular forms at level N = 4, the subscript r denotes the irreducible representations of the inhomogeneous finite modular group Γ 4 ∼ = S 4 . Here Y (6) 3I and Y (6) 3II denote the two linearly independent weight 8 modular forms in the doublet representation 2, and we adopt a similar notation for Y

The framework
We formulate our models in the framework of modular invariant approach with N = 1 global supersymmetry [2]. The field content consists of a set of chiral matter superfields Φ I and a modulus superfield τ , their modular transforms under SL(2, Z) are given by: where −k I is called the modular weight of the matter field Φ I , and ρ I (γ) is the unitary representation of Γ N . The Kähler potential is taken to be the minimal form following the convention of [2] K which gives rise to the kinetic terms of the matter fields and the modulus field after the modular symmetry breaking caused by vacuum expectation value (VEV) of τ . Notice that the modular invariance doesn't fix the Kähler potential in the bottom-up approach [65], and the Kähler potential could receive unsuppressed contributions from modular forms. However, generally both traditional flavor symmetry and modular symmetry are present in top-down approach such as the string derived standard-like models [66][67][68], the off-diagonal contributions to the Kähler metric are forbidden by the traditional flavor group and the minimal Kähler potential in Eq. (16) appears as the leading order term. Even including the whole modular dependence in the Kähler potential in these top-down models, the resulting phenomenological predictions don't differ from those which have been obtained by using just the standard Kähler potential Eq. (16). The superpotential W(Φ I , τ ) can be expanded in power series of the involved supermultiplets Φ I , where Y I 1 ...In is a modular multiplet of weight k Y as introduced in previous section. Modular invariance requires that each term of the W(Φ I , τ ) satisfies the following conditions: In order to improve the prediction power of the modular invariance approach, we include the generalized CP symmetry further. It is known that the complex modulus τ transforms under the action of gCP as [60,61,[69][70][71] up to modular transformations. The gCP transforms a generic chiral superfield Φ into the hermitian conjugate superfield with x = (t, x) and x P = (t, − x), the CP transformation matrix X r is a unitary matrix acting on flavor space. The chiral superfield Φ(x) is assigned to an irreducible unitary representation ρ r of the finite modular group, then the form of the matrix X r is strongly constrained due to the presence of modular symmetry. Firstly applying a gCP transformation followed by a modular transformation and subsequently an inverse CP transformation, the complex modulus τ and the matter field Φ transform as follows where −k denotes the modular weight of Φ. The closure of the modular transformations and gCP transformations entails that the following consistency condition has to be satisfied [60,63], where u(γ) is an outer automorphism of the modular group, Here χ(γ) is called character and it is a homomorphism of SL(2, Z) into {+1, −1}. From the relations S 4 = (ST ) 3 = 1 satisfied by the modular generators S and T , it is easy to know that only two possible values of the character are allowed [72] χ(S) = χ(T ) = 1 , or χ(S) = χ(T ) = −1 .
Consequently two possible gCP symmetries can be defined in the context of modular invariance. 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 the consistency condition of Eq. (22) on the generators S and T . For the first kind of gCP associated with the trivial character χ(S) = χ(T ) = 1, the consistency condition becomes The second kind of gCP associated with the nontrivial character χ(S) = χ(T ) = −1 has been studied in [46,63], the explicit form of X r depends on both the modular weight −k and representation assignment r of the matter field, and obviously it would be reduced the first gCP for (−1) −k ρ r (S 2 ) = 1. Notice that ρ r (S 2 ) = ±1 because of S 4 = 1. As a result, it is only relevant for the case of (−1) −k ρ r (S 2 ) = −1 which implies odd k for the inhomogeneous finite modular group Γ N , the gCP transformation X r is determined by [46,63] X r ρ * r (S)X −1 r = −ρ r (S −1 ), which can be satisfied if and only if the level N is even, the dimension of the representation ρ r is even together with the vanishing trace of ρ r (S) and ρ r (T ). If one intends to impose the second gCP in a model, the three generations of matter fields should be assigned to the direct sum of one-dimensional and two-dimensional representations of the finite modular group Γ N or Γ N , and second gCP acts nontrivially on the two matter fields in the doublet representation while the gCP transformation of the other matter field in the singlet representation can only be the first one. Moreover, the minus sign in Eq. (26) implies that the fermion mass matrix would be block diagonal and consequently some mixing angles would be constrained to be vanishing if the second gCP is implemented. Because none of quark or lepton mixing angles are vanishing in spite of some very small quark mixing angles, we shall not consider the second gCP symmetry and focus on the first gCP in the present work. As regards the CP transformation of the modular forms, it has been shown that the integral weight modular forms are in the irreducible representations of Γ N fulfilling (−1) −k ρ r (S 2 ) = 1 [44] so that only the first gCP acts on the modular forms and they transform in the same way as the matter fields under gCP, i.e., Y r (τ ) if the basis of the modular space is properly chosen [63].
The explicit form of the gCP transformation X r is determined by the consistency condition in Eq. (22) up to an overall phase for any given irreducible representation r. For the concerned first gCP and the finite modular group Γ 4 ∼ = S 4 with the basis listed in table 9, solving the consistency conditions of Eq. (25), we find that the generalized CP transformation X r is in common with the representation matrix of S, which is a combination of the modular symmetry transformation S and the canonical CP transformation. Modular invariance requires that the action is invariant under the modular transformation S, thus the gCP transformation in Eq. (27) is essentially the canonical CP transformation. Furthermore, for the level 4 modular forms built from Y 2 (τ ) and Y 3 (τ ) up to weight 8, it is straightforward to check that they transform under gCP as follows which is consistent with the general results of [63]. As given in Appendix A, all the Clebsch-Gordan (CG) coefficients in our working basis are real, thus the gCP symmetry would constrain all the coupling constants to be real. In the modular invariant theory with gCP symmetry, both modular and CP symmetries are uniquely broken by the VEV of the modulus τ . In particular, all CP violation phases arises from non-vanishing real part of τ . In the following, we shall perform a systematical classification of the Yukawa superpotential according to the transformation properties of the matter fields under the Γ 4 ∼ = S 4 modular symmetry. We assume the Higgs doublets H u and H d are S 4 trivial singlet 1 and their modular weights k Hu,H d are vanishing. Notice that k Hu,H d can always be set to zero by redefining the modular weights of matter fields.

Classifying the Yukawa couplings
The modular invariance approach is formulated in the framework of supersymmetry, and we adopt the gauge symmetry, lepton superfields, quark superfields and Higgs multiplets of the minimal supersymmetric standard model. We consider both scenarios that neutrinos are Dirac or Majorana particles, and the neutrino masses are generated by the type I seesaw mechanisms if neutrinos are Majorana particles. It is known that at least two right-handed neutrinos are necessary to generate the non-vanishing solar and atmospheric neutrino mass squared differences. We have considered both cases with two and three right-handed neutrinos. For simplicity, we denote the left-handed lepton and quark doublets as F and the right-handed lepton and quark singlets as F c , i.e. F c ∈ {u c , d c , E c , N c } and F ∈ {Q, L}. The three generations of matter fields can be assigned to transform as a triplet F (c) ∼ 3 i under S 4 modular group, the direct sum of doublet and singlet F (c) ∼ 2 ⊕ 1 i , or the direct sum of the three singlets F (c) ∼ 1 i 1 ⊕ 1 i 2 ⊕ 1 i 3 . If only two right-handed neutrinos are introduced, they can transform as a doublet or two singlets under S 4 , as discussed in Appendix C. Therefore, there will be many possible S 4 modular invariant models for quarks and leptons. The main purpose of this paper is to classify all of these possible fermion mass superpotentials. In the following, we will consider modular forms of weight less than 10, and the analytical results reached can be easily extended to much higher weight modular forms analogously.
Before going into the concrete discussion below, let us make an explanation of the notations. We denote the S 4 singlet and triplet representations as 1 ≡ 1 0 , 1 ≡ 1 1 , 3 ≡ 3 0 , 3 ≡ 3 1 . We use i, j, k, l to represent the indices of the singlet or the triplet representations, and they can only take the values 0 or 1, i.e., i, j, k, l ∈ {0, 1}. The lowercase letter a, b are used to label the component of the modular multiplets, and they can only take the value 1, 2 and 3, i.e. a, b ∈ {1, 2, 3}. For simplicity of the formula, we introduce the superfluous notations Y 2,3 which are set to zero. Moreover, we use the capital letters A, B, C to describe the degeneracy of the modular multiplets. For instance, there are two weight 6 modular forms Y (6) 3I and Y (6) 3II in the triplet representation 3. Furthermore, we introduce the operations <> and ≺ and they are defined as < i >= i (mod 2) and ≺ i = i (mod 3) which take values in the range of {0, 1} and {1, 2, 3} respectively. Notice that we define ≺ i = 3 if i is divisible by 3. It is remarkable that the general analytical expression of the fermion mass matrix can be read out for each possible representation assignment of the matter fields.
In this section, we will investigate the Yukawa superpotential for the fermion masses, which can be generally written as where all independent S 4 contractions should be considered and different singlet combinations are associated with different coefficients. The function f (Y ) is some modular form multiplet fixed by the weight and representation assignments of the matter fields F and F c . In the following, we give the concrete form of the Yukawa superpotential and the corresponding fermion mass matrix for different S 4 transformation properties of F and F c .
Let us first consider the case that both left-handed and right-handed fermions transform as triplets under S 4 . The modular weights of F c and F are denoted as k F c and k F respectively. The general Yukawa supepotential for this assignment is given by Here we have assumed that the modular form multiplets Y in all S 4 irreducible representations are present. As shown in table 1, certain modular multiplets at some specific modular weights are not allowed and thus the corresponding terms should be dropped. The fermion mass matrix can be read out from this superpotential: where repeated indices are implicitly summed over. If F and F c are quark and charged lepton fields, the case of k F c + k F = 0 is not viable, since it gives rise to three degenerate mass eigenvalues. On the other hand, if W F describe the neutrino Dirac coupling under the assumption of Majorana neutrinos, the vanishing modular weight k F c + k F = 0 is allowed.
In this case, the three generations of left-handed fermions F transform as a triplet of S 4 and the right-handed fields F c are assigned to be singlets of S 4 . The modular weight of F and F c are denoted by k F and k F c 1,2,3 respectively. Notice that permutating the assignments of the three right-handed fermions F c amount to multiplying certain permutation matrix from the right side of the mass matrix, consequently the results for the charged fermion masses and mixing matrix are left invariant. The superpotential for this assignment can be written as which leads to the following fermion mass matrix If two right-handed fields are assigned to have the same modular weight and representation assignment and they couple with a unique modular multiplet, two rows of the mass matrix would be proportional such that one mass eigenvalue would be vanishing.
In some specific cases, the mass matrix can also gives zero eigenvalue. For example, from Appendix B we can see the modular forms Y and Y (4) 3 respectively. As a consequence, the fermion mass matrix for the assignment (2,4,8) and (< i 1 + j >, < i 2 + j >, < i 3 + j >) = (0, 0, 0) will have zero mass eigenvalue as well.
Without loss of generality, we assign the first two right-handed fermions F c D = (F c 1 , F c 2 ) T to transform as a doublet under S 4 and the third one F c 3 is the singlet. The modular weights of these fields are k F c D , k F c 3 and k F . The general superpotential is of the following form: The mass matrix which can be read out from this superpotential is In some cases, the rank of M F is less than three due to the structure of the modular forms. For instance, the rank of the mass matrix is two for the assignment (k F c D +k F , k F c 3 +k F ) = (2, 4).
We interchange the representation assignments of the left-handed and the right-handed fields discussed in above. The left-handed fermions are assigned to the direct sum of a doublet F D = (F 1 , F 2 ) ∼ 2 and a singlet F 3 ∼ 1 j , while the right-handed fermions transform as a triplet under S 4 . The modular weights of these fields are denoted as k F c , k F D and k F 3 . Then we can straightforwardly read out the Yukawa superpotential for this kind of assignment, The resulting mass matrix is given by which is the transpose of the mass matrix in Eq. (35) with the indices i and j exchanged.
In this case, both left-handed and right-handed fields are assigned to the direct sum of S 4 doublet and singlet. We denote Then the Yukawa superpotential is given by which leads to The above mass matrix can be divided into four parts Let us firstly consider the (33) entry M 33 which involves the modular forms in the singlet representations of S 4 . From table 1, we see that there are only four singlet modular forms Y 1 , Y Notice that odd weight k F 3 + k F c 3 = 1, 3, 5, . . . can also lead to vanishing M 33 , but the rank of M F would be less than three so that at least one mass eigenvalue is zero. Then we proceed to consider the M 3D block consisted of the (31) and (32) . ., some mixing angles or masses are vanishing 1 . As regards the M D3 block consisted of the (13) and (23) entries, it would be vanishing if the modular weight k F c D + k F 3 is non-positive or odd. However, odd k F c D + k F 3 leads to vanishing fermion masses or mixing angles. Although either M 3D or M D3 can be vanishing, they can not be vanishing simultaneously otherwise some masses or mixing angles are constrained to be zero.
is odd in this case. As a result, both up and down quark (charged lepton and neutrino) mass matrices are block diagonal simultaneously such that some mixing angles are vanishing.
Analogous to previous cases, we find the Yukawa superpotential takes the following form, The fermion mass matrix is determined to be In above, we don't consider the singlet assignment for the left-handed fields F , because generally more free coupling constants would be necessary in the resulting quark and lepton models, and we are mainly concerned with the models with small number of free parameters in the present work.

Classifying the Majorana mass terms
In this subsection, we explore the superpotential for the Majorana mass terms, which can be written as where Λ is the characteristic scale of flavor dynamics, and all independent invariant singlets should be included. The function f (Y ) refers to the modular multiplets to ensure modular invariance, and it is fixed by the modular weight and representation of F c .
As shown in the Appendix A, the contraction 3 i × 3 i → 3 is antisymmetric. Thus the triplet modular forms transforming as 3 don't contribute to the Majorana mass terms of F c . The superpotential W F c reads as The Majorana mass matrix of F c is symmetric The contraction 2 × 2 → 1 is antisymmetric and consequently it has no contribution to mass matrix. The general superpotential for the Majorana mass is which gives rise to the following mass matrix, In the same fashion, we can read out the most general Majorana mass terms for the singlet assignment of F c , and the mass matrix is If the three generations of F c all transform as singlets under S 4 , the Lagrangian would be less constrained by modular symmetry and consequently more free parameters would be introduced in the Yukawa coupling and the Majorana mass term.

Equivalence of different assignments
The possible quark models with S 4 modular symmetry can be obtained by combining the possible forms of the up quarks and down quarks Yukawa couplings discussed in previous section 3. 3. Similarly the possible lepton models can be obtained for Dirac neutrinos, and the Majorana mass terms of the right-handed neutrinos should be considered as well if neutrinos are Majorana particles. In the present work, the left-handed quark and lepton fields are assumed to transform as a triplet or the direct sum of doublet and singlet under S 4 , all the three possible assignments: triplet, doublet plus singlet and three singlets for the righthanded quark and lepton fields would be considered. It is notable that different assignments can lead to the same predictions for fermion masses and mixing matrix. For instance, if both left-handed leptons F and right-handed charged leptons F c transform as S 4 triplets F c ∼ 3 i and F ∼ 3 j , the mass matrix can be read from Eq. (31) for any given modular weights. It Table 2: Transformation of the fermion mass matrix under changing the representation assignments of matter fields, and P is the diagonal matrix diag{1, −1, 1}. In the case that both F and F c are assigned to direct sum of doublet and singlet under S 4 , the couplings α l 1 associated with the operators (F c D F D ) 1 l should be transformed into −α l 1 , and the coupling α 1 associated with the operator (F c D F c D ) 1 should also change to −α 1 .
can be easily seen that the representation assignment F ∼ 3 <i+1> , F c ∼ 3 <j+1> gives the same charged lepton mass matrix.
Two different kinds of representation assignments can also give mass matrices related by phase transformations, as summarized in table 2. As an example, let us consider the case a ∼ 1 ja with a = 1, 2, 3, the general form of the charged lepton mass matrix is given by Eq. (43). If we change the representation assignment the charged lepton mass matrix would turn into Analogously changing the representation of the right-handed neutrinos, the light neutrino mass matrix would change as Dirac neutrinos : The phase matrix diag{1, −1, 1} can be absorbed into the lepton fields, the lepton masses and mixing parameters are left invariant. As a consequence, without loss of generality, we can take F ∼ 3 for the triplet assignment of the left-handed fields and F ∼ 2 ⊕ 1 for the doublet plus singlet assignment.
Since the signal neutrinoless double beta decay has not been observed, the nature of neutrinos is still unknown. We shall consider both Majorana and Dirac neutrinos in this work. The light neutrino masses are generated by the type-I seesaw mechanism for Majorana neutrinos, and the light neutrino mass matrix is given by the seesaw formula where M D and M N c are Dirac mass matrix and the Majorana mass matrix of the right-handed neutrinos respectively. For Dirac neutrinos, additional symmetry is generally necessary to forbid the right-handed neutrino Majorana mass term and it is usually taken to be U (1) L lepton number. In the context of modular invariance approach, the Majorana mass terms of the right-handed neutrino can be naturally forbidden by taking the modular weights of right-handed neutrinos N c to be negative integers, because there are no modular forms of negative weight.

Phenomenologically viable models and numerical results
From the general analytical expressions of the mass matrix for different representation of matter fields, we can straightforwardly obtain the possible lepton and quark models based on S 4 modular symmetry. In this work, we are interested in the models with small number of free parameters: lepton models with less than 9 free parameters and quark models with less than 11 free parameters. For each model, we perform a conventional χ 2 analysis and we use the well-known package TMinuit to numerically search for minimum of the χ 2 function and determine the best values of the input parameters. Then we evaluate masses and mixing parameters of quarks and leptons at the best fit points, and determine whether they are within the experimentally allowed 3σ regions. The overall scale factor of the mass matrix can be adjusted to reproduce any one of the mass eigenvalues. For instance, the overall factors of the charged lepton, up type quark and down type quark mass matrices are fixed by the measured values of the electron, top quark and down quark masses respectively. The overall scale of the neutrino mass matrix is determined by the solar neutrino mass square difference ∆m 2 21 . We scan over the parameter space of the models, the ratios of coupling coefficients are taken as random numbers whose absolute values freely vary in the range of [0, 10 5 ]. Moreover, the vacuum expectation value (VEV) of the complex modulus τ is also treated as a free parameter to optimize the agreement between predictions and experimental data. Since each point of τ in the complex upper-half plane can be mapped into the fundamental domain D given in Eq. (7) by a modular transformation, thus it is sufficient to limit the modulus VEV τ in the fundamental domain D. Under the CP transformation τ → −τ * , gCP invariance implies that the fermion mass matrix becomes [63]. Therefore at the CP dual point τ → −τ * , the predictions for fermion masses and mixing angles are left unchanged while the signs of all CP violation phases are flipped.
We use the fermion mass ratios, mixing angles and CP violation phases to construct the χ 2 function, the experimental data of the leptons and quarks are summarized in table 3. The data of the lepton mixing parameters are taken from the latest global fit of NuFIT v5.0 including the atmospheric neutrino data from Super-Kamiokande [73] and the neutrino mass spectrum is taken to be normal ordering for illustration. The charged lepton mass ratios are taken from [74], and the quark mixing parameters and mass ratios are adopted from [75], and they are calculated at the GUT scale M GUT = 2×10 16 GeV in a minimal SUSY breaking scenario, with SUSY breaking scale M SUSY = 1 TeV and tan β = 7.5,η b = 0.09375. The leptonic Dirac CP phase δ l CP has not been accurately measured, therefore we don't include the contribution of δ l CP in the χ 2 function. If all observables at the best fit point of a model are compatible with the experimental data at 3σ level, this model would be regarded as phenomenologically viable. In the following, we report the fitting results of the viable models with the minimal number of free parameters, and all numerical results are shown with six significant digits. Notice that Re τ and all CP violation phases flipped their signs while all other observables and free parameters are unchanged at the CP dual point.

Lepton models
For Majorana neutrinos, the minimal phenomenologically viable models only depend on five free parameters besides the complex modulus τ , and we find nine such models labelled as L1∼L9. Notice that the model L2 was firstly presented in [60]. The S 4 representation and modular weights of the lepton fields in each models are listed in  0.03888 ± 0.00062 Table 3: The central values and the 1σ errors of the mass ratios and mixing angles and CP violation phases in lepton and quark sectors. The charged lepton mass ratios averaged over tan β are taken from [2,74], and we adopt the values of the lepton mixing parameters from NuFIT v5.0 with Super-Kamiokanda atmospheric data for normal ordering [76]. The data of quark mass ratios and mixing parameters are taken from [75] with the SUSY breaking scale M SUSY = 1 TeV and tan β = 7.5,η b = 0.09375.
values of the coupling constants and the corresponding predictions for the lepton masses and mixing parameters are summarized in table 5. Although we have considered the minimal seesaw model with two right-handed neutrinos, three right-handed neutrinos are involved in these minimal models. It turns out that more free parameters are needed to accommodate the experimental data in the modular models with two right-handed neutrinos. In most modular symmetry models, both left-handed leptons L and right-handed neutrinos N c are assumed to transform as a triplet under the finite modular group while the right-handed charged leptons E c are singlets. Our models L2, L3, L4 and L5 belong to this category. It is notable that we find new possible assignments here. All the lepton fields L, E c and N c are S 4 triplet in the model L1. Both L and N c transform as triplet 3 or 3 under S 4 while the right-handed charged leptons are in the reducible representation 2 ⊕ 1 in the models L6 and L7. Furthermore, both L and E c are assigned to the direct sum of the doublet and singlet of S 4 in the models L8 and L9. From table 5, we can see that all these models can accommodate the experimental data very well, the atmospheric mixing angle θ 23 is predicted to be in the second octant. The Dirac CP violation phase δ l CP is determined to be sizable in these models, and it distributes in the range of [1.27π, 1.65π]. The upcoming generation of long-baseline neutrino oscillation experiments such as DUNE [77][78][79][80] and Hyper-Kamiokande [81] can significantly improve the sensitivity to θ 23 and δ l CP . It is expected that a 5σ discovery of CP violation can be reached after ten years of data taking over 50% of the parameter space. Thus our predictions for θ 23 and δ l CP can be tested in near future. The neutrino mass scale can be probed from direct kinematic searches, neutrinoless double decay and cosmology. The cosmological observation is sensitive to the sum of light neutrino masses m i , and the most stringent bound is m i < 0.12 eV at 95% confidence level from the Planck Collaboration [82]. All the minimal models satisfy this bound except L8 and L9 which give m i 121 meV very close to the upper limit. Notice that the cosmological bound on the neutrino masses significantly depend on the data sets that need to be combined in order to break the degeneracies of the many cosmological parameters [82]. Combining the Planck lensing with the baryon acoustic oscillation (BAO) data and the acoustic scale measured by the CMB, the neutrino mass is constrained to be i m i < 600 meV [82]. The limits also become weaker when one departs from the framework of ΛCDM plus neutrino mass to frameworks with more cosmological parameters. The direct kinematic searches provide the most model independent approach to test the neutrino mass, and the neutrino mass extracted from ordinary beta decay is m β = i |U ei | 2 m 2 i = cos 2 θ l 12 cos 2 θ l 13 m 2 1 + sin 2 θ l 12 cos 2 θ l 13 m 2 2 + sin 2 θ l 13 m 2 3 , where U is the lepton mixing matrix. From the values of lepton mixing angles and neutrino masses, we can determine the effective mass m β , as shown in table 5. The predictions for m β are around 15 meV, and they are much below the upper limit m β < 1.1 eV given by KATRIN [83]. It is expected that KATRIN can advance the sensitivity on m β by one order of magnitude down to 0.2 eV after 5 years, and the next generation experiments such as Project 8 may be able to reach the 50 meV level [84]. Therefore a positive signal of KATRIN or Project 8 in near future could rule out our models. It is known that the neutrinoless double beta (0νββ) decays of even-even nuclei are important to test the Majorana nature of neutrinos, they can provide valuable information on the neutrino mass spectrum and the CP violation phases. The amplitude of the neutrinoless double beta decay is proportional to the effective Majorana mass m ββ which is given by ei m i = cos 2 θ l 12 cos 2 θ l 13 m 1 + sin 2 θ l 12 cos 2 θ l 13 e iα 21 m 2 + sin 2 θ l 13 e i(α 31 −2δ l CP ) m 3 . (55) The strongest bound on m ββ is set by the KamLAND-Zen experiment m ββ < (61 − 165) meV [85], where the largest uncertainty arises from the computation of the associated nuclear matrix element. There are many 0νββ decay experiments in the plan and construction, which aim to improve the current bounds on m ββ . The future large scale 0νββ decay experiments have the potential of measuring the decay half-life exceeding 10 28 years. For instance, the SNO+ Phases II is expected to reach a sensitivity of 19 − 46 meV [86]. The LEGEND experiment intends to achieve a sensitivity of 15-50 meV by operating 1000 kg of detectors for 10 years [87]. The nEXO is the successor of EXO-200, and its projected m ββ sensitivity is 5.7 − 17.7 meV after 10 years of data taking [88]. Using the master formula of Eq. (55), we can determine the values of the effective Majorana neutrino mass m ββ at the best fitting points, as given in table 5. We see that the latest bound of KamLAND-Zen experiment is well satisfied and the predictions are within the reach of future tonne-scale 0νββ experiments except the models L5 and L6 which are experimentally very challenging because of the quite low values of m ββ . It is remark that these minimal viable models only use five real parameters together with the complex modulus τ to describe 12 observables: 3 charged lepton masses, 3 neutrino masses and 3 lepton mixing angles and 3 CP violating phases. Thus the values of the free parameters are strongly constrained by the experimental data and the different observables should be correlated with each other. For example, the light neutrino mass matrix only depends on the modulus τ up to an overall scale in the model L1 while there are four real couplings in the charged lepton superpotential with W e = α e (E c L) 1 Y 3 ) 1 . It is notable that the hierarchical masses of charged leptons can be reproduced although the four coupling constants α e , β e , γ e and γ e are of the same order of magnitude, as can be seen from table 5. Thus the charged lepton masses are also dictated by modular symmetry, and the hierarchical mass eigenvalues arise from the departure of τ from the self-dual fixed point τ = i [30,32,59]. Furthermore, we take the models L1 and L6 as examples, and we comprehensively scan the parameter space of these two models. The lepton masses and mixing angles are required to lie in the experimentally  preferred 3σ regions [76], we display the correlations among the free parameters and observables in figure 1 and figure 2. It is worth mentioning that the experimental data can only be accommodated in small regions of parameter space such that the predictions for the lepton mixing parameters are quite precise and their allowed regions are small as well.
Since the signal of 0νββ decay has not been observed, the possibility that neutrinos are Dirac particles can not be excluded at present. Generally additional symmetry such as lepton number conservation is necessary to forbid the Majorana mass term of the righthanded neutrinos. Modular invariance can naturally enforce Dirac neutrinos if the modular weights of the right-handed neutrinos are negative. In the same fashion, we can analyze the possible Dirac neutrino mass models with S 4 modular symmetry and gCP symmetry. We find that the phenomenologically viable models make use of at least five couplings besides the modulus τ , and four minimal models are found. The modular transformation properties of the lepton fields and the results of the χ 2 analysis are reported in table 6. Notice that the model D3 was already discussed in [40]. These models can be tested by the measurement of θ l 23 and δ l CP at future long baseline neutrino oscillation experiments [77][78][79][80][81]. The effective neutrino mass m β in beta decay is predicted to be an order of magnitude below the expected sensitivity of the KATRIN experiment [83].

Leptogenesis
Since we impose gCP as symmetry on the model, all couplings in the superpotential are constrained to be real in our working basis. As a consequence, all CP violations uniquely arise from the modulus vacuum expectation value. Thus the CP violation in leptogenesis is naturally related to the CP violation phases δ l CP , α 21 and α 31 in the lepton mixing matrix. Early studies of leptogenesis in the context of modular symmetry models without gCP symmetry can be found [20,39,89]. During the final preparations of this paper, a preprint discussing leptogenesis in a A 4 modular model with gCP appeared on the arXiv [90]. In this section, we shall study whether the measured value of the baryon asymmetry of the universe Y B0 = (0.870300 ± 0.011288) × 10 −10 [82] can be correctly generated through leptogenesis in the minimal S 4 modular invariant models found in previous section, where the subscript 0 implies "at present time".
The right-handed neutrino masses depend on the overall scale Λ in our model. In the present work, we assume that the right-handed neutrinos are heavy with masses above 10 12 GeV, thus we work in the framework of unflavored thermal leptogenesis. The modular invariance is formulated in the supersymmetric context, as shown in section 3.1, therefore we should consider supersymmetric leptogenesis. The out-of-equilibrium decays of the lightest right-handed neutrinos or sneutrinos in the early universe produce lepton asymmetries. We    Table 6: Summary of the lepton models based on S 4 modular symmetry and gCP symmetry, where neutrinos are Dirac particles. The integer k should be greater than 2 so that the modular weight of N c is negative and modular invariance forbids the Majorana mass term of right-handed neutrinos. Notice that the Higgs fields are invariant under S 4 with zero modular weight. The best fit values of the input parameters are also included. We give the predictions for neutrino mixing angles, and Dirac CP violating phase, and the neutrino masses, and the effective neutrino masses m β probed by direct kinematic search in beta decay. Note that the transformation τ → −τ * leaves all observables unchanged except shifting the sign of the CP phase δ l CP . Figure 1: The correlations among the input parameters, lepton mixing angles, CP violation phases and neutrino masses in the model L1. The lepton masses and mixing angles are required to lie in the experimentally preferred 3σ regions [76]. Notice that the transformation τ → −τ * leaves all observables unchanged except shifting the signs of the CP phases δ l CP , α 21 and α 31 , consequently we don't show the CP conjugate region. In the last panel, the red (blue) dashed lines indicate the most general allowed regions for inverted ordering (normal ordering) neutrino mass spectrum respectively, where the neutrino oscillation parameter are varied over their 3σ ranges. The present upper limit m ββ < (61 − 165) meV from KamLAND-Zen [85] is shown by horizontal grey band. The vertical grey exclusion band denotes the current bound coming from the cosmological data of i m i < 0.120eV at 95% confidence level obtained by the Planck collaboration [82].
denote the decay asymmetries for the decay of heavy neutrino into Higgs and lepton, neutrino into Higgsino and slepton, sneutrino into Higgsino and lepton, and sneutrino into Higgs and slepton as ε 1 , ε 1˜ , ε1 and ε1 respectively which are defined by [91,92] In the minimal supersymmetric standard model, all the above decay asymmetries are equal ε 1 = ε 1 = ε 1 = ε 1 [91,92]. In the basis where the Majorana mass matrix of the righthanded neutrinos is diagonal and real, the lepton asymmetry parameter ε 1 is given by [91,93]  where λ ν is the neutrino Yukawa coupling matrix in the convention (λ ν ) ij N c i (L j · H u ) and the loop function g reads Non-vanishing asymmetry parameter ε 1 requires that the off-diagonal entries of the product λ ν λ † ν are complex and different from zero. For the models L1, L4, L5, L6, the product λ ν λ † ν is proportional to the unity matrix. Consequently the lepton asymmetry ε 1 is vanishing at LO and a net baryon asymmetry can not be generated. The masses of the three right-handed neutrinos are degenerate in the models L2, L3 and L7, the baryon asymmetry is generated in the regime of resonant leptogenesis. Hence we take the model L9 as an example in the following. The lepton number asymmetry is partially converted into a non-zero baryon number asymmetry by the fast sphaleron interactions in the thermal bath in the Early Universe. For all temperature ranges, the produced baryon asymmetry normalized to the entropy density can be computed from the B − L asymmetryŶ ∆ as follows The B−L asymmetryŶ ∆ can be computed by solving the following Boltzmann equations [92,94] dY where z = M 1 /T with T being the temperature. K 1 (z) and K 2 (z) are the modified Bessel functions of the second kind. Y N 1 and Y N 1 denotes the density of the lightest right-handed neutrino N 1 with mass M 1 and its supersymmetric partners N 1 . The notations Y eq N 1 , Y eq N 1 and Y eq are corresponding equilibrium number densities and they take the following form with g * = 228.75 being the number of degrees of freedom. Moreover, the washout parameters K α and K are defined as where v u = v sin β, m * sin 2 β × 1.58 × 10 −3 eV .
At the best fit point of model L9, we find the value of K = 1.293 × 10 6 1 which implies a strong washout. In the strong washout regime, the functions f 1 (z) and f 2 (z) can be approximated as [95] where M h = 125 GeV is the mass of Higgs boson. The free parameters are fixed at their best fit values shown in table 5, notice that only the combination (α D v u ) 2 /(β N Λ) can be determined by the data of lepton masses and mixing. Numerically solving the Boltzmann equations, we find that the observed baryon asymmetry can be produced for the following values of the flavor scale: Accordingly the right-handed neutrino masses are determined to be M 1 1.985 × 10 15 GeV, M 2 6.723 × 10 15 GeV and M 3 6.834 × 10 15 GeV. The VEV of τ is the unique source of modular symmetry and gCP symmetry breaking in this model. All CP violations at both low energy and high energy should significantly depend on τ . We plot the contour region of Y B in the plane Im τ versus Re τ , where we fix all the coupling constants at their best fit values and Λ at the value in Eq. (67). The green area indicates the 3σ allowed region by the experimental data of lepton masses and mixing angles in the same plane. We see that there exists a small region of τ where both the flavor structure of the lepton and the baryon asymmetry of the universe can be explained. At the boundary of the fundamental domain D and the imaginary axis, certain residual gCP symmetry is preserved such that the lepton asymmetry ε 1 vanishes and no matter-antimatter asymmetry can be generated. Hence the VEV τ should deviate from the CP conserved points in order to obtain non-trivial CP violation in neutrino oscillation as well as a net baryon asymmetry. Furthermore, we plot the correlation between the Y B and the Dirac CP phase δ l CP in figure 3. Here the modulus vacuum expectation value τ is treated as a random complex number in the fundamental domain, the charged lepton masses and the neutrino mass squared differences and all three mixing angles are required to be within their experimentally preferred 3σ ranges [76]. Imposing the observed baryon asymmetry Y B , the allowed region of δ l CP would shrink considerably and it is around 1.28π.

Quark models
In the same fashion as we have done in the lepton sector, we can easily find out all possible quark models from the general results of section 3.2, and subsequently we perform a χ 2 analysis for each model to determine whether it can accommodate the precisely measured values of six quark masses m u , m c , m t , m d , m s , m b , and three quark mixing angles θ q 12 , θ q 13 , θ q 23 , and a CP violation phase δ q CP in the quark sector. We are concerned with the models which can reproduce the data with the smallest number of free parameters. It turns out that the minimal viable models labelled as Q1∼Q10 make use of seven real coupling constants in addition to the modulus τ , thus one prediction can be reached. The transformation of quark fields under S 4 and their modular weights are listed in table 7. It can be seen that the lefthanded quark fields Q are assigned to a triplet or doublet plus singlet of S 4 modular group, while the right-handed quark fields u c and d c are singlets or the direct sum of doublet and singlet. We present the best fit values of the input parameters and the predictions for quark masses and CKM mixing parameters in table 8. Here we have omitted these models with a high χ 2 min or a large number of parameters. It is known that the quark masses and mixing parameters have been precisely measured and their errors are quite small. The hierarchical patterns of quark masses and CKM matrix generally require more free parameters in a concrete model and it is quite difficult to explain the quark data with few (less than nine) parameters.

Toward unified description of quarks and leptons
The flavor structures of quarks and leptons are drastically different from each other, and it is not known at present whether the quark and lepton sectors are dictated by the same fundamental principle or not. In the previous two sections we have discussed individually the possible lepton and quark models with the smallest number of free parameters . In the following, we shall investigate whether quarks and leptons can be simultaneously described by the S 4 modular symmetry and gCP. In this scenario, both lepton and quark mass matrices would depend on a common complex modulus τ , and all the CP violation phases in lepton and quark sectors arise from the modulus VEV τ . The quark-lepton unification has been studied in the context of A 4 [11,14,28,31], T [45], S 4 [47] and A 5 [49] modular symmetries. The most predictive model contains fifteen parameters including the real and imaginary part of the modulus τ [47], as far as we know. The unification description of quark and lepton mixing can also be achieved in the paradigm of traditional flavor symmetry combined with gCP, the resulting lepton and quark mixing matrices can be predicted in terms of only four rotation angles if the flavor group and gCP are spontaneously broken down to Z 2 × CP by certain flavons in the neutrino, charged lepton, up quark and down quark sectors [96][97][98][99]. However, the fermion masses are not constrained in this approach, additional symmetries and fields are necessary to realize the required residual symmetry.
By comprehensively scanning the possible quark-lepton models based on S 4 modular symmetry and gCP, we find that the minimal models use fifteen independent parameters including Reτ and Imτ to explain twenty-two observables: six quark masses m u,c,t,d,s,b , three quark mixing angles θ q 12,13,23 , a quark CP violation phase δ q CP , three charged lepton masses m e,µ,τ , three neutrino masses m 1,2,3 , three lepton mixing angles θ l 12,13,23 and three leptonic CP violation phases δ l CP , α 21,31 . In the following, we will present a benchmark model which contains five real couplings in the lepton sector and eight free couplings in the quark sector. The right-handed neutrinos N c are S 4 triplet 3, while all the other lepton fields L, E c and quark fields Q, u c , d c are assigned to 2 ⊕ 1. The transformation properties of the matter fields are given by where k is an arbitrary integer. We see that the lepton sector is exactly the aforementioned lepton model L8. From the general results in sections 3.2 and 3.3, we can read out the fermion mass matrices as follows, The best fit values of the input parameters for this unified model is determined to be The masses and mixing parameters of leptons and quarks are predicted to be sin 2 θ l 12 = 0.339580 , sin 2 θ l 13 = 0.0215912 , sin 2 θ l 23 = 0.615854 , δ l CP /π = 1.33145 , α 21 /π = 1.29710 , α 31 which are compatible with experimental data at 3σ level. All the coupling constants as well as complex modulus τ are treated as random numbers, and the 3σ bounds of the mass ratios and mixing angles of both quarks and leptons are imposed. The values of τ compatible with experimental data are shown in figure 5, the light blue and red areas represent the regions favored by the quark and lepton data respectively. We see that there indeed exists a small overlap region of τ indicated by black in which the flavor structure of quarks and leptons can be accommodated simultaneously. Moreover, we display the correlations among neutrino masses and mixing parameters in figure 5. Since the common τ region of quark and lepton sectors is very small, the allowed values of all observables scatter in quite narrow ranges around their best fit values.

Conclusions and discussions
Modular invariance is an promising framework to address the flavor puzzle of the standard model. In recent years, much effort has gone into the study of lepton models based on inhomogeneous and homogeneous finite modular groups. In the present work, we have performed a systematical analysis of the possible lepton and quark models with S 4 modular symmetry. Aiming at the minimal and predictive models, we impose the gCP symmetry so that all coupling constants are constrained to be real in our working basis and the vacuum expectation of the modulus is the unique source of modular and CP symmetry breaking. In the known S 4 modular symmetry models [33-40, 53, 54], usually the three generations of left-handed lepton fields and right-handed charged lepton are assumed to transform as triplet and singlet under S 4 . Besides the singlet representations 1, 1 and triplet representations 3 and 3 , the S 4 group has a doublet irreducible representation 2. The presence of doublet representation not only introduces new features in the modular invariant lepton models, but also provides a new expedient to describe the quark sector. We give the most general analytical expressions of the modular invariant Yukawa superpotential of charged fermions and the Majorana mass terms of right-handed neutrinos. We have analyzed both scenarios that neutrinos are Majorana particles and Dirac particles. Under the assumption of Majorana neutrinos, the light neutrino masses are generated by the type I seesaw mechanism, and the conventional seesaw models with three right-handed neutrinos and the minimal seesaw models with two right-handed neutrinos are analyzed.
We have comprehensively searched for the S 4 modular invariant lepton and quark models with the lowest possible number of free parameters. After heavy numerical analysis, we find that the minimal lepton models make use of five real couplings together with the modulus τ to describe the charged lepton masses, neutrino masses, lepton mixing angles and CP violation phases. Thirteen minimal lepton models are obtained, including nine Majorana neutrino models and four Dirac neutrino models, the classification of the matter fields under modular symmetry is summarized in table 4 and table 6. Notice that the models L2 and D3 were already presented in Ref. [60] and Ref. [40] respectively while all others are new. The experimental data from neutrino oscillation, neutrinoless double decay, tritium beta decay and cosmology on neutrino mass sum can be well accommodated, as shown in table 5 and  table 6. Moreover, the predictions of these models are expected to be tested at forthcoming experiments with higher sensitivities. In most modular symmetry models, the right-handed are assumed to be singlets of modular group so that at least one parameter is introduced for each charged lepton and the hierarchies among electron, muon and tau masses can be reproduced by tuning the coupling constants. From table 4, we see that other assignments such as the triplet and double plus singlet can also be in agreement with the experimental data. The model L1 is particularly interesting, all the lepton fields L, E c and N c transform as triplet 3 under S 4 , the light neutrino mass matrix only depends on the modulus τ and an overall scale, the four coupling constants in the charged lepton mass matrices are of the same order of magnitude and the charged lepton mass hierarchies arise from the deviation from the fixed point τ = i.
Because gCP symmetry enforces all coupling constant to be real and the complex phases in the mass matrices originate from the modular forms in our models, the CP violation in leptogenesis is strongly correlated with the Dirac and Majorana CP violation phases. As an example, we have studied the leptogenesis in the model L9. By numerically solving the Boltzmann equations, we find that the baryon asymmetry of the universe, lepton masses and mixing angles can be correctly obtained in a small region of τ . The allowed range of the Dirac CP phase δ l CP would shrink significantly if the measured value of the baryon asymmetry Y B is taken into account.
As regards the quark models with S 4 modular symmetry, at least seven real coupling constants are necessary to describe the hierarchies patterns of quark masses and mixing angles, and ten minimal quark models are found, as listed in table 7. Furthermore, we investigate whether the S 4 modular symmetry can address both the lepton and quark flavor problems, and then the single complex modulus τ would be shared by the quark and lepton sectors. A typical quark-lepton unification model is presented, the lepton sector is the model L8 which contains five couplings, and the quark sector uses eight parameters. The value of τ is dominantly fixed by the experimental data of lepton masses and mixing angles. The allowed range of τ and the allowed values of quark and lepton masses and mixing parameters are very narrow, and this model can be tested in future neutrino experiments. Table 9: The representation matrices of the generatorsŜ,T ,Û as well as S, T in the five irreducible representations of S 4 in the chosen basis, where ω = e 2πi/3 is the cube root of unit.
denote 1 ≡ 1 0 , 1 ≡ 1 1 , 3 ≡ 3 0 , 3 ≡ 3 1 for the singlet or triplet representations. We shall use α i to denote the elements of first representation and β i stands for the elements of the second representation of the tensor product R 1 ⊗ R 2 , where R 1 and R 2 are two irreducible representations of S 4 .
C Classifying the minimal seesaw models with S 4 modular symmetry If the light neutrino masses originate from the type-I seesaw mechanism, the non-zero solar and atmospheric neutrino mass squared differences requires at least two right-handed neutrinos. The two right-handed neutrino models are the so-called minimal seesaw models. In the following, we shall systematically classify the neutrino superpotential for both doublet and singlet assignments of right-handed neutrinos.
• N c D = (N c 1 , N c 2 ) ∼ 2 The modular weight of N c is denoted as k N c . Since the S 4 contraction 2 ⊗ 2 → 1 is antisymmetric with respect to the two components of the doublet, the modular forms in the 1 representation can not appear in the Majorana mass terms of the right-handed neutrinos. The most general form of the superpotential for the heavy neutrino masses is given by We proceed to consider the neutrino Dirac Yukawa couplings. If the left-handed leptons transform as a triplet 3 j under S 4 with modular weight k L , the superpotential for the Dirac neutrino masses is of the form (C.10) The Dirac neutrino mass matrix can be read off as follows Under the representation transformation L : 3 j → 3 j+1 , the mass matrix M D changes into M D → diag{1, −1}M D . As a consequence, the light neutrino mass matrix given by seesaw formula is left invariant if the free coupling α is changed into −α. The transformation of the charged lepton mass matrix under L : 3 j → 3 j+1 can be read from table 2, then we know that the predictions for lepton masses and mixing matrix are preserved. Therefore it is sufficient to only consider the case of L ∼ 3 for the triplet assignment.
For the doublet plus singlet assignment of the left-handed leptons: L D ≡ (L 1 , L 2 ) ∼ 2 and L 3 ∼ 1 j whose modular weights are denoted as k L D and k L 3 respectively, the superpotential of the neutrino Yukawa coupling is 2B,≺−a H u . (C.12) Then Dirac neutrino mass matrix is given by If we change the representation L 3 : 1 j → 1 j+1 as well as the couplings α l 1 → −α l 1 , the mass matrix M D becomes M D → diag{1, −1}M D diag{1, −1, 1} and the light neutrino mass matrix transforms as M ν → diag{1, −1, 1}M ν diag{1, −1, 1}. Taking into account the charged lepton sector, the phase diag{1, −1, 1} can be absorbed by field redefinition. Thus the field L 3 can be assigned to the trivial singlet 1 of S 4 without loss of generality.
• N c 1 ∼ 1 i 1 , N c 2 ∼ 1 i 2 In this case, the superpotential of the right-handed neutrino mass terms is For the triplet assignment of left-handed lepton fields L ∼ 3 j , the superpotential of the neutrino Dirac coupling takes the following form,  The mass matrix M D is found to be In the same fashion as section 4, we have numerically analyzed the possible minimal seesaw models with S 4 modular symmetry, we find that at least eight parameters including Re τ and Im τ should be used to accommodate the experimental data of leptons. Since the resulting models contain one more free parameter than the minimal models listed in table 5 and table 6, we don't give concrete examples here.