# Resonant Reheating

Basabendu Barman  Nicolás Bernal  Yong Xu

###### Abstract

We investigate a novel reheating scenario proceeding through $s$-channel
inflaton annihilation, mediated by a massive scalar. If the inflaton $\phi$
oscillates around the minimum of a monomial potential $\propto\phi^{n}$, we
reveal the emergence of resonance phenomena originating from the dynamic
evolution of the inflaton mass for $n>2$. Consequently, a resonance appears in
both the radiation and the temperature evolution during the reheating process.
By solving the coupled Boltzmann equations, we present solutions for radiation
and temperature. We find non-trivial temperature characteristics during
reheating, depending on the value of $n$ and the masses of the inflaton and
mediator. Some phenomenological aspects of the model are explored. As a
concrete example, we show that the same mediator participates in the genesis
of dark matter, modifying the standard freeze-in dynamics. In addition, we
demonstrate that the resonant reheating scenario could be tested by next-
generation low- and high-frequency gravitational wave detectors.

MITP-24-045

##  1 Introduction

Cosmic inflation is an elegant paradigm to solve several problems in cosmology
[1, 2, 3, 4]. In the conventional cosmological narrative, the initially
(quasi) exponential expansion leads the Universe to be cold at the end of
inflation. Understanding the transition from this state to a hot, thermal,
radiation-dominated Universe, laying the foundation for Big Bang
nucleosynthesis (BBN), is vital. This heating process not only elucidates the
cosmic origins of the matter we observe, but also accounts for relics such as
photons, neutrinos, and primordial gravitational waves (GW), along with exotic
particles beyond the Standard Model (SM), for instance dark matter (DM).

The basic idea of reheating after inflation is as follows: it occurs due to
the production of particles by the oscillating scalar field $\phi$ [5, 6]. In
the most basic inflationary models, this field, responsible for driving
inflation, subsequently undergoes oscillations near the minimum of its
effective potential $V(\phi)\propto\phi^{n}$, thus inducing particle
production. Reheating can be performed through interactions between the
inflaton and some light degrees of freedom. Such interactions, which result in
energy transfer from the inflaton to the SM fields, can proceed through
inflaton decays, for example, as discussed in Refs. [7, 8, 9] or inflaton
scattering [10, 11].111The effect of time-dependent inflaton decay on
reheating dynamics has been discussed in, for example, Refs. [12, 13, 14, 15,
16].

In the case of reheating through inflaton annihilations, contact interactions
between a pair of inflatons and a pair of SM particles are typically assumed.
Inflaton annihilations through a mediator have been studied in the framework
of gravity, where the graviton plays the role of mediator. However, it has
been shown that rather steep inflaton potentials with $n>9$ are needed for
successful reheating [17, 18, 19, 20, 21]. Interestingly, this bound can be
relaxed to $n>4$ if one introduces a non-minimal coupling between gravity and
a pair of inflatons [22], or even to $n\geq 2$ if gravity couples non-
minimally to a single inflaton [23].

In this work, however, we investigate a previously unexplored novel
possibility of reheating through scalar-mediated $s$-channel scattering. In
particular, we consider a model in which a pair of inflatons $(\phi)$
annihilate into a pair of daughter particles through the $s$-channel exchange
of a massive spin-0 mediator $(S)$. The daughter particles could be Higgs
bosons in the SM or particles beyond the SM like right-handed neutrinos. If
the inflaton oscillates in a potential steeper than that of the quadratic, the
inflaton mass is not a constant but decreases with time. Consequently, the
existence of higher modes can give enough energy to the inflaton so that
mediators are produced resonantly. Once these mediator fields are produced,
they quickly decay into SM particles that make up the thermal bath. Because of
the resonant production of the mediator, the evolution of the SM temperature
features non-trivial behavior (typically a bump) during reheating, which is
very different from the standard cases characterized by a simple power law. We
find that in the presence of a resonance, the maximum temperature during
reheating is controlled by the location of the resonance relative to the end
of the reheating. In particular, when the resonance occurs close to the end of
reheating, the reheating temperature would be the maximum temperature and can
be larger than that where the resonance is absent.

Interestingly, the mediator $S$ could also play a fundamental role in the
production of particles beyond the SM, such as the DM. In particular, if the
DM has very suppressed couplings with the SM so that it never reaches thermal
equilibrium, it can be produced in the early Universe through the FIMP
mechanism [24, 25, 26, 27, 28, 29]. The dynamics of reheating determines the
abundance of DM if it is produced during reheating, and in that case it is
called an ultraviolet (UV) FIMP [30, 31, 32]. We further explore the impact of
such a resonant reheating on DM production from the visible sector or
inflaton, mediated by $S$ or graviton. We find that new and interesting
freeze-in behaviors show up mainly due to the presence of the resonance
feature during reheating.

Primordial gravitational waves (GW) serve as a robust prediction for
inflationary models. In the course of inflation, quantum fluctuations
inevitably lead to a spectrum of tensor-metric perturbations on scales beyond
the Hubble horizon. In a conventional post-inflationary scenario, these tensor
modes transition to subhorizon scales during the radiation-dominated epoch.
However, it is widely recognized that the presence of nonstandard cosmological
conditions, such as an early stiff era preceding radiation domination,
disrupts this scale invariance [33]. Under such circumstances, the GW spectrum
experiences a significant blue tilt within the frequency range corresponding
to the modes that cross the horizon during the stiff period [34, 35, 36, 37,
38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]. Such a blue-tilted
spectrum turns out to be well within the reach of several GW detector
facilities. Our scenario turns out to be testable on several next-generation
GW detectors.

The paper is organized as follows. In Section 2, we offer the model setup and
explore the features of resonance reheating with both numerical results and
analytical approximations. In Section 3, we investigate the implications for
the production of DM due to the presence of resonance. Section 4 is devoted to
discussions of the testability of the present framework with the primordial GW
spectrum. We sum up our findings in Section 5.

##  2 Post-inflationary reheating dynamics

After the end of cosmic inflation, inflatons oscillate around the minimum of
their potential, transferring energy to SM particles that subsequently
thermalize, forming the SM thermal bath. We consider the post-inflationary
oscillation of the inflaton $\phi$ at the bottom of a monomial potential
$V(\phi)$ of the form

$V(\phi)=\lambda\,\frac{\phi^{n}}{M_{P}^{n-4}}\,,$ (2.1)

where $\lambda$ is a dimensionless coupling and $M_{P}\simeq 2.4\times
10^{18}$ GeV the reduced Planck mass. This type of potential could arise from
$\alpha$-attractor $T$\- or $E$-models [51, 52], or the Starobinsky
inflationary model [1, 53, 54, 55]. Defining the energy density and pressure
of $\phi$ as $\rho_{\phi}\equiv\frac{1}{2}\,\dot{\phi}^{2}+V(\phi)$ and
$p_{\phi}\equiv\frac{1}{2}\,\dot{\phi}^{2}-V(\phi)$, the background equation-
of-state parameter is defined as $w\equiv p_{\phi}/\rho_{\phi}=(n-2)/(n+2)$
[56]. For example, for $n=4$, $w=1/3$ corresponding to a radiation-like
background.

The evolution of the energy densities of the inflaton $\rho_{\phi}$ and the SM
radiation $\rho_{R}$ can be tracked using the Boltzmann equations222We also
assume that the thermalization of the SM plasma occurs instantaneously. We
note, however, that the SM particles do not necessarily thermalize suddenly,
and thus they could be initially distributed with smaller occupation numbers
and harder momenta [57, 58, 59, 60, 61].

$\displaystyle\frac{d\rho_{\phi}}{dt}+\frac{6\,n}{2+n}\,H\,\rho_{\phi}$
$\displaystyle=-\frac{2\,n}{2+n}\,\Gamma^{2\to 2}\,\rho_{\phi}$ (2.2)
$\displaystyle\frac{d\rho_{R}}{dt}+4\,H\,\rho_{R}$
$\displaystyle=+\frac{2\,n}{2+n}\,\Gamma^{2\to 2}\,\rho_{\phi}\,,$ (2.3)

where

$H=\sqrt{\frac{\rho_{R}+\rho_{\phi}}{3\,M_{P}^{2}}}\,,$ (2.4)

is the Hubble expansion rate, and $\Gamma^{2\to 2}$ the 2-to-2 annihilation
rate of a pair of inflaton condensate into a pair of SM states. The detailed
expression of $\Gamma^{2\to 2}$ will be discussed in a moment.

During reheating, that is, when $a_{I}\leq a\leq a_{\text{rh}}$, where $a$ is
the cosmic scale factor, the expansion $H$ dominates over the interaction
$\Gamma^{2\to 2}$, and then

$\rho_{\phi}(a)\simeq\rho_{\phi}(a_{I})\left(\frac{a_{I}}{a}\right)^{\frac{6\,n}{2+n}}.$
(2.5)

Here, $a_{I}$ and $a_{\text{rh}}$ correspond to the scale factor at the end of
inflation (that is, at the beginning of reheating) and at the end of
reheating, respectively. Since the Hubble rate during reheating is dominated
by the inflaton energy density, it follows that

$H(a)\simeq
H(a_{\text{rh}})\times\begin{dcases}\left(\frac{a_{\text{rh}}}{a}\right)^{\frac{3\,n}{n+2}}&\text{
for }a_{I}\leq a\leq a_{\text{rh}}\,,\\\
\left(\frac{a_{\text{rh}}}{a}\right)^{2}&\text{ for }a_{\text{rh}}\leq
a\,,\end{dcases}$ (2.6)

taking into account that after the end of reheating the Hubble expansion rate
is dominated by free SM radiation.

The effective mass $m_{\phi}$ of the inflaton can be obtained from the second
derivative of Eq. (2.1), which reads [62]

$m_{\phi}^{2}(a)\equiv\frac{d^{2}V}{d\phi^{2}}=n\,(n-1)\,\lambda\,M_{P}^{2}\left(\frac{\phi}{M_{P}}\right)^{n-2}\simeq
n\,(n-1)\,\lambda^{\frac{2}{n}}\,M_{P}^{2}\left(\frac{\rho_{\phi}(a)}{M_{P}^{4}}\right)^{\frac{n-2}{n}}.$
(2.7)

In the last step, we have considered $\rho_{\phi}\simeq V(\phi)$ and utilized
the relation between the inflaton field $\phi$ and the potential given in Eq.
(2.1). The inflaton mass at a scale factor $a$ can be related to its value
$m_{I}\equiv m_{\phi}(a_{I})$ at the beginning of the reheating by

$m_{\phi}(a)\simeq m_{I}\left(\frac{a_{I}}{a}\right)^{\frac{3(n-2)}{n+2}}.$
(2.8)

It is interesting to note that for $n\neq 2$, $m_{\phi}$ has a field
dependence that will be inherited by the interaction rate $\Gamma^{2\to 2}$.
Additionally, for $n=2$, radiation cannot exceed the inflaton energy density,
and therefore reheating is not viable if it proceeds through annihilations
with $n<5/2$ [11].333Reheating through annihilations and $n=2$ could become
viable in the case of an inflaton nonminimally coupled to gravity [23]. In the
following, we will therefore focus on the case $n\geq 4$.

###  2.1 Bosonic reheating

As advocated in the beginning, here we are interested in the 2-to-2 scattering
of the inflaton $\phi\phi\to S\to hh$, mediated by a real singlet scalar $S$
in the $s$ channel.444We note that if $\phi\phi\to S\to hh$ is present, one
can not avoid the annihilation $\phi\phi\to SS$ through $t$ and $u$ channels.
Subsequently, the mediators $S$ decay into a pair of SM particles $S\to hh$.
However, in the present setup, this channel is subdominant, as explored in
Appendix A. This can arise from an effective Lagrangian density interaction of
the form

$\mathcal{L}\supset\frac{1}{2}\,\widetilde{\mu}_{\phi}\,\phi^{2}\,S+\frac{1}{2}\,\mu_{h}\,S\,h^{2},$
(2.9)

where $\phi$ is the inflaton field, $S$ is the mediator and $h$ a SM-like
field.555We assume quartic interaction terms, e.g.,
$\phi^{2}\,S^{2},\,S^{2}\,h^{2},\,\phi^{2}\,h^{2}$ or trilinear interaction
between the inflaton and a pair of SM Higgs $\phi\,h^{2}$, are absent. We
remain agnostic about the UV completion of our effective toy model. Note that
both couplings have a mass dimension. Since the inflaton oscillates in a
monomial potential, one must take care of the Fourier modes of the
annihilating inflaton field. In that case, after averaging over several
oscillations, the scattering rate for $\phi\phi\to S\to hh$ turns out to be

$\Gamma^{\phi\phi\to
hh}\simeq\frac{\rho_{\phi}}{m_{\phi}}\,\frac{\mu_{\phi}^{2}\,\mu_{h}^{2}}{32\pi\,m_{\phi}^{2}\left[(4m_{\phi}^{2}-m_{s}^{2})^{2}+\Gamma_{s}^{2}\,m_{s}^{2}\right]}\sqrt{1-\frac{m_{h}^{2}}{m_{\phi}^{2}}}\,\Theta(m_{\phi}-m_{h})\,,$
(2.10)

where $S$ has a mass $m_{s}$ and a total decay width

$\Gamma_{s}=\frac{1}{8\pi}\,\frac{\mu_{h}^{2}}{m_{s}}\sqrt{1-\frac{4m_{h}^{2}}{m_{s}^{2}}}\,\Theta(m_{s}-2m_{h})+\frac{1}{8\pi}\,\frac{\mu_{\phi}^{2}}{m_{s}}\sqrt{1-\frac{4\,m_{\phi}^{2}}{m_{s}^{2}}}\,\Theta(m_{s}-2m_{\phi})\,,$
(2.11)

with $\Theta$ being the Heaviside step function. Here, the effective coupling
$\mu_{\phi}$ is proportional to $\widetilde{\mu}_{\phi}$ and encodes the
summation over several oscillations of Fourier modes [63, 64, 65, 10, 62, 66].
It is interesting to note that the rate in Eq. (2.10) features a pole when
$m_{\phi}(a)=m_{s}/2$. As the inflaton mass varies during reheating, this
resonance could be reached. It occurs at the scale factor $a_{p}$ given by

$a_{p}\equiv a_{I}\left(\frac{2\,m_{I}}{m_{s}}\right)^{\frac{2+n}{3(n-2)}}.$
(2.12)

We emphasize that here we are interested in the case where $n\geq 4$.
Additionally, we will focus on the case where the resonance is reached during
the reheating era.

This system of Boltzmann equations (2.2) and (2.3) can be solved analytically.
Taking the initial conditions at $a=a_{I}$ to be $\rho_{R}(a_{I})=0$ and
$\rho_{\phi}(a_{I})=3\,H_{I}^{2}\,M_{P}^{2}$, with $H_{I}$ being the
inflationary scale,666The recent BICEP/Keck bound on tensor-to-scalar ratio
implies that $H_{I}\leq 2.0\times 10^{-5}~{}M_{P}$ [67]. This limit can also
be used to derive an upper bound on the reheating temperature:
$T_{\text{rh}}\lesssim 2.5\times 10^{15}$ GeV under the assumption that the
reheating is instantaneous [68]. the SM energy density can be approximated by

$\displaystyle\rho_{R}(a)$
$\displaystyle\simeq\frac{H_{I}^{3}\,M_{P}^{4}\,\mu_{h}^{2}\,\mu_{\phi}^{2}}{m_{I}^{7}}$
$\displaystyle\qquad\times\begin{dcases}\frac{9}{512\pi}\frac{n}{8n-17}\left(\frac{a}{a_{I}}\right)^{\frac{6(2n-7)}{2+n}}\left[1-\left(\frac{a_{I}}{a}\right)^{\frac{2(8n-17)}{n+2}}\right]&\text{for
}a_{I}\leq a\ll a_{p}\,,\\\
\frac{3\pi}{2}\frac{n}{n-2}\left(\frac{a_{I}}{a}\right)^{4}\left(\frac{2m_{I}}{m_{s}}\right)^{\frac{n-4}{3(n-2)}}\frac{m_{I}^{5}}{m_{s}^{3}\,\mu_{h}^{2}}&\text{for
}a_{p}\lesssim a\lesssim a_{t}\,,\\\
\frac{9}{32\pi}\frac{n}{2n-5}\left(\frac{a_{I}}{a}\right)^{\frac{18}{n+2}}\left[1-\left(\frac{a_{I}}{a}\right)^{\frac{2(2n-5)}{2+n}}\right]\left(\frac{m_{I}}{m_{s}}\right)^{4}&\text{for
}a_{t}\lesssim a\leq a_{\text{rh}}\,.\end{dcases}$ (2.13)

Finally, after the end of the reheating, when $a>a_{\text{rh}}$, the SM
radiation scales as free radiation and therefore $\rho_{R}(a)\propto a^{-4}$.
Some comments on the analytical solution in Eq. (2.1) are in order: $i)$
before the resonance, that is, for $a_{I}\leq a\ll a_{p}$, the terms
proportional to $\Gamma_{s}$ and $m_{s}$ in the denominator can be ignored,
$ii)$ near the resonance, the so-called narrow width approximation

$\frac{1}{(s-m^{2})^{2}+m^{2}\,\Gamma^{2}}\to\frac{\pi}{m\,\Gamma}\,\delta\left(s-m^{2}\right)$
(2.14)

is valid. It is worth mentioning here that although the resonance occurs at a
particular epoch $a=a_{p}$, it persists over a finite time. And $iii)$ after
the resonance, one can safely take $\Gamma_{s}=m_{I}=0$ in the denominator.

By equating the second and third solutions of Eq. (2.1), one finds that the
transition after the pole occurs at the scale factor $a_{t}$ defined as

$a_{t}\equiv
a_{I}\left[\frac{32\pi^{2}}{3}\,\frac{2n-5}{n-2}\left(\frac{m_{I}}{\mu_{h}}\right)^{2}\left(\frac{m_{s}}{2m_{I}}\right)^{\frac{2}{3}\frac{n-1}{n-2}}\right]^{\frac{2+n}{2(2n-5)}}.$
(2.15)

Also, the end of reheating, defined as the onset of the radiation-dominated
era, can be found by equating the third solution of Eq. (2.1) with Eq. (2.5),
and the corresponding scale factor $a_{\text{rh}}$ is

$a_{\text{rh}}\simeq
a_{I}\left(\frac{32\pi}{3}\frac{2n-5}{n}\frac{m_{I}^{3}\,m_{s}^{4}}{H_{I}\,M_{P}^{2}\,\mu_{h}^{2}\,\mu_{\phi}^{2}}\right)^{\frac{2+n}{6(n-3)}}.$
(2.16)

Finally, the SM temperature $T$ can be extracted from the SM radiation energy
density $\rho_{R}$ by the fact that

$\rho_{R}(T)=\frac{\pi^{2}}{30}\,g_{\star}(T)\,T^{4},$ (2.17)

where $g_{\star}(T)$ corresponds to the number of relativistic degrees of
freedom contributing to $\rho_{R}$. At $a=a_{p}$ the SM bath reaches a maximum
temperature $T_{\text{max}}$ [69] given by

$T_{\text{max}}\simeq\left[\frac{45}{\pi\,g_{\star}}\frac{n}{n-2}\left(\frac{m_{s}}{2m_{I}}\right)^{\frac{n+4}{n-2}}\frac{H_{I}^{3}\,M_{P}^{4}\,\mu_{\phi}^{2}}{m_{I}^{2}\,m_{s}^{3}}\right]^{\frac{1}{4}}\,.$
(2.18)

Furthermore, at the end of reheating, when $a=a_{\text{rh}}$, the SM
temperature $T_{\text{rh}}$ is

$T_{\text{rh}}\simeq\left[\frac{90}{\pi^{2}\,g_{\star}}H_{I}^{2}\,M_{P}^{2}\left(\frac{3}{32\pi}\frac{n}{2n-5}\frac{H_{I}\,M_{P}^{2}\,\mu_{h}^{2}\,\mu_{\phi}^{2}}{m_{I}^{3}\,m_{s}^{4}}\right)^{\frac{n}{n-3}}\right]^{\frac{1}{4}}\,.$
(2.19)

In general, from Eq. (2.1) the scaling of the SM temperature can be extracted,
and it is given by

$T(a)\propto\begin{dcases}a^{+\frac{3(2n-7)}{2(n+2)}}&\text{for }a_{I}\leq
a\ll a_{p}\,,\\\ a^{-1}&\text{for }a_{p}\leq a\ll a_{t}\,,\\\
a^{-\frac{9}{2(n+2)}}&\text{for }a_{t}\leq a\ll a_{\text{rh}}\,.\end{dcases}$
(2.20)

We note that before the pole the SM temperature continuously increases with a
slope $T(a)\propto a^{+\frac{3(2n-7)}{2(n+2)}}$, reaching the maximum
$T=T_{\text{max}}$ at $a=a_{p}$. Interestingly, near the resonance, the
temperature decreases as $T(a)\propto a^{-1}$, independent of the value of
$n$. In this case, the entropy injection and the Hubble expansion are
counterbalanced, so that the SM radiation scales as free radiation. Well after
the pole, the SM temperature continues to decrease, but more slowly, as
$T(a)\propto a^{-\frac{9}{2(n+2)}}$ until the end of reheating at
$a=a_{\text{rh}}$. After reheating, when the Universe becomes dominated by SM
radiation, the SM entropy is conserved and therefore $T(a)\propto a^{-1}$.
Before moving on, we clarify that in the present setup there are six free
parameters that dictate the reheating dynamics, and can be conveniently chosen
to be

$n,\,m_{I},\,H_{I},\,m_{s},\,\mu_{h},\text{ and }\mu_{\phi}\,,$ (2.21)

or

$n,\,m_{I},\,H_{I},\,\frac{a_{p}}{a_{I}},\,\frac{a_{t}}{a_{I}},\text{ and
}\frac{a_{\text{rh}}}{a_{I}}.$ (2.22)

However, one can connect one set of parameters to the other using Eqs. (2.12),
(2.15) and (2.16). For simplicity, in all numerical evaluations of this
section, we will fix $m_{I}=4\times 10^{13}$ GeV and $H_{I}=2\times 10^{13}$
GeV, as inspired by $\alpha$-attractor models [51, 52].

The fully numerical solution of the system of Boltzmann equations (2.2) and
(2.3), using the varying inflaton mass in Eq. (2.7) together with the
scattering rate in Eq. (2.10), is presented in Fig. 1 for $n=4$,
$m_{s}=10^{8}$ GeV, $\mu_{\phi}=10^{-2}$ GeV and $\mu_{h}=10^{6}$ GeV. The
left panel shows the energy densities for inflaton (blue) and SM radiation
(black), while the right panel shows the evolution of the SM bath temperature
$T$ as a function of the scale factor. The red dotted vertical lines
correspond to $a=a_{p}$, $a=a_{t}$, and $a=a_{\text{rh}}$, while the red
dotted horizontal lines (in the right panel) correspond to $T=T_{\text{max}}$,
$T=T_{t}$, and $T=T_{\text{rh}}$. Furthermore, the three dotted black lines in
the left panel show the analytical solutions in Eq. (2.1), which are in very
good agreement with the numerical solution.

Figure 1: Bosonic reheating. Left: Evolution of the energy densities for the
inflaton $\rho_{\phi}$ (solid blue) and the SM radiation $\rho_{R}$ (solid
black) as a function of the scale factor $a$. The three black dotted lines
correspond to the analytical solutions in Eq. (2.1). Right: Evolution of SM
temperature as a function of $a$. In both panels, $n=4$, $m_{s}=10^{8}$ GeV,
$\mu_{\phi}=10^{-2}$ GeV, $\mu_{h}=10^{6}$ GeV, $m_{I}=4\times 10^{13}$ GeV
and $H_{I}=2\times 10^{13}$ GeV were assumed. Here ‘before’, ‘pole’ and
‘after’ correspond to the epochs from the end of inflation until the pole
($a_{I}\leq a\leq a_{p}$), from the pole until the transition ($a_{p}\leq
a\leq a_{t}$), and from the transition to the end of reheating ($a_{t}\leq
a\leq a_{\text{rh}}$), respectively; cf. Eq. (2.1). The radiation-dominated
era starts from $a=a_{\text{rh}}$.

Figure 2: Bosonic reheating. Evolution of SM temperature $T$ as a function of
the scale factor $a$, for $n=4$ (solid), 6 (dashed), or 8 (dashed-dotted).
Left: $a_{p}/a_{I}=10^{5}$, $a_{t}/a_{I}=10^{10}$,
$a_{\text{rh}}/a_{I}=10^{15}$. Right: $m_{s}=10^{8}$ GeV, $\mu_{\phi}=10^{-2}$
GeV, and $\mu_{h}=10^{6}$ GeV. In both panels, $m_{I}=4\times 10^{13}$ GeV and
$H_{I}=2\times 10^{13}$ GeV.

In Fig. 2 the dependence of the SM temperature on $n$ is explored. We show the
evolution of $T$ as a function of the scale factor for $n=4$ (solid line),
$n=6$ (dashed line), or $n=8$ (dashed-dotted line). In the left frame, we fix
the scale factors $a_{p}/a_{I}=10^{5}$, $a_{t}/a_{I}=10^{10}$,
$a_{\text{rh}}/a_{I}=10^{15}$ (we therefore change the mediator mass and the
couplings) so that the different slopes of the SM temperature are clear. Their
behavior fits well with Eq. (2.20). In the right frame, we fix $m_{s}=10^{8}$
GeV, $\mu_{\phi}=10^{-2}$ GeV, and $\mu_{h}=10^{6}$ GeV. In this panel, it is
clear that, with an increase of $n$, the pole occurs earlier (as expected from
Eq. (2.12)), and consequently $T_{\text{max}}$ is further enhanced due to the
resonance.

In Fig. 3 we show the evolution of the inflaton (blue) and SM radiation
(black) energy densities for $n=4$ in the left panel, while the right panel
shows the variation of the SM bath temperature $T$ as a function of the scale
factor. For a mediator that is always heavier than the inflaton (solid black
lines), cf. $m_{s}\gg m_{I}$, no resonance develops during reheating.
Essentially, if the mediator is heavy, it can be integrated out, and the
reheating dynamics corresponds to the scenario in which the inflaton
annihilates through a contact interaction [10, 70, 11]. It corresponds to the
third case explored in Eqs. (2.1) and (2.20). Similarly, if the mediator is
always lighter than the inflaton (dashed black lines), the resonance cannot
develop either. However, it is interesting to note that in this case, with a
light mediator, the temperature increases monotonically until the end of the
reheating and consequently the reheating temperature
$T_{\text{rh}}=T_{\text{max}}$. It corresponds to the first case in Eqs. (2.1)
and (2.20). Note that the result presented in Fig. 1 corresponds to a
transition between the limit cases of the heavy mediator and the light
mediator.

Figure 3: Bosonic reheating. Left: Evolution of the energy densities for
inflaton $\rho_{\phi}$ (solid blue) and SM radiation $\rho_{R}$ (black) as a
function of the scale factor $a$. Right: Evolution of SM temperature as a
function of $a$. In both panels, $n=4$ and $a_{\text{rh}}/a_{I}=10^{10}$.
Solid black lines correspond to $m_{s}\gg m_{I}$ (heavy mediator), while
dashed black lines correspond to $m_{s}\ll m_{\phi}(a_{\text{rh}})$ (light
mediator). In both panels $m_{I}=4\times 10^{13}$ GeV and $H_{I}=2\times
10^{13}$ GeV.

###  2.2 Fermionic reheating

A successful reheating could also occur through inflaton annihilation into a
pair of vector-like fermions $\psi$ or e.g. right-handed neutrinos with mass
$m_{\psi}$, mediated by $S$ in the $s$ channel: $\phi\phi\to S\to\psi\psi$.
This can arise from the effective Lagrangian density

$\mathcal{L}\supset\frac{1}{2}\,\widetilde{\mu}_{\phi}\,\phi^{2}\,S+y_{f}\,S\,\bar{\psi}\,\psi\,.$
(2.23)

The corresponding annihilation rate reads

$\Gamma^{\phi\phi\to\psi\psi}\simeq\frac{\rho_{\phi}}{m_{\phi}}\frac{1}{8\pi}\frac{\mu_{\phi}^{2}\,y_{f}^{2}}{(4m_{\phi}^{2}-m_{s}^{2})^{2}+m_{s}^{2}\Gamma_{s}^{2}}\left(1-\frac{m_{\psi}^{2}}{m_{\phi}^{2}}\right)^{\frac{3}{2}},$
(2.24)

and the total decay width for $S$

$\Gamma_{s}=\frac{y_{f}^{2}\,m_{s}}{8\pi}\left(1-\frac{4m_{f}^{2}}{m_{s}^{2}}\right)^{\frac{3}{2}}\Theta(m_{s}-2m_{f})+\frac{1}{8\pi}\,\frac{\mu_{\phi}^{2}}{m_{s}}\sqrt{1-\frac{4\,m_{\phi}^{2}}{m_{s}^{2}}}\,\Theta(m_{s}-2m_{\phi})\,.$
(2.25)

The system of Boltzmann equations (2.2) and (2.3) can be analytically solved
following the same procedure as in the previous section. Consequently, the
solution for the SM radiation energy is given by

$\rho_{R}(a)\simeq\frac{9}{128\pi}\,\frac{H_{I}^{3}\,M_{P}^{4}\,y_{f}^{2}\,\mu_{\phi}^{2}}{m_{I}^{5}}\times\begin{dcases}\frac{n}{5n-11}\left(\frac{a}{a_{I}}\right)^{\frac{6(n-5)}{n+2}}\left[1-\left(\frac{a_{I}}{a}\right)^{\frac{2(5n-11)}{2+n}}\right]&\text{for
}a_{I}\leq a<a_{p},\\\
\frac{n}{n-1}\left(\frac{2m_{I}}{m_{s}}\right)^{4}\left(\frac{a_{I}}{a}\right)^{4}\left[1-\left(\frac{a_{I}}{a}\right)^{\frac{2(n-1)}{n+2}}\right]&\text{for
}a_{p}<a\,,\end{dcases}$ (2.26)

where $a_{p}$ is again the scale factor where the pole is reached (i.e.
$m_{\phi}(a_{p})=m_{s}/2$), given in Eq. (2.12), and corresponds to a maximal
temperature of the SM plasma

$T_{\text{max}}\simeq\left[\frac{135}{64\pi^{3}\,g_{\star}}\,\frac{n}{n-1}\,\frac{H_{I}^{3}\,M_{P}^{4}\,y_{f}^{2}\,\mu_{\phi}^{2}}{m_{I}^{5}}\left(\frac{2m_{I}}{m_{s}}\right)^{\frac{8(n-4)}{3(n-2)}}\right]^{\frac{1}{4}}\,.$
(2.27)

The analytical solution for the SM energy density in the fermionic case in Eq.
(2.26) can be divided into two cases (instead of three as in the bosonic case;
cf. Eq. (2.1)), as the scaling of $\rho_{R}$ during and after the pole is the
same, $\rho_{R}(a)\propto a^{-4}$.

The end of reheating, defined as the moment at which the Universe transitions
to a radiation-dominated epoch, happens at a scale factor

$a_{\text{rh}}\simeq
a_{I}\left(\frac{8\pi}{3}\,\frac{n-1}{n}\,\frac{m_{I}\,m_{s}^{4}}{H_{I}\,M_{P}^{2}\,y_{f}^{2}\,\mu_{\phi}^{2}}\right)^{\frac{n+2}{2(n-4)}},$
(2.28)

corresponding to the temperature

$T_{\text{rh}}\simeq\left[\frac{90}{\pi^{2}\,g_{\star}}\,H_{I}^{2}\,M_{P}^{2}\left(\frac{3}{8\pi}\,\frac{n}{n-1}\,\frac{H_{I}\,M_{P}^{2}\,y_{f}^{2}\,\mu_{\phi}^{2}}{m_{I}\,m_{s}^{4}}\right)^{\frac{3n}{n-4}}\right]^{\frac{1}{4}}\,.$
(2.29)

Finally, we observe from Eq. (2.26) that the SM temperature scales as

$T(a)\propto\begin{dcases}a^{\frac{3(n-5)}{2(n+2)}}&\text{for
}a_{I}<a<a_{p}\,,\\\ a^{-1}&\text{for }a_{p}<a\,.\end{dcases}$ (2.30)

The fully numerical solution of the system of Boltzmann equations for the
fermionic case is presented in Fig. 4 for $n=6$, $m_{s}=10^{8}$ GeV,
$\mu_{\phi}=10^{-2}$ GeV and $y_{f}=10^{-2}$. The left panel shows the energy
densities for the inflaton (blue) and SM radiation (black), while the right
panel shows the evolution of the SM bath temperature $T$ as a function of the
scale factor. The red dotted vertical lines correspond to $a=a_{p}$ and
$a=a_{\text{rh}}$, while the red dotted horizontal lines (in the right panel)
correspond to $T=T_{\text{max}}$ and $T=T_{\text{rh}}$. Furthermore, the two
dotted black lines in the left panel show the analytical solutions in Eq.
(2.26), which is in very good agreement with the numerical solution.

Figure 4: Fermionic reheating. Left: Evolution of the energy densities for the
inflaton $\rho_{\phi}$ (solid blue) and the SM radiation $\rho_{R}$ (solid
black) as a function of the scale factor $a$. The two black dotted lines
correspond to the analytical solutions in Eq. (2.26). Right: Evolution of the
SM temperature as a function of $a$. In both panels, $n=6$, $m_{s}=10^{8}$
GeV, $\mu_{\phi}=10^{-2}$ GeV, $y_{f}=10^{-3}$, $m_{I}=4\times 10^{13}$ GeV
and $H_{I}=2\times 10^{13}$ GeV were assumed.

##  3 Dark matter freeze-in during resonant reheating

The DM number can then be tracked by solving the Boltzmann equation for its
number density $n_{\text{dm}}$

$\frac{dn_{\text{dm}}}{dt}+3\,H\,n_{\text{dm}}=\gamma\,,$ (3.1)

where $\gamma$ is the DM production reaction rate density. As the SM entropy
is not conserved during reheating due to the annihilation of the inflaton in
SM particles, it is convenient to introduce a comoving number density
$N_{\text{dm}}\equiv n_{\text{dm}}\,a^{3}$, and therefore Eq. (3.1) can be
rewritten as

$\frac{dN_{\text{dm}}}{da}=\frac{a^{2}\,\gamma}{H}\,,$ (3.2)

which has to be numerically solved together with Eqs. (2.2) and (2.3), taking
the initial condition $N_{\text{dm}}(a_{I})=0$. To fit the whole observed DM
relic density, it is required that

$Y_{0}\,m_{\text{dm}}=\Omega
h^{2}\,\frac{1}{s_{0}}\,\frac{\rho_{c}}{h^{2}}\simeq 4.3\times
10^{-10}~{}\text{GeV},$ (3.3)

where $Y_{0}\equiv Y(T_{0})$ and $Y(T)\equiv n_{\text{dm}}(T)/s(T)$, with $s$
being the SM entropy density defined by

$s(T)=\frac{2\pi^{2}}{45}\,g_{\star s}(T)\,T^{3},$ (3.4)

and $g_{\star s}(T)$ being the number of relativistic degrees of freedom
contributing to the SM entropy. Furthermore, $\rho_{c}\simeq 1.05\times
10^{-5}\,h^{2}$ GeV/cm3 is the critical energy density, $s_{0}\simeq
2.69\times 10^{3}$ cm-3 the present entropy density [71], and $\Omega
h^{2}\simeq 0.12$ the observed abundance of DM relics [72].

In the early Universe, DM can be produced from 2-to-2 scattering of SM thermal
bath particles or inflatons, mediated by $S$.777We note that the inflaton
annihilation into a pair of $S$ particles than subsequently decay into DM
states is also viable, however it is subdominant in our scenario as discussed
in Appendix A. Additionally, it can also be mediated by the exchange of
gravitons. In the following subsections, the DM production will be discussed
in detail.

###  3.1 From the thermal bath

In this section, we investigate the impact of resonant reheating, presented in
the previous section, on dark matter (DM) yield. We consider the simplest
scenario, where the DM is a real singlet scalar $\varphi$ of mass
$m_{\text{dm}}$, with an interaction Lagrangian that reads

$\mathcal{L}_{\rm
dm}\supset\frac{1}{2}\,m_{\text{dm}}^{2}\,\varphi^{2}+\frac{1}{2}\,\mu_{\text{dm}}\,S\,\varphi^{2}\,,$
(3.5)

assuming that the only mediator between the SM, the inflaton and DM is $S$,
and again ignoring quartic couplings. Here, we consider any explicit coupling
between the DM and the inflaton, and quartic interactions between the mediator
$S$ or the SM Higgs doublet and the DM are absent by construction. We again
remain agnostic about the UV completion of such scenarios and focus mainly on
the aftermath of resonant reheating on DM production. The only interaction
that the DM can have with the visible sector is mediated by $S$, as can be
seen in Eq. (3.5). Therefore, DM can be produced in the early Universe from
the bath via the scattering of SM Higgs, mediated by $S$. This process can
occur during (inflaton domination) and after (radiation domination) reheating;
however, the effect of resonant reheating can only be seen in the former
scenario. Now, if the interaction strength between the DM and the visible
sector is sufficiently weak to ensure that the DM never thermalizes with the
radiation bath, then the production of DM takes place through the freeze-in
paradigm. This is the approach that follows.

As mentioned above, the DM yield is governed by the 2-to-2 scattering of the
bath particles (Higgs), mediated by the scalar $S$. The corresponding cross-
section for this process reads

$\sigma(s)_{hh\to\varphi\varphi}\simeq\frac{1}{8\,\pi\,s}\frac{\left(\mu_{\text{dm}}\,\mu_{h}\right)^{2}}{(s-m_{s}^{2})^{2}+\Gamma_{s}^{2}\,m_{s}^{2}}\sqrt{1-\frac{4m_{\text{dm}}^{2}}{s}}\,,$
(3.6)

where $\Gamma_{s}$ is given by (following Eqs. (2.9) and (3.5))

$\Gamma_{s}=\frac{1}{8\pi}\,\frac{\mu_{h}^{2}}{m_{s}}\sqrt{1-\frac{4m_{h}^{2}}{m_{s}^{2}}}+\frac{1}{8\pi}\,\frac{\mu_{\phi}^{2}}{m_{s}}\sqrt{1-\frac{4\,m_{\phi}^{2}}{m_{s}^{2}}}+\frac{1}{8\pi}\,\frac{\mu_{\text{dm}}^{2}}{m_{s}}\sqrt{1-\frac{4\,m_{\text{dm}}^{2}}{m_{s}^{2}}}\,,$
(3.7)

where the last term corresponds to $S$ decay into a pair of DM particles. In
the left panel of Fig. 5 the production cross section is shown as a function
of the center-of-mass energy $s$, for $m_{s}=10^{8}$ GeV,
$m_{\text{dm}}=10^{6}$ GeV, $\mu_{h}=10^{6}$ GeV and $\mu_{\phi}=10^{-2}$ GeV.
The resonance at $\sqrt{s}\simeq m_{s}$ is clearly visible as a peak in the
cross section.

Figure 5: Cross section (left) and production rate density (right), for
$m_{s}=10^{8}$ GeV, $m_{\text{dm}}=10^{6}$ GeV, $\mu_{h}=10^{6}$ GeV and
$\mu_{\phi}=10^{-2}$ GeV. The black dotted lines in the right panel correspond
to the analytical solution in Eq. (3.9).

One can further define the interaction rate density as

$\gamma=\gamma_{hh\to\varphi\varphi}=\frac{T}{2\pi^{4}}\,\int_{4m_{\text{dm}}^{2}}^{\infty}ds\,s^{3/2}\,\sigma\left(s\right)_{hh\to\varphi\varphi}\,K_{1}\left(\frac{\sqrt{s}}{T}\right),$
(3.8)

which takes the approximate form for $m_{\text{dm}}\ll m_{s}$

$\gamma_{hh\to\varphi\varphi}\simeq\frac{T}{2\pi^{4}}\,\frac{1}{8\,\pi}\frac{\left(\mu_{\text{dm}}\,\mu_{h}\right)^{2}}{1}\times\begin{dcases}\frac{4\,m_{\text{dm}}^{2}\,T}{m_{s}^{4}}\,K_{1}^{2}\left(\frac{m_{\text{dm}}}{T}\right)&\text{
for }T\ll\frac{m_{s}}{2},\\\
\frac{\pi}{\Gamma_{s}}\,K_{1}\left(\frac{m_{s}}{T}\right)\sqrt{1-\frac{4\,m_{\text{dm}}^{2}}{m_{s}^{2}}}&\text{
for }T\gg\frac{m_{s}}{2}\,,\end{dcases}$ (3.9)

or in the opposite case where $m_{s}\ll m_{\text{dm}}$

$\gamma_{hh\to\varphi\varphi}\simeq\frac{T}{2\pi^{4}}\,\frac{1}{8\,\pi}\frac{\left(\mu_{\text{dm}}\,\mu_{h}\right)^{2}}{1}\times\begin{dcases}\frac{\pi\,T^{2}}{8\,m_{\text{dm}}^{3}}\,e^{-\frac{2m_{\text{dm}}}{T}}&\text{
for }T\ll m_{\text{dm}},\\\ \frac{T}{6\,m_{\text{dm}}^{2}}&\text{ for }T\gg
m_{\text{dm}}\,,\end{dcases}$ (3.10)

with $K_{1}(x)$ being the modified Bessel function of the first kind. Here, we
have ignored the masses of the initial states, a reasonable approximation
during reheating, which occurs at a very high temperature, much before the
electroweak symmetry is broken. The right panel of Fig. 5 shows the production
rate density as a function of temperature $T$, for $m_{s}=10^{8}$ GeV,
$m_{\text{dm}}=10^{6}$ GeV, $\mu_{h}=10^{6}$ GeV and $\mu_{\phi}=10^{-2}$ GeV.
The black dotted lines correspond to the analytical solution in Eq. (3.9), in
good agreement with the fully numerical solution. Interestingly, the sharp
peak in the cross section appears as a smooth bump in the rate density.

In Fig. 6 examples of the evolution of the comoving DM abundance as a function
of the scale factor $a$ are shown for $n=4$, $a_{p}/a_{I}=10^{5}$,
$a_{t}/a_{I}=10^{10}$, $a_{\text{rh}}/a_{I}=10^{15}$, and
$m_{\text{dm}}=10^{-4}$ GeV and $\mu_{\text{dm}}=2\times 10^{-3}$ GeV (top
left), $m_{\text{dm}}=10^{0}$ GeV and $\mu_{\text{dm}}=5\times 10^{-2}$ GeV
(top right), or $m_{\text{dm}}=10^{4}$ GeV and $\mu_{\text{dm}}=3\times
10^{-1}$ GeV (bottom). In this section, we take $m_{I}=10^{12}$ GeV and
$H_{I}=5\times 10^{8}$ GeV, to avoid gravitational overproduction, as will be
seen in Section 3.3. The vertical red dotted lines correspond to the scale
factors for which $T=T_{\text{max}}$, $T=T_{t}$, or $T=T_{\text{rh}}$, while
the dashed red lines correspond to $T=m_{\text{dm}}$. All there benchmark
points fit well the observed relic abundance, represented by a horizontal blue
band. The detailed evolution of these benchmark points will be described in
the following.

Figure 6: DM production from the thermal bath. Evolution of the comoving
abundance of DM as a function of the scale factor $a$, for $n=4$,
$a_{p}/a_{I}=10^{5}$, $a_{t}/a_{I}=10^{10}$, $a_{\text{rh}}/a_{I}=10^{15}$,
and $m_{\text{dm}}=10^{-4}$ GeV and $\mu_{\text{dm}}=2\times 10^{-3}$ GeV (top
left), $m_{\text{dm}}=10^{0}$ GeV and $\mu_{\text{dm}}=5\times 10^{-2}$ GeV
(top right), or $m_{\text{dm}}=10^{4}$ GeV and $\mu_{\text{dm}}=3\times
10^{-1}$ GeV (bottom). In all panels $m_{I}=10^{12}$ GeV and $H_{I}=5\times
10^{8}$ GeV. The horizontal blue band corresponds to the abundance of DM
measured by Planck.

For light DM particles, with masses below $T_{\text{rh}}$ (and $m_{s}$), the
heavy mediator $S$ can be safely integrated out so that the cross section in
Eq. (3.6) effectively scales as
$\sigma(s)\propto\mu_{h}^{2}\mu_{\text{dm}}^{2}/(s\,m_{s}^{4})$, and therefore
the interaction rate density
$\gamma(T)\propto\mu_{h}^{2}\,\mu_{\text{dm}}^{2}\,T^{4}/m_{s}^{4}$; cf. the
first line of Eq. (3.9) in the limit $m_{\text{dm}}\ll T$. In this case, DM is
a standard IR FIMP, produced at low temperatures, mainly when $T\sim
m_{\text{dm}}$, as shown in the upper left panel of Fig. 6. The approximate
analytical solution of Eq. (3.2) reads

$Y_{0}\simeq\frac{135}{64\pi^{7}\,g_{\star
s}}\sqrt{\frac{10}{g_{\star}}}\,\frac{\mu_{\text{dm}}^{2}\,\mu_{h}^{2}\,M_{P}}{m_{\text{dm}}\,m_{s}^{4}}$
(3.11)

and, as expected, is independent on the reheating dynamics. Interestingly, to
fit the total observed relic abundance of DM, cf. Eq. (3.3), the required
inflaton coupling to DM is independent of the mass of DM. This can be seen in
Fig. 7, where the thick black line shows the parameter space in the plane
$[m_{\text{dm}},\,\mu_{\text{dm}}]$ that fits the entire observed abundance of
DM.

Figure 7: DM production from the thermal bath. Parameter space that fits the
entire observed DM abundance (black thick line), for $n=4$,
$a_{p}/a_{I}=10^{5}$, $a_{t}/a_{I}=10^{10}$, $a_{\text{rh}}/a_{I}=10^{15}$,
$m_{I}=10^{12}$ GeV and $H_{I}=5\times 10^{8}$ GeV. The black dotted lines
correspond to the analytical solutions described in the text.

The second case corresponds to intermediate masses of DM in the range
$T_{\text{rh}}\lesssim m_{\text{dm}}\lesssim T_{t}$. This case resembles the
previous one, with the difference that in the range $a_{t}<a<a_{\text{rh}}$
the temperature scales as $T(a)\propto a^{-\frac{9}{2(n+2)}}$, cf. Eq. (2.20).
DM is constantly produced in the range $T_{t}\gtrsim T\gtrsim m_{\text{dm}}$,
and stops when the temperature is of the order of its mass, due to kinematical
reasons. After production stops, during the range $m_{\text{dm}}\gtrsim
T\gtrsim T_{\text{rh}}$, the injection of entropy due to the decay of the
inflaton into the SM states dilutes the abundance of DM. In this regime, the
approximate analytical solution of Eq. (3.2) reads

$Y_{0}\simeq\frac{135}{8\pi^{8}\,g_{\star}}\sqrt{\frac{10}{g_{\star}}}\,\frac{\Gamma\left[\frac{6(n-2)}{n+2}\right]}{\Gamma\left[\frac{7n-10}{n+2}\right]}\frac{\mu_{\text{dm}}^{2}\,\mu_{h}^{2}\,M_{P}}{m_{s}^{4}\,T_{\text{rh}}}\left(\frac{T_{\text{rh}}}{m_{\text{dm}}}\right)^{\frac{4(n-2)}{3}},$
(3.12)

showing a strong dependence on the reheating dynamics, as expected in the case
of a UV FIMP. To fit the observed DM abundance, the inflaton coupling to DM
has to scale as

$\mu_{\text{dm}}\propto m_{\text{dm}}^{\frac{4n-11}{6}},$ (3.13)

as shown in Fig. 7 in the range $T_{\text{rh}}\lesssim m_{\text{dm}}\lesssim
T_{t}$.

The third case corresponds to heavy DM, with mass in the range $T_{t}\lesssim
m_{\text{dm}}\lesssim m_{s}$. The interaction rate density can be approximated
as
$\gamma(T)\propto\mu_{\text{dm}}^{2}\,m_{s}^{2}\,(T/m_{s})^{\frac{3}{2}}\,e^{-\frac{m_{s}}{T}}$,
cf. the second line of Eq. (3.9) in the limit $T\ll m_{s}$. Several comments
are in order: $i)$ the exponential factor stops DM production for temperatures
below the mass of $S$, $ii)$ since in this period the SM radiation behaves as
free radiation (cf. Eq. (2.20)), there is effectively no entropy dilution and,
therefore, the DM yield is constant, and $iii)$ there is, however, a dilution
of the abundance of DM in the era $T_{t}>T>T_{\text{rh}}$. An example of the
evolution of the DM yield for this case can be seen in the lower panel of Fig.
6. An approximate analytical solution of Eq. (3.2) is

$Y_{0}\simeq\frac{135}{4\sqrt{2}\,\pi^{\frac{11}{2}}\,g_{\star
s}}\sqrt{\frac{10}{g_{\star}}}\,\frac{\mu_{\text{dm}}^{2}\,m_{s}^{2}\,M_{P}}{T_{\text{rh}}^{5}}\left(\frac{a_{p}}{a_{\text{rh}}}\,\frac{T_{\text{max}}}{m_{s}}\right)^{\frac{6(n+1)}{n+2}}.$
(3.14)

We note that as the rate density is independent of the DM mass, it is required
that

$\mu_{\text{dm}}\propto m_{\text{dm}}^{-1/2}\,,$ (3.15)

to fit the observed DM abundance, as shown in Fig. 7

Finally, we also note that when $m_{\text{dm}}>m_{s}/2$, the DM production
becomes Boltzmann suppressed, requiring a sudden increase of
$\mu_{\text{dm}}$, as shown in Fig. 7. Overall, it is interesting to see that
this behavior significantly differs from the freeze-in during standard
reheating scenarios; see, e.g. Refs. [73, 74, 75, 76].

###  3.2 From inflaton annihilations

In our setup, another possible DM production channel occurs through the
inflaton annihilation via an $s$-channel exchange of the mediator $S$, i.e.
$\phi\phi\to S\to\varphi\varphi$. The production rate density is given by888We
note that a factor $1+\omega\equiv\frac{2n}{2+n}$ shall be included on the
right handed side of the Boltzmann equation for the DM number density.

$\gamma=\gamma_{\phi\phi\to\varphi\varphi}=+\frac{2\,n}{n+2}\,\frac{\Gamma_{\phi\phi\to\varphi\varphi}\,\rho_{\phi}}{m_{\phi}}\,,$
(3.16)

with

$\Gamma_{\phi\phi\to\varphi\varphi}\simeq\frac{\rho_{\phi}}{m_{\phi}}\,\frac{\mu_{\phi}^{2}\,\mu_{\text{dm}}^{2}}{32\pi\,m_{\phi}^{2}\left[(4m_{\phi}^{2}-m_{s}^{2})^{2}+\Gamma_{s}^{2}\,m_{s}^{2}\right]}\sqrt{1-\frac{m_{\text{dm}}^{2}}{m_{\phi}^{2}}}\,\Theta(m_{\phi}-m_{\text{dm}})\,.$
(3.17)

Figure 8: DM production from inflaton annihilations. Evolution of the comoving
DM abundance $N$ (left) and $Y$ (right) as a function of the scale factor $a$,
for $n=4$, $a_{p}/a_{I}=10^{5}$, $a_{t}/a_{I}=10^{10}$,
$a_{\text{rh}}/a_{I}=10^{15}$, and $m_{\text{dm}}=10^{-5}$ GeV and
$\mu_{\text{dm}}=2\times 10^{2}$ GeV (top), $m_{\text{dm}}=10^{-1}$ GeV and
$\mu_{\text{dm}}=25$ GeV (second from top), $m_{\text{dm}}=10^{4}$ GeV and
$\mu_{\text{dm}}=2$ GeV (third from top), or $m_{\text{dm}}=10^{9}$ GeV and
$\mu_{\text{dm}}=7\times 10^{5}$ GeV (bottom). In all panels $m_{I}=10^{12}$
GeV and $H_{I}=5\times 10^{8}$ GeV. The horizontal blue band corresponds to
the abundance of DM measured by Planck.

Examples of the evolution of the comoving DM abundance $N$ (left) and $Y$
(right) as a function of the scale factor $a$ are shown in Fig. 8 for $n=4$,
$a_{p}/a_{I}=10^{5}$, $a_{t}/a_{I}=10^{10}$, $a_{\text{rh}}/a_{I}=10^{15}$,
and $m_{\text{dm}}=10^{-5}$ GeV and $\mu_{\text{dm}}=2\times 10^{2}$ GeV
(top), $m_{\text{dm}}=10^{-1}$ GeV and $\mu_{\text{dm}}=25$ GeV (second from
top), $m_{\text{dm}}=10^{4}$ GeV and $\mu_{\text{dm}}=2$ GeV (third from top),
or $m_{\text{dm}}=10^{9}$ GeV and $\mu_{\text{dm}}=7\times 10^{5}$ GeV
(bottom). The vertical red dotted lines correspond to the scale factors for
which $T=T_{\text{max}}$, $T=T_{t}$, or $T=T_{\text{rh}}$, while the dashed
red lines correspond to $m_{\phi}=m_{\text{dm}}$. All four benchmark points
fit well the observed relic abundance, represented by a horizontal blue band
(right). The detailed evolution of these four benchmark points will be
described in the following.

In this scenario, the abundance of DM can be analytically estimated in
different regimes, depending on its mass with respect to the varying inflaton
mass. Light DM, with mass $m_{\text{dm}}<m_{\phi}(T_{\text{rh}})$, is produced
mainly at $T\simeq T_{\text{rh}}$ (top panel of Fig. 8). This corresponds to a
UV freeze-in during radiation dominance, with an interaction driven by a heavy
mediator. An approximate analytical solution of Eq. (3.2) is

$Y_{0}\simeq\frac{135}{64\pi^{3}\,g_{\star
s}}\,\frac{n}{n-3}\,\frac{\mu_{\text{dm}}^{2}\,\mu_{\phi}^{2}\,H_{I}^{3}\,M_{P}^{4}}{m_{I}^{4}\,m_{s}^{4}\,T_{\text{rh}}^{3}}\left(\frac{a_{I}}{a_{\text{rh}}}\right)^{\frac{3(8-n)}{2+n}}.$
(3.18)

To fit the total observed relic abundance of DM, cf. Eq. (3.3), the required
inflaton coupling to DM has to scale as

$\mu_{\text{dm}}\propto m_{\text{dm}}^{-\frac{1}{2}},$ (3.19)

as shown in Fig. 9 in the plane $[m_{\text{dm}},\,\mu_{\text{dm}}]$.

Figure 9: DM production from inflaton annihilations. Parameter space that fits
the entire observed DM abundance (black thick line), for $n=4$,
$a_{p}/a_{I}=10^{5}$, $a_{t}/a_{I}=10^{10}$, $a_{\text{rh}}/a_{I}=10^{15}$,
$m_{I}=10^{12}$ GeV and $H_{I}=5\times 10^{8}$ GeV. The black dotted lines
correspond to the analytical solutions described in the text.

The second case corresponds to DM with mass
$m_{\phi}(T_{\text{rh}})<m_{\text{dm}}<m_{\phi}(T_{t})$. The bulk of DM
density starts to be produced at $T=T_{t}$ and stops when
$m_{\phi}=m_{\text{dm}}$, due to kinematics, at a scale factor $a_{\text{dm}}$
given by

$a_{\text{dm}}\equiv
a_{I}\left(\frac{m_{I}}{m_{\text{dm}}}\right)^{\frac{n+2}{3(n-2)}},$ (3.20)

as shown in the second from top panel in Fig. 8. Between
$m_{\phi}=m_{\text{dm}}$ and $T=T_{\text{rh}}$ large entropy injection dilutes
the DM abundance; the final DM yield is estimated as

$Y_{0}\simeq\frac{135}{64\pi^{3}\,g_{\star
s}}\,\frac{n}{n-3}\,\frac{\mu_{\phi}^{2}\,\mu_{\text{dm}}^{2}\,M_{P}^{4}\,H_{I}^{3}}{m_{I}^{4}\,m_{s}^{4}\,T_{\text{rh}}^{3}}\left(\frac{a_{I}}{a_{\text{rh}}}\right)^{3}\left(\frac{m_{I}}{m_{\text{dm}}}\right)^{\frac{2(n-3)}{n-2}},$
(3.21)

which implies that, to fit the observed DM abundance, the inflaton coupling to
DM has to scale as

$\mu_{\text{dm}}\propto m_{\text{dm}}^{\frac{n-4}{2(n-2)}}.$ (3.22)

For the case $n=4$, $\mu_{\text{dm}}$ is independent of the DM mass, as shown
in Fig. 9.

The third case corresponds to DM with mass
$m_{\phi}(T_{t})<m_{\text{dm}}<m_{\phi}(T_{p})=2\,m_{s}$ produced mainly at
$m_{\text{dm}}=2\,m_{s}$, when $T=T_{\text{max}}$, which is a clear example of
a UV freeze-in during reheating; cf. the third panel from the top in Fig. 8.
Here, the interaction rate density in Eq. (3.17) is dominated by the pole. The
DM abundance is then diluted during the period $T_{t}>T>T_{\text{rh}}$, the
final yield being

$Y_{0}\simeq\frac{135}{16\pi^{2}\,g_{\star
s}}\,\frac{n}{n-2}\,\frac{\mu_{\phi}^{2}\,\mu_{\text{dm}}^{2}\,M_{P}^{4}\,H_{I}^{3}}{m_{I}^{2}\,m_{s}^{5}\,T_{\text{rh}}^{3}\,\Gamma_{s}}\left(\frac{a_{I}}{a_{\text{rh}}}\right)^{3}\left(\frac{m_{s}}{2m_{I}}\right)^{\frac{2}{n-2}},$
(3.23)

implying that

$\mu_{\text{dm}}\propto m_{\text{dm}}^{-\frac{1}{2}}$ (3.24)

to fit the observed DM abundance, as seen in Fig. 9.

Finally, the fourth case corresponds to the heavy DM with mass
$2\,m_{s}<m_{\text{dm}}<m_{I}$, which is produced mainly around
$m_{\text{dm}}=m_{\phi}$, when the scale factor is given by Eq. (3.20). The
abundance of DM suffers two periods of dilution, from $m_{\phi}=m_{\text{dm}}$
to $T=T_{\text{max}}$, and then during the interval $T_{t}>T>T_{\text{rh}}$;
see the lower panel of Fig. 8. The final DM yield is therefore

$Y_{0}\simeq\frac{135}{1024\pi^{3}\,g_{\star
s}}\,\frac{n}{3n-7}\left(\frac{a_{I}}{a_{\text{rh}}}\right)^{3}\frac{H_{I}^{3}\,M_{P}^{4}\,\mu_{\text{dm}}^{2}\,\mu_{\phi}^{2}}{m_{I}^{8}\,T_{\text{rh}}^{3}}\left(\frac{m_{I}}{m_{\text{dm}}}\right)^{\frac{2(3n-7)}{n-2}},$
(3.25)

implying that

$\mu_{\text{dm}}\propto m_{\text{dm}}^{\frac{5n-12}{2(n-2)}}$ (3.26)

to fit the observed DM abundance, as seen in Fig. 9. We note that DM particles
heavier than the inflaton ($m_{\text{dm}}>m_{I}$) are kinematically forbidden
in this channel.

###  3.3 Gravitational production

Gravity is the only interaction that is guaranteed to mediate between the DM
and the visible sector, and therefore it has to be taken into account. DM can
indeed be produced by the scattering of inflatons or SM particles, mediated by
the $s$-channel exchange of gravitons. Below we briefly discuss the DM
gravitational production.

#### From the thermal bath

DM is unavoidably produced from the UV freeze-in mechanism via 2-to-2
annihilation of the SM particles mediated by $s$-channel exchange of gravitons
[77, 78, 79]. The interaction rate density for such a process reads

$\gamma=\alpha\,\frac{T^{8}}{M_{P}^{4}}$ (3.27)

with $\alpha\simeq 1.9\times 10^{-4}$ for real scalar DM [80, 68, 17]. In the
cosmological scenario presented here, DM is produced mainly at
$T=T_{\text{max}}$ and then diluted during the period $T_{t}>T>T_{\text{rh}}$.
An approximate analytical solution of Eq. (3.2) is

$Y_{0}\simeq\frac{45\,\alpha}{4\pi^{2}\,g_{\star
s}}\,\frac{n+2}{n+5}\,\frac{T_{\text{max}}^{8}}{H_{I}\,M_{P}^{4}\,T_{\text{rh}}^{3}}\left(\frac{a_{I}}{a_{\text{rh}}}\right)^{3}\left(\frac{a_{p}}{a_{I}}\right)^{\frac{6(n+1)}{n+2}}.$
(3.28)

#### From inflaton annihilations

During reheating, the whole observed DM abundance can be generated via 2-to-2
annihilations of inflatons, mediated by the $s$-channel exchange of gravitons.
The interaction rate density for DM production reads [81, 82, 83]

$\gamma=\frac{2\,n}{n+2}\,\frac{\rho_{\phi}^{2}}{256\pi\,M_{P}^{4}}\,f\left(\frac{m_{\text{dm}}}{m_{\phi}}\right),$
(3.29)

with

$f(x)=\left(x^{2}+2\right)^{2}\sqrt{1-x^{2}}$ (3.30)

for real scalar DM [68]. The bulk of the DM relic density is promptly produced
at the beginning of the reheating era ($a=a_{I}$), and suffers two periods of
dilution: $a_{I}>a>a_{p}$ and $a_{t}>a>a_{\text{rh}}$. The final DM yield can
be analytically estimated as

$Y_{0}\simeq\frac{135}{256\pi^{3}\,g_{\star
s}}\,\frac{2n}{n-1}\left(\frac{H_{I}}{T_{\text{rh}}}\right)^{3}\left(\frac{a_{I}}{a_{\text{rh}}}\right)^{3}\left[1-\left(\frac{a_{I}}{a_{\text{rh}}}\right)^{\frac{6(n-1)}{n+2}}\right].$
(3.31)

Gravitation production from inflaton annihilation typically dominates over the
annihilation from thermal-bath particles and can even overtake production
through the exchange of $S$ for very large values of $H_{I}$ and
$m_{\text{dm}}$.

Figure 10: Parameter space that fits the whole observed DM abundance, for
$n=4$, $a_{p}/a_{I}=10^{5}$, $a_{t}/a_{I}=10^{10}$,
$a_{\text{rh}}/a_{I}=10^{15}$, $m_{I}=10^{12}$ GeV and $H_{I}=5\times 10^{8}$
GeV, assuming production through annihilations from SM particles in the
thermal bath (dashed blue), inflatons (dotted blue), gravitational production
(sharp cut at $m_{\text{dm}}\simeq 2\times 10^{8}$ GeV) and the combination
(thick black).

Figure 10 is an example of the combination of all previously described DM
production mechanisms, for $n=4$, $a_{p}/a_{I}=10^{5}$, $a_{t}/a_{I}=10^{10}$,
$a_{\text{rh}}/a_{I}=10^{15}$, $m_{I}=10^{12}$ GeV, and $H_{I}=5\times 10^{8}$
GeV. It shows the production through annihilations of SM particles in the
thermal bath (dashed blue; see Fig. 7), inflatons (dotted blue; see Fig. 9),
gravitational production (sharp cut at $m_{\text{dm}}\simeq 2\times 10^{8}$
GeV), and the combination of all of them (thick black). For this particular
choice of parameters, the production from the SM thermal bath dominates for
low masses. The production from inflaton scatterings is comparable near the
pole and dominates in a small region around $m_{\text{dm}}\sim 10^{8}$ GeV.
For higher masses the gravitational production from inflaton annihilations
dominates and overcloses the Universe for even higher DM masses. In summary,
in this section, we find that, within the framework of resonance reheating,
new features for DM UV freeze-in show up primarily due to the resonant
temperature behavior near the pole.

##  4 Primordial gravitational waves

The resonant reheating scenario described so far can have a possible signature
at the next-generation gravitational wave (GW) detectors, courtesy to the
blue-tilted primordial GWs with an inflationary origin. In this section we
discuss such a prospect.

The spectrum of GWs is described in terms of the fraction of their energy
density per logarithmic wavenumber (frequency) $k$ interval

$\Omega_{\text{GW}}(a,k)=\frac{1}{\rho_{c}}\,\frac{d\rho_{\text{GW}}}{d\log
k}\,,$ (4.1)

normalized to the critical density $\rho_{c}=3\,H^{2}\,M_{P}^{2}$. This
expression can be recasted as [84, 85, 86, 41]

$\Omega_{\text{GW}}(a\,,k)=\frac{1}{12}\left(\frac{k}{a\,H(a)}\right)^{2}\mathcal{P}_{T,\text{prim}}\,\mathcal{T}(a,k)\,,$
(4.2)

where the primordial tensor power spectrum $\mathcal{P}_{T,\text{prim}}$ can
be parametrized as [86, 41]

$\mathcal{P}_{T,\text{prim}}=r\,\mathcal{P}_{\zeta}(k_{\star})\left(\frac{k}{k_{\star}}\right)^{n_{T}}\,,$
(4.3)

with $k_{\star}=0.05~{}\text{Mpc}^{-1}$ being the Planck pivot scale,
$\mathcal{P}_{\zeta}(k_{\star})\simeq 2.1\times 10^{-9}$ [87] corresponding
power spectrum of the scalar perturbation at the pivot scale, $r$ is the
tensor-to-scalar ratio and $n_{T}$ represents the tensor spectral index. In
single-field inflationary models, $n_{T}\simeq-r/8$. Considering the current
bound on $r<0.036$ [67], we set $n_{T}=0$ throughout our analysis, which
corresponds to a scale-invariant primordial spectrum. In Eq. (4.2) the
transfer function $\mathcal{T}(a,\,k)$ connects primordial mode functions with
mode functions at some later time as [88, 86]

$\mathcal{T}(a,\,k)=\frac{1}{2}\left(\frac{a_{\text{hc}}}{a}\right)^{2},$
(4.4)

where the prefactor 1/2 appears due to oscillation-averaging the tensor mode
functions [86, 89, 90] and $a_{\text{hc}}$ is scale factor at the horizon
crossing defined as $a_{\text{hc}}\,H(a_{\text{hc}})=k$. Different $k$ modes
might reenter the horizon at different epochs, i.e., during radiation
domination or during reheating.

As is evident from Eq. (4.4), the transfer function characterizes the
expansion history between the moment of horizon crossing $a=a_{\text{hc}}$ of
a given mode $k$ and some later moment $a>a_{\text{hc}}$. From Eq. (4.2), one
can see that the spectral GW energy density at present (that is, at $a=a_{0}$)
reads

$\Omega_{\text{GW}}(k)\equiv\Omega_{\text{GW}}(a_{0},k)=\frac{1}{24}\left(\frac{k}{a_{0}\,H_{0}}\right)^{2}\mathcal{P}_{T,\text{prim}}\left(\frac{a_{\text{hc}}}{a_{0}}\right)^{2},$
(4.5)

where $H_{0}$ is the Hubble expansion rated as measured today. Depending on
the epoch of horizon re-entry for the perturbations, one can obtain the full
GW spectral energy density as

$\displaystyle\Omega_{\text{GW}}(a_{\text{hc}})$
$\displaystyle\simeq\Omega_{\gamma}^{(0)}\,\frac{\mathcal{P}_{T,\text{prim}}}{24}\,\frac{g_{\star}(T_{\text{hc}})}{2}\left(\frac{g_{\star
s}(T_{0})}{g_{\star s}(T_{\text{hc}})}\right)^{\frac{4}{3}}$
$\displaystyle\qquad\qquad\times\begin{dcases}\frac{g_{\star}(T_{\text{rh}})}{g_{\star}(T_{\text{hc}})}\left(\frac{g_{\star
s}(T_{\text{hc}})}{g_{\star
s}(T_{\text{rh}})}\right)^{\frac{4}{3}}\left(\frac{a_{\text{rh}}}{a_{\text{hc}}}\right)^{\frac{2(n-4)}{n+2}}&\text{for
}a_{I}<a_{\text{hc}}\leq a_{\text{rh}}\,,\\\ 1&\text{for }a_{\text{rh}}\leq
a_{\text{hc}}\leq a_{\rm eq}\,,\end{dcases}$ (4.6)

where $\Omega_{\gamma}^{(0)}\equiv\rho_{\gamma,0}/\rho_{c}=2.47\times
10^{-5}\,h^{-2}$ [86, 72] is the fraction of the energy density of photons at
the present epoch and $T_{\text{hc}}$ corresponds to the temperature at which
the corresponding mode reenters the horizon. The first line holds for the
modes that re-enter the horizon during reheating, while the second line
applies for the modes re-entering after reheating, during the radiation-
domination epoch. Note that, for $a_{\text{rh}}\leq a_{\text{hc}}\leq
a_{\text{eq}}$, we have assumed conservation of entropy from the moment of
horizon crossing until today, implying $T\propto a^{-1}\,g_{\star s}^{-1/3}$,
with $a=a_{\text{eq}}$ being the scale factor at the matter-radiation
equivalence at $T\equiv T_{\text{eq}}\simeq 0.7$ eV. The corresponding GW
frequency $f$ can be expressed as

$\displaystyle f(a_{\text{hc}})\equiv\frac{k}{2\pi\,a_{0}}=$
$\displaystyle\frac{1}{6}\sqrt{\frac{g_{\star}(T_{\text{rh}})}{10}}\left(\frac{g_{\star
s}(T_{0})}{g_{\star
s}(T_{\text{rh}})}\right)^{\frac{1}{3}}\frac{T_{0}\,T_{\text{rh}}}{M_{P}}$
$\displaystyle\quad\quad\times\begin{dcases}\left(\frac{a_{\text{rh}}}{a_{\text{hc}}}\right)^{\frac{2(n-1)}{n+2}}&\text{for
}a_{I}<a_{\text{hc}}\leq a_{\text{rh}}\,,\\\
\frac{a_{\text{rh}}}{a_{\text{hc}}}&\text{for }a_{\text{rh}}\leq
a_{\text{hc}}\leq a_{\rm eq}\,,\end{dcases}$ (4.7)

with $T_{0}\simeq 2.3\times 10^{-13}$ GeV [71]. Combining Eqs. (4) and (4) we
can rewrite the GR spectrum as a function of the frequency

$\displaystyle\Omega_{\text{GW}}(f)$
$\displaystyle\simeq\Omega_{\gamma}^{(0)}\,\frac{\mathcal{P}_{T,\text{prim}}}{24}\,\frac{g_{\star}(T_{\text{hc}})}{2}\left(\frac{g_{\star
s}(T_{0})}{g_{\star s}(T_{\text{hc}})}\right)^{\frac{4}{3}}$
$\displaystyle\times\begin{dcases}\frac{g_{\star}(T_{\text{rh}})}{g_{\star}(T_{\text{hc}})}\left(\frac{g_{\star
s}(T_{\text{hc}})}{g_{\star
s}(T_{\text{rh}})}\right)^{\frac{4}{3}}\left(\frac{f}{f_{\text{rh}}}\right)^{\frac{n-4}{n-1}}&\text{for
}f_{\text{rh}}\leq f<f_{\text{max}}\,,\\\ 1&\text{for }f_{\text{eq}}\leq f\leq
f_{\text{rh}}\,,\end{dcases}$ (4.8)

where $f_{\text{eq}}\equiv f(a_{\text{eq}})$ and

$f(a_{\text{rh}})\equiv\frac{1}{6}\sqrt{\frac{g_{\star}(T_{\text{rh}})}{10}}\left(\frac{g_{\star
s}(T_{0})}{g_{\star
s}(T_{\text{rh}})}\right)^{\frac{1}{3}}\frac{T_{0}T_{\text{rh}}}{M_{P}}\,.$
(4.9)

Finally, the GW spectrum presents a cutoff in frequency corresponding to
$a_{\text{hc}}=a_{I}$,

$f_{\text{max}}\equiv\frac{H_{I}}{2\,\pi}\,\frac{a_{I}}{a_{0}}=\frac{H_{I}}{2\,\pi}\left(\frac{g_{\star
s}(T_{0})}{g_{\star
s}(T_{\text{rh}})}\right)^{\frac{1}{3}}\frac{a_{I}}{a_{\text{rh}}}\,\frac{T_{0}}{T_{\text{rh}}}\,.$
(4.10)

From Eq. (4) one can see that modes that reenter the horizon after reheating
(that is $f<f_{\text{rh}}$) keep the same dependence with the frequency as the
original primordial spectrum $\mathcal{P}_{T,\text{prim}}$. However, if they
reenter during reheating ($f_{\text{rh}}<f<f_{\text{max}}$),
$\Omega_{\text{GW}}(f)$ gets a boost $\propto f^{\frac{n-4}{n-1}}$, which
implies that the spectrum becomes blue tilted with respect to the primordial
one. Interestingly, this only happens for $n>4$, as for $n=4$ the Universe
scales as free radiation during (and after) reheating.

Figure 11: Gravitational wave spectral energy density
$\Omega_{\text{GW}}h^{2}$ as function of the frequency $f$, and $n=4$, $n=6$
or $n=8$. The left panel corresponds to $a_{\text{rh}}/a_{I}=10^{3}$, while
the right panel corresponds to $a_{\text{rh}}/a_{I}=10^{10}$. In the figure
$H_{I}=2\times 10^{13}$ GeV and $r=0.035$. The yellow lines correspond to
projected sensitives of different detectors.

In Fig. 11, we show the GW spectral energy density at present epoch, as a
function of the frequency $f$ for $a_{\text{rh}}/a_{I}=10^{3}$ (left) and
$a_{\text{rh}}/a_{I}=10^{10}$ (right), and three choices of the inflaton
potential: $n=4$ (solid black), $n=6$ (dashed black) and $n=8$ (dashed-dotted
black). We also used $H_{I}=2\times 10^{13}$ GeV and the maximum value for the
tensor-to-scalar ratio $r=0.035$ [67]. We overlaid in yellow projected
sensitivity curves from the Big Bang Observer (BBO) [91, 92], ultimate DECIGO
[93, 94], LISA [95], $\mu$Ares [96], the cosmic explorer (CE) [97], the
Einstein Telescope (ET) [98, 99, 100, 101], and from resonant cavities [102],
see Ref. [103] for a recent review. In all cases, the GW energy density is
safe from the constraint on $\Delta N_{\text{eff}}$ [87], which quantifies the
extra radiation-like degrees of freedom around BBN. Note that for a larger
$a_{\text{rh}}/a_{I}$ (right panel) the GW spectrum receives a larger boost
because in that case the reheating lasts longer, resulting in a longer stiff
epoch that contributes to the blue tilt. As a consequence, the spectrum falls
within the sensitivity of not only the high-frequency detectors like resonant
cavity but also low-frequency detectors like BBO and DECIGO. We emphasize that
although such a spectrum is not unique to resonant reheating scenario (and can
be generic to any stiff background in the early Universe), any potential
signal of primordial GW either in the low or in the high frequencies can not
rule out the possibility of resonant production of radiation bath.

##  5 Conclusions

In this study, we investigate a reheating scenario where the transfer of
energy from the inflaton $\phi$ to the Standard Model (SM) particles proceeds
through an $s$-channel mediated scattering. We consider the inflaton to
oscillate around a generic monomial potential $\propto\phi^{n}$, with the
inflaton mass $m_{\phi}$ being field dependent and therefore varying over time
if $n>2$. If the mediator mass $m_{s}$ is of the order of the inflaton mass
$m_{\phi}$, a resonance can develop during reheating once the inflaton mass is
$m_{\phi}(\phi)\simeq m_{s}/2$. We validate these characteristics by solving a
system of coupled Boltzmann equations for the inflaton and SM energy densities
(Eqs. (2.2) and (2.3)) and observe a characteristic bump occurring both in the
radiation and in the temperature profiles, as depicted in Fig. 1. Furthermore,
we derive several analytical approximations, detailed in Eq. (2.1), which
nicely fits the numerical results. Our investigation reveals that resonance
reheating exhibits more nuanced temperature scaling compared to conventional
scenarios presented in the existing literature, which is typically
characterized by a featureless power law, as given in Eq. (2.20). In scenarios
where the mediator is always heavier than the inflaton, precluding the
development of resonance, our model successfully reproduces the results of the
existing literature. Alternatively, when the mediator is always lighter than
the inflaton, the SM temperature rises throughout the reheating phase, with
the reheating temperature serving as the maximum temperature achieved; see
Fig. 3.

Some implications for particle phenomenology and cosmology during resonant
reheating were studied. First, the mediator could play a pivotal role in
producing not only the SM states but also particles beyond the SM, such as
dark matter (DM). If the interaction rates between the DM and the visible
sector are very suppressed, DM could be produced during reheating through the
annihilation of a pair of SM or inflaton particles. Similarly to radiation, DM
production also has a very rich freeze-in phenomenology, as illustrated in
Figs. 6, 8, and 10.

Second, we also proposed the prospect of finding a signature of this model via
the primordial gravitational wave (GW) spectrum. For $n>4$, the inflationary
GW spectrum experiences a significant blue tilt within the frequency range
corresponding to the modes that cross the horizon during the stiff period. The
resulting spectrum could be well within the reach of several next-generation
low- and high-frequency GW detectors, as shown in Fig. 11.

## Acknowledgments

NB received funding from the Spanish FEDER / MCIU-AEI under the grant
FPA2017-84543-P. YX acknowledges the support from the Cluster of Excellence
“Precision Physics, Fundamental Interactions, and Structure of Matter”
(PRISMA+ EXC 2118/1) funded by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) within the German Excellence Strategy (Project No.
390831469).

##  Appendix A Annihilation to mediators

In our setup, the annihilation $\phi\phi\to ss$ through a $t$\- and
$u$-channel exchange of a $\phi$ is also kinematically open. We assume that
$s$ particles have a sufficiently large decay width and rapidly decay into SM
radiation once they are produced. The total production rate considering the
$t$ and $u$ channels is

$\Gamma^{\phi\phi\to
ss}\simeq\frac{\rho_{\phi}}{m_{\phi}}\,\frac{1}{16\pi\,m_{\phi}^{2}}\left(\frac{\mu_{\phi}^{2}}{m_{s}^{2}-2m_{\phi}^{2}}\right)^{2}\sqrt{1-\frac{m_{s}^{2}}{m_{\phi}^{2}}}\,.$
(A.1)

The solution for the radiation energy density is

$\rho_{R}(a)\simeq\frac{9}{64\pi}\,\frac{n}{8n-17}\left(\frac{a}{a_{I}}\right)^{\frac{6(2n-7)}{2+n}}\left[1-\left(\frac{a_{I}}{a}\right)^{\frac{2(8n-17)}{2+n}}\right]\frac{H_{I}^{3}\,M_{P}^{4}\,\mu_{\phi}^{4}}{M_{I}^{7}}$
(A.2)

in the limit $m_{s}\ll m_{\phi}$. Note that if $\mu_{h}=\mu_{\phi}$, Eq. (A.2)
and the first line of Eq. (2.1) are the same (up to a factor 8), making this
channel potentially important. However, if one demands an extended reheating
period, a long-lived inflaton is desirable, and hence a suppressed
$\mu_{\phi}$. Under these circumstances, the inflaton annihilation into
mediators is suppressed. We have checked that in all the benchmarks used in
this work, this channel is completely subdominant.

## References

  * [1] A.A. Starobinsky, _A New Type of Isotropic Cosmological Models Without Singularity_ , _Phys. Lett. B_ 91 (1980) 99.
  * [2] A.H. Guth, _The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems_ , _Phys. Rev. D_ 23 (1981) 347.
  * [3] A.D. Linde, _A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems_ , _Phys. Lett. B_ 108 (1982) 389.
  * [4] A. Albrecht and P.J. Steinhardt, _Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking_ , _Phys. Rev. Lett._ 48 (1982) 1220.
  * [5] L.F. Abbott, E. Farhi and M.B. Wise, _Particle Production in the New Inflationary Cosmology_ , _Phys. Lett. B_ 117 (1982) 29.
  * [6] L. Kofman, A.D. Linde and A.A. Starobinsky, _Reheating after inflation_ , _Phys. Rev. Lett._ 73 (1994) 3195 [hep-th/9405187].
  * [7] L. Kofman, A.D. Linde and A.A. Starobinsky, _Towards the theory of reheating after inflation_ , _Phys. Rev. D_ 56 (1997) 3258 [hep-ph/9704452].
  * [8] J.F. Dufaux, G.N. Felder, L. Kofman, M. Peloso and D. Podolsky, _Preheating with trilinear interactions: Tachyonic resonance_ , _JCAP_ 07 (2006) 006 [hep-ph/0602144].
  * [9] K.D. Lozanov, _Lectures on Reheating after Inflation_ , 1907.04402.
  * [10] M.A.G. García, K. Kaneta, Y. Mambrini and K.A. Olive, _Inflaton Oscillations and Post-Inflationary Reheating_ , _JCAP_ 04 (2021) 012 [2012.10756].
  * [11] N. Bernal, S. Cléry, Y. Mambrini and Y. Xu, _Probing reheating with graviton bremsstrahlung_ , _JCAP_ 01 (2024) 065 [2311.12694].
  * [12] M.A.G. García, K. Kaneta, Y. Mambrini and K.A. Olive, _Reheating and Post-inflationary Production of Dark Matter_ , _Phys. Rev. D_ 101 (2020) 123507 [2004.08404].
  * [13] R.T. Co, E. González and K. Harigaya, _Increasing Temperature toward the Completion of Reheating_ , _JCAP_ 11 (2020) 038 [2007.04328].
  * [14] A. Ahmed, B. Grzadkowski and A. Socha, _Implications of time-dependent inflaton decay on reheating and dark matter production_ , _Phys. Lett. B_ 831 (2022) 137201 [2111.06065].
  * [15] P. Arias, N. Bernal, J.K. Osiński and L. Roszkowski, _Dark matter axions in the early universe with a period of increasing temperature_ , _JCAP_ 05 (2023) 028 [2207.07677].
  * [16] B. Barman, N. Bernal, Y. Xu and Ó. Zapata, _Ultraviolet freeze-in with a time-dependent inflaton decay_ , _JCAP_ 07 (2022) 019 [2202.12906].
  * [17] S. Cléry, Y. Mambrini, K.A. Olive and S. Verner, _Gravitational portals in the early Universe_ , _Phys. Rev. D_ 105 (2022) 075005 [2112.15214].
  * [18] M.R. Haque and D. Maity, _Gravitational reheating_ , _Phys. Rev. D_ 107 (2023) 043531 [2201.02348].
  * [19] R.T. Co, Y. Mambrini and K.A. Olive, _Inflationary gravitational leptogenesis_ , _Phys. Rev. D_ 106 (2022) 075006 [2205.01689].
  * [20] B. Barman, S. Cléry, R.T. Co, Y. Mambrini and K.A. Olive, _Gravity as a portal to reheating, leptogenesis and dark matter_ , _JHEP_ 12 (2022) 072 [2210.05716].
  * [21] M.R. Haque, D. Maity and R. Mondal, _Gravitational neutrino reheating_ , _Phys. Rev. D_ 109 (2024) 063543 [2311.07684].
  * [22] S. Cléry, Y. Mambrini, K.A. Olive, A. Shkerin and S. Verner, _Gravitational portals with nonminimal couplings_ , _Phys. Rev. D_ 105 (2022) 095042 [2203.02004].
  * [23] B. Barman, N. Bernal and J. Rubio, _Rescuing Gravitational-Reheating in Chaotic Inflation_ , 2310.06039.
  * [24] J. McDonald, _Thermally generated gauge singlet scalars as selfinteracting dark matter_ , _Phys. Rev. Lett._ 88 (2002) 091304 [hep-ph/0106249].
  * [25] K.-Y. Choi and L. Roszkowski, _E-WIMPs_ , _AIP Conf. Proc._ 805 (2005) 30 [hep-ph/0511003].
  * [26] A. Kusenko, _Sterile neutrinos, dark matter, and the pulsar velocities in models with a Higgs singlet_ , _Phys. Rev. Lett._ 97 (2006) 241301 [hep-ph/0609081].
  * [27] K. Petraki and A. Kusenko, _Dark-matter sterile neutrinos in models with a gauge singlet in the Higgs sector_ , _Phys. Rev. D_ 77 (2008) 065014 [0711.4646].
  * [28] L.J. Hall, K. Jedamzik, J. March-Russell and S.M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080 [0911.1120].
  * [29] N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen, _The Dawn of FIMP Dark Matter: A Review of Models and Constraints_ , _Int. J. Mod. Phys. A_ 32 (2017) 1730023 [1706.07442].
  * [30] F. Elahi, C. Kolda and J. Unwin, _UltraViolet Freeze-in_ , _JHEP_ 03 (2015) 048 [1410.6157].
  * [31] M.A.G. García, Y. Mambrini, K.A. Olive and M. Peloso, _Enhancement of the Dark Matter Abundance Before Reheating: Applications to Gravitino Dark Matter_ , _Phys. Rev. D_ 96 (2017) 103510 [1709.01549].
  * [32] N. Bernal, F. Elahi, C. Maldonado and J. Unwin, _Ultraviolet Freeze-in and Non-Standard Cosmologies_ , _JCAP_ 11 (2019) 026 [1909.07992].
  * [33] R. Allahverdi et al., _The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe_ , _Open J.Astrophys._ 4 (2021)  [2006.16182].
  * [34] M. Giovannini, _Gravitational waves constraints on postinflationary phases stiffer than radiation_ , _Phys. Rev. D_ 58 (1998) 083504 [hep-ph/9806329].
  * [35] A. Riazuelo and J.-P. Uzan, _Quintessence and gravitational waves_ , _Phys. Rev. D_ 62 (2000) 083506 [astro-ph/0004156].
  * [36] N. Seto and J. Yokoyama, _Probing the equation of state of the early universe with a space laser interferometer_ , _J. Phys. Soc. Jap._ 72 (2003) 3082 [gr-qc/0305096].
  * [37] L.A. Boyle and A. Buonanno, _Relating gravitational wave constraints from primordial nucleosynthesis, pulsar timing, laser interferometers, and the CMB: Implications for the early Universe_ , _Phys. Rev. D_ 78 (2008) 043531 [0708.2279].
  * [38] A. Stewart and R. Brandenberger, _Observational Constraints on Theories with a Blue Spectrum of Tensor Modes_ , _JCAP_ 08 (2008) 012 [0711.4602].
  * [39] B. Li and P.R. Shapiro, _Precision cosmology and the stiff-amplified gravitational-wave background from inflation: NANOGrav, Advanced LIGO-Virgo and the Hubble tension_ , _JCAP_ 10 (2021) 024 [2107.12229].
  * [40] M. Artymowski, O. Czerwińska, Z. Lalak and M. Lewicki, _Gravitational wave signals and cosmological consequences of gravitational reheating_ , _JCAP_ 04 (2018) 046 [1711.08473].
  * [41] C. Caprini and D.G. Figueroa, _Cosmological Backgrounds of Gravitational Waves_ , _Class. Quant. Grav._ 35 (2018) 163001 [1801.04268].
  * [42] N. Bernal and F. Hajkarim, _Primordial Gravitational Waves in Nonstandard Cosmologies_ , _Phys. Rev. D_ 100 (2019) 063502 [1905.10410].
  * [43] D.G. Figueroa and E.H. Tanin, _Ability of LIGO and LISA to probe the equation of state of the early Universe_ , _JCAP_ 08 (2019) 011 [1905.11960].
  * [44] T. Opferkuch, P. Schwaller and B.A. Stefanek, _Ricci Reheating_ , _JCAP_ 07 (2019) 016 [1905.06823].
  * [45] N. Bernal, A. Ghoshal, F. Hajkarim and G. Lambiase, _Primordial Gravitational Wave Signals in Modified Cosmologies_ , _JCAP_ 11 (2020) 051 [2008.04959].
  * [46] Y. Gouttenoire, G. Servant and P. Simakachorn, _Kination cosmology from scalar fields and gravitational-wave signatures_ , 2111.01150.
  * [47] R. Caldwell et al., _Detection of early-universe gravitational-wave signatures and fundamental physics_ , _Gen. Rel. Grav._ 54 (2022) 156 [2203.07972].
  * [48] A. Chakraborty, M.R. Haque, D. Maity and R. Mondal, _Inflaton phenomenology via reheating in light of primordial gravitational waves and the latest BICEP/Keck data_ , _Phys. Rev. D_ 108 (2023) 023515 [2304.13637].
  * [49] B. Barman, A. Ghoshal, B. Grzadkowski and A. Socha, _Measuring inflaton couplings via primordial gravitational waves_ , _JHEP_ 07 (2023) 231 [2305.00027].
  * [50] B. Barman, S. Jyoti Das, M.R. Haque and Y. Mambrini, _A Gravitational Ménage à Trois: leptogenesis, primordial gravitational wave & PBH-induced reheating_, 2403.05626.
  * [51] R. Kallosh and A. Linde, _Universality Class in Conformal Inflation_ , _JCAP_ 07 (2013) 002 [1306.5220].
  * [52] R. Kallosh, A. Linde and D. Roest, _Superconformal Inflationary $\alpha$-Attractors_, _JHEP_ 11 (2013) 198 [1311.0472].
  * [53] A.A. Starobinsky, _Nonsingular Model of the Universe with the Quantum Ggravitational de Sitter Stage and its Observational Consequences_ , in _Second Seminar on Quantum Gravity_ , 1981.
  * [54] A.A. Starobinsky, _The Perturbation Spectrum Evolving from a Nonsingular Initially De-Sitter Cosmology and the Microwave Background Anisotropy_ , _Sov. Astron. Lett._ 9 (1983) 302.
  * [55] L.A. Kofman, A.D. Linde and A.A. Starobinsky, _Inflationary Universe Generated by the Combined Action of a Scalar Field and Gravitational Vacuum Polarization_ , _Phys. Lett. B_ 157 (1985) 361.
  * [56] M.S. Turner, _Coherent Scalar Field Oscillations in an Expanding Universe_ , _Phys. Rev. D_ 28 (1983) 1243.
  * [57] J. Ellis, M.A.G. García, D.V. Nanopoulos, K.A. Olive and M. Peloso, _Post-Inflationary Gravitino Production Revisited_ , _JCAP_ 03 (2016) 008 [1512.05701].
  * [58] K. Harigaya and K. Mukaida, _Thermalization after/during Reheating_ , _JHEP_ 05 (2014) 006 [1312.3097].
  * [59] M.A.G. García and M.A. Amin, _Prethermalization production of dark matter_ , _Phys. Rev. D_ 98 (2018) 103504 [1806.01865].
  * [60] M. Drees and B. Najjari, _Energy spectrum of thermalizing high energy decay products in the early universe_ , _JCAP_ 10 (2021) 009 [2105.01935].
  * [61] M. Drees and B. Najjari, _Multi-species thermalization cascade of energetic particles in the early universe_ , _JCAP_ 08 (2023) 037 [2205.07741].
  * [62] N. Bernal and Y. Xu, _WIMPs during reheating_ , _JCAP_ 12 (2022) 017 [2209.07546].
  * [63] Y. Shtanov, J.H. Traschen and R.H. Brandenberger, _Universe reheating after inflation_ , _Phys. Rev. D_ 51 (1995) 5438 [hep-ph/9407247].
  * [64] K. Ichikawa, T. Suyama, T. Takahashi and M. Yamaguchi, _Primordial Curvature Fluctuation and Its Non-Gaussianity in Models with Modulated Reheating_ , _Phys. Rev. D_ 78 (2008) 063545 [0807.3988].
  * [65] K. Kainulainen, S. Nurmi, T. Tenkanen, K. Tuominen and V. Vaskonen, _Isocurvature Constraints on Portal Couplings_ , _JCAP_ 06 (2016) 022 [1601.07733].
  * [66] B. Barman, N. Bernal, Y. Xu and Ó. Zapata, _Bremsstrahlung-induced gravitational waves in monomial potentials during reheating_ , _Phys. Rev. D_ 108 (2023) 083524 [2305.16388].
  * [67] BICEP, Keck collaboration, _Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season_ , _Phys. Rev. Lett._ 127 (2021) 151301 [2110.00483].
  * [68] B. Barman and N. Bernal, _Gravitational SIMPs_ , _JCAP_ 06 (2021) 011 [2104.10699].
  * [69] G.F. Giudice, E.W. Kolb and A. Riotto, _Largest temperature of the radiation era and its cosmological implications_ , _Phys. Rev. D_ 64 (2001) 023508 [hep-ph/0005123].
  * [70] Y. Xu, _Constraining axion and ALP dark matter from misalignment during reheating_ , _Phys. Rev. D_ 108 (2023) 083536 [2308.15322].
  * [71] Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * [72] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [73] J. Silva-Malpartida, N. Bernal, J. Jones-Pérez and R.A. Lineros, _From WIMPs to FIMPs with low reheating temperatures_ , _JCAP_ 09 (2023) 015 [2306.14943].
  * [74] M. Becker, E. Copello, J. Harz, J. Lang and Y. Xu, _Confronting dark matter freeze-in during reheating with constraints from inflation_ , _JCAP_ 01 (2024) 053 [2306.17238].
  * [75] C. Cosme, F. Costa and O. Lebedev, _Freeze-in at stronger coupling_ , 2306.13061.
  * [76] C. Cosme, F. Costa and O. Lebedev, _Temperature evolution in the Early Universe and freeze-in at stronger coupling_ , 2402.04743.
  * [77] M. Garny, M. Sandora and M.S. Sloth, _Planckian Interacting Massive Particles as Dark Matter_ , _Phys. Rev. Lett._ 116 (2016) 101302 [1511.03278].
  * [78] M. Garny, A. Palessandro, M. Sandora and M.S. Sloth, _Theory and Phenomenology of Planckian Interacting Massive Particles as Dark Matter_ , _JCAP_ 02 (2018) 027 [1709.09688].
  * [79] Y. Tang and Y.-L. Wu, _On Thermal Gravitational Contribution to Particle Production and Dark Matter_ , _Phys. Lett. B_ 774 (2017) 676 [1708.05138].
  * [80] N. Bernal, M. Dutra, Y. Mambrini, K. Olive, M. Peloso and M. Pierre, _Spin-2 Portal Dark Matter_ , _Phys. Rev. D_ 97 (2018) 115020 [1803.01866].
  * [81] Y. Mambrini and K.A. Olive, _Gravitational Production of Dark Matter during Reheating_ , _Phys. Rev. D_ 103 (2021) 115009 [2102.06214].
  * [82] N. Bernal and C.S. Fong, _Dark matter and leptogenesis from gravitational production_ , _JCAP_ 06 (2021) 028 [2103.06896].
  * [83] N. Bernal, J. Harz, M. Mojahed and Y. Xu, _Dark Matter Production in Large Field Inflation_ , 240?.soon.
  * [84] M. Maggiore, _Gravitational wave experiments and early universe cosmology_ , _Phys. Rept._ 331 (2000) 283 [gr-qc/9909001].
  * [85] Y. Watanabe and E. Komatsu, _Improved Calculation of the Primordial Gravitational Wave Spectrum in the Standard Model_ , _Phys. Rev. D_ 73 (2006) 123515 [astro-ph/0604176].
  * [86] K. Saikawa and S. Shirai, _Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model_ , _JCAP_ 05 (2018) 035 [1803.01038].
  * [87] Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [88] L.A. Boyle and P.J. Steinhardt, _Probing the early universe with inflationary gravitational waves_ , _Phys. Rev. D_ 77 (2008) 063504 [astro-ph/0512014].
  * [89] D.G. Figueroa and E.H. Tanin, _Inconsistency of an inflationary sector coupled only to Einstein gravity_ , _JCAP_ 10 (2019) 050 [1811.04093].
  * [90] G. Choi, R. Jinno and T.T. Yanagida, _Probing PeV scale SUSY breaking with satellite galaxies and primordial gravitational waves_ , _Phys. Rev. D_ 104 (2021) 095018 [2107.12804].
  * [91] J. Crowder and N.J. Cornish, _Beyond LISA: Exploring future gravitational wave missions_ , _Phys. Rev. D_ 72 (2005) 083005 [gr-qc/0506015].
  * [92] V. Corbin and N.J. Cornish, _Detecting the cosmic gravitational wave background with the big bang observer_ , _Class. Quant. Grav._ 23 (2006) 2435 [gr-qc/0512039].
  * [93] N. Seto, S. Kawamura and T. Nakamura, _Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space_ , _Phys. Rev. Lett._ 87 (2001) 221103 [astro-ph/0108011].
  * [94] H. Kudoh, A. Taruya, T. Hiramatsu and Y. Himemoto, _Detecting a gravitational-wave background with next-generation space interferometers_ , _Phys. Rev. D_ 73 (2006) 064006 [gr-qc/0511145].
  * [95] LISA collaboration, _Laser Interferometer Space Antenna_ , 1702.00786.
  * [96] A. Sesana et al., _Unveiling the gravitational universe at $\mu$-Hz frequencies_, _Exper. Astron._ 51 (2021) 1333 [1908.11391].
  * [97] D. Reitze et al., _Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO_ , _Bull. Am. Astron. Soc._ 51 (2019) 035 [1907.04833].
  * [98] S. Hild et al., _Sensitivity Studies for Third-Generation Gravitational Wave Observatories_ , _Class. Quant. Grav._ 28 (2011) 094013 [1012.0908].
  * [99] M. Punturo et al., _The Einstein Telescope: A third-generation gravitational wave observatory_ , _Class. Quant. Grav._ 27 (2010) 194002.
  * [100] B. Sathyaprakash et al., _Scientific Objectives of Einstein Telescope_ , _Class. Quant. Grav._ 29 (2012) 124013 [1206.0331].
  * [101] M. Maggiore et al., _Science Case for the Einstein Telescope_ , _JCAP_ 03 (2020) 050 [1912.02622].
  * [102] N. Herman, L. Lehoucq and A. Fúzfa, _Electromagnetic antennas for the resonant detection of the stochastic gravitational wave background_ , _Phys. Rev. D_ 108 (2023) 124009 [2203.15668].
  * [103] N. Aggarwal et al., _Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies_ , _Living Rev. Rel._ 24 (2021) 4 [2011.12414].