aainstitutetext: Department of Physics, Hokkaido University, Sapporo 060-0810,
Japanbbinstitutetext: Waseda Institute for Advanced Study, Waseda University,
1-21-1 Nishi Waseda, Shinjuku, Tokyo 169-0051, Japan

# Moduli trapping mechanism in modular flavor symmetric models

Shota Kikuchi a Tatsuo Kobayashi a Kaito Nasu b Yusuke Yamada

###### Abstract

We discuss how the moduli in modular flavor symmetric models dynamically
select enhanced symmetry points at which the residual modular symmetry renders
extra matter fields massless. The moduli dynamics non-perturbatively produces
the extra matter particles, which gives (time-dependent) effective potential
that traps the moduli to enhanced symmetry points. We show analytic estimates
of particle production rate consistent with numerical results, and the
dynamics of moduli based on the analytic estimates.

††preprint: EPHOU-23-014

##  1 Introduction

In effective theory of superstring, moduli fields, light fields associated
e.g. with higher-dimensional gravitational degrees of freedom, play crucial
roles in constructing realistic models of particle physics as well as
cosmology.

Modular symmetry is the geometrical symmetry of compact space, where moduli
transform non-trivially. Four-dimensional effective field theory derived from
superstring theory also has the modular symmetry Ferrara:1989bc ;
Ferrara:1989qb ; Lerche:1989cs ; Lauer:1989ax ; Lauer:1990tm ;
Kobayashi:2018rad ; Kobayashi:2018bff ; Ohki:2020bpo ; Kikuchi:2020frp ;
Kikuchi:2020nxn ; Kikuchi:2021ogn ; Almumin:2021fbk ; Baur:2019kwi ;
Baur:2019iai ; Nilles:2020kgo ; Nilles:2020gvu . In addition, the modular
symmetry includes finite groups such as $S_{3}$, $A_{4}$, $S_{4}$, and $A_{5}$
deAdelhartToorop:2011re . These discrete groups have been used in flavor model
building in the bottom-up approach Altarelli:2010gt ; Ishimori:2010au ;
Kobayashi:2022moq ; Hernandez:2012ra ; King:2013eh ; King:2014nza ;
Tanimoto:2015nfa ; King:2017guk ; Petcov:2017ggy . Inspired by these aspects,
modular flavor symmetric models have been receiving attention as an origin of
flavor structure in the standard model of particle physics Feruglio:2017spp ;
Kobayashi:2018vbk ; Penedo:2018nmg ; Criado:2018thu ; Kobayashi:2018scp ;
Novichkov:2018ovf ; Novichkov:2018nkm ; deAnda:2018ecu ; Kobayashi:2018wkl ;
Okada:2018yrn ; Novichkov:2018yse ; Ding:2019xna ; Nomura:2019jxj ;
Novichkov:2019sqv ; Okada:2019uoy ; deMedeirosVarzielas:2019cyj ;
Nomura:2019yft ; Kobayashi:2019rzp ; Liu:2019khw ; Okada:2019xqk ;
Kobayashi:2019mna ; Ding:2019zxk ; King:2019vhv ; Nomura:2019lnr ;
Criado:2019tzk ; Kobayashi:2019xvz ; Asaka:2019vev ; Ding:2019gof ;
Zhang:2019ngf ; Wang:2019ovr ; Kobayashi:2019uyt ; Nomura:2019xsb ;
Kobayashi:2019gtp ; Lu:2019vgm ; Wang:2019xbo ; King:2020qaj ; Abbas:2020qzc ;
Chen:2020udk ; Okada:2020oxh ; Okada:2020dmb ; Ding:2020yen ; Ding:2020msi ;
Okada:2020rjb ; Novichkov:2020eep ; Liu:2020akv ; deMedeirosVarzielas:2020kji
; Asaka:2020tmo ; Okada:2020ukr ; Wang:2020lxk ; Okada:2020brs ; Yao:2020qyy ;
Li:2021buv 111See for more references Ref. Kobayashi:2023zzc ..

In the modular flavor symmetric models, Yukawa couplings correspond to modular
forms of complex structure moduli $\tau$ , and the structure of Yukawa
couplings is determined by the vacuum expectation values (VEVs) of complex
structure moduli, which are complex scalar fields in 4D effective theory.
Therefore, moduli stabilization plays a crucial role in such flavor symmetry
models. Indeed, moduli stabilization has been studied in modular flavor
symmetric models Kobayashi:2019xvz ; Kobayashi:2019uyt ; Abe:2020vmv ;
Ishiguro:2020tmo ; Novichkov:2022wvg ; Ishiguro:2022pde ; Leedom:2022zdm ;
Knapp-Perez:2023nty .

In particular, residual $Z_{N}$ symmetries remain at fixed points
Novichkov:2018ovf . Each modular form behaves like $(\tau-\tau_{*})^{m}$
around a fixed point $\tau_{*}$, depending on the $Z_{N}$ charge $m$ of the
modular form. Thus, the modular forms are suppressed by
$\varepsilon=\tau-\tau_{*}$ at nearby fixed point $\tau_{*}$. This behavior is
important to realize the hierarchy of quark and lepton masses and mixing
angles Feruglio:2021dte ; Novichkov:2021evw ; Petcov:2022fjf ; Kikuchi:2023cap
; Abe:2023ilq ; Kikuchi:2023jap ; Abe:2023qmr ; Petcov:2023vws ; Abe:2023dvr .

One of the ways to stabilize the moduli is the use of fluxes of $p$-form
fields along compactified dimensions that yield effective potential to complex
structure moduli. If somehow the chosen flux potential yields preferable VEV
to the complex structure moduli, it would explain the flavor structure in the
standard model by chance. Another possibility is due to non-perturbative
effects. On the other hand, there might be some dynamical origin why the
complex structure moduli may take some special value in the moduli space: In
Kofman:2004yc , it was pointed out that particle production takes place when
moduli dynamically crosses the special point called the enhanced symmetry
point (ESP) at which particles coupled to the moduli become massless. The
produced particle gives effective potential to the moduli such that the moduli
are attracted to the ESP and trapped around it. Such mechanism may explain why
the VEV of moduli chooses some particular value from the landscape of moduli
space.

In this work, we consider the dynamical moduli trapping mechanism in modular
flavor symmetric models. First, with simple toy models we briefly review and
develop a numerical approach to discuss the moduli dynamics back-reacted by
the production of spectator fields. Despite its straightforwardness, we find a
common difficulty due to the computational costs in the simulations.
Furthermore, the modular flavor symmetric models contain extra complications
due to the structure of Yukawa couplings as well as moduli being coupled to
matters via Planck suppressed operators. Therefore, we will develop a semi-
analytic approach where only the first particle production event is taken into
account, and contributions from all the momentum modes are included. Using the
analytic result, we show that the moduli trapping mechanism works even for
(complicated) modular flavor symmetric models.

This paper is organized as follows. In Sec. 2, we review and discuss a fully
numerical approach, which is in principle applicable to any dynamical systems.
In Sec. 3, we review the modular flavor symmetry , and then specify the system
we consider in this work and derive the equations of motion that we need to
solve. We then analytically evaluate the amount of particle production when
moduli crosses the ESP in Sec. 4. We also numerically examine the moduli
dynamics using the analytic expressions we derived. Finally, we summarize our
results and discuss implications to more realistic models in Sec. 5. In
Appendix A, we give a brief review on particle productions in time-dependent
background. Modular forms, which we use, are shown in Appendix B. In Appendix
C, we review field dynamics in expanding Universe.

Throughout this paper, we take a natural unit convention $\hbar=1$, $c=1$ but
we write the Planck scale $M_{\rm pl}\sim 2.4\times 10^{18}$GeV explicitly to
clarify the hierarchy between various scales.

##  2 A toy model

We briefly review some ingredients that we will use in modular flavor
symmetric models. For simplicity, we discuss the following system consists of
two massless real scalar fields,

$\displaystyle S=-\frac{1}{2}\int
d^{4}x\sqrt{-g}\left[(\partial\phi)^{2}+(\partial\chi)^{2}+\lambda\phi^{2}\chi^{2}\right].$
(1)

We assume that $\phi$ behaves as a classical homogeneous field $\phi=\phi(t)$,
the background to be the Friedman-Robertson-Walker (FRW) spacetime
$ds^{2}=-dt^{2}+a^{2}(t)d{\bm{x}}^{2}$ and $\chi$ to be a quantum field, and
then the coupling term $\lambda\phi(t)^{2}\chi^{2}$ behaves as an effective
time-dependent mass term of $\chi$. We quantize $\chi$ as follows: First, the
action of $\chi$ is

$S=\frac{1}{2}\int
dtd^{3}xa^{3}\left[\dot{\chi}^{2}-a^{-2}(\partial_{i}\chi)^{2}-\lambda\phi^{2}\chi^{2}\right],$
(2)

which, by introducing $\tilde{\chi}=a^{\frac{3}{2}}\chi$, can be rewritten as

$\displaystyle S=$ $\displaystyle\frac{1}{2}\int
dtd^{3}x\left[\left(\dot{\tilde{\chi}}-\frac{3}{2}H\tilde{\chi}\right)^{2}-a^{-2}(\partial_{i}\tilde{\chi})^{2}-\lambda\phi^{2}\tilde{\chi}^{2}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{2}\int
dtd^{3}x\left[\dot{\tilde{\chi}}^{2}-a^{-2}(\partial_{i}\tilde{\chi})^{2}-\left(\lambda\phi^{2}-\frac{3}{2}\dot{H}-\frac{9}{4}H^{2}\right)\tilde{\chi}^{2}\right].$
(3)

Then, the quantum field $\hat{\chi}$ is written by

$\hat{\chi}(t,{\bm{x}})=\frac{1}{a^{\frac{3}{2}}(t)}\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\left[\hat{a}_{\bm{k}}e^{{\rm
i}{\bm{k}}\cdot{\bm{x}}}f_{k}(t)+\hat{a}^{\dagger}_{\bm{k}}e^{-{\rm
i}{\bm{k}}\cdot{\bm{x}}}f^{*}_{k}(t)\right],$ (4)

where the mode function $f_{k}(t)$ satisfies

$\ddot{f}_{k}(t)+\omega_{k}^{2}(t)f_{k}(t)=0,$ (5)

and

$\omega_{k}^{2}(t)=\frac{k^{2}}{a^{2}(t)}+\lambda\phi^{2}(t)-\frac{3}{2}\dot{H}(t)-\frac{9}{4}H^{2}(t).$
(6)

The mode function $f_{k}(t)$ satisfies the normalization condition

$f_{k}\dot{f}_{k}^{*}(t)-f^{*}_{k}(t)\dot{f}_{k}(t)={\rm i},$ (7)

which implies the canonical commutation relation

$[\hat{a}_{\bm{k}},\hat{a}^{\dagger}_{{\bm{k}}^{\prime}}]=\delta^{3}({\bm{k}}-{\bm{k}}^{\prime})\
\text{or equivalently}\
[\hat{\chi}(t,{\bm{x}}),\dot{\hat{\chi}}(t,{\bm{y}})]={\rm
i}\delta^{3}({\bm{x}}-{\bm{y}}),$ (8)

for all $t$. In general, it is impossible to solve (5) analytically, but we
are able to introduce a formal adiabatic solution

$f_{k}(t)=\frac{1}{\sqrt{2\omega_{k}(t)}}\left(\alpha_{k}(t)e^{-{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}+\beta_{k}(t)e^{{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}\right).$ (9)

The auxiliary functions $\alpha_{k}(t),\beta_{k}(t)$ satisfy

$\displaystyle\dot{\alpha}_{k}(t)=$
$\displaystyle\frac{\dot{\omega}_{k}}{2\omega_{k}}\beta_{k}e^{+2{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}},$ (10)
$\displaystyle\dot{\beta}_{k}(t)=$
$\displaystyle\frac{\dot{\omega}_{k}}{2\omega_{k}}\alpha_{k}e^{-2{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}},$ (11)

and $|\alpha_{k}(t)|^{2}-|\beta_{k}(t)|^{2}=1$ for all $t$, which is
equivalent to (7). Note that the choice of the adiabatic solution is not
unique, and we have chosen the zeroth order adiabatic solution with
$V_{k}(t)=0$ in Habib:1999cs ; Dabrowski:2014ica ; Dabrowski:2016tsx . We
choose the initial condition $\alpha_{k}(t)\to 1$ and $\beta_{k}(t)\to 0$ as
$t\to-\infty$ which realizes the adiabatic past vacuum.222The infinite past
may be replaced by some finite $t$, which does not change as long as the
initial time is sufficiently far from the time of the first particle
production event. In general, $\beta_{k}(t)$ becomes non-vanishing due to the
time-dependence of background fields, which can be physically understood as
“particle production from vacuum”. We briefly review it in Appendix A. The
production of particles back-reacts to the dynamics of the background field.
Throughout this work, we assume that the backreaction affects only the modulus
$\phi$, and the background spacetime is intact.

Let us consider the dynamics of $\phi$ back-reacted by $\chi$-particles, whose
E.O.M. is given by

$\ddot{\phi}(t)+3H\dot{\phi}+\lambda\langle\chi^{2}\rangle_{\rm
ren}\phi(t)=0,$ (12)

where $\langle\hat{\chi}^{2}(x)\rangle_{\rm ren}$ is an expectation value of a
renormalized $\hat{\chi}^{2}$ operator with the adiabatic vacuum state. The
vacuum expectation value of $\hat{\chi}^{2}$ without renormalization is
explicitly given by

$\displaystyle\langle\chi^{2}\rangle=$ $\displaystyle\lim_{y\to x}\langle
0|_{\rm in}\hat{\chi}(x)\hat{\chi}(y)|0\rangle_{\rm in}$ $\displaystyle=$
$\displaystyle\lim_{y\to
x}\int\frac{d^{3}kd^{3}k^{\prime}}{(2\pi)^{3}a^{3}(t)}\delta^{3}({\bf k}-{\bf
k}^{\prime})\left(e^{{\rm i}{\bf k}\cdot{\bf x}-{\rm i}{\bf
k}^{\prime}\cdot{\bf y}}f_{k}(x^{0})f_{k^{\prime}}^{*}(y^{0})\right)$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}a^{3}(t)}|f_{k}(t)|^{2}$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}a^{3}(t)}\frac{1}{2\omega_{k}(t)}\left(|\alpha_{k}(t)|^{2}+|\beta_{k}(t)|^{2}+\alpha_{k}(t)\beta_{k}^{*}e^{2{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}+\alpha_{k}^{*}(t)\beta_{k}e^{-2{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}\right)$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}a^{3}(t)}\frac{1}{2\omega_{k}(t)}\left(1+2|\beta_{k}(t)|^{2}+\alpha_{k}(t)\beta_{k}^{*}e^{2{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}+\alpha_{k}^{*}(t)\beta_{k}e^{-2{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}\right).$ (13)

We note that $\beta_{k}(t)$ is vanishing unless particle production occurs,
which cannot be removed by local counter terms. Therefore, we may identify the
$\beta_{k}$-independent term to be a “vacuum” contribution. We first evaluate
the “vacuum” contribution as

$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}a^{3}(t)}\frac{1}{2\omega_{k}(t)}\to$
$\displaystyle\int\frac{d^{d}k}{(2\pi)^{d}a^{d}(t)}\frac{\mu^{3-d}}{2\omega_{k}(t)}$
$\displaystyle=$
$\displaystyle\frac{2\pi^{\frac{d}{2}}}{(2\pi)^{d}\Gamma\left(\frac{d}{2}\right)}\int
dk\frac{\mu^{3-d}k^{d-1}}{2\sqrt{k^{2}+M_{\rm eff}^{2}}}$ $\displaystyle=$
$\displaystyle\frac{M_{\rm eff}^{2}}{8\pi^{2}\epsilon}-\frac{M_{\rm
eff}^{2}}{16\pi^{2}}\left[1-\gamma_{E}-\log\left(\frac{M_{\rm
eff}^{2}}{4\pi\mu^{2}}\right)\right]+\mathcal{O}(\epsilon),$ (14)

where we have used the dimensional regularization $d=3+\epsilon$ with
$\epsilon\to 0$, $M_{\rm
eff}^{2}\equiv\lambda\phi^{2}-\frac{3}{2}\dot{H}-\frac{9}{4}H^{2}$ is the
effective mass,333The effective mass $M_{\rm eff}^{2}$ may become tachyonic,
which leads to imaginary part to the effective potential. Such instability
simply shows the existence of tachyonic modes that cannot be naively
integrated out. Therefore, one has to treat such modes separately if exist. In
the following, we consider only the cases where such instability shows up. and
$\mu$ is renormalization scale parameter. We here applied dimensional
regularization such that the invariance under trivial scale transformation
$a\to ca$, ${\bm{k}}\to c{\bm{k}}$ ($\bm{x}\to c^{-1}\bm{x}$) with a constant
$c$ holds. Thus, with appropriate counter terms, we can formally write the
E.O.M. of $\phi$ as

$\ddot{\phi}(t)+\delta V(\phi)+F(t)=0,$ (15)

where $\delta V(\phi)$ denotes quantum corrected potential terms of $\phi$ and

$\displaystyle F(t)\equiv$
$\displaystyle\lambda\phi(t)\int\frac{d^{3}k}{(2\pi)^{3}a^{3}(t)}\frac{1}{2\omega_{k}(t)}\left(2|\beta_{k}(t)|^{2}+\alpha_{k}(t)\beta_{k}^{*}e^{2{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}+\alpha_{k}^{*}(t)\beta_{k}e^{-2{\rm
i}\int^{t}\omega_{k}(t^{\prime})dt^{\prime}}\right)$ $\displaystyle\approx$
$\displaystyle\lambda\phi(t)\int\frac{d^{3}k}{(2\pi)^{3}a^{3}(t)}\frac{1}{\omega_{k}(t)}n_{k}(t),$
(16)

where $n_{k}(t)\equiv|\beta_{k}(t)|^{2}$ and we have dropped fast-oscillatory
terms since it would be averaged to be zero.444The neglected terms contain
memory of the past, which makes the E.O.M. of $\phi$ an integro-differential
equation. We note that $|\beta_{k}(t)|$ typically decays faster than any
powers of $k$, and therefore, there would be no UV divergence associated with
it.555For adiabatic vacuum states $\beta(-\infty)=0$, $\beta_{k}(t)$ becomes
non-zero due to Stokes phenomena/particle production, which is not a local but
a global property with respect to time. Therefore, it is reasonable that the
terms associated with such global (non-local) effects do not lead to UV
divergences which are local effects due to short wave length modes.

We would like to give some comments about solving the dynamics of the system
under consideration:

  1. 1.

In the following, we assume that the quantum-corrected effective potential
$\delta V(\phi)=0$ throughout this work. Such a situation is effectively
realized e.g. if the theory is supersymmetric. Strictly speaking, the time-
dependence of backgrounds break supersymmetry spontaneously and there would be
some potential. We will leave the effect of the quantum corrected potential
for future work.

  2. 2.

We use Emarkov-Milne equation, which directly yields the time-dependent
particle number density: We use $f_{k}(t)=\xi_{k}(t)e^{-{\rm
i}\lambda_{k}(t)}$ as an Ansatz for the mode function. From the normalization
condition (7) the function $\lambda_{k}(t)$ satisfies

$\lambda_{k}(t)=\frac{1}{2}\int^{t}dt^{\prime}\xi_{k}^{-2}(t^{\prime})$ (17)

and the mode equation reads

$\ddot{\xi}_{k}(t)+\omega_{k}^{2}(t)\xi_{k}(t)=\frac{1}{4}\xi^{-3}_{k}(t),$
(18)

with the initial values

$\xi_{k}(t_{0})\to\frac{1}{\sqrt{2\omega_{k}(t)}},\quad\dot{\xi}_{k}(t_{0})\to
0.$ (19)

In terms of $\xi_{k}(t)$, the particle number density $n_{k}(t)$ is simply
given by

$n_{k}(t)=\frac{\xi^{2}_{k}(t)}{2\omega_{k}(t)}\left(\left(\frac{1}{2\xi_{k}^{2}(t)}-\omega_{k}(t)\right)^{2}+\frac{\dot{\xi}_{k}^{2}(t)}{\xi_{k}^{2}(t)}\right),$
(20)

which is much easier to evaluate, since it can be solved simply as a set of
differential equations. It is crucial to neglect the oscillatory terms in
$F(t)$ to avoid an integro-differential equation.

  3. 3.

As done in Kofman:2004yc  we discretize the momentum integration appearing in
$F(t)$, which allows us to perform numerical simulations. However, even with
such approximation, it is still a hard problem unless we reduce the number of
$k$-modes in numerical simulations. To do so, it would be useful to rescale
the momentum $k$ by some reference scale. Assuming that the effective mass
$M_{\rm eff}^{2}$ is dominated by the modulus coupling $\lambda\phi^{2}$, the
most relevant scale turns out to be
$v=\lambda^{\frac{1}{4}}|\dot{\phi}(t_{0})|^{\frac{1}{2}}$ where $t_{0}$ is
the time when $\phi$ first crosses $\phi=0$ at which $M_{\rm
eff}\approx\lambda^{\frac{1}{2}}|\phi|=0$.666We will explain the reason why
$v$ is a reference scale below. Therefore, we rescale dimensionful quantities
by $v$ like $\tilde{k}=k/v$,

$\tilde{\omega}_{k}=\omega_{k}/v=\sqrt{\tilde{k}^{2}+\frac{M_{\rm
eff}^{2}(t)}{v^{2}}}.$ (21)

Noting that

$n_{k}(t)=\frac{\tilde{\omega}_{k}\tilde{\xi}^{2}_{k}}{2}\left(\left(\frac{1}{2\tilde{\omega}_{k}\tilde{\xi}_{k}^{2}}-1\right)^{2}+\frac{\dot{\tilde{\xi}}^{2}}{\tilde{\omega}^{2}_{k}\tilde{\xi}_{k}^{2}}\right),$
(22)

where $\tilde{\xi}_{k}=\sqrt{v}\xi_{k}$, we rewrite $F(t)$ as

$F(t)=\lambda\phi(t)v^{2}\int\frac{d^{3}\tilde{k}}{(2\pi)^{3}}\frac{n_{k}(t)}{\tilde{\omega}_{k}(t)}=\frac{\lambda\phi(t)v^{2}}{2\pi^{2}}\int
d\tilde{k}\frac{\tilde{k}^{2}n_{k}(t)}{\tilde{\omega}_{k}(t)}.$ (23)

Let us discuss how to perform the discretization of the integral effectively.
We note that the particle number density after the first crossing of
$m^{2}(t)=0$ is approximately given by

$n_{k}=\exp\left(-\pi\frac{k^{2}}{\sqrt{\lambda}|\dot{\phi}(t_{0})|}\right)=\exp\left(-\pi\tilde{k}^{2}\right).$
(24)

This is why we have taken $v$ as the reference scale. The mode
$\pi\tilde{k}^{2}\sim\mathcal{O}(10)$ does not contribute to the integration
and here we take $\tilde{k}=N$ to be the effective cut-off of the momentum
integration where $N=\mathcal{O}(1)$ is a positive integer. Thus, $F(t)$ is
approximately given by

$F(t)\approx\frac{\lambda\phi(t)v^{2}}{2\pi^{2}}\int^{N}_{0}d\tilde{k}\frac{\tilde{k}^{2}n_{k}(t)}{\tilde{\omega}_{k}(t)}\to\frac{\lambda\phi(t)v^{2}}{2\pi^{2}a^{3}(t)}\sum_{j=1}^{n}\left(\frac{N}{n}\right)^{3}\frac{j^{2}n_{k_{j}}(t)}{\tilde{\omega}_{k_{j}}(t)},$
(25)

where $n$ is a positive integer characterizing the lattice spacing and

$\tilde{k}_{j}=\frac{Nj}{n}.$ (26)

We may recover the continuous case by taking $N,n\to\infty$. Thus, we have
reduced the problem to be a set of $n$ differential equations.777We will not
include the zero mode since it does not contribute to the effective potential
within our approximation.

Using above approximation, we find a set of equations:

$\displaystyle\tilde{\xi}_{k_{i}}^{\prime\prime}(\tilde{t})+\tilde{\omega}_{k_{i}}^{2}\tilde{\xi}_{k_{i}}(\tilde{t})=\frac{1}{4}\tilde{\xi}_{k_{i}}^{-3}(\tilde{t})$
(27)
$\displaystyle\tilde{\phi}^{\prime\prime}(\tilde{t})+3\tilde{H}\tilde{\phi}(\tilde{t})+\frac{\lambda\tilde{\phi}(\tilde{t})}{2\pi^{2}a^{3}(\tilde{t})}\sum_{j=1}^{n}\left(\frac{N}{n}\right)^{3}\frac{j^{2}n_{k_{j}}(\tilde{t})}{\tilde{\omega}_{k_{j}}(\tilde{t})}=0,$
(28)

where $i=1,2,\cdots,n$, $\tilde{t}=vt$, $\tilde{\phi}=\phi/v$, $\tilde{H}=H/v$
and the prime denotes derivative with respect to $\tilde{t}$. We show
numerical solutions in Fig 1 with parameters and initial conditions
$\lambda=1,\tilde{\phi}(0)=10,\tilde{\phi}^{\prime}(0)=-1$ in Minkowski
spacetime $a(t)=1$. As quoted in Kofman:2004yc , in this model, parametric
resonance occurs after the second and subsequent zero crossings, which
extremely enhances the number density of low $k$ modes and strengthens the
trapping effect. In numerical simulations, we have observed that the numerical
solutions with different numbers of the mode number $n$ are qualitatively
similar to each other but quantitatively different.

Figure 1: The moduli dynamics in flat spacetime. We have numerically solved
the E.O.M and mode equations with the parameters
$\lambda=1,\tilde{\phi}(0)=10,\tilde{\phi}^{\prime}(0)=-1,N=2,n=10$.

As another illustration, we consider the background

$a(t)=\left(\frac{t}{t_{0}}\right)^{\beta}=\left(\frac{\tilde{t}}{\tilde{t}_{0}}\right)^{\beta},$
(29)

where $\beta$ is a positive constant, and we have normalized the scale factor
such that $a(\tilde{t}_{0})=1$ at the initial time $\tilde{t}=\tilde{t}_{0}$.
Then,

$\tilde{H}=\frac{\beta}{tv}=\frac{\beta}{\tilde{t}}.$ (30)

In Figs 2,3 we show the numerical solution for the matter dominated universe
$\beta=\frac{2}{3}$ and the radiation dominated universe $\beta=1/2$,
respectively. These figures show that the trapping mechanism works even in the
expanding background. Note however that, if the initial Hubble parameter is
large enough $(\tilde{t}_{0}\sim{\cal O}(1))$, it is possible to slow the
trapping by the dilution of particles as well as Hubble friction.

Figure 2: The moduli trapping effect in the matter dominated Universe. The
parameters are chosen as $\tilde{t}_{0}=25$, $N=3$, $n=10$. Figure 3: The
moduli trapping effect in the radiation dominated Universe. The parameters are
chosen as $\tilde{t}_{0}=55$, $N=3$, $n=9$.

Let us consider a more involved case that instead of $\lambda\phi^{2}$ the
mass term is given by

$M_{\chi}^{2}(\phi)=\mu^{2}(\sin\pi\phi+A\sin 3\pi\phi)^{2},$ (31)

where $A$ is a real parameter. In this case, there is a true symmetric point
$\phi=n\in\mathbb{Z}$ where $\chi$ becomes massless, but there is a fake
symmetric point as shown in Fig. 4. The E.O.M. of $\phi$ turns to

$\ddot{\phi}+2\pi\mu^{2}(\sin\pi\phi+A\sin 3\pi\phi)(\cos\pi\phi+3A\cos
3\pi\phi)\langle\chi^{2}\rangle=0.$ (32)

Figure 4: The behavior of the $\phi$-dependent mass (31) for
$\mu^{2}=1,A=0.7$. $\phi=n\in\mathbb{Z}$ are the true symmetric points and
there are fake points, where the mass becomes small but non-zero.

We show two numerical solutions with different parameters in Minkowski
spacetime $a(t)=1$ in Figs 5, 6. We note that the false vacua are
$\phi=\frac{2n+1}{2}$, ($n\in\mathbb{Z}$). In the former case, the modulus is
immediately trapped to the true symmetric vacuum $\phi=1$, whereas the latter
shows that the modulus is trapped both at the false and the true symmetric
vacua, but finally reaches the true symmetric vacuum. We expect that modulus
trapping at the true vacuum is not a generic property since the effective
potential vanishes even at the false vacua, which can be seen from the E.O.M.
of $\phi$. Nevertheless it is true that the particle production takes place
more efficiently near the true vacua since $\chi$ becomes much lighter at the
point. Therefore, if the true and false vacua are sufficiently separated, we
expect the modulus to be trapped near the enhanced symmetry point.

Figure 5: A numerical solution to (32) with
$\mu=10,n=10,A=0.5,\phi(0)=1.5,\dot{\phi}(0)=-1$. Figure 6: A numerical
solution to (32) with $\mu=7,n=10,A=0.5,\phi(0)=1.8,\dot{\phi}(0)=-1$.

Before closing this section, we finally comment on some technical issues
remaining in our numerical approach shown here. We have tried our numerical
simulations with different set of parameters or initial conditions several
times. Although the behavior qualitatively converges, the results are
quantitatively unstable under the change of parameters such as mode numbers.
Such a behavior may be understood from the fact that parametric resonance
become important when the modulus $\phi$ is almost trapped to the ESP, and the
resonance is sensitive to the parameters such as momentum. As we will discuss
later, the modular flavor symmetric models contain more issues in performing
fully numerical simulations. Therefore, we will propose a semi-analytic
approach to capture the moduli dynamics. Nevertheless, the numerical approach
we have taken in this section would still be useful for some class of models.

##  3 Modular flavor symmetric models

We show the modular flavor symmetric model considered in the following
sections. In this section, we first briefly review modular flavor symmetry and
then derive the E.O.M. of complex structure moduli as well as a spectator
scalar $\chi$ while we will study in detail the $\chi$-particle production in
the next section.

###  3.1 Modular flavor symmetry

The $SL(2,Z)$ group,

$\displaystyle\gamma=\begin{pmatrix}a&b\\\ c&d\end{pmatrix},$ (33)

where $a,b,c,d$ are integer and $ad-bc=1$, is generated by two generators,
$S$, and $T$,

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

They satisfy the following algebraic relations:

$\displaystyle S^{4}=1,\qquad(ST)^{3}=1.$ (35)

This group is referred to the homogeneous modular group, $\Gamma=SL(2,Z)$.

Under the modular symmetry, the modulus $\tau$ transforms as

$\displaystyle\gamma(\tau)=\frac{a\tau+b}{c\tau+d}.$ (36)

Note that $\tau$ is invariant under $S^{2}$. That is, the generators satisfy

$\displaystyle S^{2}=1,\qquad(ST)^{3}=1,$ (37)

on the modulus $\tau$. This group is referred to the inhomogeneous modular
symmetry, $\bar{\Gamma}=PSL(2,Z)=SL(2,Z)/Z_{2}$. Also we define the congruence
subgroup $\Gamma(N)$,

$\displaystyle\Gamma(N)=\left\\{\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\in\Gamma~{}\left|~{}\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}\right.({\rm
mod}~{}N)\right\\},$ (38)

which includes $T^{N}$. Similarly, we can define $\bar{\Gamma}(N)$.

The modular forms $f(\tau)_{i}$ are the holomorphic functions of $\tau$, which
transform under the modular symmetry,

$\displaystyle f(\gamma(\tau))_{i}=(c\tau+d)^{k}\rho_{ij}(\gamma)f(\tau)_{j},$
(39)

where $k$ is the modular weight and $\rho(\gamma)_{ij}$ is a unitary matrix.
Suppose that

$\displaystyle f(\gamma(\tau))_{i}=(c\tau+d)^{k}(\gamma)f(\tau)_{i},$ (40)

for $\gamma\in\bar{\Gamma}(N)$. Then, the matrix $\rho(\gamma)_{ij}$
represents the quotient, $\Gamma_{N}=\bar{\Gamma}/\bar{\Gamma}(N)$, where
$T^{N}=1$. Interestingly, these quotients $\Gamma_{N}$ with $N=2,3,4,5$ are
isomorphic to $S_{3}$, $A_{4}$, $S_{4}$, $A_{5}$, respectively.

There are fixed points on $\tau$, i.e. $\tau={\rm i},\omega=e^{2\pi{\rm
i}/3},{\rm i}\infty$, where residual symmetries remain. That is,
$Z_{2}^{S},Z_{3}^{ST},Z_{N}^{T}$ symmetries remain at $\tau={\rm
i},\omega=e^{2\pi{\rm i}/3},{\rm i}\infty$, respectively. The modular form
$f(\tau)$ has definite charges under these residual symmetries and their
behaviors are determined by their charges.

As illustrating examples, we study modular $A_{4}$ symmetric models in this
paper. Modular flavor symmetric superpotential can be written by

$\displaystyle W=Y_{ijk}(\tau)\Phi_{i}\Phi_{j}\Phi_{k},$ (41)

where $\Phi_{i}$ are the chiral superfields with modular weights $k_{i}$ and
they have some representations under $A_{4}$. The Yukawa couplings
$Y_{ijk}(\tau)$ are modular forms. The superpotential must be invariant under
the modular symmetry including the $A_{4}$ symmetry. The Kähler potential of
chiral matter fields can be written by

$\displaystyle K=\sum_{i}\frac{1}{(2{\rm Im}\tau)^{k_{i}}}|\Phi_{i}|.$ (42)

The fundamental modular forms of $A_{4}$ have the modular weight 2 and they
correspond to the $A_{4}$ triplet $(Y_{1}(\tau),Y_{2}(\tau),Y_{3}(\tau))$
Feruglio:2017spp . Their explicit forms are shown in Appendix B. The modular
forms of higher modular weights can be obtained by their tensor products. For
example, we use the modular forms of weight 8, because all of three $A_{4}$
singlets, ${\bf 1},{\bf 1}^{\prime}$, and ${\bf 1}^{\prime\prime}$ appear as
modular forms when the weight is 8, i.e. $Y^{(8)}_{\bf 1}(\tau),Y^{(8)}_{{\bf
1}^{\prime}}(\tau)$, and $Y^{(8)}_{{\bf 1}^{\prime\prime}}(\tau)$, which are
shown in Appendix B. The three $A_{4}$ singlets, ${\bf 1},{\bf 1}^{\prime}$,
and ${\bf 1}^{\prime\prime}$ have the $Z_{3}^{ST}$ charges, 0, 1, 2,
respectively. The corresponding modular forms have additional contributions
from the automorphic factors $(c\tau+d)^{k}$.

As a simple illustrating model, in the following sections we consider a single
complex scalar spectator field that has the mass term

$\displaystyle|M_{\chi}(\tau,\bar{\tau})|^{2}|\chi|^{2},$ (43)

where $M_{\chi}(\tau,\bar{\tau})|^{2}$ is given by modular forms. Such a
complex scalar field can appear in supersymmetric models with

$\displaystyle K=$ $\displaystyle\frac{1}{(2{\rm Im}\tau)^{k}}|\Phi|^{2},$
(44) $\displaystyle W=$ $\displaystyle mY(\tau)\Phi S,$ (45)

where $\Phi$ is a chiral superfield whose scalar component is $\chi$, $k$ is
the modular weight of $\Phi$, $m$ a mass parameter, and $Y(\tau)$ a
holomorphic function of $\tau$ given by modular forms. $S$ is an additional
chiral superfield which is a singlet under modular symmetry.888Since $S$ is a
singlet, it may have a heavy mass term independent of the above one. If $S$
has a heavy mass, it can decouple from the theory even at ESPs where
$Y(\tau)\to 0$. Although there are fermions in such a model, which are as
light as the scalar $\chi$, we only discuss the complex scalar $\chi$ only.
The fermionic particle production can be similarly discussed along the line
e.g. of Greene:2000ew ; Peloso:2000hy ; Adshead:2015kza  but would be more
involved than the bosonic case we consider here because of the chiral
structure.

We would like to comment on the relation to the standard model. In the
supersymmetric standard model, mass terms of scalars and fermions through
Yukawa couplings are not generated until the electroweak symmetry is
spontaneously broken. Therefore, we do not expect standard model particles to
contribute the moduli trapping.999Scalars would have mass terms by
supersymmetry breaking and they generally depend on moduli. Therefore, the
moduli dynamics may generate scalars in the supersymmetric standard model.
Instead, we can identify $\Phi$ as a superfield that obtains its mass term
through the GUT symmetry breaking101010Another candidate would be right-handed
(s)neutrinos, where the superpotential can be written by
$W=mY(\tau)\Phi^{2}$.. Therefore, the mass parameter $m$ can be sufficiently
large (but smaller than $M_{\rm pl}$).

###  3.2 Dynamics of moduli and spectator scalars

We consider quantization of a spectator complex scalar field $\chi$ that
couples to the complex structure modulus $\tau$, and the action of $\chi$ is
given by

$S=-\int d^{4}x\sqrt{-g}\left[\frac{1}{(2{\rm
Im}\tau)^{k}}|\partial\chi|^{2}+|M_{\chi}(\tau,\bar{\tau})|^{2}|\chi|^{2}\right].$
(46)

We introduce quasi-canonical real fields

$\tau=\sqrt{2}\theta/M_{\rm pl}+{\rm i}e^{\sqrt{2}\phi/M_{\rm pl}},$ (47)

which will become a choice that makes $\phi$ canonical. We assume that the
“bare mass” $M^{2}(\tau,\bar{\tau})$ is given by modular forms as

$|M_{\chi}(\tau,\bar{\tau})|^{2}=m^{2}\left|Y^{(8)}_{\bm{1}}(\tau)\right|^{2}\
{\rm or}\ m^{2}\left|Y^{(8)}_{\bm{1}^{\prime}}(\tau)\right|^{2},$ (48)

where $m$ denotes some mass scale such as the VEV of a GUT Higgs field, and
the background spacetime is $ds^{2}=-dt^{2}+a^{2}(t)d{\bm{x}}^{2}$.
Introducing $\chi=\tilde{\chi}_{1}+{\rm i}\tilde{\chi}_{2}$, we rewrite the
action as

$\displaystyle S=$ $\displaystyle\int
dtd^{3}x\left[\frac{1}{2^{k}}a^{3}e^{-\sqrt{2}k\phi/M_{\rm
pl}}\left(\dot{\tilde{\chi}}_{i}^{2}-\frac{1}{a^{2}}({\bm{\nabla}}\tilde{\chi}_{i})^{2}\right)-a^{3}|M_{\chi}(\tau,\bar{\tau})|^{2}\tilde{\chi}^{2}_{i}\right],$
(49)

which can be made canonical by introducing a new basis
$\tilde{\chi}_{i}=\alpha\chi_{i}$ where

$\alpha\equiv\frac{2^{\frac{k-1}{2}}}{a^{\frac{3}{2}}}e^{\frac{\sqrt{2}}{2}k\phi/M_{\rm
pl}}.$ (50)

With the new basis, we find

$\displaystyle S=\sum_{i=1}^{2}\int
dtd^{3}x\left[\frac{1}{2}\left(\dot{\chi}_{i}^{2}-\frac{1}{a^{2}}({\bm{\nabla}}\chi_{i})^{2}\right)-\frac{1}{2}M_{\rm
eff}^{2}\chi_{i}^{2}\right],$ (51)

where we have assumed that the background field depends only on $t$, and used
integration by parts in the second equality. Thus, the effective mass of the
canonically normalized spectator field is

$\displaystyle M_{\rm eff}^{2}=$ $\displaystyle 2^{k}e^{\sqrt{2}k\phi/M_{\rm
pl}}|M_{\chi}|^{2}-2\frac{\dot{\alpha}^{2}}{\alpha^{2}}+\frac{\ddot{\alpha}}{\alpha}$
$\displaystyle=$ $\displaystyle 2^{k}e^{\sqrt{2}k\phi/M_{\rm
pl}}|M_{\chi}|^{2}-\frac{k^{2}}{2M_{\rm
pl}^{2}}\dot{\phi}^{2}+\frac{3\sqrt{2}k}{2M_{\rm
pl}}\dot{\phi}H+\frac{\sqrt{2}k}{2M_{\rm
pl}}\ddot{\phi}-\frac{9}{4}H^{2}-\frac{3}{2}\dot{H}$ $\displaystyle\approx$
$\displaystyle 2^{k}e^{\sqrt{2}k\phi/M_{\rm pl}}|M_{\chi}|^{2}.$ (52)

The time derivative terms of $\phi$ in the second line may be interpreted as
the higher-dimensional curvature induced mass as ${\rm Im}\tau$ is related to
the area of the extra dimensional torus. We will neglect such terms as well as
Hubble induced terms in the following.111111Since the mass becomes zero at the
critical point and the low $k$-mode can be tachyonic, it is generally non-
trivial if such approximation is allowed. We may neglect such effect by
assuming that the spectator scalar has a small mass that is comparable to the
derivative terms so that the tachyonic instability disappears. On the other
hand, the instability near the critical point may enhance the particle
production. We also notice that the effective mass is modular invariant as it
should be.

Let us consider the effective action of the complex structure modulus $\tau$
given by

$\displaystyle S$ $\displaystyle=-\int d^{4}x\sqrt{-g}\left[\frac{M_{\rm
pl}^{2}}{(2{\rm
Im}\tau)^{2}}|\partial\tau|^{2}+V(\tau,\bar{\tau})+a^{-3}M_{\rm
eff}^{2}\sum_{i=1}^{2}\langle\chi_{i}^{2}\rangle\right]$ $\displaystyle=-\int
d^{4}x\sqrt{-g}\left[\frac{1}{2}(\partial\phi)^{2}+\frac{1}{2}e^{-2\sqrt{2}\phi/M_{\rm
pl}}(\partial\theta)^{2}+V(\tau,\bar{\tau})+a^{-3}M_{\rm
eff}^{2}\sum_{i=1}^{2}\langle\chi_{i}^{2}\rangle\right].$ (53)

Note that the appearance of non-covariant expression $+a^{-3}M_{\rm
eff}^{2}\sum_{i=1}^{2}\langle\chi_{i}^{2}\rangle$ is the result of making
$\chi_{i}$ canonical.121212One could use the original variable, but would be
more involved since both the mass and the kinetic term of $\chi$ appear in the
E.O.M. of $\phi$ through the non-minimal coupling in the kinetic term. The
equations of motions of $\phi$, $\theta$ are given by

$\displaystyle-\Box\phi-\sqrt{2}e^{-2\sqrt{2}\phi/M_{\rm
pl}}(\partial\theta)^{2}-\frac{2\sqrt{2}e^{\sqrt{2}\phi/M_{\rm pl}}}{M_{\rm
pl}}{\rm Im}(\partial_{\tau}V)+\frac{1}{a^{3}}\partial_{\phi}M_{\rm
eff}^{2}\sum_{i=1}^{2}\langle\chi_{i}^{2}\rangle=0,$ (54)
$\displaystyle-\Box\theta+\frac{2\sqrt{2}}{M_{\rm
pl}}\partial_{\mu}\phi\partial^{\mu}\theta+\frac{2\sqrt{2}e^{2\sqrt{2}\phi/M_{\rm
pl}}}{M_{\rm pl}}{\rm Re}(\partial_{\tau}V)+e^{2\sqrt{2}\phi/M_{\rm
pl}}\frac{1}{a^{3}}\partial_{\theta}M_{\rm
eff}^{2}\sum_{i=1}^{2}\langle\chi_{i}^{2}\rangle=0.$ (55)

It is straightforward to generalize our previous discussion to this case by
using

$\omega_{p}^{2}=\frac{p^{2}}{a^{2}}+M_{\rm
eff}^{2}\approx\frac{p^{2}}{a^{2}}+2^{k}e^{\sqrt{2}k\phi/M_{\rm
pl}}|M_{\chi}|^{2}.$ (56)

The vacuum expectation value of $\langle\hat{\chi}_{1}^{2}\rangle$ is given by

$\displaystyle\langle\hat{\chi}_{1}^{2}\rangle\approx\frac{a^{3}M_{\rm
eff}^{2}}{8\pi^{2}\epsilon}-\frac{a^{3}M_{\rm
eff}^{2}}{16\pi^{2}}\left[1-\gamma_{E}-\log\left(\frac{M_{\rm
eff}^{2}}{4\pi\mu^{2}}\right)\right]+\int\frac{d^{3}p}{(2\pi)^{3}}\frac{n_{p}(t)}{\omega_{p}(t)}.$
(57)

Thus, the renormalized effective potential is found to be

$\displaystyle V_{\rm eff}=$ $\displaystyle\frac{M_{\rm
eff}^{4}}{8\pi^{2}\epsilon}-\frac{M_{\rm
eff}^{4}}{16\pi^{2}}\left[1-\gamma_{E}-\log\left(\frac{M_{\rm
eff}^{2}}{4\pi\mu^{2}}\right)\right]+\delta_{\rm CT}+\frac{M_{\rm
eff}^{2}}{a^{3}}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{n_{p}(t)}{\omega_{p}(t)},$
(58)

where $\delta_{\rm CT}$ denotes possible local counter-terms eliminating
divergent pieces.131313We point out that the “effective potential” contains
the derivative of the background field $\phi$ in general, which behaves as
higher-derivative terms. Such terms should be removed by appropriate counter-
terms in order to avoid (possibly) unphysical ghost degrees of freedom. As the
previous sections, we assume the cancellation of all but the last term in
(58).

In analysing the dynamics of the moduli fields $\phi,\theta$, there are
several technical issues in this model in addition to the ones discussed
within toy models: Since the coupling between $\chi$ and moduli are suppressed
by Planck scale which causes hierarchically small or large numbers, and the
moduli dependence is quite complicated, the numerical costs become more than
that in the toy models. More specifically, the effective field theory
description requires the field velocity and the Hubble parameter to be much
smaller than the Planck scale. Furthermore, if the moduli trapping occurs
after inflation, the Hubble scale needs to be less than about ${\cal
O}(10^{13}{\rm GeV})$ by the constraint on the tensor-to-scalar ratio.
Therefore, it would be natural to assume the ratio between the initial field
velocity and Planck scale to be below $\mathcal{O}(10^{-5})$, which appear in
couplings if we normalize all dimensionful parameters by the scale of initial
field velocity. The complication of the mass term given by modular forms shown
later further makes the numerical simulations difficult. Therefore, we will
develop semi-analytic approaches to discuss the moduli dynamics within modular
flavor symmetric models in the next section.

##  4 Particle production at ESPs in modular flavor symmetric models

In order to overcome some technical difficulties within our setup, we consider
a semi-analytic approach, where we analytically estimate the number density of
the spectator scalar particle produced at the first particle production time
at which $M^{2}_{\chi}(t_{0})\approx 0$. Using the estimate, we analytically
evaluate the effective potential arising after the first particle production
event, and numerically solve the E.O.M. of moduli fields with the estimated
effective potential, which does not contain the difficulties mentioned before.
Although this approach cannot capture the subsequent particle production
events, we expect such events just strengthen the trapping potential, which
does not change the dynamics of moduli qualitatively.

We first show analytic formulas for particle number density produced at the
crossing of the ESP. To do so, we need to know the behavior of the effective
mass near the ESPs. In modular symmetric models, the ESPs are $\tau={\rm
i},e^{2{\rm i}\pi/3}(\equiv\omega),+{\rm i}\infty$ where some modular forms
vanish. In the following discussion, we focus on $\tau=\omega$ as a
representative case, but similar analysis can be done for any other critical
points in a similar way.

In the following subsections, we first discuss the behavior of the effective
mass of the spectator scalar field near the ESP $\tau\sim\omega$, which
enables us to analytically estimate the amount of particle number density
produced after the crossing with the ESP.

With the aid of the analytic expression of the effective mass around
$\tau\sim\omega$, we discuss one-dimensional dynamics of $\theta$ or $\phi$
where one of them is fixed to a constant value. We give analytic formulas for
the particle number density as well as the effective forces of the produced
particle on the modulus.

###  4.1 The behavior of the effective mass near the critical point

The modular forms can be classified as representations of subgroup. The
singlet modular forms of weight 8,
$Y^{(8)}_{{\bm{1}},{\bm{1}^{\prime}},{\bm{1}^{\prime\prime}}}$ transform as
Novichkov:2021evw

$\displaystyle\left(\begin{array}[]{c}Y^{(8)}_{\bm{1}}(-(\tau+1)^{-1})\\\
Y^{(8)}_{{\bm{1}^{\prime}}}(-(\tau+1)^{-1})\\\
Y^{(8)}_{{\bm{1}^{\prime\prime}}}(-(\tau+1)^{-1})\end{array}\right)=(-(\tau+1))^{8}\left(\begin{array}[]{ccc}1&0&0\\\
0&\omega&0\\\
0&0&\omega^{2}\end{array}\right)\left(\begin{array}[]{c}Y^{(8)}_{{\bm{1}}}(\tau)\\\
Y^{(8)}_{{\bm{1}^{\prime}}}(\tau)\\\
Y^{(8)}_{{\bm{1}^{\prime\prime}}}(\tau)\end{array}\right),$ (68)

which can be equivalently written as

$\displaystyle\left(\begin{array}[]{c}Y^{(8)}_{\bm{1}}(\omega^{2}u)\\\
Y^{(8)}_{{\bm{1}^{\prime}}}(\omega^{2}u)\\\
Y^{(8)}_{{\bm{1}^{\prime\prime}}}(\omega^{2}u)\end{array}\right)=\left(\frac{1-\omega^{2}u}{1-u}\right)^{8}\left(\begin{array}[]{ccc}\omega^{-8}&0&0\\\
0&\omega^{-7}&0\\\
0&0&\omega^{-6}\end{array}\right)\left(\begin{array}[]{c}Y^{(8)}_{{\bm{1}}}(u)\\\
Y^{(8)}_{{\bm{1}^{\prime}}}(u)\\\
Y^{(8)}_{{\bm{1}^{\prime\prime}}}(u)\end{array}\right),$ (78)

where

$u\equiv\frac{\tau-\omega}{\tau-\omega^{2}},$ (79)

is a variable for the deviation from the symmetric point $\tau=\omega$.
Therefore, we obtain

$\displaystyle\left(\begin{array}[]{c}\tilde{Y}^{(8)}_{\bm{1}}(\omega^{2}u)\\\
\tilde{Y}^{(8)}_{{\bm{1}^{\prime}}}(\omega^{2}u)\\\
\tilde{Y}^{(8)}_{{\bm{1}^{\prime\prime}}}(\omega^{2}u)\end{array}\right)=\left(\begin{array}[]{ccc}\omega^{-8}&0&0\\\
0&\omega^{-7}&0\\\
0&0&\omega^{-6}\end{array}\right)\left(\begin{array}[]{c}\tilde{Y}^{(8)}_{{\bm{1}}}(u)\\\
\tilde{Y}^{(8)}_{{\bm{1}^{\prime}}}(u)\\\
\tilde{Y}^{(8)}_{{\bm{1}^{\prime\prime}}}(u)\end{array}\right),$ (89)

where $\tilde{Y}^{(8)}_{\bm{r}}(u)\equiv(1-u)^{-8}Y^{(8)}_{\bm{r}}(u)$.
Expanding both sides with respect to $\omega$ yields

$\left(\omega^{2l}-\omega^{q_{\bm{r}}-8}\right)\frac{d^{l}\tilde{Y}^{(8)}_{\bm{r}}(u)}{du^{l}}\biggr{|}_{u\to
0}=0,$ (90)

where $q_{\bm{r}}=0,1,2$ for
${\bm{1}},{\bm{1}^{\prime}},{\bm{1}^{\prime\prime}}$, respectively. This
relation implies that
$\frac{d^{l}\tilde{Y}^{(8)}_{\bm{r}}(u)}{du^{l}}\biggr{|}_{u\to 0}=0$ unless
$2l=q_{\bm{r}}-8\ (\text{mod 3})$. Thus,
$\tilde{Y}^{(8)}_{\bm{1},\bm{1}^{\prime}}(0)=0$, and
$\frac{d\tilde{Y}^{(8)}_{\bm{1},\bm{1}^{\prime\prime}}}{du}(0)=0$. Noting that

$\displaystyle\frac{\partial u}{\partial\tau}=\frac{(1-u)^{2}}{\sqrt{3}{\rm
i}},$ (91)

we find

$\displaystyle
Y^{(8)}_{\bm{r}}(\tau)|_{\tau\to\omega}=\tilde{Y}^{(8)}_{\bm{r}}(u)|_{u\to
0},$ (92)
$\displaystyle\frac{dY^{(8)}_{\bm{r}}}{d\tau}\Biggr{|}_{\tau\to\omega}=\frac{1}{\sqrt{3}{\rm
i}}\left[-8\tilde{Y}^{(8)}_{\bm{r}}+\frac{d\tilde{Y}^{(8)}_{\bm{r}}}{du}\right]_{u\to
0},$ (93)
$\displaystyle\frac{d^{2}Y^{(8)}_{\bm{r}}}{d\tau^{2}}\Biggr{|}_{\tau\to\omega}=\left[-24\tilde{Y}^{(8)}_{\bm{r}}+6\frac{d\tilde{Y}_{\bm{r}}^{(8)}}{du}-\frac{1}{3}\frac{d^{2}\tilde{Y}_{\bm{r}}^{(8)}}{du^{2}}\right]_{u\to
0}.$ (94)

Note also that

$\displaystyle
u=\frac{\tau-\omega}{\omega-\omega^{2}}+\mathcal{O}((\tau-\omega)^{2})=\frac{\tau-\omega}{\sqrt{3}{\rm
i}}+\mathcal{O}((\tau-\omega)^{2}).$ (95)

The leading order terms of each singlet
${\bm{1}},{\bm{1}^{\prime}},{\bm{1}^{\prime\prime}}$ are as follows:

$\displaystyle Y^{(8)}_{\bm{1}}(\tau)=$
$\displaystyle-\frac{1}{6}\frac{d^{2}\tilde{Y}_{\bm{1}}^{(8)}(0)}{du^{2}}(\tau-\omega)^{2}+\mathcal{O}((\tau-\omega)^{3}),$
(96) $\displaystyle Y^{(8)}_{\bm{1^{\prime}}}(\tau)=$
$\displaystyle\frac{1}{\sqrt{3}{\rm
i}}\frac{d\tilde{Y}^{(8)}_{\bm{1}^{\prime}}(0)}{du}(\tau-\omega)+\mathcal{O}((\tau-\omega)^{2}),$
(97) $\displaystyle Y^{(8)}_{\bm{1^{\prime\prime}}}(\tau)=$
$\displaystyle\tilde{Y}^{(8)}_{\bm{1^{\prime\prime}}}(0)+\mathcal{O}((\tau-\omega)^{1}),$
(98)

and the last one shows that $\bm{1}^{\prime\prime}$ cannot produce the
particles at ESP $\tau=\omega$, because $\chi$ is still massive.141414More
precisely speaking, for $\bm{r}=\bm{1}^{\prime\prime}$, particles are produced
if $\tau=\omega$ is the local minimum of $Y_{\bm{1}^{\prime\prime}}^{(8)}$,
but as the leading term is non vanishing, the particle production would be
less than that of $\bm{r}=\bm{1}\ \text{or}\ \bm{1}^{\prime}$. Thus, the
effective mass can be approximated as

$\displaystyle M_{\rm eff}^{2}=$
$\displaystyle\frac{2^{6}|C|^{2}}{9}m^{2}e^{8\sqrt{2}\phi/M_{\rm
pl}}|\tau-\omega|^{4}+\cdots$ $\displaystyle=$
$\displaystyle\frac{9|C|^{2}}{4}m^{2}\left[\left(\frac{\sqrt{2}\theta}{M_{\rm
pl}}+\frac{1}{2}\right)^{2}+\left(e^{\sqrt{2}\phi/M_{\rm
pl}}-\frac{\sqrt{3}}{2}\right)^{2}\right]^{2}+\cdots(\text{for
${\bm{r}=\bm{1}}$}),$ (99) $\displaystyle M_{\rm eff}^{2}=$
$\displaystyle\frac{2^{8}|D|^{2}}{3}m^{2}e^{8\sqrt{2}\phi/M_{\rm
pl}}|\tau-\omega|^{2}+\cdots$ $\displaystyle=$ $\displaystyle
27|D|^{2}m^{2}\left[\left(\frac{\sqrt{2}\theta}{M_{\rm
pl}}+\frac{1}{2}\right)^{2}+\left(e^{\sqrt{2}\phi/M_{\rm
pl}}-\frac{\sqrt{3}}{2}\right)^{2}\right]+\cdots(\text{for
${\bm{r}=\bm{1}^{\prime}}$}),$ (100)

where $C\equiv\frac{d^{2}\tilde{Y}^{(8)}_{\bm{1}}(0)}{du^{2}}$,
$D\equiv\frac{d\tilde{Y}^{(8)}_{\bm{1}^{\prime}}(0)}{du}$, and ellipses denote
the terms higher order in $(\tau-\omega)$.

Before going to the details of particle production, we would like give a few
comments on the behavior of the effective mass for each representation. We
notice that depending on the representation $\bm{r}$, the behavior of the mass
term near the ESP $\tau=\omega$ changes. Since the moduli dependence comes
with powers of $\chi/M_{\rm pl}$ or $\theta/M_{\rm pl}$, we can say the
particle production in $\bm{r}=\bm{1}$ case should be much smaller than that
in $\bm{r}=\bm{1}^{\prime}$. Therefore, the moduli trapping effect becomes
more significant for $\bm{r}={\bm{1}}^{\prime}$. Such observation is based on
the fact that moduli are gravitationally coupled to matter fields, which are
completely independent of particle production dynamics discussed below.

###  4.2 1D dynamics: $\phi$-fixed

In this subsection, we discuss the particle production due to the $\theta$
dynamics while $\phi$ is fixed at the ESP $\phi=\phi_{0}$ where
$e^{\sqrt{2}\phi_{0}/M_{\rm pl}}=\frac{\sqrt{3}}{2}$, and we parametrize
$\theta$ as

$\theta=v(t-t_{0})-\frac{1}{2\sqrt{2}}M_{\rm pl},$ (101)

where $t_{0}$ is the time at which $\theta$ crosses the enhanced symmetry
point and $v>0$ is the velocity of $\theta$ at $t=t_{0}$. This parametrization
is a good approximation if Hubble parameter is smaller than $\sqrt{v}$. (See
appendix C.) Then, near $t\sim t_{0}$ the effective mass can be approximated
by

$\displaystyle M_{\rm eff}^{2}\approx$ $\displaystyle
9|C|^{2}m^{2}\left(\frac{v(t-t_{0})}{M_{\rm pl}}\right)^{4}\ (\text{for
${\bm{r}=\bm{1}}$}),$ (102) $\displaystyle M_{\rm eff}^{2}\approx$
$\displaystyle 54|D|^{2}m^{2}\left(\frac{v(t-t_{0})}{M_{\rm pl}}\right)^{2}\
(\text{for ${\bm{r}=\bm{1}^{\prime}}$}),$ (103)

and the effective frequency for ${\bm{k}}$-mode is

$\omega_{k}^{2}\approx\frac{k^{2}}{a^{2}(t_{0})}+M_{\rm eff}^{2},$ (104)

where we have approximated the scale factor by its value at $t=t_{0}$. We will
show the analytic estimate of the particle production taking place around
$t\sim t_{0}$ separately for both $\bm{r}=\bm{1}$ and $\bm{r}=\bm{1}^{\prime}$
respectively.

$\bm{r}=\bm{1}$ case

In order to estimate the particle production rate, we need to find the turning
points at which $\omega_{k}^{2}(t)=0$ in a complex $t$-plane. (See
Dumlu:2010ua ; Dabrowski:2014ica ; Dabrowski:2016tsx ; Enomoto:2020xlf ;
Taya:2020dco ; Hashiba:2021npn ; Yamada:2021kqw  for review.) With (102) the
turning points are found to be

$t_{1}^{\pm}-t_{0}=e^{\pm{\rm i}\pi/4}\frac{M_{\rm
pl}}{v}\left(\frac{k^{2}}{9|C|^{2}m^{2}a^{2}(t_{0})}\right)^{\frac{1}{4}},\quad
t_{2}^{\pm}-t_{0}=e^{\pm 3{\rm i}\pi/4}\frac{M_{\rm
pl}}{v}\left(\frac{k^{2}}{9|C|^{2}m^{2}a^{2}(t_{0})}\right)^{\frac{1}{4}},$
(105)

with which the effective frequency can be rewritten as

$\displaystyle\omega_{k}^{2}\approx
R\left(-(t_{*}-t_{0})^{4}+(t-t_{0})^{4}\right),$ (106)

where $t_{*}$ can be any of turning points and

$R\equiv\frac{9|C|^{2}m^{2}v^{4}}{M_{\rm pl}^{4}}.$ (107)

Notice that there are two pairs of turning points $(t_{i}^{+},t_{i}^{-})$
$(i=1,2)$ which are complex conjugate to each other, and Stokes lines
connecting them crosses the real $t$-axis. The amount of particle production
can be approximately given by a simple formula (see e.g. Dumlu:2010ua ;
Dabrowski:2014ica ; Dabrowski:2016tsx ; Enomoto:2020xlf ; Taya:2020dco ;
Hashiba:2021npn ; Yamada:2021kqw  for details151515See also Froman:1970toy .)

$n_{k}(t)=|\beta_{k}(t_{i})|^{2}\approx\left|e^{F_{1}}+e^{{\rm
i}\theta_{12}}e^{F_{2}}\right|^{2},$ (108)

where we have assumed that $t$ is sufficiently later than $t_{0}$, and

$\displaystyle F_{i}\equiv$ $\displaystyle\exp\left({\rm
i}\int_{t_{i}^{-}}^{t_{i}^{+}}\omega_{k}(t^{\prime})dt^{\prime}\right),$ (109)
$\displaystyle\theta_{12}\equiv$ $\displaystyle\exp\left(2{\rm
i}\int_{t_{2}^{+}}^{t_{1}^{+}}\omega_{k}(t^{\prime})dt^{\prime}\right).$ (110)

More explicitly we find

$\displaystyle F_{i}=$
$\displaystyle\exp\left[-\frac{\Gamma^{2}(\frac{1}{4})}{6\sqrt{\pi}R^{\frac{1}{4}}}\left(\frac{k}{a(t_{0})}\right)^{\frac{3}{2}}\right],$
(111) $\displaystyle\theta_{12}=$
$\displaystyle\frac{\Gamma^{2}(\frac{1}{4})}{3\sqrt{\pi}R^{\frac{1}{4}}}\left(\frac{k}{a(t_{0})}\right)^{\frac{3}{2}}.$
(112)

Therefore, the particle number density after crossing two Stokes lines can be
approximately estimated as

$\displaystyle\Delta n_{k}=2e^{-\gamma_{k}}(1+\cos\gamma_{k}),$ (113)

where

$\gamma_{k}\equiv\frac{\Gamma^{2}(\frac{1}{4})}{3\sqrt{\pi}R^{\frac{1}{4}}}\left(\frac{k}{a(t_{0})}\right)^{\frac{3}{2}}.$
(114)

We have numerically checked the validity of this formula, and find a very good
agreement of the numerical result and the analytic formula as shown in Fig. 7.
Thus, we can estimate the total particle production as well as the effective
potential on the basis of the approximate formula (113).

Figure 7: Comparison between the analytic formula (113) (the blue line) and
numerical results (red dots). The momentum $p$ is in units of the velocity
$\sqrt{|v|}$ at the particle production time $t_{0}$. In this example, we have
taken $vt_{0}=50$, $vt_{\rm ini}=15$, $\beta=\frac{1}{2}$, $\sqrt{|v|}/M_{\rm
pl}=5\times 10^{-4}$, $m/\sqrt{|v|}=5\times 10^{2}$. We have also checked that
the agreement is quite well for various sets of parameters.

With the approximate particle number formula for a $k$-mode (113), we are able
to estimate the total number density to be161616The oscillatory part
$e^{\gamma_{k}}\cos\gamma_{k}$ turns out to be vanishing in the momentum
integration.

$\displaystyle\Delta N=$ $\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\Delta
n_{k}=\frac{2\times
4\pi}{8\pi^{3}}\left(\frac{\Gamma^{2}(\frac{1}{4})}{3\sqrt{\pi}R^{\frac{1}{4}}a^{\frac{3}{2}}(t_{0})}\right)^{-2}\times\frac{2}{3}=\frac{6\sqrt{R}a^{3}(t_{0})}{\pi\Gamma^{4}(\frac{1}{4})}.$
(115)

Even though we have an approximate formula of the produce particle number
density, it is still difficult to evaluate the effective force given by

$F_{\theta}(t)=\frac{{(\rm Im}\tau)^{2}\partial_{\theta}(M_{\rm
eff}^{2})}{a^{3}}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{\Delta
n_{p}}{\omega_{p}(t)}.$ (116)

Since the effective support of the integrand is localized to small $k$, we may
approximate $\omega_{p}\approx M_{\rm eff}$ in the integrand, which allows us
to approximate $F_{\theta}(t)$ as

$F_{\theta}(t)\approx\frac{3\partial_{\theta}(M_{\rm eff}^{2})}{4a^{3}M_{\rm
eff}}\Delta
N\Theta(t-t_{0})=\frac{9\sqrt{R}a^{3}(t_{0})\partial_{\theta}(M_{\rm
eff}^{2})}{2\pi\Gamma^{4}(\frac{1}{4})M_{\rm eff}a^{3}(t)}\Theta(t-t_{0}),$
(117)

where we have introduced the Heaviside theta function $\Theta(t-t_{0})$.

$\bm{r}=\bm{1}^{\prime}$ case

Next, we consider the case $\bm{r}=\bm{1}^{\prime}$ where the effective mass
is given by (103), and find a pair of turning points

$\displaystyle t^{\pm}_{1}-t_{0}=\pm{\rm i}\frac{M_{\rm
pl}k}{3\sqrt{6}|D|mv}.$ (118)

Similarly to the previous case, the particle number produced at $t=t_{0}$ can
be evaluated by

$\displaystyle n_{k}(t_{0})\approx$ $\displaystyle\exp\left(2{\rm
i}\int_{t_{1}^{-}}^{t_{1}^{+}}\omega_{k}(t^{\prime})dt^{\prime}\right)$
$\displaystyle=$ $\displaystyle\exp\left[-\frac{\pi k^{2}M_{\rm
pl}}{3\sqrt{6}|D|mv}\right].$ (119)

We show the comparison of our analytic formula (119) and numerical results in
Fig. 8. Again, we have found an excellent agreement between them.

Figure 8: Comparison between the analytic formula (119) (the blue line) and
numerical results (red dots). The momentum $p$ is in units of the velocity
$\sqrt{|v|}$ at the particle production time $t_{0}$. Here, we have taken
$vt_{0}=70$, $vt_{\rm ini}=5$, $\beta=\frac{1}{2}$, $\sqrt{|v|}/M_{\rm
pl}=10^{-4}$, $m/\sqrt{|v|}=10$. We have also checked that the agreement is
quite well for various sets of parameters.

Accordingly, the total particle density is found to be

$\Delta
N=\int\frac{d^{3}k}{(2\pi)^{3}}n_{k}(t_{0})=\frac{1}{2\pi^{2}}\times\frac{\sqrt{\pi}}{4}\left(\frac{3\sqrt{6}|D|mv}{\pi
M_{\rm
pl}}\right)^{\frac{3}{2}}=\frac{1}{8\pi^{3}}\left(\frac{3\sqrt{6}|D|mv}{M_{\rm
pl}}\right)^{\frac{3}{2}},$ (120)

and the effective force is approximately given by

$\displaystyle F_{\theta}(t)=\frac{3a^{3}(t_{0})\partial_{\theta}(M_{\rm
eff}^{2})}{4a^{3}(t)}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{n_{p}(t)}{\omega_{p}(t)}\approx\frac{3a^{3}(t_{0})(\partial_{\theta}M_{\rm
eff}^{2})(t)}{32\pi^{3}a^{3}(t)M_{\rm
eff}(t)}\left(\frac{3\sqrt{6}|D|mv}{M_{\rm
pl}}\right)^{\frac{3}{2}}\Theta(t-t_{0}).$ (121)

###  4.3 1D dynamics: $\theta$-fixed

We consider the case that $\theta$ is fixed at $\theta=-\frac{M_{\rm
pl}}{2\sqrt{2}}$ and parametrize $\phi$ as

$\phi=\frac{M_{\rm
pl}\log\left(\frac{\sqrt{3}}{2}\right)}{\sqrt{2}}+v(t-t_{0}),$ (122)

with which the effective mass can be rewritten as

$\displaystyle M_{\rm eff}^{2}=$
$\displaystyle\frac{81|C|^{2}}{64}m^{2}\left(e^{\frac{\sqrt{2}v(t-t_{0})}{M_{\rm
pl}}}-1\right)^{4}+\cdots$ $\displaystyle\approx$
$\displaystyle\frac{81|C|^{2}}{16}m^{2}\left(\frac{v(t-t_{0})}{M_{\rm
pl}}\right)^{4}\quad(\text{for ${\bm{r}=\bm{1}}$}),$ (123) $\displaystyle
M_{\rm eff}^{2}=$
$\displaystyle\frac{81|D|^{2}m^{2}}{4}\left(e^{\frac{\sqrt{2}v(t-t_{0})}{M_{\rm
pl}}}-1\right)^{2}+\cdots$ $\displaystyle\approx$
$\displaystyle\frac{81|D|^{2}m^{2}}{2}\left(\frac{v(t-t_{0})}{M_{\rm
pl}}\right)^{2}\quad(\text{for ${\bm{r}=\bm{1}^{\prime}}$}),$ (124)

where we have expanded the terms higher order in $v(t-t_{0})/M_{\rm pl}$ since
it is generally small. We notice that these expressions are the same as the
$\phi$-fixed case (102), (103) by replacing $|C|^{2}\to\frac{9}{16}|C|^{2}$
and $|D|^{2}\to\frac{3}{4}|D|^{2}$, respectively. Thus, the particle number
density of the $k$-mode is given by

$\displaystyle\Delta
n_{k}=\Biggl{\\{}\begin{array}[]{c}2e^{-\frac{2}{\sqrt{3}}\gamma_{k}}\left(1+\cos\frac{2}{\sqrt{3}}\gamma_{k}\right)\quad(\text{for
${\bm{r}=\bm{1}}$})\\\ \exp\left[-\frac{2\pi k^{2}M_{\rm
pl}}{9\sqrt{2}|D|mv}\right]\quad(\text{for
${\bm{r}=\bm{1}^{\prime}}$}).\end{array}$ (127)

Therefore, the effective force on $\phi$ due to particle production can be
written as171717Recall that $R\propto|C|^{2}$ and also that $F_{\phi}$ does
not have an extra factor $({\rm Im}\tau)^{2}$ while $F_{\theta}$ does.

$\displaystyle
F_{\phi}(t)\approx\frac{18\sqrt{R}a^{3}(t_{0})(\partial_{\phi}M_{\rm
eff}^{2})(t)}{\pi\Gamma^{4}(\frac{1}{4})M_{\rm
eff}(t)a^{3}(t)}\Theta(t-t_{0})\quad(\text{for ${\bm{r}=\bm{1}}$}),$ (128)

or

$\displaystyle F_{\phi}(t)\approx\frac{a^{3}(t_{0})(\partial_{\phi}M_{\rm
eff}^{2})(t)}{8\pi^{3}a^{3}(t)M_{\rm
eff}(t)}\left(\frac{9|D|mv}{\sqrt{2}M_{\rm
pl}}\right)^{\frac{3}{2}}\Theta(t-t_{0})\quad(\text{for
${\bm{r}=\bm{1}^{\prime}}$}).$ (129)

We note that if there are $N$ fields that has the same Yukawa coupling, the
effective force would be multiplied by $N$, which enhances the effect. We have
checked the agreement between analytic formulas (127) and numerical results
shown in Fig. 9.

Figure 9: Comparison between the analytic formula (127) (blue lines) and
numerical results (red dots). The left panel shows the ${\bm{r}=\bm{1}}$-case
and the right the ${\bm{r}=\bm{1}^{\prime}}$-case. We see a good agreement for
both cases.

We show a numerical example of the moduli dynamics by using our analytic
approximation in Fig 10. In the example, we have used the
${\bm{1}^{\prime}}$-model. We have used the formulas (146),(147) to choose
appropriate initial conditions such that $\phi$ crosses the ESP at
$t_{0}$.181818Another technical note: We have modified the denominator of
(129) as $M_{\rm eff}\to\sqrt{10^{-6}v+M_{\rm eff}^{2}}$ such that the
singular behavior at the critical point is avoided. We have checked that the
result is not affected by the IR cut-off.  As a cross-check of our numerical
solution, we show the behavior of the field velocity $\dot{\phi}$ in Fig 11.
As is clear from this simulation, the moduli trapping works despite Planck
suppressed couplings between the matter field $\chi$ and the moduli. Thus, the
moduli fields seem to prefer the ESP if they cross such a point along their
time-evolution.

We emphasize that here only the first particle production event is taken into
account, but as we see, there would be secondary and more particle production
events when the modulus crosses the ESP, which further strengthen the trapping
effect as quoted before.

Figure 10: A numerical solution of the E.O.M. of $\phi$ with (left) or without
(right) particle production effect. Here we have taken the following
parameters: The expansion parameter $\beta=\frac{1}{2}$, the initial time
$\sqrt{v}t_{\rm ini}=5$, the time of particle production $\sqrt{v}t_{0}=50$,
the Yukawa coupling parameter $m/\sqrt{v}=100$, and the species number to be
$N=1$. The field velocity at the particle production is
$\sqrt{v}=10^{-4}M_{\rm pl}$. The blue solid curve is the trajectory of $\phi$
and the red dashed line is the critical point ${\rm
Im}\tau=\frac{\sqrt{3}}{2}$.

Figure 11: The time dependence of the field velocity $\dot{\phi}$ with (left)
or without (right) particle production effect. We have taken $v$ as the unit
of the velocity. The parameters are the same as that in Fig. 10. We see that
in both cases the velocity is $1$ at $vt=vt_{0}=50$.

##  5 Summary and discussion

In this work, we have studied the moduli trapping mechanism due to the
particle production near ESPs within modular flavor symmetric models. We have
reviewed and developed a general method for numerical simulations of models
that the classical background fields and quantum fields interact with each
other. Despite generality and simplicity, we have found that such an approach
is not suitable for the application to modular flavor symmetric models because
of the complexity of the couplings between background fields and quantum
fields. Therefore, we have developed a semi-analytic approach where we
analytically evaluate the particle production near the first crossing with the
ESP, which yields the effective potential arising from produced particles.
Although this approach looses the effect of subsequent particle production
events, we have found the expected behavior, namely that the moduli in the
modular flavor symmetric models can be trapped around the ESP at which there
is a residual discrete symmetry.

Although we have studied moduli trapping effects due to the scalar field
$\chi$, we could discuss effects due to spinor and vector fields in a similar
way. As illustrating models, we have used the $A_{4}$ modular forms of weight
8, which have suppressed values around $\tau=\omega$. Similar results would be
obtained for modular forms of generic weights and other finite groups, if they
have suppression behavior around $\tau=\omega$ and $\tau={\rm i}$.
Furthermore, concrete modular flavor symmetric models include many fields,
whose masses are determined by similar modular forms. Thus several modes would
be produced around fixed points, which enlarges the effect of moduli trapping.

There are several issues that should be addressed in future work. One is to
embed this mechanism into more realistic models of particle physics. In
particular, it would be interesting to study the moduli dynamics within
magnetized orbifold models where the standard model flavor structure including
flavor mixings as well mass hierarchies has been realized Abe:2012fj ;
Abe:2014vza ; Fujimoto:2016zjs ; Buchmuller:2017vho ; Buchmuller:2017vut ;
Kikuchi:2021yog ; Hoshiya:2022qvr . It is also important to notice that the
moduli trapping due to particle production never completes the moduli
stabilization since the effective potential disappears as the particle numbers
are diluted by expansion of the Universe, which therefore requires the
introduction of additional potential terms. In particular, we have not taken
into account the 1-loop effective potential such as the Coleman-Weinberg type
potential that we have shown. Although it becomes small in supersymmetric
models, it would be important to investigate the effect of such potential and
the fate of moduli in the late Universe.

Another interesting issue would be the cosmological implications of moduli
dynamics. It has been known that the particle production during inflation may
give imprints on the curvature perturbation spectrum that can be seen from
observations of cosmic microwave background Green:2009ds ; Barnaby:2009mc ;
Pearce:2016qtn ; Flauger:2016idt ; Pearce:2017bdc . In particular, moduli in
modular flavor symmetric models may become inflaton directions. (See e.g.
Kobayashi:2016mzg ; Abe:2023ylh .) Furthermore, the presence of the
dynamically changing CP phase for the standard model matter fields may realize
baryogenesis. In such a case, the method we have applied in this work would be
useful to discuss the dynamics of matter and moduli simultaneously. We leave
these interesting questions and model buildings for future work.

## Acknowledgement

This work was supported by JSPS KAKENHI Grant Numbers JP22KJ0047 (SK) and
JP23K03375 (TK), JST SPRING Grant Number JPMJSP2119 (KN), and Waseda
University Grant for Special Research Projects (Project number: 2023C-584)
(YY).

##  Appendix A Review of particle production in time-dependent backgrounds

We give a short review of particle production in time-dependent backgrounds.
Let us consider a scalar field having a time-dependent mass,

$\displaystyle\hat{\chi}(x)=\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\left[\hat{a}_{\bm{k}}e^{+{\rm
i}{\bm{k}}\cdot{\bm{x}}}f_{k}(t)+\hat{a}_{\bm{k}}^{\dagger}e^{-{\rm
i}{\bm{k}}\cdot{\bm{x}}}f_{k}^{*}(t)\right],$ (130)

where $f(t)$ satisfies a mode equation (5) with $a(t)=1$ and the normalization
condition (7). We have “formally” introduced creation and annihilation
operators, which satisfies the canonical commutation relation (8). This is yet
insufficient to give any meanings to the vacuum state that is annihilated by
$\hat{a}_{\bm{k}}$ since we have not determined $f_{k}(t)$. More precisely
speaking, unless the boundary condition of $f_{k}(t)$ is specified, the above
expansion has no physical meaning. In Minkowski spacetime, the annihilation
operators are introduced as coefficients of positive frequency modes
$\sim\frac{1}{\sqrt{2\omega_{k}}}e^{-{\rm i}\omega_{k}t}$. In time-dependent
backgrounds, it is impossible to define a “global” positive frequency mode,
but locally (in time) it is possible to find an approximate solution. For
instance, the formal solution (9) becomes asymptotically a positive frequency
mode if $\alpha_{k}(t)\to 1$ and $\beta_{k}(t)\to 0$ as $t\to-\infty$.191919We
assume that the background fields asymptote to be time-independent in the
asymptotic past and future. Indeed, under this condition, the scalar operator
becomes

$\displaystyle\hat{\chi}(x)\underset{t\to-\infty}{\to}\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\left[\hat{a}_{\bm{k}}e^{-{\rm
i}\omega_{k}t+{\rm i}{\bm{k}}\cdot{\bm{x}}}+\hat{a}_{\bm{k}}^{\dagger}e^{+{\rm
i}\omega_{k}t-{\rm i}{\bm{k}}\cdot{\bm{x}}}\right],$ (131)

and the vacuum state that satisfies $\hat{a}_{\bm{k}}|0\rangle_{\rm in}=0$ is
understood as a vacuum state in the asymptotic past. Thus we have defined a
past (adiabatic) vacuum.

However, the “past” vacuum state is not a “vacuum” for a future observer. In
general, $\beta_{k}(t)$ becomes non-zero due to Stokes phenomena202020Review
of Stokes phenomena can be found e.g. in Dumlu:2010ua ; Dabrowski:2014ica ;
Dabrowski:2016tsx ; Enomoto:2020xlf ; Taya:2020dco ; Hashiba:2021npn ;
Yamada:2021kqw . For now, we just assume Stokes phenomena occur at some time
and we are considering the time sufficiently far from the event., $t$ being
asymptotic future time. In such a case, it is not appropriate to call
$f_{k}(t)$ as a “positive frequency mode”. Notice that the scalar operator can
be expanded as

$\displaystyle\hat{\chi}(x)$
$\displaystyle\underset{t\to+\infty}{\to}\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\left[(\alpha_{k}(\infty)\hat{a}_{\bm{k}}+\beta_{k}^{*}(\infty)\hat{a}^{\dagger}_{\bm{k}})e^{-{\rm
i}\omega_{k}t+{\rm
i}{\bm{k}}\cdot{\bm{x}}}+(\alpha_{k}^{*}(\infty)\hat{a}_{\bm{k}}^{\dagger}+\beta_{k}(\infty)\hat{a}_{\bm{k}})e^{+{\rm
i}\omega_{k}t-{\rm i}{\bm{k}}\cdot{\bm{x}}}\right]$
$\displaystyle\equiv\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\left[\hat{b}_{\bm{k}}e^{-{\rm
i}\omega_{k}t+{\rm i}{\bm{k}}\cdot{\bm{x}}}+\hat{b}^{\dagger}_{\bm{k}}e^{+{\rm
i}\omega_{k}t-{\rm i}{\bm{k}}\cdot{\bm{x}}}\right].$ (132)

The new set of creation and annihilation operators
$\hat{b}_{\bm{k}},\hat{b}_{\bm{k}}^{\dagger}$ defines a future vacuum state
$|0\rangle_{\rm out}$. Now, we find that the future creation and annihilation
operators are given by linear combinations of the past creation annihilation
operators. Note that there is no momentum exchange and the linear combination
is diagonal with respect to momenta $\bm{k}$.

How can we find particle production from vacuum? We define the future particle
number density as

$\displaystyle\hat{N}_{k}^{f}\equiv\hat{b}_{\bm{k}}^{\dagger}\hat{b}_{\bm{k}}=(\alpha_{k}^{*}(\infty)\hat{a}_{\bm{k}}^{\dagger}+\beta_{k}(\infty)\hat{a}_{\bm{k}})(\alpha_{k}(\infty)\hat{a}_{\bm{k}}+\beta_{k}^{*}(\infty)\hat{a}^{\dagger}_{\bm{k}}).$
(133)

We would like to know how much “future” particles are contained in the past
vacuum state $|0\rangle_{\rm in}$, which can be explicitly evaluated as

$\displaystyle{}_{\rm in}\langle 0|\hat{N}_{k}^{f}|0\rangle_{\rm
in}=|\beta_{k}(\infty)|^{2}\langle\hat{a}_{\bm{k}}\hat{a}_{\bm{k}}^{\dagger}\rangle_{\rm
in}=|\beta_{k}(\infty)|^{2}\frac{V}{(2\pi)^{3}}$ (134)

where we have used
$\lim_{\bm{k}\to\bm{0}}\delta^{3}({\bm{k}})=\frac{V}{(2\pi)^{3}}$ and $V$
being spatial volume. Thus, we have found that the number density for (future)
${\bm{k}}$-mode is given by

$\displaystyle n_{k}=|\beta_{k}(\infty)|^{2}.$ (135)

It turns out that the past “vacuum” state is not a “vacuum” for a future
observer viewpoint. The “particle production from vacuum” would be more
appropriate to be understood as the ambiguity of energy, which does not allow
us to globally define positive and negative modes.

We emphasize that the definition of “particle” is quite ambiguous except for
constant backgrounds. In particular, when time-dependence is not turned off,
we have to introduce adiabatic solutions such as WKB solutions, but there are
infinitely many choices of adiabatic solutions. There is no clear answer to
the question “which solution should we take?” but it is known that the optimal
definition of the “adiabatic particle number” is related to the Stokes
phenomena Dabrowski:2014ica ; Dabrowski:2016tsx .212121The optimal truncation
of the adiabatic series was originally studied in Berry1 ; Berry2 ; Berry3 .

##  Appendix B Modular forms

The modular forms consist of the log derivatives of Dedekind eta function,
which is given by

$\eta(\tau)=q^{\frac{1}{24}}\prod_{n=1}^{\infty}(1-q^{n}),\quad q=e^{2\pi{\rm
i}\tau}.$ (136)

The $A_{4}$ modular forms of weight 2 are Feruglio:2017spp

$\displaystyle Y_{1}(\tau)=\frac{{\rm
i}}{2\pi}\left(\frac{\eta^{\prime}(\tau/3)}{\eta(\tau/3)}+\frac{\eta^{\prime}((\tau+1)/3)}{\eta((\tau+1)/3)}+\frac{\eta^{\prime}((\tau+2)/3)}{\eta((\tau+2)/3)}-\frac{27\eta^{\prime}(3\tau)}{\eta(3\tau)}\right),$
(137) $\displaystyle Y_{2}(\tau)=\frac{-{\rm
i}}{\pi}\left(\frac{\eta^{\prime}(\tau/3)}{\eta(\tau/3)}+\omega^{2}\frac{\eta^{\prime}((\tau+1)/3)}{\eta((\tau+1)/3)}+\omega\frac{\eta^{\prime}((\tau+2)/3)}{\eta((\tau+2)/3)}\right),$
(138) $\displaystyle Y_{3}(\tau)=\frac{-{\rm
i}}{\pi}\left(\frac{\eta^{\prime}(\tau/3)}{\eta(\tau/3)}+\omega\frac{\eta^{\prime}((\tau+1)/3)}{\eta((\tau+1)/3)}+\omega^{2}\frac{\eta^{\prime}((\tau+2)/3)}{\eta((\tau+2)/3)}\right),$
(139)

where $\omega=e^{{\rm i}\frac{2\pi}{3}}$. They correspond to the $A_{4}$
triplet. By noting the fact that any positive integers can be written as
$3m,3m-1,3m-2$, ($m\in\mathbb{N}$) and they consist of the logarithmic
derivatives of the Dedekind eta function, one can easily rewrite these
functions as

$\displaystyle Y_{1}(\tau)=$ $\displaystyle
1+\sum_{m=1}^{\infty}\left[\frac{9mq^{m}-18m^{2m}+9mq^{3m}}{(1-q^{m})^{3}}+\frac{3(3m-2)q^{3m-2}}{1-q^{3m-2}}+\frac{3(3m-1)q^{3m-1}}{1-q^{3m-1}}-\frac{27q^{3m}}{1-q^{3m}}\right],$
(140) $\displaystyle Y_{2}(\tau)=$
$\displaystyle-6q^{-\frac{2}{3}}\sum_{m=1}^{\infty}\left[\frac{(3m-2)q^{m}}{1-q^{3m-2}}+\frac{(3m-1)q^{2}m}{1-q^{3m-1}}\right],$
(141) $\displaystyle Y_{3}(\tau)=$
$\displaystyle-6q^{-\frac{1}{3}}\sum_{m=1}^{\infty}\left[\frac{(3m-1)q^{m}}{(1-q^{3m-1}}+\frac{(3m-2)q^{2m-1}}{1-q^{3m-2}}\right].$
(142)

We will approximate the infinite sum by some (sufficiently large) finite one.
In our numerical simulations, we truncated the series up to some finite order.

By tensor products of $Y_{1}(\tau),Y_{2}(\tau)$, and $Y_{3}(\tau)$, we can
write modular forms of higher weights. The modular forms of weight 8
corresponding to three $A_{4}$ singlets, ${\bf 1},{\bf 1}^{\prime}$, and ${\bf
1}^{\prime\prime}$ can be written by Zhang:2019ngf

$\displaystyle Y^{(8)}_{\bf 1}(\tau)=(Y_{1}^{2}+2Y_{2}Y_{3})^{2},\quad
Y^{(8)}_{{\bf
1}^{\prime}}(\tau)=(Y_{1}^{2}+2Y_{2}Y_{3})(Y_{3}^{2}+2Y_{1}Y_{2}),\quad
Y^{(8)}_{{\bf 1}^{\prime\prime}}(\tau)=(Y_{3}^{2}+2Y_{1}Y_{2})^{2}.$ (143)

##  Appendix C Field dynamics without potential in expanding Universe

Here we discuss the free field dynamics in the expanding Universe. The E.O.M.
of a massless free scalar field in the FRW background is given by

$\ddot{\phi}+3H\dot{\phi}=0,$ (144)

which can also be written as

$\frac{1}{a^{3}}\frac{d}{dt}\left(a^{3}\dot{\phi}\right)=0.$ (145)

Therefore, the first integration yields

$\displaystyle\dot{\phi}(t)=\frac{a^{3}(t_{r})}{a^{3}(t)}\dot{\phi}(t_{r})=\left(\frac{t_{r}}{t}\right)^{3\beta}\dot{\phi}(t_{r}),$
(146)

where $t_{r}$ is a reference time. The secondary time integration leads to

$\displaystyle\phi(t)=\phi(t_{r})+\frac{1}{1-3\beta}\left[\left(\frac{t_{r}}{t}\right)^{3\beta-1}-1\right]t_{r}\dot{\phi}(t_{r}),$
(147)

assuming $\beta\neq\frac{1}{3}$. Near $t=t_{r}$, we may expand the above
expression as

$\displaystyle\phi(t)=\phi(t_{r})+\dot{\phi}(t_{r})(t-t_{r})-\frac{3\beta\dot{\phi}(t_{r})(t-t_{r})^{2}}{2t_{r}}+\mathcal{O}((t-t_{r})^{3}).$
(148)

Now since $H(t_{r})\sim 1/t_{r}$, we can neglect the third and higher order
terms as long as Hubble parameter at $t=t_{r}$ is sufficiently smaller than
$\sqrt{|\dot{\phi}(t_{r})|}$.

The relation between the initial velocity $v_{\rm ini}=\dot{\phi}(t_{\rm
ini})$ and $v=\dot{\phi}(t_{0})$ at some reference time $t_{0}>t_{\rm ini}$ is

$v=\left(\frac{t_{\rm ini}}{t_{0}}\right)^{3\beta}v_{\rm ini}.$ (149)

This would be useful to relate the initial velocity with the field velocity at
the crossing with enhanced symmetry points particularly in estimating the
particle production rate.

## References

  * (1) S. Ferrara, D. Lust, A. D. Shapere and S. Theisen, _Modular Invariance in Supersymmetric Field Theories_ , _Phys. Lett. B_ 225 (1989) 363.
  * (2) S. Ferrara, . D. Lust and S. Theisen, _Target Space Modular Invariance and Low-Energy Couplings in Orbifold Compactifications_ , _Phys. Lett. B_ 233 (1989) 147–152.
  * (3) W. Lerche, D. Lust and N. P. Warner, _Duality Symmetries in $N=2$ Landau-ginzburg Models_, _Phys. Lett. B_ 231 (1989) 417–424.
  * (4) J. Lauer, J. Mas and H. P. Nilles, _Duality and the Role of Nonperturbative Effects on the World Sheet_ , _Phys. Lett. B_ 226 (1989) 251–256.
  * (5) J. Lauer, J. Mas and H. P. Nilles, _Twisted sector representations of discrete background symmetries for two-dimensional orbifolds_ , _Nucl. Phys. B_ 351 (1991) 353–424.
  * (6) T. Kobayashi, S. Nagamoto, S. Takada, S. Tamba and T. H. Tatsuishi, _Modular symmetry and non-Abelian discrete flavor symmetries in string compactification_ , _Phys. Rev. D_ 97 (2018) 116002, [1804.06644].
  * (7) T. Kobayashi and S. Tamba, _Modular forms of finite modular subgroups from magnetized D-brane models_ , _Phys. Rev. D_ 99 (2019) 046001, [1811.11384].
  * (8) H. Ohki, S. Uemura and R. Watanabe, _Modular flavor symmetry on a magnetized torus_ , _Phys. Rev. D_ 102 (2020) 085008, [2003.04174].
  * (9) S. Kikuchi, T. Kobayashi, S. Takada, T. H. Tatsuishi and H. Uchida, _Revisiting modular symmetry in magnetized torus and orbifold compactifications_ , _Phys. Rev. D_ 102 (2020) 105010, [2005.12642].
  * (10) S. Kikuchi, T. Kobayashi, H. Otsuka, S. Takada and H. Uchida, _Modular symmetry by orbifolding magnetized $T^{2}\times T^{2}$: realization of double cover of $\Gamma_{N}$_, _JHEP_ 11 (2020) 101, [2007.06188].
  * (11) S. Kikuchi, T. Kobayashi and H. Uchida, _Modular flavor symmetries of three-generation modes on magnetized toroidal orbifolds_ , _Phys. Rev. D_ 104 (2021) 065008, [2101.00826].
  * (12) Y. Almumin, M.-C. Chen, V. Knapp-Pérez, S. Ramos-Sánchez, M. Ratz and S. Shukla, _Metaplectic Flavor Symmetries from Magnetized Tori_ , _JHEP_ 05 (2021) 078, [2102.11286].
  * (13) A. Baur, H. P. Nilles, A. Trautner and P. K. S. Vaudrevange, _Unification of Flavor, CP, and Modular Symmetries_ , _Phys. Lett. B_ 795 (2019) 7–14, [1901.03251].
  * (14) A. Baur, H. P. Nilles, A. Trautner and P. K. S. Vaudrevange, _A String Theory of Flavor and CP_ , _Nucl. Phys. B_ 947 (2019) 114737, [1908.00805].
  * (15) H. P. Nilles, S. Ramos-Sanchez and P. K. S. Vaudrevange, _Lessons from eclectic flavor symmetries_ , _Nucl. Phys. B_ 957 (2020) 115098, [2004.05200].
  * (16) H. P. Nilles, S. Ramos–Sánchez and P. K. S. Vaudrevange, _Eclectic flavor scheme from ten-dimensional string theory - II detailed technical analysis_ , _Nucl. Phys. B_ 966 (2021) 115367, [2010.13798].
  * (17) R. de Adelhart Toorop, F. Feruglio and C. Hagedorn, _Finite Modular Groups and Lepton Mixing_ , _Nucl. Phys. B_ 858 (2012) 437–467, [1112.1340].
  * (18) G. Altarelli and F. Feruglio, _Discrete Flavor Symmetries and Models of Neutrino Mixing_ , _Rev. Mod. Phys._ 82 (2010) 2701–2729, [1002.0211].
  * (19) H. Ishimori, T. Kobayashi, H. Ohki, Y. Shimizu, H. Okada and M. Tanimoto, _Non-Abelian Discrete Symmetries in Particle Physics_ , _Prog. Theor. Phys. Suppl._ 183 (2010) 1–163, [1003.3552].
  * (20) T. Kobayashi, H. Ohki, H. Okada, Y. Shimizu and M. Tanimoto, _An Introduction to Non-Abelian Discrete Symmetries for Particle Physicists_.  1, 2022, 10.1007/978-3-662-64679-3.
  * (21) D. Hernandez and A. Y. Smirnov, _Lepton mixing and discrete symmetries_ , _Phys. Rev. D_ 86 (2012) 053014, [1204.0445].
  * (22) S. F. King and C. Luhn, _Neutrino Mass and Mixing with Discrete Symmetry_ , _Rept. Prog. Phys._ 76 (2013) 056201, [1301.1340].
  * (23) S. F. King, A. Merle, S. Morisi, Y. Shimizu and M. Tanimoto, _Neutrino Mass and Mixing: from Theory to Experiment_ , _New J. Phys._ 16 (2014) 045018, [1402.4271].
  * (24) M. Tanimoto, _Neutrinos and flavor symmetries_ , _AIP Conf. Proc._ 1666 (2015) 120002.
  * (25) S. F. King, _Unified Models of Neutrinos, Flavour and CP Violation_ , _Prog. Part. Nucl. Phys._ 94 (2017) 217–256, [1701.04413].
  * (26) S. T. Petcov, _Discrete Flavour Symmetries, Neutrino Mixing and Leptonic CP Violation_ , _Eur. Phys. J. C_ 78 (2018) 709, [1711.10806].
  * (27) F. Feruglio, _Are neutrino masses modular forms?_ , pp. 227–266.  2019\.  1706.08749.
  * (28) T. Kobayashi, K. Tanaka and T. H. Tatsuishi, _Neutrino mixing from finite modular groups_ , _Phys. Rev. D_ 98 (2018) 016004, [1803.10391].
  * (29) J. T. Penedo and S. T. Petcov, _Lepton Masses and Mixing from Modular $S_{4}$ Symmetry_, _Nucl. Phys. B_ 939 (2019) 292–307, [1806.11040].
  * (30) J. C. Criado and F. Feruglio, _Modular Invariance Faces Precision Neutrino Data_ , _SciPost Phys._ 5 (2018) 042, [1807.01125].
  * (31) T. Kobayashi, N. Omoto, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, _Modular A 4 invariance and neutrino mixing_, _JHEP_ 11 (2018) 196, [1808.03012].
  * (32) P. P. Novichkov, J. T. Penedo, S. T. Petcov and A. V. Titov, _Modular S 4 models of lepton masses and mixing_, _JHEP_ 04 (2019) 005, [1811.04933].
  * (33) P. P. Novichkov, J. T. Penedo, S. T. Petcov and A. V. Titov, _Modular A 5 symmetry for flavour model building_, _JHEP_ 04 (2019) 174, [1812.02158].
  * (34) F. J. de Anda, S. F. King and E. Perdomo, _$SU(5)$ grand unified theory with $A_{4}$ modular symmetry_, _Phys. Rev. D_ 101 (2020) 015028, [1812.05620].
  * (35) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto, T. H. Tatsuishi and H. Uchida, _Finite modular subgroups for fermion mass matrices and baryon/lepton number violation_ , _Phys. Lett. B_ 794 (2019) 114–121, [1812.11072].
  * (36) H. Okada and M. Tanimoto, _CP violation of quarks in $A_{4}$ modular invariance_, _Phys. Lett. B_ 791 (2019) 54–61, [1812.09677].
  * (37) P. P. Novichkov, S. T. Petcov and M. Tanimoto, _Trimaximal Neutrino Mixing from Modular A4 Invariance with Residual Symmetries_ , _Phys. Lett. B_ 793 (2019) 247–258, [1812.11289].
  * (38) G.-J. Ding, S. F. King and X.-G. Liu, _Neutrino mass and mixing with $A_{5}$ modular symmetry_, _Phys. Rev. D_ 100 (2019) 115005, [1903.12588].
  * (39) T. Nomura and H. Okada, _A modular $A_{4}$ symmetric model of dark matter and neutrino_, _Phys. Lett. B_ 797 (2019) 134799, [1904.03937].
  * (40) P. P. Novichkov, J. T. Penedo, S. T. Petcov and A. V. Titov, _Generalised CP Symmetry in Modular-Invariant Models of Flavour_ , _JHEP_ 07 (2019) 165, [1905.11970].
  * (41) H. Okada and M. Tanimoto, _Towards unification of quark and lepton flavors in $A_{4}$ modular invariance_, _Eur. Phys. J. C_ 81 (2021) 52, [1905.13421].
  * (42) I. de Medeiros Varzielas, S. F. King and Y.-L. Zhou, _Multiple modular symmetries as the origin of flavor_ , _Phys. Rev. D_ 101 (2020) 055033, [1906.02208].
  * (43) T. Nomura and H. Okada, _A two loop induced neutrino mass model with modular $A_{4}$ symmetry_, _Nucl. Phys. B_ 966 (2021) 115372, [1906.03927].
  * (44) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, _Modular $S_{3}$-invariant flavor model in SU(5) grand unified theory_, _PTEP_ 2020 (2020) 053B05, [1906.10341].
  * (45) X.-G. Liu and G.-J. Ding, _Neutrino Masses and Mixing from Double Covering of Finite Modular Groups_ , _JHEP_ 08 (2019) 134, [1907.01488].
  * (46) H. Okada and Y. Orikasa, _Modular $S_{3}$ symmetric radiative seesaw model_, _Phys. Rev. D_ 100 (2019) 115037, [1907.04716].
  * (47) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, _New $A_{4}$ lepton flavor model from $S_{4}$ modular symmetry_, _JHEP_ 02 (2020) 097, [1907.09141].
  * (48) G.-J. Ding, S. F. King and X.-G. Liu, _Modular A 4 symmetry models of neutrinos and charged leptons_, _JHEP_ 09 (2019) 074, [1907.11714].
  * (49) S. F. King and Y.-L. Zhou, _Trimaximal TM 1 mixing with two modular $S_{4}$ groups_, _Phys. Rev. D_ 101 (2020) 015001, [1908.02770].
  * (50) T. Nomura, H. Okada and O. Popov, _A modular $A_{4}$ symmetric scotogenic model_, _Phys. Lett. B_ 803 (2020) 135294, [1908.07457].
  * (51) J. C. Criado, F. Feruglio and S. J. D. King, _Modular Invariant Models of Lepton Masses at Levels 4 and 5_ , _JHEP_ 02 (2020) 001, [1908.11867].
  * (52) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, _$A_{4}$ lepton flavor model and modulus stabilization from $S_{4}$ modular symmetry_, _Phys. Rev. D_ 100 (2019) 115045, [1909.05139].
  * (53) T. Asaka, Y. Heo, T. H. Tatsuishi and T. Yoshida, _Modular $A_{4}$ invariance and leptogenesis_, _JHEP_ 01 (2020) 144, [1909.06520].
  * (54) G.-J. Ding, S. F. King, X.-G. Liu and J.-N. Lu, _Modular S 4 and A4 symmetries and their fixed points: new predictive examples of lepton mixing_, _JHEP_ 12 (2019) 030, [1910.03460].
  * (55) D. Zhang, _A modular $A_{4}$ symmetry realization of two-zero textures of the Majorana neutrino mass matrix_, _Nucl. Phys. B_ 952 (2020) 114935, [1910.07869].
  * (56) X. Wang and S. Zhou, _The minimal seesaw model with a modular S 4 symmetry_, _JHEP_ 05 (2020) 017, [1910.09473].
  * (57) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto, T. H. Tatsuishi and H. Uchida, _$CP$ violation in modular invariant flavor models_, _Phys. Rev. D_ 101 (2020) 055046, [1910.11553].
  * (58) T. Nomura, H. Okada and S. Patra, _An inverse seesaw model with $A_{4}$ -modular symmetry_, _Nucl. Phys. B_ 967 (2021) 115395, [1912.00379].
  * (59) T. Kobayashi, T. Nomura and T. Shimomura, _Type II seesaw models with modular $A_{4}$ symmetry_, _Phys. Rev. D_ 102 (2020) 035019, [1912.00637].
  * (60) J.-N. Lu, X.-G. Liu and G.-J. Ding, _Modular symmetry origin of texture zeros and quark lepton unification_ , _Phys. Rev. D_ 101 (2020) 115020, [1912.07573].
  * (61) X. Wang, _Lepton flavor mixing and CP violation in the minimal type-(I+II) seesaw model with a modular $A_{4}$ symmetry_, _Nucl. Phys. B_ 957 (2020) 115105, [1912.13284].
  * (62) S. J. D. King and S. F. King, _Fermion mass hierarchies from modular symmetry_ , _JHEP_ 09 (2020) 043, [2002.00969].
  * (63) M. Abbas, _Fermion masses and mixing in modular A4 Symmetry_ , _Phys. Rev. D_ 103 (2021) 056016, [2002.01929].
  * (64) P. Chen, G.-J. Ding, J.-N. Lu and J. W. F. Valle, _Predictions from warped flavor dynamics based on the $T^{\prime}$ family group_, _Phys. Rev. D_ 102 (2020) 095014, [2003.02734].
  * (65) H. Okada and Y. Shoji, _Dirac dark matter in a radiative neutrino model_ , _Phys. Dark Univ._ 31 (2021) 100742, [2003.11396].
  * (66) H. Okada and Y. Shoji, _A radiative seesaw model with three Higgs doublets in modular $A_{4}$ symmetry_, _Nucl. Phys. B_ 961 (2020) 115216, [2003.13219].
  * (67) G.-J. Ding and F. Feruglio, _Testing Moduli and Flavon Dynamics with Neutrino Oscillations_ , _JHEP_ 06 (2020) 134, [2003.13448].
  * (68) G.-J. Ding, S. F. King, C.-C. Li and Y.-L. Zhou, _Modular Invariant Models of Leptons at Level 7_ , _JHEP_ 08 (2020) 164, [2004.12662].
  * (69) H. Okada and M. Tanimoto, _Quark and lepton flavors with common modulus $\tau$ in A4 modular symmetry_, _Phys. Dark Univ._ 40 (2023) 101204, [2005.00775].
  * (70) P. P. Novichkov, J. T. Penedo and S. T. Petcov, _Double cover of modular $S_{4}$ for flavour model building_, _Nucl. Phys. B_ 963 (2021) 115301, [2006.03058].
  * (71) X.-G. Liu, C.-Y. Yao and G.-J. Ding, _Modular invariant quark and lepton models in double covering of $S_{4}$ modular group_, _Phys. Rev. D_ 103 (2021) 056013, [2006.10722].
  * (72) I. de Medeiros Varzielas, M. Levy and Y.-L. Zhou, _Symmetries and stabilisers in modular invariant flavour models_ , _JHEP_ 11 (2020) 085, [2008.05329].
  * (73) T. Asaka, Y. Heo and T. Yoshida, _Lepton flavor model with modular $A_{4}$ symmetry in large volume limit_, _Phys. Lett. B_ 811 (2020) 135956, [2009.12120].
  * (74) H. Okada and M. Tanimoto, _Modular invariant flavor model of $A_{4}$ and hierarchical structures at nearby fixed points_, _Phys. Rev. D_ 103 (2021) 015005, [2009.14242].
  * (75) X. Wang, B. Yu and S. Zhou, _Double covering of the modular $A_{5}$ group and lepton flavor mixing in the minimal seesaw model_, _Phys. Rev. D_ 103 (2021) 076005, [2010.10159].
  * (76) H. Okada and M. Tanimoto, _Spontaneous CP violation by modulus $\tau$ in $A_{4}$ model of lepton flavors_, _JHEP_ 03 (2021) 010, [2012.01688].
  * (77) C.-Y. Yao, J.-N. Lu and G.-J. Ding, _Modular Invariant $A_{4}$ Models for Quarks and Leptons with Generalized CP Symmetry_, _JHEP_ 05 (2021) 102, [2012.13390].
  * (78) C.-C. Li, X.-G. Liu and G.-J. Ding, _Modular symmetry at level 6 and a new route towards finite modular groups_ , _JHEP_ 10 (2021) 238, [2108.02181].
  * (79) T. Kobayashi and M. Tanimoto, _Modular flavor symmetric models_ , 7, 2023\.  2307.03384.
  * (80) H. Abe, T. Kobayashi, S. Uemura and J. Yamamoto, _Loop Fayet-Iliopoulos terms in $T^{2}/Z_{2}$ models: Instability and moduli stabilization_, _Phys. Rev. D_ 102 (2020) 045005, [2003.03512].
  * (81) K. Ishiguro, T. Kobayashi and H. Otsuka, _Landscape of Modular Symmetric Flavor Models_ , _JHEP_ 03 (2021) 161, [2011.09154].
  * (82) P. P. Novichkov, J. T. Penedo and S. T. Petcov, _Modular flavour symmetries and modulus stabilisation_ , _JHEP_ 03 (2022) 149, [2201.02020].
  * (83) K. Ishiguro, H. Okada and H. Otsuka, _Residual flavor symmetry breaking in the landscape of modular flavor models_ , _JHEP_ 09 (2022) 072, [2206.04313].
  * (84) J. M. Leedom, N. Righi and A. Westphal, _Heterotic de Sitter beyond modular symmetry_ , _JHEP_ 02 (2023) 209, [2212.03876].
  * (85) V. Knapp-Perez, X.-G. Liu, H. P. Nilles, S. Ramos-Sanchez and M. Ratz, _Matter matters in moduli fixing and modular flavor symmetries_ , 2304.14437.
  * (86) F. Feruglio, V. Gherardi, A. Romanino and A. Titov, _Modular invariant dynamics and fermion mass hierarchies around $\tau=i$_, _JHEP_ 05 (2021) 242, [2101.08718].
  * (87) P. P. Novichkov, J. T. Penedo and S. T. Petcov, _Fermion mass hierarchies, large lepton mixing and residual modular symmetries_ , _JHEP_ 04 (2021) 206, [2102.07488].
  * (88) S. T. Petcov and M. Tanimoto, _$A_{4}$ modular flavour model of quark mass hierarchies close to the fixed point $\tau=\omega$_, _Eur. Phys. J. C_ 83 (2023) 579, [2212.13336].
  * (89) S. Kikuchi, T. Kobayashi, K. Nasu, S. Takada and H. Uchida, _Quark hierarchical structures in modular symmetric flavor models at level 6_ , _Phys. Rev. D_ 107 (2023) 055014, [2301.03737].
  * (90) Y. Abe, T. Higaki, J. Kawamura and T. Kobayashi, _Quark masses and CKM hierarchies from $S_{4}^{\prime}$ modular flavor symmetry_, 2301.07439.
  * (91) S. Kikuchi, T. Kobayashi, K. Nasu, S. Takada and H. Uchida, _Quark mass hierarchies and CP violation in A 4 × A4 × A4 modular symmetric flavor models_, _JHEP_ 07 (2023) 134, [2302.03326].
  * (92) Y. Abe, T. Higaki, J. Kawamura and T. Kobayashi, _Quark and lepton hierarchies from S4’ modular flavor symmetry_ , _Phys. Lett. B_ 842 (2023) 137977, [2302.11183].
  * (93) S. T. Petcov and M. Tanimoto, _$A_{4}$ Modular Flavour Model of Quark Mass Hierarchies close to the Fixed Point $\tau=i\infty$_, 2306.05730.
  * (94) Y. Abe, T. Higaki, J. Kawamura and T. Kobayashi, _Fermion Hierarchies in $SU(5)$ Grand Unification from $\Gamma_{6}^{\prime}$ Modular Flavor Symmetry_, 2307.01419.
  * (95) L. Kofman, A. D. Linde, X. Liu, A. Maloney, L. McAllister and E. Silverstein, _Beauty is attractive: Moduli trapping at enhanced symmetry points_ , _JHEP_ 05 (2004) 030, [hep-th/0403001].
  * (96) S. Habib, C. Molina-Paris and E. Mottola, _Energy momentum tensor of particles created in an expanding universe_ , _Phys. Rev. D_ 61 (2000) 024010, [gr-qc/9906120].
  * (97) R. Dabrowski and G. V. Dunne, _Superadiabatic particle number in Schwinger and de Sitter particle production_ , _Phys. Rev. D_ 90 (2014) 025021, [1405.0302].
  * (98) R. Dabrowski and G. V. Dunne, _Time dependence of adiabatic particle number_ , _Phys. Rev. D_ 94 (2016) 065005, [1606.00902].
  * (99) P. B. Greene and L. Kofman, _On the theory of fermionic preheating_ , _Phys. Rev. D_ 62 (2000) 123516, [hep-ph/0003018].
  * (100) M. Peloso and L. Sorbo, _Preheating of massive fermions after inflation: Analytical results_ , _JHEP_ 05 (2000) 016, [hep-ph/0003045].
  * (101) P. Adshead and E. I. Sfakianakis, _Fermion production during and after axion inflation_ , _JCAP_ 11 (2015) 021, [1508.00891].
  * (102) C. K. Dumlu and G. V. Dunne, _The Stokes Phenomenon and Schwinger Vacuum Pair Production in Time-Dependent Laser Pulses_ , _Phys. Rev. Lett._ 104 (2010) 250402, [1004.2509].
  * (103) S. Enomoto and T. Matsuda, _The exact WKB for cosmological particle production_ , _JHEP_ 03 (2021) 090, [2010.14835].
  * (104) H. Taya, T. Fujimori, T. Misumi, M. Nitta and N. Sakai, _Exact WKB analysis of the vacuum pair production by time-dependent electric fields_ , _JHEP_ 03 (2021) 082, [2010.16080].
  * (105) S. Hashiba and Y. Yamada, _Stokes phenomenon and gravitational particle production– How to evaluate it in practice_ , _JCAP_ 05 (2021) 022, [2101.07634].
  * (106) Y. Yamada, _Superadiabatic basis in cosmological particle production: application to preheating_ , _JCAP_ 09 (2021) 009, [2106.06111].
  * (107) N. Fröman and O. Dammert, _Tunneling and super-barrier transmission through a system of two real potential barriers_ , _Nucl. Phys. A_ 147 (1970) 627–649.
  * (108) H. Abe, T. Kobayashi, H. Ohki, A. Oikawa and K. Sumita, _Phenomenological aspects of 10D SYM theory with magnetized extra dimensions_ , _Nucl. Phys. B_ 870 (2013) 30–54, [1211.4317].
  * (109) H. Abe, T. Kobayashi, K. Sumita and Y. Tatsuta, _Gaussian Froggatt-Nielsen mechanism on magnetized orbifolds_ , _Phys. Rev. D_ 90 (2014) 105006, [1405.5012].
  * (110) Y. Fujimoto, T. Kobayashi, K. Nishiwaki, M. Sakamoto and Y. Tatsuta, _Comprehensive analysis of Yukawa hierarchies on $T^{2}/Z_{N}$ with magnetic fluxes_, _Phys. Rev. D_ 94 (2016) 035031, [1605.00140].
  * (111) W. Buchmuller and J. Schweizer, _Flavor mixings in flux compactifications_ , _Phys. Rev. D_ 95 (2017) 075024, [1701.06935].
  * (112) W. Buchmuller and K. M. Patel, _Flavor physics without flavor symmetries_ , _Phys. Rev. D_ 97 (2018) 075019, [1712.06862].
  * (113) S. Kikuchi, T. Kobayashi, Y. Ogawa and H. Uchida, _Yukawa textures in modular symmetric vacuum of magnetized orbifold models_ , _PTEP_ 2022 (2022) 033B10, [2112.01680].
  * (114) K. Hoshiya, S. Kikuchi, T. Kobayashi and H. Uchida, _Quark and lepton flavor structure in magnetized orbifold models at residual modular symmetric points_ , _Phys. Rev. D_ 106 (2022) 115003, [2209.07249].
  * (115) D. Green, B. Horn, L. Senatore and E. Silverstein, _Trapped Inflation_ , _Phys. Rev. D_ 80 (2009) 063533, [0902.1006].
  * (116) N. Barnaby, Z. Huang, L. Kofman and D. Pogosyan, _Cosmological Fluctuations from Infra-Red Cascading During Inflation_ , _Phys. Rev. D_ 80 (2009) 043501, [0902.0615].
  * (117) L. Pearce, M. Peloso and L. Sorbo, _The phenomenology of trapped inflation_ , _JCAP_ 11 (2016) 058, [1603.08021].
  * (118) R. Flauger, M. Mirbabayi, L. Senatore and E. Silverstein, _Productive Interactions: heavy particles and non-Gaussianity_ , _JCAP_ 10 (2017) 058, [1606.00513].
  * (119) L. Pearce, M. Peloso and L. Sorbo, _Resonant particle production during inflation: a full analytical study_ , _JCAP_ 05 (2017) 054, [1702.07661].
  * (120) T. Kobayashi, D. Nitta and Y. Urakawa, _Modular invariant inflation_ , _JCAP_ 08 (2016) 014, [1604.02995].
  * (121) Y. Abe, T. Higaki, F. Kaneko, T. Kobayashi and H. Otsuka, _Moduli inflation from modular flavor symmetries_ , _JHEP_ 06 (2023) 187, [2303.02947].
  * (122) M. V. Berry, _Stokes’ phenomenon; smoothing a victorian discontinuity_ , _Publications Mathématiques de l’Institut des Hautes Scientifiques_ 68 (1988) 211–221.
  * (123) M. V. Berry, _Waves near Stokes lines_ , _Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences_ 427 (1990) 265–280.
  * (124) R. Lim and M. V. Berry, _Superadiabatic tracking of quantum evolution_ , _Journal of Physics A: Mathematical and General_ 24 (1990) 3255–3264.