# Weyl Starobinsky inflation

Antonio De Felice1, Ryodai Kawaguchi2, Kotaro Mizui2, and Shinji Tsujikawa2
1Center for Gravitational Physics and Quantum Information, Yukawa Institute
for Theoretical Physics, Kyoto University, 606-8502, Kyoto, Japan
2Department of Physics, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo
169-8555, Japan

###### Abstract

In Starobinsky inflation with a Weyl squared Lagrangian $-\alpha C^{2}$, where
$\alpha$ is a coupling constant, we study the linear stability of cosmological
perturbations on a spatially flat Friedmann-Lemaître-Robertson-Walker
background. In this theory, there are two dynamical vector modes propagating
as ghosts for $\alpha>0$, whose condition is required to avoid tachyonic
instabilities of vector perturbations during inflation. The tensor sector has
four propagating degrees of freedom, among which two of them correspond to
ghost modes. However, tensor perturbations approach constants after the Hubble
radius crossing during inflation, and hence the classical instabilities are
absent. In the scalar sector, the Weyl curvature gives rise to a ghost mode
coupled to the scalaron arising from the squared Ricci scalar. We show that
two gauge-invariant gravitational potentials, which are both dynamical in our
theory, are subject to exponential growth after the Hubble radius crossing.
There are particular gauge-invariant combinations like the curvature
perturbations whose growth is suppressed, but it is not possible to remove the
instability of other propagating degrees of freedom present in the perturbed
metric. This violent and purely classical instability present in the scalar
sector makes the background unviable. Furthermore, the presence of such
classical instability makes the quantization of the modes irrelevant, and the
homogeneous inflationary background is spoiled by the Weyl curvature term.

††preprint: YITP-23-106, WUCG-23-09

##  I Introduction

The inflationary paradigm [1, 2, 3, 4] can successfully resolve several
shortcomings in Big Bang cosmology, e.g., the horizon and flatness problems.
Moreover, it can explain the origin of large-scale structures in the Universe
by stretching quantum fluctuations over super-Hubble scales during the
accelerated expansion [5, 6, 7, 8, 9]. The spectra of scalar and tensor
perturbations predicted in standard slow-roll inflation are consistent with
the observed cosmic microwave background (CMB) temperature anisotropies. After
the data release of WMAP [10] and Planck [11] satellites, we have been able to
distinguish between many different models of inflation. In particular, the
first model advocated by Starobinsky [5] is still one of the best-fit models
to the Planck CMB data [12].

In the Starobinsky model, inflation is driven by a quadratic Ricci scalar term
$\beta R^{2}$, where $\beta$ is a positive coupling constant. The period of
cosmic acceleration ends when the $\beta R^{2}$ term drops below the Ricci
scalar $R$ [13, 14, 15, 16, 17]. The quadratic curvature scalar gives rise to
a new propagating degree of freedom (DOF) dubbed “scalaron” [5] with a mass
squared $m_{S}^{2}=1/(6\beta)$ on the Minkowski background [18, 19, 20, 21,
22]. Indeed, $f(R)$ gravity given by the Lagrangian $f(R)=R+\beta R^{2}$ is
equivalent to Brans-Dicke theory [23] with a scalaron potential arising from
the gravitational sector [24, 25, 26]. From the observed amplitude of CMB
temperature fluctuations, the scalaron mass is constrained to be $m_{s}\simeq
10^{-5}M_{\rm pl}$, where $M_{\rm pl}$ is the reduced Planck mass [15, 27, 28,
29]. The Starobinsky model predicts the scalar spectral index $n_{s}\simeq
1-2/N$ and the tensor-to-scalar ratio $r\simeq 12/N^{2}$ [30, 29, 31], where
$N$ is the number of e-foldings counted backward from the end of inflation. On
scales relevant to the CMB observations ($N\simeq 55\sim 60$), the theoretical
predictions of $n_{s}$ and $r$ are well consistent with the Planck data
combined with other data [12].

From the viewpoint of an ultraviolet completion of gravity, there are also
other quadratic curvature contributions to the Lagrangian constructed from
scalar products of the Riemann tensor $R_{\mu\nu\rho\sigma}$ and the Ricci
tensor $R_{\mu\nu}$ [32]. Given that the Gauss-Bonnet curvature invariant
${\cal
G}=R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$ is a
topological term that does not affect the field equations of motion [33], the
general quadratic-order Lagrangian can be expressed in the form $L_{2}=-\alpha
C^{2}+\beta R^{2}$, where $C^{2}=2R_{\mu\nu}R^{\mu\nu}-2R^{2}/3+{\cal G}$ is a
squared Weyl curvature. This quadratic theory of gravity, which was originally
advocated by Stelle [34], is renormalizable and also asymptotically free [35].
However, the Weyl curvature generally gives rise to ghost DOFs associated with
derivative terms higher than second order in the field equations of motion
[36].

Albeit the appearance of ghosts in Weyl gravity with the Lagrangian $-\alpha
C^{2}$, the perturbative expansion about the Minkowski vacuum shows that all
the linear perturbations in scalar, vector, and tensor sectors propagate with
the speed of light [37, 38, 39]. This means that, in the absence of additional
matter sources, the perturbations are not subject to Laplacian instabilities.
Hence, on Minkowski, the Weyl ghosts can be of “soft” types [40], i.e., even
in the presence of the ghosts, the classical perturbations do not grow on the
given background. However, this situation should be different by introducing
some matter fields or by taking into account the $\beta R^{2}$ term on curved
backgrounds. The latter corresponds to Stelle’s quadratic curvature theory
mentioned above, in which case the scalaron field arising from the $\beta
R^{2}$ term can be gravitationally coupled to the Weyl ghost.

In this paper, we would like to address the stability of linear perturbations
on the inflationary background realized by Stelle’s theory. We note that there
are some related works in which the dynamics of cosmological perturbations
during inflation were discussed in the presence of the Weyl curvature term
[41, 42, 43, 44, 45]. Most of those papers assumed the existence of a
canonical scalar field besides the Weyl and Einstein-Hilbert terms. Since the
squared Weyl curvature is conformally invariant, the Lagrangian $L=R-\alpha
C^{2}+\beta R^{2}$ of Stelle’s theory can be transformed to that in the
Einstein frame with kinetic and potential energies of the scalaron field as
well as the Weyl term [46, 47, 48, 49].

The analysis of tensor perturbations on a spatially-flat Friedmann-Lemaître-
Robertson-Walker (FLRW) background [41] showed that the tensor ghosts can be
soft during inflation in that they do not grow by linear instabilities. In the
vector sector, there are two dynamical propagating DOFs arising from the Weyl
curvature term $-\alpha C^{2}$ [42]. The ghosts do not appear if $\alpha<0$,
but vector perturbations are subject to tachyonic instabilities. For
$\alpha>0$ the ghosts are present, but vector perturbations decay during
inflation. Thus, despite the presence of ghosts for $\alpha>0$, both tensor
and vector perturbations are not prone to classical Laplacian or tachyonic
instabilities. These results were already recognized in Refs. [41, 42]
according to the analysis in the Einstein frame, but we will study whether a
similar property holds in the Jordan frame. Indeed, for $\alpha>0$, the
classical instabilities are absent for both the tensor and vector sectors.

In the scalar sector, the analysis of Ref. [42] in the Einstein frame of
Stelle’s theory showed that gravitational potentials in a Newtonian gauge
exhibit rapid growth after the Hubble radius during inflation. On the other
hand, the same paper also found that the curvature perturbation in a comoving
gauge remains bounded. In Ref. [44], it was claimed that the former growth of
gravitational potentials is a gauge artifact and that scalar perturbations are
not prone to real instabilities. So far, it is not yet clear whether the
instability in the scalar sector induced by the Weyl term corresponds to a
real, physical one. To clarify this issue, we need to scrutinize whether the
instability of scalar perturbations generally persists or not independent of
the gauge choices.

In this paper, we will study the evolution of cosmological perturbations
during inflation in the Jordan frame of Stelle’s theory by paying particular
attention to the linear stability of the scalar modes. For this purpose, we
choose several different physical gauges and analytically derive the closed
fourth-order differential equations for the gravitational potentials as well
as other gauge-invariant variables like curvature perturbations. We will
explicitly show that two dynamical propagating DOFs arise from the $\beta
R^{2}$ term (i.e., scalaron) and the Weyl curvature, one of which always
behaves as a ghost mode. Therefore, in general, there are four independent
initial conditions necessary to uniquely specify the classical evolution of
the scalar sector, and that is the reason why the system can be described in
terms of a closed fourth-order differential equation for one single scalar
mode, or evidently, by two second-order differential equations for two
independent fields.

We will show that the two gravitational potentials $\Psi$ and $\Phi$, which
are both propagating DOFs, exponentially grow after the Hubble radius
crossing. This instability of $\Psi$ and $\Phi$ occurs independently of the
gauge choice made to study their dynamics. Among other relevant gauge-
invariant variables, we also find that the curvature perturbation is a
specific variable approaching a constant in the large-scale limit. However,
the exponential increase of at least one dynamical scalar DOF appearing in the
perturbed line element does not allow the FLRW spacetime to be a stable
cosmological background. Thus, the inflationary background is violated by this
real, physical, and classical instability of scalar perturbations induced by
the Weyl ghost coupled to the scalaron. Therefore, we conclude the propagating
ghost DOF in the scalar sector is not of the soft type. Well before the end of
inflation, the cosmological background is spoiled by the classical instability
and it changes to a highly inhomogeneous Universe. In such a context, we
believe that this lack of a homogeneous background makes the quantization of
perturbations irrelevant. If the ghost modes were soft, then the quantization
procedure would acquire relevance and the results could be interesting.
However, this is not the case for inflation in quadratic gravity with the Weyl
term.

This paper is organized as follows. In Sec. V, we briefly review the
background dynamics of inflation realized in Stelle’s theory. In Sec. III, we
revisit how vector perturbations propagate as truly dynamical DOFs and find
that two ghost modes propagate under the condition $\alpha>0$ for the absence
of tachyonic instabilities. In Sec. IV, we show that, despite the appearance
of ghosts arising from the Weyl term, the four dynamical DOFs of tensor
perturbations approach constants after the Hubble radius crossing. In Sec. V,
we study the evolution of scalar perturbations by choosing several different
gauges and show that, independently of the gauge choices, two modes propagate
and the FLRW background is spoiled by the presence of instabilities of at
least one of the dynamical DOF present in the perturbed line element. Although
the classical instability itself is present for any non-zero initial
conditions of scalar modes, we confirm in Sec. VI its presence by numerically
integrating the perturbation equations of motion with proper initial
conditions. Sec. VII is devoted to conclusions.

##  II Inflation in quadratic gravity

The action in quadratic gravity contains scalar products of two contractions
of the Riemann tensor $R_{\mu\nu\rho\sigma}$, Ricci tensor $R_{\mu\nu}$, and
Ricci scalar $R$. On using the property that the Gauss-Bonnet term ${\cal
G}=R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$ is
topological, the Riemann products $R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$
can be eliminated from the action. Taking the Einstein-Hilbert term $M_{\rm
Pl}^{2}R/2$ into account, the action of quadratic gravity can be expressed in
the form [34]

${\cal S}=\frac{M_{\rm pl}^{2}}{2}\int{\rm d}^{4}x\sqrt{-g}\left(R-\alpha
C^{2}+\beta R^{2}\right)\,,$ (1)

where $g$ is a determinant of the metric tensor $g_{\mu\nu}$, $\alpha$ and
$\beta$ are constants, and $C^{2}$ is the Weyl tensor squared given by

$C^{2}=2R_{\mu\nu}R^{\mu\nu}-\frac{2}{3}R^{2}+{\cal G}\,.$ (2)

Up to boundary terms, the action (1) can be expressed as

${\cal S}=\frac{M_{\rm pl}^{2}}{2}\int{\rm d}^{4}x\sqrt{-g}\left[R-2\alpha
R_{\mu\nu}R^{\mu\nu}+\left(\frac{2}{3}\alpha+\beta\right)R^{2}\right]\,.$ (3)

We consider a spatially flat FLRW background described by the line element

${\rm d}s^{2}=-N^{2}(t){\rm d}t^{2}+a^{2}(t)\delta_{ij}{\rm d}x^{i}{\rm
d}x^{j}\,,$ (4)

where $a(t)$ is a time-dependent scale factor, and $N(t)$ is a lapse function.
Varying the action (3) with respect to $N(t)$ and $a(t)$, respectively, and
setting $N(t)=1$ at the end, it follows that

$\displaystyle
H^{2}+6\beta\left(6H^{2}\dot{H}-\dot{H}^{2}+2H\ddot{H}\right)=0\,,$ (5)
$\displaystyle
2\dot{H}+3H^{2}+6\beta\left(18H^{2}\dot{H}+9\dot{H}^{2}+12H\ddot{H}+2\dddot{H}\right)=0\,,$
(6)

where $H=\dot{a}/a$ is the Hubble expansion rate, with a dot being the
derivative with respect to $t$. At the background level, the Weyl curvature
term does not contribute to the field equations of motion. This is an outcome
of the conformal invariance of the Weyl curvature term, whose components
vanish for the conformally flat background. Taking the time derivative of (5)
and combining it with Eq. (5), we obtain the same equation as (6). This means
that there is only a single independent equation of motion, Eq. (5), governing
the background dynamics.

During inflation, the Hubble expansion rate is nearly constant, and hence the
last two terms in the parenthesis of Eq. (5) are suppressed relative to the
term $6H^{2}\dot{H}$. Then, there is the approximate relation

$\dot{H}\simeq-\frac{1}{36\beta}=-\frac{m_{S}^{2}}{6}\,,$ (7)

where $m_{S}^{2}$ corresponds to the mass squared of a scalaron field given by
[18, 19, 20, 21, 22]

$m_{S}^{2}=\frac{1}{6\beta}\,.$ (8)

Provided that

$\beta>0\,,$ (9)

there is no tachyonic instability arising from the negative value of
$m_{S}^{2}$. Under the condition (9), the Hubble parameter (7) also decreases
during inflation. We will impose the condition (9) throughout the discussion
below.

From Eq. (7), we obtain the following integrated solutions:

$\displaystyle H(t)$ $\displaystyle\simeq$ $\displaystyle
H_{i}-\frac{m_{S}^{2}}{6}(t-t_{i})\,,$ (10) $\displaystyle a(t)$
$\displaystyle\simeq$ $\displaystyle
a_{i}\exp\left[H_{i}(t-t_{i})-\frac{m_{S}^{2}}{12}(t-t_{i})^{2}\right]\,,$
(11)

where $H_{i}$ and $a_{i}$ are the values of $H$ and $a$ at the onset of
inflation, respectively. We introduce the slow-roll parameter $\epsilon$, as

$\epsilon\equiv-\frac{\dot{H}}{H^{2}}\simeq\frac{m_{S}^{2}}{6H^{2}}\,.$ (12)

The end of inflation is characterized by the Hubble parameter $H_{f}$ when
$\epsilon$ becomes equivalent to 1 so that $H_{f}\simeq m_{S}/\sqrt{6}$. As we
will see below, $H_{i}$ is larger than $m_{S}$. Then, by using Eq. (10), we
can approximately estimate the time $t_{f}$ at the end of inflation, as
$t_{f}\simeq t_{i}+6H_{i}/m_{S}^{2}$. The number of e-foldings acquired during
inflation is given by

$N\equiv\int_{t_{i}}^{t_{f}}H\,{\rm d}t\simeq
H_{i}(t_{f}-t_{i})-\frac{m_{S}^{2}}{12}(t_{f}-t_{i})^{2}\simeq\frac{3H_{i}^{2}}{m_{S}^{2}}\,.$
(13)

Taking $N=60$ as a typical minimal e-folding number required to address the
horizon and flatness problems, we obtain the initial Hubble parameter
$H_{i}\simeq 4.5m_{S}$. This value translates to

$\beta H_{i}^{2}\simeq 3.4\,,$ (14)

and hence $\beta H^{2}$ is of order 1 during inflation. From the viewpoint of
ultraviolet completion of gravity, it is natural to consider the value of
$|\alpha|$ same order as $\beta$.

We note that quadratic gravity given by the action (1) corresponds to the
$f(R)=R+\beta R^{2}$ gravity with the Weyl squared term $-\alpha C^{2}$. Under
a conformal transformation of the metric tensor $g_{\mu\nu}$, the theory can
be transformed to a metric frame described by Einstein’s gravity in the
presence of a scalaron field with the potential and the Weyl squared term [46,
47, 48, 49]. In this Einstein frame, the dynamics of cosmological
perturbations during inflation were carried out in Refs. [41, 42, 43, 44, 45]
without necessarily relating the scalaron potential with the one arising from
the original Lagrangian $f(R)=R+\beta R^{2}$. In this paper, we will perform
all the analysis in the physical Jordan frame. To study the dynamics of
perturbations during inflation, we do not need to take into account additional
matter sources to the Jordan-frame action (1).

Around the background (4) with $N(t)=1$, we consider metric perturbations
which depend on the cosmic time $t$ and spatial coordinates $x^{i}$. The
perturbed line element is given by

${\rm d}s^{2}=-\left(1+2A\right){\rm
d}t^{2}+2a(t)\left(\partial_{i}B+V_{i}\right){\rm d}t\,{\rm
d}x^{i}+a^{2}(t)\left[(1+2\psi)\delta_{ij}+2\partial_{i}\partial_{j}E+2\partial_{j}F_{i}+h_{ij}\right]{\rm
d}x^{i}{\rm d}x^{j}\,,$ (15)

where we used the notation $\partial_{i}=\partial/\partial x^{i}$, and the
Latin indices represent spatial coordinates. The four quantities $A$, $B$,
$\psi$, $E$ correspond to scalar perturbations, while $V_{i}$ and $F_{i}$ are
vector perturbations satisfying the divergence-free conditions
$\delta^{ij}\partial_{j}V_{i}=0$ and $\delta^{ij}\partial_{j}F_{i}=0$. The
intrinsic tensor perturbation is given by $h_{ij}$, which satisfies the
traceless and transverse conditions ${h_{i}}^{i}=0$ and
$\delta^{ik}\partial_{k}h_{ij}=0$.

##  III Vector perturbations

In quadratic gravity given by the action (3), we first study the dynamics of
vector perturbations during inflation. Since there is the residual gauge DOF,
we choose the following gauge condition

$F_{i}=0\,,$ (16)

where $i=1,2,3$. Then, the perturbed line element in the vector sector is
given by

${\rm d}s^{2}=-{\rm d}t^{2}+2a(t)V_{i}\,{\rm d}t{\rm
d}x^{i}+a^{2}(t)\delta_{ij}{\rm d}x^{i}{\rm d}x^{j}\,.$ (17)

For practical computations, it is convenient to choose the vector-field
configuration in the form

$V_{i}=[V_{1}(t,z),V_{2}(t,z),0]\,,$ (18)

which satisfies the divergence-free condition $\partial^{i}V_{i}=0$. Expanding
the action (3) up to quadratic order in vector perturbations, integrating it
by parts, and using the background Eq. (5), the second-order action yields

${\cal S}_{v}^{(2)}=-\frac{M_{\rm Pl}^{2}}{2}\alpha\sum_{i=1,2}\int{\rm
d}^{4}x\,a\left\\{\dot{U}_{i}^{2}-\frac{U_{i}^{\prime
2}}{a^{2}}-\left[\frac{1}{2\alpha}+\frac{6\beta}{\alpha}(2H^{2}+\dot{H})\right]U_{i}^{2}\right\\}\,,$
(19)

where a prime represents the derivative with respect to $z$, and $U_{i}\equiv
V_{i}^{\prime}$. The Weyl curvature term gives rise to two dynamical vector
perturbations $U_{1}$ and $U_{2}$. For $\alpha>0$, it is clear from the action
(19) that both $U_{1}$ and $U_{2}$ behave as ghosts. On the other hand, the
vector ghosts are absent if $\alpha<0$.

Varying the action (19) with respect to $U_{i}$ (with $i=1,2$), we obtain
their equations of motion in real space. Then, we perform the Fourier
transformation

$U_{i}=\frac{1}{(2\pi)^{1/2}}\int{\rm d}k\,\tilde{U}_{i}(t,k)e^{ikz}\,,$ (20)

where $k$ is a comoving wavenumber, and $\tilde{U}_{i}$ is a function of $t$
and $k$. Omitting the tilde from $\tilde{U}_{i}(t,k)$ in the following, we
obtain the vector perturbation equations of motion in Fourier space, as

$\ddot{U}_{i}+H\dot{U}_{i}+\left[\frac{k^{2}}{a^{2}}+m_{W}^{2}+\frac{6\beta}{\alpha}(2H^{2}+\dot{H})\right]U_{i}=0\,,$
(21)

where

$m_{W}^{2}\equiv\frac{1}{2\alpha}\,,$ (22)

is the mass squared arising from the Weyl curvature [36, 37, 38, 39]. We note
that the Lagrangian $\beta R^{2}$ contributes to the vector mass through the
term $(6\beta/\alpha)(2H^{2}+\dot{H})$. During inflation ($|\dot{H}|\ll
H^{2}$), the effective mass squared of vector perturbations is approximately
given by

$m_{\rm eff}^{2}\simeq m_{W}^{2}+\frac{12\beta}{\alpha}H^{2}\,.$ (23)

Since we are considering a positive coupling $\beta$, we have $m_{\rm
eff}^{2}>0$ if $\alpha>0$ and $m_{\rm eff}^{2}<0$ if $\alpha<0$.

For the Weyl coupling constant in the range $|\alpha|\lesssim\beta$, the
effective mass squared (23) is at least of order $H^{2}$. Then, after a
squared physical wavenumber $(k/a)^{2}$ drops below $|m_{\rm eff}^{2}|$ during
inflation, Eq. (21) approximately reduces to

$\ddot{U}_{i}+H\dot{U}_{i}+m_{\rm eff}^{2}\,U_{i}\simeq 0\,.$ (24)

Since $H$ and $m_{\rm eff}^{2}$ can be approximated as constants during
inflation, we obtain the following solution

$U_{i}=c_{1}\exp\left(\frac{-H+\sqrt{H^{2}-4m_{\rm
eff}^{2}}}{2}t\right)+c_{2}\exp\left(\frac{-H-\sqrt{H^{2}-4m_{\rm
eff}^{2}}}{2}t\right)\,,$ (25)

where $c_{1}$ and $c_{2}$ are integration constants.

For $\alpha>0$ (i.e., $m_{\rm eff}^{2}>0$), if $H$ is initially in the range
$H>2m_{\rm eff}$, the amplitude of vector perturbations first decreases in
proportion to $U_{i}\propto\exp[(-H+\sqrt{H^{2}-4m_{\rm eff}^{2}})t/2]$. After
$H$ decreases below $2m_{\rm eff}$, the vector perturbation starts to
oscillate with the decreasing amplitude ($|U_{i}|\propto e^{-Ht/2}$). This
means that, even though the two vector ghosts are present for $\alpha>0$,
vector perturbations decay exponentially during inflation.

For $\alpha<0$, the negative mass squared $m_{\rm eff}^{2}$ leads to the
tachyonic growth of $U_{i}$. During inflation, the first term on the right-
hand side of Eq. (25) corresponds to a growing mode. The increase of $U_{i}$
is particularly prominent after $H^{2}$ drops below the order of $-m_{\rm
eff}^{2}$, during which $U_{i}\propto\exp(\sqrt{|m_{\rm eff}^{2}|}\,t)$. Then,
for $\alpha<0$, the FLRW background is destroyed by the growth of vector
perturbations. To avoid such an instability problem, we will focus on the
coupling in the range

$\alpha>0\,,$ (26)

in the following.

We recall that, in the above discussion, we have focused on the case
$|\alpha|\lesssim\beta$. In the coupling range $|\alpha|\gg\beta$ with $\beta
H_{i}^{2}$ of order 1, $|m_{\rm eff}^{2}|$ is much smaller than $H^{2}$ during
inflation. Taking the limit $|m_{\rm eff}^{2}|/H^{2}\to 0$ in Eq. (25), the
rapid growth of $U_{i}$ is absent. In Sec. V, however, we will show that
scalar perturbations are subject to exponential instabilities in the coupling
range $|\alpha|\gg\beta$.

##  IV Tensor perturbations

We proceed to study the dynamics of tensor perturbations $h_{ij}$ with the
perturbed line element given by

${\rm d}s^{2}=-{\rm d}t^{2}+a^{2}(t)\left(\delta_{ij}+h_{ij}\right){\rm
d}x^{i}{\rm d}x^{j}\,,$ (27)

with the traceless and transverse conditions ${h_{i}}^{i}=0$ and
$\partial^{i}h_{ij}=0$. Without loss of generality, we can consider
gravitational waves propagating along the $z$ direction, whose nonvanishing
components are

$h_{11}=-h_{22}=\frac{h_{1}(t,z)}{\sqrt{2}}\,,\qquad
h_{12}=h_{21}=\frac{h_{2}(t,z)}{\sqrt{2}}\,,$ (28)

where $h_{1}$ and $h_{2}$ are functions of $t$ and $z$. These components of
$h_{ij}$ automatically satisfy the traceless and transverse conditions
mentioned above.

Expanding the action (3) up to second order in $h_{1}$ and $h_{2}$ and
integrating the action by parts, we obtain the following quadratic-order
action

${\cal S}_{t}^{(2)}=\frac{M_{\rm Pl}^{2}}{4}\int{\rm
d}^{4}x\,a^{3}\sum_{i=1,2}\left[-\alpha\ddot{h}_{i}^{2}-\frac{\alpha}{a^{4}}h_{i}^{\prime\prime
2}+\frac{2\alpha}{a^{2}}\dot{h}_{i}^{\prime
2}+\left\\{\frac{1}{2}+\left(\alpha+6\beta\right)\left(2H^{2}+\dot{H}\right)\right\\}\dot{h}_{i}^{2}-\left\\{\frac{1}{2}+6\beta\left(2H^{2}+\dot{H}\right)\right\\}\frac{h_{i}^{\prime
2}}{a^{2}}\right]\,,$ (29)

where we also used the background Eq. (5). The presence of the Weyl curvature
term gives rise to the fourth-order differential equation for $h_{i}$. We
perform the Fourier transformation of $h_{i}$, as

$h_{i}=\frac{1}{(2\pi)^{1/2}}\int{\rm d}k\,\tilde{h}_{i}(t,k)e^{ikz}\,,$ (30)

where $\tilde{h}_{i}$ depends on $t$ and the wavenumber $k$. In Fourier space,
the second-order action $\tilde{{\cal S}}_{t}^{(2)}=\int{\rm d}^{4}x\,L_{t}$
can be obtained under the replacements
$\ddot{h}_{i}^{2}\to\ddot{\tilde{h}}_{i}^{2}$, $h_{i}^{\prime\prime 2}\to
k^{4}\tilde{h}_{i}^{2}$, $\dot{h}_{i}^{\prime 2}\to
k^{2}\dot{\tilde{h}}_{i}^{2}$, $\dot{h}_{i}^{2}\to\dot{\tilde{h}}_{i}^{2}$,
and $h_{i}^{\prime 2}\to k^{2}h_{i}^{2}$ in Eq. (29). Omitting the tilde from
$\tilde{h}_{i}$, the second-order Lagrangian density in Fourier space is given
by

$L_{t}=\frac{a^{3}M_{\rm
Pl}^{2}}{4}\sum_{i=1,2}\left[-\alpha\ddot{h}_{i}^{2}+\left\\{\frac{1}{2}+\left(\alpha+6\beta\right)\left(2H^{2}+\dot{H}\right)+\frac{2\alpha
k^{2}}{a^{2}}\right\\}\dot{h}_{i}^{2}-\frac{k^{2}}{a^{2}}\left\\{\frac{1}{2}+6\beta\left(2H^{2}+\dot{H}\right)+\frac{\alpha
k^{2}}{a^{2}}\right\\}h_{i}^{2}\right].$ (31)

To understand the appearance of ghost DOFs, we introduce Lagrangian
multipliers $\chi_{i}$ (with $i=1,2$) such that

$\bar{L}_{t}=L_{t}+\frac{a^{3}M_{\rm
Pl}^{2}\alpha}{4}\sum_{i=1,2}\left(\ddot{h}_{i}+c_{1}{\dot{h}}_{i}+c_{2}h_{i}-c_{3}\chi_{i}\right)^{2}\,,$
(32)

where $c_{i}$’s are time-dependent coefficients. We note that the coefficient
$a^{3}M_{\rm Pl}^{2}\alpha/4$ in front of $\ddot{h}_{i}^{2}$ has been
introduced to cancel the first term in Eq. (31). We fix $c_{i}$’s in Eq. (32)
to obtain the Lagrangian density containing a canonical form of the kinetic
terms $\dot{h}_{i}$ and $\dot{\chi}_{i}$, without the product
$\dot{h}_{i}\chi_{i}$. For this purpose, we choose

$\displaystyle c_{1}=3H\,,\qquad
c_{2}=\frac{k^{2}}{a^{2}}+\frac{4\alpha(H^{2}-\dot{H})+12\beta(2H^{2}+\dot{H})-3}{4\alpha}\,,\qquad
c_{3}=\frac{2}{\alpha}\,.$ (33)

After the integration by parts, the Lagrangian density (32) reduces to

$\displaystyle\bar{L}_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,a^{3}M_{\rm
Pl}^{2}\sum_{i=1,2}\Biggl{(}{\dot{h}}_{i}^{2}+2{\dot{h}}_{i}{\dot{\chi}}_{i}+\biggl{[}3\beta(H^{2}-\dot{H})(2H^{2}+\dot{H})-\frac{5H^{2}}{4}+\frac{\dot{H}}{2}-\frac{k^{2}}{a^{2}}+\frac{9\\{1-4\beta(2H^{2}+\dot{H})\\}^{2}}{32\alpha}$
(34)
$\displaystyle\qquad\qquad\qquad+\frac{\alpha}{12}\left\\{6H^{2}(H^{2}-2\dot{H})-\frac{1}{\beta}(H^{2}-\dot{H})-12\dot{H}\left(\frac{2k^{2}}{a^{2}}-5\dot{H}\right)\right\\}\biggr{]}h_{i}^{2}+\frac{2}{\alpha}\chi_{i}^{2}$
$\displaystyle\qquad\qquad\qquad-\left\\{2H^{2}-2\dot{H}+\frac{2k^{2}}{a^{2}}+\frac{3[4\beta(2H^{2}+\dot{H})-1]}{2\alpha}\right\\}h_{i}\chi_{i}\Biggr{)}\,.$

From this expression, we find that there are four dynamical perturbations
$h_{1}$, $h_{2}$, $\chi_{1}$, $\chi_{2}$ in the tensor sector, in agreement
with the analysis of Refs. [41, 38, 42]. Terms containing the product of time
derivatives in Eq. (34) can be expressed in the form

$(\bar{L}_{t})_{K}=\sum_{i=1,2}\dot{\vec{\psi}}_{i}{\bm{K}}\dot{\vec{\psi}}_{i}^{\,{\rm
T}}\,,$ (35)

where $\vec{\psi}_{i}=(h_{i},\chi_{i})$, and ${\bm{K}}$ is a $2\times 2$
symmetric matrix whose components are

$K_{11}=\frac{1}{2}\,a^{3}M_{\rm Pl}^{2}\,,\qquad
K_{12}=K_{21}=\frac{1}{2}\,a^{3}M_{\rm Pl}^{2}\,,\qquad K_{22}=0\,.$ (36)

The absence of ghosts requires the following conditions

$K_{11}=\frac{1}{2}\,a^{3}M_{\rm Pl}^{2}>0\,,\qquad{\rm and}\qquad-
K_{12}^{2}=-\frac{1}{4}\,a^{6}M_{\rm Pl}^{4}>0\,.$ (37)

While the former is satisfied, the latter is always violated. Hence the two
ghosts are present, besides the other two no-ghost modes.

In the following, we will study the propagation of tensor perturbations during
inflation. Varying the Lagrangian density (32) with respect to $\chi_{i}$, it
follows that

$\chi_{i}=\frac{1}{2}\alpha\ddot{h}_{i}+\frac{3}{2}\alpha
H\dot{h}_{i}+\frac{1}{2}\left[\frac{\alpha
k^{2}}{a^{2}}+\alpha(H^{2}-\dot{H})+3\beta(2H^{2}+\dot{H})-\frac{3}{4}\right]h_{i}\,,$
(38)

where we used the coefficients (33). We also vary Eq. (34) with respect to
$h_{i}$ and eliminate $\chi_{i}$ and their time derivatives by exploiting Eq.
(38). Then, we obtain the following fourth-order differential equation

$\displaystyle\ddddot{h_{i}}$
$\displaystyle=-6H\dddot{h_{i}}-\left[\frac{2k^{2}}{a^{2}}+\frac{\left(4\alpha+6\beta\right){\dot{H}}}{\alpha}+\frac{1+\left(22\alpha+24\beta\right)H^{2}}{2\alpha}\right]\ddot{h_{i}}$
$\displaystyle\quad-\left[\frac{2Hk^{2}}{a^{2}}+\frac{\left(\alpha+6\beta\right){\dot{H}}^{2}}{2\alpha
H}+\frac{4(\alpha+6\beta)H{\dot{H}}}{\alpha}+\frac{H\\{72\beta(\alpha+6\beta)H^{2}-\alpha+12\beta\\}}{12\alpha\beta}\right]\dot{h_{i}}$
$\displaystyle\quad-\left[\frac{k^{4}}{a^{4}}+\left(\frac{1+24\beta
H^{2}}{2\alpha}+\frac{6\beta{\dot{H}}}{\alpha}\right)\frac{k^{2}}{a^{2}}\right]h_{i}\,.$
(39)

We note that the same equation also follows by directly varying the original
Lagrangian density (31) with respect to $h_{i}$. Equation (39) governs the
dynamics of tensor perturbations.

We first solve Eq. (39) in the high momentum regime, namely, for the modes
deep inside the Hubble radius ($k/a\gg H$). After the Hubble radius crossing
during inflation, the perturbations enter the super-Hubble region $k/a<H$. The
evolution of $h_{i}$ in the latter large-scale regime will be discussed later.
Keeping only the most dominant terms for sub-Hubble perturbations and
expressing the Fourier components as
$h_{i}(t,k)=\tilde{h}_{i}(t)e^{-i\int\omega(t,k){\rm d}t}$, where
$\tilde{h}_{i}$ and $\omega$ depend on $t$ and $k=|\vec{k}|$, the WKB
approximation gives the relation $\dot{h}_{i}\simeq-i\omega h_{i}$. In this
WKB regime, we also have the inequality $|\dot{\omega}|\ll\omega^{2}$. Then,
Eq. (39) approximately reduces to

$\omega^{4}+6iH\omega^{3}-\frac{2k^{2}\omega^{2}}{a^{2}}-\frac{2iHk^{2}\omega}{a^{2}}+\frac{k^{4}}{a^{4}}\simeq
0\,.$ (40)

We search for solutions of the kind $\omega=c_{t}\,k/a$, where $c_{t}$ is the
tensor propagation speed. Substituting this dispersion relation into Eq. (40)
and taking the small-scale limit $k/(aH)\gg 1$, it follows that

$\left(c_{t}^{2}-1\right)^{2}\,\frac{k^{4}}{a^{4}}\,\left[1+\mathcal{O}\left(\frac{aH}{k}\right)\right]=0\,.$
(41)

At leading order in the expansion of the small parameter $aH/k$, we obtain

$c_{t}^{2}=1\,,$ (42)

for all the four dynamical modes $h_{1}$, $h_{2}$, $\chi_{1}$, $\chi_{2}$.
Since the propagation speeds are luminal, there are no classical Laplacian
instabilities in the tensor sector for perturbations deep inside the Hubble
radius.

Let us consider the evolution of super-Hubble tensor modes after the Hubble
radius crossing, i.e., $k/a\ll H$. Since $|\dot{H}|\ll H^{2}$ during slow-roll
inflation, Eq. (39) approximately reduces to

$\ddddot{h_{i}}+6H\dddot{h_{i}}+12\lambda_{1}H^{2}\ddot{h_{i}}+36\lambda_{2}H^{3}\dot{h_{i}}\simeq
0\,,$ (43)

where

$\lambda_{1}=\frac{(22\alpha+24\beta)H^{2}+1}{24\alpha
H^{2}}\,,\qquad\lambda_{2}=\frac{72\beta(\alpha+6\beta)H^{2}-\alpha+12\beta}{432\alpha\beta
H^{2}}\,.$ (44)

The difference between $\lambda_{1}$ and $\lambda_{2}$ is given by

$\lambda_{2}-\lambda_{1}=-\frac{3}{4}-\frac{1}{72\alpha
H^{2}}-\frac{1}{432\beta H^{2}}\,.$ (45)

For the coupling constants $\alpha$ and $\beta$ with $\alpha H^{2}\gtrsim{\cal
O}(1)$ and $\beta H^{2}\gtrsim{\cal O}(1)$ during inflation, we can neglect
the last two terms in Eq. (45) relative to $-3/4$. On using the approximate
relation $\lambda_{2}\simeq\lambda_{1}-3/4$ in this case and assuming that $H$
is constant during inflation, the solution to Eq. (43) can be expressed as

$h_{i}=A_{i}+B_{i}\,e^{-3Ht}+C_{i}\,e^{-[3-\sqrt{45-48\lambda_{1}}]Ht/2}+D_{i}\,e^{-[3+\sqrt{45-48\lambda_{1}}]Ht/2}\,,$
(46)

where $A_{i}$, $B_{i}$, $C_{i}$, and $D_{i}$ are integration constants, and

$\lambda_{1}\simeq\frac{11}{12}+\frac{\beta}{\alpha}>\frac{11}{12}\,.$ (47)

In the last inequality, we exploited the fact that both $\alpha$ and $\beta$
are positive111In the limit $\alpha\gg\beta$ and $\beta
H^{2}\gtrsim\mathcal{O}(1)$, we find that
$h_{i}=A_{i}+B_{i}e^{-3Ht}+C_{i}e^{-Ht}+D_{i}e^{-2Ht}$, so that $h_{i}$ tends
to $A_{i}$ also in this case.. For $11/12<\lambda_{1}\leq 15/16$, the last
three terms in Eq. (46) decrease exponentially. For $\lambda_{1}>15/16$, the
last two terms in Eq. (46) exhibit damped oscillations with a decreasing
amplitude proportional to $e^{-3Ht/2}$. This means that $h_{i}$ approaches the
constant value $A_{i}$.

If we consider the small Weyl coupling constant $\alpha\ll\beta$ with $\beta
H^{2}\gtrsim{\cal O}(1)$, then we have
$\lambda_{1}\simeq\lambda_{2}\simeq\beta/\alpha\gg 1$. In the limit that
$\beta/\alpha\to\infty$, the solution to Eq. (43) is given by

$h_{i}=A_{i}+B_{i}\,e^{-3Ht}+C_{i}\,e^{-(3+i\Omega)Ht/2}+D_{i}\,e^{-(3-i\Omega)Ht/2}\,,$
(48)

where $\Omega$ is a constant. In this case, the amplitude of $h_{i}$ decreases
in time as well and finally reaches a constant, $A_{i}$.

We have thus shown that, despite the appearance of ghosts, the tensor
perturbation is not subject to any classical instability. In other words, the
higher-order derivatives of $h_{i}$ appearing in the action (29) hardly affect
the standard conservation property of $h_{i}$, after the Hubble radius
crossing.

##  V Scalar perturbations

Let us next study the evolution of scalar perturbations for the perturbed line
element given by

${\rm d}s^{2}=-\left(1+2A\right){\rm d}t^{2}+2a(t)\partial_{i}B{\rm d}t{\rm
d}x^{i}+a^{2}(t)\left[(1+2\psi)\delta_{ij}+2\partial_{i}\partial_{j}E\right]{\rm
d}x^{i}{\rm d}x^{j}\,.$ (49)

We consider an infinitesimal gauge transformation

$\tilde{t}=t+\xi^{0}\,,\qquad\tilde{x}^{i}=x^{i}+\delta^{ij}\partial_{j}\xi\,,$
(50)

from one coordinate $x^{\mu}=(t,x^{i})$ to another coordinate
$\tilde{x}^{\mu}=(\tilde{t},\tilde{x}^{i})$, where $\xi^{0}$ and $\xi$ are
scalar quantities. Then, the four scalar perturbations $A$, $B$, $\psi$, and
$E$ transform, respectively, as [50, 51, 52, 53]

$\tilde{A}=A-\dot{\xi}^{0}\,,\qquad\tilde{B}=B+\frac{\xi^{0}}{a}-a\dot{\xi}\,,\qquad\tilde{\psi}=\psi-H\xi^{0}\,,\qquad\tilde{E}=E-\xi\,.$
(51)

The gauge-invariant gravitational potentials are defined by [50]

$\Psi=A+\frac{{\rm d}}{{\rm
d}t}\left[a(B-a\dot{E})\right]\,,\qquad\Phi=\psi+aH\left(B-a\dot{E}\right)\,.$
(52)

We can also construct the following gauge-invariant variables

$\mathcal{A}\equiv A-\frac{{\rm d}}{{\rm
d}t}\\!\left(\frac{\psi}{H}\right)\,,\qquad\mathcal{B}\equiv
aB+\frac{\psi}{H}-a^{2}\dot{E}\,.$ (53)

While ${\cal B}$ is related to $\Phi$ as ${\cal B}=\Phi/H$,
$\mathcal{A}=\Psi-({\rm d}/{{\rm d}t})\\!(\Phi/H)$ is not proportional to
$\Psi$.

It is known that $f(R)$ gravity with nonlinear functions of $R$ gives rise to
a scalar DOF $\phi={\rm d}f/{\rm d}R$ [5, 31]. The quadratic action (1)
contains the function $f(R)=R+\beta R^{2}$, in which case $\phi=1+2\beta R$.
Then, the perturbation of the new scalar DOF is equivalent to
$\delta\phi=2\beta\delta R$. We can construct a gauge-invariant quantity
$\zeta=\psi-H\delta\phi/\dot{\phi}$ [54, 55], or, equivalently222Note that the
field $\zeta$ is not well defined on an exact de Sitter space where $H$ is
constant. In the background for this model, the Hubble expansion rate varies
due to the slow-roll evolution of the scalaron DOF, so we can still introduce
$\zeta$.,

$\zeta=\psi-\frac{H}{\dot{R}}\delta R\,.$ (54)

There is also the following combination analogous to the Mukhanov-Sasaki
variable [56, 57]:

$\delta R_{f}=\delta R-\frac{\dot{R}}{H}\psi\,.$ (55)

which is related to $\zeta$, as $\delta R_{f}=-\dot{R}\zeta/H$.

We expand the action (3) up to quadratic order in scalar perturbations without
fixing gauges and then derive the field equations of motion by varying the
second-order action with respect to $A$, $B$, $\psi$, and $E$. These
perturbation equations of motion are written in a gauge-ready form [58, 59,
60], in that they are ready for the reader to choose a particular gauge. To
fix the spatial part of the gauge transformation vector $\xi^{\mu}$, we choose
the gauge

$E=0\,.$ (56)

In this case, the perturbation of the Ricci scalar is given by

$\delta
R=6\left(\ddot{\psi}+4H\dot{\psi}\right)-\frac{4}{a^{2}}\partial^{2}\psi-6H\dot{A}-12(2H^{2}+\dot{H})A-\frac{2}{a^{2}}\partial^{2}A-\frac{2}{a}\partial^{2}\left(\dot{B}+3HB\right).$
(57)

For the temporal part of $\xi^{\mu}$, we can consider several different gauge
choices, including (A) Newtonian gauge ($B=0$), (B) flat gauge ($\psi=0$), and
(C) unitary gauge ($\delta R=0$). The physical results, including the
stability conditions and the evolution of scalar perturbations, are
independent of the gauge choices.

We first study conditions for the absence of ghosts and then proceed to
address the dynamics of gauge-invariant perturbations for several different
gauge choices.

###  V.1 No-ghost conditions

In the flat gauge with $\psi=0$ and $E=0$, the gauge-invariant variables in
Eq. (53) reduce to ${\cal A}=A$ and ${\cal B}=aB$, respectively. We expand the
action (3) up to quadratic order in two perturbations ${\cal A}$ and ${\cal
B}$. In Fourier space, the second-order perturbed scalar Lagrangian density
$L_{s}$ contains products of the time derivatives of ${\cal A}$ and ${\cal B}$
in the form

$L_{s}\supset K_{11}\dot{\cal A}^{2}+K_{22}\dot{\cal B}^{2}+2K_{12}\dot{\cal
A}\dot{\cal B}\,,$ (58)

with the coefficients

$K_{11}=18a^{3}\beta H^{2}M_{\rm Pl}^{2}\,,\qquad
K_{22}=-\frac{2(\alpha-3\beta)M_{\rm Pl}^{2}k^{4}}{3a}\,,\qquad
K_{12}=-6a\beta HM_{\rm Pl}^{2}k^{2}\,.$ (59)

Thus, both the fields ${\cal A}$ and ${\cal B}$ propagate as dynamical
perturbations. The absence of ghosts requires the following two conditions

$\displaystyle K_{11}$ $\displaystyle=$ $\displaystyle 18a^{3}\beta
H^{2}M_{\rm Pl}^{2}>0\,,$ (60) $\displaystyle K_{11}K_{22}-K_{12}^{2}$
$\displaystyle=$ $\displaystyle-12a^{2}\alpha\beta H^{2}M_{\rm
Pl}^{4}k^{4}>0\,.$ (61)

The kinetic term $K_{11}\dot{\cal A}^{2}$ corresponds to that of the scalaron
perturbation $\mathcal{A}$, which does not behave as a ghost for $\beta>0$.
Under the no-ghost condition $\beta>0$ of the scalaron, the second equality
(61) is always violated for $\alpha>0$. Hence the other propagating DOF
corresponds to a ghost mode, which is induced by the presence of the Weyl
curvature term. To see this property more explicitly, we define the following
fields that make the kinetic matrix diagonal

${\cal A}_{1}\equiv{\cal A}+\frac{K_{12}}{K_{11}}{\cal B}\,,\qquad{\cal
B}_{1}\equiv{\cal B}\,.$ (62)

Then, the products of the time derivative of the new fields can be expressed
in the form

$L_{s}\supset K_{11}\dot{\cal
A}_{1}^{2}+\frac{K_{11}K_{22}-K_{12}^{2}}{K_{11}}\dot{\cal
B}_{1}^{2}=a^{3}\left(18\beta H^{2}M_{\rm Pl}^{2}\dot{\cal
A}_{1}^{2}-\frac{2M_{\rm Pl}^{2}\alpha k^{4}}{3a^{4}}\dot{\cal
B}_{1}^{2}\right)\,.$ (63)

For positive values of $\alpha$ and $\beta$, we also introduce the canonically
normalized fields

${\cal A}_{2}\equiv 6M_{\rm Pl}H\sqrt{\beta}{\cal A}_{1}\,,\qquad{\cal
B}_{2}\equiv M_{\rm Pl}\sqrt{\frac{4\alpha}{3}}\frac{k^{2}}{a^{2}}{\cal
B}_{1}\,,$ (64)

so that the kinetic products of the Lagrangian density can be expressed as

$(L_{s})_{K}=a^{3}\left(\frac{1}{2}\dot{{\cal
A}}_{2}^{2}-\frac{1}{2}\dot{{\cal B}}_{2}^{2}\right)\,.$ (65)

From this expression, it is clear that ${\cal A}_{2}$ and ${\cal B}_{2}$ are
the canonically normalized perturbations corresponding to the no-ghost
scalaron field and the Weyl scalar ghost, respectively.

In the above discussion we have chosen the flat gauge, but independent of the
gauge choices, the scalar sector contains one ghost and the other no-ghost
DOFs. In summary, for positive values of $\alpha$ and $\beta$, there are two
vector ghosts, two tensor ghosts, and one scalar ghost among the total eight
propagating DOFs333The recent analysis of black hole perturbations in Weyl
gravity without the $\beta R^{2}$ term [61] shows that there are seven
dynamical DOFs on a static and spherically symmetric background. Adding the
$\beta R^{2}$ term gives rise to one scalar DOF, so the total dynamical DOFs
match each other on two different backgrounds.. In the following, we study the
evolution of scalar perturbations in detail by choosing three different
gauges.

###  V.2 Newtonian gauge

The Newtonian gauge corresponds to setting $B=0$ and $E=0$. In this case the
gauge-invariant gravitational potentials in Eq. (52) reduce to $\Psi=A$ and
$\Phi=\psi$, so the perturbed line element is given by

${\rm d}s^{2}=-\left(1+2\Psi\right){\rm
d}t^{2}+a^{2}(t)\left(1+2\Phi\right)\delta_{ij}{\rm d}x^{i}{\rm d}x^{j}\,.$
(66)

In the Newtonian gauge, the coordinate transformation vector $\xi^{\mu}$ is
fixed on the FLRW background without any singularity. For this gauge choice,
both $\Psi$ and $\Phi$ play the role of two dynamical perturbations. From Eqs.
(54) and (57), the relation between $\zeta$ and the gravitational potentials
is

$\zeta=\Phi-\frac{H}{\dot{R}}\left[6(\ddot{\Phi}+4H\dot{\Phi})-\frac{4}{a^{2}}\partial^{2}\Phi-6H\dot{\Psi}-12(2H^{2}+\dot{H})\Psi-\frac{2}{a^{2}}\partial^{2}\Psi\right]\,.$
(67)

Varying the second-order action of scalar perturbations with respect to $A$,
$B$, $\psi$, and $E$, we obtain the four perturbation equations of motion,
respectively, for which we express in the form

${\cal E}_{A}=0\,,\qquad{\cal E}_{B}=0\,,\qquad{\cal E}_{\psi}=0\,,\qquad{\cal
E}_{E}=0\,,$ (68)

After setting $B=0=E$ in the end, two of the above equations are independent,
but the other two equations can be also used to obtain the closed differential
equations for $\Psi$ and $\Phi$. We also perform the Fourier transformation

$\Psi=\frac{1}{(2\pi)^{3}}\int{\rm
d}^{3}k\,\tilde{\Psi}(t,k)e^{i{\bm{k}}\cdot{\bm{r}}}\,,\qquad\Phi=\frac{1}{(2\pi)^{3}}\int{\rm
d}^{3}k\,\tilde{\Phi}(t,k)e^{i{\bm{k}}\cdot{\bm{r}}}\,,$ (69)

and work in the Fourier space by omitting the tilde from the perturbed
quantities.

For example, the fourth-order differential equation of $\Psi$ can be derived
by the following procedure. We first solve the two equations ${\cal E}_{A}=0$
and ${\cal E}_{B}=0$ for $\dddot{\Phi}$ and $\ddot{\Phi}$. Taking the time
derivative of the latter and combining it with the former, we can eliminate
the term $\dddot{\Phi}$ to obtain the other equation containing $\ddot{\Phi}$.
Then, we can solve for $\dot{\Phi}$ and $\ddot{\Phi}$ by combining the two
equations containing $\ddot{\Phi}$. Performing a similar procedure further, we
can express $\dot{\Phi}$ and $\Phi$ in terms of the derivatives of $\Psi$ up
to third order. Taking the time derivative of the $\Phi$ equation and
eliminating the $\dot{\Phi}$ term, we obtain the fourth-order differential
equation of $\Psi$ in the form

$\ddddot{\Psi}+\mu_{1}H\dddot{\Psi}+\mu_{2}H^{2}\ddot{\Psi}+\mu_{3}H^{3}\dot{\Psi}+\mu_{4}H^{4}\Psi=0\,,$
(70)

where $\mu_{1,2,3,4}$ are time-dependent dimensionless functions containing
the $k$ dependence. Due to its complexities, we do not write the explicit
forms of these coefficients.

Similarly, we can derive the fourth-order differential equation of $\Phi$ in
the form

$\displaystyle\ddddot{\Phi}+\nu_{1}H\dddot{\Phi}+\nu_{2}H^{2}\ddot{\Phi}+\nu_{3}H^{3}\dot{\Phi}+\nu_{4}H^{4}\Phi=0\,,$
(71)

where the functions $\nu_{i}$ ($i=1,2,3,4$) are not the same as $\mu_{i}$,
respectively.

For the modes deep inside the Hubble radius ($k\gg aH$), the coefficients in
Eqs. (70) and (71) reduce to

$\mu_{1}\simeq\nu_{1}\simeq
6\,,\qquad\mu_{2}\simeq\nu_{2}\simeq\frac{2k^{2}}{(aH)^{2}}\,,\qquad\mu_{3}\simeq\nu_{3}\simeq\frac{2k^{2}}{(aH)^{2}}\,,\qquad\mu_{4}\simeq\nu_{4}\simeq\frac{k^{4}}{(aH)^{4}}\,.$
(72)

In the WKB regime, we substitute the solutions
$\Psi=\Psi_{0}e^{-i\int\omega{\rm d}t}$ and $\Phi=\Phi_{0}e^{-i\int\omega{\rm
d}t}$ into Eqs. (70) and (71) with the coefficients (72), where $\Psi_{0}$ and
$\Phi_{0}$ are constants. This process leads to the same relation as Eq. (40).
Writing the dispersion relation as $\omega=c_{s}k/a$ and taking the limit
$k/a\gg H$, we obtain the squared propagation speeds

$c_{s}^{2}=1\,,$ (73)

for both $\Psi$ and $\Phi$. This means that, for the modes deep inside the
Hubble radius, the Laplacian instabilities are absent for the two
gravitational potentials.

For super-Hubble modes, we take the limit $k/(aH)\ll 1$ in the coefficients
$\mu_{i}$ and $\nu_{i}$. Moreover, we also take the limit where the slow-roll
parameter $\epsilon=-\dot{H}/H^{2}$ goes to 0. Then, the coefficients are
simplified to

$\displaystyle\mu_{1}$ $\displaystyle\simeq$
$\displaystyle\frac{864\beta^{2}(\alpha+\beta)H^{4}+18\beta(\alpha+4\beta)H^{2}-\alpha+3\beta}{6\beta
H^{2}[48\beta(\alpha+3\beta)H^{2}+\alpha+12\beta]}\,,$ (74)
$\displaystyle\mu_{2}$ $\displaystyle\simeq$
$\displaystyle-\frac{288(\alpha-12\beta)\beta^{2}(\alpha+3\beta)H^{4}-6\beta(7\alpha^{2}+36\alpha\beta+216\beta^{2})H^{2}+\alpha^{2}-6\alpha\beta-36\beta^{2}}{6\alpha\beta
H^{2}[48\beta(\alpha+3\beta)H^{2}+\alpha+12\beta]}\,,$ (75)
$\displaystyle\mu_{3}$ $\displaystyle\simeq$
$\displaystyle-\frac{1728\beta^{2}(\alpha^{2}-15\alpha\beta-24\beta^{2})H^{6}-12\beta(17\alpha^{2}+120\alpha\beta+216\beta^{2})H^{4}-\alpha(7\alpha+36\beta)H^{2}+\alpha-3\beta}{12\alpha\beta
H^{4}[48\beta(\alpha+3\beta)H^{2}+\alpha+12\beta]}\,,$ (76)
$\displaystyle\mu_{4}$ $\displaystyle\simeq$
$\displaystyle\frac{41472\beta^{4}(\alpha+\beta)H^{6}-96\beta^{2}(\alpha^{2}+33\alpha\beta+36\beta^{2})H^{4}+6\beta(\alpha^{2}-38\alpha\beta-72\beta^{2})H^{2}-\alpha(\alpha-3\beta)}{24\alpha\beta^{2}H^{4}[48\beta(\alpha+3\beta)H^{2}+\alpha+12\beta]}\,,$
(77)

and

$\displaystyle\nu_{1}$ $\displaystyle\simeq$
$\displaystyle\frac{36\beta(\alpha+\beta)H^{2}+2\alpha+3\beta}{12\beta(\alpha+3\beta)H^{2}}\,,$
(78) $\displaystyle\nu_{2}$ $\displaystyle\simeq$
$\displaystyle-\frac{2\beta(\alpha-12\beta)H^{2}-\alpha-\beta}{2\alpha\beta
H^{2}}\,,$ (79) $\displaystyle\nu_{3}$ $\displaystyle\simeq$
$\displaystyle-\frac{144\beta^{2}(\alpha^{2}-15\alpha\beta-24\beta^{2})H^{4}+4\beta(2\alpha^{2}-15\alpha\beta-36\beta^{2})H^{2}-2\alpha^{2}-7\alpha\beta-6\beta^{2}}{48\alpha\beta^{2}(\alpha+3\beta)H^{4}}\,,$
(80) $\displaystyle\nu_{4}$ $\displaystyle\simeq$
$\displaystyle\frac{864\beta^{2}(\alpha+\beta)H^{4}+72\beta(\alpha+\beta)H^{2}+2\alpha+3\beta}{24\alpha\beta(\alpha+3\beta)H^{4}}\,.$
(81)

Let us first consider the case in which the inequality

$\beta H^{2}\gg 1$ (82)

is satisfied during inflation. Then, Eqs. (74)-(77) and Eqs. (78)-(81)
approximately reduce to

$\mu_{1}\simeq\nu_{1}\simeq\frac{3(\alpha+\beta)}{\alpha+3\beta}\,,\qquad\mu_{2}\simeq\nu_{2}\simeq-\frac{\alpha-12\beta}{\alpha}\,,\qquad\mu_{3}\simeq\nu_{3}\simeq-\frac{3(\alpha^{2}-15\alpha\beta-24\beta^{2})}{\alpha(\alpha+3\beta)}\,,\qquad\mu_{4}\simeq\nu_{4}\simeq\frac{36\beta(\alpha+\beta)}{\alpha(\alpha+3\beta)}\,.$
(83)

On using these approximate coefficients and exploiting the approximation that
$H$ is constant during inflation, the solutions to Eqs. (70) and (71) in the
super-Hubble regime $k/(aH)\ll 1$ are given by

$\displaystyle\Psi=C_{1}e^{-Ht}+C_{2}e^{-\frac{3(\alpha+\beta)}{\alpha+3\beta}Ht}+C_{3}e^{(1+\sqrt{1-48\beta/\alpha})Ht/2}+C_{4}e^{(1-\sqrt{1-48\beta/\alpha})Ht/2}\,,$
(84)
$\displaystyle\Phi=D_{1}e^{-Ht}+D_{2}e^{-\frac{3(\alpha+\beta)}{\alpha+3\beta}Ht}+D_{3}e^{(1+\sqrt{1-48\beta/\alpha})Ht/2}+D_{4}e^{(1-\sqrt{1-48\beta/\alpha})Ht/2}\,,$
(85)

where $C_{i}$’s and $D_{i}$’s are constants. For $\alpha>0$ and $\beta>0$, the
third and fourth terms in Eqs. (84) and (85) grow in time, while the first and
second terms decay. Depending on the values of $\alpha$ and $\beta$, the
amplitudes of $\Psi$ and $\Phi$ increase as

$\displaystyle\\{|\Psi|,|\Phi|\\}\propto e^{Ht/2}\qquad{\rm
for}\quad\alpha<48\beta\,,$ (86) $\displaystyle\\{|\Psi|,|\Phi|\\}\propto
e^{(1+\sqrt{1-48\beta/\alpha})Ht/2}\qquad{\rm for}\quad\alpha>48\beta\,.$ (87)

In the coupling range (86), $\Psi$ and $\Phi$ exhibit oscillations with the
growing amplitudes. In the other coupling regime (87), $\Psi$ and $|\Phi|$
increase even faster than $e^{Ht/2}$. In the large Weyl coupling limit
$\alpha\gg 48\beta$, the gravitational potentials grow rapidly in proportion
to $e^{Ht}$.

The above results are valid for $\beta H^{2}$ exceeding the order 1 during
inflation. In particular, for the couplings in the ranges $\alpha\ll\beta$ and
$\beta H^{2}\gtrsim{\cal O}(1)$, the solutions to $\Psi$ and $\Phi$ correspond
to the limits $\beta/\alpha\to\infty$ in Eqs. (84) and (85). This means that
even a small Weyl coupling constant $\alpha$ induces the exponential growth of
gravitational potentials. Then, the homogeneous FLRW background is violated by
the rapid growth of $\Psi$ and $\Phi$ in the perturbed metric (66). In the
above analytic estimation we used the approximation $\beta H^{2}\gg 1$, but in
Sec. VI.3, we will show that, even for $\beta H^{2}={\cal O}(1)$, both the
amplitudes of $\Psi$ and $\Phi$ increase exponentially. The fact that $\Psi$
and $\Phi$ are subject to exponential growth does not depend on the gauge
choices either. Since all scalar perturbations are always determined through
two independent modes, at least one of them needs to be unstable. It should be
also pointed out that this instability is purely classical.

Besides $\Psi$ and $\Phi$, there are also other gauge-invariant scalar
perturbations. Let us consider the evolution of the curvature perturbation
$\zeta$ defined in Eq. (54). In the Newtonian gauge, we will describe a method
to find a closed differential equation for $\zeta$. First of all, we rewrite
the expression of $\zeta$ in Eq. (67), in terms of the gravitational
potentials in Fourier space, as

$\zeta=\Phi-\frac{H}{\dot{R}}\left[6(\ddot{\Phi}+4H\dot{\Phi})+\frac{4k^{2}}{a^{2}}\Phi-6H\dot{\Psi}-12(2H^{2}+\dot{H})\Psi+\frac{2k^{2}}{a^{2}}\Psi\right]\,.$
(88)

We can think of this relation as an equation that sets, on shell, $\zeta$ as a
function of the other fields. Therefore, we can proceed by adding a term to
the second-order scalar Lagrangian density $L_{s}$, as follows

$\bar{L}_{s}=L_{s}+b_{1}(t)\left\\{\Phi-\frac{H}{\dot{R}}\left[6(\ddot{\Phi}+4H\dot{\Phi})+\frac{4k^{2}}{a^{2}}\Phi-6H\dot{\Psi}-12(2H^{2}+\dot{H})\Psi+\frac{2k^{2}}{a^{2}}\Psi\right]-\zeta\right\\}^{2}\,,$
(89)

where $b_{1}(t)$ is a function of $t$. It is clear that at this level $\zeta$
is just a Lagrange multiplier, and its equation of motion, algebraic for
$\zeta$ itself, makes $\bar{L}_{s}$ reduce to the original Lagrangian density
$L_{s}$. We choose the coefficient $b_{1}(t)$ to cancel the term in
$\dot{\Psi}^{2}$, and, by doing so, also the term in $\ddot{\Phi}^{2}$ cancels
out from the Lagrangian.

After a few integrations by parts, we see that the field $\Psi$ can be set to
be a Lagrangian multiplier, and as such, it is integrated out from the
Lagrangian by using its equation of motion. By doing so, we arrive at an
equivalent Lagrangian density,
$\bar{L}_{s}=\bar{L}_{s}(\dot{\zeta},\dot{\Phi},\zeta,\Phi)$, which then
depends on two propagating DOFs, $\Phi$ and $\zeta$, as expected. The reduced
Lagrangian density contains the products of kinetic terms of the form

$\bar{L}_{s}\supset\bar{K}_{11}\dot{\zeta}^{2}+\bar{K}_{22}\dot{\Phi}^{2}+2\bar{K}_{12}\dot{\zeta}\dot{\Phi}\,,$
(90)

For positive values of $\alpha$ and $\beta$, we have $\bar{K}_{11}>0$ and
$\bar{K}_{11}\bar{K}_{22}-\bar{K}_{12}^{2}<0$ in the slow-roll limit. Hence
there is one ghost mode besides the other no-ghost mode. This property agrees
with the no-ghost conditions derived in Sec. V.1 for the flat gauge.

To derive the closed-form perturbation equation of $\zeta$, we proceed as
follows. From the equation of motion for the field $\Phi$, which we write in
the form ${\cal E}_{\Phi}=0$, we find an expression for $\ddot{\Phi}$ that can
be inserted into the equation for $\zeta$, i.e., ${\cal E}_{\zeta}=0$. Now, we
take the time derivative of this last equation to obtain $\dot{\cal
E}_{\zeta}=0$. We can still substitute this new equation into the expression
of $\ddot{\Phi}$, previously found, and solve it for $\dot{\Phi}$. We repeat
the step on considering now the equation $\ddot{\cal E}_{\zeta}=0$, and after
replacing it with the two expressions for $\ddot{\Phi}$ and $\dot{\Phi}$, we
can solve it with respect to $\Phi$ itself. At this point, we replace all
these $\Phi$-related expressions into the equation of motion ${\cal
E}_{\zeta}=0$. Then, we obtain the fourth-order differential equation

$\ddddot{\zeta}+\lambda_{1}H\dddot{\zeta}+\lambda_{2}H^{2}\ddot{\zeta}+\lambda_{3}H^{3}\dot{\zeta}+\lambda_{4}H^{4}\zeta=0\,,$
(91)

where $\lambda_{i}$’s are time-dependent dimensionless coefficients.

Taking the sub-Hubble limit $k/(aH)\gg 1$ with $\epsilon\to 0$, the
coefficients in Eq. (91) reduce to

$\lambda_{1}\simeq\frac{168\beta H^{2}+1}{12\beta
H^{2}}\,,\qquad\lambda_{2}\simeq\frac{2k^{2}}{(aH)^{2}}\,,\qquad\lambda_{3}\simeq\frac{(120\beta
H^{2}+1)k^{2}}{12\beta
H^{4}a^{2}}\,,\qquad\lambda_{4}\simeq\frac{k^{4}}{(aH)^{4}}\,.$ (92)

Using the solution $\zeta=\zeta_{0}e^{-i\int\omega{\rm d}t}$ under the WKB
approximation $\omega=c_{s}k/a\gg H$, the leading-order dispersion relation is
given by $\omega^{4}-2k^{2}\omega^{2}/a^{2}+k^{4}/a^{4}\simeq 0$. Hence the
curvature perturbation propagates with the luminal speed for the modes deep
inside the Hubble radius.

For super-Hubble perturbations, we take the limits $k/(aH)\to 0$ and
$\epsilon\to 0$ for the coefficients $\lambda_{i}$. For the coupling in the
range $\beta H^{2}\gg 1$, we approximately have

$\lambda_{1}\simeq\frac{8(8\alpha+3\beta)}{5\alpha+3\beta}\,,\qquad\lambda_{2}\simeq\frac{157\alpha^{2}+111\alpha\beta+36\beta^{2}}{\alpha(5\alpha+3\beta)}\,,\qquad\lambda_{3}\simeq\frac{2(\alpha+6\beta)(49\alpha+15\beta)}{\alpha(5\alpha+3\beta)}\,,\qquad\lambda_{4}\simeq
0\,.$ (93)

In this regime, we can integrate Eq. (91) to give

$\zeta=c_{1}+c_{2}e^{-\frac{(49\alpha+15\beta)Ht}{5\alpha+3\beta}}+c_{3}e^{(-3+\sqrt{1-48\beta/\alpha})Ht/2}+c_{4}e^{-(3+\sqrt{1-48\beta/\alpha})Ht/2}\,.$
(94)

For positive values of $\alpha$ and $\beta$, the last three terms in Eq. (94)
decay in time. Hence the curvature perturbation approaches a constant $c_{1}$
after the Hubble radius crossing. In the above estimation, we have used the
approximation $\beta H^{2}\gg 1$, but we have numerically confirmed that
$\zeta$ approaches a constant even for $\beta H^{2}={\cal O}(1)$.

Despite the exponential increase of two gravitational potentials during
inflation, there is a particular gauge-invariant combination $\zeta$ that does
not grow in the large-scale limit. Unlike the standard single-field slow-roll
inflation, however, we have two propagating DOFs in the scalar sector. In the
description of Lagrange multipliers explained above, the two dynamical DOFs
correspond to the perturbations $\Phi$ and $\zeta$. Even though $\zeta$ is not
enhanced after the Hubble radius crossing, the other dynamical field $\Phi$ is
subject to exponential growth. Thus, the analysis in the Newtonian gauge shows
that the Weyl curvature term violates the homogeneous inflationary background.

###  V.3 Flat gauge

In the flat gauge with the gauge conditions $\psi=0$ and $E=0$, we have ${\cal
A}=A$ and ${\cal B}=aB$ in Eq. (53). Then, the perturbed line element is given
by

${\rm d}s^{2}=-\left(1+2{\cal A}\right){\rm d}t^{2}+2\partial_{i}{\cal B}{\rm
d}t{\rm d}x^{i}+a^{2}(t)\delta_{ij}{\rm d}x^{i}{\rm d}x^{j}\,.$ (95)

On the expanding cosmological background ($H\neq 0$), the coordinate
transformation vector $\xi^{\mu}$ is always regular for the flat gauge. The
gauge-invariant variables ${\cal A}$ and ${\cal B}$ are related to $\Psi$ and
$\Phi$ according to

${\cal A}=\Psi-\frac{{\rm d}}{{\rm
d}t}\\!\left(\frac{\Phi}{H}\right)\,,\qquad{\cal B}=\frac{\Phi}{H}\,.$ (96)

While ${\cal B}$ is directly proportional to $\Phi$, ${\cal A}$ corresponds to
a combination of $\Psi$ and $\Phi$.

After setting $\psi=0=E$ in the perturbation equations of motion, the
dynamical system in the flat gauge has two propagating DOFs ${\cal A}$ and
${\cal B}$ (or $\Phi$). To derive the closed differential equation for ${\cal
A}$, we solve the two equations ${\cal E}_{A}=0$ and ${\cal E}_{B}=0$ for
$\ddot{\cal B}$ and $\dot{\cal B}$. Following a similar procedure to that
performed in the Newtonian gauge, we can express the terms $\dot{\cal B}$ and
${\cal B}$ in terms of ${\cal A}$ and its derivatives. Taking the time
derivative of the ${\cal B}$ equation and combining it with the $\dot{\cal B}$
equation, we obtain the fourth-order differential equation of ${\cal A}$ in
the form

$\ddddot{\cal A}+\tau_{1}H\dddot{\cal A}+\tau_{2}H^{2}\ddot{\cal
A}+\tau_{3}H^{3}\dot{\cal A}+\tau_{4}H^{4}{\cal A}=0\,,$ (97)

where $\tau_{i}$’s are time-dependent functions.

Taking the sub-Hubble limit $k\gg aH$ with $\epsilon\to 0$, the coefficients
in Eq. (97) reduce to

$\tau_{1}\simeq
10\,,\qquad\tau_{2}\simeq\frac{2k^{2}}{(aH)^{2}}\,,\qquad\tau_{3}\simeq\frac{6k^{2}}{(aH)^{2}}\,,\qquad\tau_{4}\simeq\frac{k^{4}}{(aH)^{4}}\,.$
(98)

On using the WKB approximation, it follows that the field ${\cal A}$
propagates with the speed of light.

For super-Hubble modes $k\ll aH$ with the coupling $\beta H^{2}\gg 1$, taking
the slow-roll limit $\epsilon\to 0$ gives

$\tau_{1}\simeq
6\,,\qquad\tau_{2}\simeq\frac{11\alpha+12\beta}{\alpha}\,,\qquad\tau_{3}\simeq\frac{6(\alpha+6\beta)}{\alpha}\,,\qquad\tau_{4}\simeq-\frac{\alpha+6\beta}{3\alpha\beta
H^{2}}\,.$ (99)

In this regime, Eq. (97) can be integrated to give

$\displaystyle{\cal A}$ $\displaystyle=c_{1}e^{-(9\alpha\beta
H+\sqrt{\Delta_{1}-12\alpha\beta
H\sqrt{\Delta_{2}}})t/(6\alpha\beta)}+c_{2}e^{-(9\alpha\beta
H-\sqrt{\Delta_{1}-12\alpha\beta H\sqrt{\Delta_{2}}})t/(6\alpha\beta)}$
$\displaystyle~{}~{}~{}+c_{3}e^{-(9\alpha\beta
H+\sqrt{\Delta_{1}+12\alpha\beta
H\sqrt{\Delta_{2}}})t/(6\alpha\beta)}+c_{4}e^{-(9\alpha\beta
H-\sqrt{\Delta_{1}+12\alpha\beta H\sqrt{\Delta_{2}}})t/(6\alpha\beta)}\,,$
(100)

where

$\displaystyle\Delta_{1}$ $\displaystyle\equiv$ $\displaystyle
9H^{2}\alpha\beta^{2}(5\alpha-24\beta)\,,$ (101) $\displaystyle\Delta_{2}$
$\displaystyle\equiv$ $\displaystyle 3\beta(\alpha+6\beta)\left[3\beta
H^{2}(\alpha+6\beta)+\alpha\right]\,.$ (102)

The amplitudes of the first two terms in Eq. (100) decrease in proportion to
$e^{-3Ht/2}$, while the third term decreases as $\propto e^{-3Ht}$. Taking the
limit $\beta H^{2}\gg 1$ for the last term of Eq. (100), the leading-order
solution is a constant term $c_{4}$. Picking up the next-to-leading
correction, we obtain the following solution

$\mathcal{A}\simeq c_{4}\left(1+\frac{t}{18\beta H}\right)\,.$ (103)

For the number of e-foldings $N$ of order 10, the correction induced by the
time-dependent terms in Eq. (103) is suppressed compared to the leading-order
constant term.

Following a similar procedure performed for the perturbation ${\cal A}$, we
can also derive the fourth-order differential equation of ${\cal B}$ in the
form

$\ddddot{\cal B}+\eta_{1}H\dddot{\cal B}+\eta_{2}H^{2}\ddot{\cal
B}+\eta_{3}H^{3}\dot{\cal B}+\eta_{4}H^{4}{\cal B}=0\,,$ (104)

where $\eta_{i}$’s are time-dependent coefficients which are not identical to
$\nu_{i}$’s in Eq. (71). If we use the variable $\Phi$ instead of ${\cal B}$,
the coefficients of the fourth-order differential equation of $\Phi$ exactly
coincide with those derived in the Newtonian gauge.

Using the WKB approximation for the modes deep inside the Hubble radius, the
perturbation ${\cal B}$ obeys

$\ddddot{\cal B}+\frac{2k^{2}}{a^{2}}\ddot{\cal B}+\frac{k^{4}}{a^{4}}{\cal
B}\simeq 0\,,$ (105)

so that ${\cal B}$ propagates with the speed of light. Taking the super-Hubble
limit ($k\ll aH$) with $\epsilon\to 0$, the coefficients in Eq. (104) reduce
to

$\displaystyle\eta_{1}$ $\displaystyle\simeq$
$\displaystyle\frac{36(\alpha+\beta)\beta
H^{2}+2\alpha+3\beta}{12(\alpha+3\beta)\beta H^{2}}\,,$ (106)
$\displaystyle\eta_{2}$ $\displaystyle\simeq$
$\displaystyle\frac{2(12\beta-\alpha)H^{2}+1}{2\alpha H^{2}}\,,$ (107)
$\displaystyle\eta_{3}$ $\displaystyle\simeq$
$\displaystyle-\frac{72\beta(\alpha^{2}-15\alpha\beta-24\beta^{2})H^{4}-2(\alpha^{2}+42\alpha\beta+36\beta^{2})H^{2}-2\alpha-3\beta}{24\alpha\beta(\alpha+3\beta)H^{4}}\,,$
(108) $\displaystyle\eta_{4}$ $\displaystyle\simeq$
$\displaystyle\frac{5184\beta^{3}(\alpha+\beta)H^{4}+12\alpha\beta(\alpha+9\beta)H^{2}+\alpha(2\alpha+3\beta)}{144\alpha\beta^{2}(\alpha+3\beta)H^{4}}\,.$
(109)

For $\beta H^{2}\gg 1$, the coefficients $\eta_{1,2,3,4}$ approximately reduce
to the values $\nu_{1,2,3,4}$ given in Eq. (83), respectively. The field
${\cal B}$ is subject to exponential growth during inflation analogous to
$\Phi$, see Eqs. (86)-(87). This means that, even though ${\cal A}$ does not
grow significantly, the other metric perturbation ${\cal B}$ in the line
element (95) increases rapidly to violate the FLRW background.

In the flat gauge, the perturbation $\zeta$ is given by

$\zeta=-\frac{H}{\dot{R}}\,\delta
R=-\frac{H}{\dot{R}}\left[\frac{2k^{2}\dot{\mathcal{B}}}{a^{2}}-6H\dot{\mathcal{A}}+\left(\frac{2k^{2}}{a^{2}}-24H^{2}-12\dot{H}\right)\mathcal{A}+\frac{4\mathcal{B}H\,k^{2}}{a^{2}}\right]\,.$
(110)

To derive the fourth-order differential equation of $\zeta$, we build the
following Lagrangian density

$\bar{L}_{s}=L_{s}+b_{2}(t)\left\\{\zeta+\frac{H}{\dot{R}}\left[\frac{2k^{2}\dot{\mathcal{B}}}{a^{2}}-6H\dot{\mathcal{A}}+\left(\frac{2k^{2}}{a^{2}}-24H^{2}-12\dot{H}\right)\mathcal{A}+\frac{4\mathcal{B}H\,k^{2}}{a^{2}}\right]\right\\}^{2}\,,$
(111)

which is equivalent to the original Lagrangian density $L_{s}$. The function
$b_{2}(t)$ needs to be chosen so that the kinetic term of $\mathcal{A}$
vanishes identically. After a few integrations by parts, the field
$\mathcal{A}$ becomes a Lagrange multiplier which can be integrated out,
leaving $\zeta$ and $\mathcal{B}$ as two dynamical DOFs. After varying
$\bar{L}_{s}$ with respect to $\zeta$ and $\mathcal{B}$, we can proceed along
the same lines as finding the equations of motion for the fields $\mathcal{A}$
and ${\cal B}$. This leads to the closed differential equation of $\zeta$ in
the form (91), where $\lambda_{i}$’s are exactly the same as those derived in
the Newtonian gauge. Hence $\zeta$ approaches a constant after the Hubble
radius crossing. However, the fact that the other dynamical perturbation
${\cal B}$ grows exponentially means that the instability of the FLRW
background cannot be avoided. We also note that the fourth-order differential
equation of $\Psi$ exactly coincides with the one obtained in the Newtonian
gauge. Hence the two gravitational potentials $\Psi$ and $\Phi$ are unstable
in the flat gauge as well, by reflecting the fact that both $\Psi$ and $\Phi$
contain the dependence of $aB$.

###  V.4 Unitary gauge

Let us finally discuss the evolution of scalar perturbations in the unitary
gauge with $\delta R=0$. Since the curvature perturbation $\zeta$ is
equivalent to $\psi$, the gauge condition in Fourier space translates to

$\delta
R=6\left(\ddot{\zeta}+4H\dot{\zeta}\right)+\frac{4k^{2}}{a^{2}}\zeta-6H\dot{A}-12\left(2H^{2}+\dot{H}\right)A+\frac{2k^{2}}{a^{2}}A+\frac{2k^{2}}{a}\left(\dot{B}+3HB\right)=0\,.$
(112)

On using this condition together with the perturbation equations of motion, we
can derive the fourth-order differential equation of $\zeta$. We first solve
the perturbation equation ${\cal E}_{A}=0$ for $\ddot{B}$ and take the time
derivative of Eq. (112) to obtain the first derivative $\dot{B}$. Using Eq.
(112) to solve for $\dot{B}$, one can express $B$ and its time derivatives in
terms of $\zeta,A$ and their time derivatives. In this way, all the
$B$-dependent quantities can be eliminated from the perturbation equations of
motion. The next step is to remove the $A$-dependent terms. On using the two
equations ${\cal E}_{\psi}=0$ and ${\cal E}_{E}=0$, we can solve for
$\dddot{A}$ and $\ddot{A}$. Then, following a similar procedure as before, it
is possible to express $A$ in terms of the derivatives of $\zeta$ up to third
order. Taking the time derivative of this equation and eliminating the
$\dot{A}$ term, we obtain the fourth-order differential equations of $\zeta$
with the exactly same coefficients as $\lambda_{1,2,3,4}$ in Eq. (91). Then,
the constancy of $\zeta$ after the Hubble radius crossing also holds in the
unitary gauge. This result is consistent with the analysis of Ref. [42] in the
Einstein frame. Similarly, we obtain the same fourth-order differential
equation of ${\cal A}$ as Eq. (97), so the solution in the super-Hubble regime
is given by Eq. (103).

The closed differential equations of $\Psi$ and $\Phi$ can be also obtained by
using the following relations

$A=\Psi-\dot{\cal B}_{u}\,,\qquad\zeta=\Phi-H{\cal B}_{u}\,,$ (113)

where

${\cal B}_{u}\equiv aB\,.$ (114)

The gauge condition (112) and the perturbation equations of motion can be now
expressed in terms of the gauge-invariant variables $\Psi$, $\Phi$, ${\cal
B}_{u}$ and their time derivatives. Indeed, the above change of variables
automatically removes the ${\cal B}_{u}$-dependent terms from the two
perturbation equations ${\cal E}_{A}=0$ and ${\cal E}_{B}=0$. Combining these
two, it is straightforward to derive the fourth-order differential equations
of $\Psi$ and $\Phi$. Again, we find that they are identical to Eqs. (70) and
(71) derived in the Newtonian gauge, respectively, with the completely same
coefficients. Hence the same instabilities of $\Psi$ and $\Phi$ are present
after the Hubble radius crossing, while the growth of $\zeta$ and ${\cal A}$
is suppressed.

In terms of the gauge-invariant variables, the perturbed line element in the
unitary gauge can be expressed as

${\rm d}s^{2}=-\left[1+2{\cal A}+2\frac{{\rm d}}{{\rm
d}t}\left(\frac{\zeta}{H}\right)\right]{\rm
d}t^{2}+\frac{2}{H}\left(\partial_{i}\Phi-\partial_{i}\zeta\right){\rm d}t{\rm
d}x^{i}+a^{2}(t)(1+2\zeta)\delta_{ij}{\rm d}x^{i}{\rm d}x^{j}\,.$ (115)

Due to the suppressed growth of ${\cal A}$ and $\zeta$, metric perturbations
in the $g_{00}$ and $g_{ij}$ components are not subject to instabilities.
However, the $\partial_{i}\Phi$ term in the $g_{0i}$ component exhibits an
exponential increase after the Hubble radius crossing. Since the
$\partial_{i}\zeta$ term approaches a constant in the super-Hubble regime, the
dominance of $\partial_{i}\Phi$ over $\partial_{i}\zeta$ in $g_{0i}$ leads to
the instability of the FLRW background. Indeed, we numerically confirmed that
the gauge-invariant perturbation ${\cal B}_{u}=(\Phi-\zeta)/H$ grows
exponentially in the super-Hubble regime due to the enhancement of $\Phi$. We
have thus analytically shown that, for any physical gauge choices, the Weyl
curvature makes the inflationary FLRW background unstable.

##  VI Numerical simulations with the discussion of initial conditions

In this section, we will numerically confirm the instability of the FLRW
background in Weyl gravity with the $\beta R^{2}$ term. For this purpose, we
first discuss the choice of initial conditions of perturbations and then
proceed to the numerical analysis.

###  VI.1 Fourth-order system and initial conditions

We have learned so far that each of the considered perturbation fields obeys a
fourth-order differential equation, which, under the WKB approximation, is
solved as the solutions describing waves propagating with the speed of light.
We will show that it is indeed possible, starting from the reduced action of a
single scalar field $v$ possessing the term ${\ddot{v}}^{2}$, to find an
equivalent Lagrangian density of two scalar fields with second-order equations
of motion. The discussion in this section can be applied to any dynamical
perturbation $v$ with some gauge choices, but in Sec. VI.2 we will consider
the flat gauge for concreteness.

For the modes deep inside the Hubble radius, we should expect the field $v$ in
Fourier space to satisfy a fourth-order equation of motion, which can be
derived by the following approximate Lagrangian

$L_{v}\simeq
Q\left({\ddot{v}}^{2}-2\frac{k^{2}}{a^{2}}\,{\dot{v}}^{2}+\frac{k^{4}}{a^{4}}\,v^{2}\right)\,.$
(116)

Here and in the following, we assume that the function $Q=Q(t,k^{2})$ can be
either positive or negative. In Sec. VI.2, we will see that it is possible to
obtain a Lagrangian density444The Lagrangian, in this case, will take the
following general form $L_{v}=Q\,(\ddot{v}^{2}-2Q_{1}\dot{v}^{2}+Q_{2}v^{2})$,
where $Q_{i}$, ($i\in\\{1,2,3\\}$), are functions of time and $k^{2}$.
reducing to the form (116) in the high-$k$ regime with $v$ related to the
perturbations ${\cal A}$ and ${\cal B}$. We note that, in the action (116), we
are assuming the high momentum regime in which the WKB approximation holds for
the dynamics of $v$.

Then, we can introduce an auxiliary field, $w$, as

$L_{v}\simeq
Q\left[{\ddot{v}}^{2}-2\frac{k^{2}}{a^{2}}\,{\dot{v}}^{2}+\frac{k^{4}}{a^{4}}\,v^{2}-(b_{2}w+\ddot{v})^{2}\right]\,,$
(117)

where $b_{2}(t,k^{2})$ is a general function, so far undetermined. By
integrating out the field $w$, we find once more the original Lagrangian
density. Therefore, the two Lagrangian densities (116) and (117) are
equivalent to each other, both of which lead to the same dynamics. It is also
clear that the term in $\ddot{v}^{2}$ vanishes in Eq. (117) for the new
Lagrangian density. At this level, we can introduce the quantity

$v=v_{2}+\frac{a^{2}b_{2}}{2k^{2}}\,w\,,$ (118)

which is meant to diagonalize the kinetic matrix. Then, we perform the other
field redefinitions

$v_{2}=\frac{a}{k}\,v_{3}\,,\qquad w=\frac{k}{a}\,w_{3}\,,$ (119)

as to make the kinetic terms only background dependent, i.e., independent of
the wavenumber $k$.

Having assumed that $Q\neq 0$, we can further introduce the following field
redefinitions

$v_{3}=\frac{a^{3/2}}{2\sqrt{|Q|}}\,v_{4}\,,\qquad
w_{3}=\frac{a^{3/2}}{b_{2}\sqrt{|Q|}}\,w_{4}\,,$ (120)

to obtain canonical kinetic terms for the fields $v_{4}$ and $w_{4}$. At this
level, the Lagrangian for the modes deep inside the Hubble radius reduces to

$L_{v}\simeq-\mathrm{sign}(Q)\left[\frac{a^{3}}{2}\,(\dot{v}_{4}^{2}-\dot{w}_{4}^{2})-\frac{k^{2}a}{4}\,(v_{4}^{2}-3w_{4}^{2}+2v_{4}w_{4})\right].$
(121)

This Lagrangian still leads to the dynamics of perturbations propagating with
the speed of light. However, the effective mass matrix
$C=\tfrac{1}{4}\,k^{2}a\,\begin{pmatrix}1&1\\\ 1&-3\end{pmatrix}$ cannot be
diagonalized by any real (finite and constant-in-time) Lorentz transformation
that would leave instead the kinetic matrix in the canonical form. Therefore,
the two modes $v_{4}$ and $w_{4}$ are not completely decoupled.

Nonetheless, we can use the WKB approximation and look for solutions of the
kind $v_{4}\propto e^{-i\int\omega{\rm d}t}$ and $w_{4}\propto
e^{-i\int\omega{\rm d}t}$, for which $\ddot{v}_{4}\propto-\omega^{2}v_{4}$ and
$\ddot{w}_{4}\propto-\omega^{2}w_{4}$. In this case, it follows that the
dispersion relation $\omega=k/a$ needs to hold, which means that both modes
propagate with the luminal speeds. Furthermore, we obtain the relation
$w_{4}=v_{4}$, and hence the dynamics of the mode $w_{4}-v_{4}$ is set to
vanish. Still, proper initial conditions need to be imposed on the field
$v_{4}$ (or $w_{4}$).

To better understand the behavior of solutions or the choice of initial
conditions, we use the perturbation equations of motion in the high-$k$ regime
expressed in terms of the conformal time $\eta=\int a^{-1}{\rm d}t$,

$\frac{{\rm d}^{2}v_{4}}{{\rm
d}\eta^{2}}\simeq-\frac{k^{2}}{2}\,(v_{4}+w_{4})\,,\qquad\frac{{\rm
d}^{2}w_{4}}{{\rm d}\eta^{2}}\simeq\frac{k^{2}}{2}\,(v_{4}-3w_{4})\,.$ (122)

Then, we obtain the following general solutions

$\displaystyle w_{4}$
$\displaystyle=c_{1}\sin(k\eta)+c_{2}\cos(k\eta)+c_{3}\,k\eta\sin(k\eta)+c_{4}\,k\eta\cos(k\eta)\,,$
(123) $\displaystyle v_{4}$
$\displaystyle=c_{1}\sin(k\eta)+c_{2}\cos(k\eta)+c_{3}\,[k\eta\sin(k\eta)+4\cos(k\eta)]+c_{4}\,[k\eta\cos(k\eta)-4\sin(k\eta)]\,,$
(124)

where $c_{1}$, $c_{2}$, $c_{3}$, and $c_{4}$ are integration constants. The
choice of exact plane-wave initial conditions corresponds to $c_{3}=0=c_{4}$
so that $v_{4}=w_{4}$. As in the usual normalization scheme of the Bunch-
Davies vacuum555From a purely classical point of view, the initial condition
representing a plane wave consists of setting ${\rm d}v_{4}/{\rm
d}\eta=-ikv_{4}$, and so on for all other higher derivatives. The
normalization of $v_{4}$ itself is not established as the field $Cv_{4}$ is
still a solution of the equations of motion ($C$ is a constant). In this case,
the physical quantity to consider is $v_{4}(\eta)/v_{4}(\eta_{i})$., we choose
a positive frequency solution $v_{4}=w_{4}={\cal C}_{0}e^{-ik\eta}$ and impose
the conditions $\sqrt{|v_{4}^{2}||(a^{3}\dot{v}_{4})^{2}|}=1/2$ and
$\sqrt{|w_{4}^{2}||(a^{3}\dot{w}_{4})^{2}|}=1/2$. This fixes the coefficient
${\cal C}_{0}$ to be $1/(a\sqrt{2k})$, and hence

$v_{4}=w_{4}=\frac{1}{a\sqrt{2k}}\,e^{-i\int\frac{k}{a}{\rm d}t}\,.$ (125)

In terms of the original perturbation $v$, we have

$v=\frac{a^{5/2}}{2k\sqrt{|Q|}}\left(v_{4}+w_{4}\right)=\frac{a^{3/2}}{\sqrt{2k^{3}|Q|}}e^{-i\int\frac{k}{a}{\rm
d}t}\,.$ (126)

We will choose this as the initial condition of $v$ for the modes deep inside
the Hubble radius. Note that we will not discuss the quantization of
perturbations in our theory. We have already shown analytically the presence
of violent classical instabilities for the modes after the Hubble radius
crossing. In other words, the scalar ghost is not of the soft type in our
theory. The classical instabilities induced by the ghost make the quantization
of perturbations irrelevant. Only for a stable classical background, it would
be worth investigating the quantization procedure.

###  VI.2 Initial conditions for the perturbations ${\cal A}$ and ${\cal B}$

In Sec. VI.1, we assumed the existence of the Lagrangian density (116) leading
to the closed fourth-order differential equation. In this section, we will
prove its existence by considering perturbations in the flat gauge. For this
purpose, we will proceed as follows. In the flat gauge, the kinetic Lagrangian
for the fields $\mathcal{A}$ and $\mathcal{B}$ was already discussed in Eq.
(58). We remind the reader about the field redefinitions that were introduced
to obtain canonical kinetic terms. The field redefinitions given in Eqs. (62)
and (64) allow us to obtain the Lagrangian density of two canonically
normalized fields ${\cal A}_{2}$ and ${\cal B}_{2}$ in the form (65). The
original perturbations ${\cal A}$ and ${\cal B}$ are related to ${\cal A}_{2}$
and ${\cal B}_{2}$, as

${\cal A}=\frac{1}{6M_{\rm Pl}H}\left(\frac{{\cal
A}_{2}}{\sqrt{\beta}}+\frac{\sqrt{3}{\cal
B}_{2}}{\sqrt{\alpha}}\right)\,,\qquad{\cal B}=\frac{\sqrt{3}a^{2}{\cal
B}_{2}}{2M_{\rm Pl}k^{2}\sqrt{\alpha}}\,.$ (127)

From the canonical expression (65), we wish to find the Lagrangian density of
the form (116). We perform the following field redefinitions

$\mathcal{A}_{2}=\frac{1}{\sqrt{2}}\,(\mathcal{Y}-\mathcal{Z})\,,\qquad\mathcal{B}_{2}=\frac{1}{\sqrt{2}}\,(\mathcal{Y}+\mathcal{Z})\,,$
(128)

so that the kinetic terms are simplified to

$L_{s}\ni-a^{3}\dot{\mathcal{Y}}\dot{\mathcal{Z}}+\dots=\mathcal{Z}\,\frac{{\rm
d}}{{\rm d}t}(a^{3}\dot{\mathcal{Y}})+\dots\,,$ (129)

up to total derivatives. Varying the full Lagrangian density $L_{s}$ with
respect to $\mathcal{Z}$, it follows that $\mathcal{Z}$ can be expressed in
terms of $\mathcal{Y}$ and its first and second-time derivatives. Substituting
$\mathcal{Z}$ into the second expression of Eq. (129), we find that the
Lagrangian density can be expressed in the form

$L_{s}=Q\,(\ddot{\mathcal{Y}}^{2}-2Q_{1}\,\dot{\mathcal{Y}}^{2}+Q_{2}\,\mathcal{Y}^{2})\,,$
(130)

where $Q$, $Q_{1}$, and $Q_{2}$ depend on $t$ and $k^{2}$. It should be
noticed that, up to this point, we have not made any approximation for
particular wavenumbers. However, if we look for the behavior of the three
quantities $Q$, $Q_{1}$, $Q_{2}$ in the high-$k$ regime, we find

$Q\simeq\frac{27\alpha\beta
a^{7}H^{2}}{\left(\alpha-3\beta\right)\left(\sqrt{\alpha}-\sqrt{3\beta}\right)^{2}k^{4}}\,,\qquad
Q_{1}\simeq\frac{k^{2}}{a^{2}}\,,\qquad Q_{2}\simeq\frac{k^{4}}{a^{4}}\,,$
(131)

which means that666For the special case where $\alpha=3\beta$, by looking at
Eq. (58), we can see that the field $\mathcal{B}$ can be easily integrated out
in terms of $\mathcal{A}$ and its first and second time derivatives.
Alternatively, we can still use the procedure described here, but now
$Q\simeq-a^{3}/(6H^{2})$ and $\mathcal{Z}\simeq
ik\mathcal{Y}/(aH)\gg\mathcal{Y}$.
$\mathrm{sign}(Q)=\mathrm{sign}(\alpha-3\beta)$. The behavior of $Q_{1}$ and
$Q_{2}$ makes sure that the propagation of the mode ${\cal Y}$ is luminal in
the WKB approximation scheme.

On using the discussion given in Sec. VI.1, ${\cal Y}$ plays the role of the
field $v$ with $Q$ given in Eq. (131). Then, we can choose the initial
condition of ${\cal Y}$, as

${\cal Y}=\frac{a^{3/2}}{\sqrt{2k^{3}|Q|}}e^{-i\int\frac{k}{a}{\rm
d}t}\simeq\frac{k^{1/2}\sqrt{|\alpha-3\beta|}\,|\sqrt{\alpha}-\sqrt{3\beta}|}{3\sqrt{6\alpha\beta}\,a^{2}H}\,e^{ik(e^{-N}-1)/(a_{i}H_{i})}\,,$
(132)

where $N=\ln(a/a_{i})$ is the e-folding number, and we set $N=0$ at the
initial time. Substituting this solution into the relation between ${\cal Z}$
and ${\cal Y}$, $\dot{{\cal Y}}$, $\ddot{{\cal Y}}$, we find

${\cal
Z}=\frac{\sqrt{\alpha}+\sqrt{3\beta}}{\sqrt{\alpha}-\sqrt{3\beta}}\,{\cal
Y}-\frac{6i\sqrt{3\alpha\beta}}{(\sqrt{\alpha}-\sqrt{3\beta})^{2}}\frac{aH}{k}{\cal
Y}+{\cal O}\left(\frac{a^{2}H^{2}}{k^{2}}\right){\cal Y}\,,$ (133)

which is valid for the modes deep inside the Hubble radius777Alternatively, we
can write the Lagrangian density in the form $L_{s}=\mathcal{Y}\,({\rm d}/{\rm
d}t)(a^{3}\dot{\mathcal{Z}})+\cdots$ and vary the action with respect to
${\cal Y}$ and use the equation of motion for ${\cal Y}$ to express $L_{s}$
with respect to ${\cal Z}$ and its derivatives. In the high-$k$ regime, the
Lagrangian density is given by
$L_{s}\simeq\tilde{Q}\,[\ddot{\mathcal{Z}}^{2}-(2k^{2}/a^{2})\,\dot{\mathcal{Z}}^{2}+(k^{4}/a^{4})\,\mathcal{Z}^{2}]$,
where
$\tilde{Q}=27\alpha\beta\,a^{7}H^{2}/[(\alpha-3\beta)(\sqrt{\alpha}+\sqrt{3\beta})^{2}k^{4}]$.
The WKB solution to ${\cal Z}$ derived by this procedure is consistent with
the leading-order relation of Eq. (133)..

On using Eqs. (127), (128), and (133), the WKB solutions to ${\cal A}$ and
${\cal B}$ are

$\displaystyle{\cal A}$ $\displaystyle=$
$\displaystyle\frac{1}{6\sqrt{2}M_{\rm Pl}H}\left[\frac{1}{\sqrt{\beta}}({\cal
Y}-{\cal Z})+\sqrt{\frac{3}{\alpha}}({\cal Y}+{\cal
Z})\right]=\frac{\sqrt{6}ai}{2(\sqrt{\alpha}-\sqrt{3\beta})M_{\rm Pl}k}{\cal
Y}\left[1+{\cal O}\left(\frac{aH}{k}\right)\right]\,,$ (134)
$\displaystyle{\cal B}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{6}a^{2}}{4M_{\rm Pl}k^{2}\sqrt{\alpha}}\left({\cal
Y}+{\cal Z}\right)=\frac{\sqrt{6}a^{2}}{2(\sqrt{\alpha}-\sqrt{3\beta})M_{\rm
Pl}k^{2}}\,{\cal Y}\left[1+{\cal O}\left(\frac{aH}{k}\right)\right]\,,$ (135)

where we used the relation (133). Recall that ${\cal Y}$ is given by Eq.
(132). The ratio between the leading-order terms to ${\cal A}$ and ${\cal B}$
is ${\cal B}/{\cal A}=-ia/k$, so the amplitude $|H{\cal B}/{\cal A}|$ is of
order $aH/k\ll 1$ for sub-Hubble modes. We will use Eqs. (134) and (135) as
the initial conditions of perturbations for the modes deep inside the Hubble
radius.

###  VI.3 Numerical integration

In the flat gauge, we numerically integrate the closed fourth-order
differential equations of ${\cal A}$ and ${\cal B}$ together with those of
$\Psi$, $\Phi$, and $\zeta$. For this purpose, we introduce the following
quantities

$h=\frac{H}{H_{i}}\,,\qquad K=\frac{k}{aH_{i}}\,,\qquad\bar{\alpha}=\alpha
H_{i}^{2}\,,\qquad\bar{\beta}=\beta H_{i}^{2}\,,\qquad\chi=H_{i}{\cal B}\,,$
(136)

where $H_{i}$ is the Hubble parameter at the onset of integration (with scale
factor $a_{i}$). We also define the following dimensionless variables

${\cal A}_{k}=\frac{k^{3/2}}{\sqrt{2\pi^{2}}}{\cal A}\,,\qquad{\cal
B}_{k}=\frac{k^{3/2}}{\sqrt{2\pi^{2}}}{\cal B}\,,\qquad\chi_{k}=H_{i}{\cal
B}_{k}\,,\qquad\Psi_{k}=\frac{k^{3/2}}{\sqrt{2\pi^{2}}}\Psi\,,\qquad\Phi_{k}=\frac{k^{3/2}}{\sqrt{2\pi^{2}}}\Phi\,,\qquad\zeta_{k}=\frac{k^{3/2}}{\sqrt{2\pi^{2}}}\zeta\,.$
(137)

Then, the two gravitational potentials can be expressed as

$\Psi_{k}={\cal
A}_{k}+h\chi_{k}^{\prime}\,,\qquad\mathrm{and}\qquad\Phi_{k}=h\chi_{k}\,,$
(138)

where a prime in this section denotes the differentiation with respect to the
number of e-foldings $N=\ln(a/a_{i})$. Together with solving the perturbation
equations of motion, we integrate the following background equations of motion

$h^{\prime}=-h\epsilon\,,\qquad\epsilon^{\prime}=\frac{3}{2}\epsilon^{2}-3\epsilon+\frac{1}{12\bar{\beta}h^{2}}\,.$
(139)

The initial value of $\epsilon$ is chosen to realize the sufficient number of
e-foldings ($N>70$) during inflation. On using Eqs. (134) and (135), we choose
the initial conditions of ${\cal A}_{k}$ and $\chi_{k}$ at $N=0$, as888All the
perturbations labeled by $k$, for instance, $\chi_{k}$, $\zeta_{k}$, etc.,
satisfy the same closed fourth-order differential equation as their
unnormalized counterparts, $\chi$, $\zeta$, etc., with different initial
conditions only by the factor $k^{3/2}/\sqrt{2\pi^{2}}$.

$\displaystyle{\cal
A}_{k}(0)=iK_{i}\frac{|\bar{\alpha}-3\bar{\beta}|^{1/2}\,{\rm
sign}\bigl{(}\sqrt{\bar{\alpha}}-\sqrt{3\bar{\beta}}\bigr{)}}{6\sqrt{2}\pi\sqrt{\bar{\alpha}\bar{\beta}}}\frac{H_{i}}{M_{\rm
Pl}}\,,\qquad\frac{{\rm d}^{n}{\cal A}_{k}}{{\rm
d}N^{n}}(0)\simeq(-iK_{i})^{n}\,{\cal A}_{k}(0)\,,$ (140)
$\displaystyle\chi_{k}(0)=\frac{|\bar{\alpha}-3\bar{\beta}|^{1/2}\,{\rm
sign}\bigl{(}\sqrt{\bar{\alpha}}-\sqrt{3\bar{\beta}}\bigr{)}}{6\sqrt{2}\pi\sqrt{\bar{\alpha}\bar{\beta}}}\frac{H_{i}}{M_{\rm
Pl}}\,,\qquad\frac{{\rm d}^{n}\chi_{k}}{{\rm
d}N^{n}}(0)\simeq(-iK_{i})^{n}\chi_{k}(0)\,.$ (141)

where $K_{i}=K(0)=k/(a_{i}H_{i})$. Notice that there is a simple relation
${\cal A}_{k}(0)=iK_{i}\,\chi_{k}(0)$ for the leading-order solution. The
ratio $H_{i}/M_{\rm Pl}$ and the couplings $\alpha$, $\beta$ determine the
initial amplitude of $\chi_{k}$. The typical Hubble scale for Starobinsky
inflation is $H_{i}/M_{\rm Pl}={\cal O}(10^{-5})$, so that
$|\chi_{k}(N=0)|={\cal O}(10^{-6})$ for the couplings $\alpha$ and $\beta$
whose orders are similar to each other. The initial conditions of
gravitational potentials are $\Psi_{k}(0)={\cal
A}_{k}(0)+\chi_{k}^{\prime}(0)$ and $\Phi_{k}(0)=\chi_{k}(0)$.

From Eq. (110), the curvature perturbation $\zeta_{k}$ can be expressed, in
terms of $\mathcal{A}_{k}$ and $\chi_{k}$, as

$\zeta_{k}=\frac{4\bar{\beta}\left[h\left(2\chi_{k}+\chi_{k}^{\prime}\right)+\mathcal{A}_{k}\right]}{1-6\bar{\beta}\epsilon\left(\epsilon-2\right)h^{2}}K^{2}-\frac{12\bar{\beta}h^{2}\left[2\mathcal{A}_{k}\,(2-\epsilon)+\mathcal{A}_{k}^{\prime}\right]}{1-6\bar{\beta}\epsilon\left(\epsilon-2\right)h^{2}}\,.$
(142)

From Eq. (140), we find that ${\cal A}_{k}(0)$ is proportional to
$K_{i}H_{i}/M_{\rm Pl}$ and hence $\zeta_{k}(0)$ contains a large term
$K_{i}^{3}(H_{i}/M_{\rm Pl})$ for the sub-Hubble modes $K_{i}\gg 1$. Because
of the relation ${\cal A}_{k}(0)=iK_{i}\,\chi_{k}(0)$, this term exactly
cancels the other contribution $\chi_{k}^{\prime}(0)K_{i}^{2}$ for the
leading-order initial condition $\chi_{k}^{\prime}(0)=-iK_{i}\chi_{k}(0)$.
This cancellation implies that we need to take into account terms of order
${\cal O}(a^{2}H^{2}/k^{2})$ in Eq. (133) to estimate the initial value of
$\zeta_{k}(0)$ correctly. We also note that Eq. (142) contains the $N$
derivatives of ${\cal A}_{k}$ and $\chi_{k}$. Since the amplitudes of ${\cal
A}_{k}$ and $\chi_{k}$ change in time, we take the $N$ derivatives of these
fields without neglecting their time dependence and incorporate them as
corrections to the seconds of Eqs. (140) and (141). These precise
manipulations show that the leading-order contributions to the first and
second terms in Eq. (142) cancel each other with respect to the large $K_{i}$
expansion. Then, the initial value of $\zeta_{k}(0)$ is typically of order
$\bar{\beta}\chi_{k}(0)K_{i}$. After deriving the precise numerical value of
$\zeta_{k}(0)$ from Eq. (142) without using the approximation, the $N$
derivatives of $\zeta_{k}$ at $N=0$ can be estimated as $({\rm
d}^{n}\zeta_{k}/{\rm d}N^{n})(0)\simeq(-iK_{i})^{n}\,\zeta(0)$.

Figure 1: Exponential increase of $|\Phi_{k}|=|h\chi_{k}|$ and $|\Psi_{k}|$
for the couplings $\bar{\alpha}=1$ and $\bar{\beta}=6$. In the regime $k\ll
aH$, the growth of both $|\Phi_{k}|$ and $|\Psi_{k}|$ can be well-fitted by an
exponential function proportional to $e^{N/2}$. The initial conditions are
chosen to be $h(0)=1$, $\epsilon(0)=0.0047$, $H_{i}/M_{\rm Pl}=10^{-5}$, and
$\Phi_{k}(0)=\chi_{k}(0)\simeq-6.31\times 10^{-7}$ (which is determined by Eq.
(141)) for the sub-Hubble mode $K(0)=100$. We also show two exponential
functions proportional to $e^{N/2}$ as a dashed pink line and a dash-dotted
light blue line, which fit well with the numerical solutions of $|\Phi_{k}|$
and $|\Psi_{k}|$, respectively. The initial conditions of ${\cal A}_{k}$ and
$\zeta_{k}$ are known from Eqs. (140) and (142), as ${\cal
A}_{k}(0)\simeq-6.34\times 10^{-7}-6.32\times 10^{-5}i$ and
$\zeta_{k}(0)\simeq 1.3\times 10^{-4}-5.9\times 10^{-4}i$. To show the
evolution of $\zeta_{k}$, we have solved its own closed fourth-order
differential equation. We see that $|\zeta_{k}|$ approaches a constant after
the Hubble radius crossing, while ${\cal A}_{k}$ exhibits very mild growth.

As we will show in the numerical calculation below, the perturbation
$\chi_{k}$ is unstable, but the decrease of the $K^{2}$ term (proportional to
$a^{-2}$) in Eq. (142) suppresses the growth of $\chi_{k}K^{2}$. Therefore,
after the Hubble radius crossing, $\zeta_{k}$ depends mostly on
$\mathcal{A}_{k}$ and its $N$-derivative. Hence, if the growth of ${\cal
A}_{k}$ is insignificant, this is also the case for $\zeta_{k}$. Furthermore,
this shows that, in the super-Hubble regime, $\zeta_{k}$ is related only to
${\cal A}_{k}$ and vice-versa. Thus, the whole scalar sector, which consists
of two independent dynamical DOFs, cannot be described by $\zeta_{k}$ and
${\cal A}_{k}$ alone in the regime $k/a\ll H$.

In Fig. 1, we plot the evolution of $|\Phi_{k}|=|h\chi_{k}|$, $|\Psi_{k}|$,
$|\mathcal{A}_{k}|$, and $|\zeta_{k}|$ for $\bar{\alpha}=1$, $\bar{\beta}=6$,
and $K(0)=100$ with the slow-roll parameter $\epsilon(0)=0.0047$. The initial
conditions for $\chi_{k}$ are instead chosen as to fulfill Eq. (141), with
$H_{i}=10^{-5}M_{\rm Pl}$. We solved the fourth-order differential equations
of $\chi_{k}$ and ${\cal A}_{k}$ and computed $\Phi_{k}$, $\Psi_{k}$ and
$\zeta_{k}$ by exploiting the relations (138) and (142). We also performed the
direct integration of the fourth-order differential equations of $\Phi_{k}$,
$\Psi_{k}$, and $\zeta_{k}$ by implementing the initial conditions for each of
them in terms of those given for $\chi_{k}(0)$ and found that the results are
in perfect agreement with those computed from $\chi_{k}$ and ${\cal A}_{k}$.

In Fig. 1, we observe that both the amplitudes of $\Phi_{k}=H{\cal B}_{k}$ and
$\Psi_{k}$ grow in proportion to $e^{Ht/2}\simeq e^{N/2}$. We recall that we
used the approximation $\beta H^{2}\gg 1$ to derive the analytic solutions
(86), but the numerical results show that this estimation is valid even for
$\beta H^{2}$ of order 1. In Fig. 1, the gravitational potentials exceed the
order 1 around the e-foldings $N=30\sim 35$ after the onset of inflation.
Thus, the exponential growth of ${\cal B}$ in the perturbed line element (95)
invalidates the FLRW background. Due to the uncertainty principle, the initial
value of the perturbation $\chi_{k}$ has a nonvanishing value related to the
energy scale $H_{i}/M_{\rm Pl}$ during inflation. Since $H_{i}/M_{\rm Pl}$
should not be much smaller than $10^{-5}$, the gravitational potentials reach
order 1 after the amplification of $e^{N/2}$ times with $N>30$. As we
estimated analytically in Sec. V, the perturbation $\zeta_{k}$ approaches a
constant after the Hubble radius crossing, while ${\cal A}_{k}$ shows very
mild growth. We recall that $\zeta_{k}$ and ${\cal A}_{k}$ are related to each
other in the super-Hubble regime and that they are not sufficient to describe
the dynamics of scalar perturbations with two propagating DOFs. Indeed, we
cannot eliminate the instability of the other dynamical perturbation ${\cal
B}_{k}=\Phi_{k}/H$.

Figure 2: The same as Fig. 1, but for the Weyl coupling $\bar{\alpha}=1/400$,
$\chi_{k}(0)\simeq-1.3\times 10^{-5}$, $\mathcal{A}_{k}(0)\simeq-1.3\times
10^{-5}-0.0013i$, and $\zeta_{k}(0)\simeq 0.0028-2.86\times 10^{-5}i$. Also
for this small coupling constant $\bar{\alpha}$, the gravitational potentials
grow as $e^{N/2}$. Figure 3: The same as Fig. 1, except for
$\bar{\alpha}=400$, $\chi_{k}(0)\simeq 1.5\times 10^{-7}$,
$\mathcal{A}_{k}(0)\simeq 1.5\times 10^{-7}+1.5\times 10^{-5}i$,
$\zeta_{k}(0)\simeq-3.2\times 10^{-5}-0.0025i$, and for the Weyl coupling
$\bar{\alpha}=400$. For this large coupling constant $\bar{\alpha}$, the
gravitational potentials grow as fast as $e^{3N/4}$.

In Fig. 2, we show the evolution of gauge-invariant perturbations by keeping
the same initial conditions and model parameters as those in Fig. 1, except
for $\bar{\alpha}$ which is set to a smaller value ($\bar{\alpha}=1/400$) and
for the initial values of the fields $\chi_{k}$, ${\cal A}_{k}$, and
$\zeta_{k}$. Even with this small value of the Weyl coupling, both
$|\Phi_{k}|$ and $|\Psi_{k}|$ increase in proportion to $e^{Ht/2}\simeq
e^{N/2}$ after the Hubble radius crossing, while the growth of $|\zeta_{k}|$
and $|{\cal A}|$ is suppressed. Indeed, this behavior is expected according to
the analytic estimations of $\Phi_{k}=H{\cal B}_{k}$, $\Psi_{k}$, ${\cal
A}_{k}$, and $\zeta_{k}$. For small values of $\alpha$, the last two terms in
Eqs. (84) and (85) grow in proportion to $e^{Ht/2}$ with oscillations999In
Fig. 2, the oscillations of $|\Psi_{k}|$ and $|\Phi_{k}|$ are not clearly seen
at large $N$, but we confirmed that they are present by enlarging the figure..
Note that Starobinsky inflation without the Weyl term ($\alpha=0$) cannot be
recovered by simply taking the limit $\alpha\to 0$ in our theory. For
$\alpha=0$ there is only a single scalar DOF arising from the $\beta R^{2}$
term, in which case the exponential growth of gravitational potentials is
absent. However, as we saw for the vector and tensor perturbations, the
squared mass of the extra modes, when the Weyl-squared term is present, is
typically of order $\alpha^{-1}$. Then, in the limit $\alpha\to 0$, the extra
modes become very massive and we should expect them to acquire a mass larger
than an ultraviolet cutoff scale of order $M_{\rm Pl}$. In this case, they may
be integrated out from the theory. On the other hand, the instability of
scalar perturbations persists for the Weyl coupling constant in the range
$\alpha\gtrsim M_{\rm Pl}^{-2}$.

In Fig. 3, the evolution of gauge-invariant perturbations is plotted for the
large Weyl coupling $\bar{\alpha}=400$ with $\bar{\beta}=6$. In this case, the
condition $\alpha>48\beta$ is satisfied and hence the two gravitational
potentials grow as Eq. (87) without oscillations. Indeed, our numerical
results in Fig. 3 demonstrate that, after the Hubble radius crossing, the
growth of $|\Phi_{k}|$ and $|\Psi_{k}|$ occurs faster than in the case
$\bar{\alpha}=1$. The suppressed growth of the other perturbations ${\cal
A}_{k}$ and $\zeta_{k}$ also agrees with the analytic estimation in the regime
$\beta H^{2}\gg 1$.

In general, out of two exponentially growing modes $\Phi_{k}$ and $\Psi_{k}$,
one can find some linear combinations of them, like $\zeta_{k}$, whose growth
is suppressed in the super-Hubble regime. In this theory, however, two scalar
perturbations determine the stability of the background and not only one.
Therefore, the stability of one linear combination is not sufficient for
guaranteeing the stability of the whole dynamical system. In the flat gauge,
the instability of the field ${\cal B}$ appearing in the perturbed metric (95)
is enough to make the whole background unstable in the super-Hubble regime. We
also solved the perturbation equations in the Newtonian gauge and obtained the
same numerical solutions for the gauge-invariant fields $\Psi_{k}$,
$\Phi_{k}$, $\zeta_{k}$, and ${\cal A}_{k}$ as those in the flat gauge. In the
Newtonian gauge, the violation of the FLRW background occurs by the growth of
two gravitational potentials $\Psi$ and $\Phi$ in the perturbed line element
(66). In the unitary gauge, we numerically confirmed the same exponential
growth of $\Phi$ and the constancy of $\zeta$ after the Hubble radius
crossing. In this case, the gauge-invariant combination ${\cal
B}_{u}=(\Phi-\zeta)/H$ appearing in the perturbed metric (115) is subject to
the exponential growth, thereby invalidating the FLRW background.

##  VII Conclusions

In this paper, we studied the dynamics of cosmological perturbations during
inflation in quadratic gravity containing the Weyl term $-\alpha C^{2}$
besides the Ricci squared term $\beta R^{2}$ in the action. Although the Weyl
curvature does not affect the background inflationary dynamics driven by the
$\beta R^{2}$ term, the evolution of perturbations is modified by the presence
of derivatives higher than second order. Since these higher-order derivatives
can give rise to ghosts, it is of interest to explore whether or not the
ghosts can lead to instabilities of the FLRW background.

As we discussed in Sec. II, geometric inflation is realized by the $\beta
R^{2}$ term with $\beta>0$, where the coupling constant $\beta$ is related to
the mass squared $m_{S}^{2}$ of a new scalar DOF (scalaron) as
$\beta=1/(6m_{S}^{2})$. To realize the number of e-foldings larger than 60, we
require that the Hubble parameter $H_{i}$ at the onset of inflation is in the
range $\beta H_{i}^{2}\gtrsim{\cal O}(1)$. If we transform the action (1) to
that in the Einstein frame, the quadratic gravity can be interpreted as the
conformally invariant Weyl theory in the presence of a canonical scalaron
field with the potential. Unlike the past related works [41, 42, 43, 44, 45],
we have carried out all the analysis in the physical Jordan frame.

In Sec. III, we showed that the Weyl term gives rise to two dynamical vector
DOFs propagating with the speed of light. For the Weyl coupling $\alpha>0$ the
two ghosts are present with the positive mass squared $m_{W}^{2}=1/(2\alpha)$,
while, for $\alpha<0$, there are no ghosts. In the latter case, however, the
negative value of $m_{W}^{2}$ leads to the tachyonic instability of vector
perturbations for $|\alpha|$ at most of order $\beta$. To avoid such an
instability which violates the inflationary FLRW background, we demand the
condition $\alpha>0$ at the expense of admitting the existence of ghosts.

In Sec. IV, we derived the second-order action of tensor perturbations and
introduced Lagrange multiplier fields $\chi_{i}$ ($i=1,2$) associated with
higher-order time derivatives. There are four dynamical DOFs in the tensor
sector, two of which behave as ghosts. Using the WKB approximation for the
modes deep inside the Hubble radius ($k/a\gg H$), the speed of tensor
perturbations is equivalent to 1 and hence the Laplacian instability is
absent. In the super-Hubble regime ($k/a\ll H$), tensor perturbations $h_{i}$
obey the fourth-order differential Eq. (43). Provided that the couplings
$\alpha$ and $\beta$ are in the ranges $\alpha>0$ and $\beta H^{2}\gtrsim 1$,
we showed that $h_{i}$’s approach constants after the Hubble radius crossing.
This means that, despite the presence of ghosts, tensor perturbations are
subject to neither Laplacian nor tachyonic instabilities.

In Sec. V, we studied the stability and evolution of scalar perturbations by
choosing several different gauge conditions. There are two dynamical
propagating DOFs in the scalar sector arising from the Lagrangians $-\alpha
C^{2}$ and $\beta R^{2}$. For $\alpha>0$ and $\beta>0$, the scalaron is not a
ghost, but the other dynamical mode behaves as a ghost. To study the dynamics
of perturbations, we also introduced several gauge-invariant perturbations
such as those defined in Eqs. (52)-(54). We chose the Newtonian, flat, and
unitary gauges and derived the closed differential equations for $\Psi$,
$\Phi$, ${\cal A}$, ${\cal B}=\Phi/H$, and $\zeta$. We found that the
coefficients of these differential equations are uniquely fixed independent of
the gauge choices. We showed that, after the Hubble radius crossing, both
$\Psi$ and $\Phi$ grow exponentially, while ${\cal A}$ and $\zeta$ approach
constants.

In the Newtonian gauge given by the perturbed line element (66), the
exponential growth of $\Psi$ and $\Phi$ occurs in the $g_{00}$ and $g_{ii}$
metric components. This violates the stability of the FLRW background after
the perturbations cross the Hubble radius during inflation. For the flat-gauge
line element (95) the growth of ${\cal A}$ is suppressed, but the exponential
increase of ${\cal B}=\Phi/H$ occurs accompanied by the enhancement of $\Phi$.
We have also numerically confirmed this behavior for the gauge-invariant
perturbations in the numerical simulations of Figs. 1-3 performed in Sec. VI.
In the unitary gauge, the perturbed line element (115) also contains the
instability mode $\Phi$ in the $g_{0i}$ metric component. We stress that these
instabilities are the physical ones arising from the gravitational interaction
between the scalaron and the other ghost DOF.

We have thus shown that the inflationary FLRW background realized by the
$\beta R^{2}$ term is violated by the presence of the Weyl term. In other
words, the Universe becomes highly inhomogeneous during inflation, being
incompatible with the observations of CMB temperature anisotropies. This
instability of scalar perturbations is present for the wide coupling range
$\alpha\gtrsim M_{\rm Pl}^{-2}$ in which the mass term $1/\sqrt{\alpha}$
associated with the Weyl term does not exceed the ultraviolet scale of order
$M_{\rm Pl}$. Unless the scalar ghost arising from the Weyl term is suitably
eliminated as a physical propagating DOF and the classical instability of the
background is removed, the quadratic curvature theory with $\alpha\neq 0$ is
excluded as a viable model of inflation (or at most, the coupling $\alpha$
must be so small that the mass of the extra modes becomes larger than the
cutoff of the theory). Related to the ghost issue, there is an approach of
“fakeon” where the ghost does not appear as a physical state after quantizing
it as a fake DOF [62, 63, 64]. There are also some approaches to the ghost
problem in quantum field theory by keeping its physical status intact [65, 66,
67, 68]. In such approaches to the ghost problem, it will be of interest to
study the stability of cosmological perturbations and resulting observational
consequences in detail.

## Acknowledgements

The work of ADF was supported by the Japan Society for the Promotion of
Science Grants-in-Aid for Scientific Research No. 20K03969 and by grant
PID2020-118159GB-C41 funded by MCIN/AEI/10.13039/501100011033. ST was
supported by the Grant-in-Aid for Scientific Research Fund of the JSPS No.
22K03642 and Waseda University Special Research Project No. 2023C-473.

## References

  * Starobinsky [1980] A. A. Starobinsky,  _Phys. Lett. B_ 91, 99 (1980).
  * Sato [1981] K. Sato,  _Mon. Not. Roy. Astron. Soc._ 195, 467 (1981).
  * Kazanas [1980] D. Kazanas,  _Astrophys. J. Lett._ 241, L59 (1980).
  * Guth [1981] A. H. Guth,  _Phys. Rev. D_ 23, 347 (1981).
  * Starobinsky [1979] A. A. Starobinsky,  _JETP Lett._ 30, 682 (1979).
  * Starobinsky [1982] A. A. Starobinsky,  _Phys. Lett. B_ 117, 175 (1982).
  * Hawking [1982] S. W. Hawking,  _Phys. Lett. B_ 115, 295 (1982).
  * Guth and Pi [1982] A. H. Guth and S. Y. Pi,  _Phys. Rev. Lett._ 49, 1110 (1982).
  * Bardeen  _et al._ [1983] J. M. Bardeen, P. J. Steinhardt, and M. S. Turner,  _Phys. Rev. D_ 28, 679 (1983).
  * Spergel  _et al._ [2003] D. N. Spergel _et al._ (WMAP),  _Astrophys. J. Suppl._ 148, 175 (2003), arXiv:astro-ph/0302209 .
  * Ade  _et al._ [2014] P. A. R. Ade _et al._ (Planck),  _Astron. Astrophys._ 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Akrami  _et al._ [2020] Y. Akrami _et al._ (Planck),  _Astron. Astrophys._ 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Starobinsky [1985] A. A. Starobinsky,  _JETP Lett._ 42, 152 (1985).
  * Vilenkin [1985] A. Vilenkin,  _Phys. Rev. D_ 32, 2511 (1985).
  * Mijic  _et al._ [1986] M. B. Mijic, M. S. Morris, and W.-M. Suen,  _Phys. Rev. D_ 34, 2934 (1986).
  * Tsujikawa  _et al._ [1999] S. Tsujikawa, K.-i. Maeda, and T. Torii,  _Phys. Rev. D_ 60, 123505 (1999), arXiv:hep-ph/9906501 .
  * Starobinsky  _et al._ [2001] A. A. Starobinsky, S. Tsujikawa, and J. Yokoyama,  _Nucl. Phys. B_ 610, 383 (2001), arXiv:astro-ph/0107555 .
  * Olmo [2005] G. J. Olmo,  _Phys. Rev. D_ 72, 083505 (2005), arXiv:gr-qc/0505135 .
  * Faraoni [2005] V. Faraoni,  _Phys. Rev. D_ 72, 124005 (2005), arXiv:gr-qc/0511094 .
  * Chiba  _et al._ [2007] T. Chiba, T. L. Smith, and A. L. Erickcek,  _Phys. Rev. D_ 75, 124014 (2007), arXiv:astro-ph/0611867 .
  * Navarro and Van Acoleyen [2005] I. Navarro and K. Van Acoleyen,  _Phys. Lett. B_ 622, 1 (2005), arXiv:gr-qc/0506096 .
  * Amendola and Tsujikawa [2008] L. Amendola and S. Tsujikawa,  _Phys. Lett. B_ 660, 125 (2008), arXiv:0705.0396 [astro-ph] .
  * Brans and Dicke [1961] C. Brans and R. H. Dicke,  _Phys. Rev._ 124, 925 (1961).
  * O’Hanlon [1972] J. O’Hanlon,  _Phys. Rev. Lett._ 29, 137 (1972).
  * Chiba [2003] T. Chiba,  _Phys. Lett. B_ 575, 1 (2003), arXiv:astro-ph/0307338 .
  * Sotiriou and Faraoni [2010] T. P. Sotiriou and V. Faraoni,  _Rev. Mod. Phys._ 82, 451 (2010), arXiv:0805.1726 [gr-qc] .
  * Mukhanov  _et al._ [1987] V. F. Mukhanov, L. A. Kofman, and D. Y. Pogosian,  _Phys. Lett. B_ 193, 427 (1987).
  * Salopek  _et al._ [1989] D. S. Salopek, J. R. Bond, and J. M. Bardeen,  _Phys. Rev. D_ 40, 1753 (1989).
  * Hwang and Noh [2001] J.-c. Hwang and H. Noh,  _Phys. Lett. B_ 506, 13 (2001), arXiv:astro-ph/0102423 .
  * Mukhanov and Chibisov [1981] V. F. Mukhanov and G. V. Chibisov,  _JETP Lett._ 33, 532 (1981).
  * De Felice and Tsujikawa [2010] A. De Felice and S. Tsujikawa,  _Living Rev. Rel._ 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Salvio [2018] A. Salvio,  _Front. in Phys._ 6, 77 (2018), arXiv:1804.09944 [hep-th] .
  * Lovelock [1971] D. Lovelock,  _J. Math. Phys._ 12, 498 (1971).
  * Stelle [1977] K. S. Stelle,  _Phys. Rev. D_ 16, 953 (1977).
  * Fradkin and Tseytlin [1982] E. S. Fradkin and A. A. Tseytlin,  _Nucl. Phys. B_ 201, 469 (1982).
  * Stelle [1978] K. S. Stelle,  _Gen. Rel. Grav._ 9, 353 (1978).
  * Hindawi  _et al._ [1996] A. Hindawi, B. A. Ovrut, and D. Waldram,  _Phys. Rev. D_ 53, 5583 (1996), arXiv:hep-th/9509142 .
  * Bogdanos  _et al._ [2010] C. Bogdanos, S. Capozziello, M. De Laurentis, and S. Nesseris,  _Astropart. Phys._ 34, 236 (2010), arXiv:0911.3094 [gr-qc] .
  * Hinterbichler and Saravani [2016] K. Hinterbichler and M. Saravani,  _Phys. Rev. D_ 93, 065006 (2016), arXiv:1508.02401 [hep-th] .
  * Smilga [2005] A. V. Smilga,  _Nucl. Phys. B_ 706, 598 (2005), arXiv:hep-th/0407231 .
  * Clunan and Sasaki [2010] T. Clunan and M. Sasaki,  _Class. Quant. Grav._ 27, 165014 (2010), arXiv:0907.3868 [hep-th] .
  * Deruelle  _et al._ [2011] N. Deruelle, M. Sasaki, Y. Sendouda, and A. Youssef,  _JCAP_ 03, 040 (2011), arXiv:1012.5202 [gr-qc] .
  * Myung and Moon [2015] Y. S. Myung and T. Moon,  _Mod. Phys. Lett. A_ 30, 1550172 (2015), arXiv:1501.01749 [gr-qc] .
  * Ivanov and Tokareva [2016] M. M. Ivanov and A. A. Tokareva,  _JCAP_ 12, 018 (2016), arXiv:1610.05330 [hep-th] .
  * Anselmi  _et al._ [2021] D. Anselmi, F. Fruzza, and M. Piva,  _Class. Quant. Grav._ 38, 225011 (2021), arXiv:2103.01653 [hep-th] .
  * Whitt [1984] B. Whitt,  _Phys. Lett. B_ 145, 176 (1984).
  * Barrow and Cotsakis [1988] J. D. Barrow and S. Cotsakis,  _Phys. Lett. B_ 214, 515 (1988).
  * Maeda [1989] K.-i. Maeda,  _Phys. Rev. D_ 39, 3159 (1989).
  * Gottlober  _et al._ [1990] S. Gottlober, H. J. Schmidt, and A. A. Starobinsky,  _Class. Quant. Grav._ 7, 893 (1990).
  * Bardeen [1980] J. M. Bardeen,  _Phys. Rev. D_ 22, 1882 (1980).
  * Kodama and Sasaki [1984] H. Kodama and M. Sasaki,  _Prog. Theor. Phys. Suppl._ 78, 1 (1984).
  * Mukhanov  _et al._ [1992] V. F. Mukhanov, H. A. Feldman, and R. H. Brandenberger,  _Phys. Rept._ 215, 203 (1992).
  * Bassett  _et al._ [2006] B. A. Bassett, S. Tsujikawa, and D. Wands,  _Rev. Mod. Phys._ 78, 537 (2006), arXiv:astro-ph/0507632 .
  * Lukash [1980] V. N. Lukash,  _Sov. Phys. JETP_ 52, 807 (1980).
  * Lyth [1985] D. H. Lyth,  _Phys. Rev. D_ 31, 1792 (1985).
  * Mukhanov [1985] V. F. Mukhanov,  _JETP Lett._ 41, 493 (1985).
  * Sasaki [1986] M. Sasaki,  _Prog. Theor. Phys._ 76, 1036 (1986).
  * Hwang and Noh [2002] J.-c. Hwang and H.-r. Noh,  _Phys. Rev. D_ 65, 023512 (2002), arXiv:astro-ph/0102005 .
  * Heisenberg  _et al._ [2018] L. Heisenberg, R. Kase, and S. Tsujikawa,  _Phys. Rev. D_ 98, 123504 (2018), arXiv:1807.07202 [gr-qc] .
  * Kase and Tsujikawa [2019] R. Kase and S. Tsujikawa,  _Int. J. Mod. Phys. D_ 28, 1942005 (2019), arXiv:1809.08735 [gr-qc] .
  * De Felice and Tsujikawa [2023] A. De Felice and S. Tsujikawa, arXiv:2307.06490 [gr-qc] .
  * Anselmi [2017] D. Anselmi,  _JHEP_ 06, 086 (2017), arXiv:1704.07728 [hep-th] .
  * Anselmi [2018] D. Anselmi,  _JHEP_ 02, 141 (2018), arXiv:1801.00915 [hep-th] .
  * Anselmi  _et al._ [2020] D. Anselmi, E. Bianchi, and M. Piva,  _JHEP_ 07, 211 (2020), arXiv:2005.10293 [hep-th] .
  * Bender and Mannheim [2008] C. M. Bender and P. D. Mannheim,  _Phys. Rev. Lett._ 100, 110402 (2008), arXiv:0706.0207 [hep-th] .
  * Donoghue and Menezes [2019] J. F. Donoghue and G. Menezes,  _Phys. Rev. D_ 100, 105006 (2019), arXiv:1908.02416 [hep-th] .
  * Salvio [2021] A. Salvio,  _Int. J. Mod. Phys. A_ 36, 2130006 (2021), arXiv:2012.11608 [hep-th] .
  * Kubo  _et al._ [2022] J. Kubo, J. Kuntz, J. Rezacek, and P. Saake,  _JCAP_ 11, 049 (2022), arXiv:2207.14329 [astro-ph.CO] .