# Dark matter production in Weyl $R^{2}$ inflation

Qing-Yang Wanga Yong Tanga,b,c,d Yue-Liang Wua,b,c,e aUniversity of Chinese
Academy of Sciences (UCAS), Beijing 100049, China
bSchool of Fundamental Physics and Mathematical Sciences,
Hangzhou Institute for Advanced Study, UCAS, Hangzhou 310024, China
cInternational Center for Theoretical Physics Asia-Pacific, Beijing/Hangzhou,
China
dNational Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China
eInstitute of Theoretical Physics, Chinese Academy of Sciences, Beijing
100190, China

###### Abstract

Dark matter and inflation are two key elements to understand the origin of
cosmic structures in modern cosmology, and yet their exact physical models
remain largely uncertain. The Weyl scaling invariant theory of gravity may
provide a feasible scheme to solve these two puzzles jointly, which contains a
massive gauge boson playing the role of dark matter candidate, and allows the
quadratic scalar curvature term, namely $R^{2}$, to realize a viable
inflationary mechanism in agreement with current observations. We ponder on
the production of dark matter in the Weyl $R^{2}$ model, including the
contribution from the non-perturbative production due to the quantum
fluctuations from inflationary vacuum and perturbative ones from scattering.
We demonstrate that there are generally three parameter regions for viable
dark matter production: (1) If the reheating temperature is larger than
$10^{3}~{}\mathrm{GeV}$, the Weyl gauge boson as dark matter can be produced
abundantly with mass larger than the inflation scale $\sim
10^{13}~{}\mathrm{GeV}$. (2) Small mass region with $3\times
10^{-13}~{}\mathrm{GeV}$ for a higher reheating temperature. (3) Annihilation
channel becomes important in the case of higher reheating temperature, which
enables the Weyl gauge boson with mass up to $4\times 10^{16}~{}\mathrm{GeV}$
to be produced through freeze-in.

##  I Introduction

Dark matter is an effective paradigm to explain the various observations in
galactic rotation curve, cosmic large scale structures, and the gravitational
lensing of clusters Rubin and Ford (1970); Allen  _et al._ (2003); Refregier
(2003); Tyson  _et al._ (1998). Studies of the primordial nucleosynthesis
Fields  _et al._ (2020) and the power spectrum of cosmic microwave background
(CMB) Aghanim  _et al._ (2020) also provide strong evidence that there is
abundant missing matter whose present total amount is about five times that of
visible matter described by the standard model of particle physics (SM). As
for the physical essence of dark matter, however, it still remains obscure.
All experiments so far have not identified whether dark matter interacts with
SM particles except gravity Lin (2019); Schumann (2019); Billard  _et al._
(2021).

In this paper, we ponder on a new kind of vector dark matter Tang and Wu
(2020a, b) from the Weyl scaling invariant theory of gravity Wu (2016, 2018).
Weyl symmetry with various forms has wide applications in cosmological physics
Weyl (1929); Mannheim and Kazanas (1989); Romero  _et al._ (2012); Bars  _et
al._ (2014); Scholz (2018). It can naturally introduce a local scaling
symmetry associated with a gauge boson (called Weyl gauge boson or WGB), which
acquires a mass from the symmetry breaking. The consideration of Weyl gravity
is also motivated by the reasonableness of introducing a quadratic scalar
curvature term $R^{2}$ into theory Ghilencea (2019); Ferreira  _et al._
(2019); Tang and Wu (2020c). It can provide a successful inflationary
mechanism with observational conformity, albeit it has a different form under
the Weyl geometry Beltran Jimenez  _et al._ (2016); Oda (2020); Ghilencea
(2021); Cai  _et al._ (2021). Our previous study has discussed Weyl $R^{2}$
inflation Tang and Wu (2020c) and briefly the production of dark matter.
However, the intrinsic connection between the inflaton and WGB has not been
fully explored. The purposes of this paper are twofold. First, we use the
latest cosmological observations to update the preferred parameter space for
Weyl $R^{2}$ inflation. Second, we identify the viable dark matter mass range
by considering all the production channels.

There are several contributions for the production of WGB. If dark matter is
already massive during inflation, there will be an important non-perturbative
production mechanism Parker (1969, 1971); Ford (1987); Lyth and Roberts
(1998); Chung  _et al._ (2001); Baer  _et al._ (2015); Graham  _et al._
(2016); Ema  _et al._ (2018, 2019); Ahmed  _et al._ (2020); Chung  _et al._
(2019); Kolb and Long (2021) that is due to the quantum fluctuations in
rapidly expanding background Ema  _et al._ (2017, 2016); Li  _et al._ (2020,
2019); Herring and Boyanovsky (2020); Babichev  _et al._ (2019); Hashiba and
Yokoyama (2019); Ling and Long (2021). Other important contributions are also
present after inflation, such as producing from inflaton decay or the
annihilation from particles in thermal bath after reheating Garny  _et al._
(2016); Tang and Wu (2016, 2017); Garny  _et al._ (2018); Chen and Kang
(2018); Bernal  _et al._ (2018); Aoki  _et al._ (2022); Clery  _et al._
(2021); Frangipane  _et al._ (2021); Redi and Tesi (2021); Barman  _et al._
(2022); Mambrini and Olive (2021); Haque and Maity (2021); Clery  _et al._
(2022). The previous study Tang and Wu (2020c) has only estimated the
annihilation channel through scalars, but not taken into account the
gravitational annihilation, non-perturbative production and the decay process.
Our goal here is to figure out how these production mechanisms determine the
relic abundance of WGB as dark matter, and identify the possible mass ranges.

The paper is organized as follows. In Sec. II, we introduce the formalism of
Weyl $R^{2}$ gravity with inflaton and massive WGB. Then in Sec. III, we
discuss the constraints on the Weyl $R^{2}$ inflation from the latest
observations. Later in Sec. IV, we focus on the WGB as dark matter. We first
discuss the equation of motion and quantization of the Weyl gauge field in a
cosmological background, and elaborate the method to calculate the energy
density of WGB produced from the inflationary fluctuations. Then based on the
inflation and reheating process, we compute the production rate and final
relic abundance of WGB from all possible channels, and show the constraints on
the mass and reheating temperature in detail. Finally, we summarize and give
our conclusion.

##  II Formalism and model

We use the following conventions in this paper: the metric
$\eta_{\mu\nu}=(-1,+1,+1,+1)$, and natural unit $\hbar=c=1$, $M_{P}\equiv
1/\sqrt{8\pi G}=2.435\times 10^{18}~{}\mathrm{GeV}=1$.

The complete Lagrangian that is motivated by gauge theory of gravity can be
found in Refs. Wu (2016, 2018, 2021a, 2021b). For the interest of dark matter
production and inflation, we can use the metric formalism and illustrate by
starting with the following simplest global scaling-invariant Lagrangian
$\mathcal{L}$ with linear and quadratic scalar curvature

$\frac{\mathcal{L}}{\sqrt{-g}}=\frac{1}{2}\phi^{2}R+\frac{\beta}{12}R^{2}-\frac{\lambda}{2}\partial^{\mu}\phi\partial_{\mu}\phi,$
(1)

where $g$ is the determinant of metric tensor $g_{\mu\nu}$, Ricci scalar $R$
are defined as usual and $\lambda$ is a nonzero real number. A scalar field
$\phi$ is introduced here to maintain the global scaling symmetry with the
following transformation rules,

$\displaystyle\mathrm{metric}:{}$ $\displaystyle g_{\mu\nu}\rightarrow
g^{\prime}_{\mu\nu}=f^{2}g_{\mu\nu},$ (2) $\displaystyle\mathrm{scalar}:{}$
$\displaystyle\phi\rightarrow\phi^{\prime}=f^{-1}\phi,$
$\displaystyle\mathrm{Ricci~{}scalar}:{}$ $\displaystyle R\rightarrow
R^{\prime}=f^{-2}R,$

where $f$ is a constant scale factor. Note that one can introduce a scalar
potential $\phi^{4}$ without spoiling the symmetry, but for our purpose in
this paper, such a term is not necessary and we neglect here.

For the global scaling transformation, the Christoffel connection

$\Gamma^{\rho}_{\mu\nu}=\frac{1}{2}g^{\rho\sigma}(\partial_{\mu}g_{\sigma\nu}+\partial_{\nu}g_{\mu\sigma}-\partial_{\sigma}g_{\mu\nu}),$
(3)

is invariant. When we extend the global symmetry into the local one, namely
$f\rightarrow f(x)$, the connection and $\mathcal{L}$ would not be invariant.
To make the theory locally scaling invariant, we construct the following
connection

$\hat{\Gamma}^{\rho}_{\mu\nu}=\Gamma^{\rho}_{\mu\nu}+(W_{\mu}g^{\rho}_{\nu}+W_{\nu}g^{\rho}_{\mu}-W^{\rho}g_{\mu\nu}),$
(4)

by introducing a vector $W_{\mu}\equiv g_{W}w_{\mu}$, called Weyl gauge field
($g_{W}$ is the gauge coupling), with the following transformation rule,

$W_{\mu}\rightarrow W^{\prime}_{\mu}=W_{\mu}-\partial_{\mu}\ln f(x).$ (5)

$\hat{\Gamma}^{\rho}_{\mu\nu}$ can be also obtained by replacing
$\partial_{\mu}g_{\rho\sigma}\rightarrow(\partial_{\mu}+2W_{\mu})g_{\rho\sigma}$
in Eq. 3. Then the corresponding Ricci scalar $\hat{R}$ should be defined
through $\hat{\Gamma}^{\rho}_{\mu\nu}$ and Lagrangian in Eq. (1) is modified
to a new form

$\frac{\mathcal{L}}{\sqrt{-g}}=\frac{1}{2}\phi^{2}\hat{R}+\frac{\beta}{12}\hat{R}^{2}-\frac{\lambda}{2}D^{\mu}\phi
D_{\mu}\phi-\frac{1}{4g_{W}^{2}}F_{\mu\nu}F^{\mu\nu},$ (6)

where the derivative in the kinetic term of $\phi$ has been replaced to the
covariant form, $D_{\mu}=\partial_{\mu}-W_{\mu}$, and the invariant field
strength of $W_{\mu}$ is defined as
$F_{\mu\nu}\equiv\partial_{\mu}W_{\nu}-\partial_{\nu}W_{\mu}$. This model is
also discussed in Ghilencea (2019); Ferreira  _et al._ (2019); Tang and Wu
(2020c).

To make the theoretical formalism more transparent for analyzing inflation and
dark matter, we shall utilize the equivalence between $f(R)$ theory and the
scalar-tensor theory Barrow and Cotsakis (1988); De Felice and Tsujikawa
(2010). We define

$F(\hat{R})=\frac{1}{2}\phi^{2}\hat{R}+\frac{\beta}{12}\hat{R}^{2},$ (7)

and introduce an auxiliary scalar field $\chi$. Then Eq. (6) can be written as
an equivalent form

$\displaystyle\frac{\mathcal{L}}{\sqrt{-g}}=$
$\displaystyle\frac{1}{2}[F^{\prime}(\chi^{2})(\hat{R}-\chi^{2})+F(\chi^{2})]-\frac{\lambda}{2}D^{\mu}\phi
D_{\mu}\phi-\frac{1}{4g_{W}^{2}}F_{\mu\nu}F^{\mu\nu}$ (8) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\phi^{2}+\frac{\beta}{3}\chi^{2}\right)\hat{R}-\frac{1}{12}\beta\chi^{4}-\frac{\lambda}{2}D^{\mu}\phi
D_{\mu}\phi-\frac{1}{4g_{W}^{2}}F_{\mu\nu}F^{\mu\nu},$

where $F^{\prime}(\chi^{2})\equiv\left|\frac{\partial
F(\hat{R})}{\partial\hat{R}}\right|_{\hat{R}=\chi^{2}}$. The equivalence
relation $\chi^{2}=\hat{R}$ can be easily derived by using the equation of
motion from $\delta\mathcal{L}/\delta\chi=0$.

Note that the local scaling symmetry of Eq. (8) allows us to choose convenient
$f(x)$ in Eq. (2) and (5) to fix the gauge and simplify the Lagrangian. To
compare with the Einstein gravity, we adopt
$f(x)=\sqrt{\phi^{2}+\frac{\beta}{3}\chi^{2}}$ and get
$\frac{1}{2}(\phi^{2}+\frac{\beta}{3}\chi^{2})\hat{R}=\frac{1}{2}\hat{R}$, or
equivalently, $\phi^{2}+\frac{\beta}{3}\chi^{2}=1$. Putting things together,
we can obtain the following Lagrangian

$\displaystyle\frac{\mathcal{L}}{\sqrt{-g}}=$
$\displaystyle\frac{1}{2}R-\frac{\lambda}{2}D^{\mu}\phi
D_{\mu}\phi-\frac{3}{4\beta}(1-\phi^{2})^{2}-\frac{1}{4g_{W}^{2}}F_{\mu\nu}F^{\mu\nu}-3W^{\mu}W_{\mu}$
(9) $\displaystyle=$
$\displaystyle\frac{1}{2}R-\frac{3\lambda}{6+\lambda\phi^{2}}\partial^{\mu}\phi\partial_{\mu}\phi-\frac{3}{4\beta}(1-\phi^{2})^{2}-\frac{1}{4g_{W}^{2}}F_{\mu\nu}F^{\mu\nu}$
$\displaystyle-\frac{1}{2}(6+\lambda\phi^{2})\left[W_{\mu}-\frac{1}{2}\partial_{\mu}\ln|6+\lambda\phi^{2}|\right]^{2}.$

where we have used the relation between $\hat{R}$ and $R$,

$\hat{R}=R-6W_{\mu}W^{\mu}-\frac{6}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}W^{\mu}),$
(10)

and neglect the surface term at infinite due to the total derivative term
above.

We define a new field variable $\Phi$,

$\Phi\equiv\pm\begin{cases}\sqrt{6}\ln\left[\sqrt{\frac{+\lambda\phi^{2}}{6}}+\sqrt{\frac{+\lambda\phi^{2}+6}{6}}\right]~{}\mathrm{for}~{}\lambda>0,\\\
\sqrt{6}\ln\left[\sqrt{\frac{-\lambda\phi^{2}}{6}}+\sqrt{\frac{-\lambda\phi^{2}-6}{6}}\right]~{}\mathrm{for}~{}\lambda<0,\end{cases}$
(11)

to get the canonical kinetic term, and make transformation for Weyl gauge
field

$\tilde{W}_{\mu}\equiv
W_{\mu}-\frac{1}{2}\partial_{\mu}\ln|6+\lambda\phi^{2}|\equiv
g_{W}\tilde{w}_{\mu}.$ (12)

Then the Lagrangian of Weyl $R^{2}$ model can be rewritten as a simple form

$\frac{\mathcal{L}}{\sqrt{-g}}=\frac{1}{2}R-\frac{1}{2}\partial^{\mu}\Phi\partial_{\mu}\Phi-\frac{1}{4g_{W}^{2}}\tilde{F}_{\mu\nu}\tilde{F}^{\mu\nu}-\frac{1}{2}m^{2}(\Phi)\tilde{W}^{\mu}\tilde{W}_{\mu}-V(\Phi),$
(13)

where we have defined the mass term

$m^{2}(\Phi)=\begin{cases}+6\cosh^{2}\left(\frac{\Phi}{\sqrt{6}}\right)~{}\mathrm{for}~{}\lambda>0,\\\
-6\sinh^{2}\left(\frac{\Phi}{\sqrt{6}}\right)~{}\mathrm{for}~{}\lambda<0,\\\
\end{cases}$ (14)

and scalar potential

$V(\Phi)=\begin{cases}\frac{3}{4\beta}\left[1-\frac{6}{\lambda}\sinh^{2}\left(\frac{\Phi}{\sqrt{6}}\right)\right]^{2}~{}\mathrm{for}~{}\lambda>0,\\\
\frac{3}{4\beta}\left[1+\frac{6}{\lambda}\cosh^{2}\left(\frac{\Phi}{\sqrt{6}}\right)\right]^{2}~{}\mathrm{for}~{}\lambda<0.\end{cases}$
(15)

Now it is clear that this theory can be reformulated as the Einstein gravity
with a massive Weyl gauge boson $\tilde{W}_{\mu}$ (WGB) and a self-interacting
scalar $\Phi$. In the following sections, we shall show that $\Phi$ can be the
inflaton field while $\tilde{W}_{\mu}$ is a dark matter candidate due to the
$Z_{2}$ symmetry, $\tilde{W}_{\mu}\rightarrow-\tilde{W}_{\mu}$.

##  III Inflation in Weyl $R^{2}$ Gravity

Now we concentrate on the inflationary mechanism in this model. The physics of
$R^{2}$ term has been extensively studied as the well-known Starobinsky model
Starobinsky (1980); Maeda (1988); Cembranos (2009); Bezrukov and Gorbunov
(2012); He  _et al._ (2018); Enckell  _et al._ (2019); Pi  _et al._ (2018);
Antoniadis  _et al._ (2018). It was proposed originally to solve the
cosmological initial singularity problem, and later demonstrated to be an
elegant mechanism to give rise to an inflationary universe. In this section,
we shall show the predictions in our model are different from those in
Starobinsky model.

Figure 1: The inflationary potentials of the Weyl $R^{2}$ model with
$\lambda=10^{2},10^{3},10^{4}$. It is a hilltop-like potential whose typical
vacumm expectation value is larger than the Planck scale.

The potential $V(\Phi)$ in Fig. 1 for a hilltop-like shape regardless of the
sign of $\lambda$. The parameter $\beta$ decides the height of the hilltop,
and $\lambda$ determines the position of valley. Fig. 1 shows the potentials
with $\beta=1.56\times 10^{10}$ in Planck unit and several values of
$\lambda>0$ ($\lambda<0$ is similar because $\cosh^{2}x=1+\sinh^{2}x$).

At the potential minimum $\Phi$ can outstrip the Planck scale when $|\lambda|$
is not too small. It is a general picture that for large $|\lambda|$ the
potential can be flat and lead to an inflationary universe when $\Phi$ slow-
rolls on the hilltop, until it falls into the valley.

Figure 2: The predictions of spectral index $n_{s}$ and tensor-to-scalar ratio
$r$ for the Weyl $R^{2}$ model with $\lambda$ from 300 to 3000, and e-folding
number $N=50,60$, in comparison with the recent constraints shown by the
BICEP/$Keck$ Collaboration Ade  _et al._ (2021).

To compare with the latest observational constraints by BICEP/$Keck$
Collaboration Ade  _et al._ (2021), we calculate the spectral index $n_{s}$
and tensor-to-scalar ratio $r$ for the primordial perturbations. These values
are related to the slow-roll parameters $\epsilon$ and $\eta$ through

$\displaystyle\epsilon=\frac{1}{2}\left[\frac{V^{\prime}(\Phi)}{V}\right]^{2},~{}\eta=\frac{V^{\prime\prime}(\Phi)}{V},~{}n_{s}=1-6\epsilon+2\eta,~{}r=16\epsilon.$
(16)

We numerically compute $n_{s}$ and $r$ for this model with $\lambda$ from 300
to 3000, and show the results in Fig. 2, where $N$ is defined as the e-folding
number of the inflationary multiples from the horizon scale as large as the
present one to the end of inflation Liddle and Leach (2003),
$N\equiv\ln(a_{e}/a_{i})\sim(50,60)$. The figure indicates that, as $\lambda$
increases, $n_{s}$ and $r$ will approach fixed values predicted in Starobinsky
model for every e-folds number $N$. The physical reasons has been explained in
detail in Ref. Tang and Wu (2020c) in which the authors show the equivalence
between this model and Starobinsky inflation as $|\lambda|\rightarrow+\infty$.
There is another noteworthy point in Fig. 2, the typical tensor-to-scalar
ratio in this model is greater than $\mathcal{O}(10^{-3})$, which implies that
the Weyl $R^{2}$ model may be tested by the next generation experiment of CMB
B-mode polarization Abazajian  _et al._ (2022). Current observation can give
constraints on $\lambda\gtrsim 500$ for $N=50$ and $\lambda\gtrsim 250$ for
$N=60$.

On the other hand, the observation of CMB perturbation also provide a
constraint on the Hubble parameter during inflation (denoted by
$H_{\mathrm{inf}}$) and the mass of inflaton $m_{\Phi}$ in this model. It is
known that the amplitude of scalar spectra $\Delta^{2}_{s}\sim
V/24\pi^{2}\epsilon$, and the observation of CMB has confirmed that
$\Delta^{2}_{s}\sim 2.1\times 10^{-9}$. Thus this result constrains the
parameter $\beta$ in Eq. (15), which determines the height of the potential
$V$ at the hilltop. In this way, we find the Hubble parameter during inflation
can be limited to $H_{\mathrm{inf}}\sim\sqrt{V/3}\sim 10^{13}~{}\mathrm{GeV}$.
Similarly, the mass of inflaton is also constrained to $m_{\Phi}\sim 2\times
10^{13}~{}\mathrm{GeV}$ defined at the minimum of the scalar potential.

##  IV Weyl gauge boson as dark matter

Now we discuss the physics of Weyl gauge boson. As we have shown, the
Lagrangian has $Z_{2}$ symmetry, $\tilde{W}_{\mu}\rightarrow-\tilde{W}_{\mu}$,
which enables the stability of $\tilde{W}_{\mu}$. The mass term
$-\frac{1}{2}m^{2}\tilde{W}^{\mu}\tilde{W}_{\mu}$ in Eq. (13), however, shows
that negative sign would appear for $\lambda<0$. This does not necessarily
mean an inconsistent theory because one could introduce another scalar field
$\psi$ with $D_{\mu}\psi D^{\mu}\psi$ that can contribute to the mass of
$\tilde{W}_{\mu}$ when $\psi\neq 0$ at its potential minimum. For
illustration, we shall only focus on the $\lambda>0$ case in the rest of this
paper.

Since the magnitude of $\Phi$ is evolving in the inflation and subsequent
reheating process, the mass term in Eq. (14) is also changing in this period,
until $\Phi$ completely stops at its potential minimum $\Phi_{0}$. We can
expand the mass term at $\Phi=\Phi_{0}$ as

$\displaystyle m^{2}=$ $\displaystyle
6\cosh^{2}\left(\frac{\varphi+\Phi_{0}}{\sqrt{6}}\right)$ (17)
$\displaystyle=$ $\displaystyle
6\cosh^{2}\left(\frac{\Phi_{0}}{\sqrt{6}}\right)+\sqrt{6}\sinh\left(\frac{2\Phi_{0}}{\sqrt{6}}\right)\varphi+\cosh^{2}\left(\frac{\Phi_{0}}{\sqrt{6}}\right)\varphi^{2}+\mathcal{O}(\varphi^{3}).$

It implies that the term $-\frac{1}{2}m^{2}\tilde{W}^{\mu}\tilde{W}_{\mu}$ not
only describes the mass of WGB, but also $\varphi-\tilde{W}$ interactions that
may lead to the decay or annihilation of inflaton into WGB after inflation.
The eventual mass of WGB is determined by the parameters $\lambda$ and $g_{W}$

$m_{W}=\sqrt{6}|g_{W}|\cosh[\Phi_{0}(\lambda)/\sqrt{6}].$ (18)

The couplings between $\varphi$ and WGB also depend on these two parameters
and can be read from above. In the rest of this section, we shall discuss the
production mechanism for WGB as a dark matter candidate in detail by including
all of the above couplings and interactions.

###  IV.1 Weyl gauge field in a cosmological background

The production of vector dark matter includes both contributions from
perturbative scattering and non-perturbative mechanism from vacuum
fluctuations due to the quantum effect in curved background Graham  _et al._
(2016); Ema  _et al._ (2018, 2019); Ahmed  _et al._ (2020). Let us concentrate
on the part of WGB in Eq. (13)

$S=\int
d^{4}x\sqrt{-g}\left[-\frac{1}{4g_{W}^{2}}\tilde{F}_{\mu\nu}\tilde{F}^{\mu\nu}-\frac{1}{2}m^{2}\tilde{W}^{\mu}\tilde{W}_{\mu}\right].$
(19)

The metric is set to the FLRW form

$g_{\mu\nu}dx^{\mu}dx^{\mu}=a^{2}(\tau)(-d\tau^{2}+d\vec{x}^{2}),$ (20)

where $a$ is the scale factor of the universe, and $\tau$ is the conformal
time, which has the following relation with the physical time $d\tau=dt/a$. We
can rewrite Eq. (19) in components

$S=\int
d^{3}xd\tau\frac{1}{2}\left(|\vec{\tilde{w}}^{\prime}-\nabla\tilde{w}_{0}|^{2}-|\nabla\times\vec{\tilde{w}}|^{2}+a^{2}\tilde{m}^{2}\tilde{w}_{0}^{2}-a^{2}\tilde{m}^{2}|\vec{\tilde{w}}|^{2}\right),$
(21)

where ${}^{\prime}\equiv\frac{\partial}{\partial\tau}$,
$\tilde{m}\equiv|g_{W}|m(\Phi)$, and $\vec{\tilde{w}}$ denotes the spatial
components.

It is obvious that $\tilde{w}_{0}$ in Eq. (21) has no kinetic term, which is
commonly known as that only three of the four components of massive vectors
are independent and physical modes. Hence, we can decouple the part of $W_{0}$
from Eq. (21) with Fourier transform, and only retain the physical modes,
including two transverse modes and a longitudinal mode

$\vec{\tilde{w}}(\tau,\vec{k})\equiv\sum_{j=\pm,L}\vec{\epsilon}_{j}(\vec{k})W_{j}(\tau,\vec{k}),$
(22)

where $\vec{\epsilon}_{j}(\vec{k})$ is the orthonormal polarization vector.
Their actions are derived as follows

$S_{\pm}=\int\frac{d^{3}kd\tau}{(2\pi)^{3}}\frac{1}{2}\left(|W_{\pm}^{\prime}|^{2}-(k^{2}+a^{2}\tilde{m}^{2})|W_{\pm}|^{2}\right),$
(23)
$S_{L}=\int\frac{d^{3}kd\tau}{(2\pi)^{3}}\frac{1}{2}\left(\frac{a^{2}\tilde{m}^{2}}{k^{2}+a^{2}\tilde{m}^{2}}|W_{L}^{\prime}|^{2}-a^{2}\tilde{m}^{2}|W_{L}|^{2}\right),$
(24)

where $k$ is the magnitude of comoving momentum $\vec{k}$, and
$\vec{k}\cdot\vec{\tilde{w}}=kW_{L}$ is used in the derivation.

For the sake of simplicity, we redefine a new field $\mathcal{W}_{L}$ for the
longitudinal mode

$\mathcal{W}_{L}\equiv\frac{a\tilde{m}}{\sqrt{k^{2}+a^{2}\tilde{m}^{2}}}W_{L}.$
(25)

Then the equation of motion for each physical mode can be written as
oscillator equations

$W_{\pm}^{\prime\prime}+\omega_{\pm}^{2}W_{\pm}=0,$ (26)
$\mathcal{W}_{L}^{\prime\prime}+\omega_{L}^{2}\mathcal{W}_{L}=0,$ (27)

with frequency variables

$\omega_{\pm}^{2}=k^{2}+a^{2}\tilde{m}^{2},$ (28)
$\displaystyle\omega_{L}^{2}=k^{2}+a^{2}\tilde{m}^{2}-\frac{k^{2}}{k^{2}+a^{2}\tilde{m}^{2}}\frac{a^{\prime\prime}}{a}+3\frac{k^{2}\tilde{m}^{2}a^{\prime
2}}{(k^{2}+a^{2}\tilde{m}^{2})^{2}}$ (29)
$\displaystyle-\frac{k^{2}}{k^{2}+a^{2}\tilde{m}^{2}}\frac{\tilde{m}^{\prime\prime}}{\tilde{m}}+3\frac{k^{2}a^{2}\tilde{m}^{\prime
2}}{(k^{2}+a^{2}\tilde{m}^{2})^{2}}-\frac{a^{\prime}}{a}\frac{\tilde{m}^{\prime}}{\tilde{m}}\frac{2k^{2}(k^{2}-2a^{2}\tilde{m}^{2})}{(k^{2}+a^{2}\tilde{m}^{2})^{2}}.$

The mode functions of WGB can be determined by solving Eq. (26) and (27).

Next, we show the relation between the energy density and the mode functions.
The energy-momentum tensor of field in a curved background is

$T_{\mu\nu}=\frac{-2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}S)}{\delta g^{\mu\nu}},$
(30)

and the energy density is defined as the vacuum expectation of $T_{00}$

$\rho=\langle 0|T_{00}|0\rangle.$ (31)

Therefore, to calculate the energy density of Weyl gauge field, we need to
define the vacuum state $|0\rangle$ of the theory. In other words, a quantized
theory is required. The first step of quantization is to give the Heisenberg
representation for Weyl gauge field

$\displaystyle\hat{W}_{\pm}(\tau,\vec{x})=\int\frac{d^{3}k}{(2\pi)^{3/2}}$
$\displaystyle\left[\hat{a}_{\vec{k}}\vec{\epsilon}_{\pm}(\vec{k})W_{\pm}(\tau,\vec{k})e^{i\vec{k}\cdot\vec{x}}+\hat{a}_{\vec{k}}^{\dagger}\vec{\epsilon}_{\pm}^{*}(\vec{k})W_{\pm}^{*}(\tau,\vec{k})e^{-i\vec{k}\cdot\vec{x}}\right],$
(32)
$\displaystyle\hat{\mathcal{W}}_{L}(\tau,\vec{x})=\int\frac{d^{3}k}{(2\pi)^{3/2}}$
$\displaystyle\left[\hat{b}_{\vec{k}}\vec{\epsilon}_{L}(\vec{k})\mathcal{W}_{L}(\tau,\vec{k})e^{i\vec{k}\cdot\vec{x}}+\hat{b}_{\vec{k}}^{\dagger}\vec{\epsilon}_{L}^{*}(\vec{k})\mathcal{W}_{L}^{*}(\tau,\vec{k})e^{-i\vec{k}\cdot\vec{x}}\right].$

Then the quantization demands that the ladder operators should satisfy the
following commutative relations

$\left[\hat{a}_{\vec{k}},\hat{a}_{\vec{k}}^{\dagger}\right]=\delta_{jj^{\prime}}\delta(\vec{k}-\vec{k}^{\prime}),~{}\left[\hat{b}_{\vec{k}},\hat{b}_{\vec{k}}^{\dagger}\right]=\delta(\vec{k}-\vec{k}^{\prime}),$
(33)

where $j$ is the index that refers to the two components of transverse mode.
The commutative relations can also be represented as a canonical form

$\displaystyle\left[\hat{W}_{\pm}(\tau,\vec{x}),\hat{W}_{\pm}^{\prime}(\tau,\vec{y})\right]=$
$\displaystyle i\delta(\vec{x}-\vec{y}),$ (34)
$\displaystyle\left[\hat{\mathcal{W}}_{L}(\tau,\vec{x}),\hat{\mathcal{W}}_{L}^{\prime}(\tau,\vec{y})\right]=$
$\displaystyle i\delta(\vec{x}-\vec{y}),$

which imply the following normalization condition

$W_{\pm}W_{\pm}^{\prime*}-W_{\pm}^{*}W_{\pm}^{\prime}=\mathcal{W}_{L}\mathcal{W}_{L}^{\prime*}-\mathcal{W}_{L}^{*}\mathcal{W}_{L}^{\prime}=i.$
(35)

Now we can define the vacuum $|0\rangle$ as a state that complies with
$\hat{a}_{\vec{k}}|0\rangle=\hat{b}_{\vec{k}}|0\rangle=0$. Then the non-
perturbative production of particles from vacuum can be calculated by
$\rho=\langle 0|T_{00}|0\rangle$, described by the time-dependent mode
functions. Finally, with Eq. (23), (24), and (30)-(33), the energy density of
Weyl gauge field in a cosmological background is obtained as

$\rho_{\pm}=\int\frac{k^{2}dk}{4\pi^{2}a^{4}}\left[|W^{\prime}_{\pm}|^{2}+(k^{2}+a^{2}\tilde{m}^{2})|W_{\pm}|^{2}-\omega_{\pm}\right],$
(36) $\displaystyle\rho_{L}=\int\frac{k^{2}dk}{4\pi^{2}a^{4}}\bigg{[}$
$\displaystyle|\mathcal{W}^{\prime}_{L}|^{2}+(k^{2}+a^{2}\tilde{m}^{2})|\mathcal{W}_{L}|^{2}+\bigg{(}\frac{(\frac{a^{\prime}}{a}+\frac{\tilde{m}^{\prime}}{\tilde{m}})k^{2}}{k^{2}+a^{2}\tilde{m}^{2}}\bigg{)}^{2}|\mathcal{W}_{L}|^{2}$
(37)
$\displaystyle-\frac{(\frac{a^{\prime}}{a}+\frac{\tilde{m}^{\prime}}{\tilde{m}})k^{2}}{k^{2}+a^{2}\tilde{m}^{2}}(\mathcal{W}^{\prime}_{L}\mathcal{W}^{*}_{L}+\mathcal{W}^{\prime*}_{L}\mathcal{W}_{L})-\omega_{L}\bigg{]},$

where the zero-point energy at flat spacetime has been subtracted by the term
$-\omega$. The total energy density is the sum of two transverse modes and one
longitudinal mode, $\rho_{tot}=2\rho_{\pm}+\rho_{L}$.

Next we numerically solve the equations of motion (26), (27) in the
inflationary background and calculate the energy density of WGB. Directly
solving these equations might be inaccurate in some cases. For instance, when
$\rho$ is many orders of magnitude smaller than the zero-point energy, the
numerical calculation of $\rho$ will meet the situation that subtracting two
very close numbers, which will make the numerical result unstable.
Fortunately, this intractable case usually has an adiabatic vacuum, which
means the vacuum evolves slowly. Then, a strategy called adiabatic
approximation can be used to deal with the problems Kofman  _et al._ (1997).
In the adiabatic approximation, the mode functions can be written as the
following form

$\displaystyle
W_{\pm}(\tau)=\mathcal{A}_{\pm}^{k}(\tau)v_{\pm}^{k}(\tau)+\mathcal{B}_{\pm}^{k}(\tau)v_{\pm}^{k*}(\tau),$
(38)
$\displaystyle\mathcal{W}_{L}(\tau)=\mathcal{A}_{L}^{k}(\tau)v_{L}^{k}(\tau)+\mathcal{B}_{L}^{k}(\tau)v_{L}^{k*}(\tau),$
$\displaystyle
v_{\pm(L)}^{k}=\frac{1}{\sqrt{2\omega_{\pm(L)}}}\exp\left[-i\int^{\tau}\omega_{\pm(L)}d\tau\right],$

where $\mathcal{A}^{k}(\tau)$ and $\mathcal{B}^{k}(\tau)$ are some functions
for the $k$-mode, and satisfy the normalization condition

$|\mathcal{A}^{k}|^{2}-|\mathcal{B}^{k}|^{2}=1.$ (39)

Substituting Eq. (38) into Eq. (32), we can redefine a new set of ladder
operators by the Bogolyubov transformation

$\tilde{a}_{\vec{k}}=\hat{a}_{\vec{k}}\mathcal{A}^{k}+\hat{a}^{\dagger}_{\vec{k}}\mathcal{B}^{k*},~{}\tilde{b}_{\vec{k}}=\hat{b}_{\vec{k}}\mathcal{A}^{k}+\hat{b}^{\dagger}_{\vec{k}}\mathcal{B}^{k*},$
(40)

which satisfy the similar commutative relations as Eq. (33)

$\left[\tilde{a}_{\vec{k}},\tilde{a}_{\vec{k}}^{\dagger}\right]=\delta_{jj^{\prime}}\delta(\vec{k}-\vec{k}^{\prime}),~{}\left[\tilde{b}_{\vec{k}},\tilde{b}_{\vec{k}}^{\dagger}\right]=\delta(\vec{k}-\vec{k}^{\prime}).$
(41)

Then the energy density of the Weyl gauge field can be written as

$\rho_{\pm/L}=\int\frac{k^{2}dk}{2\pi^{2}a^{4}}|\mathcal{B}^{k}_{\pm/L}|^{2}\omega_{\pm/L}.$
(42)

Comparing to Eq. (36) and (37), it can be seen the zero-point energy has
disappeared in this expression, which means the aforementioned numerical
problem is alleviated with the adiabatic approximation. In addition, the
function $\mathcal{B}^{k}$ obeys

$\mathcal{A}^{\prime}v=\frac{\omega_{k}^{\prime}}{2\omega_{k}}v^{*}\mathcal{B},~{}\mathcal{B}^{\prime}v^{*}=\frac{\omega_{k}^{\prime}}{2\omega_{k}}v\mathcal{A},$
(43)

which are identical to Eq. (26) and (27).

Finally, it is essential to note that the adiabatic approximation can only be
used to the case with positive $\omega^{2}$. This is the prerequisite for the
existence of adiabatic vacuum. For the transverse mode Eq. (28), this
condition is completely satisfied. But for the longitudinal mode Eq. (29), the
tachyonic situation ($\omega_{L}^{2}<0$) is possible, and it is common in the
case of $\tilde{m}\ll H_{\mathrm{inf}}$. Consequently, if the tachyonic
situation appears in the longitudinal mode, it would violate the adiabatic
condition. Then the mode functions and energy density can only be solved with
original formulas (Eq. (27) and (37)).

###  IV.2 Production of Weyl gauge boson

Now we are in a position to apply the above formalism in the cosmological
background of inflation. It can be seen from Eq. (26)-(29) that, to solve the
mode functions we also need the evolution of scale factor $a$ over conformal
time $\tau$, which is determined by the inflation and reheating in our model.

####  IV.2.1 Inflation and reheating

The evolution of cosmic scale $a$ is described by the Friedmann equation

$\mathcal{H}^{2}\equiv(\frac{a^{\prime}}{a})^{2}=\frac{a^{2}\rho}{3M_{P}^{2}},$
(44)

where ${}^{\prime}\equiv\frac{\partial}{\partial\tau}$. To solve this
equation, we need to discuss how the matter energy density evolving in our
model. To be more specific, we need reheating process to transfer the energy
of the inflaton $\Phi$ into other relativistic particles, for example through
decay into SM particles Kofman  _et al._ (1997); Albrecht  _et al._ (1982);
Kofman  _et al._ (1994); Bassett  _et al._ (2006); Lozanov (2019). Here,
instead we first consider a new scalar $\psi$ which couples to the inflaton
and SM particles simultaneously. Then in this way, the inflaton can first
decay into $\psi$ before eventually transfer its energy to the SM particles.
To maintain the scaling symmetry, we illustrate with the following toy model,

$\displaystyle\mathcal{L}_{\mathrm{reh}}=$ $\displaystyle-
g_{\psi}^{2}\phi^{2}\psi^{2}-\frac{1}{2}D^{\mu}\psi
D_{\mu}\psi+\mathcal{L}(\psi,SM)$ (45) $\displaystyle=$
$\displaystyle-\frac{6g_{\psi}^{2}}{\lambda}\sinh^{2}\left(\frac{\Phi}{\sqrt{6}}\right)\psi^{2}-\frac{1}{2}D^{\mu}\psi
D_{\mu}\psi+\mathcal{L}(\psi,SM)$ $\displaystyle\sim$
$\displaystyle-\frac{g_{\psi}^{2}}{\lambda}\Phi^{2}\psi^{2}-\frac{1}{2}D^{\mu}\psi
D_{\mu}\psi+\mathcal{L}(\psi,SM),$

where $\mathcal{L}(\psi,SM)$ denotes the coupling between $\psi$ and SM
particles, which we do not specify here. Here we have neglected higher order
interactions in the last step. We reiterate that the definition
$\Phi\equiv\varphi+\Phi_{0}$ has been introduced in the previous section. Then
the first term in Eq. (45) can be expanded as

$\displaystyle-\frac{g_{\psi}^{2}}{\lambda}\Phi^{2}\psi^{2}$
$\displaystyle=-\frac{g_{\psi}^{2}}{\lambda}(\varphi+\Phi_{0})^{2}\psi^{2}$
(46)
$\displaystyle=-\frac{g_{\psi}^{2}}{\lambda}\Phi_{0}^{2}\psi^{2}-2\frac{g_{\psi}^{2}}{\lambda}\Phi_{0}\varphi\psi^{2}-\frac{g_{\psi}^{2}}{\lambda}\varphi^{2}\psi^{2}.$

From left are the mass term of $\psi$, the three-scalar vertex, and the four-
scalar vertex. The feasibility of $\varphi$’s decay into $\psi$ requires the
mass of $\psi$ must be less than the half that of $\varphi$. Note that the
mass of inflaton is limited to $m_{\Phi}\sim 2\times 10^{13}~{}\mathrm{GeV}$
in our model, thus it is clear that $m_{\psi}\equiv
g_{\psi}\Phi_{0}/\sqrt{\lambda}<10^{13}~{}\mathrm{GeV}$. With a typical VEV
$\Phi_{0}\sim 10M_{P}$, the coupling is limited to
$g_{\psi}/\sqrt{\lambda}\lesssim 4\times 10^{-7}$, which gives an upper limit
to the reheating temperature $T_{r}$. Note that the major factor determining
$T_{r}$ is the $\varphi\rightarrow\psi\psi$ decay process caused by the three-
scalar vertex in Eq. (46), then the probability of inflaton decay is easily
derived as

$\Gamma(\varphi\rightarrow\psi\psi)=\frac{4g_{\psi}^{4}\Phi_{0}^{2}}{8\pi\lambda^{2}}\frac{\sqrt{(m_{\Phi}/2)^{2}-m_{\psi}^{2}}}{m_{\Phi}^{2}}=\frac{g_{\psi}^{2}m_{\psi}^{2}}{2\pi\lambda}\frac{\sqrt{(m_{\Phi}/2)^{2}-m_{\psi}^{2}}}{m_{\Phi}^{2}}\lesssim
10^{-20}M_{P}.$ (47)

As a consequence, the reheating temperature is estimated as $T_{r}\sim
0.4\sqrt{\Gamma M_{P}}\lesssim 10^{8}~{}\mathrm{GeV}$, which is far less than
the inflation energy scale.

With a specific reheating mechanism, the evolution of inflaton field obeys the
following equation

$\Phi^{\prime\prime}+2\mathcal{H}\Phi^{\prime}+a^{2}\frac{\partial
V}{\partial\Phi}+a\Gamma\Phi^{\prime}=0.$ (48)

Then we can write down the energy density of the inflaton

$\rho_{\Phi}=\frac{\Phi^{\prime 2}}{2a^{2}}+V(\Phi),$ (49)

and the evolving equation for $\psi$ and SM particles as radiation

$\rho_{r}^{\prime}+4\mathcal{H}\rho_{r}-\frac{\Gamma\Phi^{\prime 2}}{a}=0.$
(50)

Note that the total energy density in the Friedmann equation is the sum of
$\rho_{\Phi}$ and $\rho_{r}$. So combining Eq. (44) with (48)-(50) and (15),
we can obtain a numerical solution of $a(\tau)$, which describes the evolution
of cosmic scale within the period of inflation, reheating, and the ensuing
radiation-dominated epoch.

####  IV.2.2 Non-perturbative production

We take $\Gamma(\varphi\rightarrow\psi\psi)=10^{-20}M_{P}$ as an example in
this subsection, which is the maximum rate allowed in the toy model we
constructed above. We shall keep in mind that generally the reheating
temperature is free as long as it satisfies the bound from BBN and energy
density limit after inflation.

We solve the numerical result of $a(\tau)$, and substitute it into Eq.
(26)-(29), the mode functions for each eventual mass $m_{W}$ (defined as Eq.
(18) or $\tilde{m}|_{\Phi=\Phi_{0}}$) and co-moving momentum $k$ can be
derived with an initial condition. Because all the $k$-modes we are interested
in once had the wavelength much smaller than the horizon scale in the far past
of the inflation, it is then expected the term $k^{2}$ is dominant in Eq. (28)
and (29) under the limit of $\tau\rightarrow-\infty$. Therefore, the solution
to the equation of harmonic oscillator is suggested to be an initial condition
for both the transverse mode and the longitudinal mode, that is the well-known
Bunch-Davies vacuum state Bunch and Davies (1978)

$\lim_{\tau\rightarrow-\infty}W_{\pm}(\mathrm{or}~{}\mathcal{W}_{L})=\frac{1}{\sqrt{2k}}e^{-ik\tau}.$
(51)

With this initial condition, we illustrate with the eventual mass
$m_{W}=10^{10}~{}\mathrm{GeV}$ in Fig. 3, where $f_{k}$ is defined as a
function that is proportional to the energy density Eq. (36) and (37)

$\displaystyle f_{k,\pm}\equiv$
$\displaystyle\left[|W^{\prime}_{\pm}|^{2}+(k^{2}+a^{2}\tilde{m}^{2})|W_{\pm}|^{2}\right]/2\omega_{\pm}-\frac{1}{2}=|\mathcal{B}^{k}_{\pm}|^{2},$
(52) $\displaystyle f_{k,L}\equiv$
$\displaystyle\bigg{[}|\mathcal{W}^{\prime}_{L}|^{2}+(k^{2}+a^{2}\tilde{m}^{2})|\mathcal{W}_{L}|^{2}+\bigg{(}\frac{(\frac{a^{\prime}}{a}+\frac{\tilde{m}^{\prime}}{\tilde{m}})k^{2}}{k^{2}+a^{2}\tilde{m}^{2}}\bigg{)}^{2}|\mathcal{W}_{L}|^{2}$
$\displaystyle-\frac{(\frac{a^{\prime}}{a}+\frac{\tilde{m}^{\prime}}{\tilde{m}})k^{2}}{k^{2}+a^{2}\tilde{m}^{2}}(\mathcal{W}^{\prime}_{L}\mathcal{W}^{*}_{L}+\mathcal{W}^{\prime*}_{L}\mathcal{W}_{L})\bigg{]}/2\omega_{L}-\frac{1}{2}=|\mathcal{B}^{k}_{L}|^{2},$

and $a_{\mathrm{end}}$ denotes the scale factor when the inflation ends
(defined as the moment that the slow-roll parameter $\epsilon=1$). Moreover,
$k_{s}$ denotes the dominate mode of production, which is approximately equal
to the product of $m_{W}$ and the scale factor $a$ for the moment that the
Hubble parameter $H=m_{W}$ (if $m_{W}>H_{\mathrm{inf}}$, the dominant mode is
about $k\sim a_{\mathrm{end}}H_{\mathrm{inf}}$).

Figure 3: Evolution of $|f_{k}|$ and $\omega^{2}$ for the transverse mode (red
line, solved by the adiabatic approximation) and longitudinal mode (blue line,
solved directly in Eq. (27)) with $m_{W}=10^{10}~{}\mathrm{GeV}$,
$\Gamma=10^{-20}M_{P}$, and $k=0.5,1,2~{}k_{s}$ ($k_{s}$ denotes the dominate
mode of production). It depicts an increasing abundance of dark matters in the
co-moving volume until the Hubble parameter becomes lower than the mass of
WGB. Moreover, the existence of tachyonic enhancement makes the creation of
the longitudinal mode far more productive than the transverse mode.

There are two interesting aspects in Fig. 3. Firstly, the production of the
longitudinal mode is more abundantly than the transverse one, which has been
shown by Graham  _et al._ (2016); Ema  _et al._ (2019). This behavior can be
understood as follows. When $a\ll a_{\mathrm{end}}\lesssim k/m_{W}$, the
frequency $\omega_{\pm}$ of the transverse mode is almost constant, which
leads to a non-increasing mode function as plane wave. However, for the
longitudinal mode, we can show that $\omega_{L}^{2}\sim
k^{2}-a^{\prime\prime}/a$ with $a\ll a_{\mathrm{end}}\lesssim k/m_{W}$, here
we ignore the impact of changing mass for convenience. Since the universe is
approximately regarded as the de Sitter space during the inflationary epoch,
that is $a\sim-1/H_{\mathrm{inf}}\tau$, we have $\omega_{L}^{2}\sim
k^{2}-2a^{2}H_{\mathrm{inf}}^{2}$, which decreases as the universe expands.
Substituting it into Eq. (27), then we obtain an increasing mode function

$\mathcal{W}_{L}|_{a\ll
k/m_{W}}\sim\frac{e^{-ik\tau}}{\sqrt{2k}}\left(1+\frac{iaH_{\mathrm{inf}}}{k}\right).$
(53)

It indicates that when $a\ll k/H_{\mathrm{inf}}$, the mode function of the
longitudinal mode is also a plane wave like the transverse one. But once the
scale factor becomes large enough to make the negative $\omega_{L}^{2}$
emerge, the imaginary part of Eq. (53) will surge exponentially. This
situation is usually called the tachyonic enhancement, which leads to the high
yield of the longitudinal mode compared to the transverse one.

The second noticeable aspect is that the function $f_{k}$ for the dominant
mode $k_{s}$ gradually approaches a stationary value after the Hubble
parameter $H$ becoming smaller than $m_{W}$. This is because the term
$a^{2}m_{W}^{2}$ in Eq. (28) and (29) turns into a leading term after that
moment, which leads to a steady solution with $f_{k}\propto\rho
a^{3}\sim\mathrm{constant}$ (see Ahmed  _et al._ (2020)). And the particles
become non-relativistic due to the expansion of the universe. We can choose an
arbitrary moment that satisfies $H\ll m_{W}$ to evaluate the spectrum of
$f_{k}(k)$ around the dominant mode $k_{s}$, and calculate the integral (36),
(37) or (42) to obtain the energy density $\rho$ of WGB at that moment. In
practice, we consider a moment that is later than the end of reheating as the
conserved quantity $\rho/s$ can be defined for the WGB, where
$s=\frac{2\pi^{2}}{45}g_{*S}(T)T^{3}$ is the total entropy density for the
reheated particles, and $g_{*S}(T)$ denotes the effective number of degrees of
freedom in entropy with temperature $T$. This quantity measures the relic
abundance of non-perturbative produced particles.

Figure 4: (Left) The relation between the relic abundance and the mass of non-
perturbative produced WGB. The feasible mass of WGB should be restricted to
the range that $\rho/s\leq\rho_{0}/s_{0}$. (Right) The maximal abundance
decreases proportionately with the total decay rate of inflaton. It shows that
if $\Gamma<10^{-30}M_{P}$ or $T_{r}<10^{3}~{}\mathrm{GeV}$, the non-
perturbative produced WGB would only be part of dark matter.

Fig. 4 shows how relic abundance varies with $m_{W}$ and decay constant
$\Gamma$, where $m_{W}$ is determined by the parameters $\lambda$ and $g_{W}$,
and the parameter $\beta$ has been fixed according to the constraint of
observation on $H_{\mathrm{inf}}$ discussed in Sec. II. We can see that there
are three cases in the relationship between $\rho/s$ and $m_{W}$, which are
summarized as

$\displaystyle\rho/s\propto\begin{cases}\sqrt{m_{W}},~{}&\mathrm{for}~{}m_{W}<H_{\mathrm{reh}},\\\
\mathrm{almost~{}constant},~{}&\mathrm{for}~{}H_{\mathrm{reh}}<m_{W}<H_{\mathrm{inf}},\\\
\mathrm{exponential~{}decrease},~{}&\mathrm{for}~{}m_{W}>H_{\mathrm{inf}},\end{cases}$
(54)

where $H_{\mathrm{reh}}$ is the Hubble parameter at the end of reheating. For
the current universe, $T_{0}=2.73~{}\mathrm{K}=9.65\times 10^{-32}M_{P}$,
$g_{*S0}=3.91$, and the present energy density for the dark matters is
$\rho_{0}=2.73\times 10^{-121}M_{P}^{4}$, then we can derive
$\rho_{0}/s_{0}=1.8\times 10^{-28}M_{P}$. This observational result gives a
constraint to the mass of WGB. It can be clearly identified in Fig. 4 that the
relic abundance for the WGB produced in a non-perturbative pattern will be
larger than that of dark matters in the current universe if the mass of WGB is
in the interval of $3\times 10^{-13}~{}\mathrm{GeV}$ to $\sim 5\times
10^{13}~{}\mathrm{GeV}$ (The upper limit is in the range of
$10^{13}\mathrm{-}10^{14}~{}\mathrm{GeV}$, which is slightly dependent on the
decay rate $\Gamma$). Therefore, the observation constrains two possible
ranges of mass, $m_{W}\lesssim 3\times 10^{-13}~{}\mathrm{GeV}$ or
$m_{W}\gtrsim 5\times 10^{13}~{}\mathrm{GeV}$, such that the relic density is
not larger than the observation value.

When $m_{W}<3\times 10^{-13}~{}\mathrm{GeV}$, according to the
$\varphi-\tilde{W}$ interaction discussed in Sec. II, there will be other
production channels for WGB through the perturbative decay
$\varphi\rightarrow\tilde{W}\tilde{W}$ and the annihilation
$\varphi\varphi\rightarrow\tilde{W}\tilde{W}$. Since the contribution of decay
is much larger than that of annihilation in our case, we only consider the
decay. Based on Eq. (17) and (18), the decay rate of
$\varphi\rightarrow\tilde{W}\tilde{W}$ can be easily calculated as

$\displaystyle\Gamma(\varphi\rightarrow\tilde{W}\tilde{W})$
$\displaystyle=\frac{3\sinh^{2}(2\Phi_{0}/\sqrt{6})g_{W}^{4}}{\pi}\frac{\sqrt{(m_{\Phi}/2)^{2}-m_{W}^{2}}}{m_{\Phi}^{2}}\left(3+\frac{m_{\Phi}^{4}}{4m_{W}^{4}}-\frac{m_{\Phi}^{2}}{m_{W}^{2}}\right)$
(55)
$\displaystyle\sim\frac{3m_{\Phi}^{3}}{288\pi}\frac{\sinh^{2}(2\Phi_{0}/\sqrt{6})}{\cosh^{4}(\Phi_{0}/\sqrt{6})},~{}\mathrm{when}~{}m_{W}\ll
m_{\Phi}.$

Note that the mass of inflaton is around $m_{\Phi}\sim 2\times
10^{13}~{}\mathrm{GeV}\sim 10^{-5}M_{P}$ in this model, and the typical field
value $\Phi_{0}\sim 10M_{P}$, thus the decay rate can be estimated as
$\Gamma(\varphi\rightarrow\tilde{W}\tilde{W})\sim 10^{-17}M_{P}$. Comparing to
Eq. (47), the decay rate of $\varphi\rightarrow\tilde{W}\tilde{W}$ is far
larger than that of $\varphi\rightarrow\psi\psi$. This situation will lead to
an inadequate reheating to SM, and result in additional dark radiation in our
universe. In such a case, we would need other reheating mechanism that can
transfer most of energy in inflaton field into other particles, and
subsequently SM particles.

For the large mass range $m_{W}\gtrsim 5\times 10^{13}~{}\mathrm{GeV}$, there
is a $-\psi^{2}W^{\mu}W_{\mu}$ term due to the covariant derivative
$D_{\mu}=\partial_{\mu}-W_{\mu}$ in Eq. (45). This implies that the
annihilation $\psi\psi\rightarrow\tilde{W}\tilde{W}$ is allowed in the large
mass case. However, if the reheating temperature is too low, then the
annihilation rate will be negligible, because the number density $n$ of $\psi$
with enough energy to participate in the production of WGB is exponentially
reduced with decreasing temperature $T$, $n\propto\exp(-E/T)$. Therefore, with
a low-temperature reheating mechanism, the non-perturbative production can
reasonably provide enough abundance for WGB as dark matter with the mass
$m_{W}\sim 5\times 10^{13}~{}\mathrm{GeV}$.

Another visible aspect in Fig. 4 is that the decay rate $\Gamma$ strongly
affects the non-perturbative production. Note that $\Gamma$ denotes the total
rate of all decay channels of inflaton, including decay to SM particles and to
dark matters. Hence if $m_{W}<m_{\Phi}/2$, the
$\varphi\rightarrow\tilde{W}\tilde{W}$ channel is open, and the total decay
rate is never lower than $\Gamma(\varphi\rightarrow\tilde{W}\tilde{W})$. Then
if $m_{W}>m_{\Phi}/2$, the total decay rate is equivalent to
$\Gamma(\varphi\rightarrow\psi\psi)$ as the reheating rate. It is apparent on
the right of Fig. 4 that the maximal abundance of non-perturbative produced
WGB decreases proportionately with $\Gamma$, and if $\Gamma<10^{-30}M_{P}$,
the non-perturbative pattern is unable to produce sufficient WGB as dark
matter. Since we have shown $\Gamma(\varphi\rightarrow\tilde{W}\tilde{W})\sim
10^{-17}M_{P}$ in this model (for $m_{W}<m_{\Phi}/2$), the
$\Gamma<10^{-30}M_{P}$ case only appears when $m_{W}>m_{\Phi}/2$, which
corresponds to the reheating temperature $T_{r}\sim 0.4\sqrt{\Gamma
M_{P}}<10^{3}~{}\mathrm{GeV}$. As a result, for this range of $T_{r}$, only
heavy WGB can be a dark matter candidate, but with insufficient abundance.

####  IV.2.3 High-temperature reheating

The preceding discussions consider a low-temperature reheating mechanism that
is constructed for simplicity. In this section, we investigate the production
of WGB with a high-temperature of reheating.

Let us focus on the small mass range first. The perturbative pattern of
production is primarily the decay of inflaton
$\varphi\rightarrow\tilde{W}\tilde{W}$. The annihilation
$\psi\psi\rightarrow\tilde{W}\tilde{W}$ is also present ($\psi$ denotes an
arbitrary scalar in the thermal equilibrium after reheating), but it is
negligible, because the coupling of this process is too small according to the
relation $g_{W}\propto m_{W}$ under our definition. Hence we only need to
consider the decay. In previous sections, we have derived that the decay rate
of $\varphi\rightarrow\tilde{W}\tilde{W}$ is $\Gamma_{W}\sim 10^{-17}M_{P}$.
If this decay rate is far larger than that to the SM particles,
$\Gamma_{W}\gg\Gamma_{SM}$, the WGB can contribute as dark radiation.
Therefore, the upper limit of abundance of dark radiation constrains the lower
limit of $\Gamma_{SM}$ or $T_{r}$. Substituting
$\Gamma=\Gamma_{W}+\Gamma_{SM}$ into Eq. (48), and combining it with the
equation

$\rho_{SM/W}^{\prime}+4\mathcal{H}\rho_{SM/W}-\frac{\Gamma_{SM/W}\Phi^{\prime
2}}{a}=0,$ (56)

we can derive the energy density of dark radiation $\rho_{W}$ after reheating.
Then the amount of dark radiation at present (denoted as
$\Omega_{\mathrm{DR}}$) can be estimated. The result is, if
$\Omega_{\mathrm{DR}}<\Omega_{\mathrm{SM}}$ is expected at matter-radiation
equality, the magnitude of $\Gamma_{SM}$ cannot be lower than the level of
$\Gamma_{W}$, which restricts the reheating temperature to be $T_{r}>4\times
10^{9}~{}\mathrm{GeV}$. In conclusion, under this constraint, the small mass
WGB can be the dark matter candidate, which embodies a tiny amount of dark
radiation produced in the inflaton decay and the majority of cold component
originates from the vacuum fluctuation during inflation. Meanwhile, if it
contributes to the whole abundance of dark matters, the mass should be around
$m_{W}\sim 3\times 10^{-13}~{}\mathrm{GeV}$.

For the heavy WGB case, the decay process is forbidden kinematically. However,
the annihilation channel becomes significant, which is a kind of “freeze-in”
mechanism of production Hall  _et al._ (2010); Blennow  _et al._ (2014);
Bernal  _et al._ (2017). There are two main channels for annihilation, direct
annihilation owing to the $\psi^{2}W^{\mu}W_{\mu}$ coupling, and the other by
gravitational interaction Garny  _et al._ (2016); Tang and Wu (2016, 2017);
Garny  _et al._ (2018). The amplitude $\mathcal{A}\equiv\sum|\mathcal{M}|^{2}$
of them can be calculated as follows

$\mathcal{A}(\psi\psi\rightarrow\tilde{W}\tilde{W})=16g_{W}^{4}\left(3+\frac{\bar{\bm{s}}(\bar{\bm{s}}/4-m_{W}^{2})}{m_{W}^{4}}\right),$
(57)
$\displaystyle\mathcal{A}(\psi\psi\overset{\mathrm{grav}}{\rightarrow}\tilde{W}\tilde{W})=\kappa^{4}\left(\frac{101m_{W}^{4}m_{\psi}^{4}}{60\bar{\bm{s}}^{2}}-\frac{m_{W}^{2}m_{\psi}^{2}}{20\bar{\bm{s}}}(11m_{\psi}^{2}+m_{W}^{2})\right.$
(58)
$\displaystyle\left.+\frac{1}{240}(19m_{\psi}^{4}+76m_{W}^{2}m_{\psi}^{2})-\frac{7\bar{\bm{s}}}{240}(m_{W}^{2}+m_{\psi}^{2})+\frac{\bar{\bm{s}}^{2}}{160}\right),$

where $\bar{\bm{s}}\equiv(p_{\psi 1}+p_{\psi 2})^{2}$,
$\kappa\equiv\sqrt{32\pi G}=2$. Then the total cross section can be easily
derived as

$\sigma_{\mathrm{tot}}=\frac{1}{16\pi\bar{\bm{s}}tg_{i}^{2}}\frac{|\vec{p}_{f}|}{|\vec{p}_{i}|}\mathcal{A}_{\mathrm{tot}},$
(59)

where $t=2$ is the symmetric factor for the identical final states, $g_{i}$ is
the degree of freedom for initial state. Finally, the relic abundance of dark
matters produced in annihilation process can be calculated with the following
formula (see also in Tang and Wu (2017))

$Y=\int_{\sim 0}^{\sim
T_{r}}\frac{dT}{HTs}\left[\frac{Tg_{i}^{2}}{32\pi^{4}}\int
d\bar{\bm{s}}\sigma_{\mathrm{tot}}\sqrt{\bar{\bm{s}}}(\bar{\bm{s}}-4m_{\psi}^{2})K_{1}\left(\frac{\sqrt{\bar{\bm{s}}}}{T}\right)\right],$
(60)

where $H$ is the Hubble parameter, $s$ is the total entropy density mentioned
earlier, and $K_{1}$ is the modified Bessel function of the second kind with
order one. It is seen that for heavy particles, larger relic abundance of dark
matter would result from higher reheating temperature.

With supposing $m_{\psi}\ll m_{W}$, we can calculate the relic abundance $Y$
for each set of $T_{r}$ and $m_{W}$. Then comparing it with present abundance
of dark matter $Y_{0}\sim 5\times 10^{-9}/m_{W}$, we can get the relation
between $m_{W}$ and $T_{r}$ for $Y=Y_{0}$, as shown in Fig. 5. It shows that
the WGB in the large mass range is feasible to be the dark matter candidate
with abundance that is consistent with observation as long as the reheating
temperature $T_{r}>10^{13}~{}\mathrm{GeV}$, and the $m_{W}$ corresponding to
the proper abundance is approximately linear with $T_{r}$. Moreover, since the
instantaneous reheating limit is estimated as
$T_{\mathrm{max}}\sim(3H_{\mathrm{inf}}^{2}/g_{*S}(T_{r}))^{1/4}\sim 2\times
10^{15}~{}\mathrm{GeV}$, the WGB can reach enormous mass of up to $4\times
10^{16}~{}\mathrm{GeV}$.

Figure 5: The proper mass $m_{W}$ for the “freeze in” production (annihilation
channel) that derives an observation-consistent relic abundance is
approximately linear with the reheating temperature $T_{r}$, where the
$m_{W}<10^{14}~{}\mathrm{GeV}$ range is over-produced due to non-perturbative
production. Under the instantaneous reheating limit, the WGB can possess an
enormous mass of up to $4\times 10^{16}~{}\mathrm{GeV}$.

Now we can summarize the main results we have demonstrated in the preceding
discussions. First, if the reheating temperature is too small,
$T_{r}<10^{3}~{}\mathrm{GeV}$, the WGB can only serve as part of dark matter
with insufficient abundance. Because the non-perturbative production is
inefficient for all mass range, and the perturbative contribution is also not
productive or forbidden by the observations. Then, for
$10^{3}~{}\mathrm{GeV}<T_{r}\lesssim 4\times 10^{9}~{}\mathrm{GeV}$, only the
non-perturbative production can reasonably provide enough abundance for WGB as
dark matter, and only the large mass is allowed, $m_{W}\sim
10^{13}\mathrm{-}10^{14}~{}\mathrm{GeV}$. Moreover, for $4\times
10^{9}~{}\mathrm{GeV}\lesssim T_{r}<10^{13}~{}\mathrm{GeV}$, there are two
feasible values of mass for WGB, the small mass $m_{W}\sim 3\times
10^{-13}~{}\mathrm{GeV}$ and the large mass $m_{W}\sim
10^{14}~{}\mathrm{GeV}$. The production channel is also mainly the non-
perturbative pattern. Finally, for
$10^{13}~{}\mathrm{GeV}<T_{r}<T_{\mathrm{max}}$, the small mass case is
similar to the foregoing one, which is restricted to $m_{W}\sim 3\times
10^{-13}~{}\mathrm{GeV}$, and produced mainly in the non-perturbative channel.
But the permissible large mass case is different from the above, because the
perturbative annihilation becomes dominant under extremely high $T_{r}$, which
enables WGB with mass from $10^{14}~{}\mathrm{GeV}$ to $4\times
10^{16}~{}\mathrm{GeV}$ to be a dark matter candidate.

##  V Conclusions

We work on a scaling invariant theory of gravity with a quadratic scalar
curvature, namely the Weyl $R^{2}$ gravity. The model contains a viable
inflationary scenario that is different with the conventional Starobinsky
$R^{2}$ model. We calculate the spectral index and tensor-to-scalar ratio, and
confront them with the latest cosmological observations. Our results indicate
that this model’s predictions agree with the observations and its differences
from Starobinsky model can be tested with future CMB experiment.

The model also contains a stable gauge boson, Weyl gauge boson (WGB), which
can serve as a dark matter candidate. The breaking of local scaling symmetry
makes WGB massive in the inflationary period, which causes an inevitable non-
perturbative production for WGB, namely emerging from the inflationary quantum
fluctuations. Our investigations demonstrate that the relic abundance of non-
perturbative produced WGB is related to the mass of WGB $m_{W}$ and the decay
rate $\Gamma$ of inflaton. By comparing with the observation, we give the
possible viable mass range of WGB, $m_{W}\sim 3\times 10^{-13}~{}\mathrm{GeV}$
or $m_{W}>H_{\mathrm{inf}}$. Moreover, we also give the lower limit of decay
width of the inflaton $\Gamma$ in this model, which corresponds to the
reheating temperature $T_{r}>10^{3}~{}\mathrm{GeV}$.

We also demonstrate that there is a specific coupling between WGB and inflaton
$\Phi$, which can lead to inflaton’s decay into WGB if $m_{W}<m_{\Phi}/2$.
This provides a stringent constraint on the allowable reheating temperature in
the case of small $m_{W}$. Finally, for the case of extremely high reheating
temperature $T_{r}$, the perturbative freeze-in channel becomes important and
the dominant contribution, namely the annihilation of particles in the thermal
bath into WGB through direct coupling or gravitational interaction. It opens
new mass range of WGB up to $4\times 10^{16}~{}\mathrm{GeV}$ as dark matter
candidate.

###### Acknowledgements.

YT is supported by National Key Research and Development Program of China
(Grant No.2021YFC2201901), Natural Science Foundation of China (NSFC) under
Grants No. 11851302, the Fundamental Research Funds for the Central
Universities and Key Research Program of the Chinese Academy of Sciences,
Grant No. XDPB15. YLW is supported in part by the National Key Research and
Development Program of China under Grant No.2020YFC2201501, and NSFC under
Grants No. 11851302, No. 11851303, No. 11690022, No. 11747601, and the
Strategic Priority Research Program of the Chinese Academy of Sciences under
Grant No. XDB23030100.

## References

  * Rubin and Ford (1970) V. C. Rubin and W. K. Ford, Jr., Astrophys. J. 159, 379 (1970).
  * Allen  _et al._ (2003) S. W. Allen, A. C. Fabian, R. W. Schmidt, and H. Ebeling, Mon. Not. Roy. Astron. Soc. 342, 287 (2003), arXiv:astro-ph/0208394 .
  * Refregier (2003) A. Refregier, Ann. Rev. Astron. Astrophys. 41, 645 (2003), arXiv:astro-ph/0307212 .
  * Tyson  _et al._ (1998) J. A. Tyson, G. P. Kochanski, and I. P. Dell’Antonio, Astrophys. J. Lett. 498, L107 (1998), arXiv:astro-ph/9801193 .
  * Fields  _et al._ (2020) B. Fields, P. Molaro, and S. Sarkar (Particle Data Group Collaboration), Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Lin (2019) T. Lin, PoS 333, 009 (2019), arXiv:1904.07915 [hep-ph] .
  * Schumann (2019) M. Schumann, J. Phys. G 46, 103003 (2019), arXiv:1903.03026 [astro-ph.CO] .
  * Billard  _et al._ (2021) J. Billard _et al._ , (2021), arXiv:2104.07634 [hep-ex] .
  * Tang and Wu (2020a) Y. Tang and Y.-L. Wu, Phys. Lett. B 803, 135320 (2020a), arXiv:1904.04493 [hep-ph] .
  * Tang and Wu (2020b) Y. Tang and Y.-L. Wu, JCAP 03, 067 (2020b), arXiv:1912.07610 [hep-ph] .
  * Wu (2016) Y.-L. Wu, Phys. Rev. D93, 024012 (2016), arXiv:1506.01807 [hep-th] .
  * Wu (2018) Y.-L. Wu, Eur. Phys. J. C78, 28 (2018), arXiv:1712.04537 [hep-th] .
  * Weyl (1929) H. Weyl, Z. Phys. 56, 330 (1929).
  * Mannheim and Kazanas (1989) P. D. Mannheim and D. Kazanas, Astrophys. J. 342, 635 (1989).
  * Romero  _et al._ (2012) C. Romero, J. B. Fonseca-Neto, and M. L. Pucheu, Class. Quant. Grav. 29, 155015 (2012), arXiv:1201.1469 [gr-qc] .
  * Bars  _et al._ (2014) I. Bars, P. Steinhardt, and N. Turok, Phys. Rev. D 89, 043515 (2014), arXiv:1307.1848 [hep-th] .
  * Scholz (2018) E. Scholz, Einstein Stud. 14, 261 (2018), arXiv:1703.03187 [math.HO] .
  * Ghilencea (2019) D. M. Ghilencea, JHEP 10, 209 (2019), arXiv:1906.11572 [gr-qc] .
  * Ferreira  _et al._ (2019) P. G. Ferreira, C. T. Hill, J. Noller, and G. G. Ross, Phys. Rev. D 100, 123516 (2019), arXiv:1906.03415 [gr-qc] .
  * Tang and Wu (2020c) Y. Tang and Y.-L. Wu, Phys. Lett. B 809, 135716 (2020c), arXiv:2006.02811 [hep-ph] .
  * Beltran Jimenez  _et al._ (2016) J. Beltran Jimenez, L. Heisenberg, and T. S. Koivisto, JCAP 04, 046 (2016), arXiv:1602.07287 [hep-th] .
  * Oda (2020) I. Oda, Mod. Phys. Lett. A 35, 2050304 (2020), arXiv:2006.10867 [hep-th] .
  * Ghilencea (2021) D. M. Ghilencea, Eur. Phys. J. C 81, 510 (2021), arXiv:2007.14733 [hep-th] .
  * Cai  _et al._ (2021) R.-G. Cai, Y.-S. Hao, and S.-J. Wang, (2021), arXiv:2110.14718 [gr-qc] .
  * Parker (1969) L. Parker, Phys. Rev. 183, 1057 (1969).
  * Parker (1971) L. Parker, Phys. Rev. D 3, 346 (1971), [Erratum: Phys.Rev.D 3, 2546–2546 (1971)].
  * Ford (1987) L. H. Ford, Phys. Rev. D 35, 2955 (1987).
  * Lyth and Roberts (1998) D. H. Lyth and D. Roberts, Phys. Rev. D 57, 7120 (1998), arXiv:hep-ph/9609441 .
  * Chung  _et al._ (2001) D. J. H. Chung, P. Crotty, E. W. Kolb, and A. Riotto, Phys. Rev. D 64, 043503 (2001), arXiv:hep-ph/0104100 .
  * Baer  _et al._ (2015) H. Baer, K.-Y. Choi, J. E. Kim, and L. Roszkowski, Phys. Rept. 555, 1 (2015), arXiv:1407.0017 [hep-ph] .
  * Graham  _et al._ (2016) P. W. Graham, J. Mardon, and S. Rajendran, Phys. Rev. D 93, 103520 (2016), arXiv:1504.02102 [hep-ph] .
  * Ema  _et al._ (2018) Y. Ema, K. Nakayama, and Y. Tang, JHEP 09, 135 (2018), arXiv:1804.07471 [hep-ph] .
  * Ema  _et al._ (2019) Y. Ema, K. Nakayama, and Y. Tang, JHEP 07, 060 (2019), arXiv:1903.10973 [hep-ph] .
  * Ahmed  _et al._ (2020) A. Ahmed, B. Grzadkowski, and A. Socha, JHEP 08, 059 (2020), arXiv:2005.01766 [hep-ph] .
  * Chung  _et al._ (2019) D. J. H. Chung, E. W. Kolb, and A. J. Long, JHEP 01, 189 (2019), arXiv:1812.00211 [hep-ph] .
  * Kolb and Long (2021) E. W. Kolb and A. J. Long, JHEP 03, 283 (2021), arXiv:2009.03828 [astro-ph.CO] .
  * Ema  _et al._ (2017) Y. Ema, R. Jinno, K. Mukaida, and K. Nakayama, JCAP 02, 045 (2017), arXiv:1609.05209 [hep-ph] .
  * Ema  _et al._ (2016) Y. Ema, R. Jinno, K. Mukaida, and K. Nakayama, Phys. Rev. D 94, 063517 (2016), arXiv:1604.08898 [hep-ph] .
  * Li  _et al._ (2020) L. Li, S. Lu, Y. Wang, and S. Zhou, JHEP 07, 231 (2020), arXiv:2002.01131 [hep-ph] .
  * Li  _et al._ (2019) L. Li, T. Nakama, C. M. Sou, Y. Wang, and S. Zhou, JHEP 07, 067 (2019), arXiv:1903.08842 [astro-ph.CO] .
  * Herring and Boyanovsky (2020) N. Herring and D. Boyanovsky, Phys. Rev. D 101, 123522 (2020), arXiv:2005.00391 [astro-ph.CO] .
  * Babichev  _et al._ (2019) E. Babichev, D. Gorbunov, and S. Ramazanov, Phys. Lett. B 794, 69 (2019), arXiv:1812.03516 [hep-ph] .
  * Hashiba and Yokoyama (2019) S. Hashiba and J. Yokoyama, Phys. Rev. D 99, 043008 (2019), arXiv:1812.10032 [hep-ph] .
  * Ling and Long (2021) S. Ling and A. J. Long, Phys. Rev. D 103, 103532 (2021), arXiv:2101.11621 [astro-ph.CO] .
  * Garny  _et al._ (2016) M. Garny, M. Sandora, and M. S. Sloth, Phys. Rev. Lett. 116, 101302 (2016), arXiv:1511.03278 [hep-ph] .
  * Tang and Wu (2016) Y. Tang and Y.-L. Wu, Phys. Lett. B 758, 402 (2016), arXiv:1604.04701 [hep-ph] .
  * Tang and Wu (2017) Y. Tang and Y.-L. Wu, Phys. Lett. B 774, 676 (2017), arXiv:1708.05138 [hep-ph] .
  * Garny  _et al._ (2018) M. Garny, A. Palessandro, M. Sandora, and M. S. Sloth, JCAP 02, 027 (2018), arXiv:1709.09688 [hep-ph] .
  * Chen and Kang (2018) S.-L. Chen and Z. Kang, JCAP 05, 036 (2018), arXiv:1711.02556 [hep-ph] .
  * Bernal  _et al._ (2018) N. Bernal, M. Dutra, Y. Mambrini, K. Olive, M. Peloso, and M. Pierre, Phys. Rev. D 97, 115020 (2018), arXiv:1803.01866 [hep-ph] .
  * Aoki  _et al._ (2022) S. Aoki, H. M. Lee, A. G. Menkara, and K. Yamashita, (2022), arXiv:2202.13063 [hep-ph] .
  * Clery  _et al._ (2021) S. Clery, Y. Mambrini, K. A. Olive, and S. Verner, (2021), arXiv:2112.15214 [hep-ph] .
  * Frangipane  _et al._ (2021) E. Frangipane, S. Gori, and B. Shakya, (2021), arXiv:2110.10711 [hep-ph] .
  * Redi and Tesi (2021) M. Redi and A. Tesi, JHEP 12, 060 (2021), arXiv:2107.14801 [hep-ph] .
  * Barman  _et al._ (2022) B. Barman, N. Bernal, A. Das, and R. Roshan, JCAP 01, 047 (2022), arXiv:2108.13447 [hep-ph] .
  * Mambrini and Olive (2021) Y. Mambrini and K. A. Olive, Phys. Rev. D 103, 115009 (2021), arXiv:2102.06214 [hep-ph] .
  * Haque and Maity (2021) M. R. Haque and D. Maity, (2021), arXiv:2112.14668 [hep-ph] .
  * Clery  _et al._ (2022) S. Clery, Y. Mambrini, K. A. Olive, A. Shkerin, and S. Verner, (2022), arXiv:2203.02004 [hep-ph] .
  * Wu (2021a) Y.-L. Wu, Int. J. Mod. Phys. A 36, 2143001 (2021a), arXiv:2104.05404 [physics.gen-ph] .
  * Wu (2021b) Y.-L. Wu, Int. J. Mod. Phys. A 36, 2143002 (2021b), arXiv:2104.11078 [physics.gen-ph] .
  * Barrow and Cotsakis (1988) J. D. Barrow and S. Cotsakis, Phys. Lett. B 214, 515 (1988).
  * De Felice and Tsujikawa (2010) A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Starobinsky (1980) A. A. Starobinsky, Phys. Lett. B 91, 99 (1980).
  * Maeda (1988) K.-i. Maeda, Phys. Rev. D 37, 858 (1988).
  * Cembranos (2009) J. A. R. Cembranos, Phys. Rev. Lett. 102, 141301 (2009), arXiv:0809.1653 [hep-ph] .
  * Bezrukov and Gorbunov (2012) F. L. Bezrukov and D. S. Gorbunov, Phys. Lett. B 713, 365 (2012), arXiv:1111.4397 [hep-ph] .
  * He  _et al._ (2018) M. He, A. A. Starobinsky, and J. Yokoyama, JCAP 05, 064 (2018), arXiv:1804.00409 [astro-ph.CO] .
  * Enckell  _et al._ (2019) V.-M. Enckell, K. Enqvist, S. Rasanen, and L.-P. Wahlman, JCAP 02, 022 (2019), arXiv:1810.05536 [gr-qc] .
  * Pi  _et al._ (2018) S. Pi, Y.-l. Zhang, Q.-G. Huang, and M. Sasaki, JCAP 05, 042 (2018), arXiv:1712.09896 [astro-ph.CO] .
  * Antoniadis  _et al._ (2018) I. Antoniadis, A. Karam, A. Lykkas, and K. Tamvakis, JCAP 11, 028 (2018), arXiv:1810.10418 [gr-qc] .
  * Ade  _et al._ (2021) P. A. R. Ade _et al._ (BICEP, Keck), Phys. Rev. Lett. 127, 151301 (2021), arXiv:2110.00483 [astro-ph.CO] .
  * Liddle and Leach (2003) A. R. Liddle and S. M. Leach, Phys. Rev. D 68, 103503 (2003), arXiv:astro-ph/0305263 .
  * Abazajian  _et al._ (2022) K. Abazajian _et al._ (CMB-S4), Astrophys. J. 926, 54 (2022), arXiv:2008.12619 [astro-ph.CO] .
  * Kofman  _et al._ (1997) L. Kofman, A. D. Linde, and A. A. Starobinsky, Phys. Rev. D 56, 3258 (1997), arXiv:hep-ph/9704452 .
  * Albrecht  _et al._ (1982) A. Albrecht, P. J. Steinhardt, M. S. Turner, and F. Wilczek, Phys. Rev. Lett. 48, 1437 (1982).
  * Kofman  _et al._ (1994) L. Kofman, A. D. Linde, and A. A. Starobinsky, Phys. Rev. Lett. 73, 3195 (1994), arXiv:hep-th/9405187 .
  * Bassett  _et al._ (2006) B. A. Bassett, S. Tsujikawa, and D. Wands, Rev. Mod. Phys. 78, 537 (2006), arXiv:astro-ph/0507632 .
  * Lozanov (2019) K. D. Lozanov, (2019), arXiv:1907.04402 [astro-ph.CO] .
  * Bunch and Davies (1978) T. S. Bunch and P. C. W. Davies, Proc. Roy. Soc. Lond. A 360, 117 (1978).
  * Hall  _et al._ (2010) L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, JHEP 03, 080 (2010), arXiv:0911.1120 [hep-ph] .
  * Blennow  _et al._ (2014) M. Blennow, E. Fernandez-Martinez, and B. Zaldivar, JCAP 01, 003 (2014), arXiv:1309.7348 [hep-ph] .
  * Bernal  _et al._ (2017) N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, Int. J. Mod. Phys. A 32, 1730023 (2017), arXiv:1706.07442 [hep-ph] .