††institutetext: Institute of Theoretical Physics, School of Physics, Dalian
University of Technology, No.2 Linggong Road, Dalian, Liaoning, 116024,
P.R.China

# Inflation and Dark Matter in the $Z_{5}$ Model

XinXin Qi  Hao Sun<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We discuss the possibility of unifying dark matter physics and inflation in
the $Z_{5}$ model of the two-component dark matter. Inflation driven by the
two-component dark matter fields can be divided into two cases, singlet dark
matter inflation and mixed dark matter inflation, where both two-component
play the role of inflaton in the latter case. For dark matter, we focus on the
mixed dark matter inflation case. We show a viable parameter space that
satisfies the theoretical and dark matter relic density constraint in the case
of successful inflation. It turns out that the dark matter density is
dominated by the light component, which is consistent with the feature of the
$Z_{5}$ model of the two-component dark matter.

##  1 Introduction

The Standard Model (SM) with $\rm SU(3)_{C}\otimes SU(2)_{L}\otimes U(1)_{Y}$
gauge group has achieved unparalleled success in particle physics, thanks to
the discovery of the last missing piece in SM, the Higgs boson ATLAS:2012yve ;
CMS:2012qbp . However, when connecting SM phenomenology with that of modern
cosmology, there exist even more outstanding problems unsolved. For example,
the assumptions of inflation and dark matter, which are supported by the
precise measurement of Cosmic Microwave Background (CMB) in the Planck
Liu:2021diz  and the WMAP WMAP:2010qai  experiments.

For the issue of inflation, it works well in answering the horizon problem and
homogeneity problem in CMB. One general method of constructing an inflation
model is to introduce scalar fields coupled or uncoupled with the SM Higgs
field. Such an inflationary model can provide initial conditions for the
subsequent thermal big bang evolution of the universe. Another important
peculiarity is that it can also explain the temperature anisotropy in CMB and
large scale structures in the universe because the quantum fluctuations of the
inflation field are stretched to the cosmic scale during inflation.
Furthermore, the inflation parameters such as the spectral index, the tensor
to scalar ratio and the scalar power spectrum extracted from the experiment
Liu:2021diz ; WMAP:2010qai  can restrict the specific inflationary model when
the slow roll parameter $\epsilon\approx 1$. The simplest method is to select
the Higgs field in SM as the inflation field Bezrukov:2007ep ; Bezrukov:2010jz
without introducing extra scalar fields, and the only additional parameter
comes from the $\xi R|H|^{2}$ term. However, there are two problems in the
such model: one is that the quartic self-coupling of the Higgs field $\lambda$
becomes negative at high energy scales ($\gtrsim 10^{10}$ GeV), which destroys
the stability of the vacuum Sher:1988mj . The other is that at the scale
around $M_{Pl}/{\xi}\approx 10^{3}$ GeV Burgess:2009ea , the unitarity is
broken by requiring $\xi\sim 10^{4}$ due to the scalar power spectrum boundary
Lerner_2010 . To overcome the unitarity problem at the unitarity scale in the
Higgs inflation model, an additional dynamical degree of freedom should be
introduced Giudice:2010ka . In this case, the extra gauge singlet scalar
particle can act as the inflaton and the Higgs field in the SM will play the
role of reheating the universe. Generally speaking, people talk more about
chaotic inflation or stabinsky inflation. The former has a power law potential
of $m^{2}\phi^{2}+\lambda\phi^{4}$, and the latter has exponential potential.
We will focus on the latter case, and for more details on inflationary
cosmology, we refer to the review Linde:2014nna .

Dark Matter (DM) is another outstanding issue, which has been studied
extensively for many years. Recently, the Planck experiment Liu:2021diz  and
other astronomical observations put forward a new limit on the abundance of
DM, that is $\Omega_{DM}h^{2}\approx 0.12$. The most widely studied dark
matter scenario is the so-called weakly interacting massive particles (WIMP)
in which the mass of DM candidate particles ranges from GeV to TeV. Other DM
scenarios are also optional, including but not limited to axions, feebly
interacting massive particles (FIMP), strongly interacting massive particles
(SIMP) and forbidden dark matter Hall:2009bx ; JOUR ; Baer:2014eja ;
Hochberg:2014dra . As one of the simplest dark matter models, the gauge
singlet scalar dark matter model has been studied for a long time
Cline:2013gha ; Queiroz:2014yna . The more complex DM model with the scalar
dark matter usually carries additional discrete symmetry, which can be
regarded as a remnant of the spontaneously broken of some $U(1)$ gauge
symmetry, such as $Z_{2}$ and $Z_{N}$ Batell:2010bp ; Belanger:2014bga . On
the other hand, more studies have shown that there is no reason to assume that
the observed dark matter density is contributed just by a single new particle
Boehm:2003ha ; Ma:2006uv ; Cao:2007fy ; Hur:2007ur ; Lee:2008pc ; Zurek:2008qg
; Profumo:2009tb ; Baer:2011hx ; Esch:2014jpa . In Ref.Belanger:2020hyh , a
two-component dark matter model based on the $Z_{5}$ symmetry has been
studied. Here the $Z_{5}$ symmetry serves as the prototype of all $Z_{N}$
scenarios in which the DM particles are two complex scalars. Since the
particles carry different $Z_{5}$ charges and assume $M_{1}/2<M_{2}<2M_{1}$,
both the two DM particles are stable and can be selected as WIMP. More works
have recently attempted to combine inflation with dark matter in one
framework. In Refs. Kahlhoefer_2015 ; Lerner:2009xg ; Aravind:2015xst , a
gauge singlet scalar acts as the inflaton and then acts as a dark matter
candidate in the freeze-out scenario. Ref.Okada:2010jd  studied a scalar WIMP
dark matter candidate with non-minimal coupling to gravity, where the scalar
dark matter acts as the inflaton. Ref.Gong:2012ri  discussed inflation with
the two Higgs doublet model at low energy. In ref.Choubey:2017hsq , the
reheating after inflation was discussed in detail in the $Z_{2}$ dark matter
model Finally, a scenario linking inflation and dark matter in the model with
gauged $Z_{3}$ symmetry was discussed in ref.Choi:2020ara .

In our present paper, we connect the slow-roll inflation with the non-minimal
gravity couplings to the two-component dark matter model with $Z_{5}$
symmetry. Due to non-minimal gravity couplings, more than one scalar field
beyond the SM exist, which will lead to a pure field or multi-field inflation.
They can then play the role of DM in the universe during cooling down. We
discuss the theory constraints and inflation dynamics of the $Z_{5}$ dark
matter model.

Our paper is organized as follows. In Sec.2, we introduce the model framework.
In Sec.3, we study the inflation dynamics in the cases of singlet dark matter
inflation and mixed dark matter inflation. In Sec.4, we discuss the
constraints on the model, including unitarity, perturbativity and vacuum
stability. We discuss the DM phenomenology in the model with some fixed
parameters in Sec.5. Finally, we give our conclusion in Sec.6.

##  2 The $Z_{5}$ model of two-component dark matter for inflation

The general two-component dark matter field action under $Z_{5}$ symmetry with
the Jordan frame is given by

$\displaystyle\frac{{\cal L}_{J}}{\sqrt{-g_{J}}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}M^{2}_{Pl}R-\Big{(}\xi_{0}|H|^{2}+\xi_{1}|\Phi_{1}|^{2}+\xi_{2}|\Phi_{2}|^{2}\Big{)}R-\left|{D_{\mu}}H\right|^{2}-\left|D_{\mu}\Phi_{1}\right|^{2}-\left|D_{\mu}\Phi_{2}\right|^{2}$
(1) $\displaystyle-$ $\displaystyle{{\cal
V}_{J}}\left(H,\Phi_{1},\Phi_{2}\right)\ ,$

where the general renormalizable potential is:

$\displaystyle{{\cal V}_{J}}\left(H,\Phi_{1},\Phi_{2}\right)$
$\displaystyle=\,\,-\mu^{2}_{H}|H|^{2}+\lambda_{H}|H|^{4}+\mu_{1}^{2}|\Phi_{1}|^{2}+\lambda_{41}|\Phi_{1}|^{4}+\lambda_{S1}|H|^{2}|\Phi_{1}|^{2}$
$\displaystyle\,+\mu_{2}^{2}|\Phi_{2}|^{2}+\lambda_{42}|\Phi_{2}|{{}^{4}}+\lambda_{S2}|H|^{2}|\Phi_{2}|^{2}+\lambda_{412}|\Phi_{1}|^{2}|\Phi_{2}|^{2}$
$\displaystyle\,+\frac{1}{2}\left[\mu_{S1}\Phi^{2}_{1}\Phi_{2}^{*}+\mu_{S2}\Phi_{2}^{2}\Phi_{1}+\lambda_{31}\Phi_{1}^{3}\Phi_{2}+\lambda_{32}\Phi_{1}\Phi_{2}^{*3}+\text{h.c.}\right].$
(2)

The SM Higgs doublet and the components of the two dark matter complex
singlets are:

$H=\frac{1}{\sqrt{2}}\left(\begin{array}[]{l}0\\\
h\end{array}\right)\,,\quad\Phi_{1}=\frac{1}{\sqrt{2}}\phi_{1}e^{i\theta_{1}}\,,\quad\Phi_{2}=\frac{1}{\sqrt{2}}\phi_{2}e^{i\theta_{2}}\,.\quad$
(3)

In Eq.(1), $M_{Pl}$ represents the reduced Plank mass ($M_{Pl}\equiv
1/\sqrt{8\pi G_{N}}=2.44\times 10^{18}$ GeV), which will be set by 1 later but
will be restored if necessary. $R$ is the Ricci scalar. For the non-minimal
dimensionless couplings, we label them as $\xi_{i}$ ($i=0,1,2$) and require
their values to be positive to ensure potential stability under large field
values. $H$ is the SM Higgs doublet field and will obtain the vacuum
expectation value (vev) $v$ after spontaneously symmetry breaking. We consider
the vevs of the two-component dark matter fields to be zero under the $Z_{5}$
symmetry and the fields are labeled as $\Phi_{1}$, $\Phi_{2}$. $D_{\mu}$
stands for the covariant derivative and can be treated as the normal
derivative at the inflation epoch: $D_{\mu}\to\partial_{\mu}$. The minus sign
in the kinetic terms is needed to obtain the metric convention of $(-,+,+,+)$.
$\mathcal{V}_{J}$ is the general $Z_{5}$-invariant scalar potential
Belanger:2020hyh . We have the transform of $SM\to SM$, $\Phi_{1}\to
e^{i2\pi/5}\Phi_{1}$ and $\Phi_{2}\to e^{i4\pi/5}\Phi_{2}$ under the $Z_{5}$
symmetry. The new trilinear and quartic couplings denoted $\mu_{Si}$ and
$\lambda_{3i}~{}({i=1,2})$ respectively, can be found in the model. As we can
see, with the simultaneous disappearance of the $\mu_{Si}$ and $\lambda_{3i}$
terms, the scalar potential with $Z_{5}$ symmetry is indistinguishable from
that $Z_{2}\times Z^{{}^{\prime}}_{2}$ symmetry. In addition, by choosing the
phase of $\Phi_{1,2}$, the parameter $\lambda_{3i}$ can be real or complex.
For simplicity, we assume that all the parameters in Eq.(2) are real.

Assuming that the SM Higgs field $H$ rolls down rapidly to its vev during
inflation, as the initial condition or the effective mass is larger than the
Hubble scale, we can set $|H|=0$ and neglect all terms containing $H$. Then,
the inflation dynamics should be driven by $\Phi_{1}$ or $\Phi_{2}$, or both.

To eliminate the terms containing $R$ and $\Phi_{i}$ in Eq.(1), we perform a
conformal transformation from the Jordan frame to the Einstein frame. Note
that the parameters associated with conformal transformation will be equal to
1, making the Jordan and Einstein frame equivalent after the end of inflation.
The conformal transformation of the metric $g_{J}$ is given by Kaiser:2010ps :

$\displaystyle g^{J}_{\mu\nu}=g^{E}_{\mu\nu}/\Omega^{2}$ (4)

with

$\displaystyle\Omega^{2}\equiv 1+\xi_{1}\phi_{1}^{2}+\xi_{2}{\phi}_{2}^{2}.$
(5)

Then, we get the Einstein frame action:

$\frac{{\cal
L}_{E}}{\sqrt{-g_{E}}}=-\frac{R}{2}-\frac{1}{2}G_{ij}\partial_{\mu}\phi^{i}\partial^{\mu}\phi^{j}-{{\cal
V}_{E}}(\phi_{1},\phi_{2}),$ (6)

where:

$\displaystyle G_{ij}$
$\displaystyle=\frac{1}{\Omega^{2}}\delta_{ij}+\frac{3}{2}\frac{1}{\Omega^{4}}\frac{\partial\,\Omega^{2}}{\partial\,\phi_{i}}\frac{\partial\,\Omega^{2}}{\partial\,\phi_{j}},$
(7) $\displaystyle{\cal V}_{E}$ $\displaystyle=\frac{{\cal
V}_{J}}{\Omega^{4}}.$ (8)

Here, $G_{i,j}$ is a $2\times 2$ matrix, ${\cal V}_{E}$ denotes the scalar
potential in Einstein frame, the phase $\theta_{1}$ and $\theta_{2}$ can set
to be zero by minimizing the scalar potential at $\phi_{1}\neq 0$,
$\phi_{2}\neq 0$. In the large-field limit
$\xi_{1}\phi^{2}_{1}+\xi_{2}\phi^{2}_{2}\gg 1$, we can drop the scalar masses
and the cubic terms, then $G_{i,j}$ and ${\cal V}_{E}$ are unfolded as:

$\displaystyle G_{i,j}=$
$\displaystyle\frac{1}{1+\xi_{1}\phi^{2}_{1}+\xi_{2}\phi^{2}_{2}}\begin{pmatrix}1+\dfrac{6\xi^{2}_{1}\phi^{2}_{1}}{1+\xi_{1}\phi^{2}_{1}+\xi_{2}\phi^{2}_{2}}&\dfrac{6\xi_{1}\xi_{2}\phi_{1}\phi_{2}}{1+\xi_{1}\phi^{2}_{1}+\xi_{2}\phi^{2}_{2}}\\\
\dfrac{6\xi_{1}\xi_{2}\phi_{1}\phi_{2}}{1+\xi_{1}\phi^{2}_{1}+\xi_{2}\phi^{2}_{2}}&1+\dfrac{6\xi^{2}_{2}\phi^{2}_{2}}{1+\xi_{1}\phi^{2}_{1}+\xi_{2}\phi^{2}_{2}}\\\
\end{pmatrix}\,,$ (9) $\displaystyle{\cal V}_{E}=$
$\displaystyle\frac{\lambda_{41}\phi_{1}^{4}+\lambda_{42}\phi_{2}^{4}+\lambda_{412}\phi_{1}^{2}\phi_{2}^{2}-\lambda_{31}\phi_{1}^{3}\phi_{2}-\lambda_{32}\phi_{1}\phi_{2}^{3}}{4\left(1+\xi_{1}\phi^{2}_{1}+\xi_{2}\phi^{2}_{2}\right)^{2}}\,.$
(10)

With the non-zero mixed kinetic term $G_{1,2}$, the diagonal kinetic form can
be obtained by redefining

$\displaystyle\varphi=$
$\displaystyle\sqrt{\frac{3}{2}}\log(1+\xi_{1}\phi^{2}_{1}+\xi_{2}\phi^{2}_{2})\,,$
(11) $\displaystyle r=$ $\displaystyle\frac{\phi_{2}}{\phi_{1}}\,.$ (12)

The kinetic term and scalar potential in the Einstein frame can be rewritten
as:

$\displaystyle\mathcal{L}_{{\rm kin}}$
$\displaystyle=\frac{1}{2}\frac{(\xi_{1}^{2}+\xi_{2}^{2}r^{2})}{(\xi_{1}+\xi_{2}r^{2})^{3}}\left(1-e^{-\sqrt{\frac{2}{3}}\varphi}\right)({\partial}r)^{2}+\frac{(\xi_{1}-\xi_{2})r}{\sqrt{6}(\xi_{1}+\xi_{2}r^{2})^{2}}({\partial}\varphi)({\partial}r)$
$\displaystyle\quad+\frac{1}{2}\left(\frac{e^{\sqrt{\frac{2}{3}}{\varphi}}(1+6\xi_{1}+(1+6\xi_{2})r^{2})-6(\xi_{1}+\xi_{2}r^{2})}{6(\xi_{1}+\xi_{2}r^{2})(e^{\sqrt{\frac{2}{3}}{\varphi}}-1)}\right)({\partial}\varphi)^{2}\,,$
(13) $\displaystyle{\cal V}(r,\varphi)$ $\displaystyle={\cal
V}(r)\left(1-e^{-\sqrt{\frac{2}{3}}{\varphi}}\right)^{2}\,,$ (14)

with

${\cal
V}(r)=\frac{\lambda_{41}-\lambda_{31}r+\lambda_{412}r^{2}-\lambda_{32}r^{3}+\lambda_{42}r^{4}}{4\left(\xi_{2}r^{2}+\xi_{1}\right)^{2}}.$
(15)

The kinetic terms can be reduced further by introducing a canonical scalar
fieldChoi:2020ara :

$(\frac{dr_{c}}{dr})^{2}=\frac{(\xi^{2}_{1}+\xi^{2}_{2}r^{2})}{(\xi_{1}+\xi_{2}r^{2})^{3}},$
(16)

Then,

$\displaystyle\mathcal{L}_{{\rm kin}}$
$\displaystyle\simeq\frac{1}{2}({\partial}\varphi)^{2}+\frac{1}{2}\frac{(\xi_{1}^{2}+\xi_{2}^{2}r^{2})}{(\xi_{1}+\xi_{2}r^{2})^{3}}({\partial}r)^{2}+\frac{(\xi_{1}-\xi_{2})r}{\sqrt{6}(\xi_{1}+\xi_{2}r^{2})^{2}}({\partial}\varphi)({\partial}r)$
$\displaystyle\simeq\frac{1}{2}({\partial}\varphi)^{2}+\frac{1}{2}({\partial}r_{c})^{2}+({\partial}\varphi)({\partial}r_{c})\frac{(\xi_{1}-\xi_{2})r}{\sqrt{6}\sqrt{\xi_{1}^{2}+\xi_{2}^{2}r^{2}}\sqrt{\xi_{1}+\xi_{2}r^{2}}}\,.$
(17)

The kinetic mixing term $({\partial}\varphi)({\partial}r)$ can be canceled by
setting $\xi_{1}\approx\xi_{2}$ or suppressed by
$1/\sqrt{\xi_{1}+\xi_{2}r^{2}}$ with large non-minimal coupling. Then, the
effective action on inflation can be obtained in the form of vanishing non-
minimal couplings and kinetic mixing terms :

$\displaystyle\frac{{\cal
L}_{E}}{\sqrt{-g_{E}}}=-\frac{1}{2}{\mathcal{R}}-\frac{1}{2}({\partial}\varphi)^{2}-\frac{1}{2}({\partial}r_{c})^{2}-{\cal
V}(r_{c},\varphi)\,.$ (18)

Note that the potential in Eq.(15) has a symmetric exchange between $\Phi_{1}$
and $\Phi_{2}$, resulting in $r\rightarrow 1/r$ of nonzero values
$\lambda_{31}$ and $\lambda_{32}$, and also obtaining unequal values of non-
minimal couplings between $\xi_{1}$ and $\xi_{2}$. In the inflation epoch,
this feature can provide two distinguishable regions of $\xi_{2}/\xi_{1}$
value, i.e., $\xi_{2}\approx\xi_{1}$ and $\xi_{2}<\xi_{1}$.

Now, we discuss the stability of single-field inflation affiliated with the
potential minima(15) at a mixed point $r_{m}$. The condition for ${\cal
V}(r_{m})$ acting as an extreme minimum, requiring that:

$\displaystyle{\cal V}(r)\bigg{|}_{r=r_{m}}\geq 0\,,\quad\frac{\partial{\cal
V}(r)}{\partial r_{c}}\bigg{|}_{r=r_{m}}=0\,,\quad\frac{\partial^{2}{\cal
V}(r)}{\partial r_{c}^{2}}\bigg{|}_{r=r_{m}}\geq 0\,.$ (19)

Combining eq.(15) and eq.(16), the first and second derivatives of the
inflation potential are given respectively by:

$\displaystyle\frac{\partial{{\cal V}(r)}}{\partial r_{c}}\approx$
$\displaystyle\frac{M^{3}_{P}}{4(\xi_{1}+\xi_{2}r^{2})^{3/2}\sqrt{\xi_{1}^{2}+\xi_{2}^{2}r^{2}}}[\xi_{2}\lambda_{32}r^{4}+r^{3}(4\xi_{1}\lambda_{42}-2\xi_{2}\lambda_{412})+r^{2}(3\xi_{2}\lambda_{31}-3\xi_{1}\lambda_{32})$
$\displaystyle+r(2\xi_{1}\lambda_{412}-4\xi_{2}\lambda_{41})-\xi_{1}\lambda_{31}].$
(20)

and

$\displaystyle\frac{\partial^{2}{{\cal V}(r)}}{\partial r^{2}_{c}}\approx$
$\displaystyle\frac{M^{2}_{P}}{4(\xi_{1}+\xi_{2}r^{2})(\xi_{1}^{2}+\xi_{2}^{2}r^{2})^{2}}[r^{6}\left(2\xi_{2}^{4}\lambda_{412}-4\xi_{1}\xi_{2}^{3}\lambda_{42}\right)+r^{5}\left(\xi_{1}^{2}\xi_{2}^{2}\lambda_{32}+9\xi_{1}\xi_{2}^{3}\lambda_{32}-6\xi_{2}^{4}\lambda_{31}\right)$
$\displaystyle+r^{4}\left(8\xi_{1}^{2}\xi_{2}^{2}\lambda_{42}-10\xi_{1}\xi_{2}^{3}\lambda_{412}+12\xi_{2}^{4}\lambda_{32}\right)+r^{2}\left(12\xi_{1}^{4}\lambda_{42}-10\xi_{1}^{3}\xi_{2}\lambda_{412}+8\xi_{1}^{2}\xi_{2}^{2}\lambda_{32}\right)$
$\displaystyle+r\left(-6\xi_{1}^{4}\lambda_{32}+9\xi_{1}^{3}\xi_{2}\lambda_{31}+\xi_{1}^{2}\xi_{2}^{2}\lambda_{31}\right)+\left(2\xi_{1}^{4}\lambda_{412}-4\xi_{1}^{3}\xi_{2}\lambda_{41}\right)$
$\displaystyle+r^{3}\left(7\xi_{1}^{3}\xi_{2}\lambda_{32}-3\xi_{1}^{2}\xi_{2}^{2}\lambda_{31}-3\xi_{1}^{2}\xi_{2}^{2}\lambda_{32}+7\xi_{1}\xi_{2}^{3}\lambda_{31}\right)].$
(21)

##  3 Inflation driven by the two-component dark matter fields

We focus on inflationary dynamics driven by one dark matter field inflation or
two dark matter field. For the inflation driven by one dark matter field
$\Phi_{1}$ or $\Phi_{2}$, the eq.(2) indicates it can not work in
$r\rightarrow 0$ and $r\rightarrow\infty$ scenarios only if
$\xi_{1}\lambda_{31}=0$ and $\xi_{2}\lambda_{32}=0$ respectively. The mixed
inflation makes difference under a finite $r$ scenarios where $r$ makes the
potential minimum. The finite $r$ can be obtained by solving the quartic
equation in eq.(2). We will discuss each $r$ scenario in detail later.

###  3.1 Singlet dark matter inflation

The $r\rightarrow 0$ and $r\rightarrow\infty$ scenarios require
$\xi_{1}\lambda_{31}=0$ and $\xi_{2}\lambda_{32}=0$ respectively according to
eq.(2). For reducing parameters space, we only consider two case: A
$(\xi_{1}\neq 0,\lambda_{31}=0)$ with $r\rightarrow 0$ and B $(\xi_{2}\neq
0,\lambda_{32}=0)$ with $r\rightarrow\infty$. The minimal potential condition
of $r\rightarrow 0$ or $r\rightarrow\infty$ scenario can be found in the
second derivative potential eq.(2) and demands that:

$\xi_{1}^{4}\lambda_{412}-2\xi_{1}^{3}\xi_{2}\lambda_{41}>0$ (22)

or

$\xi_{2}^{4}\lambda_{412}-2\xi_{2}^{3}\xi_{1}\lambda_{42}>0$ (23)

###  3.2 Mixed dark matter inflation

The mixed dark matter inflation make difference when a finite $r$ makes the
inflation potential minimum and can be obtained by solving the quartic
equation in eq.(2) with the nonzero value of $\xi_{2}$ and $\lambda_{32}$.
Generally speaking, the discriminant and solutions expressions of the quartic
equation are too long and we do not present them here. For simplicity, we will
add on the conditions as follows:

$\displaystyle 4\xi_{1}\lambda_{42}-2\xi_{2}\lambda_{412}$
$\displaystyle=0\,,$ (24) $\displaystyle
3\xi_{2}\lambda_{31}-3\xi_{1}\lambda_{32}$ $\displaystyle=0\,,$ (25)
$\displaystyle 2\xi_{1}\lambda_{412}-4\xi_{2}\lambda_{41}$
$\displaystyle=0\,.$ (26)

Now the solutions of the quartic equation can be expressed as
$r_{0}=(\xi_{1}\lambda_{31}/{\xi_{2}\lambda_{32}})^{1/4}$ and we label this
scenario as case C.

We unify three cases in one frame nominally, at the inflation direction, the
Jordan-frame action can be written as follows:

$\displaystyle S$ $\displaystyle=\int
d^{4}x\,\sqrt{-g}\,\bigg{\\{}-\frac{M_{{\rm
P}}^{2}}{2}\left[1+\left(\xi_{1}+\xi_{2}r_{m}^{2}\right)\frac{\phi_{1}^{2}}{M_{{\rm
P}}^{2}}\right]\mathcal{R}$
$\displaystyle\quad-\frac{1+r_{m}^{2}}{2}g^{\mu\nu}\partial_{\mu}\phi_{1}\partial_{\nu}\phi_{1}-\frac{1}{4}(\lambda_{41}-\lambda_{31}r_{m}-\lambda_{32}r_{m}^{3}+\lambda_{42}r_{m}^{4})\phi_{1}^{4}\bigg{\\}}\,.$
(27)

After redefining $\Phi\equiv\sqrt{1+r_{m}^{2}}\,\phi_{1}$, the above action
becomes

$\displaystyle S=\int d^{4}x\,\sqrt{-g}\,\left\\{-\frac{M_{{\rm
P}}^{2}}{2}\left(1+\xi_{\Phi}\frac{\Phi^{2}}{M_{{\rm
P}}^{2}}\right)\mathcal{R}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\Phi\partial_{\nu}\Phi-\frac{1}{4}\lambda_{\Phi}\Phi^{4}\right\\}\,,$
(28)

where

$\displaystyle\xi_{\Phi}\equiv\frac{(\xi_{1}+\xi_{2}r_{m}^{2})}{1+r_{m}^{2}}\,,\qquad\lambda_{\Phi}\equiv\frac{(\lambda_{41}-\lambda_{31}r_{m}-\lambda_{32}r_{m}^{3}+\lambda_{42}r_{m}^{4})}{(1+r_{m}^{2})^{2}}\,.$
(29)

We end the part with a table[1]

| $\xi_{\Phi}$ | $\lambda_{\Phi}$ | ${\cal V}(r)$
---|---|---|---
Case A | $\xi_{1}$ | $\lambda_{41}$ | $\frac{M^{4}_{P}\lambda_{41}}{4\xi^{2}_{1}}$
Case B | $\xi_{2}$ | $\lambda_{42}$ | $\frac{M^{4}_{P}\lambda_{42}}{4\xi^{2}_{2}}$
Case C | $\frac{(\xi_{1}+\xi_{2}r_{0}^{2})}{1+r_{0}^{2}}$ | $\frac{(\lambda_{41}-\lambda_{31}r_{0}+\lambda_{412}r_{0}^{2}-\lambda_{32}r_{0}^{3}+\lambda_{42}r_{0}^{4})}{(1+r_{0}^{2})^{2}}$ | $\frac{M^{4}_{P}(\lambda_{41}-\lambda_{31}r_{0}+\lambda_{412}r_{0}^{2}-\lambda_{32}r_{0}^{3}+\lambda_{42}r_{0}^{4})}{4\left(\xi_{2}r_{0}^{2}+\xi_{1}\right)^{2}}$
Table 1: The effective parameters $\xi_{\Phi}$, $\lambda_{\Phi}$ and vacuum
energy $\mathcal{V}(r)$ in different cases. We do not expand terms with
eq.(24)-eq.(25) and $r_{0}$ just for concision

###  3.3 Inflationary observables

The potential given in eq.(14) belongs to the so-called Starobinsky
potentialsStarobinsky:1979ty ; Calmet:2016fsr , as we can see ${\cal
V}(r,\varphi)/{\cal V}(r)$ exhibits flat nature at high field values ensuring
slow roll when $r$ is determined by cases. The slow-roll parameters are given
by Baumann:2009ds :

$\displaystyle\epsilon=\frac{M_{{\rm P}}^{2}}{2}\left(\frac{{\cal
V}(r,\varphi)^{\prime}}{{\cal V}(r,\varphi)}\right)^{2}\,,\quad\eta=M_{{\rm
P}}^{2}\frac{{\cal V}(r,\varphi)^{\prime\prime}}{{\cal V}(r,\varphi)}\quad.$
(30)

where the prime denotes the derivative of the canonically normalized field
$\varphi$, which is defined by Bezrukov:2007ep :

$\displaystyle\frac{d\varphi}{d\Phi}=\frac{\sqrt{1+(1+6\xi_{\Phi})\xi_{\Phi}\Phi^{2}/{M^{2}_{P}}}}{1+\xi_{\Phi}\Phi^{2}/{M^{2}_{P}}}.$
(31)

With the help of slow-roll parameters, we can obtain the spectral index
$n_{s}$, the tensor-to-scalar ratio $r$ and the magnitude of the scalar power
spectrum $A_{s}$ at the horizon exit:

$\displaystyle n_{s}\approx 1-6\epsilon+2\eta,r\approx
16\epsilon,A_{s}=\frac{{\cal V}(r,\varphi)}{24\pi^{4}M^{2}_{P}\epsilon}.$ (32)

The inflation ends at $\epsilon\approx 1$, and the size of the universe had
expanded by $e$ times. The number of $e$-folds $N$ can be obtained as follows:

$\displaystyle
N=\frac{1}{M^{2}_{P}}\int_{\varphi_{e}}^{\varphi_{i}}d\varphi\frac{{\cal
V}(r,\varphi)}{{\cal V}^{\prime}(r,\varphi)}.$ (33)

where $\varphi_{i}$ is the value of $\varphi$ at the beginning of inflation
and $\varphi_{e}$ is the value of $\varphi$ at the end of the inflation. The
$\varphi_{e}$ can be obtained by setting slow roll inflation to the end when
$\epsilon\approx 1$ and we can compute $\varphi_{i}$ if $N$ is known. For
example, we take $N=60$, $\varphi_{e}$ can be estimated as
$M_{P}(4/3)^{1/4}/\sqrt{\xi_{\Phi}}$. Furthermore, we express $\varphi_{i}$ in
terms of $N$ as

$\displaystyle\varphi_{i}(N)\approx\sqrt{\frac{4N}{3\xi_{\Phi}}}M_{P}$ (34)

Then, in the large-field limit, $\epsilon$, $\eta$, $n_{s}$, $r$ and $A_{s}$
can be expressed in terms of the number of $e$-foldings,

$\displaystyle\epsilon$
$\displaystyle\approx\frac{3}{4N^{2}}\,,\quad\eta\approx-\frac{1}{N}\quad.$
(35) $\displaystyle n_{s}$ $\displaystyle\approx
1-\frac{2}{N}-\frac{9}{2N^{2}}\,,\quad r\approx\frac{12}{N^{2}}\,,\quad
A_{s}=(\frac{N^{2}}{72{\pi}^{2}})\frac{\lambda_{\Phi}}{\xi_{\Phi}^{2}}.$ (36)

Recent experiments have revealed the following limitations:

  * •

CMB normalization
The scalar power spectrum at the Planck pivot scale Planck:2018jri  is given
by

$\displaystyle\ln(10^{10}A_{s})=3.044\pm 0.014\quad\text{(68\% C.L. Planck
TT,TE,EE+lowE+lensing)}\,.$ (37)

Then, requiring $A_{s}\approx 2.1\times 10^{-9}$ at $N=60$, we obtain the
relation between the effective quartic and non-minimal couplings as follows:

$\displaystyle\frac{\lambda_{\Phi}}{\xi_{\Phi}^{2}}\approx 4.15\times
10^{-10}\,.$ (38)

For instance, for $\lambda_{\Phi}\simeq 0.1$ ($10^{-7}$), the non-minimal
coupling should take $\xi_{\Phi}\simeq 15500$ ($15.5$).

  * •

Spectral index and tensor-to-scalar ratio
The latest Planck data Planck:2018jri  read

$\displaystyle n_{s}$ $\displaystyle=0.9659\pm 0.0041\,$
$\displaystyle\text{(68\% C.L. Planck TT,TE,EE+lowEB+lensing)}\,,$ (39)
$\displaystyle r$ $\displaystyle<0.11\,$ $\displaystyle\text{(95\% C.L. Planck
TT,TE,EE+lowEB+lensing)}\,.$ (40)

We exhibit inflationary observables changed by the number of $e$-foldings $N$
in Fig.(1) and talk about $A_{s}$ in section 4.4.

Figure 1: The inflationary observables $\eta$, $\epsilon$, $n_{s}$, $r$
changed by $N$ according to the lastest Planck dataPlanck:2018jri .

##  4 $Z_{5}$ models and low energy constraints

###  4.1 Unitarity

The unitarity conditions come from the tree-level scalar-scalar scattering
matrix which is dominated by the quartic contact interaction. The s-wave
scattering amplitudes should lie under the perturbative unitarity limit, given
the requirement the eigenvalues of the S-matrix $M$ must be less than the
unitarity bound given by

$\displaystyle|ReM|<\frac{1}{2}.$ (41)

We extend those scattering matrix by charge and $Z_{5}$ change conservation
states and write down dominated conditions as:

$\displaystyle\frac{1}{2}\left|(2\lambda_{41}+\lambda_{412})\pm\sqrt{18\lambda^{2}_{31}+4\lambda^{2}_{41}-4\lambda_{41}\lambda_{412}+\lambda^{2}_{412}}\right|<8\pi,$
(42) $\displaystyle\frac{1}{6}\left|a_{2}\pm
2\sqrt{a_{2}^{2}-3a_{1}}\right|<8\pi.$ (43)

with

$\displaystyle
a_{1}=4(16\lambda_{41}\lambda_{42}-\lambda^{2}_{412}+24\lambda_{41}\lambda_{H}+24\lambda_{42}\lambda_{H}-2\lambda_{S1}^{2}-2\lambda_{S2}^{2}),$
(45) $\displaystyle a_{2}=4(2\lambda_{41}+2\lambda_{42}+3\lambda_{H}).$ (46)

More details can be found in Appendix A.

###  4.2 Perturbativity

We need perturbativity on the new dark matter quartic couplings,

$|3\lambda_{31}|<4\pi,|4\lambda_{41}|<4\pi|,|\lambda_{412}|<4\pi,|3\lambda_{32}|<4\pi,|4\lambda_{42}|<4\pi.$
(47)

###  4.3 Vacuum stability

The vacuum stability requires the potential under the
restriction,i.e,$\mathcal{V}>0$ for large field values. We only consider the
quartic couplings for their overwhelming contribution to $\mathcal{V}$. In
this section, we focus on the vacuum stability induced by dark matter scalars
and ignore the mixing quartic couplings with the SM Higgs, this scenario can
work when the mixing quartic couplings with the SM Higgs take positive or
small values as compared to couplings in the dark matter fields. The potential
with vacuum stability can be written as:

$\displaystyle\mathcal{V}_{DM}=\frac{1}{4}\lambda_{41}\phi_{1}^{4}+\frac{1}{4}\lambda_{42}\phi_{2}^{4}+\frac{1}{4}\lambda_{412}\phi_{1}^{2}\phi_{2}^{2}+\frac{1}{4}\lambda_{31}\phi_{1}^{3}\phi_{2}cos(3\theta_{1}+\theta_{2})+\frac{1}{4}\lambda_{32}\phi_{1}\phi_{2}^{3}cos(\theta_{1}-3\theta_{2}).$
(48)

After minimizing the potential for $3\theta_{1}+\theta_{2}$ and
$\theta_{1}-3\theta_{2}$ with $\phi_{1}\neq 0$ and $\phi_{2}\neq 0$, the above
potential has the form,

$\displaystyle\mathcal{V}_{DM}=\frac{1}{4}\lambda_{41}\phi_{1}^{4}+\frac{1}{4}\lambda_{42}\phi_{2}^{4}+\frac{1}{4}\lambda_{412}\phi_{1}^{2}\phi_{2}^{2}-\frac{1}{4}\lambda_{31}\phi_{1}^{3}\phi_{2}-\frac{1}{4}\lambda_{32}\phi_{1}\phi_{2}^{3}.$
(49)

Taking $X\equiv\frac{\phi_{1}}{\phi_{2}}$, the vacuum stability conditions
become

$\displaystyle\lambda_{41}>0,\lambda_{42}>0$ (50)

and

$\displaystyle f(X_{min})>0$ (51)

with

$\displaystyle
f(X_{min})=\frac{1}{4}\lambda_{41}X^{4}-\frac{1}{4}\lambda_{31}X^{3}+\frac{1}{4}\lambda_{412}X^{2}-\frac{1}{4}\lambda_{32}X+\frac{1}{4}\lambda_{42},$
(52)

where $X_{min}$ labels the global minimum when $f^{{}^{\prime}}(X_{min})=0$.
Then, solving $f^{{}^{\prime}}(X_{min})=0$ by the condition
$f^{{}^{\prime\prime}}(X_{min})>0$, we obtain the third condition (51) as

$\displaystyle\lambda_{41}X_{min}^{3}-\frac{3}{4}\lambda_{31}X^{2}+\frac{1}{2}\lambda_{412}X-\frac{1}{4}\lambda_{32}=0,$
(53)
$\displaystyle\lambda_{42}>\frac{1}{4}\lambda_{31}X^{3}_{min}+\frac{1}{2}\lambda_{412}X^{2}_{min}+\frac{3}{4}\lambda_{32}X_{min}.$
(54)

where

$\displaystyle
X_{min}=\Bigg{\\{}\begin{array}[]{ll}\left(P+\sqrt{P^{2}+Q^{3}}\right)^{1/3}+\left(P-\sqrt{P^{2}+Q^{3}}\right)^{1/3}\,,&D>0\\\
2\sqrt{-Q}\cos\left(\frac{1}{3}\cos^{-1}\left(\frac{P}{\sqrt{-Q^{3}}}\right)\right)\,,&D<0\end{array}$
(57)

with $D\equiv P^{2}+Q^{3}$,
$P\equiv\frac{4\lambda_{31}\lambda_{41}\lambda_{412}-8\lambda_{32}\lambda_{41}^{2}-\lambda_{31}^{3}}{64\lambda_{41}^{3}}$,
and
$Q\equiv\frac{8\lambda_{41}\lambda_{412}-3\lambda^{2}_{31}}{48\lambda_{41}^{2}}$.

###  4.4 Results

In this section, we discuss how low energy constraints and the scalar power
spectrum affect model parameters. As explained in the previous section, we
consider singlet dark matter and mixed dark matter inflation respectively. We
choose the following input parameters,

$\displaystyle\xi_{2}/\xi_{1},\xi_{1},\lambda_{31},\lambda_{32},\lambda_{412},\lambda_{41},\lambda_{42}.$
(58)

and fix $N$=62 in the analysis of scalar power spectrum.

####  4.4.1 case A and B

In the case of $r\rightarrow 0$ or $r\rightarrow\infty$ scenario, singlet dark
matter inflation takes effect under the condition
$\xi_{1}\lambda_{412}>2\xi_{2}\lambda_{41}$ or
$\xi_{2}\lambda_{412}>2\xi_{1}\lambda_{42}$. According to the symmetry of
$\xi_{1}\leftrightarrow\xi_{2}$ and $\lambda_{41}\leftrightarrow\lambda_{42}$,
Case B is equivalent to Case A, so only one of them needs to be studied and we
select the latter. To simplify parameter space further, we choose
$\lambda_{412}=3\lambda_{41}$ and $\xi_{1}\leq\xi_{2}$ so that
$\xi_{1}\lambda_{412}>2\xi_{2}\lambda_{41}$ always holds. Then, there only
three free parameters exists:$\lambda_{32},\lambda_{41},\xi_{1}$. We scan
parameter space at $\lambda_{32},\lambda_{41}\in[0.1,10]$ and
$\xi_{1}\in[10,10^{5}]$. Fig.(2) show the allowed area restricted by vacuum
stability, perturbativity, unitarity and scalar power spectrum. As we can see,
perturbativity and scalar power spectrum put more limits on $\lambda_{32}$ and
$\lambda_{41}$.

Figure 2: Parameter space for singlet dark matter inflation with
$\lambda_{412}=3\lambda_{41}$ at low-energy scale. Constraints from the model
and scalar power spectrum. The blue region represents the vacuum stability,
the dark blue and red region are favored due to perturbativity and unitarity
respectively.

####  4.4.2 case C

The mixed dark matter inflation takes place in the condition of
$\frac{\partial{{\cal V}(r)}}{\partial r_{c}}=0$ and $\frac{\partial^{2}{{\cal
V}(r)}}{\partial r^{2}_{c}}>0$ with a nonzero value of $\lambda_{31}$ and
$\lambda_{32}$. The eq.(24)-eq.(26) can be replaced by,

$\displaystyle
Y\equiv\xi_{2}/\xi_{1},\lambda_{42}=Y^{2}\lambda_{41},\lambda_{32}=Y\lambda_{31},\lambda_{412}=2Y\lambda_{41}.$
(59)

We take ${\lambda_{31},\lambda_{41},\xi_{1}}$ as free parameters and fix $Y$
with 0.1, 0.5, 1.0. The parameter space is sampled within
$\lambda_{31},\lambda_{41}\in[0.1,10]$ and $\xi_{1}\in[10,10^{5}]$. Fig.(3)
shows the viable range, bounded by vacuum stability, perturbativity, unitarity
and scalar power spectrum, and the scalar power spectrum gives a stronger
restriction on $\lambda_{31}$ and $\lambda_{41}$.

Figure 3: parameter space for mixed dark matter inflation with
$\lambda_{42}=Y^{2}\lambda_{41}$, $\lambda_{32}=Y\lambda_{31}$,
$\lambda_{412}=2Y\lambda_{41}$ at $Y=\xi_{2}/\xi_{1}=1.0,0.5,0.1$. The
parameters shown here are set at low-energy scale. The blue region represents
the vacuum stability. The dark blue and red region are favored due to
perturbativity and unitarity respectively.The histogram show the value of
$\xi$.

##  5 Dark matter

In this section, we discuss the dark matter phenomenology in the model. We
focus on the case of mixed dark matter inflation and fix $Y=1.0,0.5$ and $0.1$
respectively. In the $Z_{5}$ dark matter model, there are three $2\to 2$
processes related to dark matter, including annihilations, semi-annihilations
and conversions. The Boltzmann equations related to the two-component dark
matter are given by,

$\displaystyle\frac{dN_{1}}{dt}+3HN_{1}$ $\displaystyle=$
$\displaystyle-\sigma_{v}^{1100}\left(N_{1}^{2}-\bar{N}_{1}^{2}\right)-\sigma_{v}^{1120}\left(N_{1}^{2}-N_{2}\frac{\bar{N}_{1}^{2}}{\bar{N}_{2}}\right)-\sigma_{v}^{1122}\left(N_{1}^{2}-N_{2}^{2}\frac{\bar{N}_{1}^{2}}{\bar{N}_{2}^{2}}\right)$
(60)
$\displaystyle-\frac{1}{2}\sigma_{v}^{1112}\left(N_{1}^{2}-N_{1}N_{2}\frac{\bar{N}_{1}}{\bar{N}_{2}}\right)-\frac{1}{2}\sigma_{v}^{1222}\left(N_{1}N_{2}-N_{2}^{2}\frac{\bar{N}_{1}}{\bar{N}_{2}}\right)$
$\displaystyle-\frac{1}{2}\sigma_{v}^{1220}\left(N_{1}N_{2}-N_{2}\bar{N}_{1}\right)+\frac{1}{2}\sigma_{v}^{2210}(N_{2}^{2}-N_{1}\frac{\bar{N}_{2}^{2}}{\bar{N}_{1}})\,,$
$\displaystyle\frac{dN_{2}}{dt}+3HN_{2}$ $\displaystyle=$
$\displaystyle-\sigma_{v}^{2200}\left(N_{2}^{2}-\bar{N}_{2}^{2}\right)-\sigma_{v}^{2210}\left(N_{2}^{2}-N_{1}\frac{\bar{N}_{2}^{2}}{\bar{N}_{1}}\right)-\sigma_{v}^{2211}\left(N_{2}^{2}-N_{1}^{2}\frac{\bar{N}_{2}^{2}}{\bar{N}_{1}^{2}}\right)$
(61)
$\displaystyle-\frac{1}{2}\sigma_{v}^{2221}\left(N_{2}^{2}-N_{1}N_{2}\frac{\bar{N}_{2}}{\bar{N}_{1}}\right)-\frac{1}{2}\sigma_{v}^{1211}\left(N_{1}N_{2}-N_{1}^{2}\frac{\bar{N}_{2}}{\bar{N}_{1}}\right)$
$\displaystyle-\frac{1}{2}\sigma_{v}^{1210}\left(N_{1}N_{2}-N_{1}\bar{N}_{2}\right)+\frac{1}{2}\sigma_{v}^{1120}(N_{1}^{2}-N_{2}\frac{\bar{N}_{1}^{2}}{\bar{N}_{2}}).$

Where $N_{i}$ ($i=1,2$) are the number densities of $\phi_{i}$, and
$\bar{N}_{i}$ the equilibrium values. $\sigma_{v}^{abcd}$ corresponds to the
thermally averaged cross section of the $a+b\to c+d$ process, and the digits
$0,1,2$ above the $\sigma_{v}$ represent different particles in the model,
where $0$ is used for SM particles, $1$ represents $\phi_{1}$ or
$\phi_{1}^{\dagger}$ and $2$ for $\phi_{2}$ or $\phi_{2}^{\dagger}$. To
estimate the contribution of the three $2\to 2$ processes to the relic density
of $\phi_{1}$, one can define the following three parameters for annihilation,
semi-annihilation and conversion respectively:

$\displaystyle\chi_{an}$
$\displaystyle\equiv\frac{\sigma_{v}^{1100}}{\overline{\sigma_{v}^{1}}},\quad\chi_{se}\equiv\frac{\frac{1}{2}(\sigma_{v}^{1120}+\sigma_{v}^{1220}+\sigma_{v}^{1022})}{\overline{\sigma_{v}^{1}}},\quad\chi_{co}\equiv\frac{\sigma_{v}^{1122}+\sigma_{v}^{1112}+\sigma_{v}^{1222}}{\overline{\sigma_{v}^{1}}},$
(62)

with

$\displaystyle\overline{\sigma_{v}^{1}}\equiv\sigma_{v}^{1100}+\frac{1}{2}\sigma_{v}^{1120}+\sigma_{v}^{1122}+\sigma_{v}^{1112}+\sigma_{v}^{1222}+\frac{1}{2}\sigma_{v}^{1220}+\frac{1}{2}\sigma_{v}^{1022}.$
(63)

We calculate the DM relic density with micrOMEGAs numerically Belanger:2014vza
, where we scan the parameter space of

$\displaystyle M_{1}\subset[1500\ \mathrm{GeV},1800\
\mathrm{GeV}],M_{2}\subset[1500\rm\ GeV,3600\rm\ GeV]$ (64)

with $M_{1}$ the $\phi_{1}$ mass and $M_{2}$ the $\phi_{2}$ mass. We limit the
parameter space with relic density constraint. On the other hand, we assume
$M_{1}<M_{2}$ and $M_{2}<2M_{1}$ to guarantee both $\phi_{1}$ and $\phi_{2}$
the dark matter candidate. For different $Y$, we fix
$\lambda_{S1}=\lambda_{S2}=0.75$ and set other parameters constants
respectively, the results are given as follows.

Figure 4: Parameter space of $M_{1}-\Omega_{1}$ satisfying DM relic density
constraint in the case of $Y=0.1$, and points with different colors correspond
to the fraction of $\phi_{1}$ to the total relic density.

Figure 5: Parameter space of $M_{1}-M_{2}$ satisfying DM relic density
constraint in the case of $Y=0.1$, and points with different colors describe
the fraction of conversion in the $2\to 2$ processes.

Figure 6: Parameter space of $M_{1}-\Omega_{1}$ satisfying DM relic density
constraint in the case of $Y=0.5$, and points with different colors correspond
to the fraction of $\phi_{1}$ to the total relic density.

Figure 7: Parameter space of $M_{1}-M_{2}$ satisfying DM relic density
constraint in the case of $Y=0.5$, and points with different colors describe
the fraction of conversion in the $2\to 2$ processes.

Figure 8: Parameter space of $M_{1}-\Omega_{1}$ satisfying DM relic density
constraint in the case of $Y=1$, and points with different colors correspond
to the fraction of $\phi_{1}$ to the total relic density.

Figure 9: Parameter space of $M_{1}-M_{2}$ satisfying DM relic density
constraint in the case of $Y=1$, and points with different colors describe the
fraction of conversion in the $2\to 2$ processes.

In Fig. 5 and Fig. 5, we fix $Y=0.1$, $\lambda_{41}=1$,
$\lambda_{42}=Y^{2}\lambda_{41}$, $\lambda_{412}=0.002$,
$\lambda_{31}=0.073030$ and $\lambda_{32}=0.000024$. We give the relationship
between $M_{1}$ and $\Omega_{1}$ satisfying relic density constraint in Fig.
5, where points with different colors represent the fraction of $\phi_{1}$ in
the total relic density and $\Omega_{1}$ is the relic density of $\phi_{1}$
while $\Omega$ is the total relic density. The relic density of $\phi_{1}$ is
stringently limited within a narrow region of $(0.09,0.13)$, increasing with
the increase of $M_{1}$. In addition, the fraction of $\phi_{1}$ with
$\Omega_{1}/\Omega$ is also increasing with the increase of $M_{1}$, and for
$M_{1}>1650$ GeV, the dark matter relic density is almost composed of
$\phi_{1}$. According to Fig. 5, the viable region satisfying relic density
constraint for $M_{2}$ is [2100, 3600]GeV, and we can have a more flexible
parameter space for $M_{2}$ when $M_{1}>1650$ GeV since $\phi_{2}$ will be
less relevant in the dark matter relic density. In addition, Fig. 5 also gives
the fraction of the conversion between $\phi_{1}$ and $\phi_{2}$ among the
$2\to 2$ processes. The value of $\chi_{co}$ is at $10^{-9}$ level and can be
negligible, which means the relic density is mainly determined by the
annihilations as well as semi-annihilations.

In Fig. 7 and Fig. 7, we discuss the relationship between $M_{1}-\Omega_{1}$
and $M_{1}-M_{2}$ in the case of $Y=0.5$, where we fix other parameters with
$\lambda_{41}=1$, $\lambda_{42}=Y^{2}\lambda_{41}=0.25$, $\lambda_{412}=0.25$,
$\lambda_{31}=0.816497$ and $\lambda_{32}=0.034021$. We have a more flexible
value for $\Omega_{1}$ and $M_{2}$ compared with the case of $Y=0.1$, and the
viable value for $M_{2}$ satisfying relic density constraint can decrease to
1600 GeV. Nevertheless, the contribution of $\phi_{1}$ to the relic density is
still dominant over $\phi_{2}$, and dark matter relic density is almost given
by $\phi_{1}$ when $M_{1}$ is larger than 1700 GeV.

According to Fig. 9 and Fig. 9, we give the results of $Y=1$, where we set
$\lambda_{42}=\lambda_{41}=1.01$, $\lambda_{412}=2.02$,
$\lambda_{31}=2.332495$ and $\lambda_{32}=0.777498$. Similarly, $\phi_{1}$ is
dominant on the dark matter relic density in Fig. 9, although the conversion
of $\phi_{1}$ to $\phi_{2}$ increases up to 0.56 as in Fig. 9. Such a result
is consistent with the conclusion of ref. Belanger:2020hyh  that the dark
matter density is always dominated by the lighter component in the $Z_{5}$
model.

The elastic scattering of the dark matter off nuclei in the model arises from
the Higgs portal interaction $\lambda_{Si}$ with $i=1,2$, which can put a
stringent constraint on the parameter space. The expression of the spin-
independent (SI) cross section can be given by Belanger:2020hyh .

$\displaystyle\sigma_{i}=\frac{\lambda_{Si}^{2}}{4\pi}\frac{\mu_{R}^{2}m_{p}^{2}f_{p}^{2}}{m_{H}^{4}M_{i}^{2}}$
(65)

where $i=1,2$, $\mu_{R}$ is the reduced mass, $m_{p}$ is the proton mass,
$m_{H}$ the SM Higgs mass and $f_{p}\approx 0.3$ is the quark content of the
proton. Experiments on the direct detection of dark matter can be found in
XENON:2018voc ; PandaX-4T:2021bab , and the PandaX-4T experiments
PandaX-4T:2021bab  put the most stringent constraint on the spin-independent
dark matter.

##  6 Conclusion

In this work, we consider the dark matter and inflation in the $Z_{5}$ model
of the two-component dark matter. Two singlet scalars carrying different
$Z_{5}$ charges are introduced to the SM as the dark matter candidate. Either
one of the dark matter particles or both the two scalars can play the role of
the inflaton, which correspond to the case of singlet dark matter inflation
and mixed dark matter inflation.

We estimate the parameter space by considering all the theoretical constraints
in the case of singlet dark matter inflation and mixed dark matter inflation
respectively. We focus on the latter inflation case and then consider the dark
matter relic density constraint. In the case of successful inflation, the
viable parameter space not only satisfies correct relic abundance but is also
limited by the direct detection result. It turns out that the dark matter
density is dominated by the light component, which is consistent with the
feature of the $Z_{5}$ model of the two-component dark matter.

## Acknowledgements

Part of the work is done by the author Aigeng Yang who died of illness and
left us far away. Hao Sun is supported by the National Natural Science
Foundation of China (Grant No. 12075043, No.12147205).

##  Appendix A Unitarity constraints

With the help of mutually unmixed sets of channels with definite charge, the
full set of two-body scalar scalar scattering processes leads to $26\times 26$
matrix which can be decomposed into four block submatrices. There are three
submatrices with zero-charge channels in the initial/final states:
$S^{(1)}(10\times 10)$, $S^{(2)}(5\times 5)$, $S^{(3)}(9\times 9)$. The last
one is $S^{(4)}(6\times 6)$ with one-change channels.

The first submatrix ${\cal M}_{1}$ corresponds to the scattering whose initial
and final states are one of the following:
$(G^{0}\phi_{1},G^{0}\phi_{2},G^{0}\phi^{*}_{1},G^{0}\phi^{*}_{2},h\phi_{1},h\phi_{2},h\phi^{*}_{1},h\phi^{*}_{2},\phi_{1}\phi^{*}_{2},\phi^{*}_{1}\phi_{2})$.
From the vertex factors one can find

${\cal
M}_{1}=\left(\begin{array}[]{cccccccccc}\frac{1}{2}\lambda_{S1}&0&0&0&0&0&0&0&0&0\\\
0&\frac{1}{2}\lambda_{S2}&0&0&0&0&0&0&0&0\\\
0&0&\frac{1}{2}\lambda_{S1}&0&0&0&0&0&0&0\\\
0&0&0&\frac{1}{2}\lambda_{S2}&0&0&0&0&0&0\\\
0&0&0&0&\frac{1}{2}\lambda_{S1}&0&0&0&0&0\\\
0&0&0&0&0&\frac{1}{2}\lambda_{S2}&0&0&0&0\\\
0&0&0&0&0&0&\frac{1}{2}\lambda_{S1}&0&0&0\\\
0&0&0&0&0&0&0&\frac{1}{2}\lambda_{S2}&0&0\\\
0&0&0&0&0&0&0&0&\lambda_{412}&0\\\
0&0&0&0&0&0&0&0&0&\lambda_{412}\end{array}\right),$ (66)

We find that ${\cal M}_{1}$ has the following ten eigenvalues $e^{i}_{j}$:

$\displaystyle e^{4}_{1}=\frac{1}{2}\lambda_{S1},$ (67) $\displaystyle
e^{4}_{2}=\frac{1}{2}\lambda_{S2},$ (68) $\displaystyle
e^{2}_{3}=\lambda_{412}.$ (69)

The second submatrix ${\cal M}_{2}$ corresponds to the scattering whose
initial and final states are one of the following:
$(G^{+}G^{-},\frac{G^{0}G^{0}}{\sqrt{2}},\frac{hh}{\sqrt{2}},\phi_{1}\phi^{*}_{1},\phi_{2}\phi^{*}_{2})$,
where the $\sqrt{2}$ accounts for identical particles. Again, from the vertex
factors one can find

${\cal
M}_{2}=\left(\begin{array}[]{ccccc}4\lambda_{H}&\sqrt{2}\lambda_{H}&\sqrt{2}\lambda_{H}&\lambda_{S1}&\lambda_{S2}\\\
\sqrt{2}\lambda_{H}&3\lambda_{H}&\lambda_{H}&\frac{1}{\sqrt{2}}\lambda_{S1}&\frac{1}{\sqrt{2}}\lambda_{S2}\\\
\sqrt{2}\lambda_{H}&\lambda_{H}&3\lambda_{H}&\frac{1}{\sqrt{2}}\lambda_{S1}&\frac{1}{\sqrt{S2}}\\\
\lambda_{S1}&\frac{1}{\sqrt{2}}\lambda_{S1}&\frac{1}{\sqrt{2}}\lambda_{S1}&4\lambda_{41}&\lambda_{412}\\\
\lambda_{S2}&\frac{1}{\sqrt{2}}\lambda_{S2}&\frac{1}{\sqrt{2}}\lambda_{S2}&\lambda_{412}&4\lambda_{42}\end{array}\right).$
(70)

We find that ${\cal M}_{2}$ has the following five eigenvalues:

$\displaystyle e_{i}^{1}$ $\displaystyle=\frac{1}{2}X_{i},(i=1,2,3)$ (71)
$\displaystyle e_{4}^{2}$ $\displaystyle=2\lambda_{H}.$ (72)

Where $X_{i}$ denotes the roots of the cubic equation as follows:

$\displaystyle X^{3}-4(2\lambda_{41}+2\lambda_{42}+3\lambda_{H})X^{2}$
$\displaystyle+4(16\lambda_{41}\lambda_{42}-\lambda^{2}_{412}+24\lambda_{41}\lambda_{H}+24\lambda_{42}\lambda_{H}-2\lambda^{2}_{S1}-2\lambda^{2}_{S2})X$
$\displaystyle-16(48\lambda_{41}\lambda_{42}\lambda_{H}-3\lambda^{2}_{412}\lambda_{H}-4\lambda_{42}\lambda^{2}_{S1}+2\lambda_{412}\lambda_{S1}\lambda_{S2}-4\lambda_{41}\lambda^{2}_{S2})=0.$
(73)

We apply the Samuelson’s inequality Samine  to find the region of real roots
$X_{i}$, i.e. for the polynomial

$\displaystyle a_{n}x^{n}+a_{n-1}x^{n-1}+a_{1}x+a_{0}$ (74)

with only real roots, the roots lie in the interveral bounded by

$\displaystyle
x_{\pm}=-\frac{a_{n-1}}{na_{n}}\pm\frac{n-1}{na_{n}}\sqrt{a_{n-1}^{2}-\frac{2n}{n-1}{a_{n}}{a_{n-2}}}.$
(75)

The third submatrix ${\cal M}_{3}$ corresponds to the scattering whose initial
and dinal states are one of the following:
$(hG^{0},\frac{1}{\sqrt{2}}\phi_{1}\phi_{1},\frac{1}{\sqrt{2}}\phi^{*}_{1}\phi^{*}_{1},\frac{1}{\sqrt{2}}\phi_{2}\phi_{2},\frac{1}{\sqrt{2}}\phi^{*}_{2}\phi^{*}_{2},\phi_{1}\phi_{2},\phi^{*}_{1}\phi^{*}_{2},\phi_{1}\phi^{*}_{2},\phi^{*}_{1}\phi_{2})$.
One can find

${\cal M}_{3}=\left(\begin{array}[]{ccccccccc}2\lambda_{H}&0&0&0&0&0&0&0&0\\\
0&2\lambda_{41}&0&0&0&0&\frac{3}{\sqrt{2}}\lambda_{31}&0&0\\\
0&0&2\lambda_{41}&0&0&\frac{3}{\sqrt{2}}\lambda_{31}&0&0&0\\\
0&0&0&2\lambda_{42}&0&0&0&\frac{3}{\sqrt{2}}\lambda_{32}&0\\\
0&0&0&0&2\lambda_{42}&0&0&0&\frac{3}{\sqrt{2}}\lambda_{32}\\\
0&0&\frac{3}{\sqrt{2}}\lambda_{31}&0&0&\lambda_{412}&0&0&0\\\
0&\frac{3}{\sqrt{2}}\lambda_{31}&0&0&0&0&\lambda_{412}&0&0\\\
0&0&0&\frac{3}{\sqrt{2}}\lambda_{32}&0&0&0&\lambda_{412}&0\\\
0&0&0&0&\frac{3}{\sqrt{2}}\lambda_{32}&0&0&0&\lambda_{412}\end{array}\right).$
(76)

${\cal M}_{3}$ has the following nine eigenvalues:

$\displaystyle e_{1,2}^{4}$
$\displaystyle=\frac{1}{2}\left(2\lambda_{41}+\lambda_{412}\pm\sqrt{18\lambda_{31}^{2}+4\lambda^{2}_{41}-4\lambda_{41}\lambda_{412}+\lambda^{2}_{412}}\right),$
(77) $\displaystyle e_{3}$ $\displaystyle=2\lambda_{H}.$ (78)

The last submatrix $\mathcal{M}_{4}$ corresponds to the one-change channels
occur for two-by-two body scattering between the six charged states
$(hG^{+},G^{0}G^{+},\phi_{1}G^{+},\phi_{2}G^{+},\phi^{*}_{1}G^{+},\phi^{*}_{2}G^{+})$.
We have ${\cal M}_{4}$ and its eigenvalues as follows:

${\cal M}_{4}=\left(\begin{array}[]{cccccc}2\lambda_{H}&0&0&0&0&0\\\
0&2\lambda_{H}&0&0&0&0\\\ 0&0&\lambda_{S1}&0&0&0\\\ 0&0&0&\lambda_{S2}&0&0\\\
0&0&0&0&\lambda_{S1}&0\\\ 0&0&0&0&0&\lambda_{S2}\end{array}\right),$ (79)
$\displaystyle e^{2}_{1}$ $\displaystyle=2\lambda_{H},$ (80) $\displaystyle
e^{2}_{2}$ $\displaystyle=\lambda_{S1},$ (81) $\displaystyle e^{2}_{3}$
$\displaystyle=\lambda_{S2}.$ (82)

##  Appendix B Renormalization Group Equations

To properly connect inflation scale with low energy scale of the dark matter
physics, one should consider the RGEs of the relevant parameters. We obtain
the Renormalization Group Equations (RGEs) of the model with SARAH
Staub:2013tta  at one-loop level, and the beta functions of the relavant
quartic couplings are given by,

$\displaystyle\beta_{\lambda_{31}}$
$\displaystyle=6\lambda_{31}\Big{(}2\lambda_{41}+\lambda_{412}\Big{)}$ (83)
$\displaystyle\beta_{\lambda_{32}}$
$\displaystyle=6\lambda_{32}\Big{(}2\lambda_{42}+\lambda_{412}\Big{)}$ (84)
$\displaystyle\beta_{\lambda_{H}}$
$\displaystyle=+\frac{27}{200}g_{1}^{4}+\frac{9}{20}g_{1}^{2}g_{2}^{2}+\frac{9}{8}g_{2}^{4}-\frac{9}{5}g_{1}^{2}\lambda_{H}-9g_{2}^{2}\lambda_{H}+24\lambda_{H}^{2}+\lambda_{S1}^{2}+\lambda_{S2}^{2}$
$\displaystyle+12\lambda_{H}y_{t}^{2}-6y_{t}^{4}$ (85)
$\displaystyle\beta_{\lambda_{42}}$
$\displaystyle=20\lambda_{42}^{2}+2\lambda_{S2}^{2}+\frac{9}{2}|\lambda_{32}|^{2}+\lambda_{412}^{2}$
(86) $\displaystyle\beta_{\lambda_{412}}$
$\displaystyle=4\Big{(}2\lambda_{412}\lambda_{42}+2\lambda_{41}\lambda_{412}+\lambda_{S1}\lambda_{S2}+\lambda_{412}^{2}\Big{)}+9|\lambda_{31}|^{2}+9|\lambda_{32}|^{2}$
(87) $\displaystyle\beta_{\lambda_{41}}$
$\displaystyle=20\lambda_{41}^{2}+2\lambda_{S1}^{2}+\frac{9}{2}|\lambda_{31}|^{2}+\lambda_{412}^{2}$
(88) $\displaystyle\beta_{\lambda_{H}}$
$\displaystyle=+\frac{27}{200}g_{1}^{4}+\frac{9}{20}g_{1}^{2}g_{2}^{2}+\frac{9}{8}g_{2}^{4}-\frac{9}{5}g_{1}^{2}\lambda_{H}-9g_{2}^{2}\lambda_{H}+24\lambda_{H}^{2}+\lambda_{S1}^{2}+\lambda_{S2}^{2}$
(89) $\displaystyle+12\lambda_{H}y_{t}^{2}-6y^{4}_{t}$ (90)
$\displaystyle\beta_{\lambda_{S2}}$
$\displaystyle=+2\lambda_{412}\lambda_{S1}-\frac{9}{10}g_{1}^{2}\lambda_{S2}-\frac{9}{2}g_{2}^{2}\lambda_{S2}+8\lambda_{42}\lambda_{S2}+12\lambda_{H}\lambda_{S2}+4\lambda_{S2}^{2}+6\lambda_{S2}y_{t}^{2}$
(91) $\displaystyle\beta_{\lambda_{S1}}$
$\displaystyle=-\frac{9}{10}g_{1}^{2}\lambda_{S1}-\frac{9}{2}g_{2}^{2}\lambda_{S1}+8\lambda_{41}\lambda_{S1}+12\lambda_{H}\lambda_{S1}+4\lambda_{S1}^{2}+2\lambda_{412}\lambda_{S2}+6\lambda_{S1}y_{t}^{2}$
(92)

where $y_{t}$ is the Top Yukawa coupling. For the trilinear scalar couplings,
the RGEs are given by,

$\displaystyle\beta_{\mu_{S2}}$
$\displaystyle=4\Big{(}\lambda_{412}+\lambda_{42}\Big{)}\mu_{S2}+6\lambda_{31}\mu_{S1}^{*}+6\mu_{S1}\lambda_{32}^{*}$
(93) $\displaystyle\beta_{\mu_{S1}}$
$\displaystyle=4\lambda_{412}\mu_{S1}+4\lambda_{41}\mu_{S1}+6\lambda_{31}\mu_{S2}^{*}+6\lambda_{32}\mu_{S2}.$
(94)

## References

  * (1) ATLAS collaboration, G. Aad et al., _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Phys. Lett. B_ 716 (2012) 1–29, [1207.7214].
  * (2) CMS collaboration, S. Chatrchyan et al., _Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC_ , _Phys. Lett. B_ 716 (2012) 30–61, [1207.7235].
  * (3) Y. Liu, _Higgs inflation and its extensions and the further refining dS swampland conjecture_ , _Eur. Phys. J. C_ 81 (2021) 1122, [2112.14571].
  * (4) WMAP collaboration, E. Komatsu et al., _Seven-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation_ , _Astrophys. J. Suppl._ 192 (2011) 18, [1001.4538].
  * (5) F. L. Bezrukov and M. Shaposhnikov, _The Standard Model Higgs boson as the inflaton_ , _Phys. Lett. B_ 659 (2008) 703–706, [0710.3755].
  * (6) F. Bezrukov, A. Magnin, M. Shaposhnikov and S. Sibiryakov, _Higgs inflation: consistency and generalisations_ , _JHEP_ 01 (2011) 016, [1008.5157].
  * (7) M. Sher, _Electroweak Higgs Potentials and Vacuum Stability_ , _Phys. Rept._ 179 (1989) 273–418.
  * (8) C. P. Burgess, H. M. Lee and M. Trott, _Power-counting and the Validity of the Classical Approximation During Inflation_ , _JHEP_ 09 (2009) 103, [0902.4465].
  * (9) R. N. Lerner and J. McDonald, _Higgs inflation and naturalness_ , _Journal of Cosmology and Astroparticle Physics_ 2010 (apr, 2010) 015–015.
  * (10) G. F. Giudice and H. M. Lee, _Unitarizing Higgs Inflation_ , _Phys. Lett. B_ 694 (2011) 294–300, [1010.1417].
  * (11) A. Linde, _Inflationary Cosmology after Planck 2013_ , in _100e Ecole d’Ete de Physique: Post-Planck Cosmology_ , pp. 231–316, 2015.  1402.0526.  DOI.
  * (12) L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080, [0911.1120].
  * (13) C. E. Yaguna, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ (2011) 60.
  * (14) H. Baer, K.-Y. Choi, J. E. Kim and L. Roszkowski, _Dark matter production in the early Universe: beyond the thermal WIMP paradigm_ , _Phys. Rept._ 555 (2015) 1–60, [1407.0017].
  * (15) Y. Hochberg, E. Kuflik, T. Volansky and J. G. Wacker, _Mechanism for Thermal Relic Dark Matter of Strongly Interacting Massive Particles_ , _Phys. Rev. Lett._ 113 (2014) 171301, [1402.5143].
  * (16) J. M. Cline, K. Kainulainen, P. Scott and C. Weniger, _Update on scalar singlet dark matter_ , _Phys. Rev. D_ 88 (2013) 055025, [1306.4710].
  * (17) F. S. Queiroz and K. Sinha, _The Poker Face of the Majoron Dark Matter Model: LUX to keV Line_ , _Phys. Lett. B_ 735 (2014) 69–74, [1404.1400].
  * (18) B. Batell, _Dark Discrete Gauge Symmetries_ , _Phys. Rev. D_ 83 (2011) 035006, [1007.0045].
  * (19) G. Bélanger, K. Kannike, A. Pukhov and M. Raidal, _Minimal semi-annihilating $\mathbb{Z}_{N}$ scalar dark matter_, _JCAP_ 06 (2014) 021, [1403.4960].
  * (20) C. Boehm, P. Fayet and J. Silk, _Light and heavy dark matter particles_ , _Phys. Rev. D_ 69 (2004) 101302, [hep-ph/0311143].
  * (21) E. Ma, _Supersymmetric Model of Radiative Seesaw Majorana Neutrino Masses_ , _Annales Fond. Broglie_ 31 (2006) 285, [hep-ph/0607142].
  * (22) Q.-H. Cao, E. Ma, J. Wudka and C. P. Yuan, _Multipartite dark matter_ , 0711.3881.
  * (23) T. Hur, H.-S. Lee and S. Nasri, _A Supersymmetric U(1)-prime model with multiple dark matters_ , _Phys. Rev. D_ 77 (2008) 015008, [0710.2653].
  * (24) H.-S. Lee, _Lightest U-parity Particle (LUP) dark matter_ , _Phys. Lett. B_ 663 (2008) 255–258, [0802.0506].
  * (25) K. M. Zurek, _Multi-Component Dark Matter_ , _Phys. Rev. D_ 79 (2009) 115002, [0811.4429].
  * (26) S. Profumo, K. Sigurdson and L. Ubaldi, _Can we discover multi-component WIMP dark matter?_ , _JCAP_ 12 (2009) 016, [0907.4374].
  * (27) H. Baer, A. Lessa, S. Rajagopalan and W. Sreethawong, _Mixed axion/neutralino cold dark matter in supersymmetric models_ , _JCAP_ 06 (2011) 031, [1103.5413].
  * (28) S. Esch, M. Klasen and C. E. Yaguna, _A minimal model for two-component dark matter_ , _JHEP_ 09 (2014) 108, [1406.0617].
  * (29) G. Bélanger, A. Pukhov, C. E. Yaguna and O. Zapata, _The Z 5 model of two-component dark matter_, _JHEP_ 09 (2020) 030, [2006.14922].
  * (30) F. Kahlhoefer and J. McDonald _Journal of Cosmology and Astroparticle Physics_ 2015 (nov, 2015) 015–015.
  * (31) R. N. Lerner and J. McDonald, _Gauge singlet scalar as inflaton and thermal relic dark matter_ , _Phys. Rev. D_ 80 (2009) 123507, [0909.0520].
  * (32) A. Aravind, M. Xiao and J.-H. Yu, _Higgs Portal to Inflation and Fermionic Dark Matter_ , _Phys. Rev. D_ 93 (2016) 123513, [1512.09126].
  * (33) N. Okada and Q. Shafi, _WIMP Dark Matter Inflation with Observable Gravity Waves_ , _Phys. Rev. D_ 84 (2011) 043533, [1007.1672].
  * (34) J.-O. Gong, H. M. Lee and S. K. Kang, _Inflation and dark matter in two Higgs doublet models_ , _JHEP_ 04 (2012) 128, [1202.0288].
  * (35) S. Choubey and A. Kumar, _Inflation and Dark Matter in the Inert Doublet Model_ , _JHEP_ 11 (2017) 080, [1707.06587].
  * (36) S.-M. Choi, J. Kim, H. M. Lee and B. Zhu, _Connecting between inflation and dark matter in models with gauged Z 3 symmetry_, _JHEP_ 06 (2020) 135, [2003.11823].
  * (37) D. I. Kaiser, _Conformal Transformations with Multiple Scalar Fields_ , _Phys. Rev. D_ 81 (2010) 084044, [1003.1159].
  * (38) A. A. Starobinsky, _Spectrum of relict gravitational radiation and the early state of the universe_ , _JETP Lett._ 30 (1979) 682–685.
  * (39) X. Calmet and I. Kuntz, _Higgs Starobinsky Inflation_ , _Eur. Phys. J. C_ 76 (2016) 289, [1605.02236].
  * (40) D. Baumann, _Inflation_ , in _Theoretical Advanced Study Institute in Elementary Particle Physics: Physics of the Large and the Small_ , pp. 523–686, 2011.  0907.5424.  DOI.
  * (41) Planck collaboration, Y. Akrami et al., _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10, [1807.06211].
  * (42) G. Bélanger, F. Boudjema, A. Pukhov and A. Semenov, _micrOMEGAs4.1: two dark matter candidates_ , _Comput. Phys. Commun._ 192 (2015) 322–329, [1407.6129].
  * (43) XENON collaboration, E. Aprile et al., _Dark Matter Search Results from a One Ton-Year Exposure of XENON1T_ , _Phys. Rev. Lett._ 121 (2018) 111302, [1805.12562].
  * (44) PandaX-4T collaboration, Y. Meng et al., _Dark Matter Search Results from the PandaX-4T Commissioning Run_ , _Phys. Rev. Lett._ 127 (2021) 261802, [2107.13438].
  * (45) P. A. Samuelson, _How deviant can you be?_ , _Journal of the American Statistical Association_ 63 (1968) 1522–1525.
  * (46) F. Staub, _SARAH 4 : A tool for (not only SUSY) model builders_ , _Comput. Phys. Commun._ 185 (2014) 1773–1790, [1309.7223].


# Surface criticality of antiferromagnetic Potts model

Li-Ru Zhang  School of Science and Engineering of Mathematics and Physics,
Anhui University of Technology, Maanshan, Anhui 243002, China  Chengxiang Ding
<EMAIL_ADDRESS>School of Science and Engineering of Mathematics and
Physics, Anhui University of Technology, Maanshan, Anhui 243002, China  Youjin
Deng  Department of Modern Physics, University of Science and Technology of
China, Hefei, 230027, China  MinJiang Collaborative Center for Theoretical
Physics, College of Physics and Electronic Information Engineering, Minjiang
University, Fuzhou 350108, China  Long Zhang<EMAIL_ADDRESS>Kavli
Institute for Theoretical Sciences and CAS Center for Excellence in
Topological Quantum Computation, University of Chinese Academy of Sciences,
Beijing 100190, China

###### Abstract

We study the three-state antiferromagnetic Potts model on the simple-cubic
lattice, paying attention to the surface critical behaviors. When the nearest
neighboring interactions of the surface is tuned, we obtain a phase diagram
similar to the XY model, owing to the emergent O(2) symmetry of the bulk
critical point. For the ordinary transition, we get $y_{h1}=0.780(3)$,
$\eta_{\parallel}=1.44(1)$, and $\eta_{\perp}=0.736(6)$; for the special
transition, we get $y_{s}=0.59(1)$, $y_{h1}=1.693(2)$,
$\eta_{\parallel}=-0.391(4)$, and $\eta_{\perp}=-0.179(5)$; in the
extraordinary-log phase, the surface correlation function $C_{\parallel}(r)$
decays logarithmically, with decaying exponent $q=0.60(2)$, however, the
correlation $C_{\perp}(r)$ still decays algebraically, with critical exponent
$\eta_{\perp}=-0.442(5)$. If the ferromagnetic next nearest neighboring
surface interactions are added, we find two transition points, the first one
is a special point between the ordinary phase and the extraordinary-log phase,
the second one is a transition between the extraordinary-log phase and the
$Z_{6}$ symmetry-breaking phase, with critical exponent $y_{\rm s}=0.41(2)$.
The scaling behaviors of the second transition is very interesting, the
surface spin correlation function $C_{\parallel}(r)$ and the surface squared
staggered magnetization at this point decays logarithmically, with exponent
$q=0.37(1)$; however, the surface structure factor with the smallest wave
vector and the correlation function $C_{\perp}(r)$ satisfy power-law decaying,
with critical exponents $\eta_{\parallel}=-0.69(1)$ and
$\eta_{\perp}=-0.37(1)$, respectively.

###### pacs:

03.67.Bg, 03.65.Ud, 05.30.Rt

##  I Introduction

Phase transition and critical phenomena are hot topics in the research of
condensed matter and statistical physics. At the critical point of a
continuous transition, the system exhibits a variety of singular behaviors
characterized by algebraically decaying of correlation functions. Such type of
behaviors also appear on the surface of the system and can be different to
that of the bulk onesbinder1974 ; binder1983 ; surf1 ; surf2 ; surf3 .
Depending on the strength of the surface interactions, the surface critical
behavior can be classified as “ordinary transition”, “special transition”, and
“extraordinary transition”. Typical examples can be found in the classical
O($n$) modelsising3Dsp ; youjinOnsf ; youjinO4sf . Generally, the ordinary
transition can be found without tuning the surface interactions $J_{\rm s}$,
i.e., it is the same as the bulk ones, if it is tuned (strengthened), a phase
transition may be found, which is the special point $J_{\rm sc}$, and the
phase with $J_{\rm s}>J_{\rm sc}$ is the extraordinary phase. It is clear that
the extraordinary phases of Ising model are ordered, however, whether there is
a special transition for the surface of the Heisenberg model and what the
extraordinary phase is once controversialO3UC ; O3sf2000 ; youjinOnsf .
Similar problems also exist in XY model. Until recently, the research interest
in surface critical phenomena was renewed by the exotic surface critical
behavior in the quantum spin modelsLong2017 ; Ding2018 ; Weber2018 ; Weber2019
. A recent theoretical study pointed out that the extraordinary phase on the
surface of the Heisenberg model may be an “extraordinary-log” phase
characterized by logarithmically decaying of the surface correlation
functionMetlitski2020 , which has been verified numericallyO3sp ; similar
behavior is also found in the XY modelXYlog .

The above results show that the surface critical phenomenon has many
interesting characteristics, which are worth studying, and the related
research are continuingZhu2021 ; Weber2021 ; Toldin2021 ; Ding2021 ; Yu2021 ;
Zhu2021-2 ; Jian2021 ; Max2111 ; Max2111a . In this paper, we study the
surface critical properties of the three-state antiferromagnetic Potts model
on a simple-cubic latticeAF3-1980 ; wsk . Due to the $Z_{3}$ symmetry of the
spins and the permutation symmetry of two sets of sublattices, the ground
state of the system breaks the $Z_{6}$ symmetry, however, at the bulk critical
point, the symmetry of the order parameter is O(2), which is called “emergent
O(2) symmetry”Ding2014 ; Ding2016 , and the corresponding universality class
of the phase transition is the same as the XY model. At such “emergent O(2)”
bulk critical point, we explore the surface critical behaviors by tuning the
nearest neighboring (NN) surface interactions and the next nearest neighboring
(NNN) surface interactions, the two phase diagrams we obtained are shown in
Fig. 1. In the phase diagram (a), the extraordinary-log phase and the special
point are obtained, which is very similar to the XY modelyoujinOnsf ; XYlog ;
we can also see that the surface can not be ordered only by increasing the
strength of the surface NN interactions, this is closely related to the fact
that the two-dimensional antiferromagnetic Potts model is not ordered even the
temperature is down to zeroDing2016 . However, NNN ferromagnetic interactions
can induce $Z_{6}$ symmetry breaking phase at finite temperature for the two-
dimensional antiferromagnetic Potts modelsqNijs , this inspires us to add NNN
ferromagnetic interactions to the surface of the three-dimensional
antiferromagnetic Potts model; the results for this case are summarized in
phase diagram (b), where a new special point with very interesting scaling
behaviors is found.

The paper is arranged as follows: In Sec. II, we introduce the model and
method; in Sec. III, we present the numerical results, including the refined
bulk critical point, the critical behaviors about the ordinary transition, the
extraordinary-log phase, and the two special points; we conclude our paper in
Sec. IV.

Figure 1: Surface phase diagrams of the antiferromagnetic Potts model (1): (a)
the phase diagram with NNN interactions $J^{\prime}_{\rm s}=0$; (b) the phase
diagram with NN interactions $J_{\rm s}=1$.

##  II Model and Method

The antiferromagnetic Potts model we studied is defined on the simple-cubic
lattice

$\displaystyle\mathcal{H}=J\sum\limits_{\langle
i,j\rangle}\delta_{\sigma_{i},\sigma_{j}}+J_{s}\sum\limits_{{\langle
i,j\rangle}^{s}}\delta_{\sigma_{i},\sigma_{j}}-J_{s}^{\prime}\sum\limits_{{\langle\langle
i,j\rangle\rangle}^{\prime}}\delta_{\sigma_{i},\sigma_{j}},$ (1)

$J$, $J_{s}$ and $J_{s}^{\prime}$ are the strengths of the NN bulk
interactions, the NN surface interactions, and the NNN surface interactions,
respectively, noting that the NNN surface interactions are ferromagnetic. The
Potts spin $\sigma_{i}=1,2,3$, which can be mapped to a unit vector in the
plane

$\displaystyle\vec{\sigma}_{i}=(\cos\theta_{i},\sin\theta_{i}),$ (2)

with $\theta_{i}=2\pi\sigma_{i}/3$. Such mapping shows the symmetry of the
spin, the variables we studied are based on such vectors.

The Monte Carlo algorithm we adopt is a combination of the local update
(Metroplis algorithm) and global update (cluster algorithms)wsk . For the case
$J^{\prime}_{s}\neq 0$, the system has both antiferromagnetic and
ferromagnetic interactions, the global update is a mixture of the Swendsen-
WangSW  and WSK algorithmwsk , The high efficiency of the algorithm enables us
to perform simulations for the systems with sizes up to $L=128$. In the
simulations, the periodic boundary condition is applied along the $x$ and $y$
directions, while open boundary condition is applied along the $z$ direction.

The bulk variables we sampled include the squared (staggered) magnetization
$m_{\rm s}^{2}$, the magnetic susceptibility $\chi_{s}$, and the Binder Ratio
$Q_{\rm s}$, which are defined as

$\displaystyle m_{\rm s}^{2}$ $\displaystyle=$
$\displaystyle\langle\mathcal{M}_{\rm s}^{2}\rangle,$ (3)
$\displaystyle\chi_{\rm s}$ $\displaystyle=$ $\displaystyle
N(\langle\mathcal{M}_{\rm s}^{2}\rangle-\langle|\mathcal{M}_{\rm
s}|\rangle^{2}),$ (4) $\displaystyle Q$ $\displaystyle=$
$\displaystyle\frac{\langle\mathcal{M}_{\rm
s}^{2}\rangle^{2}}{\langle\mathcal{M}_{\rm s}^{4}\rangle},$ (5)

where $T$ is the temperature and $\mathcal{M}_{\rm s}$ is defined as

$\displaystyle\mathcal{M}_{\rm s}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum\limits_{\vec{R}}(-1)^{x+y+z}\vec{\sigma}_{\vec{R}}.$
(6)

Here $\vec{R}=(x,y,z)$ is the coordination, $N=L^{3}$ is the number of sites
of the lattice.

We also sample the correlation function $C(r)$ and correlation length $\xi$

$\displaystyle C(r)=\langle\vec{\sigma}_{i}\cdot\vec{\sigma}_{i+r}\rangle,$
(7) $\displaystyle\xi=\frac{(m^{2}_{\rm
s}/F-1)^{1/2}}{2\sqrt{\sum\limits_{i=1}^{d}\sin^{2}(\frac{k_{i}}{2})}}\;,$ (8)

where $\vec{k}$ is the “smallest wavevector” along the $x$ direction–i.e.,
$\vec{k}\equiv(2\pi/L,0,0)$; the “structure factor” $F$ is define as

$\displaystyle
F=\frac{1}{N^{2}}\Big{\langle}\big{|}\sum\limits_{\vec{R}}(-1)^{x+y+z}e^{i\vec{k}\cdot\vec{R}}\vec{\sigma}_{\vec{R}}\big{|}^{2}\Big{\rangle}.$
(9)

Generally, in a critical phase, correlation ratio $\xi/L$ assumes a universal
nonzero value in the thermodynamic limit $L\rightarrow\infty$; in a disordered
phase, correlation length $\xi$ is finite and $\xi/L$ drops to zero, while in
an ordered phase, $\xi/L$ diverges quickly since the structure factor $F$
vanishes rapidly. Therefore, like the Binder Ratio $Q_{\rm s}$, $\xi/L$ is
also very useful in locating the critical point of phase transition.

The definition of the surface variables are very similar to the bulk ones, but
the spins in Eqs. (6), (7), and (9) should be restricted to the surface ones.
For clarity, we add a subscript $1$ to the surface variables, i.e., the
surface squared magnetization, the surface magnetic susceptibility, the
surface Binder Ratio, the surface structure factor, and the surface
correlation length are written as $m_{\rm s1}$, $\chi_{\rm s1}$, $Q_{s1}$,
$F_{1}$, and $\xi_{1}$, respectively, except the surface correlation
functions, which are written as $C_{\parallel}(r)$ and $C_{\perp}(r)$; both
$C_{\parallel}(r)$ and $C_{\perp}(r)$ are calculated as Eq. (7), for
$C_{\parallel}(r)$, both site $i$ and $i+r$ are on the surface; for
$C_{\perp}(r)$, site $i$ is on the surface but site $i+r$ is in the bulk, and
the line from $i$ to $i+r$ is perpendicular to the surface.

##  III Results

###  III.1 Refine the bulk critical point

The simple-cubic antiferromagnetic Potts model has ever been studied by Monte
Carlo simulations in Refs. wsk,  and PottsWL, , whose critical exponents
coincide with the XY universality class, however, the accuracy is not very
high. A more accurate numerical results about the critical exponents can be
found in Ref. Ding2016, , while the interactions studied there are
antiferromagnetic along the $x$ and $y$ directions but ferromagnetic along the
$z$ direction, such mixed interactions does not change the universality class
of the phase transition but lead to different (bulk) critical point. In
current paper, before studying the surface critical behaviors, we refine the
buk critical point for the simple-cubic antiferromagentic model (with full
antiferromagnetic interactions).

Figure 2: Squared magnetization $m^{2}_{\rm s}$ and Binder Ratio $Q_{\rm s}$
of the simple-cubic antiferromagnetic Potts model (1).

As shown in Fig. 2, the bulk transition is clearly indicated by the bulk
squared magnetization $m^{2}_{\rm s}$ and the Binder Ratio $Q_{s}$. It can
also be shown by the bulk correlation ratio $\xi/L$, which is not shown here.
In the vicinity of the critical point, $Q_{\rm s}$ and $\xi/L$ satisfy the
finite-size scaling (FSS)fss1 ; fss2  formula

$\displaystyle Q=a_{0}+\sum\limits_{k=1}^{k_{\rm
max}}a_{k}(T-T_{c})^{k}L^{ky_{t}}+bL^{y_{1}}$ (10)

where $Q$ is $Q_{\rm s}$ or $\xi/L$, and $T_{c}$ is the critical point;
$y_{t}>0$ is the thermal critical exponent, $y_{1}<0$ is the correction-to-
scaling exponent, $a_{0}$, $a_{k}$, and $b$ are unknown parameters.

The data fitting of $Q_{\rm s}$, with fixed $y_{1}=-1$, gives
$T_{c}=1.22602(2)$ and $y_{t}=1.483(10)$; the fitting from $\xi/L$ gives
$T_{c}=1.22603(1)$ and $y_{t}=1.490(3)$; here we chose the better one
$T_{c}=1.22603(1)$ to be our final result, which will be used as the starting
point of the study of surface critical behaviors in the next subsections. We
can also see that the value of the critical exponent $y_{t}$ is consistent
with that of the XY universality classOnloop ; Lv2019 .

###  III.2 Ordinary transition

Exactly at the bulk critical point $T_{c}^{\rm bulk}=1.22603$, we perform
simulations with open boundary condition along the $z$ direction and periodic
boundary conditions along the $x$ and $y$ directions, the surface interactions
are set as $J_{\rm s}=J=1$ and $J^{\prime}_{\rm s}=0$. In this case, the
surface critical behaviors originates from the bulk correlations, which is
called the “ordinary transition”. The data of the surface squared
magnetization $m^{2}_{\rm s1}$, the surface correlation function
$C_{\parallel}$, the surface structure factor $F_{1}$, and the surface
correlation $C_{\perp}$ are fit according to the scaling formulas

$\displaystyle m^{2}_{s1}L^{2}=c+L^{2y_{\rm h1}-2}(a+bL^{y_{1}}),$ (11)
$\displaystyle C_{\parallel}(L/2)=L^{-1-\eta_{\parallel}}(a+bL^{y_{1}}),$ (12)
$\displaystyle F_{1}L^{2}=c+L^{2y_{\rm h1}-2}(a+bL^{y_{1}}),$ (13)
$\displaystyle C_{\perp}(L/2)=L^{-1-\eta_{\perp}}(a+bL^{y_{1}}),$ (14)

where $y_{\rm h1}$, $\eta_{\parallel}$, and $\eta_{\perp}$ are critical
exponent; $y_{1}<0$ is the correction-to-scaling exponents; $a$, $b$, and $c$
are unknown parameters, with $c$ the analytical part of $m^{2}_{\rm s1}L^{2}$
or $F_{1}L^{2}$, originating from the contribution of short-term correlations.
Although mathematically Eqs. (11) and (13) are equivalent to

$\displaystyle m^{2}_{s1}=c\cdot L^{-2}+L^{2y_{\rm h1}-4}(a+bL^{y_{1}}),$ (15)
$\displaystyle F_{1}=c\cdot L^{-2}+L^{2y_{\rm h1}-4}(a+bL^{y_{1}}),$ (16)

respectively, technically, in the fitting of $m^{2}_{\rm s1}$, the RHS of (15)
or (16) will be dominated by the first term and the fitting result of $y_{\rm
h1}$ is often unreliable if the value of $y_{\rm h1}$ is smaller than $1$,
which is exactly the case in ordinary transition.

The data fitting, with $y_{1}=-1$, gives $y_{\rm h1}=0.780(3)$ (from
$m^{2}_{\rm s1}$), $y_{\rm h1}=0.782(4)$ (from $F_{1}$),
$\eta_{\parallel}=1.44(1)$, and $\eta_{\perp}=0.736(6)$, these values coincide
with those of the XY modelyoujinOnsf  and satisfy the scaling laws

$\displaystyle\eta_{\parallel}=d-2y_{h1},$ (17) $\displaystyle
2\eta_{\perp}=\eta_{\parallel}+\eta,$ (18)

with $d=3$ the dimension of the system and $\eta=0.0385$ the exponent of the
bulk correlation.

###  III.3 Special transition and extraordinary phase by tuning the NN surface
interaction $J_{s}$

Figure 3: Surface squared magnetization $m^{2}_{\rm s1}$ and correlation ratio
$\xi_{1}/L$ of the antiferromagnetic Potts model (1), with $T=T_{c}^{\rm
bulk}=1.22603$ and the NNN surface interaction $J^{\prime}_{\rm s}=0$.

Exactly at the bulk critical point $T_{c}^{\rm bulk}=1.22603$, we tuning the
NN surface interaction $J_{s}$ to see whether there is a special transition,
this is confirmed by the behaviors of surface squared magnetization
$m^{2}_{\rm s1}$ and surface correlation ratio $\xi_{1}/L$, as shown in Fig.
3. This transition can also be detected by the Binder Ratio $Q_{\rm s1}$,
which is not shown here. The data of $\xi_{1}/L$ or $Q_{\rm s1}$ in the
vicinity of the special point satisfy the scaling formula (10), with $T$
replaced by $J_{\rm s}$ and the critical exponent $y_{t}$ replaced by $y_{\rm
s}$; the fitting, with $y_{1}=-1$, gives the critical point $J_{\rm
sc}=2.04119(40)$ and the critical exponent $y_{\rm s}=0.59(1)$. We can see
that the value of $y_{s}$ is consistent with that of the XY modelyoujinOnsf .

Exactly at the special transition point $J_{\rm sc}=2.04119$, we investigate
the scaling behaviors of the surface squared magnetization $m^{2}_{\rm s1}$,
the surface correlation function $C_{\parallel}$, the surface structure factor
$F_{1}$, and the surface correlation $C_{\perp}$, they also satisfy the FSS
formulas (11), (12), (13), and (14), respectively. Alternatively, $m^{2}_{\rm
s1}$ and $F_{1}$ can also be fit by Eqs. (15) and (16) respectively, because
in current case $y_{\rm h1}>1.5$, the RHS of the equation is dominated by the
second term, while the first one behaves like a correction-to-scaling term.
The data fitting gives $y_{\rm h1}=1.693(2)$ (from $m^{2}_{\rm s1}$), $y_{\rm
h1}=1.694(4)$ (from $F_{1}$), $\eta_{\parallel}=-0.391(4)$, and
$\eta_{\perp}=-0.179(5)$, these results coincide with those of the XY model
and satisfy the scaling laws (17) and (18).

Figure 4: (a) Log-log plot of $m^{2}_{\rm s1}$, $C_{\parallel}(L/2)$, and
$F_{1}$ versus $\ln(L/L_{0})$, for the antiferromagnetic Potts model (1) with
$T=T_{c}^{\rm bulk}=1.22603$, $J_{\rm s}=5$, and $J^{\prime}_{\rm s}=0$;
$L_{0}=1.54$ is taken for all the three variables. (b) Log-log plot of the
surface correlation $C_{\perp}(L/2)$.

For $J_{\rm s}>J_{\rm sc}$, we find that the system is in the extraordinary-
log phase, characterizing by the logarithmically decaying of correlation
function $C_{\parallel}(r)$ and related variablesMetlitski2020 ; O3sp ; XYlog
. Fig. 4(a) shows a typical case of $J_{s}=5$, and the data can be fit
according to

$\displaystyle m^{2}_{\rm s1}=a\cdot[\ln(L/L_{0})]^{-q},$ (19) $\displaystyle
C_{\parallel}(L/2)=a\cdot[\ln(L/L_{0})]^{-q},$ (20) $\displaystyle
F_{1}=a\cdot[\ln(L/L_{0})]^{-q^{\prime}},$ (21)

where $q$ and $q^{\prime}$ are the critical exponents, $a$ and $L_{0}$ are
nonuniversal parameters. From the fitting of $C_{\parallel}$, we get
$q=0.60(2)$ and $L_{0}=1.54(9)$; from the fitting of $m^{2}_{\rm s1}$ and
$F_{1}$, we get $q=0.61(1)$ and $q^{\prime}=1.61(1)$, respectively. We can see
that the values of $q$ and $q^{\prime}$ are consistent with those of the XY
modelXYlog  and satisfy the relation $q^{\prime}=1+q$. It should be noted that
in the fitting of $m^{2}_{\rm s1}$ and $F_{1}$, we have set the value of
$L_{0}$ fixed to 1.54, such trick is the same as that in Ref. XYlog,  for XY
model. For the surface correlation function $C_{\perp}(L)$, we find that it
still decaying algebraically, as shown in Fig. 4(b), therefore we fit it
according to Eq. (14) with $y_{1}=-1$, we get $\eta_{\perp}=-0.442(5)$.

We also study the case of $J_{s}=10$, and get $q=0.63(2)$,
$q^{\prime}=1.64(2)$, $L_{0}=1.2(1)$, and $\eta_{\perp}=-0.444(7)$. The values
of $q$ and $q^{\prime}$ coincide with the cases of $J_{s}=5$ and the XY
modelXYlog , which means the critical behaviors of the extraordinary-log phase
of the O(2) surface are universal.

In the low temperature, the order parameter of the antiferromagnetic Potts
model on the simple-cubic lattice breaks the $Z_{6}$ symmetryDing2016 ,
therefore one may ask whether the surface can be ordered with such symmetry
breaking when the surface interactions are strengthened. However, by
performing the simulations with much larger $J_{s}$, we find that the surface
is always in the extraordinary-log phase, it can not be ordered by the
strengthen of $J_{s}$. This may be related to the fact that in the limit of
$J_{s}\rightarrow\infty$, the decoupled two dimensional antiferromagnetic
Potts model can not be ordered even the temperature is zeropottssq . However,
as we will show in the next subsection, the $Z_{6}$ symmetry breaking surface
can be reached by adding the NNN ferromagnetic interactions $J^{\prime}_{s}$
in the surface.

###  III.4 New special transition induced by next nearest interactions
$J^{\prime}_{s}$

When the NNN ferromagnetic interactions $J^{\prime}_{s}$ is added to the
surface, we can find two phase transitions, as shown by the behaviors of the
surface susceptibility $\chi_{\rm s1}$ in Fig. 5(a). In the intermediate
region $J^{\prime}_{\rm sc1}<J^{\prime}_{\rm s}<J^{\prime}_{\rm sc2}$, with
$J^{\prime}_{\rm s1}\approx 0.41$ and $J^{\prime}_{\rm sc2}\approx 0.85$, the
scaling behaviors of $\chi_{\rm s1}$ is very different from that in region
$J^{\prime}_{\rm s}<J^{\prime}_{\rm sc1}$ or $J^{\prime}_{\rm
s}>J^{\prime}_{\rm sc2}$. As we will shown later, in this region, the system
is in the extraordinary-log phase. Therefore, the peak of $\chi_{\rm s1}$ at
$J^{\prime}_{\rm sc1}\approx 0.41$ indicates a special transition point like
that studied in above subsection.

Figure 5: Surface critical behaviors of the antiferromagnetic Potts model (1),
with $T=T_{c}=1.22603$ and $J_{\rm s}=1$: (a) surface susceptibility
$\chi_{\rm s1}$, (b) surface squared magnetization $m^{2}_{\rm s1}$, (c)
surface structure factor $F_{1}$, (d) surface Binder Ratio $Q_{\rm s1}$.
Figure 6: (a) Log-log plot of $F_{1}$ versus $\ln(L/L_{0})$, using the values
of $L_{0}$ fit from $m^{2}_{s1}$, listed in Table 1; (b) log-log plot of
$F_{1}$ versus $L$.

The second transition at $J^{\prime}_{\rm sc2}$, which we call “new special
point”, is revealed not only by the surface susceptibility $\chi_{s1}$, but
also by the surface squared magnetization $m^{2}_{\rm s1}$ and surface
structure factor $F_{1}$, as shown in Fig. 5(b) and (c), respectively;
however, the signature of the Binder Ratio $Q_{\rm s1}$ for this transition is
not obvious, as shown in Fig. 5(d).

We can see that the second peak of $\chi_{\rm s1}$ diverges much slower than
the first one, namely the singularity of such transition is very weak. The
transition point $J^{\prime}_{\rm sc2}\approx 0.85$ is obtained by the
position of such peak of the largest system with size $L=128$; however,
considering the finite-size effect, the critical point in the thermodynamic
limit should be a little smaller than this one. A more accurate estimation of
the critical point can be obtained from the scaling behavior of the structure
factor $F_{1}$, which satisfies the logarithmic decaying in the region
$J^{\prime}_{\rm sc1}<J^{\prime}_{\rm s}<J^{\prime}_{\rm sc2}$ but a power-law
scaling at the transition point $J^{\prime}_{\rm sc2}=0.815$; this is
demonstrated in Fig. 6.

In the region $J^{\prime}_{\rm sc1}<J^{\prime}_{\rm s}<J^{\prime}_{\rm sc2}$,
we check the FFS behaviors of $F_{1}$ for several cases, the results of data
fitting, according to Eq. (21), are listed in Table. 1. In this region,
$m^{2}_{\rm s1}$ and $C_{\parallel}$ also satisfy the logarithmic decaying,
the fitting results, according to Eq. (19) and (20), are also listed in Table.
1. These results claim that in this region the system is in an extraordinary-
log phase; the critical exponents $q$ and $q^{\prime}$ in this phase are
universal and satisfy the relation $q^{\prime}=1+q$, as that in XY modelXYlog
. In this phase, the correlation $C_{\perp}$ still satisfy the power-law
scaling formula (14), the fitting results of $\eta_{\perp}$ for several cases
are also listed in Table. 1.

An important point should be emphasized about the extraordinary-log phase is
that the surface susceptibility $\chi_{\rm s1}$ is divergent in this phase,
such property is very different from that of the ordinary phase or the
extraordinary phase with long-range order. Technically, this property can help
us to detect whether a surface is in an extraordinary-log phase. It should be
noted that, currently, it not a trivial work to distinguish an extraordinary-
log phase from a long-range ordered phase numerically without theoretical
precognition, especially when the value of the long-range order is relatively
small.

Table 1: Fitting results of $m^{2}_{\rm s1}$, $C_{\parallel}$, and $F_{1}$, where $m^{2}_{\rm s1}$ and $C_{\parallel}$ are fit according to Eqs. (19) and (20), respectively; $F_{1}$ is fit according to (21) when $J^{\prime}_{\rm s}<J^{\prime}_{sc2}$ but a power-law formula $aL^{-1-\eta_{\parallel}}$ when $J^{\prime}_{\rm s}<J^{\prime}_{sc2}$. In fitting of $C_{\parallel}$ and $m^{2}_{\rm s1}$, we let $L_{0}$ free to be fit, while in the fitting of $F_{1}$, we fix $L_{0}$ as that of $m^{2}_{\rm s1}$.  $J^{\prime}_{\rm s}$ |  $q$ ($C_{\parallel}$) | $L_{0}$ |  $q$ ($m^{2}_{\rm s1}$) | $L_{0}$ | $q^{\prime}$ | $\eta_{\perp}$
---|---|---|---|---|---|---
0.55 | 0.61(3) | 2.0(7) | 0.59(3) | 2.63(9) | 1.56(4) | -0.43(1)
0.65 | 0.56(2) | 1.36(10) | 0.60(2) | 0.89(5) | 1.59(1) | -0.43(1)
0.75 | 0.55(2) | 0.90(16) | 0.57(1) | 0.44(9) | 1.57(1) | -0.38(1)
$J^{\prime}_{\rm sc2}$ |  $q$ ($C_{\parallel}$) | $L_{0}$ |  $q$ ($m^{2}_{\rm s1}$) | $L_{0}$ | $\eta_{\parallel}$ | $\eta_{\perp}$
0.815 | 0.37(1) | 2.3(3) | 0.38(1) | 1.13(9) | $-0.69(1)$ | -0.37(1)
Figure 7: (a) $F_{1}^{-1-\eta}=F_{1}^{0.31}$ versus $J^{\prime}_{\rm s}$; the
dashed line is the critical point; (b) data collaps: $F_{1}L^{0.31}$ versus
$(J^{\prime}_{\rm s}-0.815)L^{y_{s}}$, with $y_{s}=0.41$.

At the transition point $J^{\prime}_{sc2}=0.815$, $F_{1}$ satisfy a power-law
formula $aL^{-1-\eta_{\parallel}}$, where $\eta_{\parallel}=-0.69(1)$. The
correlation function $C_{\perp}$ also satisfies the power law (14), and the
critical exponent $\eta_{\perp}=-0.37(1)$. However, the $m^{2}_{\rm s1}$ and
$C_{\parallel}$ at this point still satisfy the logarithmically decaying
formulas (19) and (20), with exponent $q=0.30(1)$, as shown in Table 1.

In the range $J^{\prime}_{\rm s}>J^{\prime}_{\rm sc2}$, the surface is in a
long-range ordered phase, the structure factor $F_{1}$ decays much faster,
which satisfies neither a logarithmically decaying formula nor a simple pow-
law formula, this is also demonstrated in Fig. 6.

In order to further explore the critical behaviors of the new special point,
we plot $F_{1}L^{-1-\eta_{\parallel}}=F_{1}L^{0.31}$ versus $J^{\prime}_{\rm
s}$ in Fig. 7(a), we can see that $F_{1}L^{-1-\eta_{\parallel}}$ plays the
role of a dimensionless variable like the Binder Ratio; Fig. 7(b) is the data
collaps plot of $F_{1}L^{0.31}$ versus $(J^{\prime}_{\rm s}-0.815)L^{y_{\rm
s}}$, with critical exponent $y_{\rm s}=0.41$. More rigorously, one can fit
the data of $F_{1}$ in the vicinity of $J^{\prime}_{\rm s}=0.815$ by the
scaling formula

$\displaystyle
F_{1}=L^{-1-\eta_{\parallel}}\big{[}a_{0}+\sum\limits_{k=1}^{k_{\rm
max}}a_{k}(J^{\prime}_{\rm s}-J^{\prime}_{\rm
sc2})^{k}L^{ky_{s}}+bL^{y_{1}}\big{]},$ (22)

which gives $\eta_{\parallel}=-0.69(1)$, $J^{\prime}_{\rm sc2}=0.815(2)$, and
$y_{\rm s}=0.41(2)$, these results are consistent with those applied in the
data collaps in Fig. 7(b), thus give a self-consistent check.

In summary, for the new special transition, it is mainly characterized by the
the following critical behaviors,

$\displaystyle J^{\prime}_{\rm sc2}=0.815(2),$ (23) $\displaystyle
y_{s}=0.41(2),$ (24) $\displaystyle F_{1}\sim aL^{-1-\eta_{\parallel}},{\rm
with}\ \eta_{\parallel}=-0.69(1),$ (25) $\displaystyle m^{2}_{\rm s1}\sim
a\cdot[\ln(L/L_{0})]^{-q},{\rm with}\ q=0.37(1),$ (26) $\displaystyle
C_{\parallel}\sim a\cdot[\ln(L/L_{0})]^{-q},{\rm with}\ q=0.38(1),$ (27)
$\displaystyle C_{\perp}\sim aL^{-1-\eta_{\perp}},{\rm
with}\quad\eta=-0.37(1).$ (28)

Figure 8: Histograms of the surface staggered magnetization of the
antiferromagnetic Potts model (1); in all the cases, the strength of the
surface NN interactions $J_{\rm s}=1$.

###  III.5 Symmetries of the surface

After studied the phase transitions, we give a short description about the
symmetries of the model. The three-state antiferromagnetic Potts model is well
known for the emergent O(2) symmetry at the bulk critical point, although the
ground state is $Z_{6}$ symmetry breaking, which is ordered by entropy. It is
an interesting question what the surface symmetry is.

As shown in Fig. 8, in the ordinary phase ($J^{\prime}_{\rm s}=0$), the
histogram visualizes the disorder properties of the surface staggered
magnetization, this confirms that the surface is disorder, and the critical
behaviors are purely induced by the bulk criticality. For the extraordinary-
log phase, we show two cases; for the case of $J^{\prime}_{\rm s}=0.65$, we
can see that symmetry of the order parameter is O(2) when the system size
reaches $L=32$, although we can still find the shadow of the $Z_{6}$ symmetry
when the system size is small ($L$=16); for the case of $J^{\prime}_{\rm
s}=0.75$, we can see that the $Z_{6}$ symmetry persists to show its effect
when the system size is $L=32$, however, the effect is relatively weak
comparing to that of $L=16$, thus we infer that in the thermodynamic limit,
the symmetry should be O(2). Such property is important, it is also helpful
for us to numerically distinguish the extraordinary-log phase from an ordered
surface with small value of long-range order.

When $J^{\prime}_{\rm s}$ is large enough, the surface is ordered and the
$Z_{6}$ symmetry is very obvious, this is also shown in Fig. 8 with
$J^{\prime}_{\rm s}=1.0$.

##  IV Conclusion and discussion

In summary, we have studied the surface critical behaviors of the three-state
antiferromanetic Potts model on the simple-cubic lattice, we obtain a phase
diagram similar to the XY model by tuning the NN surface interactions, the
universality classes of the special transition and extraordinary-log phase are
the same as the XY model. A long range-order surface which breaks the $Z_{6}$
symmetry and a new type of special transition between such ordered phase and
the extraordinary-log phase is obtained by tuning the NNN ferromagnetic
surface interactions. At such new special transition point, the scaling
behaviors are very interesting; the surface squared magnetization $m^{2}_{\rm
s1}$ and the surface correlation function $C_{\parallel}$ satisfy the
logarithmic decaying, with exponent $q=0.37(1)$, but the structure factor
$F_{1}$ and the correlation function $C_{\perp}$ still satisfy the power-law
decaying, with critical exponent $\eta_{\parallel}=-0.69(1)$ and
$\eta_{\perp}=-0.37(1)$, respectively. We also visualized the symmetries of
different phases of the surface, where the extraordinary-log phase is shown to
conserve O(2) symmetry.

We have to stress that all our conclusions about the new special transition
are based on the assumption that it is a second-order phase transition, which
does not rule out the possibility of other types of phase transitions, such as
the Berezinskii-Kosterlitz-Thouless transition.

In classical XY model, although the symmetry of the bulk critical point is
also O(2), the surface can not be ordered because of Mermin-Wagner-Hohenberg
theoremMW1966 ; H1967 , therefore such new type of special transition can not
appear in classical XY model; discrete Hamiltonian symmetry and emergent
O($n$) symmetry of bulk criticality seems a necessary condition for such new
type of special transition, such as the clock modelclock6 . In addition to the
emergent O(2) critical point, other emergent O($n$) critical point have also
been found in systems with discrete symmetry of HamiltonianDing2016 ;
Leonard2015 . It is obviously an important and interesting question to study
such new special transition and related critical properties in these models.

##  V Acknowledgment

We thank Jian-Ping Lv for valuable discussions, and Cenke Xu for the advice of
visualizing the symmetries of the sufrace by histograms. C.D. is supported by
the National Science Foundation of China under Grants Numbers 11975024 and
62175001, the Anhui Provincial Supporting Program for Excellent Young Talents
in Colleges and Universities under Grant Number gxyqZD2019023. L.Z. is
supported by the National Key R&D Program (2018YFA0305800), the National
Natural Science Foundation of China under Grants Numbers 11804337 and
12174387, CAS Strategic Priority Research Program (XDB28000000) and CAS Youth
Innovation Promotion Association. Y.D. is supported by the National Natural
Science Foundation of China under Grant No. 11625522, the Science and
Technology Committee of Shanghai under grant No. 20DZ2210100, and the National
Key R&D Program of China under Grant No. 2018YFA0306501.

## References

  * (1) K. Binder and P. C. Hohenberg, Surface effects on magnetic phase transitions, Phys. Rev. B 9, 2194 (1974).
  * (2) K. Binder, in Phase Transitions and Critical Phenomena, Vol. 8, edited by C. Domb and J. L. Lebowitz (Academic Press, London, England, 1983).
  * (3) T. C. Lubensky and M. H. Rubin, Critical phenomena in semi-infinite systems. I. $\epsilon$ expansion for positive extrapolation length, Phys. Rev. B 11, 4533 (1975).
  * (4) T. C. Lubensky and M. H. Rubin, Critical phenomena in semi-infinite systems. II. Mean-field theory, Phys. Rev. B 12, 3885 (1975).
  * (5) M. N. Barber, Scaling Relations for Critical Exponents of Surface Properties of Magnets, Phys. Rev. B 8, 407 (1973).
  * (6) M. Hasenbusch, Monte Carlo study of surface critical phenomena: The special point, Physical Review B 84, 134405 (2011).
  * (7) Y. Deng, H. W. J. Blote, and M. P. Nightingale, Surface and bulk transitions in three-dimensional O($n$) models, Phys. Rev. E 72, 016128 (2005).
  * (8) Y. Deng, Bulk and surface phase transitions in the three-dimensional O(4) spin model, Phys. Rev. E 73, 056116 (2006).
  * (9) C. S. Arnold and D. P. Pappas, Gd(0001): A Semi-Infinite Three-Dimensional Heisenberg Ferromagnet with Ordinary Surface Transition, Phys. Rev. Lett. 85, 5202 (2000).
  * (10) M. Krech, Surface scaling behavior of isotropic Heisenberg systems: Critical exponents, structure factor, and profiles, Phys. Rev. B 62, 6360 (2000).
  * (11) L. Zhang and F. Wang, Unconventional surface critical behavior induced by a quantum phase transition from the two-dimensional Affleck-Kennedy-Lieb-Tasaki phase to a Néel-ordered phase, Phys. Rev. Lett. 118, 087201 (2017).
  * (12) C. Ding, L. Zhang, and W. Guo, Engineering Surface Critical Behavior of (2 + 1)-Dimensional O(3) Quantum Critical Points, Phys. Rev. Lett. 120, 235701 (2018).
  * (13) L. Weber, F. Parisen Toldin, and S. Wessel, Nonordinary edge criticality of two-dimensional quantum critical magnets, Phys. Rev. B 98, 140403(R) (2018).
  * (14) L. Weber and S. Wessel, Nonordinary criticality at the edges of planar spin-1 Heisenberg antiferromagnets, Phys. Rev. B 100, 054437 (2019).
  * (15) M. A. Metlitski, Boundary criticality of the O(N) model in d=3 critically revisited, arXiv:2009.05119 (2020).
  * (16) F. P. Toldin, Boundary Critical Behavior of the Three-Dimensional Heisenberg Universality Class, Phys. Rev. Lett. 126, 135701 (2021).
  * (17) M. Hu, Y. Deng, and J.-P. Lv, Extraordinary-Log Surface Phase Transition in the Three-Dimensional XY Model,Phys. Rev. Lett. 127, 120603 (2021).
  * (18) W. Zhu, C. Ding, L. Zhang, and W. Guo, Surface critical behavior of coupled Haldane chains, Phys. Rev. B 103, 024412 (2021).
  * (19) L. Weber and S. Wessel, Spin versus bond correlations along dangling edges of quantum critical magnets, Phys. Rev. B 103, L020406 (2021).
  * (20) C.-M. Jian, Y. Xu, X.-C. Wu, and C. Xu, Continuous Néel-VBS quantum phase transition in non-local one-dimensional systems with SO(3) symmetry, SciPost Phys. 10, 033 (2021).
  * (21) F. P. Toldin, M. A. Metlitski, Boundary criticality of the 3d O(N) model: from normal to extraordinary, arXiv:2111.03613 (2021).
  * (22) C. Ding, W. Zhu, W.-A. Guo, L. Zhang, Special Transition and Extraordinary Phase on the Surface of a (2+1)-Dimensional Quantum Heisenberg Antiferromagnet, arXiv:2110.04762 (2021).
  * (23) X.-J. Yu, R.-Z. Huang, H.-H. Song, L. Xu, C. Ding, L. Zhang, Conformal Boundary Conditions of Symmetry-Enriched Quantum Critical Spin Chains, arXiv:2111.10945 (2021).
  * (24) W. Zhu, C. Ding, L. Zhang, W. Guo, Exotic surface behaviors induced by geometrical settings of two-dimensional dimerized quantum XXZ model, arXiv:2111.12336 (2021).
  * (25) F. P. Toldin, M. A. Metlitski, Boundary criticality of the 3d O(N) model: from normal to extraordinary, arXiv:2111.03613 (2021).
  * (26) J. Padayasi, A. Krishnan, M. A. Metlitski, I. A. Gruzberg, M. Meineri, The extraordinary boundary transition in the 3d O(N) model via conformal bootstrap, arXiv:2111.03071 (2021).
  * (27) J. R. Banavar, G. S. Grest, and D. Jasnow, Ordering and phase Transitions in antiferromagnetic potts models, Phys. Rev. Lett. 45, 1424 (1980).
  * (28) J.-S. Wang, R. H. Swendsen, and R. K. Kotecký, Three-state antiferromagneetic Potts mmodels: A monte Carlo study, 42, 2465 (1990).
  * (29) C. Yamaguchi and Y. Okabe, Three-dimensional antiferromagnetic q-state Potts models: application of the Wang-Landau algorithm, J. Phys. A: Math. Gen. 34, 8781-8794 (2001).
  * (30) C. Ding, W. Guo, and Y. Deng, Reentrance of Berezinskii-Kosterlitz-Thouless-like transitions in a three-state Potts antiferromagnetic thin film, Phys. Rev. B 90, 134420 (2014).
  * (31) Chengxiang Ding, Henk W. J. Blöte, and Youjin Deng, Emergent O($n$) symmetry in a series of three-dimensional Potts models, Phys. Rev. B 94, 104402 (2016).
  * (32) J. Salas and A. D. Sokal, The three-state square-lattice Potts antiferromagnet at zero temperature, J. Stat. Phys. 92, 729 (1998).
  * (33) M. P. M. den Nijs, M. P. Nightingale, and M. Schick, Critical fans in the antiferromagnetic three-state Potts model, Phys. Rev. B 26, 2490 (1982).
  * (34) Robert H. Swendsen and Jian-Sheng Wang, Nonuniversal critical dynamics in Monte Carlo simulations, Phys. Rev. Lett. 58, 86 (1987).
  * (35) M. P. Nightingale, in Finite-Size Scaling and Numerical Simulation of Statistical Systems, edited by V. Privaman (World Scientific, Singapore, 1990).
  * (36) M. N. Barber, in Phase Transitions and Critical Phenomena, Vol. 8, edited by C. Domb and J. L. Lebowitz (Academic Press, New York, 1983).
  * (37) Q. Liu, Y. Deng, T. M. Garoni, and H. W. J. Blöte, The O($n$) loop model on a three-dimensional lattice, Nucl. Phys. B 859, 107 (2012).
  * (38) W. Xu, Y. Sun, J.-P. Lv, and Y. Deng, High-precision Monte Carlo study of several models in the three-dimensional U(1) universality class, Phys. Rev. B 100, 064525 (2019).
  * (39) N. D. Mermin and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models, Phys. Rev. Lett. 17, 1133 (1966).
  * (40) P. C. Hohenberg, Existence of long-range order in one and two dimensions, Phys. Rev. 158, 383 (1967).
  * (41) F. Léonard and B. Delamotte, Critical Exponents Can Be Different on the Two Sides of a Transition: A Generic Mechanism, Phys. Rev. Lett. 115, 200601
  * (42) X. Zou, S. Liu, W. Guo, Surface critical properties of the three-dimensional clock model, arXiv:2204.13612.