# First-order formalism and thick branes in mimetic gravity

Qun-Ying Xie<EMAIL_ADDRESS>Qi-Ming Fu<EMAIL_ADDRESS>Tao-Tao Sui
<EMAIL_ADDRESS>Li Zhao<EMAIL_ADDRESS>Yi Zhong<EMAIL_ADDRESS>School of Information Science and Engineering, Lanzhou University, Lanzhou
730000, P. R. China  Lanzhou Center for Theoretical Physics $\&$ Research
Center of Gravitation, Lanzhou University, Lanzhou 730000, P. R. China
Institute of Physics, Shaanxi University of Technology, Hanzhong 723000, P. R.
China  Hunan Provincial Key Laboratory of High-Energy Scale Physics and
Applications, Hunan University, Changsha 410082, P. R. China

###### Abstract

In this paper, we investigate thick branes generated by a scalar field in
mimetic gravity theory. By introducing two auxiliary super-potentials, we
transform the second-order field equations of the system into a set of first-
order equations. With this first-order formalism, several types of analytical
thick brane solutions are obtained. Then, tensor and scalar perturbations are
analysed. We find that both kinds of perturbations are stable. The effective
potentials for the tensor and scalar perturbations are dual to each other. The
tensor zero mode can be localized on the brane while the scalar zero mode
cannot. Thus, the four-dimensional Newtonian potential can be recovered on the
brane.

###### keywords:

$Thick~{}brane,~{}~{}Mimetic~{}graivty$

††journal: Physics Letters B

##  1 Introduction

Modified gravity theories have obtained great development and performance in
the study of some unsolved problems in general relativity such as the dark
energy problem, the dark matter problem, the singularity problem, etc. By
isolating the conformal degree of freedom of general relativity, Chamseddine
and Mukhanov proposed a theory called mimetic gravity [1]. This theory was
studied from the view of variational principle [2]. It was shown that after
introducing the scalar potential this conformal degree of freedom becomes
dynamical and can mimic cold dark matter [1, 3] or dark energy [4, 5, 6, 7]
and can resolve the singularity problem [8] and the cosmic coincidence problem
[9]. The extension of mimetic gravity was used to investigate the inflationary
solution [10]. In Ref. [11], the authors proposed mimetic Einstein-Cartan
gravity and proved that torsion is a non-propagating field in this mimetic
gravity. Mimetic gravity theory was also extended to Horava-like theory and
applied to galactic rotation curves [12]. It was also applied to other gravity
theories such as $f(R)$ gravity [13, 14, 15, 16, 17, 18, 19], Horndeski
gravity [20, 21, 22] and Gauss-Bonnet gravity [23, 24].

On the other hand, in order to solve the gauge hierarchy problem and the
cosmological constant problem, Randall and Sundrum (RS) proposed that our
four-dimensional world could be a brane embedded in five-dimensional space-
time [25]. With the warped extra dimension, it was further found that the size
of extra dimension can be infinitely large without conflicting with Newtonian
gravitational law [26]. This charming idea has attracted substantial
researches in particle physics, cosmology, gravity theory, and other related
fields [27, 28, 29, 30, 31, 32, 33, 34, 35, 36].

Recently, one of the interesting works appeared in Ref. [37], which applied
mimetic gravity theory to the thin RSII brane model [26]. It was shown that
the mimetic scalar field can mimic the dark sectors on the brane and explain
the late time cosmic expansion in the favor of observational data, and it has
the capability to explain initial time cosmological inflation [37]. Later,
other related topics about thick branes in mimetic gravity were studied in
Refs. [40, 38, 41, 39]. Thick branes with the inner structure and the
stability of the perturbations were first investigated in Ref. [40]. Besides,
it is known that first-order formalism is a very powerful tool to obtain
analytical brane solutions [42, 43, 44]. With this formalism, the second-order
coupled field equations can be written as a set of first-order ones by
introducing one or more auxiliary super-potentials. Very recently, Bazeia et
al. used first-order formalism to find brane solutions in mimetic gravity [41,
45]. In this paper, we would like to investigate thick branes in mimetic
gravity with the help of first-order formalism by including two super-
potentials. In order to show the systematicness and effectiveness of the
first-order formalism on finding analytical brane solutions, we will utilize
polynomial, period, and mixed super-potentials. We will also investigate the
stabilities of the tensor and scalar perturbations as well as the relationship
between the localizations of the tensor and scalar zero modes.

The paper is organized as follows. In Sec. 2, we give a review of the thick
brane model and reduce the second-order field equations to the first-order
ones by introducing two auxiliary super-potentials. In Sec. 3, we obtain three
types of analytical brane solutions by considering different forms of super-
potentials. In Sec. 4 and Sec. 5, we focus on tensor and scalar perturbations,
respectively. Finally, the conclusion and discussion are given in Sec. 6.

##  2 First-order formalism for thick brane models

We consider thick branes generated by a scalar field in five-dimensional
mimetic gravity. The corresponding action is given by

$\displaystyle S\\!=\\!\int
d^{4}xdy\sqrt{-g}\left(\frac{R}{2}+L_{\phi}\right),$ (1)

where $R$ is the five-dimensional scalar curvature, the Lagrangian of the
mimetic scalar field $\phi$ is

$\displaystyle
L_{\phi}=\lambda\left[g^{MN}\partial_{M}\phi\partial_{N}\phi-U(\phi)\right]-V(\phi).$
(2)

Here $\lambda$ represents the Lagrange multiplier, $U(\phi)$ and $V(\phi)$ are
two potentials. In this paper, $x^{M}$ and $x^{\mu}$ denote respectively the
five-dimensional bulk coordinates and the four-dimensional brane ones, where
the indices $M,N,\cdots=0,1,2,3,5$ and $\mu,\nu,\cdots=0,1,2,3$.

The variations of the action (1) with respect to the background metric
$g_{MN}$, the scalar field $\phi$, and the Lagrange multiplier $\lambda$ lead
to the following field equations

$\displaystyle G_{MN}+2\lambda\partial_{M}\phi\partial_{N}\phi-L_{\phi}g_{MN}$
$\displaystyle=$ $\displaystyle 0,$ (3) $\displaystyle
2\lambda\Box^{(5)}\phi+2\nabla_{M}\lambda\nabla^{M}\phi+\lambda\frac{\partial
U}{\partial\phi}+\frac{\partial V}{\partial\phi}$ $\displaystyle=$
$\displaystyle 0,$ (4) $\displaystyle
g^{MN}\partial_{M}\phi\partial_{N}\phi-U(\phi)$ $\displaystyle=$
$\displaystyle 0.$ (5)

Here the five-dimensional d’Alembert operator is defined as
$\Box^{(5)}=g^{MN}\nabla_{M}\nabla_{N}$.

The line-element ansatz for a flat brane is given by

$\displaystyle ds^{2}\\!\\!$ $\displaystyle=\\!\\!$ $\displaystyle
e^{2A(y)}\eta_{\mu\nu}dx^{\mu}dx^{\nu}+dy^{2}$ (6) $\displaystyle=\\!\\!$
$\displaystyle a^{2}(y)\eta_{\mu\nu}dx^{\mu}dx^{\nu}+dy^{2}.$ (7)

Using the line-element (6), and considering that the scalar field is static
and depends only on the extra-dimensional coordinate, we can get the following
second-order nonlinear coupled differential field equations

$\displaystyle 6A^{\prime 2}+\lambda\Big{(}U(\phi)+\phi^{\prime
2}\Big{)}+V(\phi)$ $\displaystyle=$ $\displaystyle 0,$ (8) $\displaystyle
6A^{\prime 2}+3A^{\prime\prime}+\lambda\Big{(}U(\phi)-\phi^{\prime
2}\Big{)}+V(\phi)$ $\displaystyle=$ $\displaystyle 0,$ (9)
$\displaystyle\lambda\left(8A^{\prime}\phi^{\prime}+2\phi^{\prime\prime}+\frac{\partial
U}{\partial\phi}\right)+2\lambda^{\prime}\phi^{\prime}+\frac{\partial
V}{\partial\phi}$ $\displaystyle=$ $\displaystyle 0,$ (10)
$\displaystyle\phi^{\prime 2}-U(\phi)$ $\displaystyle=$ $\displaystyle 0.$
(11)

Here, the primes denote derivatives with respect to the extra-dimensional
coordinate $y$. It can be seen that it is not easy to analytically solve the
above second-order field equations directly. However, one can reduce them to
first-order field equations by introducing two auxiliary super-potentials [41]

$\displaystyle Q=Q(\phi),~{}~{}~{}~{}W=W(\phi),$ (12)

and providing the potential

$\displaystyle V(\phi)=Q_{\phi}W_{\phi}-\frac{2}{3}W^{2},$ (13)

where $Q_{\phi}=\frac{dQ}{d\phi}$ and $W_{\phi}=\frac{dW}{d\phi}$. The
resulting first-order field equations can be written as

$\displaystyle A^{\prime}\\!\\!$ $\displaystyle=\\!\\!$
$\displaystyle-\frac{1}{3}W(\phi),$ (14) $\displaystyle\phi^{\prime}\\!\\!$
$\displaystyle=\\!\\!$ $\displaystyle Q_{\phi},$ (15) $\displaystyle
U(\phi)\\!\\!$ $\displaystyle=\\!\\!$ $\displaystyle Q^{2}_{\phi},$ (16)
$\displaystyle\lambda(\phi)\\!\\!$ $\displaystyle=\\!\\!$
$\displaystyle-\frac{1}{2}\frac{W_{\phi}}{Q_{\phi}}.$ (17)

These equations would be helpful to give analytical brane solutions. One can
see from (16) that the scalar potential $U(\phi)$ is related with the super-
potential $Q(\phi)$, which is a function of the scalar field $\phi$. The other
potential $V(\phi)$ is determined by the two super-potentials $Q$ and $W$ (see
Eq. (13)).

Note that the first-order equations (14)-(16) can be divided into two groups.
One is Eq. (14) related with $A$ and $W$, the other is Eqs. (15) and (16)
related with $\phi,\,Q,\,U$. Therefore, these equations can be solved with
different approaches by giving different combinations from ($A,\,W$) and
$(\phi,\,Q,\,U)$. For example, we can choose $W(\phi)$ and $Q(\phi)$, or
$W(\phi)$ and $\phi(y)$.

The energy density of the brane is given by

$\displaystyle\rho(y)=T_{MN}u^{M}u^{N},$ (18)

where $u^{M}=(u^{0},0,0,0,0)$ is the velocity of a static observer. From the
condition of the velocity $g_{MN}u^{M}u^{N}=-1$, we have
$u^{0}=\text{e}^{-A}$. Thus, the energy density can be written as

$\displaystyle\rho(y)=\frac{1}{2}V(\phi(y)),$ (19)

which shows that the potential $V(\phi)$ and the profile of the scalar
$\phi(y)$ determine the distribution of the thick brane along the extra
dimension.

On the other hand, in order to localize gravity on the brane, the warp factor
$\text{e}^{2A(y)}$ should tend to zero rapidly enough as
$y\rightarrow\pm\infty$, such that the condition $\int e^{2A}dy<\infty$ is
satisfied. This will be derived in section 4. Usually, we consider the
solutions with $e^{2A(y)}|_{y\rightarrow\pm\infty}\rightarrow e^{-2k|y|}$,
which corresponds to the branes embedded in an AdS spacetime with a negative
cosmological constant. It should be pointed out that the contribution of the
cosmological constant has been included in the energy density (19) for this
case. Therefore, the net energy density should be

$\displaystyle\rho(y)=\frac{1}{2}V(\phi(y))-\Lambda.$ (20)

##  3 The thick brane solutions

Next, we will take focus on finding some specific analytical solutions of the
thick brane model by solving the first-order equations (14)-(17). Our main
motivation here is to show the systematicness and effectiveness of the first-
order formalism, which is also called the super-potential method [34, 35, 46].

###  3.1 Polynomial super-potentials

####  3.1.1 Solution I

First, supposing that one of the super-potential $Q$ has the following
polynomial form considered in Ref. [41]:

$\displaystyle Q(\phi)=k\left(v\phi-\frac{\phi^{3}}{3v}\right),$ (21)

and solving Eq. (15), we can easily get the solution of the scalar field
$\phi$

$\displaystyle\phi(y)=v\,\text{tanh}(ky).$ (22)

The solution of the potential $U$ can be read directly from Eq. (16) as

$\displaystyle U(\phi)=\frac{k^{2}}{v^{2}}\left(\phi^{2}-v^{2}\right)^{2}.$
(23)

The other super-potential $W$ is chosen as

$\displaystyle W(\phi)=\frac{3kn}{v}\phi,$ (24)

where $n$ is a non-vanishing parameter. Then, from Eq. (14), we obtain a
simple solution for the warp factor

$\displaystyle A(y)=\ln\text{sech}^{n}(ky).$ (25)

The other functions are given by

$\displaystyle\lambda(y)$ $\displaystyle=$
$\displaystyle-\frac{3n}{2v^{2}}\text{cosh}^{2}(ky),$ (26) $\displaystyle
V(\phi)$ $\displaystyle=$
$\displaystyle\frac{3nk^{2}}{v^{2}}\left[v^{2}-(2n+1)\phi^{2}\right].$ (27)

The energy density of the brane including the cosmological constant is

$\displaystyle\rho^{(\Lambda)}=\frac{V}{2}=3k^{2}n\Big{(}1-(2n+1)\tanh^{2}(ky)\Big{)}.$
(28)

From the solution (25), we can calculate the cosmological constant and hence
the net energy density

$\displaystyle\rho(y)=\rho^{(\Lambda)}-\Lambda=\frac{3}{2}k^{2}n(2n+1)\text{sech}^{2}(ky).$
(29)

The shapes of the warp factor $a(y)=\text{e}^{A(y)}$ and the energy density
$\rho(y)$ are plotted in Fig. 1, which shows that the parameter $n$ affects
the warp factor and the energy density. With the increase of $n$, the warp
factor becomes narrower while the energy density becomes larger and narrower.
The maximum of the energy density is given by
$\rho_{\text{max}}=\frac{3}{2}k^{2}n(2n+1)$ for $n>0$ or $n<-1/2$. It is
obvious that the parameter $v$ does not affect the warp factor and the energy
density, it only affects the amplitude of the scalar field $\phi$ and hence
the localization of a bulk fermion $\Psi$ when one introduces the Yukawa
coupling $\eta\bar{\Psi}\phi\Psi$ [47]. In fact, $v$ is the vacuum expectation
value of the scalar potential $U(\phi)$ given in (23).

(a) The warp factor (b) The energy density

Figure 1: The shapes of the warp factor $a(y)$ and the energy density
$\rho(y)$ of the first brane model. The parameter $n$ is set as $n=1,\,3,\,5$
for the thick red, dashed blue, thin black lines, respectively.

####  3.1.2 Solution II

Next, we also consider the same polynomial super-potential $Q$ as in previous
subsection

$\displaystyle Q(\phi)\\!\\!$ $\displaystyle=\\!\\!$ $\displaystyle
k\left(v\phi-\frac{\phi^{3}}{3v}\right).$ (30)

Therefore, we will get the same $\phi$ and $U(\phi)$ as solution I:

$\displaystyle\phi(y)$ $\displaystyle=$ $\displaystyle v\,\text{tanh}(ky),$
(31) $\displaystyle U(\phi)$ $\displaystyle=$
$\displaystyle\frac{k^{2}}{v^{2}}({\phi^{2}}-{v^{2}})^{2}.$ (32)

But different from last model, we fix $W(\phi)=Q(\phi)$, which results in the
constant solution of the Lagrange multiplier

$\displaystyle\lambda(y)=-\frac{1}{2}\frac{W_{\phi}}{Q_{\phi}}=-\frac{1}{2}.$
(33)

Then we get a different form of the warp factor via Eq. (14)

$\displaystyle
A(y)=\frac{v^{2}}{18}\text{sech}^{2}(ky)+\frac{2v^{2}}{9}\ln\text{sech}(ky).$
(34)

The scalar potential $V$ has the following $\phi^{6}$ form

$\displaystyle
V(\phi)=\frac{k^{2}}{v^{2}}\left(\phi^{2}-v^{2}\right)^{2}-\frac{2k^{2}}{9v^{2}}\left({\phi^{2}}-3v^{2}\right)^{2}\phi^{2}.$
(35)

In this case, the energy density of the brane is

$\displaystyle\rho(y)$ $\displaystyle=$
$\displaystyle\frac{k^{2}v^{2}}{108}\Big{(}27+10v^{2}+3(9+2v^{2})\text{cosh}(2ky)\Big{)}$
(36) $\displaystyle\times$ $\displaystyle\text{sech}^{6}(ky).$

In this model, the vacuum expectation value $v$ of the scalar potential
$U(\phi)$ has an explicit effect on the warp factor and the energy density.
The shapes of the scalar field, the warp factor, and the energy density do not
change with non-vanishing $v$. Note that

$\displaystyle A(|y|\rightarrow\infty)$ $\displaystyle\rightarrow$
$\displaystyle\frac{2}{9}v^{2}\left(e^{-2k|y|}-k|y|\right)$ (37)
$\displaystyle\rightarrow$ $\displaystyle-\frac{2}{9}v^{2}k|y|.$

Therefore, the five-dimensional spacetime is also asymptotic AdS.

###  3.2 Period super-potential

####  3.2.1 Solution III

Next, we try to construct another form of brane solution by giving a period
super-potential, such as

$\displaystyle W(\phi)=3kn\sin^{\frac{1}{p}}\Big{(}\frac{\phi}{v}\Big{)}.$
(38)

At the same time, the warp factor $A(y)$ is assumed to be

$\displaystyle A(y)=\ln\text{sech}^{n}(ky).$ (39)

Then the scalar field can be solved from Eq. (14) as

$\displaystyle\phi(y)=v~{}\text{arcsin}[\tanh^{p}(ky)],$ (40)

which is a kink and a double kink for $p=1$ and $p=2n+1$ with positive integer
$n$, respectively. And from Eq. (15) we know that the super-potential $Q$ is

$\displaystyle
Q(\phi)=kp^{2}v^{2}\left[\mathcal{F}(p,\phi)+\mathcal{F}(-p,\phi)\right].$
(41)

where

$\displaystyle\mathcal{F}(p,\phi)$ $\displaystyle=$
$\displaystyle(-1)^{\frac{1}{2p}}\sec^{\frac{1}{p}}\Big{(}\frac{\phi}{v}\Big{)}~{}$
$\displaystyle\times$
$\displaystyle{{}_{2}F_{1}}\left(\frac{1}{2p},\frac{1}{2p};1+\frac{1}{2p};\sec^{2}\Big{(}\frac{\phi}{v}\Big{)}\right).$

Here ${{}_{2}F_{1}}$ is the hypergeometric function. The Lagrange multiplier
and the two potentials can also be solved:

$\displaystyle\lambda(y)$ $\displaystyle=$
$\displaystyle\frac{3n\cot^{2}\Big{(}\frac{\phi}{v}\Big{)}\sin^{\frac{2}{p}}\Big{(}\frac{\phi}{v}\Big{)}}{2p^{2}v^{2}\left[\sin^{\frac{2}{p}}\Big{(}\frac{\phi}{v}\Big{)}-1\right]},$
(42) $\displaystyle V(\phi)$ $\displaystyle=$ $\displaystyle
3k^{2}n\left[1-(1+2n)\sin^{\frac{2}{p}}\Big{(}\frac{\phi}{v}\Big{)}\right],$
(43) $\displaystyle U(\phi)$ $\displaystyle=$ $\displaystyle
k^{2}p^{2}v^{2}\text{sech}^{2}\Big{(}\frac{\phi}{v}\Big{)}$ (44)
$\displaystyle\times$
$\displaystyle\left[\sin^{1+\frac{1}{p}}\Big{(}\frac{\phi}{v}\Big{)}-\sin^{1-\frac{1}{p}}\Big{(}\frac{\phi}{v}\Big{)}\right]^{2}.$
(45)

The energy density in this case is given by

$\displaystyle\rho(y)=\frac{3}{2}k^{2}n(2n+1)\text{sech}^{2}(ky),$ (46)

which is the same as Eq. (29) for the first model. In fact, from the
definition of the energy density
$\rho=T_{MN}u^{M}u^{N}=(1/2)G_{MN}u^{M}u^{N}$, we know that the energy density
will have the same configuration for the same warp factor. Here, the warp
factors in this model and the first model have the same form, and hence so do
the energy densities.

####  3.2.2 Solution IV

Similarly, for the period super-potentials

$\displaystyle Q(\phi)=W(\phi)=kv^{2}\sin\Big{(}\frac{\phi}{v}\Big{)},$ (47)

we can obtain the following solution

$\displaystyle A(y)$ $\displaystyle=$
$\displaystyle\frac{v^{2}}{3}~{}\ln\big{[}\text{sech}(ky)\big{]},$ (48)
$\displaystyle\phi(y)$ $\displaystyle=$ $\displaystyle
v~{}\text{arctan}\big{[}\sinh(ky)\big{]},$ (49) $\displaystyle\lambda(y)$
$\displaystyle=$ $\displaystyle-\frac{1}{2},$ (50) $\displaystyle V(\phi)$
$\displaystyle=$ $\displaystyle
k^{2}v^{2}\left(\cos^{2}\left(\frac{\phi}{v}\right)-\frac{2}{3}v^{2}\sin^{2}\left(\frac{\phi}{v}\right)\right),$
(51) $\displaystyle U(\phi)$ $\displaystyle=$ $\displaystyle
k^{2}v^{2}\cos^{2}\frac{\phi}{v},$ (52) $\displaystyle\rho(y)$
$\displaystyle=$
$\displaystyle\frac{1}{6}k^{2}v^{2}(2v^{2}+3)\text{sech}^{2}(ky).$ (53)

####  3.2.3 Solution V

Next, we consider different period super-potentials $W(\phi)$ and $Q(\phi)$:

$\displaystyle W(\phi)$ $\displaystyle=$ $\displaystyle
3kn\tan\Big{(}\frac{\phi}{v}\Big{)},$ (54) $\displaystyle Q(\phi)$
$\displaystyle=$
$\displaystyle\frac{1}{2}kv^{2}\sin\Big{(}\frac{2\phi}{v}\Big{)},$ (55)

The warp factor and the scalar field will have the following explicit forms

$\displaystyle A(y)$ $\displaystyle=$
$\displaystyle\ln\big{[}\text{sech}^{n}(ky)\big{]},$ (56)
$\displaystyle\phi(y)$ $\displaystyle=$ $\displaystyle v\arctan[\tanh(ky)].$
(57)

The other functions are given by

$\displaystyle\lambda(y)\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!-\frac{3n}{2v^{2}}\sec[2\arctan\tanh(ky)],$ (58)
$\displaystyle\times$ $\displaystyle(1+\tanh^{2}(ky)),$ $\displaystyle
V(\phi)\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!3k^{2}n\sec^{2}\left(\frac{\phi}{v}\right)\left[(1+n)\cos\left(\frac{2\phi}{v}\right)-n\right],$
(59) $\displaystyle U(\phi)\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!k^{2}v^{2}\cos^{2}\Big{(}\frac{2\phi}{v}\Big{)}.$
(60)

###  3.3 Mixed super-potential

####  3.3.1 Solution VI

Finally, we would like to generate the brane model by giving the pair of
super-potentials $Q(\phi)$ and $W(\phi)$ as the mixed of polynomial and
period, for example

$\displaystyle
Q(\phi)=W(\phi)=kv\Big{[}\phi+v\sin\Big{(}\frac{\phi}{v}\Big{)}\Big{]},$ (61)

which results in the constant Lagrange multiplier

$\displaystyle\lambda(y)=-\frac{1}{2}.$ (62)

Then governed by Eq. (15), the scalar field $\phi$ is determined as

$\displaystyle\phi(y)=2v\arctan(ky),$ (63)

from which one can see that the asymptotic behavior of $\phi$ is limy→±∞
$\phi(y)=\pi$. The potentials $U$ and $V$ are

$\displaystyle U(\phi)\\!\\!\\!\\!\\!$ $\displaystyle=\\!\\!\\!\\!\\!$
$\displaystyle 4k^{2}v^{2}\cos^{4}\Big{(}\frac{\phi}{2v}\Big{)},$ (64)
$\displaystyle V(\phi)\\!\\!\\!\\!\\!$ $\displaystyle=\\!\\!\\!\\!\\!$
$\displaystyle
k^{2}v^{2}\left[4\cos^{4}\Big{(}\frac{\phi}{2v}\Big{)}-\frac{2}{3}\Big{(}\phi+v\sin\Big{(}\frac{\phi}{v}\Big{)}\Big{)}^{2}\right].$
(65)

The warp factor and the energy density read as

$\displaystyle A(y)\\!\\!\\!\\!$ $\displaystyle=\\!\\!\\!\\!$
$\displaystyle-\frac{2}{3}v^{2}ky\arctan(ky),$ (66)
$\displaystyle\rho(y)\\!\\!\\!\\!$ $\displaystyle=\\!\\!\\!\\!$
$\displaystyle\frac{1}{3}k^{2}v^{2}\Bigg{[}\frac{6}{\left(k^{2}y^{2}+1\right)^{2}}+\pi^{2}v^{2}$
(67) $\displaystyle-\\!\\!\\!\\!$ $\displaystyle
v^{2}\Big{(}2\arctan(ky)+\sin\left(2\arctan(ky)\right)\Big{)}^{2}\Bigg{]}.$

The asymptotic behavior of $A(y)$ is $A(|y|\rightarrow\pm\infty\rightarrow-\pi
v^{2}k|y|/3$ .

##  4 Tensor perturbation

In this section, we consider the linear tensor fluctuation of the metric
around the background. Following the previous research works in Refs. [33, 40,
41], we perform the following coordinate transformation

$\displaystyle dz=e^{-A(y)}dy,$ (68)

to get a conformally flat metric

$\displaystyle ds^{2}\\!\\!$ $\displaystyle=\\!\\!$ $\displaystyle
e^{2A(z)}(\eta_{\mu\nu}dx^{\mu}dx^{\nu}+dz^{2}).$ (69)

To simplify the fluctuations of the metric around the background, we only
consider the transverse and traceless part of the metric fluctuation, i.e., we
consider the following tensor perturbation of the metric:

$\displaystyle ds^{2}\\!\\!$ $\displaystyle=\\!\\!$
$\displaystyle\left(e^{2A(z)}\eta_{\mu\nu}+\hat{h}_{\mu\nu}(x,z)\right)dx^{\mu}dx^{\nu}+e^{2A(z)}dz^{2}$
(70) $\displaystyle=\\!\\!$ $\displaystyle
e^{2A(z)}\left[\left(\eta_{\mu\nu}+h_{\mu\nu}(x,z)\right)dx^{\mu}dx^{\nu}+dz^{2}\right].$

Here $h_{\mu\nu}$ is the tensor perturbation of the metric, and it satisfies
the transverse and traceless conditions [34]:
${h_{\mu}}^{\mu}=\partial^{\nu}h_{\mu\nu}=0$. The non-vanishing part of the
perturbation of Einstein tensor in Eq. (3) is the $\mu\nu$ components (since
$h_{55}=0$ and $\delta G_{55}=0$) and it reads

$\displaystyle\delta G_{\mu\nu}\\!\\!\\!\\!$ $\displaystyle=\\!\\!\\!\\!$
$\displaystyle-\frac{1}{2}\Box^{(4)}h_{\mu\nu}+(6{A^{\prime}}^{2}+3A^{\prime\prime})e^{2A}h_{\mu\nu}$
(71) $\displaystyle-$ $\displaystyle
2A^{\prime}e^{2A}h_{\mu\nu}^{\prime}-\frac{1}{2}e^{2A}h^{\prime\prime}_{\mu\nu},$

where the four-dimensional d’Alembertian is defined as
$\Box^{(4)}\equiv\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}$. Using Eq. (10),
we get the perturbation equation

$\displaystyle-\frac{1}{2}\Box^{(4)}h_{\mu\nu}-2A^{\prime}e^{2A}h_{\mu\nu}^{\prime}-\frac{1}{2}e^{2A}h^{\prime\prime}_{\mu\nu}=0.$
(72)

Considering Eq. (68), we rewrite the above equation under the coordinate $z$
as

$\displaystyle\Box^{(4)}h_{\mu\nu}+3\dot{A}\dot{h}_{\mu\nu}+\ddot{h}_{\mu\nu}=0,$
(73)

where the dot represents the derivative with respect to the coordinate $z$. By
performing the following decomposition

$\displaystyle
h_{\mu\nu}(x,z)=\varepsilon_{\mu\nu}e^{ikx}h(z)e^{-\frac{3}{2}A},$ (74)

Eq. (73) leads to a kind of Schrödinger equation

$\displaystyle\Big{(}-\partial_{z}^{2}+V_{\text{T}}(z)\Big{)}h(z)=m^{2}h(z),$
(75)

where $k^{2}=-m^{2}$ with $m$ the four-dimensional mass of a graviton KK mode,
and the effective potential given by

$\displaystyle
V_{\text{T}}(z)=\frac{3}{2}\ddot{A}(z)+\frac{9}{4}\dot{A}^{2}(z).$ (76)

The effective potential of the tensor perturbation under the physical
coordinate $y$ is

$\displaystyle
V_{\text{T}}(z(y))=e^{2A(y)}\left(\frac{3}{2}A^{\prime\prime}(y)+\frac{15}{4}A^{\prime
2}(y)\right).$ (77)

The zero mode of the tensor perturbation reads as

$\displaystyle h_{0}(z)=e^{{3A(z)}/{2}}\left(c_{1}+c_{2}\int
e^{-3A(z)}dz\right).$ (78)

This general form of the zero mode was first found in $f(R)$-brane model in
Ref. [48]. One can see that the effective potentials and the zero mode of the
tensor perturbation are only determined by the warp factor $A(y)$. It is easy
to show that the tensor zero mode (78) can be localized on the brane with the
choice of $c_{2}=0$ for all brane solutions in this paper (with $n>0$ for
solutions I, III and V and $v\neq 0$ for solutions II, IV and VI). The figures
of them for solutions I, III and V are the same, and they are similar for
solutions II, IV and VI. So without loss of generality, we plot two of them
for the brane solutions I and II in Figs. 1(c) and 1(f), respectively. These
potentials have a volcano-like shape. As the parameters $n$ and $v$ increase,
the potential wells in Figs. 1(c) and 1(f) become narrower and deeper,
respectively.

(a)

(b) (c) The effective potential $V_{\text{T}}$ and the non-normalized zero
mode of the tensor perturbation for brane models I, III and V. The parameter
is set as $n=1$ (red solid thick lines), $n=3$ (blue dashed lines), and $n=5$
(black solid thin lines).

(d)

(e) (f) The effective potential $V_{\text{T}}$ and the non-normalized zero
mode of the tensor perturbation for brane model II. The parameter is set as
$v=1$ (red solid thick lines), $v=1.5$ (blue dashed lines), and $v=2$ (black
solid thin lines).

It is easy to verify that the zero modes for the above brane models are
localized around the brane. So the four-dimensional Newtonian potential can be
realized on the brane. There is no tensor tachyon mode, thus the brane is
stable against the tensor perturbation.

##  5 Scalar perturbation

At last, we come to the scalar perturbation in this section. The perturbed
metric is given by

$\displaystyle
ds^{2}=e^{2A(z)}\left[(1+2\psi)\eta_{\mu\nu}dx^{\mu}dx^{\nu}+(1+2\Phi)dz^{2}\right].$
(79)

The scalar perturbation equations can be derived as

$\displaystyle e^{2A}$
$\displaystyle\\!\\!\\!\\!\left(\frac{1}{2}U_{\phi\phi}-U_{\phi}\frac{\ddot{\phi}}{\dot{\phi}^{2}}+2\frac{\dot{A}}{\dot{\phi}}U_{\phi}-\frac{4}{3}e^{-2A}\lambda\dot{\phi}^{2}\right)\delta\phi$
(80) $\displaystyle+$
$\displaystyle\\!\\!\\!\\!\left(\frac{1}{2}\frac{e^{2A}}{\dot{\phi}}U_{\phi}+\frac{\ddot{\phi}}{\dot{\phi}}-2\dot{A}\right)\dot{\delta\phi}-\ddot{\delta\phi}=0,$

and

$\displaystyle\Phi\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\frac{\dot{\delta\phi}}{\dot{\phi}}-\frac{e^{2A}}{2\dot{\phi}^{2}}U_{\phi}\delta\phi,$
(81) $\displaystyle\Phi\\!\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!-2\psi.$ (82)

It can be seen that there is only one degree of freedom for the scalar
perturbation. Finally, by using the background equations (3)-(5), we can
replace the potential $U$ and $\lambda$ in Eq. (80) with the functions $A$ and
$\phi$ and obtain the final form of the scalar perturbation $\delta\phi$:

$\displaystyle\Bigg{(}$
$\displaystyle\\!\\!\\!\\!3\dot{A}\frac{\ddot{\phi}}{\dot{\phi}}-3\ddot{A}-2\frac{\ddot{\phi}^{2}}{\dot{\phi}^{2}}+\frac{\dddot{\phi}}{\dot{\phi}}\Bigg{)}\delta\phi$
(83) $\displaystyle+$
$\displaystyle\\!\\!\\!\\!\left(2\frac{\ddot{\phi}}{\dot{\phi}}-3\dot{A}\right)\dot{\delta\phi}-\ddot{\delta\phi}=0.$

Redefining

$\displaystyle\delta\phi(x^{\mu},z)=\overline{\delta\phi}(x^{\mu})s(z)A^{-3/2}(z)\dot{\phi}(z),$
(84)

with the four-dimensional part of $\delta\phi$ satisfying
$\Box^{(4)}\overline{\delta\phi}(x^{\mu})=0$, we get the perturbation equation
for the scalar degree of freedom $s(z)$ from Eq. (83):

$\displaystyle\left(-\partial^{2}_{z}+V_{\text{S}}(z)\right)s(z)=0,$ (85)

where the effective potential is given by

$\displaystyle
V_{\text{S}}(z)=-\frac{3}{2}\ddot{A}(z)+\frac{9}{4}\dot{A}^{2}(z)$ (86)

in the $z$ coordinate or

$\displaystyle
V_{\text{S}}(z(y))=-e^{2A}\left(\frac{3}{2}A^{\prime\prime}(y)+\frac{3}{4}A^{\prime
2}(y)\right)$ (87)

in the $y$ coordinate. Note that the corresponding effective potential of (86)
in Ref. [40] is not right since there is an error in the matter equation (45)
in that paper (the term $+2\lambda(\partial_{z}\phi)^{2}$ should be
$-2\lambda(\partial_{z}\phi)^{2}$).

Comparing (86) with (76), one can see that the effective potential of the
scalar perturbation $V_{\text{S}}$ is dual to that of the tensor perturbation
$V_{\text{T}}$: $V_{\text{S}}=V_{\text{T}}~{}(A\rightarrow-A)$, and Eqs. (85)
and (75) can be rewritten as

$\displaystyle\mathcal{P}^{\dagger}\mathcal{P}s(z)\\!\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!\\!0,$ (88)
$\displaystyle\mathcal{P}\mathcal{P}^{\dagger}h(z)\\!\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!\\!0,$ (89)

where $\mathcal{P}=\partial_{z}+\frac{3}{2}\dot{A}$. The above two equations
ensure that both the scalar and tensor perturbations are stable. The zero mode
can be obtained by replacing $A\rightarrow-A$ from the solution of the tensor
zero mode (78):

$\displaystyle s(z)=e^{{-3A(z)}/{2}}\left(c_{1}+c_{2}\int e^{3A(z)}dz\right).$
(90)

This will lead to the conclusion that only one of the tensor and scalar zero
modes can be localized on the brane.

The effective potential (87) for solution I is given by

$\displaystyle
V_{\text{S}}=\frac{3}{8}k^{2}n\text{sech}^{2}(ky)\bigg{(}4+n-n\cosh(2ky)\bigg{)},$
(91)

which is shown in Fig. 2. From this figure, it can be seen that, the value of
the effective potential at $z=0$ will increase with the parameter $n$. This
can be checked from the expression $V_{\text{S}}(0)=\frac{3k^{2}n}{2}$. The
potential has two very shallow wells, and approaches $0^{-}$ when
$y\rightarrow\pm\infty$.

Figure 2: The effective potential $V_{\text{S}}(z(y))$ for solution I. The
parameter is set as $n=1$ (red solid thick lines), $n=3$ (blue dashed lines),
and $n=5$ (black solid thin lines).

For the brane solution I in (25) with $n>0$, it is easy to show that the zero
mode (90) of the scalar perturbation cannot be localized on the brane. Thus,
there is no additional fifth force coming from the scalar perturbation. For
other brane solutions, we have also the same conclusion.

##  6 Conclusion

In this work, we investigated the super-potential method with which the
second-order equations can be reduced to the first-order ones for thick brane
models in modified gravity with Lagrange multiplier. The main step of this
method is to introduce a pair of auxiliary super-potentials, i.e., $W(\phi)$
and $Q(\phi)$. With these two super-potentials, the field equations are
rewritten as Eqs. (13)-(17). Then we try to use the method to find a series of
analytical brane solutions via some polynomial super-potentials, period super-
potentials, and mixed super-potentials. The warp factor has the same shape and
the same asymptotic behavior at the boundary of the extra dimension for all
those solutions, and all of these branes are embedded in five-dimensional AdS
spacetime. The scalar field $\phi$ is a double kink for (40) with odd integer
$p\geq 3$ or a single kink for other solutions. These shapes of the scalar
field will affect the localization properties of fermions on the brane through
the Yukawa coupling $\eta\bar{\Psi}\phi\Psi$ [47].

We also considered the tensor and scalar perturbations of the brane system. It
was shown that both equations of motion of the perturbations can be
transformed into Schrödinger-like equations. Furthermore, these equations can
be recast as the forms of (88) and (89), which show that both the
perturbations are stable. The effective potential of the tensor perturbation
is dual to that of the scalar perturbation. Therefore, only one of the tensor
and scalar zero modes can be localized on the brane. For all of our brane
solutions, the tensor zero mode can be localized on the brane while the scalar
zero mode cannot. Thus, the four-dimensional Newtonian potential can be
recovered on the brane and there is no additional fifth force contradicting
with the experiments.

## Acknowledgement

This work was supported by the National Natural Science Foundation of China
(Grant No. 11875151, No. 12047501, and No. 11705070). Yi Zhong was supported
by the Fundamental Research Funds for the Central Universities (Grant No.
531107051196).

## References

  * [1] A. H. Chamseddine and V. Mukhanov, _Mimetic Dark Matter_ , _JHEP_ 11 (2013)135, [arXiv:1308.5410].
  * [2] A. Golovnev, _On the recently proposed Mimetic Dark Matter_ , _Phys. Lett._ B 728 (2014) 39-40, [arXiv:1310.2790].
  * [3] A. H. Chamseddine, V. Mukhanov and A. Vikman, _Cosmology with Mimetic Matter_ , _JCAP_ 1406 (2014) 017, [arXiv:1403.3961].
  * [4] A. Casalino, M. Rinaldi, L. Sebastiani, and S. Vagnozzi, _Mimicking dark matter and dark energy in a mimetic model compatible with GW170817_ , _Phys. Dark Univ._ 22 (2018) 108, [arXiv:1803.02620].
  * [5] J. Matsumoto, _Unified description of dark energy and dark matter in mimetic matter model_ , [arXiv:1610.07847].
  * [6] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, _Unimodular-Mimetic Cosmology_ , _Class. Quant. Grav._ 33 (2016) 125017, [arXiv:1601.07057].
  * [7] S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, _Viable Mimetic Completion of Unified Inflation-Dark Energy Evolution in Modified Gravity_ , _Phys. Rev._ D 94 (2016) 10, 104050  [arXiv:1608.07806].
  * [8] S. Brahma, A. Golovnev, and D. H. Yeoma, _On singularity-resolution in mimetic gravity_ , _Phys. Lett._ B 782 (2018) 280, [arXiv:1803.03955].
  * [9] J. Dutta, W. Khyllep, E. N. Saridakis, N. Tamanini, and S. Vagnozzi, _Cosmological dynamics of mimetic gravity_ , _JCAP_ 1802 (2018) 041, [arXiv:1711.07290].
  * [10] S. A. H. Mansoori, A. Talebian, and H. Firouzjahi, _Mimetic Inflation_ , [arXiv:2010.13495].
  * [11] F. Izaurieta, P. Medina, N. Merino, P. Salgado, and O. Valdivia, _Mimetic Einstein-Cartan-Kibble-Sciama (ECKS) gravity_ , [arXiv:2007.07226].
  * [12] L. Sebastiani, S. Vagnozzi and R. Myrzakulov, _Mimetic gravity: a review of recent developments and applications to cosmology and astrophysics_ , _Adv. High Energy Phys._ 2017 (2017) 3156915, [arXiv:1612.08661].
  * [13] K. Nozari, and N. Sadeghnezhad, _Braneworld mimetic f(R) gravity_ , _Int. J. Geom. Meth. Mod. Phys._ 16 (2019) 03, 1950042, [ arXiv:1612.08661].
  * [14] D. Momeni, R. Myrzakulov and E. Güdekli, _Cosmological viable mimetic $f(R)$ and $f(R,T)$ theories via Noether symmetry_, _Int. J. Geom. Meth. Mod. Phys._ 12 (2015) 1550101, [arXiv:1502.00977].
  * [15] S. Nojiri and S. D. Odintsov, _Mimetic f(R) gravity: Inflation, dark energy and bounce_ , _Modern Physics Letters A_ 29 (dec, 2014) 1450211, [arXiv:1408.3561].
  * [16] S. D. Odintsov and V. K. Oikonomou, _Mimetic $F(R)$ inflation confronted with Planck and BICEP2/Keck Array data_, _Astrophys. Space Sci._ 361 (2016) 174, [arXiv:1512.09275].
  * [17] S. D. Odintsov and V. K. Oikonomou, _Viable Mimetic $F(R)$ Gravity Compatible with Planck Observations_, _Annals Phys._ 363 (2015) 503–514, [arXiv:1508.07488].
  * [18] S. D. Odintsov and V. K. Oikonomou, _Accelerating cosmologies and the phase structure of F(R) gravity with Lagrange multiplier constraints: A mimetic approach_ , _Phys. Rev._ D 93 (2016) 023517, [arXiv:1511.04559].
  * [19] V. K. Oikonomou, _Reissner-Nordström Anti-de Sitter Black Holes in Mimetic $F(R)$ Gravity_, _Universe_ 2 (2016) 10, [arXiv:1511.09117].
  * [20] F. Arroja, T. Okumura, N. Bartolo, P. Karmakar, and S. Matarrese, _Large-scale structure in mimetic Horndeski gravity_ , _JCAP_ 05(2018)050, [arXiv:1708.01850].
  * [21] F. Arroja, N. Bartolo, P. Karmakar, and S. Matarrese, _Cosmological perturbations in mimetic Horndeski gravity_ , _JCAP_ 04 (2016) 042, [arXiv:1512.09374].
  * [22] G. Cognola, R. Myrzakulov, L. Sebastiani, S. Vagnozzi and S. Zerbini, _Class. Quant. Grav._ 33 (2016) 225014 . arXiv:1601.00102.
  * [23] S. Capozziello, A. N. Makarenko, and S. D. Odintsov, _Gauss-Bonnet dark energy by Lagrange multipliers_ , _Phys. Rev._ D 87 (2013) 084037, [arXiv:1302.0093].
  * [24] A. V. Astashenok, S. D. Odintsov and V. K. Oikonomou, _Modified Gauss-Bonnet gravity with the Lagrange multiplier constraint as mimetic theory_ , _Class. Quant. Grav._ 32 (2015) 185007, [arXiv:1504.04861].
  * [25] L. Randall and R. Sundrum, _A Large mass hierarchy from a small extra dimension_ , _Phys. Rev. Lett._ 83 (1999) 3370–3373, [arXiv:hep-ph/9905221].
  * [26] L. Randall and R. Sundrum, _An Alternative to compactification_ , _Phys. Rev. Lett._ 83 (1999) 4690–4693, [arXiv:hep-th/9906064].
  * [27] H. Davoudiasl, J. L. Hewett and T. G. Rizzo, _Bulk gauge fields in the Randall-Sundrum model_ , _Phys. Lett._ B 473 (2000) 43–49, [arXiv:hep-ph/9911262].
  * [28] T. Shiromizu, K. i. Maeda and M. Sasaki, _The Einstein equation on the 3-brane world_ , _Phys. Rev._ D 62 (2000) 024012, [arXiv:gr-qc/9910076].
  * [29] T. Gherghetta and A. Pomarol, _Bulk fields and supersymmetry in a slice of AdS_ , _Nucl. Phys._ B 586 (2000) 141–162, [arXiv:hep-ph/0003129].
  * [30] T. G. Rizzo, _Introduction to Extra Dimensions_ , _AIP Conf. Proc._ 1256 (2010) 27–50, [arXiv:1003.1698].
  * [31] K. Yang, Y.-X. Liu, Y. Zhong, X.-L. Du and S.-W. Wei, _Gravity localization and mass hierarchy in scalar-tensor branes_ , _Phys. Rev._ D 86 (2012) 127502, [arXiv:1212.2735].
  * [32] K. Agashe, A. Azatov, Y. Cui, L. Randall and M. Son, _Warped Dipole Completed, with a Tower of Higgs Bosons_ , _JHEP_ 06 (2015) 196, [arXiv:1412.6468].
  * [33] C. Csaki, J. Erlich, T. J. Hollowood and Y. Shirman, _Universal aspects of gravity localized on thick branes_ , _Nucl. Phys._ B 581 (2000) 309–338, [arXiv:hep-th/0001033].
  * [34] O. DeWolfe, D. Z. Freedman, S. S. Gubser and A. Karch, _Modeling the fifth-dimension with scalars and gravity_ , _Phys. Rev._ D 62 (2000) 046008, [arXiv:hep-th/9909134].
  * [35] M. Gremm, _Four-dimensional gravity on a thick domain wall_ , _Phys. Lett._ B 478 (2000) 434–438, [arXiv:hep-th/9912060].
  * [36] Y.-X. Liu, _Introduction to Extra Dimensions and Thick Braneworlds_ , arXiv:1707.08541.
  * [37] N. Sadeghnezhad and K. Nozari, _Braneworld Mimetic Cosmology_ , _Phys. Lett._ B 769 (2017) 134–140, [arXiv:1703.06269].
  * [38] Y. Zhong, Y. Zhong, Y.-P. Zhang, W.-D. Guo, Y.-X. Liu, _Gravitational resonances in mimetic thick branes_ , _JHEP_ 04 (2019) 154, [arXiv:1812.06453].
  * [39] Q. Xiang, Y. Zhong, Q.-Y. Xie, L. zhao, _Flat and bent branes with inner structure in two-field mimetic gravity_ , [arXiv:2011.10266].
  * [40] Y. Zhong, Y. Zhong, Y.-P. Zhang, Y.-X. Liu, _Thick branes with inner structure in mimetic gravity_ , _Eur. Phys. J._ C 78 (2018) 45, [arXiv:1711.09413].
  * [41] D. Bazeia, D. A. Ferreira, D. C. Moreira, _First order formalism for thick branes in modified gravity with Lagrange multiplier_ , _EPL_ 129 (2020) 11004, [arXiv:2002.00229].
  * [42] V.I. Afonso, D. Bazeia, and L. Losano, _First-order formalism for bent brane_ , _Phys.Lett._ B 634 (2006) 526-530, [arXiv:hep-th/0601069].
  * [43] B. Janssen, P. Smyth, T. V. Riet, and B. Vercnocke, _A first-order formalism for timelike and spacelike brane solutions_ , _JHEP_ 007 (2008) 0804, [arXiv:0712.2808].
  * [44] R. Menezes, _First Order Formalism for Thick Branes in Modified Teleparallel Gravity_ , _Phys. Rev._ D 89 (2014) 125007, [arXiv:1403.5587].
  * [45] D. Bazeia, D. A. Ferreira, F. S. N. Lobo, and J. L. Rosa, _Novel modified gravity braneworld configurations with a Lagrange multiplier_ , arXiv:2011.06240.
  * [46] D. Bazeia and A. R. Gomes, _Bloch brane_ , _JHEP_ 05 (2004) 012, [arXiv:hep-th/0403141].
  * [47] Y.-X. Liu, J. Yang, Z.-H. Zhao, C.-E. Fu, and Y.-S. Duan, _Fermion Localization and Resonances on A de Sitter Thick Brane_ , _Phys. Rev._ D 80 (2009) 065019, [arXiv:0904.1785].
  * [48] Z.-Q. Cui, Z.-C. Lin, J.-J. Wan, Y.-X. Liu, and L. Zhao, _Tensor Perturbations and Thick Branes in Higher-dimensional $f(R)$ Gravity_, _JHEP_ 12 (2020) 130, [arXiv:2009.00512].