The Minimal Seesaw Model with a Modular $S_4$ Symmetry

In this paper, we incorporate the modular $S_4$ flavor symmetry into the supersymmetric version of the minimal type-I seesaw model, in which only two right-handed neutrino singlets are introduced to account for tiny Majorana neutrino masses, and explore its implications for the lepton mass spectra, flavor mixing and CP violation. The basic idea is to assign two right-handed neutrino singlets into the unique two-dimensional irreducible representation of the modular $S_4$ symmetry group. Moreover, we show that the matter-antimatter asymmetry in our Universe can be successfully explained via the resonant leptogenesis mechanism working at a relatively-low seesaw scale $\Lambda_{\rm SS} \approx 10^7~{\rm GeV}$, with which the potential problem of the gravitino overproduction can be avoided. In this connection, we emphasize that the observed matter-antimatter asymmetry may lead to a stringent constraint on the parameter space and testable predictions for low-energy observables.


Introduction
Neutrino oscillation experiments in the past two decades have firmly established that neutrinos are massive and lepton flavor mixing is significant [1,2]. To account for tiny neutrino masses, one can naturally extend the standard model (SM) by introducing three right-handed neutrinos N iR (for i = 1, 2, 3), which are singlets under the SU(2) L × U(1) Y gauge group of the SM. Therefore, the gauge-invariant Lagrangian relevant for lepton masses and flavor mixing can be written as Although tiny neutrino masses can be well accommodated in the canonical seesaw model, the flavor structures of all the relevant lepton mass matrices are left unspecified. As a consequence, the model is in general lacking of predictive power for lepton mass spectra, flavor mixing pattern and CP violation [7]. To further reduce the number of free parameters, one can consider the so-called minimal seesaw model (MSM) [8][9][10], in which only two right-handed neutrino singlets are introduced. In such a minimal scenario, the lightest neutrino appears to be massless, namely, m 1 = 0 in the case of normal neutrino mass ordering (NO) with m 1 < m 2 < m 3 or m 3 = 0 in the case of inverted neutrino mass ordering (IO) with m 3 < m 1 < m 2 . Given ∆m 2 21 ≡ m 2 2 − m 2 1 ≈ 7.39 × 10 −5 eV 2 and |∆m 2 31 | ≡ |m 2 3 − m 2 1 | ≈ 2.53 × 10 −3 eV 2 from the global-fit analysis of current neutrino oscillation data [11], we have m 1 = 0, m 2 = ∆m 2 21 ≈ 8.6 meV and m 3 = ∆m 2 31 ≈ 50.3 meV in the MSM with NO or m 3 = 0, m 1 = |∆m 2 31 | ≈ 50.3 meV and m 2 = |∆m 2 31 | + ∆m 2 21 ≈ 50.4 meV in the MSM with IO. In addition, the MSM also provides a natural framework to realize realistic and testable neutrino mass models at the TeV scale [12,13].
On the other hand, the observed pattern of leptonic flavor mixing can be explained by further imposing non-Abelian discrete flavor symmetries on the generic Lagrangian of the seesaw model. See, e.g., Refs. [14][15][16][17][18], for recent reviews on this topic. The basic strategy for the model building with discrete flavor symmetries is to suppose that there exists an overall discrete symmetry in the 1 In this paper, we shall adopt the left-right convention for all the fermion mass terms. Therefore, the Majorana mass term of ordinary left-handed neutrinos in the Lagrangian reads −ν L M ν ν C L /2 + h.c. This convention should be taken care of when one works in the supersymmetric framework, where the superpotential is usually constructed in terms of the left-handed chiral superfields. theory at some high-energy scale, which is then broken down into distinct residual symmetries in the charged-lepton and neutrino sectors at low-energy scales [19][20][21][22][23][24]. However, the breaking pattern of the discrete symmetries usually requires a number of scalar fields (i.e., flavons), which are singlets under the SM gauge symmetries, and the flavor mixing pattern is then essentially determined by the vev's of those flavons. Therefore, these neutrino mass models with discrete flavor symmetries suffer from the problems how to experimentally pin down the free parameters associated with the flavon fields and how to directly prove the existence of the flavons themselves. In this regard, the recent suggestion for modular symmetries as a solution to the flavor mixing puzzle is attractive and deserves further investigations [25]. Under one single modular symmetry, the Yukawa couplings in the lepton sector take the modular forms, in which only one complex parameter, i.e., the modulus τ , is involved. Once the value of τ is fixed, the modular symmetry is broken and the leptonic flavor mixing pattern is determined with a limited number of free coupling constants that are unconstrained by the modular symmetry. In this case, the flavon fields are no longer needed. In the literature, there have been a great number of works on the model building based on the modular group Γ N , which for a given value of N is isomorphic to the well-known non-Abelian discrete symmetry groups, e.g., Γ 2 S 3 [26][27][28][29], Γ 3 A 4 [30][31][32][33][34][35][36][37][38][39][40], Γ 4 S 4 [41][42][43] and Γ 5 A 5 [44][45][46]. Moreover, other interesting aspects of modular symmetries have also been studied, such as the combination of modular symmetries and the generalized CP symmetry [47], multiple modular symmetries [48,49], the double covering of modular groups [50], the A 4 symmetry from the modular S 4 symmetry [51,52], the modular residual symmetry [53,54] and the unification of quark and lepton flavors with modular invariance [55].
In this paper, we investigate the generic MSM with the modular S 4 symmetry and explore its implications for lepton mass spectra, flavor mixing and CP violation. The main motivation for such an investigation is to take the advantages of the minimality of the MSM and the predictive power of the modular symmetry. Two salient features of such a theoretical setup should be emphasized. First, those two right-handed neutrinos in the MSM are naturally assigned into the two-dimensional irreducible representation of the S 4 group [43,56]. Without the introduction of the flavon fields, the number of free model parameters will be kept as minimal as possible. Second, we demonstrate that the baryon number asymmetry in the present Universe, which is characterized by the baryon-to-photon density ratio η B ≡ (n b − n b )/n γ = (6.131 ± 0.041) × 10 −10 as observed by the Planck collaboration [57], can be successfully explained via the leptogenesis mechanism [58]. Different from Ref. [39], the resonant leptogenesis [59][60][61][62][63] is implemented in the MSM with the modular S 4 symmetry such that the seesaw scale can be as low as Λ SS ≈ 10 7 GeV, evading the gravitino overproduction problem in the gauge-mediated supersymmetry breaking models [64][65][66]. As the modular symmetry is intrinsically working in the supersymmetric framework [24], the supersymmetry breaking mechanism and the gravitino problem should be properly addressed in order to consistently incorporate the thermal leptogenesis mechanism.
The remaining part of this paper is organized as follows. In Sec. 2, a brief summary of the modular symmetry and the representations of the modular S 4 symmetry is given. Two simple but viable scenarios in the MSM with one single modular S 4 symmetry are then proposed in Sec. 3, where the low-energy phenomenology of lepton mass spectra, flavor mixing pattern and CP violation in these models are also discussed. In addition, we implement the resonant leptogenesis mechanism in our models to account for the matter-antimatter asymmetry and make a connection between low-and high-energy CP violation in Sec. 4. Finally, we summarize our main conclusions in Sec. 5. Some properties of the modular S 4 symmetry group are presented in Appendix A.

Modular S 4 Symmetry
The basics of modular symmetries have been expounded in previous works [24], and the properties of the modular S 4 symmetry can be found in Refs. [42,47]. In this section, we shall give a brief introduction to the modular symmetry and establish our notations. In the supersymmetric theory with a modular symmetry, the action is given by S = d 4 xd 2 θd 2 θK(τ, τ , χ, χ) + d 4 xd 2 θW(τ, χ) + d 4 xd 2 θW(τ , χ) , (2.1) where x µ (for µ = 0, 1, 2, 3), θ α and θα (for α,α = 1, 2) are the superspace coordinates, K is the Kähler potential, and W is the superpotential. In addition, τ stands for the modulus parameter and the supermultiplet is denoted by χ, for which different supermultiplets will be distinguished by its superscript χ (I) . Under the modular transformations, the action S in Eq. (2.1) must be invariant. In Eq. (2.2), γ is the element of the modular group Γ with a, b, c and d being integers satisfying ad − bc = 1, τ is an arbitrary complex number in the upper complex plane, ρ (I) (γ) denotes the representation matrix of the modular transformation γ, and k I is the weight associated with the supermultiplet χ (I) . The modular group Γ has two generators S and T , corresponding respectively to the basic transformations τ S − → −1/τ and τ T − → τ + 1, and they fulfill the identities S 2 = (ST ) 3 = 1 with 1 being the identity element. The finite modular group of the level N is then defined as Γ N ≡ Γ/Γ(N ), where Γ(N ) is the principal congruence subgroup of Γ. As we have mentioned in the previous section, Γ N with N ≤ 5 are isomorphic to the permutation symmetry groups, namely, Γ 2 S 3 , Γ 3 A 4 , Γ 4 S 4 and Γ 5 A 5 .
The invariance of S under the transformations in Eq. (2.2) demands that K(τ, τ , χ, χ) is invariant up to the Kähler transformation K(τ, τ , χ, χ) → K(τ, τ , χ, χ) + f (τ, χ) + f (τ , χ), 2 where f (τ, χ) itself is invariant under the modular transformation. At the same time, the superpotential W(τ, χ) is invariant as well and can be expanded in terms of the supermultiplets as follows W(τ, χ) = n {I 1 ,...,I n } Y I 1 ...I n (τ )χ (I 1 ) · · · χ (I n ) , (2.3) where the coefficients Y I 1 ...I n (τ ) take the modular forms, transforming under Γ N as Y I 1 ...I n (τ ) → (cτ + d) k Y ρ Y (γ)Y I 1 ...I n (τ ) , (2.4) where ρ Y is the representation matrix of Γ N and the even integer k Y is the weight of Y I 1 ...I n (τ ). In addition, k Y and ρ Y must satisfy k Y = k I 1 + · · · + k I N and ρ Y ⊗ ρ I 1 ⊗ · · · ⊗ ρ I N 1, respectively. For the symmetry group Γ 4 S 4 of our interest, there are five linearly independent modular forms of the lowest non-trivial weight k Y = 2, denoted as Y i (τ ) for i = 1, 2, · · · , 5, which form a doublet 2 and a triplet 3 under the modular S 4 symmetry transformations [41], namely, The exact expressions of Y i (τ ) (for i = 1, 2, . . . , 5) are presented in Appendix A, where one can also find how to derive the modular forms of higher weights from the products of modular forms of the weight k Y = 2. Finally, we write down the superpotential relevant for lepton masses and flavor mixing in the minimal supersymmetric standard model (MSSM) with the seesaw extension, i.e., where L, E C and N C stand for the chiral superfields containing the left-handed lepton doublets, right-handed charged-lepton singlets and right-handed neutrino singlets, respectively. In addition, H u and H d are the Higgs superfields with the hypercharges +1/2 and −1/2, and the contraction of the SU(2) L doublets has been explicitly shown with the Levi-Civita symbol ε mn (i.e., ε 12 = −ε 21 = 1 and ε 11 = ε 22 = 0). Note that f l (τ ), f D (τ ) and f R (τ ) are modular forms, which carry the flavor index and transform nontrivially under the modular flavor symmetry, and they together with the fermion superfields form the flavor singlets, each of which is labelled by the superscript i and the subscript 1. Therefore, on the right-hand side of Eq. (2.6), the summation is taken over all possible flavor singlets i together with the corresponding coefficients α i , g i and Λ i . When the modulus parameter τ is fixed, the modular symmetry is broken down and the superpotential reads where λ l and λ D turn out to be the charged-lepton and Dirac neutrino Yukawa coupling matrices, respectively, while λ R becomes the right-handed neutrino mass matrix. As a result, the flavor structures of λ l , λ D and λ R are dictated by the modular symmetry to a large extent. After the supersymmetry breaking and the spontaneous breakdown of the SU(2) L × U(1) Y gauge symmetry, one can obtain the lepton mass terms from Eq. (2.7), whose SM version has been given just below Eq. (1.1). Converting into our left-right convention for the fermion mass terms, we find the following correspondence between the lepton mass matrices and the Yukawa coupling matrices in the MSSM framework where v d = v cos β and v u = v sin β with v ≈ 246 GeV are respectively the vev of the neutral scalar component field of H d to that of H u , with tan β ≡ v u /v d being their ratio, and " * " denotes the complex conjugation.

Low-energy Phenomenology
After setting up the framework, now we propose two viable scenarios of the MSM with the modular S 4 symmetry and explore their implications for the low-energy observables. Some general remarks on the model building are helpful.
• As we have already stressed, two right-handed neutrino singlets { N C 1 , N C 2 } can naturally be assigned into the unique two-dimensional irreducible representation 2 of the S 4 symmetry group. Such an assignment will always be assumed in the present work.
• Three lepton doublets { L 1 , L 2 , L 3 } are arranged as a triplet 3 under the S 4 symmetry, while three charged-lepton singlets { E C 1 , E C 2 , E C 3 } should be assigned into the trivial 1 or nontrivial 1 singlets of S 4 . Otherwise, it will be difficult to accommodate the strong mass hierarchy of three charged leptons, namely, m e m µ m τ .
• The Higgs doublets { H u , H d } are assigned into the trivial representation 1 under the S 4 symmetry for simplicity. In this case, their modular weights k u and k d are both vanishing. Consequently, the remaining part of the MSSM irrelevant for leptonic flavor mixing need not be changed.
Given the above assignments, we are then left with the coefficients f l (τ ), f D (τ ) and f R (τ ) in Eq. (2.6), which are modular forms with the modular weights being positive and even integers. Once the weights of supermultiplets are fixed, those of modular forms can be determined from the identity k Y = k I 1 + · · · + k I N . Since the explicit expressions of f l (τ ), f D (τ ) and f R (τ ) with specific weights are known, one can then figure out the mass matrices of both charged leptons and neutrinos, from which the lepton mass spectra and flavor mixing parameters can be extracted. Due to the freedom in choosing the weights of supermultiplets, we follow the criterion that the number of free model parameters should be as small as possible. More explicitly, we first count the number of low-energy observables: three charged-lepton masses {m e , m µ , m τ }, two independent neutrino mass-squared differences {∆m 2 21 , |∆m 2 31 |} and three mixing angles {θ 12 , θ 13 , θ 23 }. Thus the number of free model parameters should be no more than eight in order to have predictive power for the other parameters, such as the CP-violating phases. To this end, we restrict the weights of f l (τ ), f D (τ ) and f R (τ ) into the value ranging from 2 to 6. After a systematic investigation, two distinct models surviving current neutrino oscillation data have been found and will be discussed in detail in the following two subsections.

Model A
In the first model, the charge assignments of the chiral superfields and the couplings under the SU(2) L gauge symmetry and the modular S 4 symmetry have been summarized in Table 1, where the corresponding modular weights are listed in the last row. Some explanations for the assignments of all the Yukawa couplings are necessary. The charged-lepton Yukawa couplings f e (τ ), f µ (τ ), f τ (τ ) are set to be the triplet 3 with a weight of 2, 3 with a weight of 4 and 3 with a weight of 6, respectively. On the other hand, two triplets f D (τ ) and f D (τ ) with the weight of 6 will be introduced. One of them is set to 3, while the other is 3 . However, there are two distinct Table 1: The charge assignments of the chiral superfields and the relevant couplings under the SU(2) L gauge symmetry and the modular S 4 symmetry in Model A, with the corresponding modular weights listed in the last row.
3 ,1 for the 3 with a weight 6 in Eq. (A.18) will be adopted for f D (τ ). Now it is straightforward to verify that the superpotential is invariant under the modular S 4 symmetry and can be decomposed into three parts W = W l + W D + W R with where the isospin indices of the SU(2) L doublets have been suppressed and the representations of the S 4 symmetry have been explicitly indicated. By using the product rules of the S 4 symmetry group collected in Appendix A, we can obtain the charged-lepton mass matrix (3.2) and the Dirac neutrino mass matrix where the complex conjugation " * " on the right-hand side should be noticed as we have explained in Eq. (2.8). Without loss of generality, we can choose α i (for i = 1, 2, 3) in Eq. (3.2) to be real and positive. Since the overall phase of any lepton mass matrix is irrelevant for lepton masses and flavor mixing, one can take g 1 in Eq. (3.3) to be real and it is convenient to parametrize the other complex parameter as g 2 /g 1 ≡ g = ge iφ g with g = | g| and φ g ≡ arg( g). In addition, the Majorana mass matrix of right-handed neutrinos is given by where Λ is real and positive parameter characterizing the absolute scale of heavy Majorana neutrino masses. With the help of the seesaw formula M ν ≈ −M D M −1 R M T D , we finally arrive at the effective neutrino mass matrix M ν and will be able to explore its implications for lepton mass spectra and flavor mixing.
Given the complex parameter τ (or equivalently two real parameters Re τ and Im τ ), one can determine three model parameters v d α 3 / √ 2, α 1 /α 3 and α 2 /α 3 from the observed charged-lepton masses m e = 0.511 MeV, m µ = 105.7 MeV and m τ = 1776.86 MeV [2]. Since the absolute scale of neutrino masses is fixed by v 2 u g 2 1 /(2Λ), we are left with only two free parameters g and φ g , which together with Re τ and Im τ will give rise to three neutrino mixing angles {θ 12 , θ 13 , θ 23 }, the ratio of two neutrino mass-squared differences ∆m 2 21 /|∆m 2 31 |, the Dirac CP-violating phase δ and the Majorana CP-violating phase σ. As the number of free parameters is two less than that of the observables, the model under consideration should be predictive. Then we proceed to explore the phenomenological implications for lepton mass spectra, flavor mixing and CP violation. To confront our model with the latest neutrino oscillation data, we perform a numerical analysis and demonstrate that the predictions are compatible with the experimental data only in the NO case at the 3σ level. Our strategy for numerical analysis is outlined below.
• First of all, the modulus parameter τ is randomly generated in the fundamental domain G, which is defined as This domain can be identified by using the basic properties of the modular forms as clearly explained in Ref. [42]. Another interesting feature of the modular forms should be noticed. If the replacement τ → −τ * is made in Eq. (A.16), one can verify that the modular forms Y i (τ ) will change to their complex-conjugate counterparts, i.e., Y i (−τ * ) = Y * i (τ ). If we further replace g with g * in M D , all the lepton mass matrices become their complex-conjugate counterparts. Under such a transformation, the theoretical predictions for all the experimental observables keep unchanged except that the signs of all CP-violating phases should be reversed.
Once the values of {Re τ, Im τ } are randomly chosen, we can extract the parameters v d α 3 / √ 2, α 1 /α 3 and α 2 /α 3 from the following identities So far all the parameters in M l have been determined. It is then easy to diagonalize the charged-lepton mass matrix via U † l M l M † l U l = Diag m 2 e , m 2 µ , m 2 τ , from which the unitary matrix U l can be obtained.
• Next the values of the other two parameters g ∈ [0, 10] and φ g ∈ [0 • , 360 • ) are randomly generated. Therefore, the effective neutrino mass matrix M ν is determined up to the overall scale parameter v 2 u g 2 1 /(2Λ), which can be fixed by Tr M ν M † ν = m 2 2 + m 2 3 = ∆m 2 21 + ∆m 2 31 in the NO case or Tr M ν M † ν = m 2 1 + m 2 2 = 2|∆m 2 32 | − ∆m 2 21 in the IO case. In practice, we also randomly choose two neutrino mass-squared differences in their allowed ranges and then M ν is numerically known. After diagonalizing M ν via U † ν M ν U * ν = Diag {m 1 , m 2 , m 3 }, we get the unitary matrix U ν . Hence, the lepton flavor mixing matrix U = U † l U ν can be calculated by using U l and U ν . In the standard parametrization [2], we have where c ij ≡ cos θ ij and s ij ≡ sin θ ij (for ij = 12, 13, 23) have been defined, δ and σ are the Dirac and Majorana CP-violating phases, respectively. Note that the lightest neutrino mass in the MSM is vanishing, so we are left with only one Majorana CP-violating phase.
To find out the allowed parameter space of {Re τ, Im τ } and {g, φ g }, we implement the globalfit results from NuFIT 4.1 [11] without including the atmospheric neutrino data from Super-     Figure 1: Allowed ranges of the model parameters {Re τ, Im τ } and {g, φ g } and the constrained ranges of low-energy observables in the NO case in Model A, where the 3σ ranges of neutrino mixing parameters and mass-squared differences from the global-fit analysis of neutrino oscillation data have been input [11].
Kamiokande. The best-fit values of three neutrino mixing angles {θ 12 , θ 13 , θ 23 }, two neutrino mass-squared differences {∆m 2 21 , ∆m 2 31 } (or {∆m 2 21 , ∆m 2 32 }), the Dirac CP-violating phase δ and their 1σ and 3σ ranges in the NO (or IO) case are summarized in Table 2. The allowed parameter space and the constraint on the mixing parameters and other low-energy observables have been shown in Fig. 1. Some comments on the numerical results are in order.
• The allowed parameter space of {Re τ, Im τ } and {g, φ g } exists only in the NO case at the 3σ level. As one can observe from the constrained range of θ 23 in the left panel in the second row of Fig. 1, only relatively small values of θ 23 survive. In particular, the maximal value of θ 23 is no more than 42.5 • , whereas the best-fit value from the global-fit analysis is 48.3 • in the NO case or 48.6 • in the IO case.
• In addition, once the neutrino mass spectrum and the mixing parameters are known, we can predict the effective neutrino mass for beta decays, i.e., The latest result from the KATRIN experiment, where the electron energy spectrum from tritium beta decays is precisely measured, indicates m β < 1.1 eV at the 90% confidence level [68,69]. With more data accumulated in KATRIN, the upper bound will be improved to m β < 0.2 eV. However, taking into account of the 3σ ranges of neutrino mixing parameters and mass-squared differences, we obtain 8.5 meV m β 9.5 meV in our model, which is far below the future sensitivity of KATRIN.
• Furthermore, three ordinary neutrinos turn out to be Majorana particles in the seesaw model, indicating that the neutrinoless double-beta decays of some even-even heavy nuclei could take place. The effective neutrino mass relevant for neutrinoless double-beta decays is defined as which can be computed with m 1 = 0 in our model. As shown in Fig. 1, we obtain 2.2 meV m ββ 4.2 meV. It is very challenging to reach such a high sensitivity in the next-generation neutrinoless double-beta decay experiments [70].
In summary, Model A is compatible with neutrino oscillation data only in the NO case at the 3σ level, but we cannot find any viable parameter space of {Re τ, Im τ } and {g, φ g } at the 1σ level. As we shall show later, the successful leptogenesis can be realized in Model A only for an extremely-high scale of heavy neutrino masses, namely, the lightest heavy Majorana neutrino mass M 1 10 11 GeV, requiring a high reheating temperature. The charge assignments of the chiral superfields and the relevant couplings under the SU(2) L gauge symmetry and the modular S 4 symmetry in Model B, with the corresponding modular weights listed in the last row.

Model B
In Model B, the charge assignments of the superfields and the couplings under the SU(2) L gauge symmetry and the modular S 4 symmetry have been given in Table 3, where the corresponding modular weights are listed in the last row. As one can see from Table 3, the charge assignments are quite similar to those in Model A. In particular, all the modular weights are exactly the same. However, these two models differ significantly in the charged-lepton sector. Now that both E C 2 and E C 3 are assigned into 1 under the S 4 symmetry group, the representations of f µ (τ ) and f τ (τ ) are accordingly changed to be 3 . Notice that f τ (τ ) is assigned as the 3 representation with a weight of 6, so again we are facing with two choices Y 3 ,2 for f τ (τ ). Notice also that Y and Y 3 ,1 for f τ (τ ) were adopted, the charged-lepton mass matrix M l would have two different rows that are proportional to each other, reducing the rank of M l to two and thus leading to m e = 0. As this is obviously in contradiction with observations, Y 3 ,2 is the unique choice for f τ (τ ). In a similar way as for Model A, we can write down the modular invariant superpotential in the lepton sector and derive the charged-lepton mass matrix and the Dirac neutrino mass matrix where the real coefficient g 1 has been factorized out from the square brackets on the right-hand side of Eq. (3.14) and the complex parameter g ≡ g 2 /g 1 = ge iφ g has been introduced as before. In Re τ   Figure 2: Allowed ranges of the model parameters {Re τ, Im τ } and {g, φ g } and the constrained ranges of low-energy observables in the NO case in Model B, where the 1σ (yellow dots) and 3σ ranges (red dots) of neutrino mixing parameters and mass-squared differences from the global-fit analysis of neutrino oscillation data have been input [11]. The best-fit values are indicated by the black stars. addition, the Majorana mass matrix of right-handed neutrinos is given by The strategy to explore the allowed parameter space and the low-energy phenomenology in Model B follows that in Model A, so we shall not repeat it here. It turns out that both NO  Figure 3: Allowed ranges of the model parameters {Re τ, Im τ } and {g, φ g } and the constrained ranges of low-energy observables in the IO case in Model B, where the 1σ (yellow dots) and 3σ ranges (red dots) of neutrino mixing parameters and mass-squared differences from the global-fit analysis of neutrino oscillation data have been input [11]. The best-fit values are indicated by the black stars. and IO cases are allowed in Model B even at the 1σ level. The allowed ranges of the model parameters in the NO and IO cases are shown in Fig. 2 and Fig. 3, respectively. Some helpful comments on the numerical results are in order.
• To determine the model parameters from neutrino oscillation data and demonstrate how well the model is consistent with observations, we construct the χ 2 -function by regarding the best-fit values q bf j of the oscillation parameters q j ∈ {sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 , ∆m 2 21 , ∆m 2 31 , δ} from the global analysis in Ref. [11] as experimental measurements, namely, (3. 16) where p i ∈ {Re τ, Im τ, g, φ g } stand for four free model parameters, q j (p i ) denote the model predictions for observables with σ j being the symmetrized 1σ uncertainties from the globalfit analysis, which has already been given in Table 2. This setup is for the NO case, while for the IO case the oscillation parameters are instead taken to be q j ∈ {sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 , ∆m 2 21 , ∆m 2 32 , δ}. The allowed ranges of model parameters can be obtained by following the standard χ 2 -fit approach. Once the model parameters are known, we can get the constraints on the observables q j , as well as the predictions for the Majorana CP-violating phase σ and the effective neutrino mass m β for beta decays and m ββ for neutrinoless double-beta decays.
• The numerical results in the NO case are shown in Fig. 2. As one can see from two plots in the first row, the whole range [−0.5, 0.5] of Re τ is allowed at the 3σ level while the phase φ g is restricted into the narrow region [0 • , 40 • ] or [320 • , 360 • ]. However, at the 1σ level, two disconnected regions −0.39 < Re τ < −0.32 and 0.32 < Re τ < 0.39 are forbidden. This is because the mixing angle θ 23 is tightly related to Re τ and the lower limit of the 1σ range of θ 23 cannot be reached in the forbidden regions of Re τ .
• The numerical results in the IO case are shown in Fig. 3. Unlike the NO case, the allowed range of Re τ is highly constrained and only the narrow region −0.04 < Re τ < 0.04 survives at the 1σ level. As one can observe from the top-left panel, the modulus parameter τ is very close to a pure imaginary number, i.e., τ ≈ i. Moreover, the phase φ g is restricted into the region around 180 • at the 1σ level. These two parameters essentially lead to very constrained regions of the CP-violating phases, namely, δ ≈ 300 • (or 60 • ) and σ ≈ −60 • (or 60 • ). Performing the χ 2 -fit analysis in the IO case, we find that the minimum χ 2 min = 2.22 is achieved at which together with the charged-lepton masses m α (for α = e, µ, τ ) determine the other parameters v d α 3 / √ 2 = 0.121 GeV, α 1 /α 3 = 13.4 and α 2 /α 3 = 2.38 × 10 −3 . The absolute neutrino mass scale is given by v 2 u g 2 1 /(2Λ) = 60.8 meV, and the neutrino mass spectrum is found to be m 3 = 0, m 1 = 49.3 meV and m 2 = 50.1 meV. In addition, three mixing angles are θ 12 = 33.8 • , θ 13 = 8.70 • and θ 23 = 48.2 • , while two CP-violating phases are δ = 301 • and σ = −57.0 • . The effective neutrino masses turn out to be m β = 28.2 meV and m ββ = 30.4 meV, which will hopefully be tested in the next-generation beta decay and neutrinoless double-beta decay experiments.
• In Appendix A, Eq. (A. 16) indicates that the modular forms Y i (τ ) (for i = 1, 2, . . . , 5) depend on τ via q = e 2πiτ . Since all the expansion coefficients are real, the complex phases of Y i (τ ) will be determined only by Re τ . Therefore there are totally two sources of CP violation in the lepton sector, namely, Re τ and φ g . From Eq. (3.18) we can see the best-fit value of Re τ is quite small in the IO case while φ g is very close to 180 • . Then one may wonder how they can lead to relatively large values of two CP-violating phases δ and σ, as shown in Fig. 3. In fact, the contributions to δ and σ from M l and M R are negligible since the imaginary parts of M l and M R are suppressed by Re τ , as we can observe from Eqs. (3.13) and (3.15). Thus CP violation should be governed by nontrivial complex phases of M D . Although both Re τ and Im g = g sin φ g are small, φ g ≈ 180 • implies that those two matrices in the square brackets on the right-hand side of Eq. (3.14) have the opposite signs, leading to a comparable real and imaginary part of M D . This is why relatively large values of CP-violating phases δ and σ come out eventually.
In summary, Model B is perfectly consistent with current neutrino oscillation data in both NO and IO cases. A generic feature of this model is that the best-fit value of Re τ is very small (i.e., Re τ = 0.0934 in the NO case and Re τ = 4.85 × 10 −3 in the IO case). Moreover, one can observe that φ g ≈ 180 • holds in the IO case, indicating that those two matrices in the square brackets on the right-hand side of Eq. (3.14) contribute destructively to the Dirac neutrino mass matrix M D but lead to significant CP-violating phases. Since the complex phases in the lepton mass matrices are completely controlled by two complex parameters τ andg, we may be able to establish a direct connection between the CP violation at low energies and that at the high-energy seesaw scale, as will be shown in the next section.

Matter-antimatter Asymmetry
As is well known, another salient feature of the canonical seesaw model with heavy right-handed neutrinos is offering a natural and attractive framework to explain the observed matter-antimatter asymmetry in our Universe via thermal leptogenesis [58,71]. The basic idea is that the CP-violating and out-of-equilibrium decays of heavy Majorana neutrinos in the early Universe generate the lepton number asymmetry, which will be subsequently converted into the baryon number asymmetry by the baryon-and lepton-number-violating sphaleron processes [72][73][74]. In the conventional scenario of thermal leptogenesis in the MSM, a strong mass hierarchy for two heavy Majorana neutrinos N 1 and N 2 is assumed, i.e., M 1 M 2 , such that the lepton number asymmetry produced in the decays of the heavier Majorana neutrino N 2 will be washed out by the lepton-number-violating processes mediated by the lighter one N 1 . Therefore, it is the lightest heavy Majorana neutrino N 1 that is responsible for the matter-antimatter asymmetry and a successful leptogenesis requires a very large mass of N 1 , e.g., M 1 (10 11 · · · 10 13 ) GeV [75,76]. The viable scenario of thermal leptogenesis with such high-scale masses of heavy Majorana neutrinos needs the reheating temperature in the post-inflation era to be high enough T rh M 1 so as to thermalize the heavy Majorana neutrino N 1 . Nevertheless, such a high reheating temperature may cause the overproduction problem of gravitinos in the gauge-mediated supersymmetric models [64][65][66]. To avoid this problem in our models with the modular S 4 symmetry, we can lower the seesaw scale down to Λ SS ∼ 10 7 GeV by implementing the resonant leptogenesis [59][60][61][62][63], for which a nearly-degenerate mass spectrum M 2 ≈ M 1 for two heavy Majorana neutrinos is necessary. In the following discussions, we shall demonstrate that the resonant leptogenesis works well only in Model B but not in Model A.
First, we have to figure out the mass spectrum of two heavy Majorana neutrinos. Given the mass matrix M R of right-handed neutrinos in Eq.
The masses of two heavy Majorana neutrinos are found to be M 1,2 = Λ |Y 1 (τ )| 2 + |Y 2 (τ )| 2 ∓ 2Im[Y * 1 (τ )Y 2 (τ )], which depend solely on the modulus parameter τ . Notice that we always assume M 1 to be the lighter one of two mass eigenvalues, so Im[Y * 1 (τ )Y 2 (τ )] > 0 should be satisfied. For Im[Y * 1 (τ )Y 2 (τ )] < 0, two mass eigenvalues have to be changed as M 1,2 = Λ |Y 1 (τ )| 2 + |Y 2 (τ )| 2 ± 2Im[Y * 1 (τ )Y 2 (τ )] and thus the diagonalizing unitary matrix in Eq. As we have already seen from the previous section, Re τ can be restricted to very small values by neutrino oscillation data. In the limit of Re τ → 0 and Im τ > 0, the modular forms Y 1 (τ ) and Y 2 (τ ) in Eq. (A.16) well approximate to where t ≡ e −(πIm τ )/2 is defined and only the terms up to O(t 2 ) are retained. Taking the best-fit value of Im τ = 1.36 (or 1.13) in the NO (or IO) case, we get t = 0.12 (or 0.17) and thus it is safe to neglect higher-order terms of O(t 4 ). With the help of Eq. (4.3), we can obtain where it is interesting to notice that the mass degeneracy parameter ∆ is proportional to |Re τ | in the limit of Re τ → 0. Given t = 0.12 (or 0.17) in the NO (or IO) case, we have ∆ ∼ O(|Re τ |).
The correlation between the mass degeneracy parameter ∆ and Re τ in the NO and IO cases of Model B is shown in the left panel of Fig. 4 and that of Fig. 5, respectively. As given by Eq. (4.4), ∆ ∼ O(|Re τ |) is expected for small values of |Re τ |. Therefore, the correlation between ∆ and |Re τ | is perfectly linear in both NO and IO cases, and it is possible to obtain ∆ ∼ 10 −5 (or ∆ ∼ 10 −7 ) for |Re τ | ∼ 10 −5 (or |Re τ | ∼ 10 −7 ) in the NO (or IO) case. We conclude that a strong mass degeneracy of heavy Majorana neutrinos is allowed in Model B and this model is also well consistent with neutrino oscillation data at low energies. However, the values of |Re τ | in Model A are restricted by neutrino oscillation data to the region 0.38 < |Re τ | < 0.48, for which the approximate expression of ∆ in Eq. (4.4) is no longer valid. Numerically, we find 1.25 < ∆ < 1.64 for 0.38 < |Re τ | < 0.48, implying that the resonant leptogenesis cannot work well in Model A.
Although for a relatively large value of ∆ in Model A the successful leptogenesis is still possible, the lightest heavy Majorana neutrino mass will be required to be greater than 10 11 GeV. Next, we calculate the lepton number asymmetries of different lepton flavors from the CPviolating decays of heavy Majorana neutrinos. Due to the interference between the tree-and one-loop-level amplitudes, the asymmetries between the rates of heavy Majorana neutrino decays N i → α + H and those of the CP-conjugate decays N i → α + H can be found [76] iα ≡ where Γ(N i → α + H) and Γ(N i → α + H) denote the decay rates of N i into leptons and antileptons for α = e, µ, τ . In the MSSM, there are three extra sources of CP asymmetries iα , i α and i α , arising from the decays of sneutrinos into leptons and Higgsinos, neutrinos into sleptons and Higgsinos, and sneutrinos into sleptons and Higgs, respectively [77]. All these four types of CP asymmetries are identical and can be expressed as [78][79][80][81][82][83][84][85][86] where x ki ≡ M 2 k /M 2 i (for k, i = 1, 2 but k = i) and Y ν ≡ U † l Y ν U * R with Y ν ≡ √ 2M D /v u being the Dirac neutrino Yukawa coupling matrix. The effective coupling matrix Y ν has been obtained by transforming into the basis where both the charged-lepton Yukawa coupling matrix and the heavy neutrino mass matrix are diagonal. In Eq. (4.6), two relevant loop functions have been defined g(x ki ) = 2(1 − x ki ) (1 − x ki ) 2 + r 2 ki , (4.8) where r ki ≡ Γ k /M i with Γ k = ( Y † ν Y ν ) kk M k /(8π) being the tree-level total decay width of N k and it serves as the regulator to remove any singularity in the limit of exact mass degeneracy where the resonant enhancement of the CP asymmetries from the one-loop self-energy corrections can be realized by the high mass degeneracy of decaying heavy Majorana neutrinos. Finally, we proceed to estimate the baryon number asymmetry via resonant leptogenesis. Since we are interested in the scenario of thermal leptogenesis working at the energy scale of Λ SS ∼ M i ≈ 10 7 GeV, the Yukawa interactions of muon and tau charged leptons turn out to be in thermal equilibrium and the lepton asymmetries of different lepton flavors from N i decays should be distinguished [87][88][89][90][91]. The generation of the lepton number asymmetry in each flavor and its subsequent depletion by the inverse decays and the lepton-number-violating scattering processes should be studied by solving the complete set of Boltzmann equations [71,76,92]. As an order-ofmagnitude estimation, the baryon number asymmetry in the present Universe, characterized by the baryon-to-photon ratio, is given by [92] η B ≈ −1.04 × 10 −2 i α iα κ iα , (4.10) where κ iα (for i = 1, 2 and α = e, µ, τ ) are the efficiency factors taking account of the washout effects on the lepton number asymmetries. In the scenario of resonant leptogenesis, both heavy Majorana neutrinos will contribute significantly to the generation and washout of lepton number asymmetries. Hence we have κ 1α ≈ κ 2α ≈ κ(K α ), and K α ≡ K 1α + K 2α with K iα ≡ P iα K i are the flavor-dependent decay parameters, where P iα = |( Y ν ) αi | 2 /( Y † ν Y ν ) ii and K i = Γ i /H(M i ) have been defined. The Hubble parameter H(T ) = 1.66 g * (T )T 2 /M pl is evaluated at the temperature T = M i , where M pl = 1.2 × 10 19 GeV is the Planck mass and g * (T ) is the number of relativistic degrees of freedom at the temperature T . Taking into account the particle content of the MSSM, one gets g * (T ) = 228.75 in the radiation-dominated epoch.
In the parameter space of our interest, where Re τ ∈ (−1.6 · · · − 0.9) × 10 −5 or Re τ ∈ (2.0 · · · 2.5) × 10 −7 in the NO or IO case in Model B is perfectly allowed by current neutrino oscillation data, we can numerically calculate the corresponding ranges of the decay parameters K α , i.e., NO : 1.79 K e 2.01 , 26.5 K µ 28.5 , 12.8 K τ 13.8 ; One can see from Eqs. (4.11) and (4.12) that except for K e in the NO case, all the decay parameters are much larger than 1. Therefore this points to the strong washout regime, where the initial lepton number asymmetries will be efficiently destroyed. Furthermore, if the right-handed neutrinos are in thermal equilibrium at T M 1 , we have [93,94] η B ≈ −1.04 × 10 −2 α ( 1α + 2α )κ(K α ) , (4.13) with where z B (K α ) = 2 + 4K 0.13 α exp(−2.5/K α ). Note that the right-hand side of Eq. (4.14) is less by a factor of one half when compared to the formula of κ(K α ) for the SM shown in Ref. [94], since the NO IO Figure 6: Illustration for the 3σ allowed regions of Re τ and φ g , which determine the sign of the baryon-to-photon ratio η B and that of the Jarlskog invariant J . In the left panel, the dark red dots correspond to η B > 0 and J < 0, while the light red ones to either η B < 0 or J > 0, in the NO case of Model B. The results in the IO case are shown in the right panel.
associated inverse-decay reactions in the supersymmetric version of MSM will double the washout rates, reducing the asymmetry by a factor of two [76].
With the help of Eq. (4.13), we are ready to compute the baryon number asymmetry η B for both NO and IO cases in Model B. For illustration, we fix the mass of N 1 at M 1 = 10 7 GeV and tan β = 10. As we have explained before, the high mass degeneracy of heavy Majorana neutrinos is required for resonant leptogenesis to work, and thus only small values of |Re τ | are viable. For this reason, we first randomly choose the values of {Re τ, Im τ, g, φ g } within their 1σ ranges allowed by neutrino oscillation data and further impose the restriction 0 < |Re τ | < 10 −4 . Once the values of {Re τ, Im τ, g, φ g } are given, the coupling matrix Y ν and the heavy Majorana neutrino masses can be immediately determined. Then, we calculate the CP asymmetries iα by using Eq. (4.6) and the baryon-to-photon ratio η B in Eq. (4.13). Finally, we compare the prediction for η B with the observed value η B = (6.131 ± 0.041) × 10 −10 , (4.15) from the latest Planck observations [57]. In the right panels of Fig. 4 and Fig. 5, we show our numerical results of η B against the mass splitting ∆ of two heavy Majorana neutrinos in the NO and IO cases of Model B, respectively. The red dots represent the theoretical predictions for η B , while the shaded band denotes the 1σ range of the currently observed value of η B in Eq. (4.15). In the chosen region of Re τ , the required range of the mass degeneracy parameter ∆ for successful resonant leptogenesis is bounded by two vertical dashed lines. From Fig. 4 and Fig. 5, we can observe that the observed baryon number asymmetry can be well accommodated in Model B, if the mass degeneracy parameter ∆ is of O(10 −5 ) in the NO case or O(10 −7 ) in the IO case. More specifically, ∆ should be lying in the range of (1.08 · · · 1.85)×10 −5 corresponding to −1.62×10 −5 < Re τ < −0.89×10 −5 in the NO case, whereas 4.18 × 10 −7 < ∆ < 5.12 × 10 −7 corresponding to 2.00 × 10 −7 < Re τ < 2.52 × 10 −7 in the IO case. Therefore, we conclude that Model B can successfully account for the baryon number asymmetry in our Universe in both NO and IO cases via resonant leptogenesis at a relatively-low energy scale Λ SS ∼ 10 7 GeV, thus avoiding the potential problem of the gravitino overproduction. Although we have fixed the masses of heavy Majorana neutrinos at 10 7 GeV, it is straightforward to see that an even lower mass scale can be realized by further reducing the mass degeneracy parameter ∆ or equivalently Re τ .
A final remark on the correlation between low-and high-energy CP violation is helpful. Since CP violation in our models is governed by two parameters Re τ and φ g , they determine both the sign of η B and that of the Jarlskog invariant J ≡ sin θ 12 cos θ 12 sin θ 23 cos θ 23 sin θ 13 cos 2 θ 13 sin δ for leptonic CP violation in neutrino oscillations [95]. In Fig. 6, the 3σ ranges of Re τ and φ g allowed by neutrino oscillation data have been presented for NO (left panel) and IO (right panel), but only the dark red dots lead to η B > 0 and J < 0 in Model B. It is evident that we need Re τ → 0 − in the NO case (or Re τ → 0 + in the IO case) to guarantee η B > 0 as well as J < 0, the latter of which is preferred by the global-fit analysis of neutrino oscillation experiments in Table 2. Moreover, the observed baryon number asymmetry also constrains the parameter space, particularly that of Re τ , as indicated in Fig. 4 and Fig. 5, so we can explore its implications for the Dirac CP-violating phase δ. Numerically, we find a very tight constraint, namely, 208 • < δ < 212 • in the NO case and 293 • < δ < 299 • in the IO case, which will be readily tested in the future long-baseline accelerator neutrino oscillation experiments [95].

Summary
The finite modular symmetry provides us with an attractive and novel way to understand lepton flavor mixing, and has recently attracted a lot of attention. In this paper, we introduce the modular S 4 symmetry to the supersymmetric version of MSM, where only two heavy right-handed neutrino singlets are added into the particle content of the SM. Within this framework, two right-handed neutrino singlets are assigned into the unique two-dimensional representation of the Γ 4 S 4 group. Two viable models are constructed to account for lepton mass spectra, flavor mixing, and the matter-antimatter asymmetry.
First, to keep our models as economical as possible, we start with the lowest non-trivial weight k Y = 2 of all the relevant modular forms, namely, f e (τ ), f µ (τ ) and f τ (τ ) for charged-lepton Yukawa couplings, f D (τ ) and f D (τ ) for the Dirac neutrino Yukawa couplings, and f R (τ ) for the right-handed neutrino mass matrix. It turns out the minimal set of the weights of these modular forms allowed by current neutrino oscillation data is (k e , k µ , k τ , k D , k R ) = (2, 4, 6, 6, 2). Given such a setup of the weights, we have found only two viable models, denoted as Model A and Model B, for which the low-energy phenomenology has been studied in detail in Sec. 3. While Model A is consistent with the global-fit results of neutrino oscillation data only at the 3σ level in the NO case, Model B works excellently in both NO and IO cases even at the 1σ level. The allowed parameter space of the model parameters, namely, the modulus parameter τ = Re τ + iIm τ and the coupling constant g = ge iφ g has been obtained. Moreover, the constrained regions of three neutrino mixing angles {θ 12 , θ 13 , θ 23 } and two CP-violating phases {δ, σ}, as well as the predictions for the effective neutrino masses m β in beta decays and m ββ in neutrinoless double-beta decays, are given. The precision measurements of these parameters will be able to verify or rule out these two economical models.
Second, as the modular symmetry is intrinsically embedded in the supersymmetric theories, we demonstrate that the resonant leptogenesis can be implemented to successfully explain the observed baryon number asymmetry in our Universe, while avoiding the potential problem of the gravitino overproduction usually encountered in the gauge-mediated supersymmetric models. We find that the mass splitting of two heavy Majorana neutrinos is solely governed by the parameter Re τ . As an excellent approximation, the mass degeneracy parameter ∆ ≡ (M 2 − M 1 )/M 1 is proportional to |Re τ | in the limit of |Re τ | 1. The successful resonant leptogenesis at a relatively-low seesaw scale can only be realized in Model B. The reason is simply that small values of Re τ in Model A are not allowed by neutrino oscillation data. For illustration, we choose M 1 = 10 7 GeV and tan β = 10, and show that the baryon number asymmetries η B = (6.131 ± 0.041) × 10 −10 can be well accommodated in Model B for −1.62 × 10 −5 < Re τ < −0.89 × 10 −5 in the NO case or 2.00 × 10 −7 < Re τ < 2.52 × 10 −7 in the IO case. Meanwhile, the observed value of η B implies a very narrow range of the Dirac CP-violating phase, namely, 208 • < δ < 212 • in the NO case and 293 • < δ < 299 • in the IO case. This is ready for test in the future long-baseline accelerator neutrino oscillation experiments.
Further investigations of the resonant leptogenesis in our models by numerically solving the complete set of Boltzmann equations will be interesting. In this connection, it is worth stressing that successful thermal leptogenesis is one of the salient features of canonical seesaw models and must be incorporated in the models with modular symmetries. As we have already seen, the baryon number asymmetry as an observable offers an independent constraint on the model parameters. Future precision data from neutrino oscillations, beta decays, neutrinoless double-beta decays and cosmological observations will hopefully be able to narrow down the parameter space and test the theoretical predictions from modular symmetries.

A The Γ 4 S 4 Symmetry Group
The permutation symmetry group S 4 has twenty-four elements and five irreducible representations, which are denoted as 1, 1 , 2, 3 and 3 . In the present work, we choose the same basis for the representation matrices of two generators S and T as in Ref. [47], namely, In this basis, we can explicitly write down the decomposition rules of the Kronecker products of any two S 4 multiplets.