\section{Introduction}
$CP$ symmetry is a fundamental discrete symmetry which plays an important role
in particle physics. In the standard model (SM), it is considered to be violated
explicitly through complex Yukawa coupling constants \cite{km} and a $\theta$
parameter in the QCD sector \cite{theta}. The former is known to explain very well
$CP$ violating phenomena in $B$ meson systems and so on \cite{pdg}.
The latter is severely constrained through an experimental search of
a neutron electric dipole moment \cite{ediploe} and causes the notorious
strong $CP$ problem \cite{strongcp}.
Peccei-Quinn (PQ) symmetry has been proposed to solve it \cite{pq}.
If we assume that the $CP$ symmetry is an original symmetry of nature,
the complex phases of Yukawa coupling constants have to be spontaneously
induced through some mechanism at high energy regions.
It may be compactification dynamics in string theory near at the Planck
scale \cite{comp-cp}.
In that case, since nonzero $\theta$ could be caused through radiative effects after
the $CP$ violation, the PQ symmetry is required to solve the strong $CP$
problem again.

As an alternative scenario for the realization of $CP$ symmetry,
we may consider it to be exact so that all coupling constants
including Yukawa couplings are real and also $\theta=0$ is kept
to some scale much lower than the Planck scale.
In that case, the $CP$ symmetry is supposed to be spontaneously broken and
this violation can be expected to be transformed to a complex phase in the CKM
matrix effectively.
If nonzero $\theta$ is not brought about in this process, it is favorable for
the strong $CP$ problem. The Nelson-Barr (NB) mechanism \cite{nb} has
been proposed as such a concrete example.
Unfortunately, radiative effects could cause nonzero $\theta$ with a magnitude
which contradicts the experimental constraints \cite{nb-theta}.\footnote{
Introduction of the PQ symmetry could solve this fault of the model.
We consider such a possibility in the extension of the model.}
However, the scenario is interesting since it can present an explanation for
the origin of the $CP$ violation at a much lower energy region than the
Planck scale. As a realization of the NB mechanism, a simple model has been
proposed in \cite{bbp}. The model is extended in \cite{patisalam} to the
lepton sector where the existence of a $CP$ violating phase in the
PMNS matrix \cite{pmns} is suggested through recent neutrino oscillation
experiments \cite{lepcp}.

Observations of the CMB fluctuation \cite{cmb,planck18} suggest the
existence of the exponential expansion of the universe called inflation.
Inflation is usually considered to be induced by some slowly rolling scalar field
called inflaton \cite{inf}.
It is a crucial problem to identify its candidate from a viewpoint of the extension
of the SM. Although the Higgs scalar has been studied as an only promising
candidate in the SM \cite{higgsinf} under an assumption that
it has a nonminimal coupling with the Ricci scalar curvature \cite{nonminimal},
several problems have been pointed out \cite{stab,unitarity,spike}.
In this situation, it is interesting to find an alternative candidate for inflaton
in a certain extension of the SM which could solve several problems in the SM.
In this sense, the model extended from a viewpoint of the $CP$
symmetry as described above could give such a promising candidate.
It contains singlet scalars which cause the spontaneous $CP$ violation
and allow the introduction of the PQ symmetry as a solution
for the strong $CP$ problem.
If they couple with the scalar curvature nonminimally, it could cause
slow-roll inflation successfully.
In this paper, we discuss such a possibility that the inflation of the universe
could be related to the $CP$ violation in the SM.
We study reheating and leptogenesis as its phenomenology caused
by extra fields introduced in the model to solve the $CP$ issues.

Remaining parts of the paper are organized as follows.
In section 2, we describe the model studied in this paper and discuss
both phases in the CKM and PMNS matrices which are derived as a result
of the spontaneous $CP$ violation. In section 3, we discuss the inflation brought
about by the singlet scalars which are related to the $CP$ issues and
the reheating. After that, we describe leptogenesis
which could show a distinguishable feature from the usual leptogenesis in the
seesaw model.  The paper is summarized in section 4.

\section{Origin of $CP$ violation}
\subsection{An extended model}
Our model is an $CP$ invariant extension of the SM with global $U(1)\times Z_4$
symmetry and several additional fields.
As fermions, we introduce
 a pair of vector-like down-type quarks $(D_L, D_R)$, a pair of vector-like charged
leptons $(E_L, E_R)$, and three right-handed singlet fermions
$N_j~(j=1,2,3)$.\footnote{Similar models with vector-like extra fermions have
been considered under different symmetry structures \cite{patisalam,mksvz}. }
We also introduce an additional doublet scalar $\eta$ and two singlet complex
scalars $\sigma$ and $S$. Their representation and charge under
the symmetry $[SU(3)_C\times SU(2)_L\times U(1)_Y]\times U(1)\times Z_4$
are summarized in Table 1.\footnote{
$Z_4$ is imposed by hand to control the couplings of the new fields
to the SM contents.}
The SM contents are assumed to have no charge of the global symmetry.
Since this global $U(1)$ has color anomaly in the same way as the KSVZ
model \cite{ksvz} for a strong $CP$ problem,
it can play the role of the PQ symmetry.
The present charge assignment for colored fermions guarantees the domain wall
number to be one ($N_{\rm DW}=1$)
so that the model can escape the domain wall problem \cite{dw,stdw}.

\begin{figure}[t]
\begin{center}
\small
\begin{tabular}{c|ccccc||c|ccccc}\hline
& $SU(3)_C$ & $SU(2)_L$ & $U(1)_Y$ & $U(1)$ & $Z_4$ & &
$SU(3)_C$ & $SU(2)_L$ & $U(1)_Y$ & $U(1)$ & $Z_4$ \\\hline
$D_L$ & {\bf 3} & {\bf 1} & $-\frac{1}{3}$ & 0 & 2 &
$D_R$ & {\bf 3} & {\bf 1} & $-\frac{1}{3}$ & 2 & 0 \\
$E_L$ & {\bf 1} & {\bf 1} & $-1$ & 0 & 2 &
$E_R$ & {\bf 1} & {\bf 1} & $-1$ & 2 & 0 \\
$\sigma$ & {\bf 1} & {\bf 1} & 0 & $-2$ & 2 &
$S$ & {\bf 1} & {\bf 1} & $0$ & 0 & 2 \\
$N_k$ & {\bf 1} & {\bf 1} & 0 & $1$ & 1 &
$\eta$ & {\bf 1} & {\bf 2} & $-\frac{1}{2}$ & $-1$ &$-1$ \\\hline
\end{tabular}
\end{center}
{\footnotesize {\bf Table 1}~~New fields added to the SM and their
representation and charge under
$[SU(3)_C\times SU(2)_L\times U(1)_Y]\times U(1)\times Z_4$.}
\end{figure}
\normalsize

The model is characterized by new Yukawa terms and
scalar potential which are invariant under the imposed symmetry
\begin{eqnarray}
-{\cal L}_Y&=&y_D\sigma\bar D_LD_R +y_E\sigma\bar E_LE_R+
\sum_{j=1}^3\Big(\frac{y_{N_j}}{2}\sigma \bar N_j^cN_j+ y_{d_j}S\bar D_L d_{R_j}
+ \tilde y_{d_j}S^\dagger\bar D_L d_{R_j} \nonumber\\
&+&y_{e_j}S\bar E_L e_{R_j} +\tilde y_{e_j}S^\dagger\bar E_L e_{R_j} +
\sum_{\alpha=1}^3 h_{\alpha j}^\ast\eta \bar\ell_\alpha N_j \Big)
+\sum_{\alpha,\beta=1}^3
\frac{y_{\alpha\beta}}{M_\ast}(\bar\ell_\alpha\phi)(\bar\ell_\beta\phi)
+ {\rm h.c.}, \nonumber \\
V&=&\lambda_1(\phi^\dagger\phi)^2+\lambda_2(\eta^\dagger\eta)^2
+\lambda_3(\phi^\dagger\phi)(\eta^\dagger\eta)+
\lambda_4(\phi^\dagger\eta)(\eta^\dagger\phi)+
\frac{\lambda_5}{2M_\ast}[\sigma(\eta^\dagger\phi)^2+{\rm h.c.}] \nonumber \\
&+&\kappa_\sigma(\sigma^\dagger\sigma)^2+\kappa_S(S^\dagger S)^2
+(\kappa_{\phi\sigma}\phi^\dagger\phi
+\kappa_{\eta\sigma}\eta^\dagger\eta)(\sigma^\dagger\sigma)
+(\kappa_{\phi S}\phi^\dagger\phi+\kappa_{\eta S}\eta^\dagger\eta) (S^\dagger S)
\nonumber \\
&+&\kappa_{\sigma S}(\sigma^\dagger\sigma)(S^\dagger S)
+ m_\phi^2\phi^\dagger\phi+ m_\eta^2\eta^\dagger\eta
+m_\sigma^2\sigma^\dagger\sigma+m_S^2S^\dagger S +V_b,
\label{model}
\end{eqnarray}
where $d_{R_j}$ and $e_{R_j}$ are the SM down-type quarks and charged leptons, respectively.
$\ell_\alpha$ is a doublet lepton and $\phi$ is an ordinary doublet Higgs scalar.
Since $CP$ invariance is assumed, parameters in Lagrangian are
considered to be all real.
In eq.~(\ref{model}), we list dominant terms up to dimension five and
$M_\ast$ is a cut-off scale of the model.
Other invariant terms are higher order and can be safely
neglected in comparison with the listed ones.
$V_b$ is composed of terms which are invariant under the global symmetry but
violate the $S$ number.

For a while, we focus on a part of field space where the field values of $\sigma$
and $S$ are much larger than both $\phi$ and $\eta$ to study the potential
composed of $\sigma$ and $S$ only.
In the present study, we assume that $V_b$ takes a following form:\footnote{
Imposed symmetry allows terms
$m_S^{\prime 2}(S^2+S^{\dagger 2})$ in $V_b$. However, we assume that their contribution
is negligible since we focus our attention on a potential valley where
$\tilde\sigma\propto \tilde S$ is satisfied and then $\cos 2\rho$ could be a constant
at the minimum of $V_b$ in that case.}
\begin{equation}
V_b=\alpha(S^4+S^{\dagger 4})+\beta\sigma^\dagger\sigma(S^2+S^{\dagger 2})
=\frac{1}{2}\tilde S^2(\alpha\tilde S^2\cos 4\rho+\beta\tilde\sigma^2\cos 2\rho),
\label{vb}
\end{equation}
where we define $\sigma=\frac{\tilde\sigma}{\sqrt 2} e^{i\theta}$ and
$S=\frac{\tilde S}{\sqrt 2}e^{i\rho}$.
Along the minimum of $V_b$ for $\rho$ which is fixed by $\frac{\partial V_b}{\partial\rho}=0$,
the potential of $\tilde\sigma$ and $\tilde S$ can be written as
\begin{equation}
V(\tilde\sigma, \tilde S)=\frac{\tilde\kappa_\sigma}{4}(\tilde \sigma^2-w^2)^2+
\frac{\tilde\kappa_S}{4}(\tilde S^2-u^2)^2+
\frac{\kappa_{\sigma S}}{4}(\tilde\sigma^2-w^2)(\tilde S^2-u^2),
\label{potv}
\end{equation}
where $\tilde\kappa_\sigma$ and $\tilde \kappa_S$ are defined as
\begin{equation}
\tilde\kappa_\sigma=\kappa_\sigma-\frac{\beta^2}{4\alpha}, \qquad
\tilde\kappa_S=\kappa_S-2\alpha,
\end{equation}
and $w$ and $u$ are the vacuum expectation values (VEVs) of
$\tilde\sigma$ and $\tilde S$.
They are supposed to be much larger than the weak scale.
They keep the gauge symmetry but break down the global symmetry $U(1)\times Z_4$
into its diagonal subgroup $Z_2$.\footnote{It guarantees the stability of the
lightest $Z_2$ odd field, which could be a dark matter (DM) candidate as discussed later.}
Since the minimum of $V_b$ can be determined by using these VEVs as
$\cos 2\rho=-\frac{\beta}{4\alpha}\frac{w^2}{u^2}$
as long as $\left|\frac{\beta}{4\alpha}\frac{w^2}{u^2}\right|\le1$ is satisfied,
the $CP$ symmetry is spontaneously broken to result in a low energy effective
model with the $CP$ violation.
On the other hand, $\theta=0$ is satisfied because of
the global $U(1)$ symmetry relevant to $\sigma$ \cite{ssbcp}.
The stability condition for the potential (\ref{potv}) can be given as
\begin{equation}
\tilde\kappa_\sigma, \tilde\kappa_S>0, \qquad
4\tilde\kappa_\sigma\tilde\kappa_S>\kappa_{\sigma S}^2.
\label{stab1}
\end{equation}

If we consider the fluctuation of $\tilde\sigma$ and $\tilde S$ around the vacua
$\langle\tilde\sigma\rangle$ and $\langle\tilde S\rangle$,
mass eigenstates are the mixture of them in general.
If we take account of the stability condition (\ref{stab1}),
mass eigenvalues can be approximately expressed as
\begin{eqnarray}
&&m^2_{\tilde S}\simeq 2\left(\tilde\kappa_S-\frac{\kappa_{\sigma S}^2}
{4\tilde\kappa_\sigma}\right)u^2 \equiv 2\hat\kappa_Su^2, \quad
m^2_{\tilde\sigma}\simeq 2\tilde\kappa_\sigma w^2
\qquad {\rm for}~~\tilde \kappa_\sigma^2 w^2 \gg \tilde\kappa_S^2u^2,  \nonumber \\
&&m_{\tilde S}^2\simeq 2\tilde\kappa_S u^2, \quad m^2_{\tilde\sigma}
\simeq 2 \left(\tilde\kappa_\sigma-\frac{\kappa_{\sigma S}^2}
{4\tilde\kappa_S}\right)w^2\equiv 2\hat\kappa_\sigma w^2
\qquad {\rm for}~~\tilde \kappa_S^2 u^2 \gg \tilde\kappa_\sigma^2 w^2.
\end{eqnarray}
Although they have a tiny subcomponent in these cases,
 a dominant component of their eigenstates
is $\tilde S$ and $\tilde\sigma$, respectively.
The mass of an orthogonal component to $\tilde S$ is found to
be $m^2_{S_\perp}=8\alpha u^2\left(1-\cos^22\rho\right)$.
Since the global $U(1)$ symmetry works as the PQ symmetry mentioned above,
and the axion decay constant is given as $f_a=w$,
the VEV $w$ should satisfy the following condition \cite{pdg,fa}:
\begin{equation}
4\times 10^8~{\rm GeV}~{^<_\sim}~ w ~{^<_\sim}~10^{11} {\rm GeV}.
\end{equation}
The NG-boson caused by the spontaneous breaking of this $U(1)$ becomes axion
\cite{axion} which is characterized by a coupling with photons \cite{lmn}
\begin{equation}
g_{a\gamma\gamma}=\frac{1.51}{10^{10}~{\rm GeV}}\left(\frac{m_a}{\rm eV}\right).
\end{equation}

In the next part, we show that the effective model after the symmetry
breaking can have
 $CP$ phases in the CKM  and PMNS matrices.
They are induced by the mass matrices for
 the down type quarks and the charged leptons
through a similar mechanism which has been discussed in \cite{bbp}
as a simple realization of the NB mechanism \cite{nb} for the
strong $CP$ problem.

\subsection{$CP$ violating phases in CKM and PMNS matrices}
Yukawa couplings of down-type quarks and charged leptons given in eq.~(\ref{model})
derive mass terms as
\begin{equation}
(\bar f_{Li}, \bar F_L){\cal M}_f
\left(\begin{array}{c} f_{R_j} \\ F_R \\ \end{array} \right)+{\rm h.c.}, \qquad
{\cal M}_f=\left(\begin{array}{cc}
m_{f_{ij}} & 0 \\ {\cal F}_{f_j} & \mu_F \\
\end{array}\right),
\label{dmass}
\end{equation}
where $f$ and $F$ represent $f=d,e$ and $F=D,E$ for down-type quarks and
charged leptons and ${\cal M}_f$ is a $4\times 4$ matrix.
 Each component of ${\cal M}_f$ is expressed as
 $m_{f_{ij}}=h_{f_{ij}}\langle\tilde\phi\rangle$,
${\cal F}_{f_j}=(y_{f_j} ue^{i\rho} + \tilde y_{f_j}
ue^{-i\rho})$ and $\mu_F=y_F w$.
This mass matrix is found to have the same form proposed in \cite{bbp}.
Since the global $U(1)$ symmetry works as the PQ symmetry and all parameters in the model
are assumed to be real, ${\rm arg}({\rm det}{\cal M}_f)=0$ is satisfied even if
radiative effects are taken into account after the spontaneous breaking of
the $CP$ symmetry \cite{patisalam}.

We consider the diagonalization of a matrix ${\cal M}_f{\cal M}_f^\dagger$
by a unitary matrix
\begin{equation}
\left(\begin{array}{cc} A_f & B_f \\ C_f& D_f \\\end{array}\right)
\left(\begin{array}{cc} m_fm_f^{\dagger} & m_f{\cal F}_f^\dagger \\
  {\cal F}_fm_f^\dagger & \mu_F^2 +{\cal F}_f{\cal F}_f^\dagger \\
\end{array}\right)
\left(\begin{array}{cc} A_f^\dagger & C_f^\dagger \\ B_f^\dagger
& D_f^\dagger \\\end{array}\right)=
\left(\begin{array}{cc} \tilde m^2_f & 0 \\ 0 & \tilde M^2_F \\\end{array}\right),
\label{mass}
\end{equation}
where $\tilde m^2_f$ is  a $3\times 3$ diagonal matrix in which generation
indices are not explicitly written. Eq.~(\ref{mass}) requires
\begin{eqnarray}
   && m_fm_f^\dagger=A_f^\dagger \tilde m^2_f A_f+
C_f^\dagger \tilde M^2_FC_f, \qquad
  {\cal F}_fm_f^\dagger=B_f^\dagger m^2A_f+ D_f^\dagger \tilde M^2_FC_f, \nonumber \\
   && \mu_F^2+{\cal F}_f{\cal F}_f^\dagger=
B_f^\dagger \tilde m^2_f B_f+ D_f^\dagger \tilde M^2_F D_f.
\end{eqnarray}
If $\mu_F^2+{\cal F}_f{\cal F}_f^\dagger$ is much larger than each
component of ${\cal F}_fm_f^\dagger$, which can be realized in the case
$u, w\gg \langle \phi\rangle$,
we find that $B_f, C_f$, and $D_f$ are approximately given as
\begin{equation}
  B_f\simeq -\frac{A_fm_f{\cal F}_f^\dagger}{\mu_F^2
+{\cal F}_f{\cal F}_f^\dagger},
 \qquad C_f\simeq\frac{{\cal F}_f m_f^\dagger}{\mu_F^2
+{\cal F}_f{\cal F}_f^\dagger},
   \qquad D_f\simeq 1.
\end{equation}
These guarantee the unitarity of the matrix $A_f$ within the
present experimental bound \cite{pdg} since
$|B_f|, |C_f|\sim \frac{|m_f|}{|{\cal F}_f|} < O(10^{-7})$ is satisfied in each component.
In such a case, it is easy to find
\begin{equation}
A_f^{-1}\tilde m_f^2A_f\simeq m_fm_f^\dagger -\frac{1}{\mu_F^2
+{\cal F}_f{\cal F}_f^\dagger}
(m_f{\cal F}_f^\dagger)({\cal F}_fm_f^\dagger), \quad
\tilde M_F^2\simeq \mu_F^2+ {\cal F}_f{\cal F}_f^\dagger.
\label{ckm}
\end{equation}
The right-hand side of the first equation corresponds to an effective squared mass matrix
of the ordinary fermions $f$. It is derived through the mixing with the extra heavy fermions $F$.
Since its second term can have complex phases in off-diagonal
components as long as $y_{f_i}\not=\tilde y_{f_i}$ is satisfied,
the matrix $A_f$ could be complex. Moreover, if $\mu_F^2~{^<_\sim}~{\cal F}_f{\cal F}_f^\dagger$
is realized, the complex phase of $A_f$ in eq.~(\ref{ckm}) could have a substantial magnitude
because the second term in the right-hand side has a comparable magnitude with the first one.
Although it can be realized for various parameter settings, we consider
a rather simple situation here,\footnote{Leptogenesis
could depend on the strength of these couplings heavily in this model as studied later. }
\begin{equation}
\langle\phi\rangle \ll w < u, \qquad  y_{f_j}\sim \tilde y_{f_j} < y_F.
\label{uw}
\end{equation}
Since the masses of vector-like fermions are expected to be of $O(10^8)$~GeV
or larger in this case, they decouple from the SM and do not contribute to
phenomena at TeV regions.

The CKM matrix is determined as $V_{\rm CKM}={O_u}^TA_d$, where $O_u$ is an
orthogonal matrix used for the diagonalization of a mass matrix for up-type quarks.
Thus, the $CP$ phase of $V_{\rm CKM}$ is caused by the one of $A_d$.
The same argument is applied to the leptonic sector and
the PMNS matrix is derived as $V_{\rm PMNS}=A_e^\dagger U_\nu$, where
$U_\nu$ is an orthogonal matrix used for the diagonalization
of a neutrino mass matrix.
The Dirac $CP$ phase in the CKM matrix and the PMNS matrix can be
induced from the same origin of $CP$ violation.
A concrete example of $A_d$ is given for a simple case in Appendix A.

\subsection{Effective model at a lower energy region}
An effective model at lower energy regions than $w$ and
$u$ can be obtained by integrating
out the heavy fields. It is reduced to the SM with a lepton sector extended as
the scotogenic neutrino mass model \cite{ma},
which is characterized by the terms invariant
under the remaining $Z_2$ symmetry
\begin{eqnarray}
-{\cal L}_{\rm scot}&=&\sum_{\alpha=1}^3
\left[\sum_{j=1}^3\left(\tilde h_{\alpha j}^\ast\bar\ell_\alpha\eta N_j
+\frac{M_{N_j}}{2}\bar N_j^cN_j \right)
+ \sum_{\beta=1}^3\frac{y_{\alpha\beta}}{M_\ast}(\bar\ell_\alpha\phi)(\bar\ell_\beta\phi)+
{\rm h.c.}\right] \nonumber \\
&+&\tilde m_\phi^2\phi^\dagger\phi+\tilde m_\eta^2\eta^\dagger\eta+
\tilde\lambda_1(\phi^\dagger\phi)^2+ \tilde\lambda_2(\eta^\dagger\eta)^2
+\tilde\lambda_3(\phi^\dagger\phi)(\eta^\dagger\eta)
+\lambda_4(\phi^\dagger\eta)(\eta^\dagger\phi) \nonumber \\
&+&\frac{\tilde\lambda_5}{2}\left[(\phi^\dagger\eta)^2+{\rm h.c.}\right],
\label{model0}
\end{eqnarray}
where neutrino Yukawa couplings $\tilde h_{\alpha j}$ are defined on the basis
for which the mass matrix of the charged leptons is diagonalized as discussed
in the previous part. Thus, they are complex now.
After the spontaneous breaking due to the VEVs of $\tilde\sigma$ and $\tilde S$,
coupling constants in eq.~(\ref{model0}) are related to the ones contained in eq.~(\ref{model}) as
\begin{eqnarray}
&&\tilde\lambda_1=\lambda_1-\frac{\kappa_{\phi\sigma}^2}{4\tilde\kappa_\sigma}
-\frac{\kappa_{\phi S}^2}{4\tilde\kappa_S}+
\frac{\kappa_{\sigma S}\kappa_{\phi\sigma}
\kappa_{\phi S}}{4\tilde\kappa_\sigma\tilde\kappa_S}, \qquad
\tilde\lambda_2=\lambda_2-\frac{\kappa_{\eta\sigma}^2}{4\tilde\kappa_\sigma}
-\frac{\kappa_{\eta S}^2}{4\tilde\kappa_S}+
\frac{\kappa_{\sigma S}\kappa_{\eta\sigma}
\kappa_{\eta S}}{4\tilde\kappa_\sigma\tilde\kappa_S}, \nonumber\\
&&\tilde\lambda_3=\lambda_3-\frac{\kappa_{\phi\sigma}\kappa_{\eta\sigma}}
{2\tilde\kappa_\sigma}
-\frac{\kappa_{\phi S}\kappa_{\eta S}}{2\tilde\kappa_S}+
\frac{\kappa_{\sigma S}\kappa_{\phi\sigma}\kappa_{\eta S}
+\kappa_{\sigma S}\kappa_{\eta\sigma}\kappa_{\phi S}}{4\tilde\kappa_\sigma\tilde\kappa_S},
\qquad \tilde\lambda_5=\lambda_5\frac{w}{M_\ast}.
\end{eqnarray}
These connection conditions should be imposed
at a certain scale $\bar M$, which is taken to be $\bar M=\tilde M_F$ in the present study.
Stability of the potential (\ref{model0}) requires the following conditions
to be satisfied
through scales $\mu<\bar M$:
\begin{equation}
\tilde\lambda_1,~\tilde\lambda_2 >0, \qquad \tilde\lambda_3,~
\tilde\lambda_3+\lambda_4-|\tilde\lambda_5|>
-2\sqrt{\tilde\lambda_1\tilde\lambda_2}.
\label{stab2}
\end{equation}
Potential stability (\ref{stab1}) and (\ref{stab2}) and perturbativity
of the model from the weak scale to the Planck scale can be examined by using
renormalization group equations (RGEs) for the coupling constants.
Relevant RGEs at $\mu>\bar M$ are given in Appendix B.
The mass parameters in eq.~(\ref{model0}) are represented as
\begin{eqnarray}
&&M_{N_j}=y_{N_j}w,\nonumber\\
&&\tilde m_\phi^2=m_\phi^2+\left(\kappa_{\phi\sigma}+
\frac{\kappa_{\phi S}\kappa_{\sigma S}}{2\tilde\kappa_S}\right)w^2
+\left(\kappa_{\phi S}+
\frac{\kappa_{\phi\sigma}\kappa_{\sigma S}}{2\tilde\kappa_\sigma}\right)u^2, \nonumber\\
&&\tilde m_\eta^2=m_\eta^2+\left(\kappa_{\eta\sigma}+
\frac{\kappa_{\eta S}\kappa_{\sigma S}}{2\tilde\kappa_S}\right)w^2
+\left(\kappa_{\eta S}+
\frac{\kappa_{\eta\sigma}\kappa_{\sigma S}}{2\tilde\kappa_\sigma}\right)u^2.
\label{cpara}
\end{eqnarray}
If $\tilde m_\eta^2>0$ is satisfied and $\eta$ has no VEV, $Z_2$ is kept as an exact
symmetry of the model. In this model, we assume both $|\tilde m_\phi|$ and $\tilde m_\eta$
have values of $O(1)$~TeV.
Since it has to be realized under the contributions from the VEVs $w$ and $u$,
parameter tunings are required.\footnote{For parameter values
assumed in this analysis, the order of required tuning is estimated as $O(10^{-10})$.}

\begin{figure}[t]
\begin{center}
\includegraphics[width=7cm]{dmrelic.pdf}
\hspace*{5mm}
\includegraphics[width=7cm]{rge-l.pdf}
\end{center}

\footnotesize{{\bf Fig.~1}~~Left:~contours of $\Omega h^2=0.12$ are plotted in the
$(\lambda_+,\tilde\lambda_3)$ plane
by a solid line for $M_{\eta_R}=0.9$ (green), 1 (red) and 1.1 (blue) in a TeV unit.
Since $\eta_R$ should be lighter than the charged components, $\lambda_4<0$ should be satisfied
which corresponds to a region above the diagonal black solid line.
Direct search bound from Xenon1T \cite{xenon} is shown by the same color
dashed line for each $M_{\eta_R}$.
Stability conditions (\ref{stab2}) restrict the allowed region to the area marked in the upper
 right quadrant (marked off by the dot-dashed lines), for which $\tilde\lambda_2=0.1$ is assumed.
Right:~an example of the running of coupling constants $\tilde\lambda_i$ for
$M=10^\mu~{\rm GeV}$. Initial values of $\tilde\lambda_3$ and $\lambda_4$ are fixed as $\tilde\lambda_3=0.445$ and
$\lambda_4=-0.545$, which are included in the allowed region of the left panel
for $M_{\eta_R}=1$~TeV.  $\tilde\lambda_3+2\sqrt{\tilde\lambda_1\tilde\lambda_2}$ and
$\lambda_++2\sqrt{\tilde\lambda_1\tilde\lambda_2}$
are also plotted by red and blue lines, respectively.}
\end{figure}

Phenomenology on neutrinos and DM could be the same
as the one which has been studied extensively in various studies \cite{scot,ks}
unless the axion is a dominant component of DM.
If the lightest neutral component of $\eta$ is DM which is identified as
its real component $\eta_R$, both DM
relic abundance and DM direct search constrain the parameters $\tilde\lambda_3$ and
$\lambda_4$ \cite{mksvz,hs}. As a reference, in Fig.~1 we show their required values
in the $(\lambda_+,\tilde\lambda_3)$ plane for the cases
$M_{\eta_R}=$0.9, 1 and 1.1 TeV where $\lambda_+=\tilde\lambda_3+\lambda_4+\tilde\lambda_5$
and $M_{\eta_R}^2=m_\phi^2+\lambda_+\langle\phi\rangle^2$.
They should be also consistent with the stability condition (\ref{stab2}).
The figure shows that these could be satisfied with rather restricted values of $\tilde\lambda_3$
and $\lambda_4$. A perturbativity requirement at $\mu>\bar M$ also constrains the model
strongly since DM relic abundance requires $\tilde\lambda_3$ and $|\lambda_4|$
to take rather large values. We have to take account of it to consider the model
at high energy regions. As an example, we show the result of RGE study in the
right panel of Fig.~1. Details of the assumed initial values are given in Appendix B.
If the initial values of $\tilde\lambda_3$ and $\lambda_4$ are chosen
at a larger $\tilde\lambda_3$ region in the left panel,
the perturbatibity condition is violated at a much lower region than the Planck scale.
We do not plot the behavior of $\kappa_i$ in this figure since their values are fixed to  be very
small and they do not change their values substantially.
If they satisfy the condition (\ref{stab1}), it can be confirmed to be kept up to the Planck scale.

Neutrino mass is forbidden at tree level due to this $Z_2$ symmetry except for
the ones generated through dimension five Weinberg operators in eq.~(\ref{model0}).
They could give a substantial contribution to the neutrino masses depending on
the cut-off scale $M_\ast$ and coupling constants $y_{\alpha\beta}$.
In the present study, however, we assume that their contribution is
negligible and the dominant contribution comes from one-loop diagrams with
$\eta$ and $N_j$ in internal lines. Its formula is given as
\begin{equation}
{\cal M}^\nu_{\alpha\beta}\simeq\sum_{j=1}^3\tilde h_{\alpha j}\tilde
h_{\beta j}\tilde\lambda_5\Lambda_j,  \qquad
\Lambda_j=\frac{\langle\phi\rangle^2}
{8\pi^2}\frac{1}{M_{N_j}}\ln\frac{M_{N_j}^2}{M_\eta^2},
\label{nmass}
\end{equation}
where $M_\eta^2=\tilde m_\eta^2+(\tilde\lambda_3+
\lambda_4)\langle\phi\rangle^2$
and $M_{N_j}\gg M_\eta$ is supposed.
As an example, one may assume a simple flavour
structure for neutrino Yukawa couplings \cite{tribi}
\begin{equation}
\tilde h_{e i}=0, ~\tilde h_{\mu i}=\tilde h_{\tau i}\equiv h_i~~(i=1,~2); \quad
\tilde h_{e 3}=\tilde h_{\mu 3}=-\tilde h_{\tau 3}\equiv h_3.
\label{tribi}
\end{equation}
This realizes a tri-bimaximal mixing which gives a simple and rather good 0-th order
approximation for the analysis of neutrino oscillation data and leptogenesis \cite{ks}.
If we impose the mass eigenvalues obtained from eq.~(\ref{nmass}) for the case
$|h_1|\ll |h_2|, |h_3|$ to satisfy the squared mass differences required by the neutrino
oscillation data, we find
\begin{equation}
|h_2^2\tilde\lambda_5|\Lambda_2=\frac{1}{2}\sqrt{\Delta m_{32}^2}, \qquad
|h_3^2\tilde\lambda_5|\Lambda_3=\frac{1}{3}\sqrt{\Delta m_{21}^2}.
\label{oscilcon}
\end{equation}
Since we have $\Lambda_{2,3}\sim7\times 10^5$~eV for $M_{2,3}\sim 10^{7}$~GeV
and $M_\eta\sim 10^3$~GeV,
the neutrino oscillation data \cite{pdg} can be explained by taking as an example
\begin{equation}
y_{N_j}\sim 10^{-2}, \quad  |h_2|\sim 6\times 10^{-3},  \quad  |h_3|\sim 2\times 10^{-3},
\quad |\tilde\lambda_5|\sim 10^{-3}.
\label{npara}
\end{equation}
Even if we impose the neutrino oscillation data, $h_1$ can take a very small value
compared with $h_{2,3}$ \cite{ks}.
It can play a crucial role for low scale leptogenesis as seen later.

\section{Inflation due to singlet scalars}
\subsection{Inflation}
It is well known that a scalar field which couples nonminimally with the Ricci scalar
can cause inflation of the universe
and the idea has been applied to the Higgs scalar in the SM \cite{higgsinf} and
its singlet scalar extensions \cite{sinfl,ext-s}.
If the singlet scalars $S$ and $\sigma$, which are related to the $CP$
issues in the SM couple with the Ricci scalar,
it can play the role of inflaton in this model.
The action relevant to the inflation is given in the Jordan frame as
\begin{eqnarray}
S_J &=& \int d^4x\sqrt{-g} \left[-\frac{1}{2}M_{\rm pl}^2R -
\xi_\sigma\sigma^\dagger\sigma R-
\xi_{S1} S^\dagger S R -\frac{\xi_{S2}}{2}(S^2+S^{\dagger 2})R\right. \nonumber \\
&+&\left.\partial^\mu \sigma^\dagger \partial_\mu \sigma
+\partial^\mu S^\dagger \partial_\mu S - V(\sigma, S) \right],
\label{inflag}
\end{eqnarray}
where $M_{\rm pl}$ is the reduced Planck mass and the coupling of $S$
is controlled by the $Z_4$ symmetry.
$V(\sigma, S)$ stands for the corresponding part in the potential (\ref{model}).
Since inflation follows very complicated dynamics if multi scalars contribute to it,
we confine our study to the inflation in a potential valley.
Moreover,  we assume $\xi_{S1}=-\xi_{S2}$ is satisfied for simplicity.
In that case, the coupling of $S$ with the Ricci scalar is reduced to $\frac{1}{2}\xi_S S_I^2R$,
where $S=\frac{1}{\sqrt 2}(S_R+iS_I)$ and $\xi_S=\xi_{S1}-\xi_{S2}$.
If $S$ is supposed to evolve along a constant $\rho$, which is determined as
a potential minimum $\frac{\partial V_b}{\partial\rho}=0$, the radial component $\tilde S$
couples with the Ricci scalar as $\frac{1}{2}\tilde\xi_S\tilde S^2R$,
where $\tilde\xi_S$ is defined as $\tilde\xi_S=\xi_S\sin^2\rho$ and
the potential $V(\sigma, S)$ is expressed by eq.(\ref{potv}).
Here we consider cases such that both $\xi_\sigma$ and $\tilde\xi_S$ are
positive only.
Stability of this potential requires the condition given in eq.~(\ref{stab1}).
We neglect the VEVs $w$ and $u$ for a while since they are
much smaller than $O(M_{\rm pl})$, which is the field values of $\tilde\sigma$ and $\tilde S$
during the inflation. We also suppose that other scalars have much smaller values
than them.

We consider the conformal transformation for a metric tensor in the Jordan frame
\begin{equation}
\tilde g_{\mu\nu}=\Omega^2g_{\mu\nu}, \qquad
\Omega^2=1+\frac{(\xi_\sigma\tilde\sigma^2+\tilde\xi_S\tilde S^2)}{M_{\rm pl}^2}.
\end{equation}
After this transformation to the Einstein frame where the Ricci scalar term takes a canonical form,
the action can be written as
\begin{eqnarray}
S_E &=& \int d^4x\sqrt{-\tilde g} \Big[-\frac{1}{2}M_{\rm pl}^2\tilde R
+\frac{1}{2}\partial^\mu\chi_\sigma \partial_\mu \chi_\sigma+\frac{1}{2}
\partial^\mu\chi_S \partial_\mu \chi_S \nonumber\\
&+&\frac{6\xi_\sigma\tilde\xi_S\frac{\tilde\sigma\tilde S}{M_{\rm pl}^2}}
{\left[(\Omega^2+\frac{6\xi_\sigma^2}{M_{\rm pl}^2}\tilde\sigma^2)
(\Omega^2+\frac{6\tilde\xi_S^2}{M_{\rm pl}^2}\tilde S^2)\right]^{1/2}}
\partial^\mu\chi_\sigma\partial_\mu\chi_S
- \frac{1}{\Omega^4}V(\tilde\sigma, \tilde S) \Big],
\label{inflag}
\end{eqnarray}
where $\chi_\sigma$ and $\chi_S$ are defined as \cite{sinfl}
\begin{equation}
\frac{\partial\chi_\sigma}{\partial\tilde\sigma}
=\frac{1}{\Omega^2}\sqrt{\Omega^2+6\xi_\sigma^2
\frac{\tilde\sigma^2}{M_{\rm pl}^2}}, \qquad
\frac{\partial\chi_S}{\partial\tilde S}
=\frac{1}{\Omega^2}\sqrt{\Omega^2+6\tilde\xi_S^2
\frac{\tilde S^2}{M_{\rm pl}^2}}.
\label{canoninf}
\end{equation}
If we introduce variables $\tilde\chi,~\varphi$ to express $\tilde\sigma$ and $\tilde S$ as
$\tilde\sigma=\tilde\chi\cos\varphi,~ \tilde S=\tilde\chi\sin\varphi$,
the potential in the Einstein frame
at the large field regions
such as $\xi_\sigma\tilde\sigma^2+\tilde\xi_S\tilde S^2\gg M_{\rm pl}^2$
can be written as
 \begin{equation}
V(\tilde\chi,\varphi)=\frac{M_{\rm pl}^4}{4}\frac{\tilde\kappa_S\sin^4\varphi+
\tilde\kappa_\sigma\cos^4\varphi+\kappa_{\sigma S}
\sin^2\varphi\cos^2\varphi}{(\xi_\sigma\cos^2\varphi+\tilde\xi_S\sin^2\varphi)^2}.
\label{vk}
\end{equation}
We find that there are three types of valley along the minimum in the $\varphi$
direction of this potential. They realize different types of inflaton. Two of them are
\begin{equation}
{\rm (i)}~\varphi=0  ~~ {\rm for}~~ 2\tilde\kappa_\sigma\tilde\xi_S<
\kappa_{\sigma S}\xi_\sigma, \qquad
{\rm (ii)}~\varphi=\frac{\pi}{2} ~~ {\rm for} ~~ 2\tilde\kappa_S\xi_\sigma<
\kappa_{\sigma S}\tilde\xi_S.
\end{equation}
In each case, a kinetic term mixing between $\chi_\sigma$ and $\chi_S$
disappears and inflaton can be identified with
$\chi_\sigma$ for (i) and $\chi_S$ for (ii),
respectively.\footnote{In different context, the inflaton dominated by $\tilde\sigma$
and the $\tilde S$ inflaton have been discussed in
\cite{smash} and \cite{patisalam}, respectively.}

Another valley which is studied in this paper is realized at
\begin{equation}
\sin^2\varphi=\frac{2\tilde\kappa_\sigma\tilde\xi_S-\kappa_{\sigma S}\xi_\sigma}
{(2\tilde\kappa_S\xi_\sigma-\kappa_{\sigma S}\tilde\xi_S)
+(2\tilde\kappa_\sigma\tilde\xi_S-\kappa_{\sigma S}\xi_\sigma)},
\label{varphi}
\end{equation}
under the conditions
\begin{equation}
 2\tilde\kappa_\sigma\tilde\xi_S>\kappa_{\sigma S}\xi_\sigma, \quad
2\tilde\kappa_S\xi_\sigma>\kappa_{\sigma S}\tilde\xi_S,
\end{equation}
where we note that these are automatically satisfied for $\kappa_{\sigma S}<0$
since eq.~(5) is imposed, and $\tilde\xi_S$ and $\xi_\sigma$ are assumed to be positive.
In this case the inflaton $\tilde\chi$ is a mixture of $\tilde\sigma$ and $\tilde S$
with a constant value of $\tilde\sigma/\tilde S$.
Although the kinetic term mixing cannot be neglected for a general $\sin\varphi$,
it can be safely neglected if we restrict it to the one in which the inflaton is dominated
by $\tilde S~(\sin^2\varphi\simeq 1)$ or $\tilde\sigma~(\sin^2\varphi\simeq 0)$.
We focus our study on the former case where
$\tilde\chi\gg\tilde\sigma$ is always satisfied during inflation.
If we additionally impose $\tilde\xi_S \gg \xi_\sigma$ on eq.~(\ref{varphi})
and assume that the relevant couplings satisfy\footnote{If we assume
$\tilde\kappa_S>|\kappa_{\sigma S}|$, $\sin\varphi$ is differently expressed as
$\sin^2\varphi=1-\frac{\tilde\kappa_S\xi_\sigma}{\tilde\kappa_\sigma\tilde\xi_S}$.
However, we do not consider such a case in this paper.}
\begin{equation}
\kappa_{\sigma S}<0, \qquad \tilde\kappa_S< |\kappa_{\sigma S}|< \tilde\kappa_\sigma,
\label{ks}
\end{equation}
$\sin\varphi$ is expressed as
$\sin^2\varphi=1+\frac{\kappa_{\sigma S}}{2\tilde\kappa_\sigma}$.
In this case, by using $\hat\kappa_S=\tilde\kappa_S-\frac{\kappa_{\sigma S}^2}
{4\tilde\kappa_\sigma}$ potential can be expressed as $V(\tilde S)=\frac{1}{4\Omega^4}
\hat\kappa_S\tilde S^4$ at the bottom of the valley and
$\cos 2\rho=-\frac{\beta}{4\alpha}\cot^2\varphi$ realizes the minimum of $V_b$ if
$\tan^2\varphi>\left|\frac{\beta}{4\alpha}\right|$ is satisfied.
Nature of the inflaton $\tilde\chi$ is fixed by the parameters $\tilde\kappa_S$,
$\tilde\kappa_\sigma$ and $\kappa_{\sigma S}$.

Squared mass of the orthogonal component to $\tilde\chi$ during inflation can be
estimated as $m_{\tilde\chi_\perp}^2=\frac{|\kappa_{\sigma S}|M_{\rm pl}^2}{2\tilde\xi_S^2}$.
Since the Hubble parameter $H_I$ satisfies
$H_I^2=\frac{\hat\kappa_SM_{\rm pl}^2}{12\tilde\xi_S^2}$ at the same period,
$H_I < m_{\tilde\chi_\perp}$ is satisfied under the condition (\ref{ks}).\footnote{
Mass of another orthogonal component $S_I$
is estimated as $m_{S_I}^2\sim \frac{8\alpha M_{\rm pl}^2}{\tilde\xi_S^2}$ for a case
$\cos 2\rho\sim 0$, which is interested in the present study
and then $96\alpha >\hat\kappa_S$ is required for $m_{S_I} >H_I$.
Although $\alpha \ll \hat\kappa_S$ is assumed from the vacuum determination,
it can be satisfied for suitable values of $\alpha$.}
 Thus, the inflation $\chi$ starts rolling along the valley within a few Hubble time
independently of an initial value of the inflaton. It justifies our analyzing the model as
a single field inflation model.
On the other hand, since $\tilde\sigma(=\sqrt\frac{|\kappa_{\sigma S}|}
{2\tilde\kappa_\sigma}\tilde\chi)>\frac{H_I}{\sqrt{2\pi}}$ is satisfied
generally,\footnote{Although it may be escapable for
$\tilde\xi_S\ll\sqrt{|\kappa_{\sigma S}|}$, it is not the case in the present model.}
the global $U(1)$ is spontaneously broken during inflation and isocurvature
fluctuation could be problematic \cite{smash}.
However, even in that case it is escapable since the axion
needs not to be a dominant component of DM in the present model.
This problem is discussed later.

The canonically normalized inflaton $\chi$ can be expressed as \cite{smash}
\begin{equation}
\Omega^2\frac{d\chi}{d\tilde S}=\sqrt{\gamma\Omega^2+6\tilde\xi_S^2
\frac{\tilde S^2}{M_{\rm pl}^2}},
\label{tchi}
\end{equation}
where $\gamma$ can be approximated along the valley as
\begin{equation}
\gamma=1-\frac{\kappa_{\sigma S}}{2\tilde\kappa_\sigma}.
\end{equation}
If we use $\gamma\simeq 1$, the potential of $\chi$ obtained through
$V(\chi)=\frac{1}{\Omega^4}V(\tilde\sigma,\tilde S)$
can be derived by using the solution of eq.~(\ref{tchi}) which is given as
\begin{equation}
\frac{\chi}{M_{\rm pl}}=-\sqrt 6~{\rm arcsinh}\left(\frac{\sqrt\frac{6}{\gamma}
\frac{\tilde\xi_S \tilde S}{M_{\rm pl}}}
{\sqrt{1+\frac{\tilde\xi_S}{M_{\rm pl}^2}\tilde S^2}}\right)+
\sqrt\frac{\gamma+6\tilde\xi_S}{\tilde\xi_S}{\rm arcsinh}
\left(\frac{\sqrt{\tilde\xi_S(1+\frac{6}{\gamma}\tilde\xi_S)}\tilde S}{M_{\rm pl}}\right).
\label{etchi}
\end{equation}
We derive the potential of $\chi$ through numerical calculation for a
typical value of $\tilde\xi_S$ by using eq.~(\ref{etchi}).
Such examples of $V(\chi)$ are shown in Fig.~2.
It can be approximated as
\begin{eqnarray}
V(\chi)=\left\{ \begin{array}{ll}
\frac{\hat\kappa_S}{4\tilde\xi_S^2}M_{\rm pl}^4 &\quad \chi>M_{\rm pl}, \\
\frac{\hat\kappa_S}{6\tilde\xi_S^2}M_{\rm pl}^2\chi^2
&\quad \frac{M_{\rm pl}}{\tilde\xi_S}
<\chi< M_{\rm pl}, \\
 \frac{\hat\kappa_S}{4}\chi^4 &\quad \chi < \frac{M_{\rm pl}}{\tilde\xi_S}.
\end{array}\right.
\label{oscilpot}
\end{eqnarray}
The inflation ends at $\chi_{\rm end}\simeq M_{\rm pl}$.
After the end of inflation, there is a substantial region where the potential
behaves as a quadratic form before it is reduced to a quartic form
at low energy regions for the case $\tilde\xi_S\gg 1$ as in the Higgs inflation.
However, such a region can be neglected for the case $\tilde\xi_S<10$,
which is the case considered in this study.
Since the inflaton oscillating in the quartic potential behaves
as radiation as shown later, the radiation domination starts
soon after the end of inflation in that case.

\begin{figure}[t]
\begin{center}
\includegraphics[width=7.5cm]{pot1.pdf}
\hspace*{3mm}
\includegraphics[width=7.5cm]{pot10.pdf}
\end{center}

\footnotesize{{\bf Fig.~2}~~Potential of the inflaton $\chi$ for
$\tilde\xi_S=1$ (left panel) and $\tilde\xi_S=10$ (right panel).
In both panels, $\tilde\xi_S/\xi_\sigma=20$ and
$\tilde\kappa_S/|\kappa_{\sigma S}|=|\kappa_{\sigma S}|/\tilde\kappa_\sigma=0.1$
are assumed and $\tilde\kappa_S$ is fixed by using
eq.~(\ref{kap}) for ${\cal N}_k=55$. As references, we also plot
approximated potential $\frac{\hat\kappa_S}{6\tilde\xi_S^2}M_{\rm pl}^2\chi^2$
and $\frac{\hat\kappa_S}{4}\chi^4$ in eq.~(\ref{oscilpot}) as $\chi^2$ and $\chi^4$.
In these plots, a Planck unit $(M_{\rm pl}=1)$ is used.}
\end{figure}

The slow-roll parameters in this model can be estimated by using
eq.~(\ref{tchi}) as \cite{inf}
\begin{equation}
\epsilon\equiv\frac{M_{\rm pl}^2}{2}\left(\frac{V^\prime}{V}\right)^2
=\frac{8M_{\rm pl}^4}{\gamma\tilde\xi_S\left(1+\frac{6}{\gamma}\tilde\xi_S\right)\tilde S^4}, \qquad
\eta\equiv M_{\rm pl}^2\frac{V^{\prime\prime}}{V}=
-\frac{8M_{\rm pl}^2}{\gamma\left(1+\frac{6}{\gamma}
\tilde\xi_S\right)\tilde S^2}.
\end{equation}
The $e$-foldings number ${\cal N}_k$ from the time when the scale
$k$ exits the horizon to the end of inflation is estimated by using eq.~(\ref{tchi}) as
\begin{eqnarray}
{\cal N}_k=\frac{1}{M_{\rm pl}^2}\int^{\chi_k}_{\chi_{\rm end}}\frac{V}{V^\prime}d\chi
=\frac{1}{8M_{\rm pl}^2}(\gamma+6\tilde\xi_S)(\tilde S_k^2-\tilde S_{\rm end}^2)
-\frac{3}{4}\ln\frac{M_{\rm pl}^2+\tilde\xi_S\tilde S_k^2}
{M_{\rm pl}^2+\tilde\xi_S\tilde S_{\rm end}^2}.
\label{efold}
\end{eqnarray}
Taking account of these,
the slow-roll parameters in this inflation scenario is found to be approximated as
$\epsilon\simeq \frac{3}{4{\cal N}_k^2}$ and $\eta\simeq -\frac{1}{{\cal N}_k}$.
The field value of inflaton during the inflation is found to be  expressed as
$\chi_k=\frac{\sqrt 6}{2}M_{\rm pl}\ln(32\tilde\xi_S{\cal N}_k)$
by using eqs.~(\ref{etchi}) and (\ref{efold}), and its potential $V_k(\equiv V(\chi_k))$ takes
a constant value as shown in eq.~(\ref{oscilpot}).
On the other hand, if we use $\epsilon=1$ at the end of inflation,
 the inflaton potential is estimated as
$V_{\rm end}(\equiv V(\chi_{\rm end}))\simeq 0.072
\frac{\hat\kappa_S}{\tilde\xi_S^2}M_{\rm pl}^4$ which is found to be a
good approximation from Fig.~2.

The spectrum of density perturbation predicted by the inflation
is known to be expressed as \cite{inf}
\begin{equation}
{\cal P}(k)=A_s\left(\frac{k}{k_\ast}\right)^{n_s-1},  \qquad
A_s=\frac{V}{24\pi^2M_{\rm pl}^4\epsilon}\Big|_{k_\ast}.
\label{power}
\end{equation}
If we use the Planck data $A_s=(2.101^{+0.031}_{-0.034})\times 10^{-9}$
at $k_\ast=0.05~{\rm Mpc}^{-1}$ \cite{planck18}, we find the
Hubble parameter during the inflation to be $H_I=1.4\times 10^{13}
\left(\frac{60}{{\cal N}_{k_\ast}}\right)$~GeV and the relation
\begin{equation}
\hat\kappa_S\simeq 4.13\times 10^{-10}\tilde\xi_S^2 \left(\frac{60}{{\cal N}_{k_\ast}}\right)^2,
\label{kap}
\end{equation}
which should be satisfied at the horizon exit time of the scale $k_\ast$.
We confine our study to the case $\tilde\xi_S<10$.

\begin{figure}[t]
\begin{center}
\includegraphics[width=8cm]{ns-r.pdf}
\end{center}

\footnotesize{{\bf Fig.~3}~~Predicted values of the scalar spectral index $n_s$
and the tensor-to-scalar ratio $r$ in this model. They are read off as the values
at intersection points of two lines with a fixed value of $\tilde\xi_S$ or ${\cal N}_k$.
The coupling constant $\hat\kappa_S$ is varied in a range
from $10^{-7}$ to $10^{-10}$.}
\end{figure}

In Fig.~3, we plot predicted values
for the scalar spectral index $n_s$ and the tensor-to-scalar ratio $r$
in the present model.
Since the quartic coupling $\hat\kappa_S$ is a free parameter of the model
under the constraint (\ref{kap}),
we vary $\hat\kappa_S$ in the range $10^{-10}\le \hat\kappa_S\le 10^{-7}$
for fixed values of $\tilde\xi_S$ or ${\cal N}_k$.
The CMB constraint (\ref{kap}) is satisfied at intersection points of
the lines with a fixed value of $\tilde\xi_S$ or ${\cal N}_k$.
The figure shows that the constraints of the observed CMB data
\cite{planck18} are satisfied for the supposed parameters.

After the end of inflation, the inflaton $\chi$ starts oscillation in the potential $V(\chi)$.
At this stage, the description by $\chi$ is no longer justified, especially, at
the small field regions. $\varphi$ is not constant in general there. $\tilde S$ and $\tilde\sigma$
should be treated independently. In the following study, however,
we confine our study to a special inflaton trajectory and estimate the reheating
phenomena by using $\chi$ to give a rough evaluation of reheating temperature
under the assumption that inflaton follows a constant $\varphi$ trajectory.\footnote{
During the first several oscillations, both $\varphi$ and $\rho$ can
be numerically confirmed to take constant values.
This suggests that single field treatment is rather good
during the first few oscillations at least.}
In this case, inflaton oscillation is described by the equation
\begin{equation}
\frac{d^2\chi}{dt^2}+3H\frac{d\chi}{dt}+V^\prime(\chi)=0.
\label{infeq}
\end{equation}
Since the amplitude of $\chi$ evolves approximately as
$\Phi(t)=\frac{\tilde\xi_S}{\sqrt{\pi\hat\kappa_S}t}$ in the quadratic
potential after the end of inflation, the inflaton $\chi$ oscillates
$\frac{1}{2\pi\sqrt{3\pi}}(\tilde\xi_S-1)$ times before the potential (\ref{oscilpot})
changes from a quadratic form to a quartic one.
This means that preheating under the quadratic potential could play
no substantial role for the case $\tilde\xi_S<10$.
In such a case, we need to consider the preheating in the quartic potential only.
The model with the quartic potential $V(\chi)=\frac{\hat\kappa_S}{4}\chi^4$
becomes conformally invariant \cite{confpre}.

If we introduce dimensionless conformal time $\tau$, which is defined
by using a scale factor $a$ as $a d\tau=\sqrt{\hat\kappa_S}\chi_{\rm end}dt$
and also a rescaled field $f=\frac{a\chi}{ \chi_{\rm end}}$,
eq.~(\ref{infeq}) can be rewritten as
\begin{equation}
\frac{d^2f}{d\tau^2}+f^3=0.
\label{infeq1}
\end{equation}
The solution of this equation which describes the inflaton oscillation is known
to be given by a Jacobi elliptic function
$f(\tau)={\rm cn} \left(\tau-\tau_i,\frac{1}{\sqrt 2}\right)$.\footnote{If we take
$\tau_i\simeq 2.44$, $f(\tau)$ can be approximated by $\cos\left(\frac{2\pi}{\tau_0}\tau\right)$,
where $\tau_0$ is expressed by using the complete elliptic integral of the first kind $K$
as $\tau_0=4K(\frac{1}{\sqrt 2})$ \cite{confpre}.}
From the Friedman equation for this inflaton oscillation, we find
\begin{equation}
a(\tau)=\frac{\chi_{\rm end}}{2\sqrt3 M_{\rm pl}}\tau, \qquad
\tau=2(3\hat\kappa_S M_{\rm pl}^2)^{1/4}\sqrt t.
\label{rd}
\end{equation}
Since $H=1/2t$ is satisfied, this oscillation era is radiation dominated.
If we take into account such a feature of the model that radiation domination
starts just after the end of inflation, the $e$-foldings number ${\cal N}_k$
can be expressed by noting a relation $k=a_kH_k$ as
\begin{equation}
{\cal N}_k=56.7-\ln\left(\frac{k}{a_0H_0}\right)+
\ln\left(\frac{V_{\rm end}^{1/4}}{10^{14}~{\rm GeV}}\right)
+2\ln\left(\frac{V_k^{1/4}}{V_{\rm end}^{1/4}}\right),
\label{efolds}
\end{equation}
where $H_k^2=\frac{V_k}{3M_{\rm pl}^2}$ and the suffix 0 stands for the present
value of each quantity. Reheating temperature dependence of ${\cal N}_k$ is weak or lost differently
from the usual case \cite{inf} where substantial matter domination is assumed
to follow the inflation era.\footnote{If reheating occurs through a perturbative process
at $\chi~{^<_\sim}~u$ where matter domination is realized, its effect on ${\cal N}_k$
could also be negligible as long as $\Gamma>H$ is satisfied at that stage where
$\Gamma$ is inflaton decay width.}
In the next part, we discuss the reheating temperature expected to be
realized in the present model.

\subsection{Preheating and reheating}
Before proceeding to the study of particle production under the background
oscillation of inflaton, we need to know the mass of the relevant particles
which is induced through the interaction with inflaton.
Such interactions are given as
\begin{eqnarray}
&&\Big[-\frac{y_D}{\sqrt 2}\frac{\kappa_{\sigma S}}{2\tilde\kappa_\sigma}
\chi\bar D_LD_R
-\frac{y_E}{\sqrt 2}\frac{\kappa_{\sigma S}}{2\tilde\kappa_\sigma}
\chi\bar E_LE_R+
\sum_{j=1}^3\Big\{
\frac{1}{\sqrt 2}\left(y_{d_j} e^{i\rho}
+ \tilde y_{d_j} e^{-i\rho}\right)\chi\bar D_L d_{R_j} \nonumber\\
&&+\frac{1}{\sqrt 2}\left(y_{e_j} e^{i\rho} +
\tilde y_{e_j} e^{-i\rho}\right)\chi\bar E_L e_{R_j}
-\frac{y_{N_j}}{2\sqrt 2}\frac{\kappa_{\sigma S}}{2\tilde\kappa_\sigma}
\chi \bar N_j^cN_j \Big\}
+ {\rm h.c.}\Big] \nonumber \\
&&+\frac{1}{2}(\kappa_{\phi S}\phi^\dagger\phi+\kappa_{\eta S}\eta^\dagger\eta)
\chi^2
-\frac{\kappa_{\sigma S}}{2\tilde\kappa_\sigma}
(\kappa_{\phi\sigma}\phi^\dagger\phi
+\kappa_{\phi\sigma}\eta^\dagger\eta)\chi^2.
\label{oscm}
\end{eqnarray}
The particles interacting with the inflaton $\chi$ have mass varying
with the oscillation of $\chi$ and their mass can be read off from
eq.~(\ref{oscm}) as
\begin{eqnarray}
&&M_{N_j}\simeq \frac{y_{N_j}}{\sqrt 2}\frac{|\kappa_{\sigma S}|}
{2\tilde\kappa_\sigma}\chi, \quad
\tilde M_{F}\simeq \frac{\chi}{\sqrt 2}\left[\sum_{j=1}^3(y_{f_j}^2+\tilde y_{f_j}^2) +
y_F^2\frac{\kappa_{\sigma S}^2}
{4\tilde\kappa_\sigma^2}\right]^{\frac{1}{2}}, \nonumber \\
&&m_\phi^2\simeq\frac{1}{2}\left(\kappa_{\phi S}+
\frac{|\kappa_{\sigma S}|}{\tilde\kappa_\sigma}\kappa_{\phi\sigma}\right)\chi^2,
\quad
m_\eta^2\simeq\frac{1}{2}\left(\kappa_{\eta S}+
\frac{|\kappa_{\sigma S}|}{\tilde\kappa_\sigma}
\kappa_{\eta\sigma}\right)\chi^2,
\label{imass1}
\end{eqnarray}
where $F=D$ or $E$ should be understood for $f=d$ or $e$, respectively.
Since the effect of nonminimal coupling is negligible during this oscillation period,
it is convenient to use the components of $\sigma$ and $S$, which are
parallel and orthogonal to the inflaton $\chi$ to describe their interactions.
If we indicate each of them as $\sigma_\parallel$, $\sigma_\perp$,
$S_\parallel$, and $S_\perp$, then their interactions are expressed as
\begin{eqnarray}
&&\frac{\tilde\kappa_S}{4}(S_\parallel^2+S_\perp^2-u^2)^2
+\frac{\kappa_{\sigma S}}{4}(S_\parallel^2+S_\perp^2-u^2)
(\sigma_\parallel^2+\sigma_\perp^2-w^2)
+\frac{\tilde\kappa_\sigma}{4}\left(\sigma_\parallel^2
+\sigma_\perp^2-w^2\right)^2.\nonumber \\
\end{eqnarray}
By combining these interactions with the composition of $\chi$,
their masses are found to be given by\footnote{It should be noted
that the mass of $\sigma$
could have another non-negligible contribution which is
induced by explicit breaking of the global $U(1)$ symmetry
brought about by the quantum gravitational effects.
We do not take account of it in the present study.}
\begin{eqnarray}
&& m^2_{S_\parallel}\simeq \left(3\hat\kappa_S+
\frac{\kappa_{\sigma S}^2}{2\tilde\kappa_\sigma}\right)\chi^2, \qquad
m_{\sigma_\parallel}^2\simeq \left(|\kappa_{\sigma S}|+
\frac{\kappa_{\sigma S}^2}{4\tilde\kappa_\sigma}\right)\chi^2, \nonumber \\
&&m^2_{S_\perp}\simeq\hat\kappa_S \chi^2, \qquad
m_{\sigma_\perp}^2\simeq
\frac{\kappa_{\sigma S}^2}{4\tilde\kappa_\sigma}\chi^2.
\label{imass2}
\end{eqnarray}

The coupling constants relevant to these masses are restricted
through the assumed inflaton composition and the realization of the
$CP$ phases in the CKM and PMNS matrices.
The discussion in the previous sections shows that such requirements
are satisfied for
\begin{equation}
\hat\kappa_S < |\kappa_{\sigma S}|<\tilde\kappa_\sigma, \qquad
\hat\kappa_S< y_{N_j}, y_{f_j}, \tilde y_{f_j}.
\label{rehpara}
\end{equation}
We assume additionally
\begin{equation}
\hat\kappa_S < g_\phi\equiv \kappa_{\phi S}+
\frac{|\kappa_{\sigma S}|}{\tilde\kappa_\sigma}\kappa_{\phi\sigma},\qquad
\hat\kappa_S < g_\eta\equiv \kappa_{\eta S}+
\frac{|\kappa_{\sigma S}|}{\tilde\kappa_\sigma}\kappa_{\eta\sigma}.
\label{pecoupl}
\end{equation}
Since the oscillation frequency of the inflaton is
$\sim \sqrt{\hat\kappa_S}\chi$, decays or annihilations of the inflaton are
kinematically forbidden except for the one to $\sigma_\perp$
as found from eqs.~(\ref{imass1}) and (\ref{imass2}).
In $\sigma_\perp$ case, the inflaton reaction rate to it is much smaller than the
Hubble parameter at this period because of the smallness of its coupling with the inflaton,
energy drain from the inflaton to $\sigma_\perp$ is ineffective to be neglected.
As a result, the energy transfer from the inflaton oscillation to excited particles
is expected to occur at the time when the inflaton crosses the zero where
the resonant particle production is possible.

Preheating under the background inflaton oscillation can generate the excitations of
$\chi$ itself and other scalars $\psi$ which couple with $\chi$
at its zero crossing \cite{pre}.
In a quartic potential case \cite{confpre}, the model becomes
conformally invariant and
the time evolution equations of $\chi_k(\simeq S_k)$
and $\psi_k$, which are the comoving modes with a momentum $k$,
 can be transformed to the simple ones by rescaling them
to the dimensionless quantities
in the  same way as eq.~(\ref{infeq1}). They are given as
\begin{eqnarray}
&& \frac{d^2}{d\tau^2}X_k+\omega_k^2X_k=0, \qquad
\omega_k^2=\bar k^2+3f(\tau)^2, \nonumber \\
&& \frac{d}{d\tau^2}F_k+\tilde\omega_k^2F_k=0, \qquad
\tilde\omega_k^2=\bar k^2+\frac{g_\psi}{\hat\kappa_S}f(\tau)^2,
\label{teq}
\end{eqnarray}
where the rescaled variables are defined as
\begin{eqnarray}
X_k=\frac{a\chi_k}{\chi_{\rm end} }, \qquad
F_k=\frac{a\psi_k}{\chi_{\rm end} }, \qquad
\bar k=\frac{ak}{\chi_{\rm end}\sqrt{\hat\kappa_S}}.
\end{eqnarray}
Function $f(\tau)$ is the solution of eq.~(\ref{infeq1}) and $g_\psi$ stands
for a coupling constant of the relevant particle
$\psi(=\sigma, S_{\perp},\phi,\eta)$ with
the inflaton $\chi$, which can read off from eqs.~(\ref{imass1}) and (\ref{imass2}).
Amplitudes $X_k$ and $F_k$ are known to show the exponential behavior
$\propto e^{\mu_k\tau}$ with a characteristic exponent $\mu_k$, which is
determined by a parameter $g_\psi/\hat\kappa_S$.
Using the solutions of eq.~(\ref{teq}), the number density
of the produced particle $\psi$ can be calculated as
\begin{equation}
n_k^\psi=\frac{\tilde\omega_k}{2\hat\kappa_S}
\left(\frac{|F_k^\prime|^2}{\tilde\omega_k^2} +|F_k|^2\right)-\frac{1}{2}.
\end{equation}

Particle production based on eq.~(\ref{teq}) at the inflaton zero crossing
has been studied in \cite{confpre} and it is shown to be characterized by
the parameter $g_\psi/\hat\kappa_S$.
We classify the relevant couplings into five groups
\begin{equation}
{\rm (A)}~\frac{g_{\sigma_\parallel}}{\hat\kappa_S}\gg 1, \quad
{\rm (B)}~\frac{g_{S_\parallel}}{\hat\kappa_S}=3, \quad
{\rm (C)}~\frac{g_{S_\perp}}{\hat\kappa_S}=1, \quad
{\rm (D)}~\frac{g_{\sigma_\perp}}{\hat\kappa_S} \ll 1,\quad
{\rm (E)}~\frac{g_\phi}{\hat\kappa_S},~\frac{g_\eta}{\hat\kappa_S}>1,
\end{equation}
where we note that couplings in (A)-(D) are fixed by the present
inflaton composition but the ones in (E) are not constrained.
Now we consider the resonant particle production in each group.
A maximum value of characteristic exponent in (D) is very small so that it plays no effective
role also in preheating.
In (B) and (C), both the fluctuations of $S_\parallel$
and $S_\perp$ are produced fast, but it stops as soon as $\langle |S_\parallel|^2\rangle$
and $\langle |S_\perp|^2\rangle$ reach a certain value such as $0.5\chi_{\rm end}^2/a^2$.
Although a maximum value $\mu_{\max}$ of the characteristic exponent of (B) is much
smaller than the one of (C) and also the resonance band of (B) is much narrower than (C),
the interaction $S_\parallel^2S_\perp^2$ accelerates the production of fluctuations
of $S_\parallel$ through rescattering and they reach the similar value \cite{smash, rescat}.
Since the backreaction of these fluctuations to the inflaton oscillation restructures
the resonance band,  the resonant particle production stops before causing
much more conversion of the inflaton oscillation energy to particle excitations.
Moreover, since the decay of excitations produced through these processes are also closed
kinematically, these could not play an efficient role in reheating.
In (A), since $\sigma_\parallel$ also couples to the inflaton directly,
the resonant production of its excitation stops at a certain stage due to the same
reason as (B) and (C). Even if the excited particles are allowed to
decay to fermions $F$
and $N_j$ kinematically, the decay width is much smaller than
the Hubble parameters to be neglected.
As a result, if the process due to (E) is not effective,
preheating cannot play any role for reheating and reheating proceeds
through perturbative processes after the amplitude of inflaton is smaller
than the VEV $u$.

Here, we have to note that there is a possibility in (E) where the energy
transfer from the inflaton oscillation to radiation proceeds through preheating
since the produced excitations can decay to relativistic particles differently from (B) and (C).
In this case, $\phi$ and $\eta$ are produced as excitations at the zero crossing
of the inflaton where an adiabaticity condition $\tilde\omega_k^\prime<\tilde\omega_k^2$
could be violated for certain values of $\bar k$.
By using the analytic solution of eq.~(\ref{teq}) derived in \cite{confpre},
the momentum distribution $n_{k}^\psi$ of the produced particle $\psi$ through
one zero crossing of the inflaton can be estimated as
\begin{equation}
n_{\bar k}^\psi= e^{2\mu_k\frac{\tau_0}{2}}=e^{-\left(\bar k/\bar k_c\right)^2},
\qquad \bar k_c^2=\sqrt\frac{g_\psi}{2\pi^2\hat\kappa_S},
\end{equation}
where $\tau_0$ is an inflaton oscillation period and $\tau_0=7.416$.
The resonance is efficient for $\bar k<\bar k_c$.
Thus, the particle number density produced during one zero crossing of the
inflaton is
\begin{equation}
\bar n^\psi =\int\frac{d^3\bar k }{(2\pi)^3}n_{\bar k}^\psi
=\int\frac{d^3\bar k}{(2\pi)^3}e^{-(\bar k/\bar k_c)^2}
=\frac{\bar k_c^3}{8\pi^{3/2}}.
\end{equation}

The energy transfer from the inflaton oscillation to relativistic particles
is caused through the decay of the produced particles $\psi(=\phi, \eta)$
and thermalization proceeds.
They can decay to light fermions through $\phi\rightarrow \bar qt$
with a top Yukawa coupling $h_t$
and $\eta \rightarrow \bar\ell N$ with neutrino Yukawa couplings $h_j$, respectively.
Here, we should note that $\eta$ can be heavier than $N_j$ at this
stage even if $\eta$
is the lightest one with $Z_2$ odd parity at the weak scale.
It is caused by the inflaton composition in the present model
as found from eq.~(\ref{imass1}).
Their decay widths in the comoving frame are given
by using the conformally rescaled unit as
\begin{eqnarray}
&&\bar\Gamma_\psi=\frac{c_\psi y_\psi^2}{8\pi}\bar m_\psi,  \qquad
\bar m_\psi=\frac{a m_\psi}{ \chi_{\rm end} \sqrt{\hat\kappa_S}}
=\sqrt\frac{g_\psi}{\hat\kappa_S}f(\tau),
\end{eqnarray}
where $\psi=\phi,~\eta$, and $c_\psi$ are internal degrees
of freedom $c_\phi=3$ and $c_\eta=1$.
The Yukawa coupling $y_\psi$ represents $y_\phi=h_t$ and $y_\eta=h_j$.
Since $\bar\Gamma_{\psi}^{-1}<\tau_0/2$ is satisfied for $g_\psi> 4\times
10^{-7}\left(\frac{\hat\kappa_S}{10^{-8}}\right)$, the produced $\psi$
decays to the light fermions completely before the next inflaton
zero crossing \cite{inst} and then it is not accumulated in such cases.
We fix $\tau=0$ at the first inflaton zero crossing so that $f(\tau)$ can be
expressed approximately as $f(\tau)=f_0\sin(cf_0\tau)$.
Transferred energy density through the $\psi$ decay during a half period of
oscillation can be estimated as\footnote{$\bar\rho$ is
defined as the energy density in the comoving frame
by using the conformally rescaled variables.}
\begin{equation}
\delta\bar\rho_r=\int^{\tau_0/2}_0d\tau\bar\Gamma_\psi\bar m_\psi\bar n_\psi
e^{-\int_0^\tau\bar\Gamma_\psi\tau^\prime}
=\frac{1}{8\pi^{3/2}(2\pi^2)^{3/4}}\left(\frac{g_\psi}{\hat\kappa_S}\right)^{5/4}
Y(f_0,\gamma_\psi),
\end{equation}
where $\gamma_\psi$ and $Y(f_0,\gamma_\psi)$ are defined by using
$c=2\pi/\tau_0$ as
\begin{equation}
\gamma_\psi=\frac{c_\psi y_\psi^2}{8\pi c}\sqrt\frac{g_\psi}{\hat\kappa_S}, \qquad
Y(f_0,\gamma_\psi)=c\gamma_\psi\int^{\tau_0/2}_0d\tau f_0^2\sin^2(cf_0\tau)
e^{-2\gamma_\psi\sin^2(\frac{cf_0\tau}{2})}.
\end{equation}

The energy density transferred to the light particles is accumulated at each
inflaton zero crossing linearly and its averaged value for $\tau$ is
estimated as
\begin{equation}
\bar\rho_r(\tau)=\frac{2\tau}{\tau_0}\delta\bar\rho_r
=6.5\times 10^{-4}\left(\frac{g_\psi}{\hat\kappa_S}\right)^{5/4}
Y(f_0,\gamma_\psi)\tau,
\label{rde}
\end{equation}
where the substantial change of $f_0$ is assumed to be negligible
during $\tau$.
Since the total energy density of the inflaton oscillation energy $\bar\rho_\chi$
and the transferred energy $\bar\rho_r$ to light particles
is conserved, reheating temperature realized through
this process can be estimated from $\bar\rho_{\chi_{\rm end}}=\bar\rho_r$.
It can be written by transferring it to the physical unit as
\begin{equation}
\frac{1}{4\hat\kappa_S}\left(\frac{\sqrt{\hat\kappa_S}
\chi_{\rm end}}{a}\right)^4=\frac{\pi^2}{30}g_\ast T_R^4,
\end{equation}
where we use $\bar\rho_{\chi_{\rm end}}=\frac{1}{4\hat\kappa_S}$ and $g_\ast=130$.
By applying eqs.~(\ref{rd}) and (\ref{rde}) to this formula,
we find\footnote{The same result can be obtained by using the relation
$H=\frac{1}{2t}$ in the radiation dominated era together
with eqs.~(\ref{rd}) and (\ref{rde}).}
\begin{equation}
T_R=5.9\times 10^{15}g_\psi^{5/4}Y(f_0,\gamma_\psi)~{\rm GeV}.
\label{nonptr}
\end{equation}
Since $h_t\gg h_j$ is satisfied, reheating temperature is expected to
be determined by
the produced $\phi$ as long as $\phi$ is dominantly produced.

If preheating cannot produce relativistic particles effectively,
the dominant energy is still kept in the inflaton oscillation.
When the oscillation amplitude of $\chi$ decreases to be $O(u)$,
the inflaton starts decaying to the light particles through the perturbative processes.
Since the mass pattern is expected under the present assumption
for the coupling constants in (\ref{rehpara}) to be
\begin{equation}
2\tilde m_\eta<m_\chi< \tilde M_D, \tilde M_E,
\label{mspec}
\end{equation}
the inflaton decay is expected to occur mainly through
$\chi\rightarrow \eta^\dagger\eta$ and
$\chi\rightarrow \phi^\dagger\phi$ at tree level.
The decay width of $\psi(=\phi,~\eta)$ is estimated as
\begin{eqnarray}
\Gamma_\psi&\simeq&
\frac{g_\psi^2}{16\pi \hat\kappa_S}m_\chi,
\label{width}
\end{eqnarray}
where $g_\psi$ is defined in eq.~(\ref{pecoupl}).
After the inflaton decays to $\eta^\dagger\eta$ and $\phi^\dagger\phi$,
the SM contents are expected to be thermalized through gauge interactions
with $\eta$ and $\phi$ immediately.
Since $\Gamma_\psi > H$ is satisfied for $g_\psi>10^{-7.1}
\left(\frac{\hat\kappa_S}{10^{-8}}\right)^{1/2}
\left(\frac{u}{10^{11}~{\rm GeV}}\right)^{1/2}$ at $\chi\simeq u$,
reheating temperature in such a case can be estimated through
$\frac{1}{4}\hat\kappa_Su^4=\frac{\pi^2}{30}g_\ast T_R^4$ as\footnote{
Although a larger value of $u$ can make the reheating
temperature much higher, its upper bound exits.
Since a larger $g_\psi$ is required in that case,
 it could violate the perturbativity of the model and cause the upper bound for it.
For example, if we consider the case with $y_F=10^{-1.2}$,
the perturbativity is violated for $g_\psi>10^{-4.4}$. As a result,
the reheating temperature due to the perturbative process
is bounded as $T_R<6.3 \times 10^{13}$ GeV.}
\begin{equation}
T_R\simeq
2.8\times10^8\left(\frac{\hat\kappa_S}{10^{-8}}\right)^{1/4}
\left(\frac{u}{10^{11}~{\rm GeV}}\right)~{\rm GeV},
\end{equation}
which is independent of $g_\psi$. However, if $\Gamma_\psi > H$ is not satisfied
because of a small $g_\psi$, the reheating temperature is expected to be determined
through $\Gamma_\psi=H$ and then becomes smaller proportionally to $g_\psi$.

\begin{figure}[t]
\begin{center}
\includegraphics[width=7.5cm]{bothtr.pdf}
\hspace*{5mm}
\includegraphics[width=7.5cm]{ptr.pdf}
\end{center}

\footnotesize{{\bf Fig.~4}~~Left : reheating temperature $T_R$ (GeV) predicted
for both the preheating (a) and the perturbative process (b).
They are plotted as a function of $g_\phi$.
In the case (a), the decay $\phi\rightarrow \bar qt$ is assumed.
Right : contribution to the reheating temperature due to
a component $\sigma$ of the inflaton in a case $\kappa_{\phi S}=0$.
The reheating temperature $T_R$ (GeV) is plotted as a function of $\kappa_{\phi\sigma}$.
Each line represents $T_R$ for several values of
$(\tilde\kappa_S, |\kappa_{\sigma S}|/\tilde\kappa_\sigma)$, where
$\tilde\kappa_\sigma$ is fixed at $10^{-4.5}$ and $10^{-5.3}$ for
$\tilde\kappa_S=10^{-8}$ and $10^{-9}$, respectively.}
\end{figure}

In the left panel of Fig.~4, for a case $\psi=\phi$,
the expected reheating temperature through both processes is plotted
as a function $g_\phi$ in a case
$\tilde\kappa_S=10^{-8}$, $\tilde\kappa_\sigma=10^{-4.5}$,
$|\kappa_{\sigma S}|/\tilde\kappa_\sigma=10^{-1.2}$, and $u=10^{11}$ GeV.
It shows that the reheating temperature is determined by the perturbative
process at $g_\phi<10^{-6}$.
We also found from the figure that the reheating at $g_\phi>10^{-6}$
proceeds through the preheating. Even if the dominant component $S$ of
the inflaton has no coupling with $\phi$ so that
$g_\phi\simeq|\kappa_{\sigma S}|\kappa_{\phi\sigma}/\tilde\kappa_\sigma$,
the preheating is caused by the component $\sigma$.
It is shown in the right panel where the reheating temperature is plotted
for $\kappa_{\phi\sigma}$ by varying
$|\kappa_{\sigma S}|/\tilde\kappa_\sigma$.\footnote{The condition
$\bar\Gamma_\psi^{-1}<\tau_0/2$ can be
confirmed for the parameters used here. }
These figures show that $T_R>2.3\times 10^8$ GeV can be realized
if $g_\phi>4\times 10^{-8}$ is satisfied.
However, since the perturbativity of the model is found to be violated
at $g_\phi>10^{-4.4}$ as mentioned in the previous footnote,
$\kappa_{\phi\sigma}<10^{-4.4}$ and $\kappa_{\phi\sigma}<10^{-2.6}$
should be satisfied, and then the reheating temperature cannot be
higher than $\sim 10^{10}$~GeV as found from the figure.
Since the decay of $\phi$ is so effective, it decays soon
after their production and much before the inflaton amplitude becomes large
during the oscillation. This makes the energy transfer in the preheating
inefficient.

In the usual leptogenesis in the seesaw scenario, the right-handed
neutrinos are supposed to be thermalized only through the neutrino
Yukawa couplings $h_j$.
In the present model, neutrino mass eigenvalues obtained from eq.~(\ref{nmass})
require $h_{2,3}=O(10^{-3})$ to explain the neutrino oscillation data
as discussed at a part of eq.~(\ref{npara}).
On the other hand, reheating temperature is found to satisfy
$T_R~{^>_\sim}~10^8$ GeV from the above discussion.
Since the decay width $\Gamma_{N_{2,3}}$ of $N_{2,3}$ and the reheating
temperature $T_R$ satisfy $\Gamma_{N_{2,3}}>H(T_R)$
and $T_R>M_{N_{2,3}}$, $N_{2,3}$ are also expected to be
in the thermal equilibrium through the inverse decay simultaneously
at the reheating period.
In the case of $N_1$, however, it depends on the magnitude of its Yukawa
coupling $h_1$ which can be much smaller than others.
We should note that $N_1$ could be effectively generated in the thermal bath,
even if $h_1$ is extremely small, through the scattering of extra fermions
which are expected in the thermal equilibrium through gauge interactions
in the case $\tilde M_D, \tilde M_E< T_R$.
It is a noticeable feature of the present model which opens a window for
low scale leptogenesis.

\section{Phenomenological signature of the model}
\subsection{Leptogenesis}
The most interesting feature of this inflation scenario is that
thermal leptogenesis could generate sufficient baryon number asymmetry even
for $M_{N_1}<10^9$ GeV without relying on resonance effect.
In the ordinary seesaw framework,
neutrino mass is generated as
$(m_\nu)_{\alpha\beta}=\frac{h_{\alpha j}h_{\beta j}\langle\phi\rangle^2}{M_{N_j}}$
through Yukawa interaction $h_{\alpha j}\bar\ell_\alpha\phi N_j$.
Baryon number asymmetry in the universe \cite{baryon} is
expected to be generated by the same interaction through
thermal leptogenesis \cite{leptg}.
If we assume the sufficient lepton asymmetry is generated
through the out-of-equilibrium decay of the lightest right-handed neutrino,
which has been in the thermal equilibrium, then
the reheating temperature $T_R$ is required to be larger than its mass
$T_R>M_{N_1}$.
Moreover, since it has to be produced sufficiently in the thermal bath,
its Yukawa coupling $h_{\alpha 1}$ should not be so small. On the other hand,
the neutrino mass formula gives a severer upper bound on $h_{\alpha 1}$
for a smaller $M_{N_1}$ under the constraints of neutrino oscillation data.
These impose a lower bound for $M_{N_1}$ such as $10^9$ GeV \cite{di}.
This condition for $M_{N_1}$ is not changed even if $T_R\gg 10^9$ GeV is satisfied.
The problem is caused by such a feature of the model that both the production and
the out-of-equilibrium decay of the right-handed neutrino have to be caused only
by the same neutrino Yukawa coupling.
It does not change in the original scotogenic model either \cite{ks}.
In that model, the right-handed neutrino mass can be much smaller than $10^9$ GeV
keeping the neutrino Yukawa couplings to be rather larger values by fixing
$|\lambda_5|$ at a smaller value in a consistent way with
the neutrino oscillation data. However,
the washout of the generated lepton number due to the inverse decay of the right-handed
neutrinos becomes so effective in that case.
As a result, successful leptogenesis cannot be realized
for a lighter right-handed neutrino than $10^8$ GeV.\footnote{Low scale
leptogenesis in the scotogenic model has been studied intensively
in \cite{lowlept1}. However, the lightest
right-handed neutrino is assumed to be in the thermal equilibrium initially there.}
It is a notable aspect in the present model that this situation
can be changed by the particles
which are introduced to explain the $CP$ issues in the SM.

We note that the interaction between the right-handed
neutrino $N_1$ and extra vector-like fermions $F$ mediated by $\tilde\sigma$
could change the situation.\footnote{The similar mechanism has been discussed
in models with a different type of inflaton \cite{hs,lowlept2}.}
The lightest right-handed neutrino $N_1$ can be effectively produced
in the thermal bath through the extra fermions scattering
$\bar D_LD_R, \bar E_LE_R\rightarrow N_1N_1$
mediated by $\tilde\sigma$ if $D_{L,R}$ and/or $E_{L,R}$ are in the thermal
equilibrium at a certain temperature $T$.  In that case,
both conditions $T>\tilde M_F, M_{N_1}$ and $\Gamma_{FF}\simeq H(T)$
are required to be satisfied,
 where $\Gamma_{FF}$ is the reaction rate of this scattering.
Mass of these fermions is determined by the VEVs $u$ and $w$ which should
be larger than the lower bound of PQ symmetry breaking scale.
Since the rough estimation of  $\Gamma_{FF}\simeq H(T)$ for
relativistic $F$ and $N_1$ gives
\begin{equation}
T\simeq 5.8\times 10^8\left(\frac{y_F}{10^{-1.2}}\right)^2
\left(\frac{y_{N_1}}{10^{-2}}\right)^2 {\rm GeV},
\label{scat}
\end{equation}
we find that $T>\tilde M_F, M_{N_1}$ could be satisfied for suitable values of
$y_F$ and $y_{N_1}$.
It is crucial that this does not depend on the magnitude of
the $N_1$ Yukawa coupling $h_1$.
If an extremely small value is assumed for $h_1$,
successful leptogenesis is allowed in a consistent way with neutrino oscillation data
even for $M_{N_1}<10^9$ GeV.

After $N_1$ is produced in the thermal bath through the scattering of
the extra fermions mediated by $\tilde\sigma$,
it is expected to decay to $\ell_\alpha\eta^\dagger$ by a strongly suppressed
Yukawa coupling.
Since its substantial decay occurs after the washout processes
are frozen out, the generated lepton number asymmetry can be efficiently
converted to the baryon number asymmetry through sphaleron processes.
This scenario can be checked by solving Boltzmann equations for $Y_{N_1}$ and
$Y_L(\equiv Y_\ell-Y_{\bar\ell})$, where $Y_\psi$ is
defined as $Y_\psi=\frac{n_\psi}{s}$
by using the $\psi$ number density $n_\psi$ and the entropy density $s$.
Boltzmann equations analyzed here are given as
\begin{eqnarray}
&&\frac{dY_{N_1}}{dz}=-\frac{z}{sH(M_{N_1})}\left(\frac{Y_{N_1}}{Y^{\rm eq}_{N_1}}
-1\right)\left[ \gamma_D^{N_1}+
\left(\frac{Y_{N_1}}{Y^{\rm eq}_{N_1}}+1\right)\sum_{F=D,E}
\gamma_{F}\right], \nonumber \\
&&\frac{dY_{L}}{dz}=-\frac{z}{sH(M_{N_1})}\left[\varepsilon
\left(\frac{Y_{N_1}}{Y^{\rm eq}_{N_1}}-1\right)
\gamma_D^{N_1} -\frac{2Y_L}{Y_\ell^{\rm eq}}
\sum_{j=1,2,3}\left(\frac{\gamma_D^{N_j}}{4} +\gamma_{N_j}\right)\right],
\label{beq}
\end{eqnarray}
where $z=\frac{M_{N_1}}{T}$ and an equilibrium value of $Y_\psi$ is
represented by $Y_\psi^{\rm eq}$.
$H(T)$ is the Hubble parameter at temperature $T$ and
the $CP$ asymmetry $\varepsilon$ for the decay of $N_1$ is expressed as
\begin{eqnarray}
\varepsilon&=&\frac{1}{8\pi}\sum_{j=2,3}
\frac{{\rm Im}[ \sum_{\alpha}(\tilde h_{\alpha 1}
\tilde h_{\alpha j}^\ast)]^2}
{\sum_{\alpha} \tilde h_{\alpha 1}\tilde h_{\alpha 1}^\ast}
F\left(\frac{M_{N_j}^2}{M_{N_1}^2}\right)
\nonumber\\
&=&\frac{1}{16 \pi}\left[4|h_2|^2F\left(\frac{y_{N_2}^2}{y_{N_1}^2}\right)
\sin 2(\theta_1-\theta_2)+|h_3|^2F\left(\frac{y_{N_3}^2}{y_{N_1}^2}\right)
\sin 2(\theta_1-\theta_3)\right],
\label{epsilon}
\end{eqnarray}
where $h_j=|h_j|e^{i\theta_j}$ and $F(x)=\sqrt x\left[1-(1+x)\ln\frac{1+x}{x}\right]$.
A reaction density for the decay
$N_j\rightarrow \ell_\alpha\eta^\dagger$ and for the lepton number violating
scattering mediated by $N_j$ is expressed by $\gamma_D^{N_j}$
and $\gamma_{N_j}$, respectively \cite{ks}.
$\gamma_{F}$ represents a reaction density for the scattering
$\bar D_LD_R,\bar E_LE_R\rightarrow N_1N_1$.
We assume that $(D_L,D_R)$ and $(E_L,E_R)$ are in the thermal equilibrium
and $Y_{N_1}=Y_L=0$ at $z=z_R(\equiv\frac{M_{N_1}}{T_R})$.

Now we fix the model parameters for numerical study of eq.~(\ref{beq}) by
taking account of the discussion in the previous part.
We consider two cases for the VEVs of the singlet scalars such that
\begin{equation}
{\rm (I)}~ w=10^9~{\rm GeV}, ~u=10^{11}~{\rm GeV},
\qquad  {\rm (II)}~ w=10^{11}~{\rm GeV},~ u=10^{13}~{\rm GeV},
\end{equation}
where the axion could be a dominant DM in case (II).
The parameters $\tilde\kappa_S, \tilde\kappa_\sigma$, and $\kappa_{\sigma S}$,
which characterize the inflaton $\chi$, are fixed to $\tilde\kappa_S=10^{-8},
\tilde\kappa_\sigma=10^{-4.5}$, and $|\kappa_{\sigma S}|=10^{-6.1}$.
These are used in the right panel of  Fig.~4.
The condition ${\cal F}_f{\cal F}_f^\dagger> \mu_F^2$ for which the
$CP$ phases in the CKM and PMNS matrices can be generated is
reformulated as $\delta\equiv\tilde M_F/\mu_F>\sqrt 2$.
If we confine our study on a case $y_f=(0,0,y)$ and $\tilde y_f=(0,\tilde y,0)$
for simplicity,\footnote{It is considered as an example in Appendix A.}
we have a relation $y^2+\tilde y^2=(\delta^2-1)\frac{w^2}{u^2}y_F^2$
among Yukawa couplings of the extra fermions.
We fix them as $\delta=\sqrt 3$, $\tilde y/y=0.5$ and $y_D=y_E=10^{-1.2}$
at the scale $\bar M$.
Parameters relevant to the neutrino mass generation are fixed as
\begin{eqnarray}
&& y_{N_2}=2\times 10^{-2},  \quad  y_{N_3}=4\times 10^{-2}
\quad {\rm for ~(I)~ and~ (II)}, \nonumber\\
&& y_{N_1}=7\times 10^{-3}, \quad |h_1|=6\times 10^{-7}
\quad |\tilde\lambda_5|=10^{-3}, \quad M_\eta=1~{\rm TeV}
\quad {\rm for~ (I)},
\nonumber \\
&&y_{N_1}= 10^{-3}, \quad |h_1|=6\times 10^{-5},
\quad |\tilde\lambda_5|=5\times 10^{-3}, \quad M_\eta=0.9~{\rm TeV}
\quad {\rm for~ (II)}.
\label{para}
\end{eqnarray}
These parameters fix the mass of relevant particles as
\begin{eqnarray}
{\rm (I)}&& M_{N_1}=7\times 10^6~{\rm GeV}, \quad M_{N_2}= 2\times 10^7~{\rm GeV}, \quad
M_{N_3}=4\times 10^7~{\rm GeV}, \nonumber \\
&&\tilde M_D= \tilde M_E=1.1\times 10^8~{\rm GeV}, \quad m_{\tilde\sigma}=8\times 10^6~{\rm GeV},
\nonumber \\
{\rm (II)}&&M_{N_1}=10^8~{\rm GeV},  \quad M_{N_2}= 2\times 10^9~{\rm GeV}, \quad
M_{N_3}=4\times 10^9~{\rm GeV}, \nonumber \\
&&\tilde M_D= \tilde M_E=1.1\times 10^{10}~{\rm GeV}. \quad m_{\tilde\sigma}=8\times 10^8~{\rm GeV}.
\end{eqnarray}
Although the mediator has a small component $\tilde S$, it can be safely
treated as $\tilde\sigma$.
For these parameters, the $CP$ asymmetry $\varepsilon$ in the $N_1$ decay takes
a value of $O(10^{-6})$ in both cases if the maximum $CP$ phase is assumed.
DM is determined by the couplings $\tilde\lambda_3$ and $\lambda_4$.
Since they are fixed so as to realize the correct DM abundance by
the neutral component of $\eta$ in the case (I), it cannot saturate the required
DM abundance for the same $\tilde\lambda_3$ and $\lambda_4$ in the case (II)
as found from Fig.~1. The axion could be a dominant component of DM
in the case (II) since $w$ is taken to be a sufficient value for it.

We give a remark on these couplings here.
It is crucial to examine whether the above parameters used in this analysis are
consistent with the potential stability conditions (\ref{stab1}), (\ref{stab2}) and
also the perturbativity of the model under constraints coming
from the requirements for the DM relic abundance and the reheating temperature.
If DM relic abundance is realized by the neutral component of $\eta$,
both $\tilde\lambda_3$ and $\lambda_4$ should take values shown in Fig.~1.
On the other hand, the reheating temperature required for sufficient
leptogenesis can be realized for $\kappa_{\phi\sigma}~{^>_\sim}~10^{-4}$
or $\kappa_{\phi S}~{^>_\sim}~10^{-7}$ as found from the analysis of
the reheating temperature.
Since they can give rather large contributions to the $\beta$-functions
of the scalar quartic couplings $\kappa_{\eta\sigma}$ and $\kappa_{\eta S}$
for example, the perturbativity up to the inflation scale could be violated.
An upper bound on $g_\phi$ has to be imposed to escape it and
it results in an upper bound on the reheating temperature discussed already.
The parameter sets used here have been confirmed to satisfy these
conditions through the RGEs study. Details of the used parameters
in the analysis are addressed in Appendix B and an example of the results of this
study is presented in the right panel of Fig.~1.

\begin{figure}[t]
\begin{center}
\includegraphics[width=7.5cm]{fbasym-w9.pdf}
\hspace*{3mm}
\includegraphics[width=7.5cm]{fbasym-w11.pdf}
\end{center}

{\footnotesize {\bf Fig.~5}~Evolution of $Y_{N_1}$ and $Y_L$ obtained as
solutions of Boltzmann equations.  Results of the case (I) is shown
in the left panel for $\kappa_{\phi\sigma}=10^{-4}$ and
$\kappa_{\phi S}=\kappa_{\eta S}=0$.
Results of the the case (II) is shown in the right panel for
$\kappa_{\phi S}=\kappa_{\eta S}=10^{-6}$
and $\kappa_{\phi\sigma}=\kappa_{\eta\sigma}=0$.
Other parameters in each case are given in the text.
Initial values for them are fixed as $Y_{N_1}=Y_L=0$ at $z=z_R$.
$\rho_{N_1}/\rho_{\rm R}$ represents a ratio of the energy density
of $N_1$ to the one of radiation. }
\end{figure}

Solutions of the Boltzmann equations in  the cases (I) and (II) are
are shown in Fig.~5.  The lightest right-handed neutrino mass in each case is
$M_{N_1}=7\times 10^6$ and $10^8$~GeV.
In both cases, the sufficient baryon number asymmetry is found to be produced.
The figure for the case (I) shows clearly that the present scenario works well.
$Y_{N_1}$ reaches a value near $Y_{N_1}^{\rm eq}$ through the scattering
of the extra fermions as expected. Substantial out-of-equilibrium decay occurs
 at $z>10$ to generate the lepton number asymmetry.
The delay of the decay due to the small $h_1$ could make the washout of
the lepton number asymmetry ineffective.
On the other hand,  we cannot definitely find a signature of the scenario
in the case (II), where the $N_1$ mass is near the bound for which the usual
leptogenesis can generate the required baryon number asymmetry in the original
scotogenic model \cite{ks}. The figure shows that additional contribution to the
$N_1$ production starts at $z\simeq 0.1$.
It is considered to be brought about by the $N_1$ inverse decay since it is
expected to become effective
around $z\sim \left(\frac{6.3\times 10^{-5}}{h_1}\right)
\left(\frac{M_{N_1}}{10^8~{\rm GeV}}\right)$.
The figure shows that it plays a main role for the $N_1$ production.
As a result, the leptogenesis in this case results in the ordinary one where the
lower mass bound of $N_1$ is $O(10^8)$ GeV.

These results show that the model with suitable parameters can
generate a sufficient amount of baryon number asymmetry through
leptogenesis even if the reheating temperature is lower than $10^9$~GeV
as long as $M_{N_1}<T_R$ GeV is satisfied.
In the present model, both the right-handed neutrino mass and
the extra fermion mass
are determined by $M_{N_j}=y_{N_j}w$ and $\tilde M_F=\delta y_Fw$.
Since $w$ is fixed as the $PQ$ symmetry breaking scale and
$\delta>1$ is imposed for the realization of the substantial $CP$ phases in the CKM
and PMNS matrices, their mass cannot be arbitrarily smaller than $10^9$ GeV
under the condition that the scattering of the extra fermions to the right-handed
neutrinos is effective.
Because of this reason, low scale leptogenesis, which can be a distinguishable feature
of the model, tends to be allowed only for the case where the PQ symmetry breaking
occurs at a neighborhood of its lower bound.
Even in that case, successful leptogenesis is expected to be realized only in the range
 $M_{N_1}>4\times10^6\left(\frac{y_F}{10^{-1.2}}\right)^{-1/2}$ GeV for the case
$T_R>10^8$ GeV since $T$ in eq.~(\ref{scat}) should satisfy $T>\tilde M_F, M_{N_1}$
for the sufficient $N_1$ production.

\begin{figure}
\begin{center}
\includegraphics[width=7.5cm]{rt-basym.pdf}
\hspace*{5mm}
\includegraphics[width=7.5cm]{l5-basym.pdf}
\end{center}
\vspace*{-5mm}
\footnotesize{{\bf Fig.~6}~~Left : baryon number asymmetry generated
at the expected reheating temperature $T_R$ (GeV) in the case (I).
It is plotted by a red solid line for
$y_D=10^{-1.2}$ and by green crosses for $y_D=10^{-1.5}$.
Right : baryon number asymmetry generated at $T_R=10^9$ GeV for various
values of $\tilde\lambda_5$ and $y_{N_1}$ in the case (I).
In both panels, other parameters are fixed at the ones given in (\ref{para}). }
\end{figure}

In Fig. 6, we show the baryon number asymmetry $Y_B$ generated  in the case (I)
varying the values of relevant parameters. In the left panel, $Y_B$ is plotted
as a function of the reheating temperature, which is fixed by the inflaton composition and
its coupling with $\phi$ and $\eta$. Two values of $y_F$ are used in this plot.
The $N_1$ production in the present scenario depends on the reheating
temperature and the couplings $y_F,~ y_{N_1}$.
A red solid line representing $Y_B$ is expected for
the parameters given in eq.(\ref{para}) for the case (I).
It becomes larger and reaches an upper bound $Y_B\simeq 2.5\times 10^{-10}$
when the reheating temperature increases to $10^{10}$~GeV.
This behavior can be understood if we take into account that
the equilibrium number density of extra fermions are suppressed by
the Boltzmann factor at lower reheating temperature and then
the $N_1$ production due to the scattering of extra fermions is suppressed.
We also plot $Y_B$ for a smaller value of $y_F$ by green crosses at some typical
$T_R$. They show that $Y_B$ takes smaller values for a smaller $y_F$
since the $N_1$ production cross section is proportional
to $y_F^2$.\footnote{Since the effect of Boltzmann suppression
caused by its mass $\tilde M_F=\delta y_Fw$ could be dominant at lower $T_R$
compared with the effect on the cross section, the smaller $y_F$ gives a larger
$Y_B$ at $T_R<10^8$ GeV in this case.}
In the right panel,  $Y_B$ is plotted by varying
$|\tilde\lambda_5|$ and $y_{N_1}$.  A red solid line represents it as a function of
$|\tilde\lambda_5|$ for a fixed $y_{N_1}=5\times 10^{-3}$.
Since the neutrino oscillation data have to
be imposed on eq.~(\ref{oscilcon}), Yukawa couplings $h_{2,3}$ are settled
by $|\tilde\lambda_5|$, $y_{N_2}$ and $y_{N_3}$.
The $CP$ asymmetry $\varepsilon$ and the washout of the generated
lepton number asymmetry are mainly determined by $h_{2,3}$ for the fixed
$y_{N_2}$ and $y_{N_3}$ as found from eq.(\ref{epsilon}).
Since a smaller $|\tilde\lambda_5|$ makes $h_{2,3}$ larger and then both $\varepsilon$
and washout larger, $Y_B$ takes a maximum value for a certain $|\tilde\lambda_5|$,
which is found in the figure.
We also plot $Y_B$ by varying $y_{N_1}$ for a fixed $|\tilde\lambda_5|$ in the same panel.
A smaller $y_{N_1}$ makes the $N_1$ production less effective for a fixed
$y_F$ and then its lower bound is expected to appear for successful leptogenesis.
It gives the lower bound of $M_{N_1}$ as $\sim 4\times 10^6$ GeV
as predicted above.

Although other parameters are fixed at the ones given in (\ref{para})
in these figures, it is useful to give remarks on their dependence here.
If $\delta$ takes a larger value, the mass of extra fermions $\tilde M_F$
becomes larger to suppress the reaction density $\gamma_F$
due to the Boltzmann factor.
As a result, the $N_1$ number density generated through
the scattering becomes smaller and the resulting $Y_B$
also becomes smaller. If $h_1$ is much smaller, $N_1$ decay delays and
the entropy produced through the decay of relic $N_1$ might dilute
the generated lepton number asymmetry.

\subsection{Dark matter and isocurvature fluctuations}
This model has two DM candidates.
One  is the lightest neutral component of $\eta$ with $Z_2$ odd parity which is
an indispensable ingredient of the model.
It is known to be a good DM candidate which does not cause any contradiction
with known experimental data as long as its mass is in the TeV range
where the coannihilation can be effective \cite{idm,highmass,ks}.
As found from Fig.~1, both the DM abundance and the DM direct
search bound can be satisfied if the couplings $\tilde\lambda_3$ and $|\lambda_4|$
take suitable values of $O(1)$.
Although these parameters could affect the perturbativity of the scalar quartic
couplings through the radiative corrections, we can safely escape such problems
in certain parameter regions.
The results obtained for the case (I) in the previous part
are derived by supposing that the required DM is $\eta_R$.

Axion is another promising candidate in the model.
However, the axion could be a dominant component of DM only for
$f_a\sim 10^{11}$~GeV although it depends on the contribution from
the axion string decay \cite{axiondm}.
We consider the case (II) as such an example.
As described before, the PQ symmetry is spontaneously broken
during the inflation since the inflaton contains the radial component of $\sigma$.
As a result, the axion appears as the phase $\theta$ of $\sigma$.
Since the axion potential is flat during the inflation, the axion gets
a quantum fluctuation $\delta A=\left(H/2\pi\right)^2$ and it
can cause isocurvature fluctuation in the CMB amplitude \cite{isocurv,isocurv2}.
A canonically normalized axion $A$ is defined by noting eq.~(\ref{canoninf}) as
\begin{equation}
\frac{\partial A}{\partial \theta}=\frac{\tilde\sigma}{\Omega^2}\sqrt{\Omega^2+
6\xi_\sigma\frac{\tilde\sigma^2}{M_{\rm pl}^2}}
\simeq\frac{\sqrt{\chi M_{\rm pl}}}{\tilde\xi_S^{1/4}}
\frac{|\kappa_{\sigma S}|}{\tilde\kappa_\sigma} \equiv\chi_{\rm iso}.
\end{equation}
Since the axion interacts with other fields very weakly, it causes the isocurvature
fluctuation as the fluctuation of its number density $n_A$.
The amplitude of its power spectrum can be expressed as
\begin{equation}
{\cal P}_i(k)=\Big\langle\left|\frac{\delta n_A}{n_A}\right|^2\Big\rangle
=\frac{H_k^2}{\pi^2\chi_{\rm iso}^2\langle \theta^2\rangle}.
\end{equation}
Since the axion is only a source of the isocurvature fluctuation in this model,
its fraction in the power spectrum is given as
\begin{equation}
\alpha=\frac{R_a^2{\cal P}_i(k)} {R_a^2{\cal P}_i(k)+{\cal P}_s(k)}
\simeq 8~\epsilon~\tilde\xi_S^{1/2}\frac{M_{\rm pl}}{\chi_k}
\frac{R_a^2}{\langle\theta^2\rangle}\left(\frac{\tilde\kappa_\sigma}{\kappa_{\sigma S}}\right)^2,
\end{equation}
where ${\cal P}_s(k)=A_s$ which is given in eq.(\ref{power}).
$R_a$ is a fraction of the axion energy density in the CDM and defined as
$R_a=\Omega_a/\Omega_{\rm CDM}$.
If we use a relation \cite{qcdaxion}
\begin{equation}
R_a=\frac{\langle\theta^2\rangle}{6\times 10^{-6}}\left(\frac{f_a}{10^{16}~{\rm GeV}}\right)^{7/6},
\end{equation}
we find
\begin{equation}
\alpha=3.25\times 10^{-5}~\tilde\xi_S^{1/2}\frac{M_{\rm pl}}{\chi_k}
\left(\frac{55}{{\cal N}_k}\right)^2
\left(\frac{f_a}{10^{10}~{\rm GeV}}\right)^{7/6}R_a
\left(\frac{\tilde\kappa_\sigma}{\kappa_{\sigma S}}\right)^2.
\end{equation}
Since the Planck data constrain $\alpha$ as $\alpha\le 0.037$ at
$k=0.05$~Mpc$^{-1}$ \cite{planck18},
we have a condition for the model to be consistent with the present observation of
the isocurvature fluctuation in the CMB as
\begin{equation}
R_a<\frac{67}{\tilde\xi_S^{1/2}}\frac{\chi_k}{M_{\rm pl}}\left(\frac{{\cal N}_k}{55}\right)^2
\left(\frac{10^{11}~{\rm GeV}}{w}\right)^{7/6}
\left(\frac{\kappa_{\sigma S}}{\tilde\kappa_\sigma}\right)^2,
\end{equation}
where $f_a=w$ is used.
In the case (I), this gives no constraint and the parameters used in the present study
to estimate the reheating temperature in Fig.~4 and the baryon number asymmetry
are consistent with the observational data.
DM can be identified with the neutral component of $\eta$.
On the other hand, in the case (II), the isocurvature condition can be satisfied
for $R_a<0.21$ if $|\kappa_{\sigma S}|/\tilde\kappa_\sigma=10^{-1.6}$ is assumed
and for $R_a<0.034$ if $|\kappa_{\sigma S}|/\tilde\kappa_\sigma=10^{-2}$ is assumed.
The isocurvature constraint forbids the axion to be a dominant DM component and
the neutral component of $\eta$ is required to play a role of the dominant DM
in this case also.

\section{Summary}
We have proposed a model which could give an explanation for
the origin of the $CP$ phases in both the CKM and PMNS matrices
and the strong $CP$ problem.
It is a simple extension of the SM with vector-like extra fermions and several
scalars. In order to control the couplings of new fields, global symmetry is imposed.
If the $CP$ symmetry is spontaneously broken in a singlet scalar sector at
an intermediate scale, it can be transformed to the CKM and PMNS matrices
through the mixing between the extra fermions and the ordinary quarks
or the charged leptons.
On the other hand, since the colored extra fermions play the same role
as the ones in the KSVZ model for the strong $CP$ problem,
the strong $CP$ problem could be solved through the PQ mechanism.
After the symmetry breaking due to the singlet scalars, the leptonic sector of
the model is reduced to the scotogenic model, which can explain the small
neutrino masses and the DM abundance due to the remnant discrete symmetry
of the imposed symmetry.
Singlet scalars introduced to explain the $CP$ issues can play a role of inflaton
if it has a nonminimal coupling with the Ricci scalar.
We suppose this coupling is of order one. In that case,
although it gives the similar prediction for the scalar spectral index and
the tensor to scalar ratio to the one of the Higgs inflation,
reheating phenomena is different from it since the radiation domination starts
just after the end of inflation.

The model has a notable phenomenological feature in addition to these.
The extra fermions which are introduced for the $CP$ issues could make the
thermal leptogenesis generate the sufficient baryon number asymmetry
even if the lightest right-handed neutrino mass is much lower than $10^9$~GeV,
which is the well-known lower bound of the right-handed neutrino mass for
successful leptogenesis in the ordinary seesaw scenario.
Although the model allows low scale leptogenesis, it is difficult to
distinguish it from other thermal leptogenesis model experimentally.
However, if we consider its supersymmetric extension, it could give a possibility
to escape the gravitino problem.
The model is constrained by the isocurvature fluctuation which is caused
by the spontaneous breaking of the PQ symmetry during the inflation.
We find that its present observation can be consistent with the model even if
the DM relic abundance is imposed on the model.
Although the relic density of axion should be a small fraction of the DM,
there is a neutral component of the inert doublet scalar as an alternative
candidate of DM in the model.
It can explain the DM abundance just as in the scotogenic model without
affecting other predictions of the model.
It is remarkable that the model has potentiality to explain various issues
in the SM although the model is rather simple.
It may deserve further study.

\section*{Appendix A~~ A simple example for $A_f$}
In this Appendix, we present a simple example which could bring
about a phase in the CKM matrix.
In this example, we assume $w=10^9$ GeV and $u = 10^{11}$ GeV,
and also the relevant Yukawa couplings $h_d$, $y_d$ and
$\tilde y_d$ to be written by using real constant parameters as
\begin{equation}
h_d=c\left(\begin{array}{ccc}
\epsilon^4 &\epsilon^3 &p_1\epsilon^3\\
\epsilon^3 &\epsilon^2 & p_2\epsilon^2 \\
\epsilon^2 & p_3 & -p_3 \\
\end{array}\right), \quad
y_d=(0, y, 0), \quad \tilde y_d=(0, 0, \tilde y).
\label{yd}
\end{equation}
As long as $\epsilon$ satisfies $\epsilon < 1$,
the down type quark mass matrix $m_d(\equiv h_d\langle\tilde\phi\rangle)$
has hierarchical mass eigenvalues.
Here, we introduce $X_{ij}$ and $Y_{ij}$ whose definition is given as
\begin{eqnarray}
&&X_{ij}=1+p_ip_j
+\frac{y^2+\tilde y^2p_ip_j
+y\tilde y(p_i+p_j)\cos 2\rho}{y^2+\tilde y^2}
\left(1-\frac{1}{\delta^2}\right), \nonumber \\
&&Y_{ij}=\frac{y\tilde y(p_i-p_j)\sin 2\rho}{y^2+\tilde y^2}
\left(1-\frac{1}{\delta^2}\right),
\end{eqnarray}
where $\delta$ is defined as $\delta=\tilde M_F/\mu_F$.
If we define $R_{ij}$ and $\theta_{ij}$ by using these quantities as
\begin{equation}
R_{ij}=\sqrt{X_{ij}^2+Y_{ij}^2}, \quad \tan\theta_{ij}=\frac{Y_{ij}}{X_{ij}},
\end{equation}
the component of eq.~(\ref{ckm}) is found to be expressed under the assumption
$\mu_D^2< {\cal F}_d{\cal F}_d^\dagger$ as
\begin{equation}
(A_d^{-1}m^2A_d)_{ij}=c^2\langle\tilde\phi\rangle^2 \epsilon_{ij}R_{ij}e^{i\theta_{ij}},
\label{diag}
\end{equation}
where $\epsilon_{ij}$ is defined as
\begin{equation}
\epsilon_{11}=\epsilon^6, \quad \epsilon_{22}=\epsilon^4, \quad
\epsilon_{33}=1, \quad \epsilon_{12}=\epsilon_{21}=\epsilon^5, \quad
\epsilon_{13}=\epsilon_{31}=\epsilon^3, \quad
\epsilon_{23}=\epsilon_{32}=\epsilon^2.
\end{equation}

By solving eq.~(\ref{diag}), we find that $A_d$ is approximately written as
\begin{equation}
A_d\simeq\left(
\begin{array}{ccc}
1 & -\lambda & \lambda^3\left(\frac{X_{23}}{|\alpha|^2X_{33}}e^{i\vartheta}
-\frac{X_{13}}{|\alpha|^3X_{33}}\right) \\
\lambda & 1 & -\lambda^2 \frac{X_{23}}{|\alpha|^2X_{33}}e^{i\vartheta}\\
 \lambda^3\frac{X_{13}}{|\alpha|^3X_{33}}&
 \lambda^2\frac{X_{23}}{|\alpha|^2X_{33}}e^{-i\vartheta} & 1 \\
\end{array}
\right),
\end{equation}
where the constants $\lambda$, $\alpha$, and $\vartheta$ are defined by
\begin{equation}
\alpha=\frac{X_{12}X_{33}-X_{13}X_{23}
e^{-i(\theta_{23}+\theta_{12}-\theta_{13})}}
{X_{22}X_{33}-X_{23}^2}, \quad
\lambda=|\alpha|\epsilon,\quad
\vartheta=\arg(\alpha) +\theta_{23}+\theta_{12}-\theta_{13}.
\end{equation}
This expression shows that $A_d$ could have a nontrivial phase which
gives the origin of the CKM phase. If the diagonalization matrix $O^L$ for a
mass matrix of the up type quarks takes an almost diagonal form,
the CKM matrix could be obtained as $V_{\rm CKM}\simeq A_d$.
As an example, if we assume $\cos\rho=\frac{\pi}{4}$ and fix other
parameters as
\begin{eqnarray}
&&y_F=10^{-1.5}, \quad \delta=\sqrt{3}, \quad y=4\times 10^{-4}, \quad
 \tilde y=2\times 10^{-4},   \nonumber \\
&&p_1=1.1, \quad p_2=-0.9, \quad p_3=1,  \quad \epsilon=0.2, \quad c=0.014,
\end{eqnarray}
then we obtain $\lambda=0.22$ and the Jarlskog invariant \cite{jar}
as $J(\equiv {\rm Im}[ A_{12}A_{13}^\ast A_{23} A_{22}^\ast])=-1.6\times 10^{-6}$.
The mass eigenvalues for the down type quarks are obtained as
\begin{eqnarray}
&&m_d=\left\{X_{11}-\frac{X_{13}^2}{X_{33}}+|\alpha|^2\left(X_{22}
-\frac{X_{23}^2}{X_{33}}-2\right)\right\}^{1/2}\epsilon^3c\langle\tilde\phi\rangle
\simeq 3.3~{\rm MeV}, \nonumber \\
&&m_s=\left(X_{22}-\frac{X_{23}^2}{X_{33}}\right)^{1/2}\epsilon^2c\langle\tilde\phi\rangle
\simeq 138~{\rm MeV}, \nonumber\\
&&m_b=X_{33}^{1/2}c\langle\tilde\phi\rangle\simeq 4.2~{\rm GeV}.
\end{eqnarray}
Although a diagonalization matrix $A_e$ for the charged lepton sector
may be considered to take the same form as $A_d$, it is not favorable
for a large $CP$ phase
in the PMNS matrix. In that case, since $A_e$ is the nearly diagonal, large flavor mixing
has to be caused only by the neutrino sector and the Dirac $CP$ phase in the PMNS
matrix becomes small as a result. A large $CP$ phase in the PMNS matrix requires
$A_e$ to have rather large off-diagonal elements also.

\section*{Appendix B~~ RGEs for coupling constants}
In order to examine both the stability and the perturbativity of the model
from the weak scale to the Planck scale, we have to know the running of the
coupling constants in eq.~(\ref{model}).
If we fix these coupling constants at an intermediate scale
and solve the RGEs to the Planck scale, we can find their values
throughout the scale. The one-loop RGEs of the relevant coupling constants
are given as

\begin{eqnarray}
&&16\pi^2\mu\frac{\partial \lambda_1}{\partial\mu}=24\lambda_1^2
+\lambda_3^2+(\lambda_3+\lambda_4)^2 + \kappa_{\phi\sigma}^2
+ \kappa_{\phi S}^2
+\frac{3}{8}\left(3g^4+g^{\prime 4}+2g^2g^{\prime 2}\right) \nonumber \\
&&\hspace*{20mm}-3\lambda_1\left(3g^2+g^{\prime 2}-4h_t^2\right)-6h_t^4,  \nonumber\\
&&16\pi^2\mu\frac{\partial \lambda_2}{\partial\mu}=24\lambda_2^2+\lambda_3^2
+(\lambda_3+\lambda_4)^2 + \kappa_{\eta\sigma}^2
+ \kappa_{\eta S}^2
+\frac{3}{8}\left(3g^4+g^{\prime 4}+2g^2g^{\prime 2}\right)
 \nonumber \\
&&\hspace*{20mm}-3\lambda_2\left(3g^2+g^{\prime 2}\right) +4\lambda_2\left(2h_2^2+3h_3^2\right)-8h_2^4-18h_3^4, \nonumber\\
&&16\pi^2\mu\frac{\partial \lambda_3}{\partial\mu}=2(\lambda_1+\lambda_2)
(6\lambda_3+2\lambda_4)
+4\lambda_3^2+2\lambda_4^2 +2\kappa_{\phi\sigma}\kappa_{\eta\sigma}+
 2\kappa_{\phi S}\kappa_{\eta S}\nonumber \\
&&\hspace*{20mm}+\frac{3}{4}\left(3g^4+g^{\prime 4}-2g^2g^{\prime 2}\right)
-3\lambda_3\left(3g^2+g^{\prime 2}\right)
+2\lambda_3\left(3h_t^2+2h_2^2+3h_3^2\right), \nonumber \\
&&16\pi^2\mu\frac{\partial \lambda_4}{\partial\mu}=4(\lambda_1+\lambda_2)\lambda_4+8\lambda_3\lambda_4+4\lambda_4^2
+3g^2g^{\prime 2}-3\lambda_4\left(3g^2+g^{\prime 2}\right) \nonumber \\
&&\hspace*{20mm}+2\lambda_4\left(3h_t^2+2h_2^2+3h_3^2\right),  \nonumber \\
&&16\pi^2\mu\frac{\partial \kappa_S}{\partial\mu}=20\kappa_S^2
+\kappa_{\sigma S}^2+2(\kappa_{\phi S}^2+\kappa_{\eta S}^2)
+4\kappa_S\left[3(y_{d}^2+\tilde y_{d}^2)+y_{e}^2+\tilde y_{e}^2\right] \nonumber \\
&&\hspace*{20mm}-2\left[3(y_{d}^2+\tilde y_{d}^2)^2)
+(y_{e}^2+\tilde y_{e}^2)^2\right], \nonumber \\
&&16\pi^2\mu\frac{\partial \kappa_\sigma}{\partial\mu}=20\kappa_\sigma^2
+\kappa_{\sigma S}^2+2(\kappa_{\phi \sigma}^2+\kappa_{\eta \sigma}^2)
+4\kappa_\sigma\left(3y_D^2+y_E^2+\frac{1}{2}y_{N_3}^2\right)-
2\left(3y_D^4+y_E^4+\frac{1}{2}y_{N_3}^4\right), \nonumber \\
&&16\pi^2\mu\frac{\partial \kappa_{\sigma S}}{\partial\mu}= 4\kappa_{\sigma S}^2+
8(\kappa_S+\kappa_{\sigma})\kappa_{\sigma S}+
2(\kappa_{\phi S}\kappa_{\phi\sigma}
+\kappa_{\eta S}\kappa_{\eta\sigma})\nonumber \\
&&\hspace*{20mm}+2\kappa_{\sigma S}\left[3(y_{d}^2
+\tilde y_{d}^2)+y_{e}^2+\tilde y_{e}^2
+3y_D^2+y_E^2+\frac{1}{2}y_{N_3}^2\right]
-4[3y_D^2(y_d^2+\tilde y_d^2)+y_E^2(y_e^2+\tilde y_e^2)], \nonumber \\
&&16\pi^2\mu\frac{\partial \kappa_{\phi\sigma}}{\partial\mu}=
4\kappa_{\phi\sigma}^2+
2\kappa_{\sigma S}\kappa_{\phi S}+2\kappa_{\eta\sigma}(2\lambda_3+\lambda_4)
+4\kappa_{\phi\sigma}(3\lambda_1+2\kappa_\sigma)\nonumber \\
&&\hspace*{20mm}+2\kappa_{\phi\sigma}\left(3y_D^2+y_E^2+\frac{1}{2}y_{N_3}^2+
3h_t^2-\frac{9}{4}g^2
-\frac{3}{4}g^{\prime 2}\right), \nonumber \\
&&16\pi^2\mu\frac{\partial \kappa_{\eta\sigma}}{\partial\mu}=4\kappa_{\eta\sigma}+
2\kappa_{\sigma S}\kappa_{\eta S}+2\kappa_{\phi\sigma}(2\lambda_3+\lambda_4)
+4\kappa_{\eta\sigma}(3\lambda_2+2\kappa_\sigma) \nonumber \\
&&\hspace*{20mm}+2\kappa_{\eta\sigma}\left(3y_D^2+y_E^2+\frac{1}{2}y_{N_3}^2 +2h_2^2+3h_3^2
 -\frac{9}{4}g^2-\frac{3}{4}g^{\prime 2}\right)-4(h_2^2+h_3^2)y_{N_3}^2 \nonumber \\
&&16\pi^2\mu\frac{\partial \kappa_{\phi S}}{\partial\mu}=4\kappa_{\phi S}^2+
2\kappa_{\sigma S}\kappa_{\phi \sigma}+2\kappa_{\eta S}(2\lambda_3+\lambda_4)
+4\kappa_{\phi S}(3\lambda_1+2\kappa_S) \nonumber \\
&&\hspace*{20mm}+2\kappa_{\phi S}\left[3(y_{d}^2+\tilde y_{d}^2)
+y_{e}^2+\tilde y_{e}^2+3h_t^2 -\frac{9}{4}g^2-\frac{3}{4}g^{\prime 2}\right],
\nonumber
\end{eqnarray}
\begin{eqnarray}
&&16\pi^2\mu\frac{\partial \kappa_{\eta S}}{\partial\mu}=4\kappa_{\eta S}^2+
2\kappa_{\sigma S}\kappa_{\eta\sigma}+2\kappa_{\phi S}(2\lambda_3+\lambda_4)+
4\kappa_{\eta S}(3\lambda_2+2\kappa_S)\nonumber\\
&&\hspace*{20mm}+2\kappa_{\eta S}\left[3(y_{d}^2+\tilde y_{d}^2)
+y_{e}^2+\tilde y_{e}^2+2h_2^2+3h_3^2
-\frac{9}{4}g^2-\frac{3}{4}g^{\prime 2}\right], \nonumber \\
&&16\pi^2\mu\frac{\partial y_{d}}{\partial\mu}=y_{d}\left[-8g_s^2
-\frac{2}{3}g^{\prime 2}+\frac{1}{2}y_D^2
+4y_{d}^2+3\tilde y_{d}^2+y_{e}^2+\tilde y_{e}^2\right], \nonumber \\\
&&16\pi^2\mu\frac{\partial\tilde y_{d}}{\partial\mu}=\tilde y_{d}
\left[-8g_s^2-\frac{2}{3}g^{\prime 2}+\frac{1}{2}y_D^2
+3y_{d}^2+4\tilde y_{d}^2+y_{e}^2+\tilde y_{e}^2\right], \nonumber \\
&&16\pi^2\mu\frac{\partial y_{e}}{\partial\mu}=y_{e}\left[-6g^{\prime 2}+\frac{1}{2}y_E^2
+3(y_{d}^2+\tilde y_{d}^2)+2y_{e}^2+\tilde y_{e}^2\right], \nonumber \\
&&16\pi^2\mu\frac{\partial \tilde y_{e}}{\partial\mu}=\tilde y_{e}\left[-6g^{\prime 2}+\frac{1}{2}y_E^2
+3(y_{d}^2+\tilde y_{d}^2)+y_{e}^2+2\tilde y_{e}^2\right], \nonumber \\
&&16\pi^2\mu\frac{\partial y_D}{\partial\mu}=y_D\left(-8g_s^2-\frac{2}{3}g^{\prime 2}
+4y_D^2+y_E^2+\frac{1}{2}y_d^2+\frac{1}{2}\tilde y_d^2+\frac{1}{2}y_{N_3}^2\right), \nonumber \\
&&16\pi^2\mu\frac{\partial y_E}{\partial\mu}=y_E\left(-6g^{\prime 2}
+3y_D^2+2y_E^2+\frac{1}{2}y_e^2+\frac{1}{2}\tilde y_e^2+\frac{1}{2}y_{N_3}^2\right), \nonumber \\
&&16\pi^2\mu\frac{\partial y_{N_3}}{\partial\mu}=y_{N_3}\left[3y_D^2+y_E^2
+\frac{3}{2}y_{N_3}^2+2(h_2^2+h_3^2)\right], \nonumber \\
&&16\pi^2\mu\frac{\partial h_2}{\partial\mu}=h_2\left(-\frac{9}{4}g^2-
\frac{3}{4}g^{\prime 2}+5h_2^2+3h_3^2\right),
 \nonumber\\
&&16\pi^2\mu\frac{\partial h_3}{\partial\mu}=h_3\left(-\frac{9}{4}g^2-
\frac{3}{4}g^{\prime 2}+2h_2^2+\frac{15}{2}h_3^2+\frac{1}{2}y_{N_3}^2\right),
 \nonumber\\
&&16\pi^2\mu\frac{\partial h_t}{\partial\mu}=h_t\left(\frac{9}{2}h_t^2-8g_s^2-\frac{9}{4}g^2
-\frac{17}{12}g^{\prime 2}\right),  \nonumber\\
&&16\pi^2\mu\frac{\partial g_s}{\partial\mu}=-\frac{19}{3}g_s^3, \qquad
16\pi^2\mu\frac{\partial g}{\partial\mu}=-3g^3, \qquad
16\pi^2\mu\frac{\partial g^\prime}{\partial\mu}=\frac{79}{9}g^{\prime 3},
\end{eqnarray}
where $g_s$, $g$, and $g^\prime$ are the gauge coupling constants of the SM.
In these equations, we assume eq.~(\ref{tribi}) for neutrino Yukawa couplings
and eq.~(\ref{yd}) for $y_f$ and $\tilde y_f$, and also only
$y_{N_3}$ is taken into account. Contributions from $V_b$ are also neglected
in these RGEs.

Initial values of the coupling constants which are used in the RGE study of the right panel
of Fig.~1 are taken to be the same ones which are used in the analysis of the leptogenesis
for the case (I).
A part of them are fixed at the weak scale as
\begin{eqnarray}
&&\tilde\lambda_1=0.13, \quad \tilde\lambda_2=0.1,  \quad \tilde\lambda_3=0.445, \quad
\lambda_4=-0.545, \quad \lambda_5=-10^{-3}, \nonumber \\
&&h_1= 6\times 10^{-7}, \quad h_2= 8.3\times 10^{-3} , \quad  h_3=3.9\times 10^{-3},
\end{eqnarray}
where $\tilde\lambda_{3}$ and $\lambda_4$ are fixed by the DM constraint shown
in the left panel of Fig.~1 and
$h_{2,3}$ are fixed by using the nautrino mass formula and the neutrino oscillation data.
Remaining ones are fixed at a scale $\bar M$ as
\begin{eqnarray}
&&\tilde\kappa_S=10^{-8},   \quad \tilde\kappa_\sigma=10^{-4.5}, \quad  |\kappa_{\sigma S}|=10^{-6.1},
\quad \kappa_{\phi\sigma}=\kappa_{\eta\sigma}=10^{-4},  \quad
\kappa_{\phi S}=\kappa_{\eta S}=10^{-7}\nonumber \\
&& y_{N_1}=7\times 10^{-3}. \quad y_{N_2}=2\times 10^{-2}, \quad y_{N_3}=4\times 10^{-2},
\nonumber\\
&&y_D=y_E=10^{-1.2}, \quad y_e=y_d=8\times 10^{-4}, \quad \tilde y_e=\tilde y_d=4\times 10^{-4}.
\end{eqnarray}
These satisfy the imposed conditions (\ref{stab1}), (\ref{uw}), (\ref{stab2}), (\ref{ks}) and (\ref{kap}).

\section*{Acknowledgements}
This work is partially supported by a Grant-in-Aid for Scientific Research (C)
from Japan Society for Promotion of Science (Grant No. 18K03644).
N.~S.~R is supported by \textit{Program 5000 Doktor} under Ministry of
Religious Affairs (MORA), Republic of Indonesia.

\newpage
\bibliographystyle{unsrt}