# Transfer Learning for Molecular Property Predictions from Small Data Sets

Thorren Kirschbaum  [  Annika Bande  [<EMAIL_ADDRESS>
###### Abstract

Machine learning has emerged as a new tool in chemistry to bypass expensive
experiments or quantum-chemical calculations, for example, in high-throughput
screening applications. However, many machine learning studies rely on small
data sets, making it difficult to efficiently implement powerful deep learning
architectures such as message passing neural networks. In this study, we
benchmark common machine learning models for the prediction of molecular
properties on small data sets, for which the best results are obtained with
the message passing neural network PaiNN, as well as SOAP molecular
descriptors concatenated to a set of simple molecular descriptors tailored to
gradient boosting with regression trees. To further improve the predictive
capabilities of PaiNN, we present a transfer learning strategy that uses large
data sets to pre-train the respective models and allows to obtain more
accurate models after fine-tuning on the original data sets. The pre-training
labels are obtained from computationally cheap ab initio or semi-empirical
models and corrected by simple linear regression on the target data set to
obtain labels that are close to those of the original data. This strategy is
tested on the Harvard Oxford Photovoltaics data set (HOPV, HOMO-LUMO-gaps),
for which excellent results are obtained, and on the Freesolv data set
(solvation energies), where this method is unsuccessful due to a complex
underlying learning task and the dissimilar methods used to obtain pre-
training and fine-tuning labels. Finally, we find that the final training
results do not improve monotonically with the size of the pre-training data
set, but pre-training with fewer data points can lead to more biased pre-
trained models and higher accuracy after fine-tuning.

###### keywords:

Machine Learning, Transfer Learning, Message Passing Neural Networks, Semi-
Empirical Methods, Density Functional Theory

Helmholtz-Zentrum Berlin] Theory of Electron Dynamics and Spectroscopy,
Helmholtz-Zentrum Berlin für Materialien und Energie GmbH, Hahn-Meitner-Platz
1, 14109 Berlin, Germany Helmholtz-Zentrum Berlin] Theory of Electron Dynamics
and Spectroscopy, Helmholtz-Zentrum Berlin für Materialien und Energie GmbH,
Hahn-Meitner-Platz 1, 14109 Berlin, Germany [LU Hannover]Institute of
Inorganic Chemistry, Leibniz University Hannover, Callinstr. 9, 30167
Hannover, Germany

##  1 Introduction

In the realm of molecular sciences, predicting the properties of molecules
stands as a pivotal pursuit across numerous domains, including
pharmaceuticals, materials science, and environmental studies. Traditionally,
this endeavor has been tackled by quantum chemistry, relying on sophisticated
theoretical frameworks and computational simulations to unravel the intricate
behaviors of molecules. However, the advent of machine learning (ML) has
ushered in a paradigm shift, offering an alternative approach that leverages
data-driven methodologies efficiently compute molecular properties.1, 2, 3, 4,
5

One notable application of ML in this domain is its integration into high
throughput screening workflows, where it replaces conventional quantum
chemistry methods.3 Through the judicious application of ML algorithms,
researchers can expedite the process of identifying promising molecules for
drug discovery or materials development, thereby accelerating the pace of
scientific innovation. However, the powerful deep neural networks (NNs) often
used in such approaches usually need large amounts of training data to finally
make accurate predictions.

Machine learning algorithms, particularly deep NNs, have exhibited remarkable
success in discerning intricate patterns and relationships within data sets of
substantial size. In recent developments, the message passing neural network
PaiNN has demonstrated highest performance in predicting molecular properties,
achieving state-of-the-art results across various datasets.6, 7 Another
prevalent approach in ML for molecular sciences entails the utilization of
tailored molecular descriptors, such as SOAP,8, 9 in conjunction with
traditional ML models like kernel ridge regression (KRR). In general, as a
rule of thumb, deep learning architectures are to be preferred when $>$
103–104 training points are available, and different approaches, such as KRR,
should be considered otherwise.10 Yet, in the context of molecular property
prediction, data sets often fall short regarding the typically required large
amounts of data for deep NNs. This mismatch often presents a formidable
obstacle, impeding the efficacy of ML in scenarios characterized by small data
sets, which are common in practical applications.11, 12, 13 Bridging this
divide necessitates innovative strategies and methodologies tailored to the
unique challenges posed by small data sets. In this paper, we present a
transfer learning strategy to address this pressing need.

In transfer learning, the knowledge learned in one setting is used to improve
performance on a related target task.14 In a first step, known as pre-
training, an NN is trained on a large set of data that is closely related to
the target task, yielding an NN with pre-optimized parameters. In the second
step, termed fine-tuning, the same model is re-trained using data of the
target task. Now, since the model is already close to the optimal set of
parameters after the pre-training, much less data is needed during fine-tuning
to obtain an accurate model for the target task. This strategy has proven
useful in various fields, including image classification15 and captioning,16
gaming strategies,17 and social network analysis.18 In the physical sciences,
some previous studies have used transfer learning to improve property
predictions of molecules and materials.19, 20 Fur such tasks, Sun et al. found
that supervised pre-training generally helps, if the pre-training labels are
closely aligned with the downstream tasks, but the success of efficient
transfer learning strongly depends on the employed data sets, ML models, and
their hyperparameters.21 Moreover, a range of studies recently investigated
the use of transfer learning to improve atomic force predictions for molecular
dynamics simulations.22, 23, 24, 25, 26, 27, 28, 29

In this study, we start by undertaking a rigorous benchmarking, evaluating
standard ML techniques in the context of small molecular data sets. In
particular, we compare the use of SOAP descriptors tailored to a standard ML
algorithm (KRR, NN, or gradient boosting) to the message passing neural
network PaiNN for regression. The comparison is done both on the Harvard
Organic Photovoltaics (HOPV) data set,30 containing photovoltaic properties of
350 organic molecules, and on the Freesolv data set,31 which contains
solvation energies of 643 small organic molecules.

Secondly, we propose a transfer learning strategy as a means of augmenting
predictive modeling with small data sets. It involves leveraging computations
from computationally cheap ab initio or semi-empirical models coupled with
simple linear regression, to generate pre-training labels for PaiNN.

##  2 Methods

In this study, we use two small common physical chemistry data sets, the
Harvard Organic Photovoltaics (HOPV) data set30 and the Freesolv data set.31
The HOPV data set contains 350 molecules and their photovoltaic data obtained
from both experiments and quantum-chemical calculations. In this study, we
focus on the HOMO-LUMO gaps obtained from PBE0/def2-SVP density functional
theory, which is a common ab initio method used for standard chemistry data
sets.32, 33, 34, 35, 36, 37, 7 Note that the version of the HOPV data set used
by us (HOPV_15_revised_2) contains a single molecule with an Se atom, which we
regarded as a structural outlier and therefore removed from the data set prior
to training. The Freesolv data set contains experimental and calculated
solvation energies of 643 small molecules,31 from which we use the values
obtained from ”alchemical” calculations based on molecular dynamics
simulations by the Mobley group.38

For benchmarking standard ML methods with HOPV and Freesolv, we employed the
PaiNN implementation of the schnetpack Python library,39 as well as the SOAP
implementation of the dscribe library40 in junction with either sklearn’s
GradientBoostingRegressor (GBoost), kernel ridge regression, or a simple feed-
forward NN. The GBoost algorithm profited from an additional set of simple
molecular descriptors concatenated to the SOAP descriptor. For NN learning, to
limit the size of the NNs, the sparse SOAP descriptors were compressed by
principal component analysis, retaining 99.999 % of the variance.9, 7 The
hyperparameters of PaiNN were optimized using an iterative grid search, and
the hyperparameters of the SOAP-ML models were optimized using Bayesian
optimization. For PaiNN, a learning rate decay and early stopping were
employed, and results were averaged over five runs. The final model
hyperparameters are listed in the appendix.

In the second part of this study, we present a method to improve the
predictions of the message passing neural network PaiNN,6 which was also
identified as the best ML model for the data sets at hand (see Results and
Discussion), via transfer learning. For this, we need (i) a (small) molecular
data set with high-accuracy data of the target property ($y_{\textrm{true}}$),
(ii) a large data set of molecular structures to be used for pre-training, and
(iii) a computationally inexpensive method from the field of quantum chemistry
to compute the target property for the pre-training data. Cheap ab initio
methods, like density functional theory in the linear density approximation
(LDA-DFT), or semi-empirical methods typically yield results which are
qualitatively correct, but quantitatively far off.41 To correct the results
obtained from these methods, we use the respective quantum-chemical (QC)
method to compute the target property on the original high-quality data set
(i), $y_{\textrm{QC}}$, and apply simple linear regression to obtain the
coefficients $a$ and $b$ that optimize the equation $(a\cdot
y_{\textrm{QC}}+b)-y_{\textrm{true}}=0$ for this data set. Subsequently, we
use the same QC method to compute the target property for the structures of
the pre-training data set (ii), again $y_{\textrm{QC}}$, and finally obtain
the pre-training properties as $y_{\textrm{pt}}=a\cdot y_{\textrm{QC}}+b$,
thought to be similar to that of the methods used in (i).

In this study, we use the molecular structures of the OE62 data set as pre-
training structures for the HOPV data set, and the molecular structures of the
QM9 data set as pre-training structures for Freesolv (ii). The OE62 data set
contains crystal structures of 62 k medium-sized organic molecules,37 and the
QM9 data set contains DFT-optimized structures of 134 k small molecules.42 To
obtain pre-training labels for these data sets (iii), HOMO-LUMO gaps and
solvation energies were calculated with LDA-DFT using the ORCA software,43 and
tight-binding density functional theory (DFTB) using the GFN2-xTB method44
implemented in the XTB software.45 Solvation energies were obtained as the
difference in single point energy of the molecule in vacuum and the molecule
with an implicit solvent model for water (conductor-like polarizable continuum
model, CPCM in ORCA,46 and generalized Born model with surface area
contributions, GBSA in XTB,47, 48, 49 respectively). In Figure 1, we display
the scatter plots of the original HOPV and Freesolv data, respectively,
against the data obtained from LDA and XTB, respectively, corrected by linear
regression. These plots clearly show that the obtained HOMO-LUMO gaps are very
close to the original HOPV data, but the solvation energies are quite
dissimilar to the original Freesolv data. Furthermore, for both data sets, the
properties obtained from LDA are closer to the original data those obtained
from XTB. The mean absolute errors of the fits are: 0.05 eV (HOPV LDA), 0.12
eV (HOPV XTB), 1.1 kcal/mol (Freesolv LDA), 1.7 kcal/mol (Freesolv XTB).

Figure 1: Scatter plots of the original data plotted against data obtained
from XTB (left) or LDA (right) plus linear regression for the data sets HOPV
(HOMO-LUMO-gaps, top) and Freesolv (solvation energies, bottom).

Furthermore, we also GBoost trained on HOPV and Freesolv, respectively, to
compute HOMO-LUMO gaps for OE62 structures and solvation energies for QM9
structures, respectively, and used these for pre-training PaiNN. The results
are then compared to the pre-training described above, in which the pre-
training labels are obtained from cheap quantum-chemical methods.

On average, HOMO-LUMO gap computations took approximately 106 s (LDA-DFT),
0.22 s (XTB), and $<$ 0.1 s (GBoost) per data point in OE62, respectively.
Solvation energy computations took approximately 70 s (LDA-DFT), 0.15 s (XTB),
and $<$ 0.1 s (GBoost) per data point in QM9, respectively, and 1.33 s (XTB)
per data point in OE62.

Pre-training on QM9 data was done for 200 epochs (100 k training, 10 k
validation, 20 k test data points). Pre-training on OE62 data was done for 250
epochs for LDA and XTB data (50 k training, 5 k validation, 6 k test data
points), and for 300 epochs for GBoost data (35 k training, 5 k validation,
5537 testing data points). For pre-training with GBoost data, the pre-training
data containing elements which did not occur in the fine-tuning data set were
omitted to keep the employed SOAP descriptors reasonably small and therefore
minimize memory requirements. Finally, pre-training was done using 35 k
training data points (SOAP+SD+GBoost/OE62 for HOPV, and SOAP+SD+GBoost/QM9 for
Freesolv) or 30 k data points (SOAP+SD+GBoost/OE62 for Freesolv),
respectively. For each method, three models were trained with different random
seeds and the best model as evaluated on the test subset was selected for
fine-tuning.

##  3 Results and Discussion

###  3.1 ML benchmarks

We start by benchmarking standard chemistry ML models for regression on the
two small, but high-quality data sets under investigation, HOPV30 and
Freesolv,31 listed under (i) in the Methods section. The results are
summarized in Table 1 as mean absolute errors (MAEs, root mean squared errors
(RMSE) are given in parenthesis). For both data sets, very low errors are
obtained with PaiNN, yielding MAEs of 0.20 eV for HOPV and 0.56 kcal/mol for
Freesolv, respectively. For HOPV, the SOAP+SD+GBoost model yields even
slightly better results with an MAE of 0.19 eV, while 0.86 kcal/mol are
obtained for regression on Freesolv. The results obtained with these
algorithms are qualitatively similar to those obtained in previous ML studies
on HOPV50, 51, 52, 53, 54, 55, 56 and Freesolv.57, 58, 59, 60, 61, 62, 63, 64,
56, 65, 66, 67

Table 1: Mean absolute errors for HOPV HOMO-LUMO gaps (in eV) and Freesolv solvation energies (in kcal/mol) from four different models. The first three models use SOAP molecular descriptors connected to either KRR, a dense NN, or GBoost. The GBoost model uses an additional set of simple molecular descriptors (SD) concatenated with the SOAP descriptors. The final model is the message passing NN PaiNN. Uncertainties are measured in RMSE (in parentheses). Best MAE results are bold. Model |  | HOPV
---
MAE (RMSE) / eV
| Freesolv
---
MAE (RMSE) / kcal/mol
SOAP+KRR | 0.31 (0.46) | 1.7 (2.2)
SOAP+NN | 0.25 (0.33) | 2.1 (3.2)
SOAP+SD+GBoost | 0.19 (0.30) | 0.86 (1.30)
PaiNN | 0.20 (0.31) | 0.56 (0.86)

The remaining models SOAP+NN and SOAP+KRR perform considerably worse on both
data sets, in direct comparison SOAP+NN yielding better results for HOPV while
SOAP+KRR yields better results for Freesolv. We also note that the MAEs from
the different models are much closer for regression on HOPV, where the worst
model (SOAP+KRR) has a 63 % higher MAE than the best model (SOAP+SD+GBoost).
In contrast, for Freesolv, the worst model (SOAP+NN) has a 375 % higher MAE
than the best model (PaiNN), indicating that the learning task is much more
complex than for HOPV (cf. also Figure 1. For the following transfer learning
with PaiNN, the best remaining model was selected to generate the respective
pre-training labels, which is SOAP+SD+GBoost for both data sets.

###  3.2 Transfer Learning

We use both cheap models from quantum chemistry (LDA-DFT and XTB-DFTB)
corrected by simple linear regression, and an ML regressor (SOAP+SD+GBoost,
from here on named GBoost) to obtain labels for pre-training the PaiNN models
(details in the Methods section). For HOPV, consisting of medium-sized organic
molecules, we use the structures of the OE62 data set for pre-training, which
contains ca. 62 k diverse structures of medium-sized organic molecules. For
Freesolv, which on the other hand consists of small organic molecules, we use
the structures of the QM9 data set for pre-training, which contains ca. 132 k
organic molecules with up to nine heavy (non-H) atoms. Tables 2 and 3
summarize the results obtained for regression on the HOPV and Freesolv data
set, respectively, after PaiNN training from scratch (cf. Table 1), or after
pre-training on data labeled with either GBoost, XTB, or LDA, respectively.
The learning curves displaying the MAE vs. the number of fine-tuning training
points are shown in Figure 2 (here, pre-training is done on the full
respective pre-training data set), and the learning curves displaying the MAE
vs. the number of pre-training points are shown in Figure 3 (here, fine-tuning
is done on the full respective fine-tuning data set), each for the HOPV set on
top and the Freesolv set in the bottom.

Table 2: MAEs for HOPV HOMO-LUMO gaps (in eV) from PaiNN obtained after training from scratch, or after pre-training on OE62 data with labels obtained from GBoost, XTB or LDA, respectively. Uncertainties are measured in RMSE (in parentheses). Best MAE results are bold. Pre-training |  | HOPV
---
MAE (RMSE) / eV
None | 0.20 (0.31)
OE62+GBoost | 0.18 (0.25)
OE62+XTB | 0.13 (0.17)
OE62+LDA | 0.13 (0.19)
Table 3: MAEs for Freesolv solvation energies (in kcal/mol) from PaiNN obtained after training from scratch, or after pre-training on OE62 data with labels obtained from GBoost, XTB or LDA, respectively. Uncertainties are measured in RMSE (in parentheses). Best MAE results are bold. Pre-training |  | Freesolv
---
MAE (RMSE) / kcal/mol
None | 0.56 (0.86)
QM9+GBoost | 0.74 (1.19)
QM9+XTB | 0.60 (0.99)
QM9+LDA | 0.64 (1.00)
Figure 2: Fine-tuning learning curves for PaiNN training on HOPV (top) and
Freesolv (bottom), after training from scratch (black) or after pre-training
on OE62 or QM9 data, respectively, with labels obtained from GBoost (red), XTB
(light blue) or LDA-DFT (dark blue). The MAE is plotted against the number of
training examples used for fine-tuning (log scale). Figure 3: Pre-training
learning curves for PaiNN training on HOPV (top) and Freesolv (bottom), after
training from scratch (black star) or after pre-training on OE62 or QM9 data,
respectively, with labels obtained from GBoost (red), XTB (light blue) or LDA-
DFT (dark blue). The mean absolute error (MAE) is plotted against the number
of training examples used for pre-training.

For HOPV (Table 2), pre-training results in an improved learning performance.
The MAE is reduced from 0.20 eV (training from scratch) to 0.18 eV after pre-
training on OE62+GBoost, and further to 0.13 eV after pre-training on OE62+XTB
or OE62+LDA. The pre-training with GBoost is less effective, probably due to
the fact that for this method, some training points were removed from OE62 to
avoid creating overly large SOAP descriptors for the ML model. However, the
GBoost model might also produce less accurate labels than the QC models
corrected by linear regression, since the GBoost MAE (0.19 eV, Table 1) is
higher than the MAEs of the LDA or XTB plus linear regression fit on HOPV
(MAEs of 0.05 eV for LDA, and 0.12 eV for XTB, respectively).

The fine-tuning learning curves for HOPV (Figure 2, top) show that this
improvement of ca. 0.07 eV over training from scratch (black curve) after pre-
training on OE62+XTB (light blue) or OE62+LDA (dark blue) is relatively
stable, independent of the number of HOPV fine-tuning training points. Only,
after pre-training on OE62+GBoost, the improvement obtained from pre-training
is much higher when fine-tuning on increasingly many HOPV data points (red
curve).

The HOPV pre-training learning curves (Figure 3, top) confirm that for all
pre-training data set sizes tested, the labels obtained from the QC models
lead to a higher accuracy as compared to the ML-generated labels.
Interestingly, for GBoost and LDA labels, increasing the size of the pre-
training data set does not monotonically improve the final training results.
Instead, there appears to be an optimal training set size at around 20 k pre-
training points for GBoost-labeled data and 30 k pre-training points for LDA-
labeled data, respectively, after which further improving the pre-training
data set size results in higher MAEs after fine-tuning. For XTB-labeled data,
a similar optimal pre-training data set size may exist outside the range of
the data available here. This indicates that a pre-trained model that is too
strongly fitted on the pre-training task may be worse for transfer learning
than a more biased model that only learned to grasp high-level concepts from a
limited amount of training data.

Regression learning on Freesolv, on the other hand, is not improved by pre-
training the PaiNN models (Table 3). Instead, the final training MAEs are
increased from 0.56 kcal/mol (training from scratch) to 0.74 kcal/mol
(GBoost), 0.60 kcal/mol (XTB), and 0.64 kcal/mol (LDA). The fine-tuning
learning curves (Figure 2, bottom) show that an improvement is achieved when
fine-tuning is done with very few ($\leq$ ca. 64) data points, QM9+LDA fine-
tuning resulting in the largest improvement. However, for larger fine-tuning
data set sizes, the results obtained after QM9+XTB or QM9+LDA pre-training are
slightly worse as compared to training from scratch, and pre-training with
QM9+GBoost results in the worst performance. Furthermore, the pre-training
learning curves (Figure 3, bottom) indicate that the loss of performance is
approximately independent from the size of the pre-training data set.
Moreover, for pre-training data obtained from GBoost and LDA, the final models
tend to even decrease in quality after being pre-trained on larger data sets.

As discussed before, the regression task on Freesolv might be more complex as
compared to regression on HOPV. As a consequence, the pre-training labels
generated by GBoost are less reliable for Freesolv solvation energies than
those generated for HOPV HOMO-LUMO gaps (cf. Table 1). Additionally, the
methods used for producing the respective pre-training labels (QC methods) are
quite similar to the method used for generating HOPV data (standard DFT),30
but quite dissimilar to the method used for generating Freesolv data
(”alchemical” calculations based on molecular dynamics simulations, see also
the scatter plots in Figure 1).38 These factors may be the reason for why the
employed transfer learning strategy does work very well for use with the HOPV
data set, but not for Freesolv.

Another hypothesis for why transfer learning might fail in certain scenarios
is that the underlying regression task is too easy, and therefore no
additional pre-training is needed.21 However, previous results for regression
on molecular energies on the QM9 data set (free energy $G$, internal energy
$U$, enthalpy $H$) using PaiNN indicated that MAEs as low as ca. 0.14 kcal/mol
can be obtained on a sufficiently large data set, which is considerably lower
than the MAEs obtained here.6 Also, these quantities ($G$, $U$, $H$) can
easily be calculated from only the respective molecular structure, but
solvation energy computations additionally need to take into account the shape
and intermolecular interactions of the respective solvation shell, thus adding
another level of complexity. Therefore, we anticipate that the regression task
on Freesolv is not too easy for transfer learning.

In an attempt to obtain pre-trained models that are more robust and therefore
allegedly more suitable for transfer learning, we added additional data points
to the respective pre-training data sets. For this, we used XTB to generate
additional labels for QM9 (HOMO-LUMO gaps) and OE62 (solvation energies).
However, we found that pre-training on a joint OE62 and QM9 data set did not
improve the final fine-tuning results as compared to pre-training on OE62 or
QM9 data only, respectively, for both HOPV and Freesolv training. This is
surprising, because NNs usually get more accurate when increasing the amount
of training data, and new dissimilar data points can help to create more
robust models. Here, instead, the introduction of more data points that are
dissimilar to the data contained in the final fine-tuning data set apparently
impedes the learning process. In particular, for Freesolv, the underlying
problem of having a complex learning task and dissimilar methods for
generating pre-training and fine-tuning labels could not be solved by
augmenting the pre-training data set. This observation is in line with
previous work on transfer learning in chemistry, which stated that the pre-
training and fine-tuning data sets should be as closely aligned as possible to
obtain the best results.19

Following this idea, we aimed to instead make the pre-training and fine-tuning
data sets more similar to each other: The HOPV data set consists almost
entirely of S-functionalized aromatic molecules, while the OE62 data set used
for pre-training contains a much more diverse set of molecules, including more
structural variety and a wider range of atomic species. Accordingly, we
filtered all molecules from the OE62 pre-training data set containing any of
the elements As, Se, Br, Te or I, and furthermore removed all molecules that
do not contain any sulfur atoms, thereby reducing the data set size from
61,489 to 41,487 molecules. Pre-training with 20 k less molecules from the
filtered OE62 data set resulted in a slight loss of performance compared to
pre-training of the full OE62+XTB data set (+0.01 eV MAE), which is
approximately the same result as obtained after non-selective filtering of the
pre-training data set (+0.02 eV MAE, cf. Figure 3). This indicates that
increased similarity between pre-training and fine-tuning structures might
slightly improve the learning performance by yielding more tailored pre-
trained models, but the effect found here is very small.

Finally, previous works have used discriminative fine-tuning to further
optimize transfer learning.29 In this method, the learning rate during fine-
tuning is adapted layer-wise, such that the first layers have a much smaller
learning rate than the final layers of the model. The idea is that the first
layers, which learn to grasp fundamental features of the molecular data, will
be sufficiently trained during pre-training, while the last layers, which
derive the final predictions, will need more fine-tuning. We implemented
discriminative fine-tuning for PaiNN, following the approach outlined by
Zaverkin et al.29 Here, the learning rate (lr) for layer n$-$1 was adjusted as
lrn-1 = lrn/5 with lr = 5e-4 being the learning rate of the last layer.
However, we found that this method did not yield any improvements over non-
discriminative fine-tuning for transfer learning with PaiNN. In contrast, we
found that in general, modifications of the learning rate and learning rate
decay had a significant impact on the learning performance, and sticking to
the parameters identified by the initial hyperparameter sweep continued to
yield the best results.

##  4 Conclusions

This study investigated the use of common ML algorithms for use with small
data sets, focusing on HOPV (350 data points) and Freesolv (643 data points).
We found that both the MPNN PaiNN, and gradient boosting with regression trees
operating on SOAP molecular descriptors concatenated to a set of simple
molecular descriptors yield accurate results, clearly outperforming the
SOAP+NN and SOAP+KRR architectures.

Furthermore, we presented a method for data-efficient ML on small physical
chemistry data sets using transfer learning. The pre-training data was
obtained from computationally cheapQC calculations corrected by simple linear
regression calibrated on the fine-tuning data set. This method utilizes the
fact that many cheap ab initio or semi-empirical techniques yield
quantitatively incorrect but qualitatively correct data, which can easily be
corrected via linear regression. For the HOPV data set, the introduced method
outperforms transfer learning with labels obtained from an ML model. Here,
pre-training with labels computed by either XTB or LDA result in very similar
learning performances, but XTB computations were much faster than LDA
computations. Accordingly, labelling the pre-training data with XTB should be
preferred.

However, the same method did not yield improved learning results for
regression on the Freesolv data set. This was due to the fact that the
underlying learning task is much more complex and the data used for pre-
training is too dissimilar to the data contained in Freesolv.

Furthermore, the ML results indicate that pre-training with a limited amount
of data points may be beneficial, since the obtained models are slightly more
biased and less strongly fitted to the pre-training task. Furthermore, we
found that augmenting the pre-training data sets with additional structures
that are less similar to the fine-tuning data does not improve the learning
performance. Finally, for the models trained here, sticking to the
hyperparameters identified by the initial hyperparameter sweeps was more
efficient than adapting the learning rates as done in discriminative fine-
tuning.

In conclusion, for chemistry ML tasks that use small data sets, an effective
transfer learning strategy should use pre-training on molecules that are
similar to the molecules of the downstream task, a labelling method that is
closely aligned to the method used for labelling the fine-tuning data, and an
optimized training set size for the pre-training.

##  5 Appendix

###  5.1 Machine Learning Model Hyperparameters

We here provide details of the ML models used for regression on both HOPV and
Freesolv (cf. Table 1).

The final SOAP+KRR models used the following hyperparameters:
HOPV: $r_{cut}^{SOAP}=5.63$, $n_{max}^{SOAP}=7$, $l_{max}^{SOAP}=7$,
$alpha_{KRR}=1.12$, kernel: cosine.
Freesolv: $r_{cut}^{SOAP}=9.53$, $n_{max}^{SOAP}=8$, $l_{max}^{SOAP}=8$,
$alpha_{KRR}=0.59$, kernel: = cosine.

In the SOAP+NN model, the neural network consists of densely connected layers
with steadily decreasing numbers of neurons, including batch normalization and
a dropout layer with $p=0.3$ after the first layer. Here, 20 % of the training
data was used for validation, and early stopping was employed. The final
models used the following hyperparameters:
HOPV: $r_{cut}^{SOAP}=7.22$, $n_{max}^{SOAP}=8$, $l_{max}^{SOAP}=8$,
$n_{layers}=2$, $lr=0.00860$, nonlinearity: swish.
Freesolv: $r_{cut}^{SOAP}=9.94$, $n_{max}^{SOAP}=7$, $l_{max}^{SOAP}=7$,
$n_{layers}=4$, $lr=0.00627$, nonlin.: swish.

In the SOAP+SD+GBoost model, the following simple molecular descriptors were
employed: total number of atoms, number of atoms for each atom type in the
data set, one-hot-encoding for each atom type (present/not present in the
respective molecule), mean and standard deviation of the carbon-carbon-
distance of neighboring carbon atoms within the molecule. Using these
descriptors with the KRR or NN did not yield any improvement for the
respective learning algorithm. The final SOAP+SD+GBoost models used the
following hyperparameters:
HOPV: $r_{cut}^{SOAP}=9.71$, $n_{max}^{SOAP}=5$, $l_{max}^{SOAP}=8$,
$lr_{GBoost}=0.059$, $n\\_estimators_{GBoost}=295$, $max\\_depth_{GBoost}=4$.
Freesolv: $r_{cut}^{SOAP}=6.85$, $n_{max}^{SOAP}=4$, $l_{max}^{SOAP}=4$,
$lr_{GBoost}=0.250$, $n\\_estimators_{GBoost}=158$, $max\\_depth_{GBoost}=4$.

The final PaiNN models used the following hyperparameters:
HOPV: $radial\\_basis=GaussianRBF$, $r_{cut}=5.0$, $n\\_rbf=20$,
$n\\_atom\\_basis=30$, $n\\_interactions=3$, $cutoff\\_fn=CosineCutoff$,
$lr=5e-4$, $lr\\_decay\\_factor=0.5$, $lr\\_decay\\_patience=5$,
$batchsize=10$.
Freesolv: $radial\\_basis=GaussianRBF$, $r_{cut}=8.0$, $n\\_rbf=30$,
$n\\_atom\\_basis=200$, $n\\_interactions=3$, $cutoff\\_fn=CosineCutoff$,
$lr=5e-4$, $lr\\_decay\\_factor=0.5$, $lr\\_decay\\_patience=5$,
$batchsize=10$.

TK and AB acknowledge support from the Helmholtz Einstein International Berlin
Research School in Data Science (HEIBRiDS). Computing resources were kindly
provided by the Freie Universität Berlin hpc cluster Curta68 and by the
Helmholtz-Zentrum Dresden-Rossendorf. The introduction was in part formulated
with the help of ChatGPT 3.5.

## References

  * Montavon et al. 2013 Montavon, G.; Rupp, M.; Gobre, V.; Vazquez-Mayagoitia, A.; Hansen, K.; Tkatchenko, A.; Müller, K.-R.; Anatole von Lilienfeld, O. Machine learning of molecular electronic properties in chemical compound space. _New J. Phys._ 2013, _15_ , 095003
  * Hansen et al. 2015 Hansen, K.; Biegler, F.; Ramakrishnan, R.; Pronobis, W.; von Lilienfeld, O. A.; Müller, K.-R.; Tkatchenko, A. Machine Learning Predictions of Molecular Properties: Accurate Many-Body Potentials and Nonlocality in Chemical Space. _J. Phys. Chem. Lett._ 2015, _6_ , 2326–2331
  * Dral 2020 Dral, P. O. Quantum Chemistry in the Age of Machine Learning. _J. Phys. Chem. Lett._ 2020, _11_ , 2336–2347
  * Fung et al. 2021 Fung, V.; Zhang, J.; Juarez, E.; Sumpter, B. G. Benchmarking graph neural networks for materials chemistry. _npj Comput. Mater._ 2021, _7_ , 1–8
  * Singh et al. 2022 Singh, K.; Münchmeyer, J.; Weber, L.; Leser, U.; Bande, A. Graph Neural Networks for Learning Molecular Excitation Spectra. _J. Chem. Theory Comput._ 2022, _18_ , 4408–4417
  * Schütt et al. 2021 Schütt, K.; Unke, O.; Gastegger, M. Equivariant message passing for the prediction of tensorial properties and molecular spectra. Proceedings of the 38th International Conference on Machine Learning. 2021; pp 9377–9388
  * Kirschbaum et al. 2023 Kirschbaum, T.; von Seggern, B.; Dzubiella, J.; Bande, A.; Noé, F. Machine Learning Frontier Orbital Energies of Nanodiamonds. _J. Chem. Theory Comput._ 2023, _19_ , 4461–4473
  * Bartók et al. 2013 Bartók, A. P.; Kondor, R.; Csányi, G. On representing chemical environments. _Phys. Rev. B_ 2013, _87_ , 184115
  * Parsaeifard et al. 2021 Parsaeifard, B.; De, D. S.; Christensen, A. S.; Faber, F. A.; Kocer, E.; De, S.; Behler, J.; Lilienfeld, O. A. v.; Goedecker, S. An assessment of the structural resolution of various fingerprints commonly used in machine learning. _Mach. Learn.: Sci. Technol._ 2021, _2_ , 015018
  * Unke et al. 2021 Unke, O. T.; Chmiela, S.; Sauceda, H. E.; Gastegger, M.; Poltavsky, I.; Schütt, K. T.; Tkatchenko, A.; Müller, K.-R. Machine Learning Force Fields. _Chem. Rev._ 2021, _121_ , 10142–10186
  * Ramakrishnan et al. 2015 Ramakrishnan, R.; Dral, P. O.; Rupp, M.; von Lilienfeld, O. A. Big Data Meets Quantum Chemistry Approximations: The ?-Machine Learning Approach. _J. Chem. Theory Comput._ 2015, _11_ , 2087–2096
  * Zhang and Ling 2018 Zhang, Y.; Ling, C. A strategy to apply machine learning to small datasets in materials science. _npj Comput. Mater._ 2018, _4_ , 1–8
  * Mazouin et al. 2022 Mazouin, B.; Schöpfer, A. A.; Lilienfeld, O. A. v. Selected machine learning of HOMO–LUMO gaps with improved data-efficiency. _Mater. Adv._ 2022, 8306–8316
  * Zhuang et al. 2020 Zhuang, F.; Qi, Z.; Duan, K.; Xi, D.; Zhu, Y.; Zhu, H.; Xiong, H.; He, Q. A comprehensive survey on transfer learning. _Proceedings of the IEEE_ 2020, _109_ , 43–76
  * Saenko et al. 2010 Saenko, K.; Kulis, B.; Fritz, M.; Darrell, T. Computer Vision-ECCV 2010, edited by K. Daniilidis, P. Maragos & N. Paragios, Vol. 6316 of Lecture Notes in Computer Science. 2010
  * Vinyals et al. 2015 Vinyals, O.; Toshev, A.; Bengio, S.; Erhan, D. Show and tell: A neural image caption generator. Proceedings of the IEEE conference on computer vision and pattern recognition. 2015; pp 3156–3164
  * Sharma et al. 2007 Sharma, M.; Holmes, M.; Santamaria, J.; Irani, A.; Isbell, J. CL, & Ram, A.(2007). Transfer learning in real-time strategy games using hybrid CBR/RL. Proceedings of the Twentieth International Joint Conference on Artificial Intelligence. 2007; pp 1041–1046
  * Tang et al. 2016 Tang, J.; Lou, T.; Kleinberg, J.; Wu, S. Transfer learning to infer social ties across heterogeneous networks. _ACM Trans. Inf. Syst._ 2016, _34_ , 1–43
  * Li et al. 2022 Li, H.; Zhao, X.; Li, S.; Wan, F.; Zhao, D.; Zeng, J. Improving molecular property prediction through a task similarity enhanced transfer learning strategy. _Iscience_ 2022, _25_ , 105231
  * Shui et al. 2022 Shui, Z.; Karls, D.; Wen, M.; Tadmor, E.; Karypis, G.; others Injecting domain knowledge from empirical interatomic potentials to neural networks for predicting material properties. _Adv. Neural Inf. Process. Syst._ 2022, _35_ , 14839–14851
  * Sun et al. 2022 Sun, R.; Dai, H.; Yu, A. W. Does GNN Pretraining Help Molecular Representation? _Adv. Neural Inf. Process. Syst._ 2022, _35_ , 12096–12109
  * Smith et al. 2019 Smith, J. S.; Nebgen, B. T.; Zubatyuk, R.; Lubbers, N.; Devereux, C.; Barros, K.; Tretiak, S.; Isayev, O.; Roitberg, A. E. Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning. _Nat. Commun._ 2019, _10_ , 2903
  * Käser et al. 2020 Käser, S.; Koner, D.; Christensen, A. S.; von Lilienfeld, O. A.; Meuwly, M. _J. Phys. Chem. A_ 2020, _124_ , 8853–8865
  * Zheng et al. 2021 Zheng, P.; Zubatyuk, R.; Wu, W.; Isayev, O.; Dral, P. O. Artificial intelligence-enhanced quantum chemical method with broad applicability. _Nat. Commun._ 2021, _12_ , 7022
  * Käser et al. 2022 Käser, S.; Richardson, J. O.; Meuwly, M. Transfer learning for affordable and high-quality tunneling splittings from instanton calculations. _J. Chem. Theory Comput._ 2022, _18_ , 6840–6850
  * Gao et al. 2022 Gao, X.; Gao, W.; Xiao, W.; Wang, Z.; Wang, C.; Xiang, L. Supervised Pretraining for Molecular Force Fields and Properties Prediction. _arXiv_ 2022, arXiv:2211.14429
  * Zhang et al. 2022 Zhang, D.; Bi, H.; Dai, F.-Z.; Jiang, W.; Zhang, L.; Wang, H. DPA-1: Pretraining of Attention-based Deep Potential Model for Molecular Simulation. _arXiv_ 2022, arXiv:2208.08236
  * Chen et al. 2023 Chen, M. S.; Lee, J.; Ye, H.-Z.; Berkelbach, T. C.; Reichman, D. R.; Markland, T. E. Data-efficient machine learning potentials from transfer learning of periodic correlated electronic structure methods: Liquid water at AFQMC, CCSD, and CCSD (T) accuracy. _J. Chem. Theory Comput._ 2023, _19_ , 4510–4519
  * Zaverkin et al. 2023 Zaverkin, V.; Holzmüller, D.; Bonfirraro, L.; Kästner, J. Transfer learning for chemically accurate interatomic neural network potentials. _Phys. Chem. Chem. Phys._ 2023, _25_ , 5383–5396
  * Lopez et al. 2016 Lopez, S. A.; Pyzer-Knapp, E. O.; Simm, G. N.; Lutzow, T.; Li, K.; Seress, L. R.; Hachmann, J.; Aspuru-Guzik, A. The Harvard organic photovoltaic dataset. _Sci. Data_ 2016, _3_ , 160086
  * Duarte Ramos Matos et al. 2017 Duarte Ramos Matos, G.; Kyu, D. Y.; Loeffler, H. H.; Chodera, J. D.; Shirts, M. R.; Mobley, D. L. Approaches for Calculating Solvation Free Energies and Enthalpies Demonstrated with an Update of the FreeSolv Database. _J. Chem. Eng. Data_ 2017, _62_ , 1559–1569
  * Adamo and Barone 1999 Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. _J. Chem. Phys._ 1999, _110_ , 6158–6170
  * Perdew et al. 1996 Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. _J. Chem. Phys._ 1996, _105_ , 9982–9985
  * Schäfer et al. 1992 Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. _J. Chem. Phys._ 1992, _97_ , 2571–2577
  * Weigend and Ahlrichs 2005 Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. _Phys. Chem. Chem. Phys._ 2005, _7_ , 3297–3305
  * Rupp et al. 2012 Rupp, M.; Tkatchenko, A.; Müller, K.-R.; von Lilienfeld, O. A. Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning. _Phys. Rev. Lett._ 2012, _108_ , 058301
  * Stuke et al. 2020 Stuke, A.; Kunkel, C.; Golze, D.; Todorović, M.; Margraf, J. T.; Reuter, K.; Rinke, P.; Oberhofer, H. Atomic structures and orbital energies of 61,489 crystal-forming organic molecules. _Sci. Data_ 2020, _7_ , 58
  * Mobley and Guthrie 2014 Mobley, D. L.; Guthrie, J. P. FreeSolv: a database of experimental and calculated hydration free energies, with input files. _J. Comput. Aided Mol. Des._ 2014, _28_ , 711–720
  * Schütt et al. 2019 Schütt, K. T.; Kessel, P.; Gastegger, M.; Nicoli, K. A.; Tkatchenko, A.; Müller, K.-R. SchNetPack: A Deep Learning Toolbox For Atomistic Systems. _J. Chem. Theory Comput._ 2019, _15_ , 448–455
  * Himanen et al. 2020 Himanen, L.; Jäger, M. O. J.; Morooka, E. V.; Federici Canova, F.; Ranawat, Y. S.; Gao, D. Z.; Rinke, P.; Foster, A. S. DScribe: Library of descriptors for machine learning in materials science. _Comput. Phys. Commun._ 2020, _247_ , 106949
  * Thiel 2014 Thiel, W. Semiempirical quantum–chemical methods. _WIREs Comput. Mol. Sci._ 2014, _4_ , 145–157
  * Ramakrishnan et al. 2014 Ramakrishnan, R.; Dral, P. O.; Rupp, M.; von Lilienfeld, O. A. Quantum chemistry structures and properties of 134 kilo molecules. _Sci. Data_ 2014, _1_ , 140022
  * Neese et al. 2020 Neese, F.; Wennmohs, F.; Becker, U.; Riplinger, C. The ORCA quantum chemistry program package. _J. Chem. Phys._ 2020, _152_ , 224108
  * Bannwarth et al. 2019 Bannwarth, C.; Ehlert, S.; Grimme, S. GFN2-xTB—An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions. _J. Chem. Theory Comput._ 2019, _15_ , 1652–1671
  * Bannwarth et al. 2021 Bannwarth, C.; Caldeweyher, E.; Ehlert, S.; Hansen, A.; Pracht, P.; Seibert, J.; Spicher, S.; Grimme, S. Extended tight-binding quantum chemistry methods. _WIREs Comput. Mol. Sci._ 2021, _11_ , e1493
  * Barone and Cossi 1998 Barone, V.; Cossi, M. Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model. _J. Phys. Chem. A_ 1998, _102_ , 1995–2001
  * Onufriev et al. 2004 Onufriev, A.; Bashford, D.; Case, D. A. Exploring protein native states and large-scale conformational changes with a modified generalized born model. _Proteins_ 2004, _55_ , 383–394
  * Sigalov et al. 2006 Sigalov, G.; Fenley, A.; Onufriev, A. Analytical electrostatics for biomolecules: Beyond the generalized Born approximation. _J. Chem. Phys._ 2006, _124_ , 124902
  * Lange and Herbert 2012 Lange, A. W.; Herbert, J. M. Improving Generalized Born Models by Exploiting Connections to Polarizable Continuum Models. I. An Improved Effective Coulomb Operator. _J. Chem. Theory Comput._ 2012, _8_ , 1999–2011
  * Paul et al. 2019 Paul, A.; Furmanchuk, A.; Liao, W.-k.; Choudhary, A.; Agrawal, A. Property Prediction of Organic Donor Molecules for Photovoltaic Applications Using Extremely Randomized Trees. _Mol. Inform._ 2019, _38_ , 1900038
  * Zhao et al. 2020 Zhao, Z.-W.; del Cueto, M.; Geng, Y.; Troisi, A. Effect of Increasing the Descriptor Set on Machine Learning Prediction of Small Molecule-Based Organic Solar Cells. _Chem. Mater._ 2020, _32_ , 7777–7787
  * Munshi et al. 2021 Munshi, J.; Chen, W.; Chien, T.; Balasubramanian, G. Transfer Learned Designer Polymers For Organic Solar Cells. _J. Chem. Inf. Model._ 2021, _61_ , 134–142
  * Eibeck et al. 2021 Eibeck, A.; Nurkowski, D.; Menon, A.; Bai, J.; Wu, J.; Zhou, L.; Mosbach, S.; Akroyd, J.; Kraft, M. Predicting Power Conversion Efficiency of Organic Photovoltaics: Models and Data Analysis. _ACS Omega_ 2021, _6_ , 23764–23775
  * Fare et al. 2022 Fare, C.; Fenner, P.; Benatan, M.; Varsi, A.; Pyzer-Knapp, E. O. A multi-fidelity machine learning approach to high throughput materials screening. _npj Comput. Mater._ 2022, _8_ , 1–9
  * Moore et al. 2022 Moore, G. J.; Bardagot, O.; Banerji, N. Deep Transfer Learning: A Fast and Accurate Tool to Predict the Energy Levels of Donor Molecules for Organic Photovoltaics. _Adv. Theory Simul._ 2022, _5_ , 2100511
  * Li et al. 2022 Li, H.; Zhao, X.; Li, S.; Wan, F.; Zhao, D.; Zeng, J. Improving molecular property prediction through a task similarity enhanced transfer learning strategy. _iScience_ 2022, _25_ , 105231
  * Hutchinson and Kobayashi 2019 Hutchinson, S. T.; Kobayashi, R. Solvent-Specific Featurization for Predicting Free Energies of Solvation through Machine Learning. _J. Chem. Inf. Model._ 2019, _59_ , 1338–1346
  * Scheen et al. 2020 Scheen, J.; Wu, W.; Mey, A. S. J. S.; Tosco, P.; Mackey, M.; Michel, J. Hybrid Alchemical Free Energy/Machine-Learning Methodology for the Computation of Hydration Free Energies. _J. Chem. Inf. Model._ 2020, _60_ , 5331–5339
  * Weinreich et al. 2021 Weinreich, J.; Browning, N. J.; von Lilienfeld, O. A. Machine learning of free energies in chemical compound space using ensemble representations: Reaching experimental uncertainty for solvation. _J. Chem. Phys._ 2021, _154_ , 134113
  * Borgis et al. 2021 Borgis, D.; Luukkonen, S.; Belloni, L.; Jeanmairet, G. Accurate prediction of hydration free energies and solvation structures using molecular density functional theory with a simple bridge functional. _J. Chem. Phys._ 2021, _155_ , 024117
  * Vermeire and Green 2021 Vermeire, F. H.; Green, W. H. Transfer learning for solvation free energies: From quantum chemistry to experiments. _Chem. Eng. J._ 2021, _418_ , 129307
  * Gao et al. 2021 Gao, P.; Yang, X.; Tang, Y.-H.; Zheng, M.; Andersen, A.; Murugesan, V.; Hollas, A.; Wang, W. Graphical Gaussian process regression model for aqueous solvation free energy prediction of organic molecules in redox flow batteries. _Phys. Chem. Chem. Phys._ 2021, _23_ , 24892–24904
  * Pathak et al. 2021 Pathak, Y.; Mehta, S.; Priyakumar, U. D. Learning Atomic Interactions through Solvation Free Energy Prediction Using Graph Neural Networks. _J. Chem. Inf. Model._ 2021, _61_ , 689–698
  * Chung et al. 2022 Chung, Y.; Vermeire, F. H.; Wu, H.; Walker, P. J.; Abraham, M. H.; Green, W. H. Group Contribution and Machine Learning Approaches to Predict Abraham Solute Parameters, Solvation Free Energy, and Solvation Enthalpy. _J. Chem. Inf. Model._ 2022, _62_ , 433–446
  * Low et al. 2022 Low, K.; Coote, M. L.; Izgorodina, E. I. Explainable Solvation Free Energy Prediction Combining Graph Neural Networks with Chemical Intuition. _J. Chem. Inf. Model._ 2022, _62_ , 5457–5470
  * Zhang et al. 2022 Zhang, D.; Xia, S.; Zhang, Y. Accurate Prediction of Aqueous Free Solvation Energies Using 3D Atomic Feature-Based Graph Neural Network with Transfer Learning. _J. Chem. Inf. Model._ 2022, _62_ , 1840–1848
  * Zhang et al. 2023 Zhang, Z.-Y.; Peng, D.; Liu, L.; Shen, L.; Fang, W.-H. Machine Learning Prediction of Hydration Free Energy with Physically Inspired Descriptors. _J. Phys. Chem. Lett._ 2023, _14_ , 1877–1884
  * Bennett et al. 2020 Bennett, L.; Melchers, B.; Proppe, B. _Curta: A General-purpose High-Performance Computer at ZEDAT, Freie Universität Berlin_ ; Freie Universität Berlin, Germany, 2020


CTP-SCU/2021029, APCTP Pre2020 - 020

# Linear seesaw model with a modular $S_{4}$ flavor symmetry

Takaaki Nomura<EMAIL_ADDRESS>College of Physics, Sichuan University,
Chengdu 610065, China  Hiroshi Okada<EMAIL_ADDRESS>Asia Pacific
Center for Theoretical Physics (APCTP) - Headquarters San 31, Hyoja-dong, Nam-
gu, Pohang 790-784, Korea  Department of Physics, Pohang University of Science
and Technology, Pohang 37673, Republic of Korea

###### Abstract

We discuss a linear seesaw model with as minimum field content as possible,
introducing a modular $S_{4}$ with the help of gauged $U(1)_{B-L}$ symmetries.
Due to rank two neutrino mass matrix, we have a vanishing neutrino mass
eigenvalue, and only the normal mass hierarchy of neutrinos is favored through
the modular $S_{4}$ symmetry. In our numerical $\Delta\chi^{2}$ analysis, we
especially find rather sharp prediction on sum of neutrino masses to be around
$60$ meV in addition to the other predictions.

##  I Introduction

Neutrino sector is theoretically unconfirmed in the standard model (SM),
because only two mass squared differences and three mixings are experimentally
found and the scale of mass is extremely minuscule compared to the other three
sectors in the SM. Thus, a lot of scientists expect that neutrino sector would
possess new physics. A gauged $B-L$(baryon number minus lepton number) $U(1)$
symmetry; $U(1)_{B-L}$, is one of the promising prescriptions to generate such
tiny neutrino masses introducing three right-handed neutrinos with rather
heavy masses ($M_{R}$), that is called canonical seesaw model Seesaw1 ;
Seesaw2 ; Seesaw3 ; Seesaw4 . Since the mass scale is sometime expected to be
the one of a grand unified theory($M_{GUT}\sim 10^{15}$ GeV) in order to be
small neutrino masses, its scale cannot be confirmed by our current
experiments, and the spontaneous $U(1)_{B-L}$ symmetry breaking scale is
naturally expected to be the same energy as the cut-off scale; $M_{R}\sim
M_{GUT}$.

In order to achieve a successful neutrino mass model within our
scale$(\sim$TeV), another heavy neutral fermions($S_{L}$) with left-handed
chirality are introduced along this line of idea 111Inverse seesaw model
Mohapatra:1986bd ; Wyler:1982dd  also requests both of $N_{R}$ and $S_{L}$ and
the neutrino mass could be realized within TeV scale. But this model may not
require GUT scale.. Now that there exist $N_{R}$ and $S_{L}$, these can be
embedded into the middle scale with $SU(3)_{C}\otimes SU(2)_{L}\otimes
SU(2)_{R}\otimes U(1)_{B-L}$ Wyler:1982dd ; Akhmedov:1995ip ; Akhmedov:1995vm
, moreover, which is included in $SO(10)$ group Malinsky:2005bi . Then, the
neutrino mass matrix would not be suppressed by the middle scale but by
$M_{GUT}$ scale when appropriate charge assignments are assigned for each
fields in a supersymmetric theory Malinsky:2005bi . Thus, we can test
phenomenologies with our current experiments, supposing the middle scale
breaking can occur at TeV scale. This type of model is called by ”linear
seesaw”. 222Notice here that our model has a different mechanism from the one
of original linear seesaw, even though our neutrino model can be realized
within TeV scale. Therefore, our neutrino mass matrix directly depends on the
scale of $U(1)_{B-L}$ breaking and ratio between vacuum expectation values
(VEVs) of two Higgs doublet model as can be seen in the main text, while we
expect higher energy scale such as GUT in order that the modular field must
break and get VEV denoted by $\tau$ in our literature.  Since these models
typically require more free parameters than the other three sectors in the SM
fermion, flavor symmetries are also introduced in these models frequently in
order to reduce the parameters and get predictions (if possible).

In 2017, attractive flavor symmetries are proposed by papers Feruglio:2017spp
; deAdelhartToorop:2011re , in which they applied modular non-Abelian discrete
flavor symmetries to quark and lepton sectors. One remarkable advantage is
that any dimensionless couplings can also be transformed as non-trivial
representations under those symmetries. Therefore, we do not need so many
scalars to find a predictive mass matrix. Another advantage is that we have a
modular weight from the modular origin that can play a role in stabilizing DM
when appropriate charge assignments are distributed to each of the fields of
models. Along the line of this idea, a vast reference has recently appeared in
the literature, e.g., $A_{4}$ Feruglio:2017spp ; Criado:2018thu ;
Kobayashi:2018scp ; Okada:2018yrn ; Nomura:2019jxj ; Okada:2019uoy ;
deAnda:2018ecu ; Novichkov:2018yse ; Nomura:2019yft ; Okada:2019mjf ;
Ding:2019zxk ; Nomura:2019lnr ; Kobayashi:2019xvz ; Asaka:2019vev ;
Zhang:2019ngf ; Gui-JunDing:2019wap ; Kobayashi:2019gtp ; Nomura:2019xsb ;
Wang:2019xbo ; Okada:2020dmb ; Okada:2020rjb ; Behera:2020lpd ; Behera:2020sfe
; Nomura:2020opk ; Nomura:2020cog ; Asaka:2020tmo ; Okada:2020ukr ;
Nagao:2020snm ; Okada:2020brs ; Yao:2020qyy ; Chen:2021zty ; Kashav:2021zir ;
Okada:2021qdf ; deMedeirosVarzielas:2021pug ; Nomura:2021yjb ;
Hutauruk:2020xtk ; Ding:2021eva ; Nagao:2021rio , $S_{3}$ Kobayashi:2018vbk ;
Kobayashi:2018wkl ; Kobayashi:2019rzp ; Okada:2019xqk ; Mishra:2020gxg ;
Du:2020ylx , $S_{4}$ Penedo:2018nmg ; Novichkov:2018ovf ; Kobayashi:2019mna ;
King:2019vhv ; Okada:2019lzv ; Criado:2019tzk ; Wang:2019ovr ; Zhao:2021jxg ;
King:2021fhl ; Ding:2021zbg ; Zhang:2021olk ; gui-jun , $A_{5}$
Novichkov:2018nkm ; Ding:2019xna ; Criado:2019tzk , double covering of $A_{5}$
Wang:2020lxk ; Yao:2020zml ; Wang:2021mkw ; Behera:2021eut , larger groups
Baur:2019kwi , multiple modular symmetries deMedeirosVarzielas:2019cyj , and
double covering of $A_{4}$ Liu:2019khw ; Chen:2020udk ; Li:2021buv , $S_{4}$
Novichkov:2020eep ; Liu:2020akv , and the other types of groups
Kikuchi:2020nxn ; Almumin:2021fbk ; Ding:2021iqp ; Feruglio:2021dte ;
Kikuchi:2021ogn ; Novichkov:2021evw  in which masses, mixing, and CP phases
for the quark and/or lepton have been predicted 333For interested readers, we
provide some literature reviews, which are useful to understand the non-
Abelian group and its applications to flavor structure Altarelli:2010gt ;
Ishimori:2010au ; Ishimori:2012zz ; Hernandez:2012ra ; King:2013eh ;
King:2014nza ; King:2017guk ; Petcov:2017ggy .. Moreover, a systematic
approach to understanding the origin of CP transformations has been discussed
in Ref. Baur:2019iai , and CP/flavor violation in models with modular symmetry
was discussed in Refs. Kobayashi:2019uyt ; Novichkov:2019sqv , and a possible
correction from Kähler potential was discussed in Ref. Chen:2019ewa .
Furthermore, systematic analysis of the fixed points (stabilizers) has been
discussed in Ref. deMedeirosVarzielas:2020kji . It would be interesting to
consider linear seesaw model with local $U(1)_{B-L}(\subset SO(10))$ under
modular symmetry since these symmetries can be originated from a string
theory. Moreover nature of modular symmetry can be used to realize linear
seesaw mechanism in addition to constraining flavor structure.

In this study, we propose a linear seesaw model under modular $S_{4}$ with the
help of $U(1)_{B-L}$ symmetry, in which we try to construct the predictive
model as minimum as possible. Due to rank two neutrino mass matrix, we have a
vanishing neutrino mass eigenvalue. Furthermore, only the normal mass
hierarchy of neutrinos is favored through the modular $S_{4}$ symmetry. In our
numerical, we perform $\Delta\chi^{2}$ analysis in the neutrino sector,
considering non-unitarity constraint.

This paper is organized as follows. In Sec. II, we review our model,
constructing renormalizable Lagrangian and mass matrices in lepton sector.
Then, we formulate the neutrino mass matrix with rank two, in which we
estimate structure of the neutrino mass matrix in expansion of modulus. And we
derive several observables in lepton sector. At the end of this section we
discuss the non-unitarity bound. In Sec. III we perform $\Delta\chi^{2}$
analysis in the lepton sector, and show some predictions through our model. In
Sec. IV, we give summary and discussion. In Appendix, we explain the modular
$S_{4}$ symmetry.

##  II Model

| Fermions | Scalars
---|---|---
|  $L_{L}$ |  $\overline{e_{R}},[\overline{\mu_{R}},\overline{\tau_{R}}]$ |  $\overline{N_{R}}$ |  $S_{L_{1}},S_{L_{2}}$ |  $H_{1}$ |  $H_{2}$ |  $\varphi$
$SU(2)_{L}$ | $\bm{2}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{2}$ | $\bm{2}$ | $\bm{1}$
$U(1)_{Y}$ | $-\frac{1}{2}$ | $1$ | $0$ | $0$ | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$
$U(1)_{B-L}$ | $-1$ | $1$ | $1$ | $0$ | $0$ | $1$ | $-1$
$S_{4}$ | ${\bf 3^{\prime}}$ | ${\bf 1},{\bf 2}$ | $\bm{3}$ | $\bm{1},\bm{1^{\prime}}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$
$-k_{I}$ | $-1$ | $-1,-3$ | $-3$ | $-1$ | $0$ | $0$ | $0$
Table 1: Lepton and boson particle contents and their charge assignments under
$SU(2)_{L}\times U(1)_{Y}\times U(1)_{B-L}\times S_{4}\times(-k_{I})$ where
$L_{L}\equiv[L_{L_{e}},L_{L_{\mu}},L_{L_{\tau}}]^{T}$ $k_{I}$ is the number of
modular weight.

###  II.1 Model review

In this section we review our model framework for linear seesaw mechanism,
introducing $B-L$ local Abelian symmetry; $U(1)_{B-L}$, and modular $A_{4}$
symmetry. As for fermion sector, we add two left-handed neutral fermions
$\\{S_{L_{1}},S_{L_{2}}\\}$ that belong to isospin singlet, where they have
zero charge under $U(1)_{B-L}$, $\\{\bm{1},\bm{1^{\prime}}\\}$ under $S_{4}$,
and $-1$ under $-k_{I}$, respectively. Also, we introduce three right-handed
neutral fermions $\overline{N_{R}}$ that belong to isospin singlet, where they
have $1$ charge under $U(1)_{B-L}$, $\bm{3}$ under $S_{4}$, and $-3$ under
$-k$, respectively. The SM left-handed leptons
$L_{L}\equiv[L_{L_{e}},L_{L_{\mu}},L_{L_{\tau}}]^{T}$ belong to $-1$ charge
under $U(1)_{B-L}$, $\bm{3^{\prime}}$ under $S_{4}$, and $-1$ under $-k$,
respectively. While the SM right-handed leptons
$\\{\overline{e_{R}},[\overline{\mu_{R}},\overline{\tau_{R}}]\\}$ belong to
$+1$ charge under $U(1)_{B-L}$, $\\{{\bf 1},{\bf 2}\\}$ under $S_{4}$, and
$\\{-1,-3\\}$ under $-k_{I}$, respectively.

As for scalar sector, we adopt two Higgs doublet $H_{1},H_{2}$ plus an isospin
singlet field $\varphi$. An isospin singlet $\varphi$ has $-1$ charge under
the $U(1)_{B-L}$. Here, $H_{1}$ is SM-like Higgs that has zero charge under
$U(1)_{B-L}$ and $-k_{I}$ while $H_{2}$ has $1$ charge under $U(1)_{B-L}$, and
zero under $-k_{I}$. We denote each of vacuum expectation values (VEVs) to be
$\langle H_{1,2}\rangle\equiv[0,v_{1,2}/\sqrt{2}]^{T}$, and
$\langle\varphi\rangle\equiv v_{\varphi}/\sqrt{2}$. We summarize our particle
content and their assignments in Table 1.

Then the valid lepton Yukawa Lagrangian is symbolized by

$\displaystyle-\mathcal{L}_{\rm
lepton}=\mathcal{L}_{M_{\ell}}+\mathcal{L}_{\rm M_{D}}+\mathcal{L}_{\rm
M^{\prime}_{D}}+\mathcal{L}_{M_{NS}},$ (II.1)

where $\mathcal{L}_{M_{\ell}}$ is charged lepton Yukawa Lagrangian.
$\mathcal{L}_{\rm M_{D}}$ is the one of [$\overline{N_{R}}L_{L}\tilde{H}_{1}$]
where $\tilde{H}\equiv i\sigma_{2}H^{*}$. $\mathcal{L}_{\rm M^{\prime}_{D}}$
is the one of [$\overline{L_{L}^{C}}S_{L}H_{2}$]. $\mathcal{L}_{M_{NS}}$ is
the one of [$\overline{N_{R}}S_{L}\varphi$]. $[\cdots]$ implies that concrete
flavor structures are manifolded. We will see each of their structures below.

The scalar potential of our model is written by

$\displaystyle V$
$\displaystyle=m_{\varphi}^{2}\varphi^{*}\varphi+m_{1}^{2}H_{1}^{\dagger}H_{1}+m_{2}^{2}H_{2}^{\dagger}H_{2}-\mu_{12}(H_{1}^{\dagger}H_{2}\varphi+h.c.)$
$\displaystyle+\lambda_{1}(H_{1}^{\dagger}H_{1})^{2}+\lambda_{2}(H_{2}^{\dagger}H_{2})^{2}+\lambda_{\varphi}(\varphi^{*}\varphi)^{2}+\lambda_{3}(H_{1}^{\dagger}H_{1})(H_{2}^{\dagger}H_{2})+\lambda_{4}(H_{1}^{\dagger}H_{2})(H_{2}^{\dagger}H_{1})$
$\displaystyle+\lambda_{\varphi
H_{1}}(H_{1}^{\dagger}H_{1})(\varphi^{*}\varphi)+\lambda_{\varphi
H_{2}}(H_{2}^{\dagger}H_{2})(\varphi^{*}\varphi),$ (II.2)

where $h.c.$ stands for Hermitian conjugate. We consider $\varphi$ develops a
VEV at much higher scale than $H_{1,2}$. Then after $\varphi$ developing a VEV
the scalar potential matches with that of two Higgs doublet model (THDM)
without $(H_{1}^{\dagger}H_{2})^{2}$ term due to $U(1)_{B-L}$ symmetry. Also
the Yukawa couplings associated with two Higgs doublet are those of type-I
THDM since only $H_{1}$ can couples to the SM fermions. In our analysis below
we thus do not discuss THDM part further and focus on neutrino mass.

###  II.2 Valid Lagrangians

Before discussing the valid Lagrangians for lepton sector, we define Yukawa
couplings under modular $S_{4}$ symmetry as follows: $S_{4}$ doublet with
$-k_{I}=2$ and $4$ are respectively denoted by $Y^{(2)}_{\bf
2}\equiv[y_{1},y_{2}]^{T}$ and $[y^{\prime}_{1},y^{\prime}_{2}]^{T}$, $S_{4}$
triplet with $-k_{I}=2$ is $Y^{(2)}_{\bf
3^{\prime}}\equiv[y_{3},y_{4},y_{5}]^{T}$, $S_{4}$ triplets with $-k_{I}=4$
are $Y^{(4)}_{\bf 3}\equiv[y^{\prime}_{3},y^{\prime}_{4},y^{\prime}_{5}]^{T}$,
and $Y^{(4)}_{\bf
3^{\prime}}\equiv[y^{\prime\prime}_{3},y^{\prime\prime}_{4},y^{\prime\prime}_{5}]^{T}$.
Each of structures is explicitly written in Appendix.

Charged lepton mass matrix:
The renormalizable Lagrangian for charged-lepton sector is given by

$\displaystyle-\mathcal{L}_{M_{\ell}}$
$\displaystyle=\alpha_{\ell}\overline{e_{R}}(y_{3}L_{L_{e}}+y_{4}\bar{L}_{L_{\tau}}+y_{5}\bar{L}_{L_{\mu}})\tilde{H}_{1}$
$\displaystyle+\beta_{\ell}\left[\frac{\sqrt{3}}{2}\overline{\mu_{R}}(y^{\prime}_{4}L_{L_{\mu}}+y^{\prime}_{5}L_{L_{\tau}})+\overline{\tau_{R}}\left(-y^{\prime}_{3}L_{L_{e}}+\frac{1}{2}(y^{\prime}_{4}L_{L_{\tau}}+y^{\prime}_{5}L_{L_{\mu}})\right)\right]\tilde{H}_{1}$
$\displaystyle+\gamma_{\ell}\left[\frac{\sqrt{3}}{2}\overline{\tau_{R}}(y^{\prime\prime}_{4}L_{L_{\mu}}+y^{\prime\prime}_{5}L_{L_{\tau}})+\overline{\mu_{R}}\left(y^{\prime\prime}_{3}L_{L_{e}}-\frac{1}{2}(y^{\prime\prime}_{4}L_{L_{\tau}}+y^{\prime\prime}_{5}L_{L_{\mu}})\right)\right]\tilde{H}_{1}+{\rm
h.c.},$ (II.3)

where $\\{\alpha_{\ell},\beta_{\ell},\gamma_{\ell}\\}$ are real parameters
without loss of generality. Then the charged-lepton mass matrix after the
spontaneous symmetry breaking is given by

$\displaystyle(M_{\ell})_{RL}=\frac{v_{1}}{\sqrt{2}}\begin{pmatrix}\alpha_{\ell}y_{3}&\alpha_{\ell}y_{5}&\alpha_{\ell}y_{4}\\\
\gamma_{\ell}y^{\prime\prime}_{3}&\frac{\sqrt{3}}{2}\beta_{\ell}y^{\prime}_{4}-\frac{1}{2}\gamma_{\ell}y^{\prime\prime}_{5}&\frac{\sqrt{3}}{2}\beta_{\ell}y^{\prime}_{5}-\frac{1}{2}\gamma_{\ell}y^{\prime\prime}_{4}\\\
-\beta_{\ell}y^{\prime}_{3}&\frac{\sqrt{3}}{2}\gamma_{\ell}y^{\prime\prime}_{4}+\frac{1}{2}\beta_{\ell}y^{\prime}_{5}&\frac{\sqrt{3}}{2}\gamma_{\ell}y^{\prime\prime}_{5}+\frac{1}{2}\beta_{\ell}y^{\prime}_{4}\\\
\end{pmatrix}.$ (II.4)

The charged-lepton mass eigenvalues are obtained by diagonalizing ${\rm
diag}[m_{e},m_{\mu},m_{\tau}]=V_{R_{\ell}}^{\dagger}M_{\ell}V_{L_{\ell}}$,
where $V_{L_{\ell},R_{\ell}}$ are unitary matrices. In our numerical analysis,
we will determine the free parameters
$\\{\alpha_{\ell},\beta_{\ell},\gamma_{\ell}\\}$ so as to fit the three
charged-lepton mass eigenstates after giving all the numerical values, by
applying the relations: Here, we fix
$\alpha_{\ell},\beta_{\ell},\gamma_{\ell}$ in order to find the experimental
three charged-lepton masses, applying the following relations:

$\displaystyle{\rm
Tr}[{M_{\ell}}^{\dagger}M_{\ell}]=|m_{e}|^{2}+|m_{\mu}|^{2}+|m_{\tau}|^{2},\quad{\rm
Det}[{M_{\ell}}^{\dagger}M_{\ell}]=|m_{e}|^{2}|m_{\mu}|^{2}|m_{\tau}|^{2},$
$\displaystyle({\rm Tr}[{M_{\ell}}^{\dagger}]M_{\ell})^{2}-{\rm
Tr}[({M_{\ell}}^{\dagger}M_{\ell})^{2}]=2(|m_{e}|^{2}|m_{\mu}|^{2}+|m_{\mu}|^{2}|m_{\tau}|^{2}+|m_{e}|^{2}|m_{\tau}|^{2}).$
(II.5)

Neutral fermion mass matrices:
At first, we construct the valid Lagrangian of Dirac mass matrix
$\mathcal{L}_{\rm M_{D}}$. This is given by

$\displaystyle\mathcal{L}_{\rm M_{D}}$
$\displaystyle=\alpha_{D}\left[y^{\prime}_{3}(\overline{N_{R_{3}}}L_{L_{\mu}}-\overline{N_{R_{2}}}L_{L_{\tau}})+y^{\prime}_{4}(\overline{N_{R_{1}}}L_{L_{\tau}}-\overline{N_{R_{3}}}L_{L_{e}})+y^{\prime}_{5}(\overline{N_{R_{2}}}L_{L_{e}}-\overline{N_{R_{1}}}L_{L_{\mu}})\right]\tilde{H}_{1}$
$\displaystyle+\beta_{D}\left[y^{\prime\prime}_{3}(\overline{N_{R_{3}}}L_{L_{\tau}}-\overline{N_{R_{2}}}L_{L_{\mu}})+y^{\prime\prime}_{4}(-\overline{N_{R_{1}}}L_{L_{\mu}}-\overline{N_{R_{2}}}L_{L_{e}})+y^{\prime\prime}_{5}(\overline{N_{R_{1}}}L_{L_{\tau}}+\overline{N_{R_{3}}}L_{L_{e}})\right]\tilde{H}_{1}$
$\displaystyle+\gamma_{D}\left[\frac{\sqrt{3}}{2}y^{\prime}_{1}(\overline{N_{R_{2}}}L_{L_{\mu}}+\overline{N_{R_{3}}}L_{L_{\tau}})+y^{\prime}_{2}\\{-\overline{N_{R_{1}}}L_{L_{e}}+\frac{1}{2}(\overline{N_{R_{2}}}L_{L_{\tau}}+\overline{N_{R_{3}}}L_{L_{\mu}})\\}\right]\tilde{H}_{1}+{\rm
h.c.},$ (II.6)

where we suppose $\alpha_{D}$ to be real while $\beta_{D},\ \gamma_{D}$ to be
complex after rephasing of fields. Similar to the charged-lepton sector we
find the Dirac mass matrix as follows:

$\displaystyle(M_{D})_{RL}=\frac{v_{1}}{\sqrt{2}}\begin{pmatrix}-\gamma_{D}y^{\prime}_{2}&-\alpha_{D}y^{\prime}_{5}-\beta_{D}y^{\prime\prime}_{4}&\alpha_{D}y^{\prime}_{4}+\beta_{D}y^{\prime\prime}_{5}\\\
\alpha_{D}y^{\prime}_{5}-\beta_{D}y^{\prime\prime}_{4}&-\beta_{D}y^{\prime\prime}_{3}+\frac{\sqrt{3}}{2}\gamma_{D}y^{\prime}_{1}&-\alpha_{D}y^{\prime}_{3}+\frac{1}{2}\gamma_{D}y^{\prime}_{2}\\\
-\alpha_{D}y^{\prime}_{4}+\beta_{D}y^{\prime\prime}_{5}&\alpha_{D}y^{\prime}_{3}+\frac{1}{2}\gamma_{D}y^{\prime}_{2}&\beta_{D}y^{\prime\prime}_{3}+\frac{\sqrt{3}}{2}\gamma_{D}y^{\prime}_{1}\\\
\end{pmatrix}.$ (II.7)

For our convenience to analyze the neutrino oscillation, we redefine
$M_{D}\equiv\frac{v_{1}}{\sqrt{2}}\tilde{M_{D}}$

Another Dirac Lagrangian is induced via $\mathcal{L}_{\rm M^{\prime}_{D}}$,
which is given by

$\displaystyle\mathcal{L}_{\rm M^{\prime}_{D}}$
$\displaystyle=\alpha^{\prime}_{D}\left[y_{3}\overline{L^{C}_{L_{e}}}+y_{4}\overline{L^{C}_{L_{\tau}}}+y_{5}\overline{L^{C}_{L_{\mu}}}\right]S_{L_{1}}H_{2}+{\rm
h.c.},$ (II.8)

Then, we find another Dirac mass matrix as follows:

$\displaystyle(M^{\prime}_{D})_{L_{L}S_{L}}=\frac{v_{2}}{\sqrt{2}}\alpha^{\prime}_{D}\begin{pmatrix}y_{3}&0\\\
y_{5}&0\\\ y_{4}&0\\\ \end{pmatrix}.$ (II.9)

Being the same as the reason for $M_{D}$, we redefine
$M^{\prime}_{D}\equiv\frac{v_{2}}{\sqrt{2}}\alpha^{\prime}_{D}\tilde{M^{\prime}_{D}}$
.

The third term of Lagrangian $\mathcal{L}_{M_{NS}}$ is given by

$\displaystyle\mathcal{L}_{M_{NS}}=$
$\displaystyle\alpha_{NS}\left[y^{\prime}_{3}\overline{N_{R_{e}}}+y^{\prime}_{4}\overline{N_{R_{\tau}}}+y^{\prime}_{5}\overline{N_{R_{\mu}}}\right]S_{L_{1}}\varphi+\beta_{NS}\left[y^{\prime\prime}_{3}\overline{N_{R_{e}}}+y^{\prime\prime}_{4}\overline{N_{R_{\tau}}}+y^{\prime\prime}_{5}\overline{N_{R_{\mu}}}\right]S_{L_{2}}\varphi+{\rm
h.c.},$ (II.10)

where $\alpha_{NS},\beta_{NS}$ are real without loss of generality. Then, we
obtain the mass matrix

$\displaystyle
M_{NS}=\frac{v_{\varphi}}{\sqrt{2}}\begin{pmatrix}y^{\prime}_{3}&y^{\prime\prime}_{3}\\\
y^{\prime}_{5}&y^{\prime\prime}_{5}\\\ y^{\prime}_{4}&y^{\prime\prime}_{4}\\\
\end{pmatrix}\begin{pmatrix}\alpha_{NS}&0\\\ 0&\beta_{NS}\\\
\end{pmatrix}\equiv\frac{v_{\varphi}}{\sqrt{2}}\tilde{M}_{NS}.$ (II.11)

In basis of $[\nu_{L},N_{R}^{C},S_{L}]^{T}$, the neutral fermion mass matrix
is given by

$\displaystyle M_{N}=\begin{pmatrix}0_{3\times 3}&M_{D}^{T}&M^{\prime}_{D}\\\
M_{D}&0_{3\times 3}&M_{NS}\\\ m^{\prime T}_{D}&M_{NS}^{T}&0_{2\times
2}\end{pmatrix}.$ (II.12)

Then, the active neutrino mass matrix is given by

$\displaystyle m_{\nu}$
$\displaystyle=M^{\prime}_{D}(M_{NS}^{T}M_{NS})^{-1}M_{NS}^{T}M_{D}+[M^{\prime}_{D}(M_{NS}^{T}M_{NS})^{-1}M_{NS}^{T}M_{D}]^{T}$
$\displaystyle=\frac{v_{1}v_{2}}{\sqrt{2}v_{\varphi}}\left(\tilde{M}^{\prime}_{D}(\tilde{M}_{NS}^{T}\tilde{M}_{NS})^{-1}\tilde{M}_{NS}^{T}\tilde{M}_{D}+[\tilde{M}^{\prime}_{D}(\tilde{M}_{NS}^{T}\tilde{M}_{NS})^{-1}\tilde{M}_{NS}^{T}\tilde{M}_{D}]^{T}\right)$
$\displaystyle=\kappa\tilde{m}_{\nu},$ (II.13)

where $\kappa\equiv\frac{v_{1}v_{2}}{\sqrt{2}v_{\varphi}}$ and we assume mass
hierarchies among $M_{D},M^{\prime}_{D}\ll M_{NS}$. 444Mass hierarchies is
dynamically achieved in refs. Wang:2015saa ; Das:2017ski .. The neutrino mass
eigenvalues are obtained as follows: $D_{\nu}=\kappa
D_{\nu}=U_{\nu}^{T}m_{\nu}U_{\nu}=\kappa U_{\nu}^{T}\tilde{m}_{\nu}U_{\nu}$,
where $U_{\nu}$ is a unitary matrix. Then, the Pontecorvo-Maki-Nakagawa-Sakata
(PMNS) matrix is given by $U_{PMNS}\equiv V^{\dagger}_{L_{\ell}}U_{\nu}$.
Notice that $m_{\nu}$ is rank two, thus the lightest neutrino mass is zero.
Here $\kappa$ is described by one experimental values and dimensionless
neutrino mass eigenstates as follows:

$\displaystyle{\rm(NH)}:\ \kappa^{2}=\frac{|\Delta m_{\rm
atm}^{2}|}{\tilde{D}_{\nu_{3}}^{2}},\quad{\rm(IH)}:\ \kappa^{2}=\frac{|\Delta
m_{\rm atm}^{2}|}{\tilde{D}_{\nu_{2}}^{2}},$ (II.14)

where $\Delta m_{\rm atm}^{2}$ is the atmospheric neutrino mass-squared
difference and NH(IH) stands for normal(inverted) ordering respectively.
Subsequently, the solar mass difference squared can be written in terms of
$\kappa$ as follows:

$\displaystyle{\rm(NH)}:\ \Delta m_{\rm
sol}^{2}={\kappa^{2}}{\tilde{D}_{\nu_{2}}^{2}},\quad{\rm(IH)}:\ \Delta m_{\rm
sol}^{2}={\kappa^{2}}({\tilde{D}_{\nu_{2}}^{2}-\tilde{D}_{\nu_{1}}^{2}}),$
(II.15)

which can be compared to the observed value. In other words, we explicitly
write the mass eigenvalues in terms of $\Delta m_{\rm atm}^{2}$ and $\Delta
m_{\rm sol}^{2}$ as follows:

$\displaystyle{\rm(NH)}:\ D_{\nu_{1}}^{2}=0,\ D_{\nu_{2}}^{2}=\Delta m_{\rm
sol}^{2},\ D_{\nu_{3}}^{2}=\Delta m_{\rm atm}^{2},$ (II.16)
$\displaystyle{\rm(IH)}:\ D_{\nu_{1}}^{2}=\Delta m_{\rm atm}^{2}-\Delta m_{\rm
sol}^{2},\ D_{\nu_{2}}^{2}=\Delta m_{\rm sol}^{2},\ D_{\nu_{3}}^{2}=0,$
(II.17)

which implies that NH is hierarchical but IH is degenerate, since $\Delta
m_{\rm sol}^{2}/\Delta m_{\rm atm}^{2}<<1$. Here, we expand
$|\tilde{m}_{\nu}|^{2}$ in terms of $q\equiv e^{2p\pi i\tau}<<1$. Then the
mass matrix is given by

$\displaystyle|\tilde{m}_{\nu}|^{2}\sim\begin{pmatrix}{\cal O}(1)&{\cal
O}(q)&0\\\ {\cal O}(q)&{\cal O}(1)&{\cal O}(q^{2})\\\ 0&{\cal
O}(q^{2})&0\end{pmatrix}.$ (II.18)

The ratio between two nonzero squared eigenvalues $R$ is estimated by

$\displaystyle R={\cal O}(q)<<1.$ (II.19)

It suggests that the neutrino mass eigenvalues tend to be hierarchical,
therefore NH is favored. In fact, we would not obtain the allowed region
within $3\sigma$ for IH in our numerical analysis. Thus, we focus on NH
hereafter.

In our model, PMNS matrix is parametrized by three mixing angle
$\theta_{ij}(i,j=1,2,3;i<j)$, one CP violating Dirac phase $\delta_{CP}$, and
one Majorana phase $\alpha_{21}$ as follows:

$U_{PMNS}=\begin{pmatrix}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta_{CP}}\\\
-s_{12}c_{23}-c_{12}s_{23}s_{13}e^{i\delta_{CP}}&c_{12}c_{23}-s_{12}s_{23}s_{13}e^{i\delta_{CP}}&s_{23}c_{13}\\\
s_{12}s_{23}-c_{12}c_{23}s_{13}e^{i\delta_{CP}}&-c_{12}s_{23}-s_{12}c_{23}s_{13}e^{i\delta_{CP}}&c_{23}c_{13}\end{pmatrix}\begin{pmatrix}1&0&0\\\
0&e^{i\frac{\alpha_{21}}{2}}&0\\\ 0&0&1\end{pmatrix},$ (II.20)

where $c_{ij}$ and $s_{ij}$ stand for $\cos\theta_{ij}$ and $\sin\theta_{ij}$
respectively. Then, these mixings are given in terms of the components of
$U_{PMNS}$ as follows:

$\displaystyle\sin^{2}\theta_{13}=|(U_{PMNS})_{13}|^{2},\quad\sin^{2}\theta_{23}=\frac{|(U_{PMNS})_{23}|^{2}}{1-|(U_{PMNS})_{13}|^{2}},\quad\sin^{2}\theta_{12}=\frac{|(U_{PMNS})_{12}|^{2}}{1-|(U_{PMNS})_{13}|^{2}}.$
(II.21)

Also we compute the Jarlskog invariant $J_{CP}$ that is derived from PMNS
matrix elements as follows:

$J_{CP}=\text{Im}[U_{e1}U_{\mu 2}U_{e2}^{*}U_{\mu
1}^{*}]=s_{23}c_{23}s_{12}c_{12}s_{13}c^{2}_{13}\sin\delta_{CP}.$ (II.22)

Majorana phase is estimated in terms of other invariant $I_{1}$ as follows:

$I_{1}=\text{Im}[U^{*}_{e1}U_{e2}]=c_{12}s_{12}c_{13}^{2}\sin\left(\frac{\alpha_{21}}{2}\right).$
(II.23)

In addition, the effective mass for the neutrinoless double beta decay is
written by

$\displaystyle{\rm(NH)}:\ \langle
m_{ee}\rangle=\kappa|\tilde{D}_{\nu_{2}}s^{2}_{12}c^{2}_{13}e^{i\alpha_{21}}+\tilde{D}_{\nu_{3}}s^{2}_{13}e^{-2i\delta_{CP}}|,$
(II.24) $\displaystyle{\rm(IH)}:\ \langle
m_{ee}\rangle=\kappa|\tilde{D}_{\nu_{1}}c^{2}_{12}c^{2}_{13}+\tilde{D}_{\nu_{2}}s^{2}_{12}c^{2}_{13}e^{i\alpha_{21}}|,$
(II.25)

where its value could be measured by KamLAND-Zen in future KamLAND-Zen:2016pfg
. We will adopt the neutrino experimental data at 3$\sigma$ interval in Nufit
5.0 Esteban:2018azc ; Nufit .

Non-unitarity:
Here, let us briefly discuss non-unitarity matrix $U^{\prime}_{PMNS}$. This is
typically parametrized by the form

$\displaystyle
U^{\prime}_{PMNS}\equiv\left(1-\frac{1}{2}F^{\dagger}F\right)U_{PMNS},$
(II.26)

where $F\equiv(M_{NS}^{T}M_{NS})^{-1}M_{NS}^{T}M_{D}$ is a hermitian matrix,
and $U^{\prime}_{PMNS}$ represents the deviation from the unitarity. Applying
global constraints Fernandez-Martinez:2016lgt , one finds Agostinho:2017wfs

$\displaystyle|FF^{\dagger}|\leq\left[\begin{array}[]{ccc}2.5\times
10^{-3}&2.4\times 10^{-5}&2.7\times 10^{-3}\\\ 2.4\times 10^{-5}&4.0\times
10^{-4}&1.2\times 10^{-3}\\\ 2.7\times 10^{-3}&1.2\times 10^{-3}&5.6\times
10^{-3}\\\ \end{array}\right].$ (II.30)

In our case,
$F\equiv(M_{NS}^{T}M_{NS})^{-1}M_{NS}^{T}M_{D}=\frac{v_{1}}{v_{\varphi}}(\tilde{M}_{NS}^{T}\tilde{M}_{NS})^{-1}\tilde{M}_{NS}^{T}\tilde{M}_{D}$.
Since $v_{\varphi}$ is freely taken to be large (while $v_{1}={\cal O}(100)$
GeV at most), we easily satisfy this bound.

##  III Numerical analysis

In this section, we carry out numerical $\Delta\chi^{2}$ analysis searching
for parameters satisfying neutrino oscillation data and non-unitarity
constraint, and show our predictions, where we use the best fit values of
charged-lepton masses. Notice here that we focus on NH since IH is disfavored
by analytical estimation as can be seen in the previous section.

In our numerical analysis, we randomly scan free parameters in following
ranges

$\displaystyle\\{\alpha_{D},|\beta_{D}|,|\gamma_{D}|,\alpha_{NS},\beta_{NS}\\}\in[10^{-5},1.0],\quad
t_{\beta}\in[10,100],$ (III.1)

where $\tau$ runs over the fundamental region, and $t_{\beta}\equiv
v_{1}/v_{2}$ and $\sqrt{v_{1}^{2}+v_{2}^{2}}=246$ GeV. Under the above
regions, we perform numerical analysis. Fig. 1 shows the correlation between
real part of $\tau$ and imaginary part of $\tau$, where the blue points are
allowed within 2, green ones within 3, and red one within 5 of
$\Delta\chi^{2}$ analysis for five accurately known observables $\Delta
m^{2}_{\rm atm},\Delta m^{2}_{\rm sol},s_{12}^{2},s_{23}^{2},s_{13}^{2}$ in
Nufit 5.0 Esteban:2018azc ; Nufit . The real part runs whole the range, but
the imaginary part is localized at the region of $[1.35-1.5]$.

Figure 1: The allowed region of modulus $\tau$ , where the blue points are
allowed within 2, green ones within 3, and red one within 5 of
$\Delta\chi^{2}$ analysis.

Fig. 2 demonstrates the correlation between the Majorana phase $\alpha_{21}$
and Dirac CP phase $\delta_{CP}$. The legend is the same as the case of Fig.1.
We find a clear feature on the region that $\alpha_{21}$ is within
$\\{0^{\circ}-120^{\circ},240^{\circ}-310^{\circ}\\}$, while $\delta_{CP}$ is
within $\\{70^{\circ}-290^{\circ}\\}$.

Figure 2: Correlation between a Majorana phase $\alpha_{21}$ and
$\delta_{CP}$, where the legend is the same as the case of Fig.1.

The upper figures in Fig. 3 show correlation between the effective mass for
the neutrinoless double beta decay $\langle m_{ee}\rangle$ and the sum of
neutrino masses $\sum m_{i}$ in unit of [eV], where the legend is the same as
the case of Fig.1. We find that $\langle m_{ee}\rangle$ is allowed within
$[0.001-0.004]$ eV. On the other hand $\sum m_{i}$ is restricted to be around
$0.06$ eV, which would be a sharp prediction of this model.

Figure 3: Predicted correlation between the effective mass for the
neutrinoless double beta decay $\langle m_{ee}\rangle$ and the sum of neutrino
masses $\sum m_{i}$ in unit of [eV], where the legend is the same as the case
of Fig.1.

Finally, we show a benchmark in table 2, where we select it so that
$\sqrt{\Delta\chi^{2}}$ is minimum. The mass matrices for dimensionless
neutrino and charged-lepton are found as

$\displaystyle\tilde{m}_{\nu}$
$\displaystyle=\left[\begin{array}[]{ccc}37.1331-17.9304i&-16.2576+55.2032i&35.2235+24.3462i\\\
-16.2576+55.2032i&-40.0223-38.1947i&-55.9229+44.0753i\\\
35.2235+24.3462i&-55.9229+44.0753i&-1.23231+74.8642i\\\ \end{array}\right],$
(III.5) $\displaystyle M_{\ell}$
$\displaystyle=\left[\begin{array}[]{ccc}0.00232+0.000467i&-0.00111+0.00314i&-0.00453+0.00274i\\\
-0.282+0.102i&-0.0524-0.469i&0.214-0.559i\\\
0.336-0.277i&0.362+0.706i&-0.160+1.23i\\\ \end{array}\right].$ (III.9)

$\tau$ | $0.113762+1.43906i$
---|---
$t_{\beta}$ | $98.6$
$[\alpha_{\ell},\gamma_{\ell},\beta_{\ell}]$ | $[1.72\times 10^{-5},\ 6.15\times 10^{-4},\ 8.06\times 10^{-4}]$
$[\alpha_{D},\alpha_{NS},\beta_{NS}]$ | $[-0.0112907,\ -0.00078512,\ 0.0203379]$
$[\beta_{D},\gamma_{D}]$ | $[-0.01432-0.00360i,-1.90\times 10^{-5}-2.34\times 10^{-6}i]$
$\Delta m^{2}_{\rm atm}$ | $2.53\times 10^{-3}{\rm eV}^{2}$
$\Delta m^{2}_{\rm sol}$ | $7.48\times 10^{-5}{\rm eV}^{2}$
$\sin^{2}\theta_{12}$ | $0.289$
$\sin^{2}\theta_{23}$ | $0.565$
$\sin^{2}\theta_{13}$ | $0.02207$
$[\delta_{CP},\ \alpha_{21}]$ | $[248^{\circ},\,80.3^{\circ}]$
$\sum m_{i}$ |  $58.5$ meV
$\langle m_{ee}\rangle$ |  $1.65$ meV
$\sqrt{\Delta\chi^{2}}$ | $1.40$
Table 2: A benchmark point of our input parameters and observables, where we
select it so that $\sqrt{\Delta\chi^{2}}$ is minimum.

##  IV Summary and discussion

We have studied a linear seesaw model with as minimum field content as
possible, introducing a modular $S_{4}$ with the help of $U(1)_{B-L}$
symmetries. Due to rank two neutrino mass matrix, we have had a vanishing
neutrino mass eigenvalue. Furthermore, only the normal mass hierarchy of
neutrinos is favored through the modular $S_{4}$ symmetry. In our numerical
$\Delta\chi^{2}$ analysis, we have found rather sharp prediction on sum of
neutrino masses to be around $60$ meV. Imaginary part of $\tau$ is restricted
at 1.35-1.5, while real part runs whole the range in the fundamental region.
The other remarks are listed below:

  1. 1.

$\alpha_{21}$ is within $\\{0^{\circ}-120^{\circ},240^{\circ}-310^{\circ}\\}$,
while $\delta_{CP}$ is within $\\{70^{\circ}-290^{\circ}\\}$.

  2. 2.

$\langle m_{ee}\rangle$ is allowed by $\\{0.001-0.004\\}$ eV.

Therefore our model indicates several predictions in neutrino sector that is
due to minimal structure with $S_{4}$ modular symmetry.

## Acknowledgments

This research was supported by an appointment to the JRG Program at the APCTP
through the Science and Technology Promotion Fund and Lottery Fund of the
Korean Government. This was also supported by the Korean Local Governments -
Gyeongsangbuk-do Province and Pohang City (H.O.). H. O. is sincerely grateful
for the KIAS member.

## Appendix

Here we review some properties of modular $S_{4}$ symmetry. In general, the
modular group $\bar{\Gamma}$ is the group of linear fractional transformation
$\gamma$ acting on the modulus $\tau$ which belongs to the upper-half complex
plane and transforms as

$\tau\longrightarrow\gamma\tau=\frac{a\tau+b}{c\tau+d}\ ,~{}~{}{\rm
where}~{}~{}a,b,c,d\in\mathbb{Z}~{}~{}{\rm and}~{}~{}ad-bc=1,~{}~{}{\rm
Im}[\tau]>0~{}.$ (IV.1)

This is isomorphic to $PSL(2,\mathbb{Z})=SL(2,\mathbb{Z})/\\{I,-I\\}$
transformation. Then modular transformation is generated by two
transformations $S$ and $T$ defined as follows;

$\displaystyle S:\tau\longrightarrow-\frac{1}{\tau}\ ,\qquad\qquad
T:\tau\longrightarrow\tau+1\ ,$ (IV.2)

and they satisfy the following algebraic relations,

$S^{2}=\mathbb{I}\ ,\qquad(ST)^{3}=\mathbb{I}\ .$ (IV.3)

Here we introduce the series of groups $\Gamma(N)~{}(N=1,2,3,\dots)$ which are
defined by

$\displaystyle\begin{aligned} \Gamma(N)=\left\\{\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in SL(2,\mathbb{Z})~{},~{}~{}\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}~{}~{}({\rm
mod}N)\right\\}\end{aligned},$ (IV.4)

and we define $\bar{\Gamma}(2)\equiv\Gamma(2)/\\{I,-I\\}$ for $N=2$. Since the
element $-I$ does not belong to $\Gamma(N)$ for $N>2$ case, we have
$\bar{\Gamma}(N)=\Gamma(N)$, that are infinite normal subgroup of
$\bar{\Gamma}$ known as principal congruence subgroups. We thus obtain finite
modular groups as the quotient groups defined by
$\Gamma_{N}\equiv\bar{\Gamma}/\bar{\Gamma}(N)$. For these finite groups
$\Gamma_{N}$, $T^{N}=\mathbb{I}$ is imposed, and the groups $\Gamma_{N}$ with
$N=2,3,4,5$ are isomorphic to $S_{3}$, $A_{4}$, $S_{4}$ and $A_{5}$,
respectively deAdelhartToorop:2011re .

Modular forms of level $N$ are holomorphic functions $f(\tau)$ which are
transformed under the action of $\Gamma(N)$ given by

$f(\gamma\tau)=(c\tau+d)^{k}f(\tau)~{},~{}~{}\gamma\in\Gamma(N)~{},$ (IV.5)

where $k$ is the so-called as the modular weight.

Here we discuss the modular symmetric theory framework without imposing
supersymmetry explicitly, considering the $S_{4}$ ($N=4$) modular group. Under
the modular transformation in Eq.(IV.1), a field $\phi^{(I)}$ is also
transformed as

$\phi^{(I)}\to(c\tau+d)^{-k_{I}}\rho^{(I)}(\gamma)\phi^{(I)},$ (IV.6)

where $-k_{I}$ is the modular weight and $\rho^{(I)}(\gamma)$ denotes an
unitary representation matrix of $\gamma\in\Gamma(4)$. Thus Lagrangian such as
Yukawa terms can be invariant if sum of modular weight from fields and modular
form in corresponding term is zero (also invariant under $S_{4}$ and gauge
symmetry).

The kinetic terms and quadratic terms of scalar fields can be written by

$\sum_{I}\frac{|\partial_{\mu}\phi^{(I)}|^{2}}{(-i\tau+i\bar{\tau})^{k_{I}}}~{},\quad\sum_{I}\frac{|\phi^{(I)}|^{2}}{(-i\tau+i\bar{\tau})^{k_{I}}}~{},$
(IV.7)

which is invariant under the modular transformation and overall factor is
eventually absorbed by a field redefinition consistently. Therefore the
Lagrangian associated with these terms should be invariant under the modular
symmetry.

The basis of modular forms with weight 2, $Y^{(2)}_{\rm 2}=(y_{1},y_{2})$ and
$Y^{(2)}_{\rm 3^{\prime}}=(y_{3},y_{4},y_{5})$, transforming as a doublet and
triplet under $S_{4}$, are found in terms of the Dedekind eta-function
$\eta(\tau)$ and its derivative Feruglio:2017spp :

$\displaystyle y_{1}(\tau)$ $\displaystyle=$
$\displaystyle\frac{i}{8}\left(8\frac{\eta^{\prime}(\tau+\frac{1}{2})}{\eta(\tau+\frac{1}{2})}+32\frac{\eta^{\prime}(4\tau)}{\eta(4\tau)}-\frac{\eta^{\prime}(\frac{\tau}{4})}{\eta(\frac{\tau}{4})}-\frac{\eta^{\prime}(\frac{\tau+1}{4})}{\eta(\frac{\tau+1}{4})}-\frac{\eta^{\prime}(\frac{\tau+2}{4})}{\eta(\frac{\tau+2}{4})}-\frac{\eta^{\prime}(\frac{\tau+3}{4})}{\eta(\frac{\tau+3}{4})}\right),$
$\displaystyle y_{2}(\tau)$ $\displaystyle=$
$\displaystyle\frac{i\sqrt{3}}{8}\left(\frac{\eta^{\prime}(\frac{\tau}{4})}{\eta(\frac{\tau}{4})}-\frac{\eta^{\prime}(\frac{\tau+1}{4})}{\eta(\frac{\tau+1}{4})}+\frac{\eta^{\prime}(\frac{\tau+2}{4})}{\eta(\frac{\tau+2}{4})}-\frac{\eta^{\prime}(\frac{\tau+3}{4})}{\eta(\frac{\tau+3}{4})}\right)$
$\displaystyle y_{3}(\tau)$ $\displaystyle=$ $\displaystyle
i\left(\frac{\eta^{\prime}(\tau+\frac{1}{2})}{\eta(\tau+\frac{1}{2})}-4\frac{\eta^{\prime}(4\tau)}{\eta(4\tau)}\right),$
$\displaystyle y_{4}(\tau)$ $\displaystyle=$
$\displaystyle\frac{i}{4\sqrt{2}}\left(-\frac{\eta^{\prime}(\frac{\tau}{4})}{\eta(\frac{\tau}{4})}+i\frac{\eta^{\prime}(\frac{\tau+1}{4})}{\eta(\frac{\tau+1}{4})}+\frac{\eta^{\prime}(\frac{\tau+2}{4})}{\eta(\frac{\tau+2}{4})}-i\frac{\eta^{\prime}(\frac{\tau+3}{4})}{\eta(\frac{\tau+3}{4})}\right),$
$\displaystyle y_{5}(\tau)$ $\displaystyle=$
$\displaystyle\frac{i}{4\sqrt{2}}\left(-\frac{\eta^{\prime}(\frac{\tau}{4})}{\eta(\frac{\tau}{4})}-i\frac{\eta^{\prime}(\frac{\tau+1}{4})}{\eta(\frac{\tau+1}{4})}+\frac{\eta^{\prime}(\frac{\tau+2}{4})}{\eta(\frac{\tau+2}{4})}+i\frac{\eta^{\prime}(\frac{\tau+3}{4})}{\eta(\frac{\tau+3}{4})}\right).$
(IV.8)

$y_{i}$’s can be expanded in terms of $q$ as follows:

$\displaystyle y_{1}$
$\displaystyle=-3\pi\left(\frac{b_{1}}{8}+3b_{5}\right),\ y_{2}=3\sqrt{3}\pi
b_{3},\ y_{3}=-\pi\left(-\frac{b_{1}}{4}+2b_{5}\right),$ $\displaystyle y_{4}$
$\displaystyle=-\pi\sqrt{2}b_{2},\ y_{5}=-4\pi\sqrt{2}b_{4},$ (IV.9)

where $b_{i}$ are given by

$\displaystyle b_{1}\sim 1,\ b_{2}\sim q,\ b_{3}\sim q^{2},\ b_{4}\sim 0,\
b_{5}\sim 0,$ (IV.10)

with $q=\exp(2\pi i\tau)$ and $|q|\ll 1$ Novichkov:2019sqv .

Then, Yukawas with higher weights are constructed by multiplication rules of
$S_{4}$, and one finds the following couplings:

$\displaystyle Y^{(4)}_{\bf 2}=\left[\begin{array}[]{c}y_{2}^{2}-y_{1}^{2}\\\
2y_{1}y_{2}\\\ \end{array}\right],\quad Y^{(4)}_{\bf
3}=\left[\begin{array}[]{c}-2y_{2}y_{3}\\\ \sqrt{3}y_{1}y_{5}+y_{2}y_{4}\\\
\sqrt{3}y_{1}y_{4}+y_{2}y_{5}\\\ \end{array}\right],\quad Y^{(4)}_{\bf
3^{\prime}}=\left[\begin{array}[]{c}2y_{1}y_{3}\\\
\sqrt{3}y_{2}y_{5}-y_{1}y_{4}\\\ \sqrt{3}y_{2}y_{4}-y_{1}y_{5}\\\
\end{array}\right].$ (IV.19)

## References

  * (1) P. Minkowski, Phys. Lett. B 67, 421 (1977);
  * (2) T. Yanagida, in _Proceedings of the Workshop on the Unified Theory and the Baryon Number in the Universe_ (O. Sawada and A. Sugamoto, eds.), KEK, Tsukuba, Japan, 1979, p. 95;
  * (3) M. Gell-Mann, P. Ramond, and R. Slansky, _Supergravity_ (P. van Nieuwenhuizen et al. eds.), North Holland, Amsterdam, 1979, p. 315; S. L. Glashow, _The future of elementary particle physics_ , in _Proceedings of the 1979 Cargèse Summer Institute on Quarks and Leptons_ (M. Levy et al. eds.), Plenum Press, New York, 1980, p. 687;
  * (4) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980).
  * (5) R. N. Mohapatra and J. W. F. Valle, Phys. Rev. D 34, 1642 (1986).
  * (6) D. Wyler and L. Wolfenstein, Nucl. Phys. B 218, 205 (1983).
  * (7) E. K. Akhmedov, M. Lindner, E. Schnapka and J. W. F. Valle, Phys. Lett. B 368, 270 (1996) [hep-ph/9507275].
  * (8) E. K. Akhmedov, M. Lindner, E. Schnapka and J. W. F. Valle, Phys. Rev. D 53, 2752 (1996) [hep-ph/9509255].
  * (9) M. Malinsky, J. C. Romao and J. W. F. Valle, Phys. Rev. Lett. 95 (2005), 161801 doi:10.1103/PhysRevLett.95.161801 [arXiv:hep-ph/0506296 [hep-ph]].
  * (10) F. Feruglio, doi:10.1142/9789813238053_0012 [arXiv:1706.08749 [hep-ph]].
  * (11) R. de Adelhart Toorop, F. Feruglio and C. Hagedorn, Nucl. Phys. B 858, 437-467 (2012) doi:10.1016/j.nuclphysb.2012.01.017 [arXiv:1112.1340 [hep-ph]].
  * (12) J. C. Criado and F. Feruglio, SciPost Phys. 5, no.5, 042 (2018) doi:10.21468/SciPostPhys.5.5.042 [arXiv:1807.01125 [hep-ph]].
  * (13) T. Kobayashi, N. Omoto, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, JHEP 11, 196 (2018) doi:10.1007/JHEP11(2018)196 [arXiv:1808.03012 [hep-ph]].
  * (14) H. Okada and M. Tanimoto, Phys. Lett. B 791, 54-61 (2019) doi:10.1016/j.physletb.2019.02.028 [arXiv:1812.09677 [hep-ph]].
  * (15) T. Nomura and H. Okada, Phys. Lett. B 797, 134799 (2019) doi:10.1016/j.physletb.2019.134799 [arXiv:1904.03937 [hep-ph]].
  * (16) H. Okada and M. Tanimoto, Eur. Phys. J. C 81, no.1, 52 (2021) doi:10.1140/epjc/s10052-021-08845-y [arXiv:1905.13421 [hep-ph]].
  * (17) F. J. de Anda, S. F. King and E. Perdomo, Phys. Rev. D 101, no.1, 015028 (2020) doi:10.1103/PhysRevD.101.015028 [arXiv:1812.05620 [hep-ph]].
  * (18) P. P. Novichkov, S. T. Petcov and M. Tanimoto, Phys. Lett. B 793, 247-258 (2019) doi:10.1016/j.physletb.2019.04.043 [arXiv:1812.11289 [hep-ph]].
  * (19) T. Nomura and H. Okada, Nucl. Phys. B 966, 115372 (2021) doi:10.1016/j.nuclphysb.2021.115372 [arXiv:1906.03927 [hep-ph]].
  * (20) H. Okada and Y. Orikasa, [arXiv:1907.13520 [hep-ph]].
  * (21) G. J. Ding, S. F. King and X. G. Liu, JHEP 09, 074 (2019) doi:10.1007/JHEP09(2019)074 [arXiv:1907.11714 [hep-ph]].
  * (22) T. Nomura, H. Okada and O. Popov, Phys. Lett. B 803, 135294 (2020) doi:10.1016/j.physletb.2020.135294 [arXiv:1908.07457 [hep-ph]].
  * (23) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, Phys. Rev. D 100, no.11, 115045 (2019) [erratum: Phys. Rev. D 101, no.3, 039904 (2020)] doi:10.1103/PhysRevD.100.115045 [arXiv:1909.05139 [hep-ph]].
  * (24) T. Asaka, Y. Heo, T. H. Tatsuishi and T. Yoshida, JHEP 01, 144 (2020) doi:10.1007/JHEP01(2020)144 [arXiv:1909.06520 [hep-ph]].
  * (25) D. Zhang, Nucl. Phys. B 952, 114935 (2020) doi:10.1016/j.nuclphysb.2020.114935 [arXiv:1910.07869 [hep-ph]].
  * (26) G. J. Ding, S. F. King, X. G. Liu and J. N. Lu, JHEP 12, 030 (2019) doi:10.1007/JHEP12(2019)030 [arXiv:1910.03460 [hep-ph]].
  * (27) T. Kobayashi, T. Nomura and T. Shimomura, Phys. Rev. D 102, no.3, 035019 (2020) doi:10.1103/PhysRevD.102.035019 [arXiv:1912.00637 [hep-ph]].
  * (28) T. Nomura, H. Okada and S. Patra, Nucl. Phys. B 967, 115395 (2021) doi:10.1016/j.nuclphysb.2021.115395 [arXiv:1912.00379 [hep-ph]].
  * (29) X. Wang, Nucl. Phys. B 957, 115105 (2020) doi:10.1016/j.nuclphysb.2020.115105 [arXiv:1912.13284 [hep-ph]].
  * (30) H. Okada and Y. Shoji, Nucl. Phys. B 961, 115216 (2020) doi:10.1016/j.nuclphysb.2020.115216 [arXiv:2003.13219 [hep-ph]].
  * (31) H. Okada and M. Tanimoto, [arXiv:2005.00775 [hep-ph]].
  * (32) M. K. Behera, S. Singirala, S. Mishra and R. Mohanta, [arXiv:2009.01806 [hep-ph]].
  * (33) M. K. Behera, S. Mishra, S. Singirala and R. Mohanta, [arXiv:2007.00545 [hep-ph]].
  * (34) T. Nomura and H. Okada, [arXiv:2007.04801 [hep-ph]].
  * (35) T. Nomura and H. Okada, [arXiv:2007.15459 [hep-ph]].
  * (36) T. Asaka, Y. Heo and T. Yoshida, Phys. Lett. B 811, 135956 (2020) doi:10.1016/j.physletb.2020.135956 [arXiv:2009.12120 [hep-ph]].
  * (37) H. Okada and M. Tanimoto, Phys. Rev. D 103, no.1, 015005 (2021) doi:10.1103/PhysRevD.103.015005 [arXiv:2009.14242 [hep-ph]].
  * (38) K. I. Nagao and H. Okada, [arXiv:2010.03348 [hep-ph]].
  * (39) H. Okada and M. Tanimoto, JHEP 03, 010 (2021) doi:10.1007/JHEP03(2021)010 [arXiv:2012.01688 [hep-ph]].
  * (40) C. Y. Yao, J. N. Lu and G. J. Ding, JHEP 05 (2021), 102 doi:10.1007/JHEP05(2021)102 [arXiv:2012.13390 [hep-ph]].
  * (41) P. Chen, G. J. Ding and S. F. King, JHEP 04 (2021), 239 doi:10.1007/JHEP04(2021)239 [arXiv:2101.12724 [hep-ph]].
  * (42) M. Kashav and S. Verma, [arXiv:2103.07207 [hep-ph]].
  * (43) H. Okada, Y. Shimizu, M. Tanimoto and T. Yoshida, [arXiv:2105.14292 [hep-ph]].
  * (44) I. de Medeiros Varzielas and J. Lourenço, [arXiv:2107.04042 [hep-ph]].
  * (45) T. Nomura, H. Okada and Y. Orikasa, [arXiv:2106.12375 [hep-ph]].
  * (46) P. T. P. Hutauruk, D. W. Kang, J. Kim and H. Okada, [arXiv:2012.11156 [hep-ph]].
  * (47) G. J. Ding, S. F. King and J. N. Lu, [arXiv:2108.09655 [hep-ph]].
  * (48) K. I. Nagao and H. Okada, [arXiv:2108.09984 [hep-ph]].
  * (49) T. Kobayashi, K. Tanaka and T. H. Tatsuishi, Phys. Rev. D 98, no.1, 016004 (2018) doi:10.1103/PhysRevD.98.016004 [arXiv:1803.10391 [hep-ph]].
  * (50) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto, T. H. Tatsuishi and H. Uchida, Phys. Lett. B 794, 114-121 (2019) doi:10.1016/j.physletb.2019.05.034 [arXiv:1812.11072 [hep-ph]].
  * (51) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, PTEP 2020, no.5, 053B05 (2020) doi:10.1093/ptep/ptaa055 [arXiv:1906.10341 [hep-ph]].
  * (52) H. Okada and Y. Orikasa, Phys. Rev. D 100, no.11, 115037 (2019) doi:10.1103/PhysRevD.100.115037 [arXiv:1907.04716 [hep-ph]].
  * (53) S. Mishra, [arXiv:2008.02095 [hep-ph]].
  * (54) X. Du and F. Wang, JHEP 02, 221 (2021) doi:10.1007/JHEP02(2021)221 [arXiv:2012.01397 [hep-ph]].
  * (55) J. T. Penedo and S. T. Petcov, Nucl. Phys. B 939, 292-307 (2019) doi:10.1016/j.nuclphysb.2018.12.016 [arXiv:1806.11040 [hep-ph]].
  * (56) P. P. Novichkov, J. T. Penedo, S. T. Petcov and A. V. Titov, JHEP 04, 005 (2019) doi:10.1007/JHEP04(2019)005 [arXiv:1811.04933 [hep-ph]].
  * (57) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, JHEP 02, 097 (2020) doi:10.1007/JHEP02(2020)097 [arXiv:1907.09141 [hep-ph]].
  * (58) S. F. King and Y. L. Zhou, Phys. Rev. D 101, no.1, 015001 (2020) doi:10.1103/PhysRevD.101.015001 [arXiv:1908.02770 [hep-ph]].
  * (59) H. Okada and Y. Orikasa, [arXiv:1908.08409 [hep-ph]].
  * (60) J. C. Criado, F. Feruglio and S. J. D. King, JHEP 02, 001 (2020) doi:10.1007/JHEP02(2020)001 [arXiv:1908.11867 [hep-ph]].
  * (61) X. Wang and S. Zhou, JHEP 05, 017 (2020) doi:10.1007/JHEP05(2020)017 [arXiv:1910.09473 [hep-ph]].
  * (62) Y. Zhao and H. H. Zhang, JHEP 03 (2021), 002 doi:10.1007/JHEP03(2021)002 [arXiv:2101.02266 [hep-ph]].
  * (63) S. F. King and Y. L. Zhou, JHEP 04 (2021), 291 doi:10.1007/JHEP04(2021)291 [arXiv:2103.02633 [hep-ph]].
  * (64) G. J. Ding, S. F. King and C. Y. Yao, [arXiv:2103.16311 [hep-ph]].
  * (65) X. Zhang and S. Zhou, [arXiv:2106.03433 [hep-ph]].
  * (66) Bu-Yao Qu, Xiang-Gan Liu, Ping-Tao Chen, Gui-Jun Ding [arXiv:2106.11659 [hep-ph]].
  * (67) P. P. Novichkov, J. T. Penedo, S. T. Petcov and A. V. Titov, JHEP 04, 174 (2019) doi:10.1007/JHEP04(2019)174 [arXiv:1812.02158 [hep-ph]].
  * (68) G. J. Ding, S. F. King and X. G. Liu, Phys. Rev. D 100, no.11, 115005 (2019) doi:10.1103/PhysRevD.100.115005 [arXiv:1903.12588 [hep-ph]].
  * (69) X. Wang, B. Yu and S. Zhou, Phys. Rev. D 103, no.7, 076005 (2021) doi:10.1103/PhysRevD.103.076005 [arXiv:2010.10159 [hep-ph]].
  * (70) C. Y. Yao, X. G. Liu and G. J. Ding, Phys. Rev. D 103, no.9, 095013 (2021) doi:10.1103/PhysRevD.103.095013 [arXiv:2011.03501 [hep-ph]].
  * (71) X. Wang and S. Zhou, [arXiv:2102.04358 [hep-ph]].
  * (72) M. K. Behera and R. Mohanta, [arXiv:2108.01059 [hep-ph]].
  * (73) A. Baur, H. P. Nilles, A. Trautner and P. K. S. Vaudrevange, Phys. Lett. B 795, 7-14 (2019) doi:10.1016/j.physletb.2019.03.066 [arXiv:1901.03251 [hep-th]].
  * (74) I. de Medeiros Varzielas, S. F. King and Y. L. Zhou, Phys. Rev. D 101, no.5, 055033 (2020) doi:10.1103/PhysRevD.101.055033 [arXiv:1906.02208 [hep-ph]].
  * (75) X. G. Liu and G. J. Ding, JHEP 08, 134 (2019) doi:10.1007/JHEP08(2019)134 [arXiv:1907.01488 [hep-ph]].
  * (76) P. Chen, G. J. Ding, J. N. Lu and J. W. F. Valle, Phys. Rev. D 102, no.9, 095014 (2020) doi:10.1103/PhysRevD.102.095014 [arXiv:2003.02734 [hep-ph]].
  * (77) C. C. Li, X. G. Liu and G. J. Ding, [arXiv:2108.02181 [hep-ph]].
  * (78) P. P. Novichkov, J. T. Penedo and S. T. Petcov, Nucl. Phys. B 963, 115301 (2021) doi:10.1016/j.nuclphysb.2020.115301 [arXiv:2006.03058 [hep-ph]].
  * (79) X. G. Liu, C. Y. Yao and G. J. Ding, Phys. Rev. D 103, no.5, 056013 (2021) doi:10.1103/PhysRevD.103.056013 [arXiv:2006.10722 [hep-ph]].
  * (80) S. Kikuchi, T. Kobayashi, H. Otsuka, S. Takada and H. Uchida, JHEP 11, 101 (2020) doi:10.1007/JHEP11(2020)101 [arXiv:2007.06188 [hep-th]].
  * (81) Y. Almumin, M. C. Chen, V. Knapp-Pérez, S. Ramos-Sánchez, M. Ratz and S. Shukla, JHEP 05 (2021), 078 doi:10.1007/JHEP05(2021)078 [arXiv:2102.11286 [hep-th]].
  * (82) G. J. Ding, F. Feruglio and X. G. Liu, SciPost Phys. 10 (2021), 133 doi:10.21468/SciPostPhys.10.6.133 [arXiv:2102.06716 [hep-ph]].
  * (83) F. Feruglio, V. Gherardi, A. Romanino and A. Titov, JHEP 05 (2021), 242 doi:10.1007/JHEP05(2021)242 [arXiv:2101.08718 [hep-ph]].
  * (84) S. Kikuchi, T. Kobayashi and H. Uchida, [arXiv:2101.00826 [hep-th]].
  * (85) P. P. Novichkov, J. T. Penedo and S. T. Petcov, JHEP 04 (2021), 206 doi:10.1007/JHEP04(2021)206 [arXiv:2102.07488 [hep-ph]].
  * (86) G. Altarelli and F. Feruglio, Rev. Mod. Phys. 82, 2701-2729 (2010) doi:10.1103/RevModPhys.82.2701 [arXiv:1002.0211 [hep-ph]].
  * (87) H. Ishimori, T. Kobayashi, H. Ohki, Y. Shimizu, H. Okada and M. Tanimoto, Prog. Theor. Phys. Suppl. 183, 1-163 (2010) doi:10.1143/PTPS.183.1 [arXiv:1003.3552 [hep-th]].
  * (88) H. Ishimori, T. Kobayashi, H. Ohki, H. Okada, Y. Shimizu and M. Tanimoto, Lect. Notes Phys. 858, 1-227 (2012) doi:10.1007/978-3-642-30805-5
  * (89) D. Hernandez and A. Y. Smirnov, Phys. Rev. D 86, 053014 (2012) doi:10.1103/PhysRevD.86.053014 [arXiv:1204.0445 [hep-ph]].
  * (90) S. F. King and C. Luhn, Rept. Prog. Phys. 76, 056201 (2013) doi:10.1088/0034-4885/76/5/056201 [arXiv:1301.1340 [hep-ph]].
  * (91) S. F. King, A. Merle, S. Morisi, Y. Shimizu and M. Tanimoto, New J. Phys. 16, 045018 (2014) doi:10.1088/1367-2630/16/4/045018 [arXiv:1402.4271 [hep-ph]].
  * (92) S. F. King, Prog. Part. Nucl. Phys. 94, 217-256 (2017) doi:10.1016/j.ppnp.2017.01.003 [arXiv:1701.04413 [hep-ph]].
  * (93) S. T. Petcov, Eur. Phys. J. C 78, no.9, 709 (2018) doi:10.1140/epjc/s10052-018-6158-5 [arXiv:1711.10806 [hep-ph]].
  * (94) A. Baur, H. P. Nilles, A. Trautner and P. K. S. Vaudrevange, Nucl. Phys. B 947, 114737 (2019) doi:10.1016/j.nuclphysb.2019.114737 [arXiv:1908.00805 [hep-th]].
  * (95) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto, T. H. Tatsuishi and H. Uchida, Phys. Rev. D 101, no.5, 055046 (2020) doi:10.1103/PhysRevD.101.055046 [arXiv:1910.11553 [hep-ph]].
  * (96) P. P. Novichkov, J. T. Penedo, S. T. Petcov and A. V. Titov, JHEP 07, 165 (2019) doi:10.1007/JHEP07(2019)165 [arXiv:1905.11970 [hep-ph]].
  * (97) M. Tanimoto and K. Yamamoto, [arXiv:2106.10919 [hep-ph]].
  * (98) M. C. Chen, S. Ramos-Sánchez and M. Ratz, Phys. Lett. B 801, 135153 (2020) doi:10.1016/j.physletb.2019.135153 [arXiv:1909.06910 [hep-ph]].
  * (99) I. de Medeiros Varzielas, M. Levy and Y. L. Zhou, JHEP 11, 085 (2020) doi:10.1007/JHEP11(2020)085 [arXiv:2008.05329 [hep-ph]].
  * (100) A. Das, T. Nomura, H. Okada and S. Roy, Phys. Rev. D 96, no. 7, 075001 (2017) [arXiv:1704.02078 [hep-ph]].
  * (101) W. Wang and Z. L. Han, Phys. Rev. D 92 (2015), 095001 [arXiv:1508.00706 [hep-ph]].
  * (102) A. Gando et al. [KamLAND-Zen Collaboration], Phys. Rev. Lett. 117, no. 8, 082503 (2016) Addendum: [Phys. Rev. Lett. 117, no. 10, 109903 (2016)] [arXiv:1605.02889 [hep-ex]].
  * (103) I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, JHEP 1901, 106 (2019) [arXiv:1811.05487 [hep-ph]].
  * (104) I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni, and T. Schwetz, NuFIT 4.1 (2019), www.nu-fit.org, (2019).
  * (105) E. Fernandez-Martinez, J. Hernandez-Garcia and J. Lopez-Pavon, JHEP 08, 033 (2016) doi:10.1007/JHEP08(2016)033 [arXiv:1605.08774 [hep-ph]].
  * (106) N. R. Agostinho, G. C. Branco, P. M. F. Pereira, M. N. Rebelo and J. I. Silva-Marcos, Eur. Phys. J. C 78, no.11, 895 (2018) doi:10.1140/epjc/s10052-018-6347-2 [arXiv:1711.06229 [hep-ph]].