11institutetext:  LIX, Ecole Polytechnique, Institut Polytechnique de Paris
Palaiseau, France
22institutetext: DEIF Wind Power Technology Austria GmbH
Klagenfurt, Austria

# An Improved Yaw Control Algorithm for Wind Turbines via Reinforcement
Learning

Alban Puech 1122 Jesse Read 11

###### Abstract

Yaw misalignment, measured as the difference between the wind direction and
the nacelle position of a wind turbine, has consequences on the power output,
the safety and the lifetime of the turbine and its wind park as a whole. We
use reinforcement learning to develop a yaw control agent to minimise yaw
misalignment and optimally reallocate yaw resources, prioritising high-speed
segments, while keeping yaw usage low. To achieve this, we carefully crafted
and tested the reward metric to trade-off yaw usage versus yaw alignment (as
proportional to power production), and created a novel simulator (environment)
based on real-world wind logs obtained from a REpower MM82 2MW turbine. The
resulting algorithm decreased the yaw misalignment by 5.5% and 11.2% on two
simulations of 2.7 hours each, compared to the conventional active yaw control
algorithm. The average net energy gain obtained was 0.31% and 0.33%
respectively, compared to the traditional yaw control algorithm. On a single
2MW turbine, this amounts to a 1.5k-2.5k euros annual gain, which sums up to
very significant profits over an entire wind park.

###### Keywords:

Wind turbine control Multi-objective reinforcement learning yaw control

##  1 Introduction

As the world tries to move away from fossil fuels, wind energy appears to be
one of the most promising renewable energy sources. Wind turbines convert
kinetic energy from the wind into electricity. Energy output mainly depends on
the wind speed, the turbine blade diameter and the generator size, but it is
now known that the layout, as well as the operation of turbines, can have a
significant impact on the power performance of wind parks [1, 24]. In
particular, the orientation of the wind turbine rotor against the wind,
achieved by the yaw mechanism, is of key importance to ensure that the maximal
amount of wind energy is extracted. To do so, one needs to minimise the yaw
misalignment $\gamma$, or yaw error, measured as the difference between the
wind direction $\phi$ and the nacelle position $\theta$:

$\gamma=\phi-\theta$ (1)

Based on a collection of 300 turbines, [17] showed that more than 50 % of wind
turbines were misaligned to the wind direction, while in [5], the loss yield
by a misalignment of only 7 degrees on a 2MW wind turbine was estimated to be
around 6000 euros/year. Yaw misalignment also has an impact on the wind park
as a whole since it increases the wake effect left by the rotor on the
downstream turbines, as shown by [7]. Finally, yaw misalignment reduces the
lifespan of the turbine, putting them under high constraint loads [22].

The yaw control algorithm that we use in the present paper as a baseline is a
conventional yaw control algorithm (CYCA) deployed in most modern fixed and
variable speed horizontal axis wind turbines. In this paper, the results of
CYCA are obtained in two different ways. When needed, we make the distinction
between the version of CYCA running on the wind turbine and which control
sequence (sequence of control action) is constructed from the turbine logs
(CYCA-L), and its simulated version CYCA-S, which control sequence is obtained
using a yaw control simulator from DEIF. It is important to note that both
CYCA versions correspond to the same algorithm, but their results differ due
to external factors (e.g. maintenance, cable untwisting, emergency stop) that
interfere with the correct execution of CYCA while running on the turbine.
CYCA takes the wind direction measured over a small duration preceding the yaw
control as an estimate for the nacelle position to reach during the next
control cycle. It typically suffers from irreducible operation delays stemming
from lags in communications and its fixed (relatively slow) yaw speed. Beyond
slow reaction time to wind direction change, the yaw resources are
intentionally limited by the controller because a too-frequent yaw usage risks
damage to key components of the mechanism such as the bearings.

To limit unnecessary yaw actuation and load on the mechanism, CYCA actuates
the yaw mechanism only after the cumulative error of yaw misalignment exceeds
some threshold. As a result, the same importance is given to yaw alignment in
low-speed segments, where the potential energy gain by reducing yaw
misalignment is low, and to higher speed segments, where most of the energy
gain can be achieved. We summarise the working principle of CYCA in Fig. 1. An
optimal yaw system control algorithm would allocate yaw resources in order to
minimise the yaw misalignment, with a priority given to the regions where most
of the energy gain can be obtained.

Figure 1: Simplified control loop of CYCA. The cumulative yaw error is updated
with the new yaw error computed from the wind direction measurement. The value
of the threshold defines the sensibility of CYCA to small or short yaw
misalignment. The length of the control cycle varies but is usually around 1
second.

Little research has been done on new yaw control strategies compared to other
turbine operations, such as pitch control. Some yaw control strategies rely on
direct or indirect measurement of the wind direction, while others rely on
maximum power point tracking. Most of them depend on wind measurements that
are sensitive to external factors. Moreover, yaw alignment is often achieved
at the expense of yaw usage. They also do not consider the power output for
the allocation of the yaw resources. To address those issues, we propose in
this paper a new data-driven yaw control strategy that bases the control
decision on prediction and pattern recognition rather than past measurements.
This new strategy, relying on a Reinforcement-Learning (RL)-trained agent,
increases the power output by reducing the overall yaw misalignment and by
intelligently re-allocating yaw resources. The contributions of this paper
are:

  * •

A novel RL-based yaw control algorithm for wind turbines that optimally
allocates yaw resources and better tracks wind direction change, keeping yaw
usage low.

  * •

A new realistic simulation yaw environment that allows training and
benchmarking of our algorithm against several baselines (simulated and real-
world conventional yaw control algorithm) using past wind conditions obtained
from turbine logs.

  * •

A new multi-objective reward function to handle short wind direction
variations without the use of deadband functions and new evaluation metrics
that take into account the yaw drive consumption changes.

We organise the remainder of the paper as follows: We will first introduce the
novel Reinforcement-Learning control algorithm (RLYCA) and discuss its key
components. We will then focus on the RL environment used for both training
and assessment of the novel strategy. Finally, we will discuss the results of
the proposed strategy and the possible improvement opportunities.

##  2 Related Work and Existing yaw control strategies

There are mainly 4 different types of active yaw control strategies. Contrary
to passive yaw control strategies, which are based on wind force to orient the
rotor, active yaw control strategies use a servomechanism to change the rotor
position. Among the active yaw strategies, three of them are deterministic:
the ones with advanced measurements, without wind direction measurement and
with indirect measurements. Recently, RL methods were proposed and output a
control action based on pattern recognition.

Yaw control with advanced measurements. Yaw control with advanced measurements
relies on remote instruments like Lidars and Sodars to get more accurate wind
direction and speed values, as in [5]. As the wind direction measure is made
upstream of the turbine, it is not affected by the wake effect and gives a
more accurate measure of the yaw misalignment, which is then used by an active
yaw control algorithm. Alternatively, a nacelle-mounted Lidar can be used to
learn a correction function that is then applied to the wind direction
measurement obtained from the wind vane [12]. Both methods were shown to help
at correcting the yaw angle, but they require expensive measurement tools.

Yaw control without wind direction measurement. Other methods, without wind
direction measurement, have been proposed recently [15, 11]. They rely on
maximum power point tracking. The difference between the actual power output
and the theoretical power output computed using the power curve gives
information on the power loss coming from the yaw misalignment, and, in turn,
gives an estimate of the yaw misalignment. A benefit of those methods is that
they no longer depend on the wind direction measurements that are often, as
discussed previously, disturbed by the operation of the wind turbine and by
the wake effect. However, they are sensitive to wind speed variation since the
theoretical output power is computed out of it and to any other external
factor such as snow or dirt that could reduce the power output of the turbine.

Yaw control with indirect measurement. The methods that gained the most
attention recently are the ones based on indirect measurements. They rely on
short term prediction of the wind direction to feed a model predictive control
(MPC) algorithm that anticipates wind direction change [10, 8, 9]. Although
the algorithm is no longer directly using the wind direction measure, its
performance is still highly dependent on the accuracy of the wind direction
prediction. In high wind variability conditions, the accuracy tends to be low
(Dzulfikri, Nuryanti and Erdani [10] showed a mean absolute error of 9.9
degrees for a 1-minute ahead prediction of the wind direction on the testing
set), which restricts the use of such methods to a relatively small space of
wind conditions. It is important to note that, to this day, none of the
proposed methods relying on wind direction prediction significantly
outperforms the traditional active yaw control strategy at reducing the yaw
misalignment.

RL-based yaw control strategies. Finally, the RL methods are an alternative to
the other data-driven methods presented previously as they do not use the wind
direction and speed measurements in a predefined, deterministic way. They
rather determine the best control action based on the current wind and
position state, learning from different simulation scenarios or online, while
running on the turbine controller. Those techniques are the most promising as
they always output the optimal control action to the best of the model
knowledge acquired from the training. They work in a black box manner and are
not impacted by the measurement inaccuracy since no control strategy needs to
be designed in advance, and since no decision relies directly on the value of
the wind direction or speed. The performance of the models then depends on the
abstraction and generalisation potential of the algorithm and the quantity and
relevance of the training data for the given task, which can be a limitation
in situations where we do not have enough data. Those strategies make easier
the task of reallocating yaw resources by prioritising segments of high
potential gain since the algorithm can learn to identify those segments.

To this day, most of the research proposing an RL-trained agent focuses on
long term control. Saenz-Aguirre [19] proposed a model that actuates the
nacelle position every 24 hours only. This present paper shows a method that
is suitable for most wind turbines that usually spend 3 to 10 % of the time
yawing.
In addition to that, the wind conditions used in the simulators often make the
control decision easier to determine. The wind direction used in the
simulations ran to test the control algorithms against CYCA often has the same
mean over any segment of the dataset and shows periodic sudden change with
some added noises of small magnitude. The wind speed is often taken to be
almost stable with very small variations around its mean [19, 20]. Instead, we
base our simulations on real-world data obtained from turbine logs. Most of
the papers proposing new yaw control techniques do not show their performances
for different levels of yaw usage and the bench-markings against current yaw
control strategies are often limited. In the present paper, we offer a clear
comparison with both CYCA-S and CYCA-L and we take into account the power
consumption change of the yaw drive in the computation of the energy gains.
The originality of this paper is that it offers a yaw control algorithm that
can directly and realistically replace CYCA, increasing the power extracted
while keeping the yaw usage reasonable.

##  3 Proposed Yaw Control Strategy

###  3.1 Goals and constraints of the strategy

Figure 2: Left: Power curve of 2MW turbine and its different regions. The
turbine operates only after cut-in and before cut-out. The rated power
corresponds to the maximum power the turbine can produce. Right: Power loss
factor as a function of the yaw misalignment with $\alpha=3$.

One of the main motivations for our new yaw control algorithm is to maximise
the wind power extracted by the turbine. Power loss due to yaw misalignment
occurs in region 2 of the power curve, shown in Fig.2. Indeed, in region 1,
that is, before cut-in, the wind is too slow and the turbine is not generating
electricity. In region 4, after cut-out, the turbine is shut down and the
rotor is braked to avoid over-speeding. Moreover, yaw misalignment does not
have a significant impact on the power production in region 3, since the
turbine is already outputting the maximum amount of power that it can output.
However, the turbine will require a larger wind speed to reach rated power in
case of yaw misalignment. In region 2, the power $P_{\gamma}$ generated by a
wind turbine with a yaw misalignment $\gamma$ is expressed in degrees as:

$P_{\gamma}=P_{*}\cos^{\alpha}\gamma$ (2)

where $P_{*}$ is the power extracted assuming perfect yaw alignment. Moreover,
$P_{*}$ is directly related to the wind speed $v$, the circular swept area of
the turbine $A$, the air density $\rho$ and the Betz factor [23] $c$:

$P_{*}=\frac{1}{2}\rho Av^{3}c$ (3)

Many values have been proposed for the cosine exponent $\alpha$. Its value
significantly impacts the calculated loss caused by yaw misalignment and it is
often assumed that the value of $\alpha$ is three, as supported by [3] and by
analysis conducted by [14]. Most recent studies like [22] showed that the wake
effect caused by upstream turbines can have an impact on the value of
$\alpha$. Dahlberg and Montgomerie [4] showed that the value of the cosine
exponent could range from 1.7 to 5.1. Thus, we will stick $\alpha=3$,
corresponding to the theoretical analysis proposed in [14]. The power loss
factor associated with a yaw misalignment in region 2 is shown in Fig. 2. The
other benefit of reducing the yaw misalignment is to prevent the turbine from
undergoing large mechanical and fatigue loads. Damiani [6] showed that
misaligned turbines are more subject to damage in case of a storm but also see
their lifespan reduced. Fleming [13] showed an increase in the blade, yaw
bearings and tower loading in a situation of yaw misalignment.

Yaw misalignment minimisation should not be achieved at the expense of yaw
usage. Indeed, yawing too often and for too long can harm the yaw mechanism.
It is thus crucial to keep the time spent yawing under some limit mainly
defined by the turbine safety guidelines, in our case, less than 10 % of the
time should be spent yawing. Hence, the goal is not to achieve perfect yaw
alignment, but rather to find the optimal trade-off between yaw correction and
yaw actuation. The main constraint is to avoid actuating the yaw systems too
often and for too long. As the loss induced by yaw misalignment is cubic,
reducing the yaw misalignment when it is large already allows us to decrease
most of the loss caused by the misalignment. Moreover, and as discussed
earlier, most of the yaw usage should be done on high-speed segments, where
yaw misalignment has a higher cost on the power output.

###  3.2 A proposed yaw control strategy via Reinforcement Learning

Reinforcement learning (RL) involves teaching an agent progressing in an
environment to take actions that maximise expected gains defined by a reward
function. The agent has limited information concerning its current state and
learns its policy based on experience from repeated episodes.

We define a control cycle as the episode of fixed duration (e.g. 3 seconds)
during which information is received by the agent which then outputs a control
action. In the case of RLYCA-A, during a control cycle, the agent takes an
action $A_{t}$ by applying its internal policy, upon observing its current
state $S_{t}$. After the action is taken, a reward $R_{t+1}$ and a new state
$S_{t+1}$ is computed. The corresponding RL task can be described by Fig. 3.

Figure 3: Simplified block Diagram of the RL task. The new state $S_{t+1}$ and
the reward $R_{t+1}$ is computed only after the agent took its action $A_{t}$.

We train the agent to learn the optimal control policy for the yaw system.
Optimal control theory and reinforcement learning are closely related as they
both refer to techniques used to find a control for a system that evolves in
time and interacts with its environment, such that a certain objective
function is maximised. However, an optimal control algorithm takes an action
that maximises the predefined objective function in a deterministic way. In
reinforcement learning, the agent does not know the reward function in advance
and has to learn it through experiences to define its policy.

In this paper, the proposed yaw control strategy is based on a Proximal Policy
Optimisation Algorithm (PPO) [21]. Unlike Deep Q-learning, used in [19, 20],
PPO can not learn from offline stored data and requires online training.

Designing the action space. At each control cycle, the agent outputs one of
the following actions:

  * •

0: Clockwise rotation for the full duration of the control cycle

  * •

1: No rotation

  * •

2: Counterclockwise rotation for the full duration of the control cycle

Defining the reward function. The goal of the agent is to learn the optimal
policy, which is the one that maximises the expected return, computed as the
sum of rewards over an episode. We define a reward function in two parts, each
accounting for one of our two goals: maximise power production, and limit the
yaw mechanism usage.
The first term $R^{(1)}$ of the reward function accounts for the power
production maximisation aspect. To maximise the power extracted during the
control cycle $t$, we need to minimise the yaw misalignment $\gamma_{t}$ as
expressed in Eq. (1). However, more importance should be given to the control
cycles where the wind speed $v_{t}$ is high (but still lower than the rated
speed), as Eq. (3) shows that the power extracted assuming perfect alignment
$P_{0}$ increases with the cube of the wind speed $v$. Hence, we define
$R^{(1)}_{t+1}$, the first negative component of the reward for the action
taken during the control cycle $t$:

$R^{(1)}_{t+1}=-\gamma_{t+1}^{2}\tilde{v}_{t+1}^{3}$ (4)

Weighting the yaw misalignment $\gamma_{t+1}$ squared by the standardised wind
speed $\tilde{v}_{t+1}$ cubed allows us to penalise more yaw misalignment when
the power output is high, and thus encourages a better yaw resource allocation
The second term of the reward function, denoted by $R^{(2)}$, accounts for the
second objective of RLYCA-A: yaw usage limitation. To keep the yaw count low
and reduce the time spent yawing, we should prevent the agent from reacting to
short wind direction variations (in time) or to small yaw misalignments (in
magnitude) that are not worth being corrected compared to larger yaw
misalignment. Indeed, since the loss increases with the cube of the yaw
misalignment, a correction of one degree will yield more gain when the yaw
misalignment is large than when it is small. We thus want our agent to react
in priority to large yaw misalignment, and to spare the yaw mechanism the rest
of the time. Hence, we reward the agent for observing a period of at least $k$
control cycles without taking any moving actions. We define $R^{(2)}_{t+1}$,
the second component of the reward for the action $a_{t}$ taken during control
cycle $t$ as

$R^{(2)}_{t+1}=w\prod_{i=t-k+1}^{t}\mathbf{1}\\{a_{i}=1\\}$ (5)

where the weight $w\in\mathbb{R^{+}}$ encodes the importance given to the
preservation of the yaw mechanism.To illustrate the principle of this loss
function, let us assume that $w=40$, that the controller has not taken a
moving action in the last $k-1$ control cycles, and that the wind conditions
remain similar during the control cycle. Assume further that the maximum angle
covered by the yaw mechanism in one control cycle is 3 degrees. If the agent
takes the correct moving action (yaws in the correct direction), the yaw error
can be corrected by a maximum of 3 degrees. The second component of the reward
function will be 0, and the reward will be given by
$R_{t+1}(moving)=-(\gamma_{t}-3)^{2}\tilde{v}_{t}^{3}$.
If the agent chooses the stationary action, the reward will be given by the
sum of the two reward components defined in Eqs. (4) and (5), where the yaw
misalignment remains the same since no correction was made and since we
assumed stable wind conditions. The reward will thus be:
$R_{t+1}(stationary)=-(\gamma_{t})^{2}\tilde{v}_{t}^{3}+40$.
Fig. 4 (left) shows the reward if a stationary action is taken (blue), and if
a moving action is taken (orange). For the moving action to give a larger
reward than the stationary action, the yaw error needs to be larger than 8.1
degrees, where the two curves intersect. This can be interpreted as the fact
that the agent will decide to move the yaw mechanism only if the yaw
misalignment is larger than this value, and will not actuate the yaw mechanism
otherwise. Fig. 4 (right) shows, for different standardised wind speed values,
the minimum yaw misalignment for the moving action to give a larger reward
than the stationary action. We designed the reward function such that, when
the wind speed is high, the agent is more willing to correct small yaw
misalignment. Note that in region 3 of the power curve, the reward function
penalises yaw misalignment in the same manner as described before, but not for
the same reason. Indeed, although there is no power loss due to yaw
misalignment in this region, the constraint loads are higher, which justifies
yaw misalignment to be equally penalised.

Figure 4: Left: Reward for the stationary action and the moving action as a
function of the yaw misalignment under fixed $\tilde{v}=1$ and $w=40$. Right:
Minimum yaw misalignment (vertical axis) starting from which the moving action
has a higher reward than the stationary action (for $w=40$); as the wind speed
increases, we encourage more and more yaw correction.

Choosing the state space. The reward being dependent on the yaw misalignment,
the wind speed and the past actions, we aim at including that information into
the state $S_{t}$ given to the agent to take its action $A_{t}$ during the
control cycle $t$. We define $j$ as the number of lagged values or ancestor
values given for each feature of the state. $S_{t}$ is thus expressed as in
Eq. (6), where $\gamma_{t}$, $\phi_{t}$ and $\tilde{v}_{t}$ denote the yaw
misalignment, the wind direction and the standardised wind speed during the
control cycle $t$, respectively.

$S_{t}=\begin{bmatrix}a_{t-1}&\gamma_{t}&\phi_{t}&\tilde{v}_{t}\\\
a_{t-2}&\gamma_{t-1}&\phi_{t-1}&\tilde{v}_{t-1}\\\
a_{t-3}&\gamma_{t-2}&\phi_{t-2}&\tilde{v}_{t-2}\\\
\vdots&\vdots&\vdots&\vdots\\\
a_{t-j}&\gamma_{t-j+1}&\phi_{t-j+1}&\tilde{v}_{t-j+1}\end{bmatrix}$ (6)

It is important to note that, in the control cycle $t$, state $S_{t}$ contains
the wind direction $\phi_{t}$ and wind speed $v_{t}$, but not $v_{t+1}$ nor
$\phi_{t+1}$ (the agent cannot know in advance the future wind direction and
speed).

##  4 Yaw system simulation environment

We develop a yaw system environment that we use with both RLYCA-A and CYCA-S.
In the case of RLYCA-A, we use the environment for both training and testing.
In the case of CYCA-S, we use the yaw system environment for testing only as
the latter algorithm does not require to be trained. Finally, CYCA-L testing
does not require any simulator since the control sequence, as well as all the
measurements needed for the bench-marking, are directly extracted from the
turbine logs.

###  4.1 Datasets

The yaw system simulation environment simulates the wind conditions recorded
at a real-life wind park. Namely, we collected two time-series data sets each
consisting of 21000 turbine log points, recorded in September and February of
2021. Each $t$-th data point ($t=1,\ldots,21000$) contains the mean wind
direction and the mean wind speed over 1 second, stacked into state $S_{t}$ as
illustrated in 6. Dataset 1 exhibits a relatively steady wind direction; of
$34.1^{\circ}\pm 9.7$ with an interquartile range of $13.7^{\circ}$. Dataset 2
has similar statistical features (wind direction of $41.4^{\circ}\pm 11.6$
with an interquartile range of $14.2^{\circ}$) but shows fast wind direction
changes of large magnitude between 10000s and 12500s and between 15000s and
20000s.

Figure 5: The first dataset (left) shows a steady wind direction while the
second one (right) shows 2 large wind direction changes. Both are used for our
experiments. For each dataset, we train RLYCA-A on its first half of a dataset
and test it on the other.

###  4.2 Working principle of the yaw system simulation environment

The working principle of the yaw system environment used for the RL agent is
described in Fig. 6. We update the new nacelle position $\theta_{t+1}$ as in
7, where $r$ denotes the yaw rate and $p$ the duration of the control cycle.

$\theta_{t+1}=\theta_{t}+p(action-1)r$ (7)

The new yaw misalignment is then computed using the new nacelle position and
wind direction. We compute the power assuming perfect yaw alignment using the
power curve shown in Fig. 2 and the power under the effect of yaw misalignment
is obtained from Eq. (2), provided that $P_{*}$ belongs to region 2 of the
power curve. Otherwise, we assume no power loss ($P_{\gamma}=P_{*}$). We
simulate a communication delay of a duration equal to the period of the
control cycle. We can assume that the yaw mechanism receives the control
signal $a_{t}$ when the next control cycle ($t+1$) is executed. In practice,
this means that we iterate in the wind time series before we compute the new
yaw misalignment.This greatly contributes to the complexity of the task: the
wind conditions can change drastically between the real-time where a given
action is decided and its application.

Figure 6: Working principle of the yaw system environment

##  5 Experiments and Discussions

In this section, we compare our novel algorithm RLYCA-A with benchmarks CYCA-L
and CYCA-S on the two datasets which we develop into a simulated environment
using open AI gym [2].

###  5.1 Implementation and parameters

RLYCA-A uses a PPO agent, implemented in PyTorch [16] and stablebaselines3
[18]. We use an Actor-Critic policy with a 2 dense layers of 64 neurons for
the value network and 2 layers of 64 neurons for the policy network; trained
over 200,000 steps. For CYCA-S we obtain the control sequence from a control
script developed at DEIF. For CYCA-L we reconstruct its control sequence from
the real-world turbine logs that saved the position of the nacelle at the same
time as the wind condition used for the experiment was saved.

The power curves used for the yaw system simulated environment is shown in
Fig. 2. Table 1 gives the parameters used for the yaw system simulated
environment.

Table 1: Parameters used for the experiment on both datasets. They are
identical to those of the turbine from which the controller logs were obtained
(REpower MM82).

parameter name | value
---|---
environment parameters
$k$ number of consecutive stationary action to get reward | 2
$j$ number of lagged values in state | 12
$w$ weight of $R^{(2)}$ | 40
$v_{r}$ rated speed | 14m/s
$p$ cycle period | 10s
$comm$ communication delay | 10s
$P_{yawdrive}$ | 18 kWh
RL agent parameters
learning rate | 0.003
number of steps before update | 2048
batch size | 64
number of epochs | 10
discount factor | 0.99
generalised Advantage Estimator factor | 0.95

The yaw alignment performances of the algorithms are compared using the
average yaw misalignment computed as stated in Eq. (1). We assess the power
extraction performance using the estimate of the power produced given by the
yaw system environment and stated in the previous subsection. Since the three
experiments are conducted in the same wind conditions given by the testing
dataset, the only variable that we use in the computation of the power output
and that differs between the three approaches is the nacelle position, which
is the direct result of the applications of the three algorithms. We compare
the yaw usage of all methods using three metrics: The proportion of the time
spent yawing over the simulated duration of the experiment, the angle covered
by the yaw mechanism, and the number of yaw actuations. Additionally, we
compute the difference in yaw mechanism power consumption $\Delta(t)$ as shown
in Eq. (8), where $\delta_{\phi_{RL}}(t)$ and $\delta_{\phi_{simu}}(t)$ denote
the change in nacelle position angle during the control cycle $t$ obtained
with RLYCA-A and with CYCA-S, respectively. Moreover, $r$ denotes the yaw rate
and $P_{yawdrive}$ denotes the power consumption of the yaw mechanism.

$\Delta(t)=\frac{\delta_{\phi_{RL}}(t)-\delta_{\phi_{simu}}(t)}{r}P_{yawdrive}$
(8)

###  5.2 Experiment in steady wind direction conditions (dataset/simulation 1)

We first compare the algorithms in steady wind direction conditions. Fig. 7
shows the nacelle position obtained for the three algorithms. Most of the yaw
misalignment decrease is achieved between 5000 and 7000 seconds (first bar
plot from top). This results in a large power output increase in the same time
interval with a cumulative gain of 1.5kw between 6000 and 6500 seconds. This
power production improvement requires more yaw actuations, which translates
into a yaw consumption increase that is mainly concentrated in the segment
where the wind direction is decreasing. The yaw usage is also increased in the
first half of the dataset but it does not improve the power output. This
phenomenon is something that should be corrected in a future version of
RLYCA-A. The power output change, net of any increase in the yaw consumption,
is negative (power loss) in the first half of the dataset because the increase
in yaw usage does not lead to a power production increase. Indeed, the power
output is increased, but the yaw consumption change is larger than this gain,
which results in a negative net power change. The region where the net power
output is positive corresponds to a high speed (and thus high power) region so
that the relative energy gain is translated into large energy gains. On the
opposite, the region where the net power change is negative is a low-speed
region, so that the absolute power loss does not affect too much the overall
performance of the algorithm. Table 2 shows the results obtained with the
three algorithms, while Table 3 summarises the results achieved by RLYCA-A.
All in all, the algorithm succeeded in allocating more resources to the high-
speed segment, resulting in a large power output increase. The relative energy
gain achieved by RLYCA-A compared to CYCA-S is 0.4 %, and the net energy gain
amounts to 0.31 %. This gain is due to a better allocation of the yaw
resources.

###  5.3 Experiment in variable wind direction conditions (dataset/simulation
2)

We now compare the algorithms in the simulation derived on the variable wind
direction dataset. Fig. 8 shows the nacelle position obtained for the three
algorithms. Most of the yaw misalignment decrease is achieved in the second
half of the dataset. RLYCA-A reacts quicker to the change in wind direction
than CYCA-S.

Figure 7: Comparison of control signals of the three algorithms on the steady
wind direction dataset. The power and yaw error comparisons shown on the bar
plots are made against CYCA-S.

The power output is increased in almost every region of the dataset. The yaw
consumption is also more important but the increase in power output
compensates for such additional power consumption everywhere except after 9000
seconds. Indeed, the additional yaw actuations after 9000 seconds do not
translate into a better yaw alignment. All in all, the relative energy gain
achieved by RLYCA-A is 0.75 %, but the net energy gain is 0.33, as the yaw
usage is doubled. The wind speed is stable on the testing set, as shown in
Fig. 5, so that, this time, the energy gain can not be explained by a better
allocation of the yaw resources. In this example, it is the overall yaw
misalignment decrease achieved by RLYCA-A (18.5 % compared to CYCA-S and 24.4
% compared to CYCA-L) that explains such a high energy gain.

Figure 8: Comparison of control signals of the three algorithms on the
variable wind direction dataset. The power and yaw error comparisons are made
against CYCA-S. Table 2: Result summary of the three algorithms. Only the
power outputs for CYCA-S and RLYCA-A are shown as the power output for CYCA-L
is impacted by external factors that are not considered in the simulation
(e.g. dirt or snow).

| Steady wind | variable wind
---|---|---
| RLYCA-A | CYCA-S | CYCA-L | RLYCA-A | CYCA-S | CYCA-L
average yaw error (deg) | 6.18 | 6.52 | 6.91 | 5.71 | 7.01 | 7.56
power output (kWh) | 1173.1 | 1168.5 |  | 741.5 | 736.0 |
angle covered (deg) | 111.0 | 53.5 | 27.3 | 261 | 112.1 | 91.8
yaw count | 34 | 15 | 5 | 73 | 27 | 14
time spent yawing (%) | 3.7 | 2.0 | 0.4 | 8.7 | 4.2 | 1.6

Table 3: Performance of RLYCA-A compared to CYCA-S. The values indicated
between parenthesis indicate the comparison with CYCA-L.

| steady wind | variable wind
---|---|---
average yaw error decrease (%) | 5.5 (11.2) | 18.5 (24.4)
energy gain (%) | 0.40 | 0.75
net energy gain (%) | 0.31 | 0.33

##  6 Conclusion and Future Work

In this paper, we proposed a reinforcement learning yaw control strategy to
optimally allocate yaw resources and track wind direction change, with the
goal to maximise the power extracted, while keeping yaw usage under a 10 %
bound. We compared our algorithm to a simulation of the current algorithm
running on a REpower MM82 and to the control signals gathered from the turbine
logs. The comparison was done on two 10000 seconds-long test datasets obtained
from past wind measurements logged by the turbine controller. The RL algorithm
achieved better wind direction change tracking, which led to a decrease in the
yaw misalignment. It also better allocates yaw resources, prioritising high
power (and thus high potential gain) segments. Compared to the simulated
baseline yaw control algorithm, our new algorithm improves by 5.5 and 18.5 %
the yaw alignment in the steady wind and variable wind conditions,
respectively. When comparing its performance to the control sequence obtained
from the turbine logs, the last metrics become 11.2 and 24.4 %. Taking into
account the increase in yaw drive consumption, the net energy gain observed on
the two datasets was 0.31 and 0.33 %. For a single 2MW wind class 2 turbine,
this amounts to an annual 1500-2500 euros gain11182 meters rotor diameter,
8-11 GWh annual power production, 80 % of the time spent in region 2, cosine-
cube power loss law, 0.09 euros per kWh., which, at the scale of a wind park,
represents substantial gains. Throughout our work, we not only showed the
relevance of data-driven methods for wind turbines control, but we more
generally demonstrated the potential of reinforcement learning to deal with
the uncertain nature of renewable energy resources and renewable energy
structure behaviour. In the future, we plan to use wind turbine digital twins
to have an even more accurate simulation of the yaw mechanism, a model of the
constraint loads, as well as a more accurate estimation of the power output.
The power output could then be directly used in the reward function. We also
plan to include data from upstream turbines to the state space, which could
help in anticipating the wind changes. Finally, the next step would be to move
to the wind farm level, as it has been shown that local yaw alignment can be
detrimental at the farm level and because new intended yaw misalignment
strategies are being discussed.

Reproducible Paper. The code and the data sets used for this paper are
available at https://github.com/albanpuech/RL-yaw-control-algorithm-for-wind-
turbines

Acknowledgements.  We would like to thank Mohamed Alami Chehbourne from LIX as
well as Walter Telsnig, Anatoliy Zabrovskiy and Martin Göldner from DEIF for
the helpful discussions and insights on the topic during the preparation of
this paper.

## References

  * [1] Al-Addous, M., Jaradat, M., Albatayneh, A., Wellmann, J., Al Hmidan, S.: The significance of wind turbines layout optimization on the predicted farm energy yield. Atmosphere 11(1) (2020)
  * [2] Brockman, G., Cheung, V., Pettersson, L., Schneider, J., Schulman, J., Tang, J., Zaremba, W.: OpenAI gym (2016)
  * [3] Burton, T., Sharpe, D., Jenkins, N., Bossanyi, E.: Wind Energy Handbook. John Wiley & Sons (2001), https://books.google.at/books?id=4UYm893y-34C
  * [4] Dahlberg, J., Montgomerie, B.: Technical report, final report part 2, wake effects and other loads. Research program of the utgrunden demonstration offshore wind farm (02 2005)
  * [5] Dalmas J., P.A.: Étude comparative de trois appareils pour la mesure d’alignement nacelle. Conference: 5th International Renewable Energy Congress (IREC) (2014)
  * [6] Damiani, et al.: Assessment of wind turbine component loads under yaw-offset conditions. Wind Energy Science 3, 173–189 (04 2018)
  * [7] Dijk, M., Wingerden, J.W., Ashuri, T., Li, Y., Rotea, M.: Yaw-misalignment and its impact on wind turbine loads and wind farm power output. Journal of Physics: Conference Series 753, 062013 (09 2016)
  * [8] Dongran, et al.: Wind direction prediction for yaw control of wind turbines. International Journal of Control, Automation and Systems 15 (08 2017)
  * [9] Dongran, S., Li, L., Yang, J., Joo, Y.H.: A model predictive control for the yaw control system of horizontal-axis wind turbines. Energy Procedia 158, 237–242 (02 2019)
  * [10] Dzulfikri, Z., Nuryanti, N., Erdani, Y.: Design and implementation of artificial neural networks to predict wind directions on controlling yaw of wind turbine prototype. Journal of Robotics and Control (JRC) 1(1), 20–26 (2019)
  * [11] Farret, F., Pfitscher, L., Bernardon, D.: Sensorless active yaw control for wind turbines. Journal of Science and Medicine in Sport 2, 1370 – 1375 vol.2 (02 2001)
  * [12] Fleming, et al.: Field-test results using a nacelle-mounted lidar for improving wind turbine power capture by reducing yaw misalignment. Journal of Physics: Conf. Ser. 524 (06 2014)
  * [13] Fleming, P., et al.: Simulation comparison of wake mitigation control strategies for a two-turbine case. Wind Energy 18(12), 2135–2143 (2015)
  * [14] Kragh, K., Hansen, M.: Potential of power gain with improved yaw alignment. Wind Energy 18(6), 979–989 (2015)
  * [15] Mademlis, C., Mesemanolis, A., Karakasis, N., Nalmpantis, T.: Active yaw control in a horizontal axis wind system without requiring wind direction measurement. IET Renewable Power Generation 10 (06 2016)
  * [16] Paszke, et al.: Pytorch: An imperative style, high-performance deep learning library. In: Advances in Neural Information Processing Systems 32, pp. 8024–8035. archive (2019)
  * [17] Pedersen, M.: Yaw misalignment and power curve analysis. EWEA Analysis of operating wind farms (2016)
  * [18] Raffin, A., Hill, A., Gleave, A., Kanervisto, A., Ernestus, M., Dormann, N.: Stable-baselines3: Reliable reinforcement learning implementations. Journal of Machine Learning Research 22(268), 1–8 (2021), http://jmlr.org/papers/v22/20-1364.html
  * [19] Saenz-Aguirre, et al.: Artificial neural network based reinforcement learning for wind turbine yaw control. Energies 12, 436 (01 2019)
  * [20] Saenz-Aguirre, A., Zulueta, E., Fernandez-Gamiz, U., Ulazia, A., Teso-Fz-Betono, D.: Performance enhancement of the artificial neural network–based reinforcement learning for wind turbine yaw control. Wind Energy 23(3), 676–690 (2020)
  * [21] Schulman, J., Wolski, F., Dhariwal, P., Radford, A., Klimov, O.: Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347 (2017)
  * [22] Urbán, A.M., Liew, J., Dellwik, E., Larsen, G.C.: The effect of wake position and yaw misalignment on power loss in wind turbines. Journal of Physics: Conf. Ser. 1222(1) (2019)
  * [23] Venkata, S., Krishnamurthy, S.: Wind energy explained: Theory, design, and application. Power and Energy Magazine, IEEE 1, 50– 51 (12 2003)
  * [24] Wan, S., Cheng, L., Sheng, X.: Effects of yaw error on wind turbine running characteristics based on the equivalent wind speed model. Energies 8(7), 6286–6301 (2015)


CTPU-PTC-23-04
EPHOU-23-006

Quark and lepton hierarchies from $S_{4}^{\prime}$ modular flavor symmetry
Yoshihiko Abe${}^{\ a}$<EMAIL_ADDRESS>Tetsutaro Higaki${}^{\ b}$
<EMAIL_ADDRESS>Junichiro Kawamura${}^{\ b,c}$
<EMAIL_ADDRESS>and Tatsuo Kobayashi${}^{\ d}$
<EMAIL_ADDRESS>

a Department of Physics, University of Wisconsin-Madison, Madison, WI 53706,
USA

b Department of Physics, Keio University, Yokohama, 223-8522, Japan

c Center for Theoretical Physics of the Universe, Institute for Basic Science
(IBS), Daejeon 34051, Korea

d Department of Physics, Hokkaido University, Sapporo 060-0810, Japan

We propose models in which the hierarchical structures of the masses and
mixing in both quark and lepton sectors are explained by the $S_{4}^{\prime}$
modular flavor symmetry near the fixed point $\tau\sim i\infty$. The model
provides the first explicit example which explains hierarchies of both quarks
and leptons. The hierarchies are realized by powers of $\epsilon=e^{2\pi
i\tau/4}=\mathcal{O}(0.01)$ and $2\,\mathrm{Im}\,\tau\sim 5$, where $\tau$
being the modulus. The small parameter $\epsilon$ plays a role of flavon in
the Froggatt-Nielson mechanism under the residual $Z_{4}^{T}$ symmetry, and
powers of $2\,\mathrm{Im}\,\tau$ in the Yukawa couplings are controlled by
modular weights via the canonical normalization. The doublet quarks are
identified to a $S_{4}^{\prime}$ triplet to explain the hierarchical structure
of the quark mixing angles, while the doublet leptons are composed of three
singlets for the large mixing angles in the lepton sector. We show that the
$S_{4}^{\prime}$ modular symmetry alone can explain the hierarchies in both
quark and lepton sectors by $\mathcal{O}(1)$ coefficients.

##  1 Introduction

Modular flavor symmetry is an intriguing way to understand the origin of the
flavor structure of the quarks and leptons in the Standard Model (SM) [1, 2,
3, 4, 5, 6, 7, 8, 9]. Yukawa coupling constants are transformed as modular
forms under the modular symmetry, and hence these are holomorphic functions of
the modulus $\tau$. The finite modular symmetries $\Gamma_{N}$,
$N\in\mathbb{N}$, are generalizations of the discrete non-Abelian flavor
symmetries [10] which have been extensively studied in the literature [11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21]. For instance, $\Gamma_{4}$ is isomorphic
to $S_{4}$ symmetry. Phenomenology of the modular flavor symmetries have been
studied in Refs. [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,
37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55,
56, 57, 58, 59, 60, 61, 62].

Recently, it has been shown that the modular flavor symmetry can explain the
hierarchical structure of the quark masses and the Cabbibo-Kobayashi-Maskawa
(CKM) mixing matrix, based on modular $A_{4}$ [63], $\Gamma_{6}$ [64],
$S_{4}^{\prime}$ [65] and $A_{4}\times A_{4}\times A_{4}$ [66] symmetries. The
lepton sector is also studied in Refs. [67, 68]. In these models, hierarchical
structures are induced where the modulus $\tau$ has its vacuum expectation
value (VEV) near a fixed point, so that a residual $\mathbbm{Z}_{n}$ symmetry
is approximately unbroken. The residual symmetry realizes the well-known
Froggatt-Nielsen (FN) mechanism [69, 70]. For instance, the $m$-th power of
flavon VEV is replaced by $e^{-2\pi m\mathrm{Im}\,\tau/N}$ when the modulus is
stabilized near $\tau\sim i\infty$ and the symmetry generated by $T$, namely
$\mathbbm{Z}_{N}^{T}$, is approximately unbroken. The residual symmetries are
used for model buildings in Refs. [25, 34, 4, 60, 68].

In this work, we extend the analysis in Ref. [65] utilizing
$\Gamma_{4}^{\prime}\simeq S_{4}^{\prime}$ symmetry to explain the quark
hierarchy at $\tau\sim i\infty$. First of all, the lepton sector is studied
under the $S_{4}^{\prime}$ symmetry, so that the mass hierarchies of the
charged leptons are explained in the same manner as those of the quarks, as
well as the neutrino oscillation data. This model provides the first explicit
example in which the hierarchies in both quark and lepton sectors are
explained by a common flavor modular symmetry. Secondly, we consider the co-
existence of non-hatted and hatted representations in a same type of fermions,
whereas we assumed that there is no such co-existence in Ref. [65]. This
generalization allows us to explain the fermion hierarchies by completely
$\mathcal{O}\left({1}\right)$ coefficients unlike the models proposed in our
previous study [65], where there are $\mathcal{O}\left({0.1}\right)$
hierarchies among the coefficients 111 We pointed out that this small
hierarchy can be interpreted by another modular $S_{3}$ flavor symmetry. .

The rest of this paper is organized as follows. The modular flavor symmetry at
$N=4$ is briefly reviewed in Sec. 2, and then the models are constructed in
Sec. 3. Finally, we conclude in Sec. 4. The modular forms are shown in
Appendix A.

##  2 Modular symmetry at $N=4$

We review the modular symmetry at level $N=4$. More detailed discussions are
found in Ref. [7]. We consider the series of groups $\Gamma(N)$,
$N\in\mathbb{N}$, called principal congruence subgroups, defined as

$\displaystyle\Gamma(N):=\left\\{\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in\mathrm{SL}({2},{\mathbb{Z}}),\quad\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\equiv\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}\mod N\right\\},$
(1)

where $\Gamma:=SL(2,\mathbb{Z})=\Gamma(1)$ is the special linear group of
$2\times 2$ matrices of integers with determinant equal to one, i.e. $ad-
bc=1$. The group $\Gamma$ acts on the complex variable $\tau$
($\mathrm{Im}\,\tau>0$) as

$\displaystyle\tau\to\frac{a\tau+b}{c\tau+d}.$ (2)

Under the actions of these generators, $\tau$ is transformed as

$\displaystyle\tau\xrightarrow{S}-\frac{1}{\tau},\quad\tau\xrightarrow{T}\tau+1,\quad\tau\xrightarrow{R}\tau.$
(3)

Since $R$ does not change $\tau$, one can consider the quotient group
$\overline{\Gamma}:=\Gamma/\mathbb{Z}_{2}^{R}$, where $\mathbb{Z}_{2}^{R}$
being $\mathbb{Z}_{2}$ symmetry generated by $R$, which has one-to-one
correspondence with the action to $\tau$. The finite modular group
$\Gamma_{N}^{\prime}$ ($\Gamma_{N}$) is defined as a quotient group
$\Gamma_{N}^{\prime}:=\Gamma/\Gamma(N)$
($\Gamma_{N}:=\overline{\Gamma}/\Gamma(N)$). The group $\Gamma_{N}^{\prime}$
is generated by the generators,

$\displaystyle S=\begin{pmatrix}0&1\\\ -1&0\end{pmatrix},\quad
T=\begin{pmatrix}1&1\\\ 0&1\end{pmatrix},\quad R=\begin{pmatrix}-1&0\\\
0&-1\end{pmatrix},$ (4)

satisfying the following relations,

$\displaystyle S^{2}=R,\quad(ST)^{3}=R^{2}=T^{N}=\mathbf{1},\quad TR=RT.$ (5)

Those for $\Gamma_{N}$ are given by setting $R=\mathbf{1}$. In a case of
$N=4$, $\Gamma_{4}$ is isomorphic to the $S_{4}$ symmetry, and
$\Gamma_{4}^{\prime}$ is isomorphic to the double covering of $S_{4}$, namely
$S_{4}^{\prime}$. Under the $S_{4}^{\prime}$ symmetry, there are 10
irreducible representations,

$\displaystyle
1,1^{\prime},2,3,3^{\prime}\quad\text{and}\quad\hat{1},\hat{1}^{\prime},\hat{2},\hat{3},\hat{3}^{\prime},$
(6)

where the non-hatted and hatted representations transform under $R$ trivially
and non-trivially, respectively, i.e. $Rr=r$ and $R\hat{r}=-\hat{r}$ for a
representation $r$. We choose the basis in which $T$ is diagonal and $S$ is
real. The representation matrices of the doublet $2$ and the triplet $3$ are
respectively given by

$\displaystyle\rho_{S}(2)=\frac{1}{2}\begin{pmatrix}-1&\sqrt{3}\\\
\sqrt{3}&1\end{pmatrix},\quad\rho_{T}(2)=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix},$ (7)

and

$\displaystyle\rho_{S}(3)=-\frac{1}{2}\begin{pmatrix}0&\sqrt{2}&\sqrt{2}\\\
\sqrt{2}&-1&1\\\ \sqrt{2}&1&-1\\\
\end{pmatrix},\quad\rho_{T}(3)=\begin{pmatrix}-1&0&0\\\ 0&-i&0\\\
0&0&i\end{pmatrix}.$ (8)

The primed or/and hatted representations are related as

$\displaystyle\rho_{S}(r)=$ $\displaystyle\
-\rho_{S}(r^{\prime})=-i\rho_{S}(\hat{r})=i\rho_{S}(\hat{r}^{\prime}),$
$\displaystyle\rho_{T}(r)=$ $\displaystyle\
-\rho_{T}(r^{\prime})=i\rho_{T}(\hat{r})=-i\rho_{T}(\hat{r}^{\prime}),$ (9)
$\displaystyle\mathbf{1}=\rho_{R}(r)=$ $\displaystyle\
\rho_{R}(r^{\prime})=-\rho_{R}(\hat{r})=-\rho_{R}(\hat{r}^{\prime}).$

Table 1:  The number of representations of the modular forms at the weight $k\leq 11$ in the $S_{4}^{\prime}$ modular symmetry. The representations for odd weights should be understood as the hatted ones.  weight | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11
---|---|---|---|---|---|---|---|---|---|---|---
$1$ | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 0
$1^{\prime}$ | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 1 | 1 | 1
$2$ | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 2 | 1 | 2 | 2
$3$ | 1 | 0 | 1 | 1 | 2 | 1 | 2 | 2 | 3 | 2 | 3
$3^{\prime}$ | 0 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | 3 | 3

The modular form $Y(\tau)$ with representation $r$ and weight $k$ transforms
under the finite modular symmetry $\Gamma_{4}^{\prime}$ as

$\displaystyle Y^{(k)}_{r}(\tau)\to(c\tau+d)^{k}\rho(r)Y^{(k)}_{r}(\tau),$
(10)

where $\rho(r)$ is a representation matrix. We assume that a matter field $f$,
with representation $r_{f}$ and weight $k_{f}$, transforms in the same manner,

$\displaystyle f\to(c\tau+d)^{k_{f}}\rho(r_{f})f.$ (11)

There are $2k+1$ independent modular forms at a weight $k$. For $k=1$, there
is a $\hat{3}$ representation,

$\displaystyle
Y^{(1)}_{\hat{3}}(\tau)=\begin{pmatrix}\sqrt{2}\epsilon(\tau)\theta(\tau)\\\
\epsilon^{2}(\tau)\\\ -\theta^{2}(\tau)\end{pmatrix}.$ (12)

The functions $\theta$ and $\epsilon$ are defined as [4]

$\displaystyle\theta(\tau):=\frac{\eta(2\tau)^{5}}{\eta(\tau)^{2}\eta(4\tau)^{2}},\quad\epsilon(\tau):=\frac{2\eta(4\tau)^{5}}{\eta(2\tau)},$
(13)

where $\eta(\tau)$ is the Dedekind eta function. Their $q$-expansions are
given by

$\displaystyle\theta(\tau)=1+2\sum_{n=1}^{\infty}q^{n^{2}},\quad\epsilon(\tau)=2q^{1/4}\sum_{n=0}^{\infty}q^{n(n+1)},$
(14)

with $q:=e^{2\pi i\tau}$. These functions are approximately
$\theta(\tau)\simeq 1$ and $\epsilon(\tau)\sim 2q^{1/4}\ll 1$, for
$2\pi\mathrm{Im}\,\tau\gg 1$, where the symmetry generated by $T$, namely
$\mathbb{Z}_{4}^{T}$, is a good symmetry. The modular forms with higher
weights can be constructed from products of $Y^{(1)}_{\hat{3}}$. The number of
representations for the weights $k\leq 11$ are shown in Table 1, and the
explicit forms of the modular forms used in the models are listed in Appendix
A. In this work, we assume that normalizations of the modular forms, which can
not be determined from the modular symmetry, do not induce additional
hierarcihal structures other than those from $\epsilon$.

The residual symmetry $\mathbb{Z}_{4}^{T}$ realizes the FN mechanism. Under
the $\mathbb{Z}_{4}^{T}$ symmetry, $q^{1/4}$ plays the role of flavon whose
charge is unity, because $q^{1/4}\xrightarrow{T}iq^{1/4}$. From Eqs. (7) and
(8), the hierarchical structures of the modular forms of $2$ and $3$ are read
as

$\displaystyle Y_{2}\sim\begin{pmatrix}1\\\ \epsilon^{2}\end{pmatrix},\quad
Y_{3}\sim\begin{pmatrix}\epsilon^{2}\\\ \epsilon^{3}\\\
\epsilon\end{pmatrix}.$ (15)

The hierarchical structures of the other representations can be read in the
same way.

##  3 Models

###  3.1 Hierarchical structures

We aim to explain the hierarchical structure of the quarks and leptons in the
SM by the $S_{4}^{\prime}$ modular flavor symmetry. In terms of
$\epsilon\sim\mathcal{O}\left({0.01}\right)$, the hierarchies of the fermion
masses may be expressed as

$\displaystyle(m_{u},m_{c},m_{t})\sim$ $\displaystyle\
(\epsilon^{3},\epsilon^{2},1),\quad(m_{d},m_{s},m_{b})\sim\epsilon^{p}(\epsilon^{2},\epsilon,1),\quad(m_{e},m_{\mu},m_{\tau})\sim\epsilon^{p}(\epsilon^{2},\epsilon,1),$
(16)

where $p=0,1$, and the CKM and PMNS matrices are given by

$\displaystyle V_{\mathrm{CKM}}\sim$ $\displaystyle\
\begin{pmatrix}1&\epsilon&\epsilon^{2}\\\ \epsilon&1&\epsilon\\\
\epsilon^{2}&\epsilon&1\end{pmatrix},\quad
V_{\mathrm{PMNS}}\sim\begin{pmatrix}1&1&1\\\ 1&1&1\\\ 1&1&1\end{pmatrix}.$
(17)

In addition to these hierarchical structures, there are hierarchical
structures from $2\mathrm{Im}\,\tau=-(4/\pi)\log\epsilon\sim 5$ due to the
canonical normalization of matter kinetic terms depending on the modular
weights as discussed later. The level $N=4$ is the minimal number to realize
these texture with up to $\epsilon^{3}$ which may be necessary to explain the
observed values with $\mathcal{O}\left({1}\right)$ coefficients.

The texture in Eq. (16) is realized if the representations of the quarks and
leptons are

$\displaystyle u^{c}=$ $\displaystyle\ 1\oplus 1\oplus\hat{1}^{\prime},\quad
d^{c}=\begin{cases}1\oplus 1\oplus 1\\\
\hat{1}^{\prime}\oplus\hat{1}^{\prime}\oplus\hat{1}^{\prime}\end{cases},\quad
Q=3,\quad e^{c}=\begin{cases}3\\\ \hat{3}^{\prime}\end{cases},\quad L=1\oplus
1\oplus 1,$ (18)

where the top (bottom) case for $d^{c}$ and $e^{c}$ realizes $p=1$ ($p=0$).
There are equivalent combinations which give the same Yukawa couplings, e.g.
$u^{c}=1^{\prime}\oplus 1^{\prime}\oplus\hat{1}$ and $Q=3^{\prime}$ for the up
Yukawa couplings. The triplets can be split into a doublet and a singlet, or
three singlets, but we do not consider these possibilities because these are
less predictive than the case of the triplets. We consider the different
texture of the quarks from the model proposed in our recent analysis [65]. In
the current paper, we allow the co-existence of the hatted and non-hatted
representations, in a same type of fermions as assigned in $u^{c}$. This
generalization allows $V_{tb}\sim\mathcal{O}\left({1}\right)$ even if $Q$ is a
triplet, unlike the model without the co-existence.

Table 2:  Assignments of the fermions and Higgs doublets under $G_{\mathrm{EW}}:=SU(2)_{L}\times U(1)_{Y}$, $S_{4}^{\prime}$ and the weight $k$.  | $u^{c}_{i}$ | $d^{c}_{i}$ | $Q$ | $e^{c}$ | $L_{i}$ | $H_{u}$ | $H_{d}$
---|---|---|---|---|---|---|---
$G_{\mathrm{EW}}$ | $1_{-2/3}$ | $1_{1/3}$ | $2_{1/6}$ | $1_{1}$ | $2_{-1/2}$ | $2_{1/2}$ | $2_{-1/2}$
$S_{4}^{\prime}$ | $(1,1,\hat{1}^{\prime})$ |  $(1,1,1)$ or $(\hat{1}^{\prime},\hat{1}^{\prime},\hat{1}^{\prime})$ | $3$ |  $3$ or $\hat{3}^{\prime}$ | $(1,1,1)$ | $1$ | $1$
$k$ | $-k_{u_{i}}$ | $-k_{d_{i}}$ | $-k_{Q}$ | $-k_{e}$ | $-k_{L_{i}}$ | $0$ | 0

We study the models with representations of the quarks and leptons shown in
Eq. (18), so that the texture in Eqs. (16) and (17) is realized. The matter
contents of the quarks and leptons and our notation of the modular weights are
summarized in Table 2. The Yukawa coupling terms in the superpotential are
schematically given by

$\displaystyle W=$ $\displaystyle\
H_{u}\left\\{\sum_{a=1}^{2}\alpha_{a}\left(QY_{3}^{({k_{u_{a}}+k_{Q}})}\right)_{1}u_{a}^{c}+\alpha_{3}\left(QY_{\hat{3}}^{(k_{u_{3}}+k_{Q})}u_{3}^{c}\right)_{1}\right\\}$
(19)
$\displaystyle\quad+H_{d}\sum_{i=1}^{3}\left\\{\beta_{i}\left(QY^{(k_{d_{i}}+k_{Q})}_{\mathbf{3}}d^{c}_{i}\right)_{1}+\gamma_{i}L_{i}\left(Y^{(k_{e}+k_{L_{i}})}_{\mathbf{3}}e^{c}\right)_{1}\right\\}+\sum_{i,j=1}^{3}\frac{c_{ij}}{\Lambda}Y_{{1}}^{(k_{L_{i}}+k_{L_{j}})}L_{i}H_{u}L_{j}H_{u}$
$\displaystyle=:$ $\displaystyle\
H_{u}QY_{u}u^{c}+H_{d}QY_{d}d^{c}+H_{d}LY_{e}e^{c}+H_{u}L^{T}C_{n}H_{u}L,$

where $\mathbf{3}=3$ ($\hat{3}$) for $p=1$ ($p=0$) in the Yukawa terms to
which $H_{d}$ couples. Here, $(\cdots)_{1}$ is the trivial singlet combination
of the product inside the parenthesis. The number of the coefficients
$\alpha_{i}$, $\beta_{i}$ and $\gamma_{i}$ will change accordingly to the
number of modular forms for a given representation $r$ and modular weight $k$,
as explicitly shown after assigning modular weights. We consider the Weinberg
operator for the Majorana neutrino masses given by the last term. We assume
that $L$ is a non-hatted singlet such that the corresponding modular form is a
trivial-singlet. The modular forms in the Weinberg operator are the non-
trivial singlets $1^{\prime}$ if $L$ are the hatted singlets, but we do not
consider this possibility for simplicity.

The Kähler potential of the chiral superfield $f$ with a weight $-k_{f}$ is
given by

$\displaystyle
K\supset\frac{f^{\dagger}f}{(-i\tau+i\overline{\tau})^{k_{f}}},$ (20)

and hence the couplings are modified by the canonical normalization as

$\displaystyle\left[Y_{f}\right]_{ij}\to\left(\sqrt{2\,\mathrm{Im}\,\tau}\right)^{k_{Y}}\left[Y_{f}\right]_{ij},\quad\left[C_{n}\right]_{ij}\to\left(\sqrt{2\,\mathrm{Im}\,\tau}\right)^{k_{L_{i}}+k_{L_{j}}}\left[C_{n}\right]_{ij},$
(21)

for $i,j=1,2,3$ and $f=u,d,e$. Here, $k_{Y}$ is the weight of the Yukawa
coupling.

The hierarchical structures of the Yukawa matrices before the canonical
normalization are given by

$\displaystyle Y_{u}\sim\begin{pmatrix}\epsilon^{2}&\epsilon^{2}&\epsilon\\\
\epsilon&\epsilon&1\\\ \epsilon^{3}&\epsilon^{3}&\epsilon^{2}\\\
\end{pmatrix},\quad Y_{d}\sim
Y_{e}^{T}\sim\begin{pmatrix}\epsilon^{2}&\epsilon^{2}&\epsilon^{2}\\\
\epsilon&\epsilon&\epsilon\\\ \epsilon^{3}&\epsilon^{3}&\epsilon^{3}\\\
\end{pmatrix},\quad C_{n}\sim\frac{1}{\Lambda}\begin{pmatrix}1&1&1\\\ 1&1&1\\\
1&1&1\end{pmatrix}.$ (22)

These structures realize hierarchical structures in Eqs. (16) and (17). The CP
phase from the modulus VEV appears only in $\epsilon(\tau)$ up to
$\mathcal{O}\left({\epsilon^{4}}\right)$, but the phases can be absorbed by
the redefinition of the fermions. Thus there should be CP phases from
$\mathcal{O}\left({1}\right)$ coefficients rather than that from $\tau$. This
fact also implies that the Yukawa couplings are independent of
$\mathrm{Re}\,\tau$ as a good approximation.

###  3.2 Model for $p=1$

Table 3:  Values of the masses and mixing angles at the fitted points in the
models in the case of $p=1$ (left) and $p=0$ (right). The second, third and
fourth columns show the model predictions, experimental central values and its
errors, respectively. See the text in detail. $s_{ij}^{Q}$ is the mixing angle
between the $i$-th and $j$-th quark in the CKM matrix and $s_{ij}^{2}$ is the
square of mixing angle in the PMNS matrix. $\delta_{\mathrm{CKM}}$ and
$\delta_{\mathrm{PMNS}}$ are the CP phases.
$R^{21}_{32}:=(m_{2}^{2}-m_{1}^{2})/(m_{3}^{2}-m_{2}^{2})$ is the ratio of the
neutrino mass squared difference. These are defined in the standard
parametrization [71].

obs. | value | center | error
---|---|---|---
$y_{u}$$/10^{-6}$ | 4.44 | 2.85 | 0.88
$y_{c}$$/10^{-3}$ | 1.481 | 1.479 | 0.052
$y_{t}$ | 0.5322 | 0.5320 | 0.0053
$y_{d}$$/10^{-5}$ | 1.94 | 1.93 | 0.21
$y_{s}$$/10^{-4}$ | 3.88 | 3.82 | 0.21
$y_{b}$$/10^{-2}$ | 2.097 | 2.100 | 0.021
$y_{e}$$/10^{-6}$ | 7.816 | 7.816 | 0.047
$y_{\mu}$$/10^{-3}$ | 1.6496 | 1.6500 | 0.0099
$y_{\tau}$$/10^{-2}$ | 2.808 | 2.805 | 0.028
$s_{12}^{Q}$ | 0.22520 | 0.22541 | 0.00072
$s_{23}^{Q}$$/10^{-2}$ | 4.007 | 3.998 | 0.064
$s_{13}^{Q}$$/10^{-3}$ | 3.43 | 3.48 | 0.13
$\delta_{\mathrm{CKM}}$ | 1.2395 | 1.2080 | 0.0540
$R^{21}_{32}$$/10^{-2}$ | 3.053 | 3.070 | 0.084
$s_{12}^{2}$ | 0.302 | 0.307 | 0.013
$s_{23}^{2}$ | 0.547 | 0.546 | 0.021
$s_{13}^{2}$$/10^{-2}$ | 2.203 | 2.200 | 0.070
$\delta_{\mathrm{PMNS}}$ | -0.85 | -2.01 | 0.63

obs. | value | center | error
---|---|---|---
$y_{u}$$/10^{-6}$ | 3.27 | 2.74 | 0.85
$y_{c}$$/10^{-3}$ | 1.418 | 1.419 | 0.050
$y_{t}$ | 0.5030 | 0.5029 | 0.0050
$y_{d}$$/10^{-5}$ | 6.93 | 7.32 | 0.81
$y_{s}$$/10^{-3}$ | 1.383 | 1.450 | 0.078
$y_{b}$$/10^{-2}$ | 7.997 | 7.976 | 0.080
$y_{e}$$/10^{-5}$ | 2.966 | 2.966 | 0.018
$y_{\mu}$$/10^{-3}$ | 6.259 | 6.261 | 0.038
$y_{\tau}$ | 0.1077 | 0.1074 | 0.0011
$s_{12}^{Q}$ | 0.22533 | 0.22541 | 0.00072
$s_{23}^{Q}$$/10^{-2}$ | 4.003 | 4.011 | 0.064
$s_{13}^{Q}$$/10^{-3}$ | 3.47 | 3.49 | 0.13
$\delta_{\mathrm{CKM}}$ | 1.2032 | 1.2080 | 0.0540
$R^{21}_{32}$$/10^{-2}$ | 3.070 | 3.070 | 0.084
$s_{12}^{2}$ | 0.307 | 0.307 | 0.013
$s_{23}^{2}$ | 0.541 | 0.546 | 0.021
$s_{13}^{2}$$/10^{-2}$ | 2.199 | 2.200 | 0.070
$\delta_{\mathrm{PMNS}}$ | -2.30 | -2.01 | 0.63

Now we consider the model for $p=1$. The modular weights are given by

$\displaystyle(k_{u_{1}},k_{u_{2}},k_{u_{3}})=(0,4,3),\quad(k_{d_{1}},k_{d_{2}},k_{d_{3}})=(0,2,4),\quad
k_{Q}=4,\quad k_{e}=4,\quad k_{L}=(0,2,4),$ (23)

so that the superpotential is given by

$\displaystyle W=$ $\displaystyle\
H_{u}\left\\{\alpha_{1}\left(QY_{3}^{(4)}\right)_{1}u_{1}^{c}+\alpha_{2}^{i_{Y}}\left(QY_{3}^{i_{Y}(8)}\right)_{1}u_{2}^{c}+\alpha_{3}^{i_{Y}}\left(QY_{\hat{3}}^{i_{Y}(7)}u_{3}^{c}\right)_{1}\right\\}$
(24)
$\displaystyle+H_{d}\left\\{\sum_{a=1}^{2}\beta_{a}\left(QY^{(2+2a)}_{{3}}\right)_{1}d^{c}_{a}+\beta_{3}^{i_{Y}}\left(QY^{i_{Y}(8)}_{{3}}\right)_{1}d^{c}_{3}\right\\}$
$\displaystyle+H_{d}\left\\{\sum_{a=1}^{2}\gamma_{a}L_{a}\left(Y^{(2+2a)}_{{3}}e^{c}\right)_{1}+\gamma_{3}^{i_{Y}}L_{3}\left(Y^{i_{Y}(8)}_{{3}}e^{c}\right)_{1}\right\\}+\sum_{i,j=1}^{3}\frac{c_{ij}}{\Lambda}Y^{(2i+2j-4)}_{{1}}L_{i}H_{u}L_{j}H_{u},$

where the summation over the modular forms $i_{Y}=1,2$ are implicit. The
modular weights of $d^{c}_{i}$, as well as $L_{i}$ should be different among
the flavors to make all of the fermions massive. For instance, the Yukawa
matrix $Y_{d}$ has only one non-zero singular value if all of the modular
weights are $4$ or $6$, where there is only one $3$ representation. Hence, the
assignment in Eq. (23) is the minimal number to have three massive down quarks
and charged leptons. We assume that the 18 coefficients are real except
$\alpha_{3}^{1}$ and $\gamma_{3}^{1}$ to realize the non-zero CP phases in the
CKM and PMNS matrices. Note that $c_{12}$ for Weinberg operator is irrelevant,
because there is no trivial singlet at $k=2$. Since the absolute size of the
neutrino mass can be fitted by $\Lambda$, we only fit to the ratio of the
neutrino mass squared difference
$R^{21}_{32}:=(m_{2}^{2}-m_{1}^{2})/(m_{3}^{2}-m_{2}^{2})$.

We fitted $\mathrm{Im}\,\tau$,
$\tan\beta:={\langle{H_{u}^{0}}\rangle}/{\langle{H_{d}^{0}}\rangle}$ and the
coefficients to explain the Yukawa couplings and the CKM angles at the GUT
scale in the MSSM with $M_{\mathrm{SUSY}}=10~{}\mathrm{TeV}$ and neglecting
the threshold corrections [72]. We assume the $1\%$ relative uncertainties for
the third generation fermions, and those shown in Table 3 of Ref. [72] are
used for the other observables. The neutrino data is taken from Ref. [71]
under the assumption of the normal ordering, as $m_{1}<m_{2}<m_{3}$ is
naturally predicted from the weight assignment.

With the weight assignment, the hierarchical structures with the powers of
$t:=2\mathrm{Im}\,\tau$ is given by

$\displaystyle(y_{u},y_{c},y_{t})\sim$ $\displaystyle\
\left(\epsilon^{3}/t^{3/2},\epsilon^{2}t^{1/2},1\right)y_{t}\sim(5\times
10^{-7},7\times 10^{-4},0.5),$ (25)
$\displaystyle(y_{d},y_{s},y_{b})\sim(y_{e},y_{\mu},y_{\tau})\sim$
$\displaystyle\ \left(\epsilon^{3}/t^{3/2},\epsilon^{2}/t^{1/2},\epsilon
t^{1/2}\right)y_{t}\sim(5\times 10^{-7},1\times 10^{-4},0.03),$
$\displaystyle(s_{12}^{Q},s_{23}^{Q},s_{13}^{Q})\sim$ $\displaystyle\
\left(\epsilon,\epsilon,\epsilon^{2}\right)\sim(0.02,0.02,0.0006),$
$\displaystyle(R^{21}_{32},s_{12}^{2},s_{23}^{2},s_{13}^{2})\sim$
$\displaystyle\
\left(1/t^{4},1/t^{2},1/t^{2},1/t^{4}\right)\sim\left(0.001,0.03,0.03,0.001\right),$

where $\mathrm{Im}\,\tau=2.83$ and $y_{t}=0.5$ are used for the numerical
estimation. Here, $y_{f}$, where $f$ being the SM charged fermions, are the
Yukawa couplings in the mass basis. Since the top Yukawa coupling $y_{t}$ is
roughly given by $\alpha_{3}^{1}t^{7/2}$, the coefficient $\alpha_{3}^{1}$ as
well as the other coefficients should be so small that the factors from powers
of $t$ are canceled. $s_{ij}^{Q}$ is the mixing angle between the $i$-th and
$j$-th quark in the CKM matrix and $s_{ij}^{2}$ is the square of mixing angle
in the PMNS matrix in the standard parametrization [71]. Among the
observables, some of the light fermion Yukawa couplings $y_{d}$, $y_{e}$,
$y_{\mu}$, the Cabbibo angle $s_{12}^{Q}$, and the neutrino mixing angles are
smaller than the observed values by $\mathcal{O}\left({10}\right)$. Most of
the differences are explained by the modular forms, e.g. the factor
$7\sqrt{2}$ in the coefficients of the element with
$\mathcal{O}\left({\epsilon^{3}}\right)$ in $Y^{2(7)}_{\hat{3}}$. Whereas,
$\left|{c_{13}/c_{33}}\right|\sim t$ will be necessary to compensate the
suppression by $\mathcal{O}\left({t}\right)$ in the neutrino mixing. Since all
of the values are predicted to be smaller than the observed values, there will
be no cancellation among the different contributions. Thus, in the current
model, the cancellation required to explain the charm mass in the model
studied in Ref. [65] is absent.

We found the parameter sets which can explain the observed data within
$1.8\sigma$. At the benchmark point, $\tan\beta=3.6666$,
$\mathrm{Im}\,\tau=2.8258$, $|\alpha_{3}^{1}|=1.2334\times 10^{-3}$,

$\displaystyle\frac{1}{|\alpha_{3}^{1}|}\begin{pmatrix}\alpha_{1}\\\
\alpha^{1}_{2}\\\ \alpha^{2}_{2}\\\ \alpha^{1}_{3}\\\
\alpha^{2}_{3}\end{pmatrix}=$ $\displaystyle\ \begin{pmatrix}-0.6863\\\
-1.7953\\\ 0.8360\\\ e^{-1.6360i}\\\
-1.1016\end{pmatrix},\quad\frac{1}{|\alpha_{3}^{1}|}\begin{pmatrix}\beta_{1}\\\
\beta_{2}\\\ \beta^{1}_{3}\\\
\beta^{2}_{3}\end{pmatrix}=\begin{pmatrix}1.6470\\\ -1.8393\\\ 0.6233\\\
-0.6245\end{pmatrix},$ (26)
$\displaystyle\frac{1}{|\alpha_{3}^{1}|}\begin{pmatrix}\gamma_{1}\\\
\gamma_{2}\\\ \gamma^{1}_{3}\\\ \gamma^{2}_{3}\end{pmatrix}=$ $\displaystyle\
\begin{pmatrix}-1.6692\\\ -2.7814\\\ 0.8032\times e^{0.9284i}\\\
2.7497\end{pmatrix},\quad\frac{1}{|\alpha_{3}^{1}|}\begin{pmatrix}c_{11}\\\
c_{22}\\\ c_{33}\\\ c_{13}\\\ c_{23}\end{pmatrix}=\begin{pmatrix}1.7273\\\
-1.9311\\\ -0.6692\\\ -5.3871\\\ -2.0488\end{pmatrix}.$

The absolute values of the coefficients are in the range $[0.62,5.4]$. As
expected from the estimation Eq. (25), $\left|{c_{13}}\right|$ is larger than
the other coefficients to realize the neutrino mixing angles. This factor
would be explained by numerical factors of modular forms of non-singlet
representations when we consider the see-saw mechanism as a UV completion of
the Weinberg operator.

###  3.3 Model for $p=0$

Now we consider the model with $p=0$. The modular weights are given by

$\displaystyle(k_{u_{1}},k_{u_{2}},k_{u_{3}})=(4,8,7),\quad(k_{d_{1}},k_{d_{2}},k_{d_{3}})=(1,3,5),\quad
k_{Q}=0,\quad k_{e}=1,\quad k_{L}=(2,4,4),$ (27)

so that the superpotential is given by

$\displaystyle W=$ $\displaystyle\
H_{u}\left\\{\alpha_{1}\left(QY_{3}^{(4)}\right)_{1}u_{1}^{c}+\alpha_{2}^{i_{Y}}\left(QY_{3}^{i_{Y}(8)}\right)_{1}u_{2}^{c}+\alpha_{3}^{i_{Y}}\left(QY_{\hat{3}}^{i_{Y}(7)}u_{3}^{c}\right)_{1}\right\\}$
(28)
$\displaystyle+H_{d}\left\\{\sum_{a=1}^{2}\beta_{a}\left(QY^{(1+2a)}_{{\hat{3}}}d^{c}_{a}\right)_{1}+\beta_{3}^{i_{Y}}\left(QY^{i_{Y}(5)}_{{\hat{3}}}d^{c}_{3}\right)_{1}\right\\}$
$\displaystyle+H_{d}\left\\{\gamma_{1}L_{1}\left(Y^{(3)}_{\hat{3}}e^{c}\right)_{1}+\sum_{b=2}^{3}\gamma_{b}^{i_{Y}}L_{b}\left(Y^{i_{Y}(5)}_{\hat{3}}e^{c}\right)_{1}\right\\}+\sum_{i,j=1}^{3}\frac{c_{ij}}{\Lambda}Y^{(8-2\delta_{i1}-2\delta_{j1})}_{{1}}L_{i}H_{u}L_{j}H_{u}.$

We assigned the same modular weights for $L_{2}$ and $L_{3}$, so that there is
no hierarchy in the second and third generations in the neutrino sector. The
three charged leptons are still massive because there are two modular forms of
the $\hat{3}$ representation at $k=5$. There are 20 coefficients in the model.
The hierarchical structure of the observables are estimated as

$\displaystyle(y_{u},y_{c},y_{t})\sim$ $\displaystyle\
\left(\epsilon^{3}/t^{3/2},\epsilon^{2}t^{1/2},1\right)y_{t}\sim(1\times
10^{-7},3\times 10^{-4},0.5),$ (29) $\displaystyle(y_{d},y_{s},y_{b})\sim$
$\displaystyle\
\left(\epsilon^{2}/t^{3},\epsilon/t^{2},1/t\right)y_{t}\sim(6\times
10^{-7},2\times 10^{-4},0.08),$ $\displaystyle(y_{e},y_{\mu},y_{\tau})\sim$
$\displaystyle\
\left(\epsilon^{2}/t^{2},\epsilon/t,1/t\right)y_{t}\sim(3\times
10^{-6},0.001,0.08),$ $\displaystyle(s_{12}^{Q},s_{23}^{Q},s_{13}^{Q})\sim$
$\displaystyle\
\left(\epsilon,\epsilon,\epsilon^{2}\right)\sim(0.02,0.02,3\times 10^{-4}),$
$\displaystyle(R^{21}_{32},s_{12}^{2},s_{23}^{2},s_{13}^{2})\sim$
$\displaystyle\ \left(1,1/t^{2},1,1/t^{2}\right)\sim(1,0.03,1,0.03),$

where $\mathrm{Im}\,\tau=3.07$. In this case, $y_{u}$, $y_{d}$, $s_{12}^{Q}$,
$s_{13}^{Q}$ and $s_{12}^{2}$ are more than 10 times smaller than the observed
values, especially $y_{d}$ is about $125$ times smaller. Again, the gaps in
the quark sector can be explained by the modular forms, e.g. the factors of
$5$ in the $\mathcal{O}\left({\epsilon^{2}}\right)$ elements in
$Y^{(3)}_{\hat{3}}$ and $Y^{2(5)}_{\hat{3}}$. Thus the observed values will be
explained by the $\mathcal{O}\left({1}\right)$ coefficients. We assume that
the neutrinos are in the normal ordering, because the model predicts
$m_{1}\lesssim m_{2},m_{3}$ and $s_{12}^{2},s_{13}^{2}\lesssim s_{23}^{2}$
which more easily fits to it rather than the inverse ordering case.

We found the benchmark point which can explain the experimental values within
$0.9\sigma$. The values are given by $\tan\beta=14.1755$,
$\mathrm{Im}\,\tau=3.0744$, $|\alpha_{3}^{1}|=8.6886\times 10^{-4}$,

$\displaystyle\frac{1}{|\alpha_{3}^{1}|}\begin{pmatrix}\alpha_{1}\\\
\alpha^{1}_{2}\\\ \alpha^{2}_{2}\\\ \alpha^{1}_{3}\\\
\alpha^{2}_{3}\end{pmatrix}=$ $\displaystyle\ \begin{pmatrix}-1.0336\\\
1.2757\\\ -2.2480\\\ e^{0.5391i}\\\
4.8979\end{pmatrix},\quad\frac{1}{|\alpha_{3}^{1}|}\begin{pmatrix}\beta_{1}\\\
\beta_{2}\\\ \beta^{1}_{3}\\\
\beta^{2}_{3}\end{pmatrix}=\begin{pmatrix}2.5696\\\ -3.3719\\\ 4.0290\\\
-0.8087\end{pmatrix},$ (30)
$\displaystyle\frac{1}{|\alpha_{3}^{1}|}\begin{pmatrix}\gamma_{1}\\\
\gamma^{1}_{2}\\\ \gamma^{2}_{2}\\\ \gamma^{1}_{3}\\\
\gamma^{2}_{3}\end{pmatrix}=$ $\displaystyle\ \begin{pmatrix}4.0267\\\
4.6090\\\ 0.8186\\\ 4.4630\times e^{-2.9326i}\\\
-0.8028\end{pmatrix},\quad\frac{1}{|\alpha_{3}^{1}|}\begin{pmatrix}c_{11}\\\
c_{22}\\\ c_{33}\\\ c_{12}\\\ c_{13}\\\
c_{23}\end{pmatrix}=\begin{pmatrix}-1.3694\\\ 1.2175\\\ 1.3739\\\ 1.4887\\\
-6.0655\\\ -1.1843\end{pmatrix}.$

The absolute values of the coefficients are within the range $[0.8,6.1]$, and
hence the ratio of the smallest and largest coefficients is $7.6$. The
coefficient $c_{13}$ is $\mathcal{O}\left({t}\right)$ larger than the other
coefficients of the Weinberg operator to obtain the large mixing angles. This
factor would be explained in the see-saw mechanism as in the $p=1$ model.

##  4 Conclusion

In this work, we constructed models to explain the hierarchical structures of
both quark and lepton sectors under the modular flavor symmetry
$\Gamma^{\prime}_{4}$, which is isomorphic to $S_{4}^{\prime}$. The
hierarchical structure is realized by the FN mechanism due to the residual
$\mathbbm{Z}_{4}^{T}$ symmetry, where the small parameter is given by
$\epsilon(\tau)$ having the charge one under $\mathbbm{Z}_{4}^{T}$. In
addition, another hierarchical structure is induced from the powers of
$2\mathrm{Im}\,\tau$, originated from the canonical normalization depending on
the modular weights. Unlike the model proposed in our previous paper [65], we
consider the case in which the non-hatted and hatted representations appear in
the same type of fermions, as in $u^{c}$ of Eq. (18). This generalization
allows us to explain the fermion hierarchies with the
$\mathcal{O}\left({1}\right)$ coefficients by the $S_{4}^{\prime}$ symmetry
alone. In fact, the ratio of the smallest and largest coefficients are 8.6
(7.6) in the first (second) model. It is interesting that the small deviations
from the hierarchical structure by $\epsilon$ in Eqs. (16) and (17) are
resolved by the powers of $2\mathrm{Im}\,\tau$ and the numerical factors in
the modular forms.

The lepton sector is studied in addition to the quark sector, where we assume
that the Majorana neutrino masses are realized from the Weinberg operator. The
weak doublet lepton $L$, including the neutrinos, should have a common charge
under the $\mathbbm{Z}^{T}_{4}$ symmetry, so that the lepton mixing angles are
$\mathcal{O}\left({1}\right)$. This indicates that the three weak doublet
leptons are the same singlet representation. The small hierarchy is induced
from the powers of $2\mathrm{Im}\,\tau$, since the weights of $L$ may be
different for three massive charged leptons. This naturally explains the
smallness of the angle $s_{13}^{2}$, while it also predicts small angle
$s_{12}^{2}$, see Eqs. (25) and (29). These small differences are simply
resolved by the difference in the $\mathcal{O}\left({1}\right)$ coefficients.

## Acknowledgment

The work of J.K. is supported in part by the Institute for Basic Science
(IBS-R018-D1). This work is supported in part by he Grant-in-Aid for
Scientific Research from the Ministry of Education, Science, Sports and
Culture (MEXT), Japan No. 22K03601 (T.H.) and 18K13534 (J.K.). The work of
Y.A. is supported by JSPS Overseas Research Fellowships.

##  Appendix A Modular forms of $\Gamma_{4}^{\prime}$

We collect the modular forms used in the model, which can be constructed from
products of the modular form at $k=1$ in Eq. (12). The products of the non-
trivial singlets are given by

$\displaystyle 1^{\prime}\otimes
1^{\prime}=\hat{1}\otimes\hat{1}^{\prime}=1,\quad\hat{1}\otimes\hat{1}=\hat{1}^{\prime}\otimes\hat{1}^{\prime}=1^{\prime},\quad
1^{\prime}\otimes\hat{1}^{\prime}=\hat{1},\quad
1^{\prime}\otimes\hat{1}=\hat{1}^{\prime}.$ (31)

The products of $2$ and $3$ are given by

$\displaystyle 2(u)\otimes 2(v)=$ $\displaystyle\
(u_{1}v_{1}+u_{2}v_{2})_{1}\oplus(u_{1}v_{2}-u_{2}v_{1})_{1^{\prime}}\oplus\begin{pmatrix}u_{2}v_{2}-u_{1}v_{1}\\\
u_{1}v_{2}+u_{2}v_{1}\end{pmatrix}_{2},$ (32) $\displaystyle 2(u)\otimes
3(\phi)=$ $\displaystyle\ \begin{pmatrix}-2u_{1}\phi_{1}\\\
u_{1}\phi_{2}-\sqrt{3}u_{2}\phi_{3}\\\ u_{1}\phi_{3}-\sqrt{3}u_{2}\phi_{2}\\\
\end{pmatrix}_{3}\oplus\begin{pmatrix}-2u_{2}\phi_{1}\\\
u_{2}\phi_{2}+\sqrt{3}u_{1}\phi_{3}\\\ u_{2}\phi_{3}+\sqrt{3}u_{1}\phi_{2}\\\
\end{pmatrix}_{3^{\prime}},$ $\displaystyle 3(\phi)\otimes 3(\psi)=$
$\displaystyle\
\left(\phi_{1}\psi_{1}+\phi_{2}\psi_{3}+\phi_{3}\psi_{2}\right)_{1}\oplus\begin{pmatrix}2\phi_{1}\psi_{1}-\phi_{2}\psi_{3}-\phi_{3}\psi_{2}\\\
\sqrt{3}\left(\phi_{2}\psi_{2}+\phi_{3}\psi_{3}\right)\end{pmatrix}_{2}$
$\displaystyle\quad\oplus\begin{pmatrix}\phi_{2}\psi_{2}-\phi_{3}\psi_{3}\\\
-\phi_{3}\psi_{1}-\phi_{1}\psi_{3}\\\
\phi_{1}\psi_{2}+\phi_{2}\psi_{1}\end{pmatrix}_{3}\oplus\begin{pmatrix}\phi_{2}\psi_{3}-\phi_{3}\psi_{2}\\\
\phi_{1}\psi_{2}-\phi_{2}\psi_{1}\\\
\phi_{3}\psi_{1}-\phi_{1}\psi_{3}\end{pmatrix}_{3^{\prime}}.$

The products of the other representations with prime and/or hat are calculated
in the same manner as the singlets in Eq. (31), with noting that primed
doublets, $2^{\prime}$ and $\hat{2}^{\prime}$ denoted by
$(u_{1},u_{2})_{\mathbf{2}^{\prime}}$, should be understood as
$(u_{2},-u_{1})_{\mathbf{2}}$, where $\mathbf{2}=2,\hat{2}$.

The modular forms used in the models are given by

$\displaystyle Y^{(3)}_{\hat{3}}=$ $\displaystyle\
\begin{pmatrix}2\sqrt{2}\epsilon\theta\left(\epsilon^{4}+\theta^{4}\right)\\\
-\epsilon^{2}(\epsilon^{4}-5\theta^{4})\\\
-\theta^{2}(5\epsilon^{4}-\theta^{4})\end{pmatrix},$ (33) $\displaystyle
Y^{(4)}_{1}=$ $\displaystyle\
\epsilon^{8}+14\epsilon^{4}\theta^{4}+\theta^{8},$ $\displaystyle\quad
Y^{(4)}_{3}=$ $\displaystyle\
\epsilon\theta\left(\epsilon^{4}-\theta^{4}\right)\begin{pmatrix}-\sqrt{2}\epsilon\theta\\\
-\epsilon^{2}\\\ \theta^{2}\end{pmatrix},$ $\displaystyle Y^{1(5)}_{\hat{3}}=$
$\displaystyle\
\epsilon\theta(\epsilon^{4}-\theta^{4})\begin{pmatrix}\epsilon^{4}-\theta^{4}\\\
2\sqrt{2}\epsilon\theta^{3}\\\ 2\sqrt{2}\epsilon^{3}\theta\end{pmatrix},$
$\displaystyle\quad Y^{2(5)}_{\hat{3}}=$ $\displaystyle\
\begin{pmatrix}16\sqrt{2}\epsilon^{5}\theta^{5}\\\
\epsilon^{2}(\epsilon^{8}+10\epsilon^{4}\theta^{4}+5\theta^{8})\\\
-\theta^{2}(5\epsilon^{8}+10\epsilon^{4}\theta^{4}+\theta^{8})\end{pmatrix},$
$\displaystyle Y^{(6)}_{1}=$ $\displaystyle\
(\epsilon^{4}+\theta^{4})(\epsilon^{8}-34\epsilon^{4}\theta^{4}+\theta^{8}),$
$\displaystyle\quad Y^{(6)}_{3}=$ $\displaystyle\
\epsilon\theta\left(\epsilon^{4}-\theta^{4}\right)\begin{pmatrix}-2\sqrt{2}\epsilon\theta(\epsilon^{4}+\theta^{4})\\\
\epsilon^{2}\left(\epsilon^{4}-5\theta^{4}\right)\\\
\theta^{2}\left(5\epsilon^{4}-\theta^{4}\right)\end{pmatrix},$ $\displaystyle
Y^{1(7)}_{\hat{3}}=$ $\displaystyle\
\begin{pmatrix}32\sqrt{2}\epsilon^{5}\theta^{5}\left(\epsilon^{4}+\theta^{4}\right)\\\
-\epsilon^{2}(\epsilon^{12}-19\epsilon^{8}\theta^{4}-45\epsilon^{4}\theta^{8}-\theta^{12})\\\
-\theta^{2}(\epsilon^{12}+45\epsilon^{8}\theta^{4}+19\epsilon^{4}\theta^{8}-\theta^{12})\end{pmatrix},$
$\displaystyle\quad Y^{2(7)}_{\hat{3}}=$ $\displaystyle\
\epsilon\theta(\epsilon^{4}-\theta^{4})\begin{pmatrix}\epsilon^{8}-\theta^{8}\\\
-\sqrt{2}\epsilon\theta^{3}(7\epsilon^{4}+\theta^{4})\\\
-\sqrt{2}\epsilon^{3}\theta(\epsilon^{4}+7\theta^{4})\\\ \end{pmatrix},$
$\displaystyle Y^{(8)}_{1}=$ $\displaystyle\
(\epsilon^{8}+14\epsilon^{4}\theta^{4}+\theta^{8})^{2},$ $\displaystyle
Y^{1(8)}_{3}=$ $\displaystyle\
\epsilon\theta(\epsilon^{4}-\theta^{4})\begin{pmatrix}16\sqrt{2}\epsilon^{5}\theta^{5}\\\
\epsilon^{2}\left(\epsilon^{8}+10\epsilon^{4}\theta^{4}+5\theta^{8}\right)\\\
-\theta^{2}\left(5\epsilon^{8}+10\epsilon^{4}\theta^{4}+\theta^{8}\right)\end{pmatrix},\quad$
$\displaystyle Y^{2(8)}_{3}=$ $\displaystyle\
\epsilon^{2}\theta^{2}(\epsilon^{4}-\theta^{4})^{2}\begin{pmatrix}\epsilon^{4}-\theta^{4}\\\
2\sqrt{2}\epsilon\theta^{3}\\\ 2\sqrt{2}\epsilon^{3}\theta\\\ \end{pmatrix}.$

These are normalized such that absolute value of the coefficient of the
largest element at $\epsilon\ll 1$ is unity.

## References

  * [1] F. Feruglio, Are neutrino masses modular forms?, pp. 227–266.  2019\.  arXiv:1706.08749.
  * [2] T. Kobayashi, K. Tanaka, and T. H. Tatsuishi, Neutrino mixing from finite modular groups, Phys. Rev. D 98 (2018), no. 1 016004, [arXiv:1803.10391].
  * [3] J. T. Penedo and S. T. Petcov, Lepton Masses and Mixing from Modular $S_{4}$ Symmetry, Nucl. Phys. B 939 (2019) 292–307, [arXiv:1806.11040].
  * [4] P. P. Novichkov, J. T. Penedo, S. T. Petcov, and A. V. Titov, Modular A5 symmetry for flavour model building, JHEP 04 (2019) 174, [arXiv:1812.02158].
  * [5] G.-J. Ding, S. F. King, and X.-G. Liu, Neutrino mass and mixing with $A_{5}$ modular symmetry, Phys. Rev. D 100 (2019), no. 11 115005, [arXiv:1903.12588].
  * [6] X.-G. Liu and G.-J. Ding, Neutrino Masses and Mixing from Double Covering of Finite Modular Groups, JHEP 08 (2019) 134, [arXiv:1907.01488].
  * [7] P. P. Novichkov, J. T. Penedo, and S. T. Petcov, Double cover of modular $S_{4}$ for flavour model building, Nucl. Phys. B 963 (2021) 115301, [arXiv:2006.03058].
  * [8] X.-G. Liu, C.-Y. Yao, and G.-J. Ding, Modular invariant quark and lepton models in double covering of $S_{4}$ modular group, Phys. Rev. D 103 (2021), no. 5 056013, [arXiv:2006.10722].
  * [9] X.-G. Liu, C.-Y. Yao, B.-Y. Qu, and G.-J. Ding, Half-integral weight modular forms and application to neutrino mass models, Phys. Rev. D 102 (2020), no. 11 115035, [arXiv:2007.13706].
  * [10] R. de Adelhart Toorop, F. Feruglio, and C. Hagedorn, Finite Modular Groups and Lepton Mixing, Nucl. Phys. B 858 (2012) 437–467, [arXiv:1112.1340].
  * [11] G. Altarelli and F. Feruglio, Discrete Flavor Symmetries and Models of Neutrino Mixing, Rev. Mod. Phys. 82 (2010) 2701–2729, [arXiv:1002.0211].
  * [12] H. Ishimori, T. Kobayashi, H. Ohki, Y. Shimizu, H. Okada, and M. Tanimoto, Non-Abelian Discrete Symmetries in Particle Physics, Prog. Theor. Phys. Suppl. 183 (2010) 1–163, [arXiv:1003.3552].
  * [13] H. Ishimori, T. Kobayashi, H. Ohki, H. Okada, Y. Shimizu, and M. Tanimoto, An introduction to non-Abelian discrete symmetries for particle physicists, vol. 858.  2012\.
  * [14] D. Hernandez and A. Y. Smirnov, Lepton mixing and discrete symmetries, Phys. Rev. D 86 (2012) 053014, [arXiv:1204.0445].
  * [15] S. F. King and C. Luhn, Neutrino Mass and Mixing with Discrete Symmetry, Rept. Prog. Phys. 76 (2013) 056201, [arXiv:1301.1340].
  * [16] S. F. King, A. Merle, S. Morisi, Y. Shimizu, and M. Tanimoto, Neutrino Mass and Mixing: from Theory to Experiment, New J. Phys. 16 (2014) 045018, [arXiv:1402.4271].
  * [17] M. Tanimoto, Neutrinos and flavor symmetries, AIP Conf. Proc. 1666 (2015), no. 1 120002.
  * [18] S. F. King, Unified Models of Neutrinos, Flavour and CP Violation, Prog. Part. Nucl. Phys. 94 (2017) 217–256, [arXiv:1701.04413].
  * [19] S. T. Petcov, Discrete Flavour Symmetries, Neutrino Mixing and Leptonic CP Violation, Eur. Phys. J. C 78 (2018), no. 9 709, [arXiv:1711.10806].
  * [20] F. Feruglio and A. Romanino, Lepton flavor symmetries, Rev. Mod. Phys. 93 (2021), no. 1 015007, [arXiv:1912.06028].
  * [21] T. Kobayashi, H. Ohki, H. Okada, Y. Shimizu, and M. Tanimoto, An Introduction to Non-Abelian Discrete Symmetries for Particle Physicists.  1, 2022.
  * [22] J. C. Criado and F. Feruglio, Modular Invariance Faces Precision Neutrino Data, SciPost Phys. 5 (2018), no. 5 042, [arXiv:1807.01125].
  * [23] T. Kobayashi, N. Omoto, Y. Shimizu, K. Takagi, M. Tanimoto, and T. H. Tatsuishi, Modular A4 invariance and neutrino mixing, JHEP 11 (2018) 196, [arXiv:1808.03012].
  * [24] G.-J. Ding, S. F. King, and X.-G. Liu, Modular A4 symmetry models of neutrinos and charged leptons, JHEP 09 (2019) 074, [arXiv:1907.11714].
  * [25] P. P. Novichkov, J. T. Penedo, S. T. Petcov, and A. V. Titov, Modular S4 models of lepton masses and mixing, JHEP 04 (2019) 005, [arXiv:1811.04933].
  * [26] T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto, and T. H. Tatsuishi, New $A_{4}$ lepton flavor model from $S_{4}$ modular symmetry, JHEP 02 (2020) 097, [arXiv:1907.09141].
  * [27] X. Wang and S. Zhou, The minimal seesaw model with a modular S4 symmetry, JHEP 05 (2020) 017, [arXiv:1910.09473].
  * [28] P. Chen, G.-J. Ding, J.-N. Lu, and J. W. F. Valle, Predictions from warped flavor dynamics based on the $T^{\prime}$ family group, Phys. Rev. D 102 (2020), no. 9 095014, [arXiv:2003.02734].
  * [29] I. de Medeiros Varzielas, S. F. King, and Y.-L. Zhou, Multiple modular symmetries as the origin of flavor, Phys. Rev. D 101 (2020), no. 5 055033, [arXiv:1906.02208].
  * [30] T. Asaka, Y. Heo, T. H. Tatsuishi, and T. Yoshida, Modular $A_{4}$ invariance and leptogenesis, JHEP 01 (2020) 144, [arXiv:1909.06520].
  * [31] T. Asaka, Y. Heo, and T. Yoshida, Lepton flavor model with modular $A_{4}$ symmetry in large volume limit, Phys. Lett. B 811 (2020) 135956, [arXiv:2009.12120].
  * [32] F. J. de Anda, S. F. King, and E. Perdomo, $SU(5)$ grand unified theory with $A_{4}$ modular symmetry, Phys. Rev. D 101 (2020), no. 1 015028, [arXiv:1812.05620].
  * [33] T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto, and T. H. Tatsuishi, Modular $S_{3}$-invariant flavor model in SU(5) grand unified theory, PTEP 2020 (2020), no. 5 053B05, [arXiv:1906.10341].
  * [34] P. P. Novichkov, S. T. Petcov, and M. Tanimoto, Trimaximal Neutrino Mixing from Modular A4 Invariance with Residual Symmetries, Phys. Lett. B 793 (2019) 247–258, [arXiv:1812.11289].
  * [35] T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto, T. H. Tatsuishi, and H. Uchida, Finite modular subgroups for fermion mass matrices and baryon/lepton number violation, Phys. Lett. B 794 (2019) 114–121, [arXiv:1812.11072].
  * [36] H. Okada and M. Tanimoto, CP violation of quarks in $A_{4}$ modular invariance, Phys. Lett. B 791 (2019) 54–61, [arXiv:1812.09677].
  * [37] H. Okada and M. Tanimoto, Towards unification of quark and lepton flavors in $A_{4}$ modular invariance, Eur. Phys. J. C 81 (2021), no. 1 52, [arXiv:1905.13421].
  * [38] T. Nomura and H. Okada, A modular $A_{4}$ symmetric model of dark matter and neutrino, Phys. Lett. B 797 (2019) 134799, [arXiv:1904.03937].
  * [39] H. Okada and Y. Orikasa, Modular $S_{3}$ symmetric radiative seesaw model, Phys. Rev. D 100 (2019), no. 11 115037, [arXiv:1907.04716].
  * [40] T. Nomura and H. Okada, A two loop induced neutrino mass model with modular $A_{4}$ symmetry, Nucl. Phys. B 966 (2021) 115372, [arXiv:1906.03927].
  * [41] T. Nomura, H. Okada, and O. Popov, A modular $A_{4}$ symmetric scotogenic model, Phys. Lett. B 803 (2020) 135294, [arXiv:1908.07457].
  * [42] J. C. Criado, F. Feruglio, and S. J. D. King, Modular Invariant Models of Lepton Masses at Levels 4 and 5, JHEP 02 (2020) 001, [arXiv:1908.11867].
  * [43] S. F. King and Y.-L. Zhou, Trimaximal TM1 mixing with two modular $S_{4}$ groups, Phys. Rev. D 101 (2020), no. 1 015001, [arXiv:1908.02770].
  * [44] G.-J. Ding, S. F. King, X.-G. Liu, and J.-N. Lu, Modular S4 and A4 symmetries and their fixed points: new predictive examples of lepton mixing, JHEP 12 (2019) 030, [arXiv:1910.03460].
  * [45] I. de Medeiros Varzielas, M. Levy, and Y.-L. Zhou, Symmetries and stabilisers in modular invariant flavour models, JHEP 11 (2020) 085, [arXiv:2008.05329].
  * [46] D. Zhang, A modular $A_{4}$ symmetry realization of two-zero textures of the Majorana neutrino mass matrix, Nucl. Phys. B 952 (2020) 114935, [arXiv:1910.07869].
  * [47] T. Nomura, H. Okada, and S. Patra, An inverse seesaw model with $A_{4}$ -modular symmetry, Nucl. Phys. B 967 (2021) 115395, [arXiv:1912.00379].
  * [48] T. Kobayashi, T. Nomura, and T. Shimomura, Type II seesaw models with modular $A_{4}$ symmetry, Phys. Rev. D 102 (2020), no. 3 035019, [arXiv:1912.00637].
  * [49] J.-N. Lu, X.-G. Liu, and G.-J. Ding, Modular symmetry origin of texture zeros and quark lepton unification, Phys. Rev. D 101 (2020), no. 11 115020, [arXiv:1912.07573].
  * [50] X. Wang, Lepton flavor mixing and CP violation in the minimal type-(I+II) seesaw model with a modular $A_{4}$ symmetry, Nucl. Phys. B 957 (2020) 115105, [arXiv:1912.13284].
  * [51] S. J. D. King and S. F. King, Fermion mass hierarchies from modular symmetry, JHEP 09 (2020) 043, [arXiv:2002.00969].
  * [52] M. Abbas, Fermion masses and mixing in modular A4 Symmetry, Phys. Rev. D 103 (2021), no. 5 056016, [arXiv:2002.01929].
  * [53] H. Okada and Y. Shoji, Dirac dark matter in a radiative neutrino model, Phys. Dark Univ. 31 (2021) 100742, [arXiv:2003.11396].
  * [54] H. Okada and Y. Shoji, A radiative seesaw model with three Higgs doublets in modular $A_{4}$ symmetry, Nucl. Phys. B 961 (2020) 115216, [arXiv:2003.13219].
  * [55] G.-J. Ding and F. Feruglio, Testing Moduli and Flavon Dynamics with Neutrino Oscillations, JHEP 06 (2020) 134, [arXiv:2003.13448].
  * [56] H. Okada and M. Tanimoto, Quark and lepton flavors with common modulus $\tau$ in $A_{4}$ modular symmetry, arXiv:2005.00775.
  * [57] H. Okada and M. Tanimoto, Modular invariant flavor model of $A_{4}$ and hierarchical structures at nearby fixed points, Phys. Rev. D 103 (2021), no. 1 015005, [arXiv:2009.14242].
  * [58] K. I. Nagao and H. Okada, Neutrino and dark matter in a gauged $U(1)_{R}$ symmetry, JCAP 05 (2021) 063, [arXiv:2008.13686].
  * [59] X. Wang, B. Yu, and S. Zhou, Double covering of the modular $A_{5}$ group and lepton flavor mixing in the minimal seesaw model, Phys. Rev. D 103 (2021), no. 7 076005, [arXiv:2010.10159].
  * [60] H. Okada and M. Tanimoto, Spontaneous CP violation by modulus $\tau$ in $A_{4}$ model of lepton flavors, JHEP 03 (2021) 010, [arXiv:2012.01688].
  * [61] C.-Y. Yao, J.-N. Lu, and G.-J. Ding, Modular Invariant $A_{4}$ Models for Quarks and Leptons with Generalized CP Symmetry, JHEP 05 (2021) 102, [arXiv:2012.13390].
  * [62] H. Kuranaga, H. Ohki, and S. Uemura, Modular origin of mass hierarchy: Froggatt-Nielsen like mechanism, JHEP 07 (2021) 068, [arXiv:2105.06237].
  * [63] S. T. Petcov and M. Tanimoto, $A_{4}$ Modular Flavour Model of Quark Mass Hierarchies close to the Fixed Point $\tau=\omega$, arXiv:2212.13336.
  * [64] S. Kikuchi, T. Kobayashi, K. Nasu, S. Takada, and H. Uchida, Quark hierarchical structures in modular symmetric flavor models at level 6, arXiv:2301.03737.
  * [65] Y. Abe, T. Higaki, J. Kawamura, and T. Kobayashi, Quark masses and CKM hierarchies from $S_{4}^{\prime}$ modular flavor symmetry, arXiv:2301.07439.
  * [66] S. Kikuchi, T. Kobayashi, K. Nasu, S. Takada, and H. Uchida, Quark mass hierarchies and CP violation in $A_{4}\times A_{4}\times A_{4}$ modular symmetric flavor models, arXiv:2302.03326.
  * [67] F. Feruglio, V. Gherardi, A. Romanino, and A. Titov, Modular invariant dynamics and fermion mass hierarchies around $\tau=i$, JHEP 05 (2021) 242, [arXiv:2101.08718].
  * [68] P. P. Novichkov, J. T. Penedo, and S. T. Petcov, Fermion mass hierarchies, large lepton mixing and residual modular symmetries, JHEP 04 (2021) 206, [arXiv:2102.07488].
  * [69] C. D. Froggatt and H. B. Nielsen, Hierarchy of Quark Masses, Cabibbo Angles and CP Violation, Nucl. Phys. B 147 (1979) 277–298.
  * [70] T. Higaki and J. Kawamura, A low-scale flavon model with a $Z_{N}$ symmetry, JHEP 03 (2020) 129, [arXiv:1911.09127].
  * [71] Particle Data Group Collaboration, R. L. Workman and Others, Review of Particle Physics, PTEP 2022 (2022) 083C01.
  * [72] S. Antusch and V. Maurer, Running quark and lepton parameters at various scales, JHEP 11 (2013) 115, [arXiv:1306.6879].