#  USTC-ICTS/PCFT-21-20 Flavor mixing and CP violation from the interplay of
$S_{4}$ modular group and gCP

Bu-Yao Qu111E-mail<EMAIL_ADDRESS>Xiang-Gan Liu222E-mail:
<EMAIL_ADDRESS>Ping-Tao Chen333E-mail<EMAIL_ADDRESS>Gui-Jun Ding444E-mail<EMAIL_ADDRESS>
Interdisciplinary Center for Theoretical Study and Department of Modern
Physics,
University of Science and Technology of China, Hefei, Anhui 230026, China

###### Abstract

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.

##  1 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
$\tau$ and the Yukawa couplings are modular forms which are holomorphic
functions of $\tau$. 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
$\Gamma=\mathrm{SL}(2,\mathbb{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 $\Gamma_{N}\equiv\bar{\Gamma}/\bar{\Gamma}(N)$ or its double covering
group $\Gamma^{\prime}_{N}\equiv\Gamma/\Gamma(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 $\Gamma_{2}\cong
S_{3}$ [4, 5, 6, 7], $\Gamma_{3}\cong A_{4}$ [2, 8, 4, 9, 10, 11, 5, 12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32],
$\Gamma_{4}\cong S_{4}$ [33, 34, 35, 36, 37, 38, 39, 21, 40], $\Gamma_{5}\cong
A_{5}$ [41, 42, 38], $\Gamma_{7}\cong\mathrm{PSL}(2,\mathbb{Z}_{7})$ [43],
$\Gamma^{\prime}_{3}\cong T^{\prime}$ [44, 45], $\Gamma^{\prime}_{4}\cong
S^{\prime}_{4}$ [46, 47] and $\Gamma^{\prime}_{5}\cong A^{\prime}_{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
[10, 6, 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
$\Gamma_{N}$ [2]. If the modular weights are general integers, the modular
forms can be organized into irreducible multiplets of the double covering
group $\Gamma^{\prime}_{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 $\mathrm{SL}(2,\mathbb{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 $\tilde{\Gamma}_{4}\cong\tilde{S}_{4}$ [55] and the
weight $k/5$ modular forms with finite metaplectic modular group
$\tilde{\Gamma}_{5}\cong A^{\prime}_{5}\times 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 $\tau$ 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
$\tau=i,-1/2+i\sqrt{3}/2,i\infty$ [58, 32, 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\geq 2$ [63]. Notice that the symplectic group coincides with the modular
group $\mathrm{SL}(2,\mathbb{Z})$ when $g=1$. In the symmetric basis for the
modular generators $S$ and $T$ with
$\rho_{\mathbf{r}}(S)=\rho^{T}_{\mathbf{r}}(S)$ and
$\rho_{\mathbf{r}}(T)=\rho^{T}_{\mathbf{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
$\mathbf{1}$ and $\mathbf{1}^{\prime}$, a doublet $\mathbf{2}$, and two
triplets $\mathbf{3}$ and $\mathbf{3}^{\prime}$. 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
$\Gamma_{4}\cong 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 $\text{Re}\tau$ and $\text{Im}\tau$. 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 $\tau$, and a benchmark model is presented. In modular
invariant models which also fulfill gCP invariance, the VEV of the modulus
$\tau$ 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 $\Gamma_{4}\cong 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
$\mathrm{SL}(2,\mathbb{Z})$ of degree two over integers:

$\mathrm{SL}(2,\mathbb{Z})=\left\\{\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\Big{|}ad-bc=1\,,\quad a,b,c,d\in\mathbb{Z}\,\right\\}\,.$ (1)

$\mathrm{SL}(2,\mathbb{Z})$ is often denoted $\Gamma$ as well. It have two
generators $S$ and $T$ with

$S=\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}\,,\qquad T=\begin{pmatrix}1&1\\\
0&1\end{pmatrix}\,,$ (2)

which satisfy the relations

$S^{4}=(ST)^{3}=\mathbb{1}\,,\quad S^{2}T=TS^{2}\,.$ (3)

The modular group $\Gamma$ has an important class of normal subgroups called
principal congruence subgroup of level $N$ which is defined as

$\Gamma(N)=\left\\{\gamma\in\mathrm{SL}(2,\mathbb{Z})~{}~{}\Big{|}~{}~{}\gamma\equiv\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\mod N\right\\}\,.$ (4)

Note that $\Gamma(1)=\mathrm{SL}(2,\mathbb{Z})$ and $T^{N}\in\Gamma(N)$. We
can obtain the finite modular group from the quotient group:

$\Gamma_{N}\equiv\bar{\Gamma}/\bar{\Gamma}(N):~{}~{}~{}~{}S^{2}=(ST)^{3}=T^{N}=1\,,~{}\quad~{}N<6\,,$
(5)

where $\bar{\Gamma}$ and $\bar{\Gamma}(N)$ are the projective groups
$\bar{\Gamma}=\Gamma/\\{\pm 1\\}$ and $\bar{\Gamma}(N)=\Gamma(N)/\\{\pm 1\\}$.
Notice that $\Gamma(N)\cong\bar{\Gamma}(N)$ for $N>2$. The group $\Gamma_{N}$
is usually called inhomogeneous finite modular group of level $N$, and
similarly the homogeneous finite modular group $\Gamma^{\prime}_{N}$ can be
defined as $\Gamma^{\prime}_{N}\equiv\Gamma/\Gamma(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\geq 6$ [64]. The group $\Gamma_{N}$ is isomorphic to
quotient group of $\Gamma^{\prime}_{N}$ over its center $\\{1,S^{2}=-1\\}$,
consequently $\Gamma^{\prime}_{N}$ is the double covering group of
$\Gamma_{N}$, and $\Gamma^{\prime}_{N}$ has twice as many elements as
$\Gamma_{N}$. Top-down constructions in string theory generally lead to
homogeneous finite modular group $\Gamma^{\prime}_{N}$ rather than the
inhomogeneous finite modular group $\Gamma_{N}$ [61, 62].

The modular group $\mathrm{SL}(2,\mathbb{Z})$ acts on the upper half plane
$\mathcal{H}=\\{\tau\in\mathbb{C}~{}|~{}\text{Im}(\tau)>0\\}$ by the linear
fractional transformation:

$\gamma\tau=\dfrac{a\tau+b}{c\tau+d}\,,\qquad\gamma=\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in\Gamma\,.$ (6)

Obviously $\gamma$ and $-\gamma$ give the same action on the modulus $\tau$,
thus each linear fractional transformation corresponds to an element of the
projective special linear group $\bar{\Gamma}$. If all the points in the orbit
of a modulus $\tau$ is identified, we obtain the coset space
$\mathcal{H}/\Gamma$ which is the so-called fundamental domain $\mathcal{D}$
of $\mathrm{SL}(2,\mathbb{Z})$:

$\mathcal{D}=\left\\{\tau\in\mathcal{H}~{}\Big{|}~{}|\tau|\geq
1\,,~{}-1/2\leq\text{Re}(\tau)\leq 1/2\right\\}\,,$ (7)

which is a hyperbolic triangle bounded by the vertical lines
$\text{Re}(\tau)=\frac{1}{2}$, $\text{Re}(\tau)=-\frac{1}{2}$ and the circle
$|\tau|=1$. Every point $\tau\in\mathcal{H}$ is equivalent to a point of
$\mathcal{D}$ via the action of $\mathrm{SL}(2,\mathbb{Z})$, and no two
distinct points inside $\mathcal{D}$ are equivalent under the action of
$\mathrm{SL}(2,\mathbb{Z})$ and two points of $\mathcal{D}$ are in the same
orbit only if they lie on the boundary of $\mathcal{D}$.

The modular form of integral weight $k$ and level $N$ is a holomorphic
function of $\tau$, and it transforms under $\Gamma(N)$ as follows

$f(\gamma\tau)=(c\tau+d)^{k}f(\tau)\,,\qquad\forall\gamma=\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in\Gamma(N)\,.$ (8)

Therefore the weight $k$ “differential form” $f(\tau)(d\tau)^{k/2}$ is
invariant under the action of every element of $\Gamma(N)$. The modular forms
of weight $k$ and level $N$ span a finite dimensional linear space
$\mathcal{M}_{k}(\Gamma(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
$\mathcal{M}(\Gamma(N))=\bigoplus^{\infty}_{k=0}\mathcal{M}_{k}(\Gamma(N))$ of
all modular forms of level $N$ forms a graded ring. Furthermore, it has been
proved that the finite dimensional space $\mathcal{M}_{2k}(\Gamma(N))$ can be
decomposed into irreducible representations of the finite modular groups
$\Gamma_{N}$ [2, 44] up to the automorphy factor $(c\tau+d)^{2k}$. That is to
say, it is always possible to choose a basis in $\mathcal{M}_{2k}(\Gamma(N))$
so that $Y^{(2k)}_{\mathbf{r}}=(f_{1}(\tau),f_{2}(\tau),\dots)^{T}$ transform
under the full modular group $\Gamma$ as

$Y^{(2k)}_{\mathbf{r}}(\gamma\tau)=(c\tau+d)^{2k}\rho_{\mathbf{r}}(\gamma)Y^{(2k)}_{\mathbf{r}}(\tau)\,,\qquad\forall\gamma=\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in\Gamma\,,$ (9)

where $\rho_{\mathbf{r}}(\gamma)$ is the irreducible representation of
quotient group $\Gamma_{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, 55, 46]:

$\displaystyle\mathcal{M}_{k}(\Gamma(4))$
$\displaystyle=\bigoplus_{a+b=2k,\,a,b\geq
0}\mathbb{C}\frac{\eta^{2b-2a}(4\tau)\eta^{5a-b}(2\tau)}{\eta^{2a}(\tau)}\,$
$\displaystyle=\bigoplus_{a+b=2k,\,a,b\geq
0}\mathbb{C}\theta^{a}_{2}(\tau)\theta^{b}_{3}(\tau)\,,$ (10)

where Dedekind eta function $\eta(\tau)$ is defined by

$\eta(\tau)=q^{1/24}\prod_{n=1}^{\infty}\left(1-q^{n}\right),\qquad q\equiv
e^{i2\pi\tau}\,,$ (11)

and the theta constants are defined as

$\theta_{2}(\tau)=\sum_{m\in\mathbb{Z}}e^{2\pi
i\tau(m+1/2)^{2}}\,,\quad\theta_{3}(\tau)=\sum_{m\in\mathbb{Z}}e^{2\pi i\tau
m^{2}}\,.$ (12)

Thus the dimension of the modular space $\mathcal{M}_{2k}(\Gamma(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 $\vartheta_{1}$ and $\vartheta_{2}$ which
are linear combinations of $\theta_{2}(\tau)$ and $\theta_{3}(\tau)$ as
follows [55],

$\displaystyle\vartheta_{1}(\tau)=\omega^{2}\theta_{3}(\tau)+(i+\omega)\theta_{2}(\tau)\,,$
$\displaystyle\vartheta_{2}(\tau)=\frac{\sqrt{2}+\sqrt{6}}{2}\theta_{3}(\tau)+e^{i\pi/4}\theta_{2}(\tau)\,,$
(13)

with $\omega=e^{2\pi i/3}$. In particular, the weight 2 modular multiplets
$Y_{\mathbf{2}}^{(2)}$ and $Y_{\mathbf{3}}^{(2)}$ can be written as [55]

$\displaystyle Y^{(2)}_{\mathbf{2}}\equiv\begin{pmatrix}Y_{1}\\\
Y_{2}\end{pmatrix}=\frac{e^{i\pi/3}}{9+6\sqrt{3}}\begin{pmatrix}2\sqrt{2}\vartheta_{1}\vartheta_{2}^{3}-\vartheta_{1}^{4}\\\
\vartheta_{2}^{4}+2\sqrt{2}\vartheta_{1}^{3}\vartheta_{2}\end{pmatrix}\,,$
$\displaystyle Y^{(2)}_{\mathbf{3}}\equiv\begin{pmatrix}Y_{3}\\\ Y_{4}\\\
Y_{5}\end{pmatrix}=\frac{e^{i\pi/3}}{3(2+\sqrt{3})}\begin{pmatrix}3\vartheta_{1}^{2}\vartheta_{2}^{2}\\\
\vartheta_{2}^{4}-\sqrt{2}\vartheta_{1}^{3}\vartheta_{2}\\\
\vartheta_{1}^{4}+\sqrt{2}\vartheta_{1}\vartheta_{2}^{3}\end{pmatrix}\,.$ (14)

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.

Modular weight $2k$ | Modular forms $Y^{(2k)}_{\mathbf{r}}$
---|---
$2k=2$ | $Y^{(2)}_{\mathbf{2}},~{}Y^{(2)}_{\mathbf{3}}$
$2k=4$ | $Y^{(4)}_{\mathbf{1}},~{}Y^{(4)}_{\mathbf{2}},~{}Y^{(4)}_{\mathbf{3}},~{}Y^{(4)}_{\mathbf{3}^{\prime}}$
$2k=6$ | $Y^{(6)}_{\mathbf{1}},~{}Y^{(6)}_{\mathbf{1}^{\prime}},~{}Y^{(6)}_{\mathbf{2}},~{}Y^{(6)}_{\mathbf{3}I},~{}Y^{(6)}_{\mathbf{3}II},~{}Y^{(6)}_{\mathbf{3}^{\prime}}$
$2k=8$ | $Y^{(8)}_{\mathbf{1}},~{}Y^{(8)}_{\mathbf{2}I},~{}Y^{(8)}_{\mathbf{2}II},~{}Y^{(8)}_{\mathbf{3}I},~{}Y^{(8)}_{\mathbf{3}II},~{}Y^{(8)}_{\mathbf{3^{\prime}}I},~{}Y^{(8)}_{\mathbf{3}^{\prime}II}$
Table 1: Summary of the even weight modular forms at level $N=4$, the
subscript $\mathbf{r}$ denotes the irreducible representations of the
inhomogeneous finite modular group $\Gamma_{4}\cong S_{4}$. Here
$Y^{(6)}_{\mathbf{3}I}$ and $Y^{(6)}_{\mathbf{3}II}$ denote the two linearly
independent weight 8 modular forms in the doublet representation $\mathbf{2}$,
and we adopt a similar notation for $Y^{(8)}_{\mathbf{2}I}$,
$Y^{(8)}_{\mathbf{2}II}$ and $Y^{(8)}_{\mathbf{3}I},~{}Y^{(8)}_{\mathbf{3}II}$
and $Y^{(8)}_{\mathbf{3^{\prime}}I},~{}Y^{(8)}_{\mathbf{3}^{\prime}II}$.

##  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 $\Gamma_{4}\cong S_{4}$ and gCP.

###  3.1  The framework

We formulate our models in the framework of modular invariant approach with
$\mathcal{N}=1$ global supersymmetry [2]. The field content consists of a set
of chiral matter superfields $\Phi_{I}$ and a modulus superfield $\tau$, their
modular transforms under $\mathrm{SL}(2,\mathbb{Z})$ are given by:

$\tau\to\gamma\tau=\frac{a\tau+b}{c\tau+d}\,,\qquad\Phi_{I}\to(c\tau+d)^{-k_{I}}\rho_{I}(\gamma)\Phi_{I}\,,$
(15)

where $-k_{I}$ is called the modular weight of the matter field $\Phi_{I}$,
and $\rho_{I}(\gamma)$ is the unitary representation of $\Gamma_{N}$. The
Kähler potential is taken to be the minimal form following the convention of
[2]

$\mathcal{K}(\Phi_{I},\bar{\Phi}_{I};\tau,\bar{\tau})=-h\Lambda^{2}\log(-i\tau+i\bar{\tau})+\sum_{I}(-i\tau+i\bar{\tau})^{-k_{I}}|\Phi_{I}|^{2}~{}~{}~{},$
(16)

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 $\tau$. 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
$\mathcal{W}(\Phi_{I},\tau)$ can be expanded in power series of the involved
supermultiplets $\Phi_{I}$,

$\mathcal{W}(\Phi_{I},\tau)=\sum_{n}Y_{I_{1}...I_{n}}(\tau)~{}\Phi_{I_{1}}...\Phi_{I_{n}}\,,$
(17)

where $Y_{I_{1}...I_{n}}$ is a modular multiplet of weight $k_{Y}$ as
introduced in previous section. Modular invariance requires that each term of
the $\mathcal{W}(\Phi_{I},\tau)$ satisfies the following conditions:

$k_{Y}=k_{I_{1}}+...+k_{I_{n}},~{}\quad~{}\rho_{Y}\otimes\rho_{I_{1}}\otimes\ldots\otimes\rho_{I_{n}}\ni\mathbf{1}\,.$
(18)

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 $\tau$ transforms under the action of gCP as [60, 61, 69, 70, 71]

$\tau\stackrel{{\scriptstyle\mathcal{CP}}}{{\longmapsto}}-\tau^{*}$ (19)

up to modular transformations. The gCP transforms a generic chiral superfield
$\Phi$ into the hermitian conjugate superfield

$\Phi(x)\stackrel{{\scriptstyle\mathcal{CP}}}{{\longmapsto}}X_{\mathbf{r}}\overline{\Phi}(x_{\mathcal{P}})\,,$
(20)

with $x=(t,\vec{x})$ and $x_{\mathcal{P}}=(t,-\vec{x})$, the CP transformation
matrix $X_{\mathbf{r}}$ is a unitary matrix acting on flavor space. The chiral
superfield $\Phi(x)$ is assigned to an irreducible unitary representation
$\rho_{\mathbf{r}}$ of the finite modular group, then the form of the matrix
$X_{\mathbf{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 $\tau$ and the matter field $\Phi$ transform as follows

$\displaystyle\hskip-28.90755pt\tau\stackrel{{\scriptstyle\mathcal{CP}}}{{\longrightarrow}}-\tau^{*}\stackrel{{\scriptstyle\gamma}}{{\longrightarrow}}-\frac{a\tau^{*}+b}{c\tau^{*}+d}\stackrel{{\scriptstyle\mathcal{CP}^{-1}}}{{\longrightarrow}}\frac{a\tau-b}{-c\tau+d}\,,$
$\displaystyle\hskip-28.90755pt\Phi(x)\stackrel{{\scriptstyle\mathcal{CP}}}{{\longrightarrow}}X_{\mathbf{r}}\overline{\Phi}(x_{\mathcal{P}})\stackrel{{\scriptstyle\gamma}}{{\longrightarrow}}(c\tau^{*}+d)^{-k}X_{\mathbf{r}}\rho^{*}_{\mathbf{r}}(\gamma)\overline{\Phi}(x_{\mathcal{P}})\stackrel{{\scriptstyle\mathcal{CP}^{-1}}}{{\longrightarrow}}(-c\tau+d)^{-k}X_{\mathbf{r}}\rho^{*}_{\mathbf{r}}(\gamma)X^{-1}_{\mathbf{r}}\Phi(x)\,,$
(21)

where $-k$ denotes the modular weight of $\Phi$. The closure of the modular
transformations and gCP transformations entails that the following consistency
condition has to be satisfied [60, 63],

$X_{\mathbf{r}}\rho^{*}_{\mathbf{r}}(\gamma)X^{-1}_{\mathbf{r}}=\chi^{-k}(\gamma)\rho_{\mathbf{r}}(u(\gamma))\,,$
(22)

where $u(\gamma)$ is an outer automorphism of the modular group,

$\gamma=\begin{pmatrix}a&b\\\ c&d\end{pmatrix}\longmapsto
u(\gamma)=\chi(\gamma)\begin{pmatrix}a&-b\\\ -c&d\end{pmatrix}\,.$ (23)

Here $\chi(\gamma)$ is called character and it is a homomorphism of
$\mathrm{SL}(2,\mathbb{Z})$ into $\left\\{+1,-1\right\\}$. 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]

$\chi(S)=\chi(T)=1\,,~{}~{}\text{or}~{}~{}\chi(S)=\chi(T)=-1\,.$ (24)

Consequently two possible gCP symmetries can be defined in the context of
modular invariance. We see that the CP transformation $X_{\mathbf{r}}$ maps
the modular group element $\gamma$ onto another element $u(\gamma)$ and the
group structure of the modular symmetry is preserved, i.e.
$u(\gamma_{1}\gamma_{2})=u(\gamma_{1})u(\gamma_{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
$\chi(S)=\chi(T)=1$, the consistency condition becomes

$X_{\mathbf{r}}\rho^{*}_{\mathbf{r}}(S)X^{-1}_{\mathbf{r}}=\rho_{\mathbf{r}}(S^{-1}),\qquad
X_{\mathbf{r}}\rho^{*}_{\mathbf{r}}(T)X^{-1}_{\mathbf{r}}=\rho_{\mathbf{r}}(T^{-1})\,,$
(25)

The second kind of gCP associated with the nontrivial character
$\chi(S)=\chi(T)=-1$ has been studied in [46, 63], the explicit form of
$X_{\mathbf{r}}$ depends on both the modular weight $-k$ and representation
assignment $\mathbf{r}$ of the matter field, and obviously it would be reduced
the first gCP for $(-1)^{-k}\rho_{\mathbf{r}}(S^{2})=1$. Notice that
$\rho_{\mathbf{r}}(S^{2})=\pm 1$ because of $S^{4}=1$. As a result, it is only
relevant for the case of $(-1)^{-k}\rho_{\mathbf{r}}(S^{2})=-1$ which implies
odd $k$ for the inhomogeneous finite modular group $\Gamma_{N}$, the gCP
transformation $X_{\mathbf{r}}$ is determined by [46, 63]

$X_{\mathbf{r}}\rho^{*}_{\mathbf{r}}(S)X^{-1}_{\mathbf{r}}=-\rho_{\mathbf{r}}(S^{-1}),\qquad
X_{\mathbf{r}}\rho^{*}_{\mathbf{r}}(T)X^{-1}_{\mathbf{r}}=-\rho_{\mathbf{r}}(T^{-1})\,,$
(26)

which can be satisfied if and only if the level $N$ is even, the dimension of
the representation $\rho_{\mathbf{r}}$ is even together with the vanishing
trace of $\rho_{\mathbf{r}}(S)$ and $\rho_{\mathbf{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 $\Gamma^{\prime}_{N}$ or
$\Gamma_{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 $\Gamma^{\prime}_{N}$ fulfilling $(-1)^{-k}\rho_{\mathbf{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_{\mathbf{r}}(\tau)\stackrel{{\scriptstyle\mathcal{CP}}}{{\longmapsto}}Y_{\mathbf{r}}(-\tau^{*})=X_{\mathbf{r}}Y^{*}(\tau)$
if the basis of the modular space is properly chosen [63].

The explicit form of the gCP transformation $X_{\mathbf{r}}$ is determined by
the consistency condition in Eq. (22) up to an overall phase for any given
irreducible representation $\mathbf{r}$. For the concerned first gCP and the
finite modular group $\Gamma_{4}\cong S_{4}$ with the basis listed in table 9,
solving the consistency conditions of Eq. (25), we find that the generalized
CP transformation $X_{\mathbf{r}}$ is in common with the representation matrix
of $S$,

$X_{\mathbf{r}}=\rho_{\mathbf{r}}(S)\,,$ (27)

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)}_{\mathbf{2}}(\tau)$ and
$Y^{(2)}_{\mathbf{3}}(\tau)$ up to weight 8, it is straightforward to check
that they transform under gCP as follows

$Y^{(k)}_{\mathbf{r}}(\tau)\stackrel{{\scriptstyle\mathcal{CP}}}{{\longmapsto}}Y^{(k)}_{\mathbf{r}}(-\tau^{*})=X_{\mathbf{r}}Y^{(k)*}_{\mathbf{r}}(\tau)\,,$
(28)

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 $\tau$. In
particular, all CP violation phases arises from non-vanishing real part of
$\tau$. In the following, we shall perform a systematical classification of
the Yukawa superpotential according to the transformation properties of the
matter fields under the $\Gamma_{4}\cong S_{4}$ modular symmetry. We assume
the Higgs doublets $H_{u}$ and $H_{d}$ are $S_{4}$ trivial singlet
$\mathbf{1}$ and their modular weights $k_{H_{u},H_{d}}$ are vanishing. Notice
that $k_{H_{u},H_{d}}$ can always be set to zero by redefining the modular
weights of matter fields.

###  3.2  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}\in\\{u^{c},d^{c},E^{c},N^{c}\\}$ and $F\in\\{Q,L\\}$. The
three generations of matter fields can be assigned to transform as a triplet
$F^{(c)}\sim\mathbf{3}^{i}$ under $S_{4}$ modular group, the direct sum of
doublet and singlet $F^{(c)}\sim\mathbf{2}\oplus\mathbf{1}^{i}$, or the direct
sum of the three singlets
$F^{(c)}\sim\mathbf{1}^{i_{1}}\oplus\mathbf{1}^{i_{2}}\oplus\mathbf{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
$\mathbf{1}\equiv\mathbf{1}^{0}$, $\mathbf{1}^{\prime}\equiv\mathbf{1}^{1}$,
$\mathbf{3}\equiv\mathbf{3}^{0}$, $\mathbf{3}^{\prime}\equiv\mathbf{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\in\\{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\in\\{1,2,3\\}$. For simplicity of the formula, we introduce
the superfluous notations $Y_{\mathbf{1},2}^{(k)}$, $Y_{\mathbf{1},3}^{(k)}$
and $Y_{\mathbf{2},3}^{(k)}$ 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)}_{\mathbf{3}I}$ and $Y^{(6)}_{\mathbf{3}II}$ in the triplet
representation $\mathbf{3}$. Furthermore, we introduce the operations $<>$ and
$\prec\succ$ and they are defined as $<i>=i~{}(\mathrm{mod}~{}2)$ and $\prec
i\succ=i~{}(\mathrm{mod}~{}3)$ which take values in the range of $\\{0,1\\}$
and $\\{1,2,3\\}$ respectively. Notice that we define $\prec i\succ=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

$\mathcal{W}_{F}=\alpha(F^{c}FH_{u/d}f(Y))_{\bf{1}}\,,$ (29)

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}$.

  * •

$F^{c}\sim\mathbf{3}^{i}\,,~{}~{}F\sim\mathbf{3}^{j}\,.$

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

$\displaystyle\mathcal{W}_{F}=$
$\displaystyle\alpha\left((F^{c}F)_{\mathbf{1}^{<i+j>}}Y_{\mathbf{1}^{<i+j>}}^{(k_{F^{c}}+k_{F})}\right)_{\mathbf{1}}H_{u/d}+\sum_{A}\beta_{A}\left((F^{c}F)_{\mathbf{2}}Y_{\mathbf{2}A}^{(k_{F^{c}}+k_{F})}\right)_{\mathbf{1}}H_{u/d}$
(30)
$\displaystyle+\sum_{l=0}^{1}\sum_{B}\gamma^{l}_{B}\left((F^{c}F)_{\mathbf{3}^{l}}Y_{\mathbf{3}^{l}B}^{(k_{F^{c}}+k_{F})}\right)_{\mathbf{1}}H_{u/d}$
$\displaystyle=$
$\displaystyle\sum_{a=1}^{3}\sum_{b=1}^{3}F^{c}_{a}F_{b}\bigg{\\{}\alpha
Y_{\mathbf{1}^{<i+j>},\prec a+b-1\succ}^{(k_{F^{c}}+k_{F})}+(-1)^{\prec
a+b+1\succ(i+j)}\sum_{A}\beta_{A}Y_{\mathbf{2}A,\prec
a+b-2\succ}^{(k_{F^{c}}+k_{F})}$
$\displaystyle+\sum_{l=0}^{1}\sum_{B}\gamma^{l}_{B}Y_{\mathbf{3}^{l}B,\prec
3-a-b\succ}^{(k_{F^{c}}+k_{F})}\left[\delta_{ab}\left(1-(-1)^{(i+j+l)}\right)-(\epsilon_{ba\prec-
b-a\succ})^{(i+j+l+1)}\right]\bigg{\\}}H_{u/d}\,.$

Here we have assumed that the modular form multiplets
$Y^{(k_{F^{c}}+k_{F})}_{\mathbf{1}}$,
$Y^{(k_{F^{c}}+k_{F})}_{\mathbf{1}^{\prime}}$,
$Y^{(k_{F^{c}}+k_{F})}_{\mathbf{2}}$, $Y^{(k_{F^{c}}+k_{F})}_{\mathbf{3}}$ and
$Y^{(k_{F^{c}}+k_{F})}_{\mathbf{3}^{\prime}}$ 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:

$\displaystyle M_{F}=$
$\displaystyle{}\alpha\left(\begin{matrix}Y_{\mathbf{1}^{<i+j>}}^{(k_{F^{c}}+k_{F})}&0&0\\\
0&0&Y_{\mathbf{1}^{<i+j>}}^{(k_{F^{c}}+k_{F})}\\\
0&Y_{\mathbf{1}^{<i+j>}}^{(k_{F^{c}}+k_{F})}&0\end{matrix}\right)v_{u/d}$ (31)
$\displaystyle+$
$\displaystyle\beta_{A}\left(\begin{matrix}0&(-1)^{i+j}Y_{\mathbf{2}A,1}^{(k_{F^{c}}+k_{F})}&Y_{\mathbf{2}A,2}^{(k_{F^{c}}+k_{F})}\\\
(-1)^{i+j}Y_{\mathbf{2}A,1}^{(k_{F^{c}}+k_{F})}&Y_{\mathbf{2}A,2}^{(k_{F^{c}}+k_{F})}&0\\\
Y_{\mathbf{2}A,2}^{(k_{F^{c}}+k_{F})}&0&(-1)^{i+j}Y_{\mathbf{2}A,1}^{(k_{F^{c}}+k_{F})}\end{matrix}\right)v_{u/d}$
$\displaystyle+$
$\displaystyle\gamma^{l}_{B}\left(\begin{matrix}[1-(-1)^{i+j+l}]Y_{\mathbf{3}^{l}B,1}^{(k_{F^{c}}+k_{F})}&(-1)^{i+j+l}Y_{\mathbf{3}^{l}B,3}^{(k_{F^{c}}+k_{F})}&-Y_{\mathbf{3}^{l}B,2}^{(k_{F^{c}}+k_{F})}\\\
-Y_{\mathbf{3}^{l}B,3}^{(k_{F^{c}}+k_{F})}&[1-(-1)^{i+j+l}]Y_{\mathbf{3}^{l}B,2}^{(k_{F^{c}}+k_{F})}&(-1)^{i+j+l}Y_{\mathbf{3}^{l}B,1}^{(k_{F^{c}}+k_{F})}\\\
(-1)^{i+j+l}Y_{\mathbf{3}^{l}B,2}^{(k_{F^{c}}+k_{F})}&-Y_{\mathbf{3}^{l}B,1}^{(k_{F^{c}}+k_{F})}&[1-(-1)^{i+j+l}]Y_{\mathbf{3}^{l}B,3}^{(k_{F^{c}}+k_{F})}\end{matrix}\right)v_{u/d}\,,$

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 $\mathcal{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.

  * •

$F^{c}\sim\mathbf{1}^{i_{1}}\oplus\mathbf{1}^{i_{2}}\oplus\mathbf{1}^{i_{3}}\,,~{}~{}F\sim\mathbf{3}^{j}\,.$

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

$\displaystyle\mathcal{W}_{F}=$
$\displaystyle{}[\alpha(F^{c}_{1}Ff_{1}(Y))_{\bf{1}}+\beta(F^{c}_{2}Ff_{2}(Y))_{\bf{1}}+\gamma(F^{c}_{3}Ff_{3}(Y))_{\bf{1}}]H_{u/d}$
(32) $\displaystyle=$
$\displaystyle\sum_{b=1}^{3}\sum_{A}\alpha_{A}F^{c}_{1}F_{b}Y_{\mathbf{3}^{<i_{1}+j>}A,\prec
2-b\succ}^{(k_{F^{c}_{1}}+k_{F})}H_{u/d}+\sum_{b=1}^{3}\sum_{B}\beta_{B}F^{c}_{2}F_{b}Y_{\mathbf{3}^{<i_{2}+j>}B,\prec
2-b\succ}^{(k_{F^{c}_{2}}+k_{F})}H_{u/d}$ $\displaystyle+$
$\displaystyle\sum_{b=1}^{3}\sum_{C}\gamma_{C}F^{c}_{3}F_{b}Y_{\mathbf{3}^{<i_{3}+j>}C,\prec
2-b\succ}^{(k_{F^{c}_{3}}+k_{F})}H_{u/d}\,.$

which leads to the following fermion mass matrix

$M_{F}=\left(\begin{matrix}\alpha_{A}Y_{\mathbf{3}^{<i_{1}+j>}A,1}^{(k_{F^{c}_{1}}+k_{F})}&\alpha_{A}Y_{\mathbf{3}^{<i_{1}+j>}A,3}^{(k_{F^{c}_{1}}+k_{F})}&\alpha_{A}Y_{\mathbf{3}^{<i_{1}+j>}A,2}^{(k_{F^{c}_{1}}+k_{F})}\\\
\beta_{B}Y_{\mathbf{3}^{<i_{2}+j>}B,1}^{(k_{F^{c}_{2}}+k_{F})}&\beta_{B}Y_{\mathbf{3}^{<i_{2}+j>}B,3}^{(k_{F^{c}_{2}}+k_{F})}&\beta_{B}Y_{\mathbf{3}^{<i_{2}+j>}B,2}^{(k_{F^{c}_{2}}+k_{F})}\\\
\gamma_{C}Y_{\mathbf{3}^{<i_{3}+j>}C,1}^{(k_{F^{c}_{3}}+k_{F})}&\gamma_{C}Y_{\mathbf{3}^{<i_{3}+j>}C,3}^{(k_{F^{c}_{3}}+k_{F})}&\gamma_{C}Y_{\mathbf{3}^{<i_{3}+j>}C,2}^{(k_{F^{c}_{3}}+k_{F})}\end{matrix}\right)v_{u/d}\,.$
(33)

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_{\mathbf{3}I}^{(8)}$ and $Y_{\mathbf{3}II}^{(8)}$ are parallel to
$Y_{\mathbf{3}}^{(2)}$ and $Y_{\mathbf{3}}^{(4)}$ respectively. As a
consequence, the fermion mass matrix for the assignment
$(k_{F^{c}_{1}}+k_{F},k_{F^{c}_{2}}+k_{F},k_{F^{c}_{3}}+k_{F})=(2,4,8)$ and
$(<i_{1}+j>,<i_{2}+j>,<i_{3}+j>)=(0,0,0)$ will have zero mass eigenvalue as
well.

  * •

$F^{c}\sim\mathbf{2}\oplus\mathbf{1}^{i}\,,~{}~{}F\sim\mathbf{3}^{j}\,.$

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:

$\displaystyle\mathcal{W}_{F}$ $\displaystyle=$
$\displaystyle{}[\alpha(F^{c}_{D}Ff_{1}(Y))_{\bf{1}}+\beta(F^{c}_{3}Ff_{2}(Y))_{\bf{1}}]H_{u/d}$
(34) $\displaystyle=$
$\displaystyle\sum_{a=1}^{2}\sum_{b=1}^{3}\sum_{l=0}^{1}\sum_{A}(-1)^{(a+1)(j+l)}\alpha^{l}_{A}F^{c}_{a}F_{b}Y_{\mathbf{3}^{l}A,\prec
2+a-b\succ}^{(k_{F^{c}_{D}}+k_{F})}H_{u/d}$
$\displaystyle~{}~{}+\sum_{b=1}^{3}\sum_{B}\beta_{B}F^{c}_{3}F_{b}Y_{\mathbf{3}^{<i+j>}B,\prec
2-b\succ}^{(k_{F^{c}_{3}}+k_{F})}H_{u/d}\,.$

The mass matrix which can be read out from this superpotential is

$M_{F}=\left(\begin{matrix}\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,2}^{(k_{F^{c}_{D}}+k_{F})}&\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,1}^{(k_{F^{c}_{D}}+k_{F})}&\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,3}^{(k_{F^{c}_{D}}+k_{F})}\\\
(-1)^{j+l}\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,3}^{(k_{F^{c}_{D}}+k_{F})}&(-1)^{j+l}\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,2}^{(k_{F^{c}_{D}}+k_{F})}&(-1)^{j+l}\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,1}^{(k_{F^{c}_{D}}+k_{F})}\\\
\beta_{B}Y_{\mathbf{3}^{<i+j>}B,1}^{(k_{F^{c}_{3}}+k_{F})}&\beta_{B}Y_{\mathbf{3}^{<i+j>}B,3}^{(k_{F^{c}_{3}}+k_{F})}&\beta_{B}Y_{\mathbf{3}^{<i+j>}B,2}^{(k_{F^{c}_{3}}+k_{F})}\end{matrix}\right)v_{u/d}\,.$
(35)

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)$.

  * •

$F^{c}\sim\mathbf{3}^{i}\,,~{}~{}F\sim\mathbf{2}\oplus\mathbf{1}^{j}\,.$

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})\sim\mathbf{2}$ and a
singlet $F_{3}\sim\mathbf{1}^{j}$, while the right-handed fermions
$F^{c}=(F^{c}_{1},F^{c}_{2},F^{c}_{3})$ 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,

$\displaystyle\mathcal{W}_{F}=$
$\displaystyle{}[\alpha(F^{c}F_{D}f_{F_{D}}(Y))_{\bf{1}}+\beta(F^{c}F_{3}f_{F_{3}}(Y))_{\bf{1}}]H_{u/d}$
(36) $\displaystyle=$
$\displaystyle\sum_{a=1}^{3}\sum_{b=1}^{2}\sum_{l=0}^{1}\sum_{A}(-1)^{(b+1)(i+l)}\alpha^{l}_{A}F_{b}F^{c}_{a}Y_{\mathbf{3}^{l}A,\prec
2+b-a\succ}^{(k_{F^{c}}+k_{F_{D}})}H_{u/d}$
$\displaystyle+\sum_{a=1}^{3}\sum_{B}\beta_{B}F_{3}F^{c}_{a}Y_{\mathbf{3}^{<i+j>}B,\prec
2-a\succ}^{(k_{F^{c}}+k_{F_{3}})}H_{u/d}\,.$

The resulting mass matrix is given by

$M_{F}=\left(\begin{matrix}\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,2}^{(k_{F^{c}}+k_{F_{D}})}&(-1)^{i+l}\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,3}^{(k_{F^{c}}+k_{F_{D}})}&\beta_{B}Y_{\mathbf{3}^{<i+j>}B,1}^{(k_{F^{c}}+k_{F_{3}})}\\\
\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,1}^{(k_{F^{c}}+k_{F_{D}})}&(-1)^{i+l}\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,2}^{(k_{F^{c}}+k_{F_{D}})}&\beta_{B}Y_{\mathbf{3}^{<i+j>}B,3}^{(k_{F^{c}}+k_{F_{3}})}\\\
\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,3}^{(k_{F^{c}}+k_{F_{D}})}&(-1)^{i+l}\alpha^{l}_{A}Y_{\mathbf{3}^{l}A,1}^{(k_{F^{c}}+k_{F_{D}})}&\beta_{B}Y_{\mathbf{3}^{<i+j>}B,2}^{(k_{F^{c}}+k_{F_{3}})}\end{matrix}\right)v_{u/d}\,,$
(37)

which is the transpose of the mass matrix in Eq. (35) with the indices $i$ and
$j$ exchanged.

  * •

$F^{c}\sim\mathbf{2}\oplus\mathbf{1}^{i}\,,~{}~{}F\sim\mathbf{2}\oplus\mathbf{1}^{j}\,.$

In this case, both left-handed and right-handed fields are assigned to the
direct sum of $S_{4}$ doublet and singlet. We denote $F_{D}=(F_{1},F_{2})$,
$F^{c}_{D}=(F^{c}_{1},F^{c}_{2})$ which transform as doublet under $S_{4}$
while $F_{3}$, $F^{c}_{3}$ are singlets. The modular weights of these fields
are $k_{F^{c}_{D}}$, $k_{F^{c}_{3}}$, $k_{F_{D}}$ and $k_{F_{3}}$. Then the
Yukawa superpotential is given by

$\displaystyle W_{F}=$
$\displaystyle{}[\alpha(F^{c}_{D}F_{D}f_{DD}(Y))_{\bf{1}}+\beta(F^{c}_{D}F_{3}f_{D3}(Y))_{\bf{1}}+\gamma(F^{c}_{3}F_{D}f_{3D}(Y))_{\bf{1}}+\delta(F^{c}_{3}F_{3}f_{33}(Y))_{\mathbf{1}}]H_{u/d}$
(38) $\displaystyle=$
$\displaystyle\Big{[}\alpha^{l}_{1}(F^{c}_{1}F_{2}+(-1)^{l}F^{c}_{2}F_{1})Y_{\mathbf{1}^{l}}^{(k_{F^{c}_{D}}+k_{F_{D}})}+\alpha_{2A}(F^{c}_{1}F_{1}Y_{\mathbf{2}A,1}^{(k_{F^{c}_{D}}+k_{F_{D}})}+F^{c}_{2}F_{2}Y_{\mathbf{2}A,2}^{(k_{F^{c}_{D}}+k_{F_{D}})})$
$\displaystyle+\beta_{B}F_{3}(F^{c}_{1}Y_{\mathbf{2}B,2}^{(k_{F^{c}_{D}}+k_{F_{3}})}+(-1)^{j}F^{c}_{2}Y_{\mathbf{2}B,1}^{(k_{F^{c}_{D}}+k_{F_{3}})})+\gamma_{C}F^{c}_{3}(F_{1}Y_{\mathbf{2}C,2}^{(k_{F^{c}_{3}}+k_{F_{D}})}$
$\displaystyle+(-1)^{i}F_{2}Y_{\mathbf{2}C,1}^{(k_{F^{c}_{3}}+k_{F_{D}})})+\delta
F^{c}_{3}F_{3}Y_{\mathbf{1}^{<i+j>}}^{(k_{F^{c}_{3}}+k_{F_{3}})}\Big{]}H_{u/d}\,,$

which leads to

$M_{F}=\left(\begin{matrix}\alpha_{2A}Y_{\mathbf{2}A,1}^{(k_{F^{c}_{D}}+k_{F_{D}})}&\alpha^{l}_{1}Y_{\mathbf{1}^{l}}^{(k_{F^{c}_{D}}+k_{F_{D}})}&\beta_{B}Y_{\mathbf{2}B,2}^{(k_{F^{c}_{D}}+k_{F_{3}})}\\\
(-1)^{l}\alpha^{l}_{1}Y_{\mathbf{1}^{l}}^{(k_{F^{c}_{D}}+k_{F_{D}})}&\alpha_{2A}Y_{\mathbf{2}A,2}^{(k_{F^{c}_{D}}+k_{F_{D}})}&(-1)^{j}\beta_{B}Y_{\mathbf{2}B,1}^{(k_{F^{c}_{D}}+k_{F_{3}})}\\\
\gamma_{C}Y_{\mathbf{2}C,2}^{(k_{F^{c}_{3}}+k_{F_{D}})}&(-1)^{i}\gamma_{C}Y_{\mathbf{2}C,1}^{(k_{F^{c}_{3}}+k_{F_{D}})}&\delta
Y_{\mathbf{1}^{<i+j>}}^{(k_{F^{c}_{3}}+k_{F_{3}})}\end{matrix}\right)v_{u/d}\,.$
(39)

The above mass matrix can be divided into four parts

$M_{F}=\left(\begin{matrix}M_{DD}&M_{D3}\\\
M_{3D}&M_{33}\end{matrix}\right)\,.$ (40)

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^{(4)}_{\mathbf{1}}$,
$Y^{(6)}_{\mathbf{1}}$, $Y^{(6)}_{\mathbf{1}^{\prime}}$ and
$Y^{(8)}_{\mathbf{1}}$ up to weight 8. Hence $M_{33}$ would be vanishing if
the following conditions are fulfilled

$\displaystyle k_{F_{3}}+k_{F^{c}_{3}}<0\,,$ (41) or
$\displaystyle{}k_{F_{3}}+k_{F^{c}_{3}}=2\,,~{}\begin{cases}\rho_{F_{3}}=\mathbf{1}\,,\rho_{F^{c}_{3}}=\mathbf{1}\\\
\rho_{F_{3}}=\mathbf{1}^{\prime}\,,\rho_{F^{c}_{3}}=\mathbf{1}^{\prime}\end{cases}$
or $\displaystyle k_{F_{3}}+k_{F^{c}_{3}}\neq
6~{}\,,\begin{cases}\rho_{F_{3}}=\mathbf{1}^{\prime}\,,\rho_{F^{c}_{3}}=\mathbf{1}\\\
\rho_{F_{3}}=\mathbf{1}\,,\rho_{F^{c}_{3}}=\mathbf{1}^{\prime}\end{cases}\,.$

Notice that odd weight $k_{F_{3}}+k_{F^{c}_{3}}=1,3,5,\ldots$ 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)$ entries, it would be vanishing if the
modular weights fulfill $k_{F^{c}_{3}}+k_{F_{D}}\leq 0$ or
$k_{F^{c}_{3}}+k_{F_{D}}=1,3,5,\ldots$. For the case of odd modular weight
$k_{F^{c}_{3}}+k_{F_{D}}=1,3,5,\ldots$, some mixing angles or masses are
vanishing555Non-zero fermion masses requires that $k_{F^{c}_{3}}+k_{F_{3}}$
and $k_{F^{c}_{D}}+k_{F_{D}}$ are even while $k_{F^{c}_{D}}+k_{F_{3}}$ 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.. 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.

  * •

$F^{c}\sim\mathbf{1}^{i_{1}}\oplus\mathbf{1}^{i_{2}}\oplus\mathbf{1}^{i_{3}}\,,~{}~{}F\sim\mathbf{2}\oplus\mathbf{1}^{j}\,.$

Analogous to previous cases, we find the Yukawa superpotential takes the
following form,

$\displaystyle\mathcal{W}_{F}$ $\displaystyle=$
$\displaystyle[\alpha(F^{c}_{1}F_{D}f_{1D}(Y))_{\bf{1}}+\beta(F^{c}_{2}F_{D}f_{2D}(Y))_{\bf{1}}+\gamma(F^{c}_{3}F_{D}f_{3D}(Y))_{\bf{1}}$
(42)
$\displaystyle+\delta_{1}(F^{c}_{1}F_{3}f_{13}(Y))_{\mathbf{1}}+\delta_{2}(F^{c}_{2}F_{3}f_{23}(Y))_{\mathbf{1}}+\delta_{3}(F^{c}_{3}F_{3}f_{33}(Y))_{\mathbf{1}}]H_{u/d}$
$\displaystyle=$
$\displaystyle\Big{[}\alpha_{A}F^{c}_{1}(F_{1}Y_{\mathbf{2}A,2}^{(k_{F^{c}_{1}}+k_{F_{D}})}+(-1)^{i_{1}}F_{2}Y_{\mathbf{2}A,1}^{(k_{F^{c}_{1}}+k_{F_{D}})})+\delta_{1}F^{c}_{1}F_{3}Y_{\mathbf{1}^{<i_{1}+j>}}^{(k_{F^{c}_{1}}+k_{F_{3}})}$
$\displaystyle+\beta_{B}F^{c}_{2}(F_{1}Y_{\mathbf{2}B,2}^{(k_{F^{c}_{2}}+k_{F_{D}})}+(-1)^{i_{2}}F_{2}Y_{\mathbf{2}B,1}^{(k_{F^{c}_{2}}+k_{F_{D}})})+\delta_{2}F^{c}_{2}F_{3}Y_{\mathbf{1}^{<i_{2}+j>}}^{(k_{F^{c}_{2}}+k_{F_{3}})}$
$\displaystyle+\gamma_{C}F^{c}_{3}(F_{1}Y_{\mathbf{2}C,2}^{(k_{F^{c}_{3}}+k_{F_{D}})}+(-1)^{i_{3}}F_{2}Y_{\mathbf{2}C,1}^{(k_{F^{c}_{3}}+k_{F_{D}})})+\delta_{3}F^{c}_{3}F_{3}Y_{\mathbf{1}^{<i_{3}+j>}}^{(k_{F^{c}_{3}}+k_{F_{3}})}\Big{]}H_{u/d}\,.$

The fermion mass matrix is determined to be

$M_{F}=\left(\begin{matrix}\alpha_{A}Y_{\mathbf{2}A,2}^{(k_{F^{c}_{1}}+k_{F_{D}})}&(-1)^{i_{1}}\alpha_{A}Y_{\mathbf{2}A,1}^{(k_{F^{c}_{1}}+k_{F_{D}})}&\delta_{1}Y_{\mathbf{1}^{<i_{1}+j>}}^{(k_{F^{c}_{1}}+k_{F_{3}})}\\\
\beta_{B}Y_{\mathbf{2}B,2}^{(k_{F^{c}_{2}}+k_{F_{D}})}&(-1)^{i_{2}}\beta_{B}Y_{\mathbf{2}B,1}^{(k_{F^{c}_{2}}+k_{F_{D}})}&\delta_{2}Y_{\mathbf{1}^{<i_{2}+j>}}^{(k_{F^{c}_{2}}+k_{F_{3}})}\\\
\gamma_{C}Y_{\mathbf{2}C,2}^{(k_{F^{c}_{3}}+k_{F_{D}})}&(-1)^{i_{3}}\gamma_{C}Y_{\mathbf{2}C,1}^{(k_{F^{c}_{3}}+k_{F_{D}})}&\delta_{3}Y_{\mathbf{1}^{<i_{3}+j>}}^{(k_{F^{c}_{3}}+k_{F_{3}})}\end{matrix}\right)v_{u/d}\,.$
(43)

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.

###  3.3  Classifying the Majorana mass terms

####

In this subsection, we explore the superpotential for the Majorana mass terms,
which can be written as

$\mathcal{W}_{F^{c}}=\Lambda(F^{c}F^{c}f(Y))_{\bf{1}}\,,$ (44)

where $\Lambda$ 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}$.

  * •

$F^{c}\sim\mathbf{3}^{k}\,.$

As shown in the Appendix A, the contraction
$\mathbf{3}^{i}\times\mathbf{3}^{i}\rightarrow\mathbf{3}$ is antisymmetric.
Thus the triplet modular forms transforming as $\mathbf{3}$ don’t contribute
to the Majorana mass terms of $F^{c}$. The superpotential
$\mathcal{W}_{F^{c}}$ reads as

$\displaystyle\mathcal{W}_{F^{c}}=$
$\displaystyle{}\left[\alpha\left((F^{c}F^{c})_{\mathbf{1}}Y_{\mathbf{1}}^{(2k_{F^{c}})}\right)_{\mathbf{1}}+\sum_{A}\beta_{A}\left((F^{c}F)_{\mathbf{2}}Y_{\mathbf{2}A}^{(2k_{F^{c}})}\right)_{\mathbf{1}}+\sum_{B}\gamma^{\prime}_{B}\left((F^{c}F^{c})_{\mathbf{3}^{\prime}}Y_{\mathbf{3}^{\prime}B}^{(2k_{F^{c}})}\right)_{\mathbf{1}}\right]\Lambda$
(45) $\displaystyle=$ $\displaystyle\sum_{a=1}^{3}\sum_{b=1}^{3}\Lambda
F^{c}_{a}F^{c}_{b}\\{\alpha Y_{\mathbf{1},\prec
a+b-1\succ}^{(2k_{F^{c}})}+\sum_{A}\beta_{A}Y_{\mathbf{2}A,\prec
a+b-2\succ}^{(2k_{F^{c}})}+\sum_{B}\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}A,\prec
3-a-b\succ}^{(2k_{F^{c}})}(3\delta_{ab}-1)\\}\,.$

The Majorana mass matrix of $F^{c}$ is symmetric

$M_{F^{c}}=\left(\begin{matrix}\alpha
Y_{\mathbf{1}}^{(2k_{F^{c}})}+2\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}B,1}^{(2k_{F^{c}})}&\beta_{A}Y_{\mathbf{2}A,1}^{(2k_{F^{c}})}-\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}B,3}^{(2k_{F^{c}})}&\beta_{A}Y_{\mathbf{2}A,2}^{(2k_{F^{c}})}-\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}B,2}^{(2k_{F^{c}})}\\\
\beta_{A}Y_{\mathbf{2}A,1}^{(2k_{F^{c}})}-\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}B,3}^{(2k_{F^{c}})}&\beta_{A}Y_{\mathbf{2}A,2}^{(2k_{F^{c}})}+2\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}B,2}^{(2k_{F^{c}})}&\alpha
Y_{\mathbf{1}}^{(2k_{F^{c}})}-\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}B,1}^{(2k_{F^{c}})}\\\
\beta_{A}Y_{\mathbf{2}A,2}^{(2k_{F^{c}})}-\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}B,2}^{(2k_{F^{c}})}&\alpha
Y_{\mathbf{1}}^{(2k_{F^{c}})}-\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}B,1}^{(2k_{F^{c}})}&\beta_{A}Y_{\mathbf{2}A,1}^{(2k_{F^{c}})}+2\gamma^{\prime}_{B}Y_{\mathbf{3}^{\prime}B,3}^{(2k_{F^{c}})}\end{matrix}\right)\Lambda\,.$
(46)

  * •

$F^{c}\sim\mathbf{2}\oplus\mathbf{1}^{i}\,.$

The contraction $\mathbf{2}\times\mathbf{2}\rightarrow\mathbf{1}^{\prime}$ is
antisymmetric and consequently it has no contribution to mass matrix. The
general superpotential for the Majorana mass is

$\displaystyle\mathcal{W}_{F^{c}}=$
$\displaystyle{}[\alpha(F^{c}_{D}F^{c}_{D}f_{DD}(Y))_{\bf{1}}+2\beta(F^{c}_{D}F^{c}_{3}f_{D3}(Y))_{\bf{1}}+\gamma(F^{c}_{3}F^{c}_{3}f_{33}(Y))_{\mathbf{1}}]\Lambda$
(47) $\displaystyle=$
$\displaystyle\Big{[}2\alpha_{1}F^{c}_{1}F^{c}_{2}Y_{\mathbf{1}^{l}}^{(2k_{F^{c}_{D}})}+\alpha_{2A}(F^{c}_{1}F^{c}_{1}Y_{\mathbf{2}A,1}^{(2k_{F^{c}_{D}})}+F^{c}_{2}F^{c}_{2}Y_{\mathbf{2}A,2}^{(2k_{F^{c}_{D}})})$
$\displaystyle+2\beta_{B}F^{c}_{3}(F^{c}_{1}Y_{\mathbf{2}B,2}^{(k_{F^{c}_{D}}+k_{F^{c}_{3}})}+(-1)^{i}F^{c}_{2}Y_{\mathbf{2}B,1}^{(k_{F^{c}_{D}}+k_{F^{c}_{3}})})+\gamma
F^{c}_{3}F^{c}_{3}Y_{\mathbf{1}}^{(2k_{F^{c}_{3}})}\Big{]}\Lambda\,,$

which gives rise to the following mass matrix,

$M_{F^{c}}=\left(\begin{matrix}\alpha_{2A}Y_{\mathbf{2}A,1}^{(2k_{F^{c}_{D}})}&\alpha_{1}Y_{\mathbf{1}}^{(2k_{F^{c}_{D}})}&\beta_{B}Y_{\mathbf{2}B,2}^{(k_{F^{c}_{D}}+k_{F^{c}_{3}})}\\\
\alpha_{1}Y_{\mathbf{1}}^{(2k_{F^{c}_{D}})}&\alpha_{2A}Y_{\mathbf{2}A,2}^{(2k_{F^{c}_{D}})}&(-1)^{i}\beta_{B}Y_{\mathbf{2}B,1}^{(k_{F^{c}_{D}}+k_{F^{c}_{3}})}\\\
\beta_{B}Y_{\mathbf{2}B,2}^{(k_{F^{c}_{3}}+k_{F^{c}_{D}})}&(-1)^{i}\beta_{B}Y_{\mathbf{2}B,1}^{(k_{F^{c}_{3}}+k_{F^{c}_{D}})}&\gamma
Y_{\mathbf{1}}^{(2k_{F^{c}_{3}})}\end{matrix}\right)\Lambda\,.$ (48)

  * •

$F^{c}\sim\mathbf{1}^{i_{1}}\oplus\mathbf{1}^{i_{2}}\oplus\mathbf{1}^{i_{3}}\,.$

In the same fashion, we can read out the most general Majorana mass terms for
the singlet assignment of $F^{c}$,

$\displaystyle\mathcal{W}_{F^{c}}$
$\displaystyle{}=\sum_{a=1}^{3}\sum_{b=1}^{3}\Lambda\alpha_{ab}F^{c}_{a}F^{c}_{b}f_{ab}(Y)$
(49)
$\displaystyle=\sum_{a=1}^{3}\Lambda\alpha_{aa}F^{c}_{a}F^{c}_{a}Y_{\mathbf{1}}^{(2k_{F^{c}_{a}})}+2\sum_{1\leq
a<b\leq
3}\Lambda\alpha_{ab}F^{c}_{a}F^{c}_{b}Y_{\mathbf{1}^{<i_{a}+i_{b}>}}^{(k_{F^{c}_{a}}+k_{F^{c}_{b}})}\,,$

and the mass matrix is

$M_{F^{c}}=\left(\begin{matrix}\alpha_{11}Y_{\mathbf{1}}^{(2k_{F^{c}_{1}})}~{}&\alpha_{12}Y_{\mathbf{1}^{<i_{1}+i_{2}>}}^{(k_{F^{c}_{1}}+k_{F^{c}_{2}})}~{}&\alpha_{13}Y_{\mathbf{1}^{<i_{1}+i_{3}>}}^{(k_{F^{c}_{1}}+k_{F^{c}_{3}})}\\\
\alpha_{12}Y_{\mathbf{1}^{<i_{1}+i_{2}>}}^{(k_{F^{c}_{1}}+k_{F^{c}_{2}})}~{}&\alpha_{22}Y_{\mathbf{1}}^{(2k_{F^{c}_{2}})}~{}&\alpha_{23}Y_{\mathbf{1}^{<i_{2}+i_{3}>}}^{(k_{F^{c}_{2}}+k_{F^{c}_{3}})}\\\
\alpha_{13}Y_{\mathbf{1}^{<i_{1}+i_{3}>}}^{(k_{F^{c}_{1}}+k_{F^{c}_{3}})}~{}&\alpha_{23}Y_{\mathbf{1}^{<i_{2}+i_{3}>}}^{(k_{F^{c}_{2}}+k_{F^{c}_{3}})}~{}&\alpha_{33}Y_{\mathbf{1}}^{(2k_{F^{c}_{3}})}\end{matrix}\right)\Lambda\,.$
(50)

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.

###  3.4 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 right-handed 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}\sim\mathbf{3}^{i}$ and $F\sim\mathbf{3}^{j}$, the mass matrix can be
read from Eq. (31) for any given modular weights. It can be easily seen that
the representation assignment
$F\sim\mathbf{3}^{<i+1>}\,,~{}F^{c}\sim\mathbf{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 $L_{D}=(L_{1},L_{2})^{T}\sim\mathbf{2}$,
$L_{3}\sim\mathbf{1}^{i_{3}}$, $E^{c}_{a}\sim\mathbf{1}^{j_{a}}$ 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

$L_{3}:\mathbf{1}^{i_{3}}\rightarrow\mathbf{1}^{<i_{3}+1>}\,,~{}E^{c}_{a}:\mathbf{1}^{j_{a}}\rightarrow\mathbf{1}^{<j_{a}+1>}\,,$
(51)

the charged lepton mass matrix would turn into

$M_{e}\rightarrow M_{e}\mathrm{diag}\\{1,-1,1\\}\,.$ (52)

Analogously changing the representation of the right-handed neutrinos, the
light neutrino mass matrix would change as

$\displaystyle\text{Dirac~{}neutrinos}:\left\\{\begin{array}[]{ll}M_{\nu}\rightarrow
M_{\nu}\mathrm{diag}\\{1,-1,1\\}&\text{for}~{}N^{c}\sim\mathbf{3}^{k}~{}\text{or}~{}N^{c}_{a}\sim\mathbf{1}^{k_{a}}\,,\\\\[7.22743pt]
M_{\nu}\rightarrow\mathrm{diag}\\{1,-1,1\\}M_{\nu}\mathrm{diag}\\{1,-1,1\\}&\text{for}~{}N^{c}\sim\mathbf{2}\oplus\mathbf{1}^{k}\,,\end{array}\right.$
(55)
$\displaystyle\text{Majorana~{}neutrinos}:M_{\nu}\rightarrow\mathrm{diag}\\{1,-1,1\\}M_{\nu}\mathrm{diag}\\{1,-1,1\\}\,.$
(56)

The phase matrix $\mathrm{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\sim\mathbf{3}$ for the
triplet assignment of the left-handed fields and
$F\sim\mathbf{2}\oplus\mathbf{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 $M_{\nu}=-M_{D}^{T}M_{N^{c}}^{-1}M_{D}$,
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.

$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}F^{c}$ $F~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ | $\mathbf{3}^{j}\rightarrow\mathbf{3}^{<j+1>}$ | $\mathbf{2}\oplus\mathbf{1}^{j}\rightarrow\mathbf{2}\oplus\mathbf{1}^{<j+1>}$ | Majorana mass matrix
---|---|---|---
$\mathbf{3}^{i}\rightarrow\mathbf{3}^{<i+1>}$ | $M_{F}\rightarrow M_{F}$ | $M_{F}\rightarrow M_{F}P$ | $M_{F^{c}}\rightarrow M_{F^{c}}$
$\mathbf{2}\oplus\mathbf{1}^{i}\rightarrow\mathbf{2}\oplus\mathbf{1}^{<i+1>}$ | $M_{F}\rightarrow PM_{F}$ | $M_{F}\rightarrow PM_{F}P$ | $M_{F^{c}}\rightarrow PM_{F^{c}}P$
$\mathbf{1}^{i_{1}}\oplus\mathbf{1}^{i_{2}}\oplus\mathbf{1}^{i_{3}}$ | $M_{F}\rightarrow M_{F}$ | $M_{F}\rightarrow M_{F}P$ | $M_{F^{c}}\rightarrow M_{F^{c}}$
$\rightarrow\mathbf{1}^{<i_{1}+1>}\oplus\mathbf{1}^{<i_{2}+1>}\oplus\mathbf{1}^{<i_{3}+1>}$

Table 2: Transformation of the fermion mass matrix under changing the
representation assignments of matter fields, and $P$ is the diagonal matrix
$\mathrm{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
$\alpha_{1}^{l}$ associated with the operators
$(F^{c}_{D}F_{D})_{\mathbf{1}^{l}}$ should be transformed into
$-\alpha^{l}_{1}$, and the coupling $\alpha_{1}$ associated with the operator
$(F^{c}_{D}F^{c}_{D})_{\mathbf{1}}$ should also change to $-\alpha_{1}$.

##  4 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 $\chi^{2}$ analysis
and we use the well-known package TMinuit to numerically search for minimum of
the $\chi^{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\sigma$ 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 $\Delta 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 $\tau$ is also treated as a free parameter to optimize the
agreement between predictions and experimental data. Since each point of
$\tau$ in the complex upper-half plane can be mapped into the fundamental
domain $\mathcal{D}$ given in Eq. (7) by a modular transformation, thus it is
sufficient to limit the modulus VEV $\langle\tau\rangle$ in the fundamental
domain $\mathcal{D}$. Under the CP transformation $\tau\rightarrow-\tau^{*}$,
gCP invariance implies that the fermion mass matrix becomes
$M_{F}(-\tau^{*})=\rho^{*}_{F^{c}}(S)M^{*}_{F}(\tau)\rho^{\dagger}_{F}(S)$ for
the charged fermion and
$M_{F^{c}}(-\tau^{*})=\rho^{*}_{F^{c}}(S)M^{*}_{F^{c}}(\tau)\rho^{\dagger}_{F^{c}}(S)$
for the Majorana mass matrix of $F^{c}$ [63]. Therefore at the CP dual point
$\tau\rightarrow-\tau^{*}$, 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 $\chi^{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_{\text{GUT}}=2\times 10^{16}$ GeV in a minimal SUSY breaking scenario, with
SUSY breaking scale $M_{\text{SUSY}}=1$ TeV and
$\tan\beta=7.5,\bar{\eta}_{b}=0.09375$. The leptonic Dirac CP phase
$\delta^{l}_{CP}$ has not been accurately measured, therefore we don’t include
the contribution of $\delta^{l}_{CP}$ in the $\chi^{2}$ function. If all
observables at the best fit point of a model are compatible with the
experimental data at $3\sigma$ 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
$\mathrm{Re}\langle\tau\rangle$ and all CP violation phases flipped their
signs while all other observables and free parameters are unchanged at the CP
dual point.

Leptons | Quarks
---|---
Observables | Central value and $1\sigma$ error | Observables | Central value and $1\sigma$ error
$m_{e}/m_{\mu}$ | $0.0048\pm 0.0002$ | $m_{u}/m_{c}$ | $(1.9286\pm 0.6017)\times 10^{-3}$
$m_{\mu}/m_{\tau}$ | $0.0565\pm 0.0045$ | $m_{c}/m_{t}$ | $(2.7247\pm 0.1200)\times 10^{-3}$
$\Delta m_{21}^{2}/10^{-5}\text{eV}^{2}$ | $7.42^{+0.21}_{-0.20}$ | $m_{d}/m_{s}$ | $(5.0528\pm 0.6192)\times 10^{-2}$
$\Delta m_{31}^{2}/10^{-3}\text{eV}^{2}$ | $2.517^{+0.026}_{-0.028}$ | $m_{s}/m_{b}$ | $(1.7684\pm 0.0975)\times 10^{-2}$
$\delta^{l}_{CP}/\pi$ | $1.0944^{+0.1500}_{-0.1333}$ | $\delta^{q}_{CP}$ | $69.213^{\circ}\pm 3.115^{\circ}$
$\sin^{2}\theta^{l}_{12}$ | $0.304^{+0.012}_{-0.012}$ | $\theta^{q}_{12}$ | $0.22736\pm 0.00073$
$\sin^{2}\theta^{l}_{13}$ | $0.02219^{+0.00062}_{-0.00063}$ | $\theta^{q}_{13}$ | $0.00338\pm 0.00012$
$\sin^{2}\theta^{l}_{23}$ | $0.573^{+0.016}_{-0.020}$ | $\theta^{q}_{23}$ | $0.03888\pm 0.00062$
Table 3: The central values and the $1\sigma$ 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\beta$ are taken from [74, 2], 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_{\text{SUSY}}=1$ TeV and $\tan\beta=7.5,\bar{\eta}_{b}=0.09375$.

###  4.1 Lepton models

For Majorana neutrinos, the minimal phenomenologically viable models only
depend on five free parameters besides the complex modulus $\tau$, and we find
nine such models labelled as L1$\sim$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 table 4. The best fit 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 $\mathbf{3}$ or
$\mathbf{3}^{\prime}$ under $S_{4}$ while the right-handed charged leptons are
in the reducible representation $\mathbf{2}\oplus\mathbf{1}^{\prime}$ 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 $\theta_{23}$ is predicted to be in the
second octant. The Dirac CP violation phase $\delta^{l}_{CP}$ is determined to
be sizable in these models, and it distributes in the range of
$[1.27\pi,1.65\pi]$. 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 $\theta_{23}$ and
$\delta^{l}_{CP}$. It is expected that a $5\sigma$ discovery of CP violation
can be reached after ten years of data taking over $50\%$ of the parameter
space. Thus our predictions for $\theta_{23}$ and $\delta^{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 $\sum m_{i}$, and the most
stringent bound is $\sum 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 $\sum m_{i}\simeq 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 $\sum_{i}m_{i}<600$ meV [82].
The limits also become weaker when one departs from the framework of
$\Lambda$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_{\beta}=\sqrt{\sum_{i}|U_{ei}|^{2}m^{2}_{i}}=\sqrt{\cos^{2}\theta^{l}_{12}\cos^{2}\theta^{l}_{13}m^{2}_{1}+\sin^{2}\theta^{l}_{12}\cos^{2}\theta^{l}_{13}m^{2}_{2}+\sin^{2}\theta^{l}_{13}m^{2}_{3}}~{}\,,$
(57)

where $U$ is the lepton mixing matrix. From the values of lepton mixing angles
and neutrino masses, we can determine the effective mass $m_{\beta}$, as shown
in table 5. The predictions for $m_{\beta}$ are around 15 meV, and they are
much below the upper limit $m_{\beta}<1.1$ eV given by KATRIN [83]. It is
expected that KATRIN can advance the sensitivity on $m_{\beta}$ 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\nu\beta\beta$) 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_{\beta\beta}$ which is given by

$\displaystyle m_{\beta\beta}$ $\displaystyle=$
$\displaystyle\Big{|}\sum_{i}U^{2}_{ei}m_{i}\Big{|}$ (58) $\displaystyle=$
$\displaystyle\Big{|}\cos^{2}\theta^{l}_{12}\cos^{2}\theta^{l}_{13}m_{1}+\sin^{2}\theta^{l}_{12}\cos^{2}\theta^{l}_{13}e^{i\alpha_{21}}m_{2}+\sin^{2}\theta^{l}_{13}e^{i(\alpha_{31}-2\delta^{l}_{CP})}m_{3}\Big{|}\,.$

The strongest bound on $m_{\beta\beta}$ is set by the KamLAND-Zen experiment
$m_{\beta\beta}<(61-165)$ meV [85], where the largest uncertainty arises from
the computation of the associated nuclear matrix element. There are many
$0\nu\beta\beta$ decay experiments in the plan and construction, which aim to
improve the current bounds on $m_{\beta\beta}$. The future large scale
$0\nu\beta\beta$ 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_{\beta\beta}$ sensitivity is $5.7-17.7$ meV after 10 years of
data taking [88]. Using the master formula of Eq. (58), we can determine the
values of the effective Majorana neutrino mass $m_{\beta\beta}$ 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\nu\beta\beta$ experiments except the models L5 and L6
which are experimentally very challenging because of the quite low values of
$m_{\beta\beta}$.

It is remark that these minimal viable models only use five real parameters
together with the complex modulus $\tau$ 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
$\tau$ up to an overall scale in the model L1 while there are four real
couplings in the charged lepton superpotential with
$\mathcal{W}_{e}=\alpha^{e}(E^{c}L)_{\mathbf{1}}Y^{(4)}_{\mathbf{1}}+\beta^{e}(E^{c}LY^{(4)}_{\mathbf{2}})_{\mathbf{1}}+\gamma^{e}(E^{c}LY^{(4)}_{\mathbf{3}})_{\mathbf{1}}+\gamma^{\prime
e}(E^{c}LY^{(4)}_{\mathbf{3}^{\prime}})_{\mathbf{1}}$. It is notable that the
hierarchical masses of charged leptons can be reproduced although the four
coupling constants $\alpha^{e}$, $\beta^{e}$, $\gamma^{e}$ and $\gamma^{\prime
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 $\langle\tau\rangle$
from the self-dual fixed point $\tau=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\sigma$ 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\nu\beta\beta$ 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 right-handed 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 $\tau$, and four minimal models are
found. The modular transformation properties of the lepton fields and the
results of the $\chi^{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 $\theta^{l}_{23}$ and $\delta^{l}_{CP}$ at future long baseline
neutrino oscillation experiments [77, 78, 79, 80, 81]. The effective neutrino
mass $m_{\beta}$ in beta decay is predicted to be an order of magnitude below
the expected sensitivity of the KATRIN experiment [83].

| L1 | L2 | L3 | L4 | L5 | L6 | L7 | L8 | L9
---|---|---|---|---|---|---|---|---|---
$\rho_{L}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}$
$\rho_{E^{c}}$ | $\mathbf{3}$ | $\mathbf{1}\oplus\mathbf{1}^{\prime}\oplus\mathbf{1}$ | $\mathbf{1}\oplus\mathbf{1}^{\prime}\oplus\mathbf{1}^{\prime}$ | $\mathbf{1}\oplus\mathbf{1}^{\prime}\oplus\mathbf{1}$ | $\mathbf{1}\oplus\mathbf{1}^{\prime}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}^{\prime}$ | $\mathbf{2}\oplus\mathbf{1}^{\prime}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}^{\prime}$
$\rho_{N^{c}}$ | $\mathbf{3}$ | $\mathbf{3}^{\prime}$ | $\mathbf{3}^{\prime}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}^{\prime}$ | $\mathbf{3}$ | $\mathbf{3}$
$k_{L}$ | $-1$ | $2$ | $2$ | $-1$ | $-1$ | $-1$ | $2$ | $1,1$ | $1,1$
$k_{E^{c}}$ | $5$ | $0,2,2$ | $0,2,4$ | $3,5,7$ | $3,5,9$ | $7,5$ | $2,2$ | $1,1$ | $1,1$
$k_{N^{c}}$ | $1$ | $0$ | 0 | $1$ | $1$ | $1$ | $0$ | $1$ | $1$

Table 4:  Summary of the representation and modular weight assignments of the
matter fields in the minimal phenomenologically viable lepton models based on
$S_{4}$ modular symmetry and gCP symmetry, the neutrinos are assumed to be
Majorana particles. Notice that the Higgs fields are invariant under $S_{4}$
with zero modular weight.

Models | L1 | Models | L2 | L3 | Models | L4 | L5
---|---|---|---|---|---|---|---
$\mathrm{Re}\langle\tau\rangle$ | $-0.187799$ | $\mathrm{Re}\langle\tau\rangle$ | 0.101211 | 0.101211 | $\mathrm{Re}\langle\tau\rangle$ | $-0.178216$ | 0.0667032
$\mathrm{Im}\langle\tau\rangle$ | 1.08926 | $\mathrm{Im}\langle\tau\rangle$ | 1.01587 | 1.01587 | $\mathrm{Im}\langle\tau\rangle$ | 1.09848 | 1.17412
$\beta^{e}/\alpha^{e}$ | $-0.997574$ | $\beta^{e}/\alpha^{e}$ | 11.1874 | 11.1873 | $\beta^{e}/\alpha^{e}$ | 1441.41 | 1012.34
$\gamma^{e}/\alpha^{e}$ | $-0.978209$ | $\gamma^{e}/\alpha^{e}$ | 0.00260186 | 0.00209056 | $\gamma^{e}_{1}/\alpha^{e}$ | 48.4687 | 7.25588
$\gamma^{\prime}{}^{e}/\alpha^{e}$ | $-0.288994$ | $\gamma^{D}/\beta^{D}$ | 0.0139677 | 0.0139677 | $\gamma^{e}_{2}/\alpha^{e}$ | $-49.8784$ | $-55.3767$
$\alpha^{e}v_{d}(\mathrm{MeV})$ | 228.251 | $\alpha^{e}v_{d}(\mathrm{MeV})$ | 42.0558 | 42.0555 | $\alpha^{e}v_{d}(\mathrm{MeV})$ | 0.384973 | 0.629924
$\frac{(\alpha^{D}v_{u})^{2}}{\beta^{N}\Lambda}(\mathrm{meV})$ | 28.4038 | $\frac{(\beta^{D}v_{u})^{2}}{\alpha^{N}\Lambda}(\mathrm{meV})$ | 9.40363 | 9.40363 | $\frac{(\alpha^{D}v_{u})^{2}}{\beta^{N}\Lambda}(\mathrm{meV})$ | 27.5524 | 22.6716
$\sin^{2}\theta^{l}_{12}$ | 0.320753 | $\sin^{2}\theta_{12}$ | 0.305480 | 0.305480 | $\sin^{2}\theta^{l}_{12}$ | 0.307702 | 0.304328
$\sin^{2}\theta^{l}_{13}$ | 0.0219057 | $\sin^{2}\theta_{13}$ | 0.0221663 | 0.0221663 | $\sin^{2}\theta^{l}_{13}$ | 0.0221664 | 0.0221541
$\sin^{2}\theta^{l}_{23}$ | 0.521416 | $\sin^{2}\theta_{23}$ | 0.485929 | 0.485929 | $\sin^{2}\theta^{l}_{23}$ | 0.503135 | 0.491327
$\delta^{l}_{CP}/\pi$ | 1.33248 | $\delta^{l}_{CP}/\pi$ | 1.64135 | 1.64135 | $\delta^{l}_{CP}/\pi$ | 1.32001 | 1.46324
$\alpha_{21}/\pi$ | 1.32139 | $\alpha_{21}/\pi$ | 0.353191 | 0.353191 | $\alpha_{21}/\pi$ | 1.31196 | 1.14629
$\alpha_{31}/\pi$ | 0.528182 | $\alpha_{31}/\pi$ | 1.25878 | 1.25878 | $\alpha_{31}/\pi$ | 0.509414 | 1.99551
$m_{1}/\mathrm{meV}$ | 14.0612 | $m_{1}/\mathrm{meV}$ | 12.1830 | 12.1830 | $m_{1}/\mathrm{meV}$ | 13.6298 | 11.1822
$m_{2}/\mathrm{meV}$ | 16.4808 | $m_{2}/\mathrm{meV}$ | 14.9106 | 14.9106 | $m_{2}/\mathrm{meV}$ | 16.1143 | 14.1047
$m_{3}/\mathrm{meV}$ | 51.7056 | $m_{3}/\mathrm{meV}$ | 51.5812 | 51.5812 | $m_{3}/\mathrm{meV}$ | 51.8358 | 51.2977
$m_{\beta}/\mathrm{meV}$ | 16.5872 | $m_{\beta}/\mathrm{meV}$ | 15.0395 | 15.0395 | $m_{\beta}/\mathrm{meV}$ | 16.2311 | 14.2324
$m_{\beta\beta}/\mathrm{meV}$ | 9.03230 | $m_{\beta\beta}/\mathrm{meV}$ | 12.0458 | 12.0458 | $m_{\beta\beta}/\mathrm{meV}$ | 8.80743 | 3.45226
$\chi^{2}_{\text{min}}$ | 9.01 | $\chi^{2}_{\text{min}}$ | 18.97 | 18.97 | $\chi^{2}_{\text{min}}$ | 12.31 | 16.68

Models | L6 | Models | L7 | Models | L8 | L9
---|---|---|---|---|---|---
$\mathrm{Re}\langle\tau\rangle$ | 0.0546503 | $\mathrm{Re}\langle\tau\rangle$ | 0.101532 | $\mathrm{Re}\langle\tau\rangle$ | $-0.482373$ | 0.193695
$\mathrm{Im}\langle\tau\rangle$ | 1.17782 | $\mathrm{Im}\langle\tau\rangle$ | 1.01583 | $\mathrm{Im}\langle\tau\rangle$ | 1.27223 | 0.991159
$\beta^{e}/\alpha^{e}_{1}$ | 0.0248284 | $\beta^{e}/\alpha^{e}$ | 37.5920 | $\beta^{e}/\alpha^{e}$ | 2834.28 | 2833.91
$\alpha^{e}_{2}/\alpha^{e}_{1}$ | $-0.901642$ | $\alpha^{\prime}{}^{e}/\alpha^{e}$ | $-1.01173$ | $\gamma^{e}/\alpha^{e}$ | 160.150 | 160.149
$\alpha^{\prime}{}^{e}/\alpha^{e}_{1}$ | $-1.17279$ | $\gamma^{D}/\beta^{D}$ | 0.0139492 | $\beta^{D}/\alpha^{D}$ | 1.07525 | $-1.07524$
$\alpha^{e}_{1}v_{d}(\mathrm{MeV})$ | 202.486 | $\alpha^{e}v_{d}(\mathrm{MeV})$ | 12.5174 | $\alpha^{e}v_{d}(\mathrm{MeV})$ | 0.459000 | 0.394985
$\frac{(\alpha^{D}v_{u})^{2}}{\beta^{N}\Lambda}(\mathrm{meV})$ | 22.5089 | $\frac{(\beta^{D}v_{u})^{2}}{\alpha^{N}\Lambda}(\mathrm{meV})$ | 9.40482 | $\frac{(\alpha^{D}v_{u})^{2}}{\beta^{N}\Lambda}(\mathrm{meV})$ | 12.2127 | 6.90400
$\sin^{2}\theta^{l}_{12}$ | 0.303992 | $\sin^{2}\theta_{12}$ | 0.305479 | $\sin^{2}\theta^{l}_{12}$ | 0.301955 | 0.301963
$\sin^{2}\theta^{l}_{13}$ | 0.0221904 | $\sin^{2}\theta_{13}$ | 0.0221666 | $\sin^{2}\theta^{l}_{13}$ | 0.022121 | 0.0221229
$\sin^{2}\theta^{l}_{23}$ | 0.573536 | $\sin^{2}\theta_{23}$ | 0.485966 | $\sin^{2}\theta^{l}_{23}$ | 0.613488 | 0.613493
$\delta^{l}_{CP}/\pi$ | 1.35611 | $\delta^{l}_{CP}/\pi$ | 1.64081 | $\delta^{l}_{CP}/\pi$ | 1.27970 | 1.27971
$\alpha_{21}/\pi$ | 1.11691 | $\alpha_{21}/\pi$ | 0.353732 | $\alpha_{21}/\pi$ | 1.26955 | 1.26956
$\alpha_{31}/\pi$ | 1.99179 | $\alpha_{31}/\pi$ | 1.25913 | $\alpha_{31}/\pi$ | 0.411336 | 0.411352
$m_{1}/\mathrm{meV}$ | 11.1010 | $m_{1}/\mathrm{meV}$ | 12.1879 | $m_{1}/\mathrm{meV}$ | 30.6674 | 30.6690
$m_{2}/\mathrm{meV}$ | 14.0403 | $m_{2}/\mathrm{meV}$ | 14.9146 | $m_{2}/\mathrm{meV}$ | 31.8495 | 31.8511
$m_{3}/\mathrm{meV}$ | 51.2817 | $m_{3}/\mathrm{meV}$ | 51.5824 | $m_{3}/\mathrm{meV}$ | 58.6517 | 58.6547
$m_{\beta}/\mathrm{meV}$ | 14.1711 | $m_{\beta}/\mathrm{meV}$ | 15.0435 | $m_{\beta}/\mathrm{meV}$ | 31.8999 | 31.9016
$m_{\beta\beta}/\mathrm{meV}$ | 3.77388 | $m_{\beta\beta}/\mathrm{meV}$ | 12.0465 | $m_{\beta\beta}/\mathrm{meV}$ | 17.6025 | 17.6035
$\chi^{2}_{\text{min}}$ | 0.00134 | $\chi^{2}_{\text{min}}$ | 18.96 | $\chi^{2}_{\text{min}}$ | 6.45 | 6.45

Table 5: The best fit values of the input parameters for the minimal lepton
models listed in table 4, where neutrinos are assumed to be Majorana
particles. We give the values of the neutrino mixing angles, Dirac and
Majorana CP violating phases and the neutrino masses at the best fitting
points. The notations $m_{\beta}$ and $m_{\beta\beta}$ denote the effective
neutrino masses measured in beta decay and neutrinoless double decay
respectively. Note that the transformation $\tau\rightarrow-\tau^{*}$ leaves
all observables unchanged except shifting the signs of the CP phases
$\delta^{l}_{CP}$, $\alpha_{21}$ and $\alpha_{31}$.  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\sigma$ regions
[76]. Notice that the transformation $\tau\rightarrow-\tau^{*}$ leaves all
observables unchanged except shifting the signs of the CP phases
$\delta^{l}_{CP}$, $\alpha_{21}$ and $\alpha_{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\sigma$ ranges. The present upper limit
$m_{\beta\beta}<(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 $\sum_{i}m_{i}<0.120$eV at 95% confidence level
obtained by the Planck collaboration [82].  Figure 2: The correlations among
the input parameters, lepton mixing angles, CP violation phases and neutrino
masses in the model L6, where we adopt the same convention as figure 1.

Models | $\rho_{L}$ | $\rho_{E^{c}}$ | $\rho_{N^{c}}$ | $k_{L}$ | $k_{E^{c}}$ | $k_{N^{c}}$
---|---|---|---|---|---|---
D1 | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{3}$ | $\mathbf{3}$ | $k,k$ | $4-k$ | $2-k\,(k>2)$
D2 | $\mathbf{3}$ | $\mathbf{2}\oplus\mathbf{1}^{\prime}$ | $\mathbf{3}^{\prime}$ | $k$ | $4-k,4-k$ | $2-k\,(k>2)$
D3 | $\mathbf{3}$ | $\mathbf{1}\oplus\mathbf{1}^{\prime}\oplus\mathbf{1}$ | $\mathbf{3}^{\prime}$ | $k$ | $2-k,4-k,4-k$ | $2-k\,(k>2)$
D4 | $\mathbf{3}$ | $\mathbf{1}\oplus\mathbf{1}^{\prime}\oplus\mathbf{1^{\prime}}$ | $\mathbf{3}^{\prime}$ | $k$ | $2-k,4-k,6-k$ | $2-k\,(k>2)$

Model | D1 | Model | D2 | Models | D3 | D4
---|---|---|---|---|---|---
$\mathrm{Re}\langle\tau\rangle$ | 0.341595 | $\mathrm{Re}\langle\tau\rangle$ | 0.106064 | $\mathrm{Re}\langle\tau\rangle$ | 0.105759 | 0.105759
$\mathrm{Im}\langle\tau\rangle$ | 1.36933 | $\mathrm{Im}\langle\tau\rangle$ | 1.00322 | $\mathrm{Im}\langle\tau\rangle$ | 1.00325 | 1.00325
$\beta^{\nu}/\alpha^{\nu}$ | 1.02513 | $\beta^{e}/\alpha^{e}_{1}$ | 37.5323 | $\beta^{e}/\alpha^{e}$ | 11.0459 | 11.0458
$\alpha^{e}_{2}/\alpha^{e}_{1}$ | $-0.984449$ | $\gamma^{\nu}/\beta^{\nu}$ | 0.00502027 | $\gamma^{e}/\alpha^{e}$ | 0.00258625 | 0.00203595
$\beta^{e}/\alpha^{e}_{1}$ | 21.8849 | $\alpha^{e}_{2}/\alpha^{e}_{1}$ | $-1.01168$ | $\gamma^{\nu}/\beta^{\nu}$ | 0.00502515 | 0.00502516
$\alpha^{e}_{1}v_{d}(\mathrm{MeV})$ | 22.8857 | $\alpha^{e}_{1}v_{d}(\mathrm{MeV})$ | 12.2227 | $\alpha^{e}v_{d}(\mathrm{MeV})$ | 41.5209 | 41.5205
$\alpha^{\nu}v_{u}(\mathrm{meV})$ | 18.5350 | $\beta^{\nu}v_{u}(\mathrm{meV})$ | 24.5106 | $\beta^{\nu}v_{u}(\mathrm{meV})$ | 24.5128 | 24.5128
$\sin^{2}\theta^{l}_{12}$ | 0.302259 | $\sin^{2}\theta_{12}$ | 0.307374 | $\sin^{2}\theta^{l}_{12}$ | 0.307295 | 0.307294
$\sin^{2}\theta^{l}_{13}$ | 0.0221697 | $\sin^{2}\theta_{13}$ | 0.0221678 | $\sin^{2}\theta^{l}_{13}$ | 0.0221678 | 0.221678
$\sin^{2}\theta^{l}_{23}$ | 0.596554 | $\sin^{2}\theta_{23}$ | 0.477773 | $\sin^{2}\theta^{l}_{23}$ | 0.477732 | 0.477732
$\delta^{l}_{CP}/\pi$ | 1.36268 | $\delta^{l}_{CP}/\pi$ | 1.56812 | $\delta^{l}_{CP}/\pi$ | 1.56839 | 1.56839
$m_{1}/\mathrm{meV}$ | 25.4036 | $m_{1}/\mathrm{meV}$ | 29.1666 | $m_{1}/\mathrm{meV}$ | 29.1665 | 29.1665
$m_{2}/\mathrm{meV}$ | 26.8187 | $m_{2}/\mathrm{meV}$ | 30.4071 | $m_{2}/\mathrm{meV}$ | 30.4070 | 30.4070
$m_{3}/\mathrm{meV}$ | 56.1239 | $m_{3}/\mathrm{meV}$ | 58.0277 | $m_{3}/\mathrm{meV}$ | 58.0364 | 58.0364
$m_{\beta}/\mathrm{meV}$ | 26.8833 | $m_{\beta}/\mathrm{meV}$ | 30.4744 | $m_{\beta}/\mathrm{meV}$ | 30.4746 | 30.4745
$\chi^{2}_{\text{min}}$ | 2.29 | $\chi^{2}_{\text{min}}$ | 22.77 | $\chi^{2}_{\text{min}}$ | 27.79 | 22.79

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_{\beta}$ probed by
direct kinematic search in beta decay. Note that the transformation
$\tau\rightarrow-\tau^{*}$ leaves all observables unchanged except shifting
the sign of the CP phase $\delta^{l}_{CP}$.

###  4.2 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 $\delta^{l}_{CP}$, $\alpha_{21}$ and $\alpha_{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\pm 0.011288)\times 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 $\Lambda$ 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 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
$\varepsilon_{1\ell}$, $\varepsilon_{1\tilde{\ell}}$,
$\varepsilon_{\tilde{1}\ell}$ and $\varepsilon_{\tilde{1}{\ell}}$ respectively
which are defined by [91, 92]

$\displaystyle\varepsilon_{1\ell}\equiv\frac{\Gamma_{N_{1}\ell}-\Gamma_{N_{1}\bar{\ell}}}{\Gamma_{N_{1}\ell}+\Gamma_{N_{1}\bar{\ell}}},~{}~{}~{}\quad\varepsilon_{1\tilde{\ell}}\equiv\frac{\Gamma_{N_{1}\widetilde{\ell}}-\Gamma_{N_{1}\widetilde{\ell}^{*}}}{\Gamma_{N_{1}\widetilde{\ell}}+\Gamma_{N_{1}\widetilde{\ell}^{*}}},$
(59)
$\displaystyle\varepsilon_{\tilde{1}\ell}\equiv\frac{\Gamma_{\widetilde{N}_{1}^{*}\ell}-\Gamma_{\widetilde{N}_{1}\bar{\ell}}}{\Gamma_{\widetilde{N}_{1}^{*}\ell}+\Gamma_{\widetilde{N}_{1}\bar{\ell}}},~{}~{}~{}\quad\varepsilon_{\tilde{1}\tilde{\ell}}\equiv\frac{\Gamma_{\widetilde{N}_{1}\tilde{\ell}}-\Gamma_{\widetilde{N}_{1}^{*}\tilde{\ell^{*}}}}{\Gamma_{\widetilde{N}_{1}\tilde{\ell}}+\Gamma_{\widetilde{N}_{1}^{*}\tilde{\ell^{*}}}}.$
(60)

In the minimal supersymmetric standard model, all the above decay asymmetries
are equal
$\varepsilon_{1\ell}=\varepsilon_{1\widetilde{\ell}}=\varepsilon_{\widetilde{1}\ell}=\varepsilon_{\widetilde{1}\widetilde{\ell}}$
[91, 92]. In the basis where the Majorana mass matrix of the right-handed
neutrinos is diagonal and real, the lepton asymmetry parameter
$\varepsilon_{1\ell}$ is given by [91, 93]

$\varepsilon_{1\ell}=\frac{1}{8\pi{\left(\lambda_{\nu}\lambda_{\nu}^{\dagger}\right)}_{11}}\sum_{i=2,3}\texttt{Im}\left\\{\left[{\left(\lambda_{\nu}\lambda_{\nu}^{\dagger}\right)}_{1i}\right]^{2}\right\\}~{}g\left(\frac{M_{i}^{2}}{M_{1}^{2}}\right)\,,$
(61)

where $\lambda_{\nu}$ is the neutrino Yukawa coupling matrix in the convention
$(\lambda_{\nu})_{ij}N^{c}_{i}(L_{j}\cdot H_{u})$ and the loop function $g$
reads

$\displaystyle
g(x)=\sqrt{x}\left[\frac{2}{1-x}-\mathrm{ln}\left(\frac{1+x}{x}\right)\right]\stackrel{{\scriptstyle
x\gg 1}}{{\longrightarrow}}-\frac{3}{\sqrt{x}}\,.$ (62)

Non-vanishing asymmetry parameter $\varepsilon_{1\ell}$ requires that the off-
diagonal entries of the product $\lambda_{\nu}\lambda_{\nu}^{\dagger}$ are
complex and different from zero. For the models L1, L4, L5, L6, the product
$\lambda_{\nu}\lambda_{\nu}^{\dagger}$ is proportional to the unity matrix.
Consequently the lepton asymmetry $\varepsilon_{1\ell}$ 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
$\hat{Y}_{\Delta}$ as follows

$\displaystyle Y_{B}=\frac{10}{31}\hat{Y}_{\Delta}\,,$ (63)

The $B-L$ asymmetry $\hat{Y}_{\Delta}$ can be computed by solving the
following Boltzmann equations [92, 94]

$\displaystyle\frac{\mathrm{d}Y_{N_{1}}}{\mathrm{d}z}$ $\displaystyle=$
$\displaystyle-\,2K\,z\,\frac{K_{1}(z)}{K_{2}(z)}\,f_{1}(z)\,(Y_{N_{1}}-Y^{\mathrm{eq}}_{N_{1}})\;,$
$\displaystyle\frac{\mathrm{d}Y_{\widetilde{N}_{1}}}{\mathrm{d}z}$
$\displaystyle=$
$\displaystyle-\,2K\,z\,\frac{K_{1}(z)}{K_{2}(z)}\,f_{1}(z)\,(Y_{\widetilde{N}_{1}}-Y^{\mathrm{eq}}_{\widetilde{N}_{1}})\;,$
$\displaystyle\frac{\mathrm{d}\hat{Y}_{\Delta}}{\mathrm{d}z}$ $\displaystyle=$
$\displaystyle-\;(\varepsilon_{1\ell}+\varepsilon_{1\widetilde{\ell}})\,K\,z\,\frac{K_{1}(z)}{K_{2}(z)}\,f_{1}(z)\,(Y_{N_{1}}-Y^{\mathrm{eq}}_{N_{1}})$
(64)
$\displaystyle-\;(\varepsilon_{\widetilde{1}\ell}+\varepsilon_{\widetilde{1}\widetilde{\ell}})\,K\,z\,\frac{K_{1}(z)}{K_{2}(z)}\,f_{1}(z)\,(Y_{\widetilde{N}_{1}}-Y^{\mathrm{eq}}_{\widetilde{N}_{1}})$
$\displaystyle-\;K\,z\,\frac{K_{1}(z)}{K_{2}(z)}\,f_{2}(z)\,\frac{\hat{Y}_{\Delta}}{\hat{Y}^{\mathrm{eq}}_{\ell}}\,(Y^{\mathrm{eq}}_{N_{1}}+Y^{\mathrm{eq}}_{\widetilde{N}_{1}})\;,$

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_{\widetilde{N}_{1}}$ denotes the density of the lightest right-handed
neutrino $N_{1}$ with mass $M_{1}$ and its supersymmetric partners
$\widetilde{N}_{1}$. The notations $Y^{\mathrm{eq}}_{N_{1}}$,
$Y^{\mathrm{eq}}_{\widetilde{N}_{1}}$ and $\hat{Y}^{\mathrm{eq}}_{\ell}$ are
corresponding equilibrium number densities and they take the following form

$\hat{Y}_{\ell}^{\mathrm{eq}}=Y^{\mathrm{eq}}_{\tilde{\ell}}+Y^{eq}_{\ell},\quad
Y^{\mathrm{eq}}_{\tilde{\ell}}\simeq
Y^{eq}_{\ell}\simeq\frac{45}{\pi^{4}g_{*}},\quad
Y_{N_{1}}^{\mathrm{eq}}(z)=Y_{\widetilde{N}_{1}}^{\mathrm{eq}}(z)\simeq\frac{45}{2\pi^{4}g_{*}}z^{2}K_{2}(z)\,,$
(65)

with $g_{*}=228.75$ being the number of degrees of freedom. Moreover, the
washout parameters $K_{\alpha}$ and $K$ are defined as

$\displaystyle
K_{\alpha}=\frac{\tilde{m}_{1\alpha}}{m^{*}},\qquad\tilde{m}_{1\alpha}\equiv\frac{|\lambda_{1\alpha}|^{2}v_{u}^{2}}{M_{1}},\qquad\,K=\sum_{\alpha}K_{\alpha}\,,$
(66)

where

$\displaystyle v_{u}=v\sin\beta,\qquad\quad m^{*}\simeq\sin^{2}\beta\times
1.58\times 10^{-3}\,\text{eV}\,.$ (67)

At the best fit point of model L9, we find the value of $K=1.293\times
10^{6}\gg 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]

$f_{1}(z)=2f_{2}(z)=\left[\frac{K_{s}}{zK}+\frac{z}{t}\ln\left(1+\frac{t}{z}\right)\right]\frac{K_{2}(z)}{K_{1}(z)}\,,$
(68)

with

$t=\frac{K}{K_{s}\ln(M_{1}/M_{h})},\qquad\quad\frac{K_{s}}{K}=\frac{9}{8\pi^{2}}\,.$
(69)

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 $(\alpha^{D}v_{u})^{2}/(\beta^{N}\Lambda)$ 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:

$\displaystyle\Lambda=3.36\times 10^{15}\mathrm{GeV}.$ (70)

Accordingly the right-handed neutrino masses are determined to be $M_{1}\simeq
1.985\times 10^{15}$ GeV, $M_{2}\simeq 6.723\times 10^{15}$ GeV and
$M_{3}\simeq 6.834\times 10^{15}$ GeV. The VEV of $\tau$ 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
$\langle\tau\rangle$. We plot the contour region of $Y_{B}$ in the plane
Im$\langle\tau\rangle$ versus Re$\langle\tau\rangle$, where we fix all the
coupling constants at their best fit values and $\Lambda$ at the value in Eq.
(70). The green area indicates the $3\sigma$ 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 $\tau$ 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 $\mathcal{D}$ and the imaginary axis,
certain residual gCP symmetry is preserved such that the lepton asymmetry
$\varepsilon_{1\ell}$ vanishes and no matter-antimatter asymmetry can be
generated. Hence the VEV $\langle\tau\rangle$ 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 $\delta^{l}_{CP}$ in figure 3. Here the modulus vacuum expectation value
$\langle\tau\rangle$ 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\sigma$ ranges [76]. Imposing the observed baryon asymmetry
$Y_{B}$, the allowed region of $\delta^{l}_{CP}$ would shrink considerably and
it is around $1.28\pi$.

Figure 3: The contour plot $Y_{B}$ in the
Re$\langle\tau\rangle$-Im$\langle\tau\rangle$ plane. The green region
represent the region for which both lepton masses and lepton mixing angles are
compatible with experimental data at the $3\sigma$ level or better.  Figure 4:
The correlation between $Y_{B}$ and $\delta_{CP}$. The horizonal light blue
band denotes the $95\%$CL region of $Y_{B}$, and the vertical light green band
represents the 3$\sigma$ range of $\delta^{l}_{CP}$.

###  4.3 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 $\chi^{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 $\theta^{q}_{12}$, $\theta^{q}_{13}$, $\theta^{q}_{23}$, and a CP
violation phase $\delta^{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$\sim$Q10 make use of seven real coupling constants in addition to the
modulus $\tau$, 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 left-handed 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 $\chi^{2}_{\text{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.

Models | Q1 | Q2 | Q3 | Q4 | Q5
---|---|---|---|---|---
$\rho_{Q}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$
$\rho_{u^{c}}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}^{\prime}$ | $\mathbf{2}\oplus\mathbf{1}^{\prime}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}^{\prime}$
$\rho_{d^{c}}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{1}\oplus\mathbf{1}^{\prime}\oplus\mathbf{1}$ | $\mathbf{1}\oplus\mathbf{1}\oplus\mathbf{1}$
$k_{Q}$ | $k$ | $k$ | $k$ | $k$ | $k$
$k_{u^{c}}$ | $6-k,4-k$ | $6-k,4-k$ | $6-k,6-k$ | $2-k,6-k$ | $2-k,8-k$
$k_{d^{c}}$ | $2-k,6-k$ | $2-k,6-k$ | $2-k,6-k$ | $2-k,4-k,8-k$ | $2-k,4-k,6-k$
Models | Q6 | Q7 | Q8 | Q9 | Q10
$\rho_{Q}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}$
$\rho_{u^{c}}$ | $\mathbf{2}\oplus\mathbf{1}^{\prime}$ | $\mathbf{2}\oplus\mathbf{1}^{\prime}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}$ | $\mathbf{2}\oplus\mathbf{1}$
$\rho_{d^{c}}$ | $\mathbf{1}\oplus\mathbf{1}\oplus\mathbf{1}$ | $\mathbf{1}^{\prime}\oplus\mathbf{1}\oplus\mathbf{1}$ | $\mathbf{1}^{\prime}\oplus\mathbf{1}^{\prime}\oplus\mathbf{1}$ | $\mathbf{1}^{\prime}\oplus\mathbf{1}\oplus\mathbf{1}$ | $\mathbf{1}^{\prime}\oplus\mathbf{1}^{\prime}\oplus\mathbf{1}$
$k_{Q}$ | $k,k$ | $k,k$ | $k,k$ | $k,k$ | $k,k$
$k_{u^{c}}$ | $-k,6-k$ | $-k,6-k$ | $-k,4-k$ | $-k,4-k$ | $-k,4-k$
$k_{d^{c}}$ | $4-k,2-k,-k$ | $4-k,4-k,2-k$ | $6-k,4-k,2-k$ | $6-k,2-k,-k$ | $6-k,2-k,-k$

Table 7:  Classification of quark fields in the minimal quark models with
$S_{4}$ modular symmetry and gCP symmetry, where $k$ can be any integer.
Notice that the Higgs fields are invariant under $S_{4}$ with zero modular
weight.

Models | Q1 | Q2 | Q3 | Models | Q4 | Q5
---|---|---|---|---|---|---
$\mathrm{Re}\langle\tau\rangle$ | $-0.436841$ | $-0.437014$ | $-0.436195$ | $\mathrm{Re}\langle\tau\rangle$ | $-0.00687942$ | 0.493588
$\mathrm{Im}\langle\tau\rangle$ | 1.81494 | 1.81474 | 1.81557 | $\mathrm{Im}\langle\tau\rangle$ | 1.00472 | 0.874580
$\beta^{u}/\alpha_{1}^{u}$ | 0.00343099 | 0.000237927 | 0.000235881 | $\beta_{1}^{u}/\alpha^{u}$ | 0.479893 | 1505.91
$\beta_{1}^{d}/\alpha^{d}$ | 0.304628 | 0.304084 | 0.304619 | $\beta^{d}/\alpha^{d}$ | 680.772 | 612.613
$\alpha_{2}^{u}/\alpha_{1}^{u}$ | 1.03883 | 1.03878 | 1.03904 | $\gamma^{d}_{1}/\alpha^{d}$ | 3.32280 | 34.3723
$\alpha^{\prime}{}^{u}/\alpha_{1}^{u}$ | 1.00004 | 1.00006 | 1.00000 | $\beta_{2}^{u}/\alpha_{1}^{u}$ | $-226.241$ | $-0.186995$
$\beta_{2}^{d}/\alpha^{d}$ | 10.9113 | 10.8935 | 10.9002 | $\gamma^{d}_{2}/\alpha^{d}$ | 39.7179 | $-5.96481$
$\alpha_{1}^{u}v_{u}(\mathrm{GeV})$ | 8.84642 | 8.84655 | 8.84579 | $\alpha^{u}v_{u}(\mathrm{GeV})$ | 0.0808452 | 0.0865103
$\alpha^{d}v_{d}(\mathrm{GeV})$ | 0.0236747 | 0.0237130 | 0.0236988 | $\alpha^{d}v_{d}(\mathrm{GeV})$ | 0.000336231 | 0.000279067
$\theta^{q}_{12}$ | 0.227433 | 0.227325 | 0.227439 | $\theta^{q}_{12}$ | 0.227402 | 0.227351
$\theta^{q}_{13}$ | 0.00338504 | 0.00337795 | 0.00340799 | $\theta^{q}_{13}$ | 0.00318930 | 0.00310537
$\theta^{q}_{23}$ | 0.0388938 | 0.0389309 | 0.0387763 | $\theta^{q}_{23}$ | 0.0386561 | 0.0399389
$\delta^{q}_{CP}/{}^{\circ}$ | 69.4363 | 69.5128 | 69.4400 | $\delta^{q}_{CP}/{}^{\circ}$ | 69.4280 | 70.1150
$m_{u}/m_{c}$ | 0.00192985 | 0.00192161 | 0.00192901 | $m_{u}/m_{c}$ | 0.00260041 | 0.00310718
$m_{c}/m_{t}$ | 0.00273544 | 0.00273426 | 0.00274461 | $m_{c}/m_{t}$ | 0.00265548 | 0.00297164
$m_{d}/m_{s}$ | 0.0458926 | 0.0458843 | 0.0459400 | $m_{d}/m_{s}$ | 0.0504604 | 0.0507499
$m_{s}/m_{b}$ | 0.0176515 | 0.0176858 | 0.0176518 | $m_{s}/m_{b}$ | 0.0177025 | 0.0176849
$\chi^{2}_{\text{min}}$ | 0.58521 | 0.588937 | 0.663495 | $\chi^{2}_{\text{min}}$ | 4.26505 | 16.2615
Models | Q7 | Q8 | Q10 | Models | Q6 | Q9
$\mathrm{Re}\langle\tau\rangle$ | $-0.495895$ | 0.307358 | 0.307719 | $\mathrm{Re}\langle\tau\rangle$ | $-0.495886$ | $-0.191783$
$\mathrm{Im}\langle\tau\rangle$ | 0.875601 | 2.21966 | 2.21896 | $\mathrm{Im}\langle\tau\rangle$ | 0.875587 | 2.21983
$\delta^{u}/\alpha^{u}$ | 0.218871 | 0.607105 | 0.606631 | $\delta^{u}/\alpha^{u}$ | 0.219074 | 0.612005
$\gamma^{u}/\alpha^{u}$ | 1687.47 | 258.331 | 257.759 | $\gamma^{u}/\alpha^{u}$ | 1692.88 | 258.949
$\beta^{d}/\alpha^{d}$ | 4.33152 | 0.461192 | 0.455737 | $\beta^{d}/\alpha^{d}$ | 102.878 | 116.077
$\gamma^{d}/\alpha^{d}$ | 457.852 | 119.140 | 119.036 | $\delta^{d}_{3}/\alpha^{d}$ | 0.167218 | 0.150554
$\delta^{d}/\alpha^{d}$ | 27.0580 | 8.82734 | 8.79643 | $\delta^{d}_{1}/\alpha^{d}$ | 5.93625 | -8.38832
$\alpha^{u}v_{u}(\mathrm{GeV})$ | 0.244549 | 0.244977 | 0.245520 | $\alpha^{u}v_{u}(\mathrm{GeV})$ | 0.243972 | 0.244414
$\alpha^{d}v_{d}(\mathrm{GeV})$ | 0.00114978 | 0.00554085 | 0.00554569 | $\alpha^{d}v_{d}(\mathrm{GeV})$ | 0.00511694 | 0.00568697
$\theta^{q}_{12}$ | 0.227357 | 0.227366 | 0.227330 | $\theta^{q}_{12}$ | 0.227368 | 0.227365
$\theta^{q}_{13}$ | 0.00333108 | 0.00332068 | 0.00332544 | $\theta^{q}_{13}$ | 0.00332646 | 0.00333852
$\theta^{q}_{23}$ | 0.0389070 | 0.0389306 | 0.0390165 | $\theta^{q}_{23}$ | 0.0388726 | 0.0389099
$\delta^{q}_{CP}/{}^{\circ}$ | 69.2159 | 69.3956 | 69.2651 | $\delta^{q}_{CP}/{}^{\circ}$ | 69.4425 | 69.0265
$m_{u}/m_{c}$ | 0.00333311 | 0.00332201 | 0.00332677 | $m_{u}/m_{c}$ | 0.00332849 | 0.00334163
$m_{c}/m_{t}$ | 0.00273137 | 0.00273614 | 0.00274221 | $m_{c}/m_{t}$ | 0.00272492 | 0.00272986
$m_{d}/m_{s}$ | 0.0505789 | 0.0505368 | 0.0499235 | $m_{d}/m_{s}$ | 0.0504933 | 0.0505272
$m_{s}/m_{b}$ | 0.0176638 | 0.0176765 | 0.0176811 | $m_{s}/m_{b}$ | 0.0176839 | 0.0176835
$\chi^{2}_{\text{min}}$ | 5.64026 | 5.64895 | 5.70768 | $\chi^{2}_{\text{min}}$ | 5.63916 | 5.65984

Table 8:  Results of fit for the quark models listed in table 7.

###  4.4 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 $\tau$, and all
the CP violation phases in lepton and quark sectors arise from the modulus VEV
$\langle\tau\rangle$. The quark-lepton unification has been studied in the
context of $A_{4}$ [11, 14, 28, 31], $T^{\prime}$ [45], $S^{\prime}_{4}$ [47]
and $A^{\prime}_{5}$ [49] modular symmetries. The most predictive model
contains fifteen parameters including the real and imaginary part of the
modulus $\tau$ [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}\times 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 $\text{Re}\tau$ and $\text{Im}\tau$ to
explain twenty-two observables: six quark masses $m_{u,c,t,d,s,b}$, three
quark mixing angles $\theta^{q}_{12,13,23}$, a quark CP violation phase
$\delta^{q}_{CP}$, three charged lepton masses $m_{e,\mu,\tau}$, three
neutrino masses $m_{1,2,3}$, three lepton mixing angles
$\theta^{l}_{12,13,23}$ and three leptonic CP violation phases
$\delta^{l}_{CP}$, $\alpha_{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 $\mathbf{3}$, while all the other lepton fields $L$,
$E^{c}$ and quark fields $Q$, $u^{c}$, $d^{c}$ are assigned to
$\mathbf{2}\oplus\mathbf{1}$. The transformation properties of the matter
fields are given by

$\displaystyle\rho_{u^{c}}=\rho_{d^{c}}=\rho_{Q}=\rho_{E^{c}}=\rho_{L}=\mathbf{2}\oplus\mathbf{1}\,,\,\rho_{N}^{c}=\mathbf{3}\,,$
$\displaystyle
k_{u^{c}_{D}}=k_{d^{c}_{D}}=2-k\,,\,k_{u^{c}_{3}}=-k\,,\,k_{d^{c}_{3}}=8-k\,,\,k_{Q_{D}}=k_{Q_{3}}=k\,,\,$
$\displaystyle k_{E^{c}_{D}}=k_{E^{c}_{3}}=k_{L_{D}}=k_{L_{3}}=k_{N^{c}}=1\,,$
(71)

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,

$\displaystyle M_{e}$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\alpha^{e}Y_{\mathbf{2},1}^{(2)}&0&\beta^{e}Y_{\mathbf{2},2}^{(2)}\\\
0&\alpha^{e}Y_{\mathbf{2},2}^{(2)}&\beta^{e}Y_{\mathbf{2},1}^{(2)}\\\
\gamma^{e}Y_{\mathbf{2},2}^{(2)}&\gamma^{e}Y_{\mathbf{2},1}^{(2)}&0\end{matrix}\right)v_{d}\,,~{}~{}~{}~{}M_{D}=\left(\begin{matrix}\alpha^{D}Y_{\mathbf{3},2}^{(2)}&\alpha^{D}Y_{\mathbf{3},3}^{(2)}&\beta^{D}Y_{\mathbf{3},1}^{(2)}\\\
\alpha^{D}Y_{\mathbf{3},1}^{(2)}&\alpha^{D}Y_{\mathbf{3},2}^{(2)}&\beta^{D}Y_{\mathbf{3},3}^{(2)}\\\
\alpha^{D}Y_{\mathbf{3},3}^{(2)}&\alpha^{D}Y_{\mathbf{3},1}^{(2)}&\beta^{D}Y_{\mathbf{3},2}^{(2)}\end{matrix}\right)v_{u}\,,$
$\displaystyle M_{N^{c}}$ $\displaystyle=$
$\displaystyle\beta^{N}\left(\begin{matrix}0&Y_{\mathbf{2},1}^{(2)}&Y_{\mathbf{2},2}^{(2)}\\\
Y_{\mathbf{2},1}^{(2)}&Y_{\mathbf{2},2}^{(2)}&0\\\
Y_{\mathbf{2},2}^{(2)}&0&Y_{\mathbf{2},1}^{(2)}\end{matrix}\right)\Lambda\,,\,~{}~{}~{}~{}M_{u}=\left(\begin{matrix}\alpha^{u}Y_{\mathbf{2},1}^{(2)}&0&\beta^{u}Y_{\mathbf{2},2}^{(2)}\\\
0&\alpha^{u}Y_{\mathbf{2},2}^{(2)}&\beta^{u}Y_{\mathbf{2},1}^{(2)}\\\
0&0&\delta^{u}\end{matrix}\right)v_{u}\,,\,$ $\displaystyle M_{d}$
$\displaystyle=$
$\displaystyle\left(\begin{matrix}\alpha^{d}Y_{\mathbf{2},1}^{(2)}&0&\beta^{d}Y_{\mathbf{2},2}^{(2)}\\\
0&\alpha^{d}Y_{\mathbf{2},2}^{(2)}&\beta^{d}Y_{\mathbf{2},1}^{(2)}\\\
\gamma^{d}_{1}Y_{\mathbf{2}I,2}^{(8)}+\gamma^{d}_{2}Y_{\mathbf{2}II,2}^{(8)}&\gamma^{d}_{1}Y_{\mathbf{2}I,1}^{(8)}+\gamma^{d}_{2}Y_{\mathbf{2}II,1}^{(8)}&\delta^{d}Y_{\mathbf{1}}^{(8)}\end{matrix}\right)v_{d}\,.$
(72)

The best fit values of the input parameters for this unified model is
determined to be

$\displaystyle\mathrm{Re}\langle\tau\rangle=-0.477059\,,\,~{}\mathrm{Im}\langle\tau\rangle=1.28145\,,\,~{}\beta^{u}/\alpha^{u}=285.975\,,\,~{}\delta^{u}/\alpha^{u}=1.09959\,,$
$\displaystyle\gamma^{d}_{1}/\alpha^{d}=2.52377\,,\,~{}\beta^{d}/\alpha^{d}=-94.0867\,,\,~{}\gamma^{d}_{2}/\alpha^{d}=6.99468\,,\,\delta^{d}/\alpha^{d}=-83.5316\,,$
$\displaystyle\beta^{e}/\alpha^{e}=2820.85\,,\,~{}\gamma^{e}/\alpha^{e}=159.393\,,\,\beta^{D}/\alpha^{D}=1.07257\,,\,~{}\alpha^{u}v_{u}=0.216464\,\mathrm{GeV}\,,$
$\displaystyle\alpha^{d}v_{d}=0.00294862\,\mathrm{GeV}\,,~{}\alpha^{e}v_{d}=0.461746\,\mathrm{MeV}\,,\,~{}\dfrac{(\alpha^{D}v_{u})^{2}}{\beta^{N}\Lambda}=11.9256\,\mathrm{meV}\,.$
(73)

The masses and mixing parameters of leptons and quarks are predicted to be

$\displaystyle\sin^{2}\theta^{l}_{12}=0.339580\,,\,~{}~{}\sin^{2}\theta^{l}_{13}=0.0215912\,,\,~{}~{}\sin^{2}\theta^{l}_{23}=0.615854\,,$
$\displaystyle\delta^{l}_{CP}/\pi=1.33145\,,\,~{}~{}\alpha_{21}/\pi=1.29710\,,\,~{}~{}\alpha_{31}/\pi=0.476112\,,$
$\displaystyle
m_{1}=29.1179\,\mathrm{meV}\,,~{}~{}\,m_{2}=30.3604\,\mathrm{meV}\,,\,~{}~{}m_{3}=57.1162\,\mathrm{meV}\,,$
$\displaystyle
m_{e}/m_{\mu}=0.00480040\,,\,~{}~{}m_{\mu}/m_{\tau}=0.0565053\,,\,~{}~{}m_{\beta\beta}=16.3936\,\mathrm{meV}\,,$
$\displaystyle
m_{\beta}=30.4062\,\mathrm{meV}\,,\,~{}~{}m_{u}/m_{c}=0.00192854\,,\,~{}~{}m_{c}/m_{t}=0.00272947\,,$
$\displaystyle
m_{d}/m_{s}=0.0459963\,,\,~{}~{}m_{s}/m_{b}=0.0178060\,,\,~{}~{}\theta^{q}_{12}=0.227380\,,$
$\displaystyle\theta^{q}_{13}=0.00311374\,,\,~{}~{}\theta^{q}_{23}=0.0394220\,,\,~{}~{}\delta^{q}_{CP}=68.6949^{\circ}\,,$
(74)

which are compatible with experimental data at $3\sigma$ level. All the
coupling constants as well as complex modulus $\tau$ are treated as random
numbers, and the $3\sigma$ bounds of the mass ratios and mixing angles of both
quarks and leptons are imposed. The values of $\tau$ 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 $\tau$ 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 $\tau$ 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.

Figure 5: The experimentally favored values of $\tau$ are displayed in the
first row in the quark-lepton unification model. The quark (lepton) masses and
mixing parameters are compatible with experimental data at $3\sigma$ level or
better in the light blue (red) area, and the common values of $\tau$ are
indicated with black. The second and the third rows are for the correlation
among the neutrino masses and mixing parameters.

##  5 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, 34, 35, 36, 37,
38, 39, 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 $\mathbf{1}$,
$\mathbf{1}^{\prime}$ and triplet representations $\mathbf{3}$ and
$\mathbf{3}^{\prime}$, the $S_{4}$ group has a doublet irreducible
representation $\mathbf{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 $\tau$ 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 $\mathbf{3}$ under
$S_{4}$, the light neutrino mass matrix only depends on the modulus $\tau$ 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 $\tau=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 $\tau$. The allowed range of the Dirac
CP phase $\delta^{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 $\tau$ 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 $\tau$ is dominantly fixed by the experimental data of lepton masses
and mixing angles. The allowed range of $\tau$ 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.

### Acknowledgements

This work is supported by the National Natural Science Foundation of China
under Grant Nos 11975224, 11835013 and 11947301.

## Appendix

##  Appendix A $S_{4}$ group

The finite modular group $\Gamma_{4}$ is isomorphic to the permutation group
$S_{4}$ which is the group of all permutations of four elements. Geometrically
$S_{4}$ is the group of orientation-preserving symmetries of the cube or
equivalently the octahedron. As shown in Eq. (5), the inhomogeneous finite
modular group $\Gamma_{4}\cong S_{4}$ can be generated by the modular
generators $S$ and $T$ satisfying the relations

$S^{2}=(ST)^{3}=T^{4}=\mathbb{1}\,.$ (A.1)

In the paradigm of traditional flavor symmetry, it is convenient to express
the $S_{4}$ group in terms of three generators $\hat{S}$, $\hat{T}$ and
$\hat{U}$ obeying the multiplication rules [100]:

$\displaystyle\hat{S}^{2}=\hat{T}^{3}=\hat{U}^{2}=(\hat{S}\hat{T})^{3}=(\hat{S}\hat{U})^{2}=(\hat{T}\hat{U})^{2}=(\hat{S}\hat{T}\hat{U})^{4}=1\,.$
(A.2)

The generators $\hat{S}$ and $\hat{T}$ alone generate the group $A_{4}$, and
the generators $\hat{T}$ and $\hat{U}$ alone generate the group $S_{3}$. The
two different choices of generators are related as follows

$S=\hat{S}\hat{U},\,T=\hat{S}\hat{T}^{2}\hat{U},\,ST=\hat{T}\,,~{}~{}\text{or}~{}~{}\hat{S}=(ST^{2})^{2},\,\hat{T}=ST,\,\hat{U}=T^{2}ST^{2}\,.$
(A.3)

The $S_{4}$ group has 24 elements and 5 irreducible representations including
two singlet representations $\mathbf{1}$, $\mathbf{1}^{\prime}$, a doublet
representation $\mathbf{2}$ and two triplet representations $\mathbf{3}$ and
$\mathbf{3}^{\prime}$. In this work, we choose the same representation basis
as that of [100], i.e. the representation matrix of the generator $\hat{T}$ is
diagonal. We summarize the representation matrices of the generators in table
9.

|  $\rho_{\mathbf{r}}(\hat{S}$) | $\rho_{\mathbf{r}}(\hat{T})$ | $\rho_{\mathbf{r}}(\hat{U})$ | $\rho_{\mathbf{r}}(S)$ | $\rho_{\mathbf{r}}(T)$
---|---|---|---|---|---
${\bf 1}$, ${\bf 1^{\prime}}$ | 1 | 1 | $\pm 1$ | $\pm 1$ | $\pm 1$
${\bf 2}$ | $\left(\begin{array}[]{cc}1&~{}0\\\ 0&~{}1\end{array}\right)$ | $\left(\begin{array}[]{cc}\omega&~{}0\\\ 0&~{}\omega^{2}\end{array}\right)$ | $\left(\begin{array}[]{cc}0&~{}1\\\ 1&~{}0\end{array}\right)$ | $\left(\begin{array}[]{cc}0&~{}1\\\ 1&~{}0\end{array}\right)$ | $\left(\begin{array}[]{cc}0{}&~{}\omega^{2}\\\ \omega{}&~{}0\end{array}\right)$
${\bf 3}$, ${\bf 3^{\prime}}$ | $\dfrac{1}{3}\left(\begin{array}[]{ccc}-1&~{}2{}&2\\\ 2&~{}-1{}&2\\\ 2&~{}2{}&-1\end{array}\right)$ | $\left(\begin{array}[]{ccc}1&~{}0{}&0\\\ 0&~{}\omega^{2}{}&0\\\ 0&~{}0{}&\omega\end{array}\right)$ | $\mp\left(\begin{array}[]{ccc}1&~{}0{}&0\\\ 0&~{}0{}&1\\\ 0&~{}1{}&0\end{array}\right)$ | $\pm\dfrac{1}{3}\left(\begin{array}[]{ccc}1&~{}-2{}&-2\\\ -2&~{}-2{}&1\\\ -2&~{}1{}&-2\end{array}\right)$ | $\pm\dfrac{1}{3}\left(\begin{array}[]{ccc}1&~{}-2\omega^{2}{}&-2\omega\\\ -2&~{}-2\omega^{2}{}&\omega\\\ -2&~{}\omega^{2}{}&-2\omega\end{array}\right)$

Table 9: The representation matrices of the generators
$\hat{S}$,$\hat{T}$,$\hat{U}$ as well as $S,T$ in the five irreducible
representations of $S_{4}$ in the chosen basis, where $\omega=e^{2\pi i/3}$ is
the cube root of unit.

We now list the Kronecker products and the corresponding Clebsch-Gordan
coefficients which are quite useful when constructing $S_{4}$ modular
invariant models. For convenience, we denote $\mathbf{1}\equiv\mathbf{1}^{0}$,
$\mathbf{1}^{\prime}\equiv\mathbf{1}^{1}$, $\mathbf{3}\equiv\mathbf{3}^{0}$,
$\mathbf{3}^{\prime}\equiv\mathbf{3}^{1}$ for the singlet or triplet
representations. We shall use $\alpha_{i}$ to denote the elements of first
representation and $\beta_{i}$ stands for the elements of the second
representation of the tensor product $\mathbf{R_{1}}\otimes\mathbf{R_{2}}$,
where $\mathbf{R_{1}}$ and $\mathbf{R_{2}}$ are two irreducible
representations of $S_{4}$.

$\displaystyle{}\mathbf{1}^{i}\otimes\mathbf{1}^{j}=\mathbf{1}^{<i+j>}\sim\alpha\beta\,,$
(A.4)
$\displaystyle\mathbf{1}^{i}\otimes\mathbf{2}=\mathbf{2}\sim\left(\begin{matrix}\alpha\beta_{1}\\\
(-1)^{i}\alpha\beta_{2}\end{matrix}\right)\,,$
$\displaystyle\mathbf{1}^{i}\otimes\mathbf{3}^{j}=\mathbf{3}^{<i+j>}\sim\left(\begin{matrix}\alpha\beta_{1}\\\
\alpha\beta_{2}\\\ \alpha\beta_{3}\end{matrix}\right)\,,$
$\displaystyle\mathbf{2}\otimes\mathbf{2}=\mathbf{1}^{0}\oplus\mathbf{1}^{1}\oplus\mathbf{2}\,,\quad\begin{cases}\mathbf{1}^{l}\sim\alpha_{1}\beta_{2}+(-1)^{l}\alpha_{2}\beta_{1}\\\
\mathbf{2}\sim\left(\begin{matrix}\alpha_{2}\beta_{2}\\\
\alpha_{1}\beta_{1}\end{matrix}\right)\end{cases}\,,$
$\displaystyle\mathbf{2}\otimes\mathbf{3}^{i}=\mathbf{3}^{0}\oplus\mathbf{3}^{1}\,,\quad\mathbf{3}^{l}\sim\left(\begin{matrix}\alpha_{1}\beta_{2}\\\
\alpha_{1}\beta_{3}\\\
\alpha_{1}\beta_{1}\end{matrix}\right)+(-1)^{i+l}\left(\begin{matrix}\alpha_{2}\beta_{3}\\\
\alpha_{2}\beta_{1}\\\ \alpha_{2}\beta_{2}\end{matrix}\right)\,,$
$\displaystyle\mathbf{3}^{i}\otimes\mathbf{3}^{j}=\mathbf{1}^{<i+j>}\oplus\mathbf{2}\oplus\mathbf{3}^{0}\oplus\mathbf{3}^{1}\,,\quad\begin{cases}\mathbf{1}^{<i+j>}\sim\alpha_{1}\beta_{1}+\alpha_{2}\beta_{3}+\alpha_{3}\beta_{2}\\\
\mathbf{2}\sim\left(\begin{matrix}\alpha_{1}\beta_{3}+\alpha_{2}\beta_{2}+\alpha_{3}\beta_{1}\\\
(-1)^{i+j}(\alpha_{1}\beta_{2}+\alpha_{2}\beta_{1}+\alpha_{3}\beta_{3})\end{matrix}\right)\\\
\mathbf{3}^{l}\sim\left(\begin{matrix}\alpha_{1}\beta_{1}-\alpha_{3}\beta_{2}\\\
\alpha_{3}\beta_{3}-\alpha_{2}\beta_{1}\\\
\alpha_{2}\beta_{2}-\alpha_{1}\beta_{3}\end{matrix}\right)-(-1)^{i+j+l}\left(\begin{matrix}\alpha_{1}\beta_{1}-\alpha_{2}\beta_{3}\\\
\alpha_{3}\beta_{3}-\alpha_{1}\beta_{2}\\\
\alpha_{2}\beta_{2}-\alpha_{3}\beta_{1}\end{matrix}\right)\end{cases}\,,$

where $i,j,l=0,1$ and we define the notation $<i>\equiv
i~{}(\mathrm{mod}~{}2)$.

##  Appendix B Modular multiplets of weight $4,6,8$ at level $N=4$

From the tensor products of lower weight modular forms with the help of the
Clebsch Gordan coefficients of $S_{4}$ in Appendix A, we can get the higher
weight modular forms. In the following we construct the weight 4, weight 6 and
weight 8 modular multiplets. By expressing the modular forms $Y_{1,2,3,4,5}$
in terms of $\vartheta_{1}$ and $\vartheta_{2}$ as shown in Eq. (2), we can
easily identify the linearly independent modular multiplets of higher weights
without examining the cumbersome non-linear constraints which relate redundant
multiplets coming from tensor products. The linear space of modular forms of
level 4 and weight $k$ has dimension $2k+1$. In the following, we shall
present the explicit form of the modular forms of weight $4,6,8$.

The weight 4 modular forms can be generated from the tensor products of
$Y^{(2)}_{\mathbf{2,3}}$. Using the Kronecker products
$\mathbf{2}\otimes\mathbf{2}=\mathbf{1}\oplus\mathbf{1^{\prime}}\oplus\mathbf{2}$,
$\mathbf{2}\otimes\mathbf{3}=\mathbf{3}\oplus\mathbf{3^{\prime}}$ and
$\mathbf{3}\otimes\mathbf{3}=\mathbf{1}\oplus\mathbf{2}\oplus\mathbf{3}\oplus\mathbf{3^{\prime}}$,
the weight 4 modular forms can be arranged into different $S_{4}$ irreducible
representations $\mathbf{1}$, $\mathbf{2}$, $\mathbf{3}$ and
$\mathbf{3^{\prime}}$. To be more specific, we have

$\displaystyle
Y^{(4)}_{\mathbf{1}}=\left(Y^{(2)}_{\mathbf{2}}Y^{(2)}_{\mathbf{2}}\right)_{\mathbf{1}}=2Y_{1}Y_{2}\,,$
$\displaystyle
Y^{(4)}_{\mathbf{2}}=\left(Y^{(2)}_{\mathbf{2}}Y^{(2)}_{\mathbf{2}}\right)_{\mathbf{2}}=\Big{(}Y^{2}_{2},~{}Y^{2}_{1}\Big{)}^{T}\,,$
$\displaystyle
Y^{(4)}_{\mathbf{3}}=\left(Y^{(2)}_{\mathbf{2}}Y^{(2)}_{\mathbf{3}}\right)_{\mathbf{3}}=\Big{(}Y_{1}Y_{4}+Y_{2}Y_{5},~{}Y_{1}Y_{5}+Y_{2}Y_{3},~{}Y_{1}Y_{3}+Y_{2}Y_{4}\Big{)}^{T},$
$\displaystyle
Y^{(4)}_{\mathbf{3^{\prime}}}=\left(Y^{(2)}_{\mathbf{2}}Y^{(2)}_{\mathbf{3}}\right)_{\mathbf{3^{\prime}}}=\Big{(}Y_{1}Y_{4}-Y_{2}Y_{5},~{}Y_{1}Y_{5}-Y_{2}Y_{3},~{}Y_{1}Y_{3}-Y_{2}Y_{4}\Big{)}^{T}\,.$
(B.5)

Then we proceed to construct weight 6 modular multiplets from the tensor
products of $Y^{(2)}_{\mathbf{2,3}}$ and
$Y^{(4)}_{\mathbf{1,2,3,3^{\prime}}}$. We find thirteen linearly independent
weight 6 modular forms which can be decomposed into two singlets $\mathbf{1}$,
$\mathbf{1}^{\prime}$, a doublet $\mathbf{2}$ and three triplets $\mathbf{3}$
and $\mathbf{3}^{\prime}$ of $S_{4}$ as follows

$\displaystyle
Y^{(6)}_{\mathbf{1}}=\left(Y^{(2)}_{\mathbf{2}}Y^{(4)}_{\mathbf{2}}\right)_{\mathbf{1}}=Y^{3}_{1}+Y^{3}_{2}\,,$
$\displaystyle
Y^{(6)}_{\mathbf{1^{\prime}}}=\left(Y^{(2)}_{\mathbf{2}}Y^{(4)}_{\mathbf{2}}\right)_{\mathbf{1^{\prime}}}=Y^{3}_{1}-Y^{3}_{2}\,,$
$\displaystyle
Y^{(6)}_{\mathbf{2}}=\left(Y^{(2)}_{\mathbf{2}}Y^{(4)}_{\mathbf{1}}\right)_{\mathbf{2}}=2Y_{1}Y_{2}\Big{(}Y_{1},~{}Y_{2}\Big{)}^{T}\,,$
$\displaystyle
Y^{(6)}_{\mathbf{3}I}=\left(Y^{(2)}_{\mathbf{3}}Y^{(4)}_{\mathbf{1}}\right)_{\mathbf{3}}=2Y_{1}Y_{2}\Big{(}Y_{3},~{}Y_{4},~{}Y_{5}\Big{)}^{T}\,,$
$\displaystyle
Y^{(6)}_{\mathbf{3}II}=\left(Y^{(2)}_{\mathbf{3}}Y^{(4)}_{\mathbf{2}}\right)_{\mathbf{3}}=\Big{(}Y_{2}^{2}Y_{4}+Y_{1}^{2}Y_{5},~{}Y_{2}^{2}Y_{5}+Y_{1}^{2}Y_{3},~{}Y_{2}^{2}Y_{3}+Y_{1}^{2}Y_{4}\Big{)}^{T}\,,$
$\displaystyle