el Ahmed0000-0002-2907-2433

dan Grzadkowski0000-0001-9980-6335

a Socha0000-0002-4924-9267

# Higgs boson induced reheating and ultraviolet frozen-in dark matter

Aqeel Ahmed  Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117
Heidelberg, Germany  Bohdan Grzadkowski  Anna Socha  Faculty of Physics,
University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland

###### Abstract

A reheating phase in the early universe is an essential part of all
inflationary models during which not only the Standard Model (SM) quanta are
produced but it can also shed light on the production of dark matter. In this
work, we explore a class of reheating models where the reheating is induced by
a cubic interaction of the inflaton $\phi$ to the SM Higgs boson ${\bm{h}}$ of
the form $g_{h\phi}M_{\rm Pl}\phi|{\bm{h}}|^{2}$ adopting the
$\alpha$-attractor T-model of inflation. Assuming inflaton as a background
field such interaction implies a $\phi$-dependent mass term of the Higgs boson
and a non-trivial phase-space suppression of the reheating efficiency. As a
consequence, the reheating is prolonged and the maximal temperature of the SM
thermal bath is reduced. In particular, due to oscillations of the inflaton
field the $\phi$-dependent Higgs boson mass results in periodic transitions
between phases of broken and unbroken electroweak gauge symmetry. The
consequences of these rapid phase transitions have been studied in detail. A
purely gravitational reheating mechanism in the presence of the inflaton
background, i.e., for $g_{h\phi}=0$, has also been investigated. It turned out
that even though it may account for the total production of SM radiation in
the absence of $g_{h\phi}$, its contribution to the reheating is subdominant
for the range of $g_{h\phi}$ considered in this work. Approximate analytical
solutions of Boltzmann equations for energy densities of the inflaton and SM
radiation have been obtained. As a dark matter candidate a massive Abelian
vector boson, $X_{\mu}$, has been considered. Various production mechanisms of
$X_{\mu}$ have been discussed including (i) purely gravitational production
from the inflaton background, (ii) gravitational freeze-in from the SM quanta,
(iii) inflaton decay through a dim-5 effective operator, and (iv) Higgs portal
freeze-in and Higgs decay through a dim-6 effective operator. Parameters that
properly describe the observed relic abundance have been determined.

###### Contents

  1. I intro
  2. II model
  3. III inflaton
  4. IV reheating
    1. IV.1 Evolution of the inflaton and SM energy densities
    2. IV.2 Higgs boson production from oscillating inflaton
    3. IV.3 Reheating through massless Higgs bosons
    4. IV.4 Reheating through massive Higgs bosons
    5. IV.5 Tachyonic resonant production of Higgs bosons
    6. IV.6 Numerical analysis
  5. V DM-effects
    1. V.1 Inflaton induced gravitational DM production
    2. V.2 Gravitational freeze-in production
    3. V.3 Inflaton decay
    4. V.4 Higgs boson decay
    5. V.5 Higgs portal freeze-in production
  6. VI summary
  7. A constraints
  8. B RDMproduction
    1. B.1 S-matrix for the SM and DM production
    2. B.2 Production of the SM and DM from the inflaton
    3. B.3 BeqnPhi

##  I Introduction

The most successful theory of the early Universe is cosmic inflation which
results in a period of exponential expansion [1, 2]. The theory of inflation
can be effectively described by a slowly-rolling single scalar field, called
the inflaton $\phi$, with an approximately flat potential. Inflation
successfully explains puzzles of the early Universe, e.g., the horizon
problem, flatness problem, and seeds for Large Scale Structures (for a review
see [3]). Many of the non-trivial features of the inflationary paradigm can be
tested in cosmological observations. For instance, the slow-roll single field
inflation naturally leads to a nearly scale-invariant power spectrum, which
matches very well with the recent observations from the cosmic microwave
background (CMB) measurements [4]. During the inflationary phase, the Universe
expands exponentially, and therefore, inflation ends with an empty (no
radiation/matter) and cold (non-thermal) Universe with total energy density
stored in the inflaton field. To populate the Universe, one needs a mechanism
that converts the inflaton energy density to the Standard Model (SM) and
possibly to the dark sector. The process of transferring energy density from
the inflaton field to the SM through perturbative decays is referred to as
reheating [5, 6, 7, 8, 9, 10].

The perturbative reheating can be realized through some interactions between
the inflaton and the SM fields. The lowest dimensional SM gauge singlet
operator is $|{\bm{h}}|^{2}\\!=\\!{\bm{h}}^{\dagger}{\bm{h}}$, where
${\bm{h}}$ is the SM Higgs doublet. Hence, in generic scalar field
inflationary models, one would expect the leading inflaton–SM interaction is
through the $\phi|{\bm{h}}|^{2}$ term. Whereas interactions of the scalar
inflaton with the SM gauge singlet operators involving fermions and gauge
bosons are higher dimensional $(\geq 5)$ and, therefore, would be suppressed
by the inflationary scale $\Lambda$. In order not to spoil the flatness of the
inflaton potential, such interactions are expected to be subdominant during
the inflationary phase. However, after the end of inflation, the perturbative
reheating predominantly follows through the inflaton-Higgs interaction 111We
note that non-perturbative effects can potentially become relevant for larger
inflaton-Higgs couplings, which lead to tachyonic resonant production of Higgs
modes, i.e., the so-called preheating regime. However, the tachyonic resonant
Higgs production is suppressed in the presence of relatively large quartic SM
Higgs coupling, see Sec. IV.5. Therefore, the dominant mechanism of energy
transfer from the inflaton field to the Higgs field remains the perturbative
decay, see also [11, 12].. Regarding the inflaton field as a classical
background field, the $\phi|{\bm{h}}|^{2}$ term is also a source of
$\phi$-dependent Higgs mass that oscillates in time due to coherent
oscillations of the inflaton field.

In this work, we aim to study in detail consequences of Higgs dynamics due to
inflaton–Higgs interaction during the reheating phase. In particular, our goal
is to analyze the implications of the $\phi$-dependent Higgs mass, which
oscillates and results in rapid transitions between phases of broken and
unbroken electroweak symmetry. We notice that this non-trivial Higgs mass
leads to the suppression of perturbative decays of the inflaton field to Higgs
boson pairs, which not only leads to elongations of the reheating period but
also suppresses the production of the SM radiation energy density. As a
result, evolution of the temperature of the SM bath is modified, which in turn
can significantly affect the freeze-in production of dark matter (DM) during
the reheating phase, see also [13, 14, 15, 16]. The reheating dynamics due to
this non-trivial $\phi$-dependent Higgs mass is referred to as the massive
reheating scenario, whereas for a comparison we consider the case where such
mass effects are neglected, hence referred to as the massless reheating
scenario. In this work we consider reheating through the SM Higgs boson,
however, the results obtained here are straightforward to generalize to any
other scalar field which interacts with the inflation field.

As an example, we employ the $\alpha$-attractor T-model of inflation [17, 18]
whose potential is approximately flat for large inflaton field values suitable
for inflation and it has a monomial shape of the form $\propto\phi^{2n}$
during the reheating phase. In this work, we consider generic $n$ which leads
to an effective equation of state $w\\!=\\!(n-1)/(n+1)$ during the reheating
phase and determines the evolution of the Universe during this period 222This
work is a continuation and significant extension of the research described in
our earlier paper [12] which was limited to the $n\\!=\\!1$ case only and
focused on implications of time-dependent inflaton decay width.. We provide
analytic and numerical results for the dynamics of the inflaton field during
the inflationary and reheating phases in the presence of inflaton–Higgs
interaction.

Furthermore, to investigate non-trivial implications of the Higgs-induced
reheating on the DM production, we consider a model with the DM candidate
being a massive vector field $X_{\mu}$ of a dark Abelian gauge symmetry
$U(1)_{X}$. The vector DM $X_{\mu}$ interacts with the SM as well as with the
inflaton field through gravity and higher dimensional operators suppressed by
the Planck mass. We study production of such DM particles during the reheating
phase. Since the DM interactions with the SM and inflaton are assumed to be
Planck mass suppressed, therefore vector DM is not in thermal equilibrium with
the SM bath. We study DM production through its gravitational interactions
with the inflaton background field [19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 12, 30, 31] as well as through the annihilation of SM particles, also
known as the gravitational freeze-in mechanism [31, 32, 33, 34, 35, 36, 37].
Moreover, DM production through direct inflaton decay is also investigated. In
the case of vector DM, the leading contribution from the inflaton decay turns
out to be triggered by a dim-5 operator suppressed by the Planck mass. Another
source of DM production is through an effective dim-6 operator suppressed by
the Planck mass squared involving the Higgs doublet and DM fields [38, 35, 37,
39]. In this case, the production can occur due to the annihilation and decay
of the Higgs bosons.

The paper is organized as follows. In Sec. II we describe details of our
model, which include the inflaton, the SM Higgs, and the vector DM. In Sec.
III, we study inflaton dynamics during the inflationary phase as well as
during the early stages of reheating, where the SM radiation energy density is
negligible compared to that of the inflaton field. The reheating dynamics is
presented in Sec. IV where we analyze non-trivial effects of the inflaton-
induced Higgs mass on the production of SM radiation quanta, including effects
of rapid electroweak phase transitions. In particular, we present the analytic
and numerical results for the SM radiation energy density and the SM bath
temperature evolution during this phase. In Sec. V we discuss implications of
the Higgs-induced reheating on the gravitational production of DM via graviton
exchange from the inflaton background and SM radiation. Moreover, we study the
DM production due to the inflaton decay and through the direct annihilation
and decays of Higgs bosons. We summarize our findings in Sec. VI. In Appendix
A, we present recent constraints on the inflationary model parameters due to
recent CMB measurements by Planck collaboration [4]. Moreover, we supplement
our results with a detailed derivation of the direct production of the SM and
DM particles in the presence of oscillating inflaton background as well as the
inflaton-induced gravitational production in Appendix B.

##  II The model

In this section, we present our model to describe the inflationary/reheating
dynamics and dark matter production. It has been assumed that the mass scale
for interactions between the SM, inflaton, and/or DM sector is set by the
Planck mass $M_{\rm Pl}$. We consider the following action where the SM,
inflaton, and vector DM interact minimally with gravity,

$S=\int\mathrm{d}^{4}x\sqrt{-g}\bigg{[}\frac{M_{\rm Pl}^{2}}{2}R+{\cal
L}_{\phi}+{\cal L}_{\rm SM}+{\cal L}_{\rm DM}+{\cal L}_{\rm int}\bigg{]},$ (1)

with ${\cal L}_{\phi}$, ${\cal L}_{\rm SM}$ and ${\cal L}_{\rm DM}$ being the
Lagrangian densities for the inflaton, SM, and DM, respectively, whereas,
${\cal L}_{\rm int}$ describes interactions among the inflaton, SM, vector DM,
and graviton. Above $M_{\rm Pl}\\!=\\!2.4\times 10^{18}\;\hbox{GeV}$ denotes
the reduced Planck mass, while $R$ is the Ricci scalar for background metric
$g_{\mu\nu}$ and $g$ denotes its determinant. We consider the background
metric as the FLRW metric, i.e.,

$ds^{2}=dt^{2}-a(t)^{2}d{\bm{x}}^{2},$ (2)

where the ${\bm{x}}$ vector denotes the three spatial coordinates and $a(t)$
is the scale factor.

The Lagrangian density for the inflaton field reads

${\cal
L}_{\rm\phi}=\frac{1}{2}\partial_{\mu}\phi\,\partial^{\mu}\phi-V(\phi),$ (3)

where $V(\phi)$ denotes the inflaton potential. Note that in the above
Lagrangian, and hereafter, the indices are raised and lowered via the
background metric $g_{\mu\nu}$. In this work, we consider the
$\alpha$-attractor T-model of inflation [17, 18], such that the inflaton
potential is

$\displaystyle V(\phi)$
$\displaystyle=\Lambda^{4}\tanh^{2n}\bigg{(}\\!\frac{|\phi|}{M}\\!\bigg{)}\simeq\begin{dcases}\Lambda^{4}&|\phi|\gg\\!M\\\
\Lambda^{4}\Big{|}\frac{\phi}{M}\Big{|}^{2n}&|\phi|\ll\\!M\end{dcases},$ (4)

where $\Lambda$ determines the scale of inflation, whereas $M$ is related to
the reduced Planck mass with the $\alpha$ parameter of the $\alpha$-attractor
T-model as $M\equiv\sqrt{6\alpha}\,M_{\rm Pl}$. The potential $V(\phi)$ has a
minimum at $\phi\\!=\\!0$ for positive values of $n$. The current experimental
data [40, 4] constrain values of the potential parameters. The inflation scale
$\Lambda$ can be limited by the amplitude of the scalar perturbations,
$A_{s}$, and the spectral tilt $n_{s}$ from inflationary observables of the
CMB. Planck 2018 data [40] sets the upper limit $\Lambda\lesssim 1.4\times
10^{16}\;\hbox{GeV}$, see Appendix A. Furthermore, the current upper bound on
the tensor to scalar power spectrum ratio, $r\lesssim 0.032$, limits the value
of the $\alpha$ parameter or $M$ from above, such that, $M\lesssim 10M_{\rm
Pl}$. Hereinafter, without loss of generality, we fix $\alpha\\!=\\!1/6$, such
that $M\\!=\\!M_{\rm Pl}$, and $\Lambda\\!=\\!3\times 10^{-3}M_{\rm Pl}$.

In Fig. 1 we have plotted $V(\phi)/\Lambda^{4}$ as a function of the $\phi$
field for $\alpha\\!=\\!1/6$ with benchmark values of
$n\\!=\\!\\{2/3,1,3/2,2\\}$. For large field values, i.e., $|\phi|\gg\\!M$,
the potential has a plateau $V(\phi)\simeq\Lambda^{4}$. This region is
suitable for cosmic inflation, as the flatness of $V(\phi)$ guarantees that
cosmic inflation lasts long enough. In particular, requiring 50 e-folds of
inflation leads to an initial condition on $\phi$ (indicated by the colored,
empty dots) at the scale factor $a_{\star}$, which gives
$\phi(a_{\star})\\!\simeq\\!\\{3.14,3.35,3.55,3.69\\}\,M_{\rm Pl}$ for
$n\\!=\\!\\{2/3,1,3/2,2\\},$ respectively. When $\phi$ rolls down to smaller
values, i.e., $|\phi|\sim\\!M$, inflation ends, and $\phi$ starts to oscillate
around the bottom of its potential. For smaller field values, $|\phi|\ll\\!M$,
i.e., during the reheating phase, the potential can be well approximated by
the power-law form $\Lambda^{4}\big{(}|\phi|/M\big{)}^{2n}$ (dashed curves).
Note that for $n\\!=\\!1$ and $n\\!=\\!2$ we reproduce the standard quadratic
and quartic inflaton potential, respectively.

Figure 1: The $V(\phi)/\Lambda^{4}$ term as a function of $\phi/M_{\rm Pl}$
for different values of $n$ and $\alpha\\!=\\!1/6$. Solid curves present the
full form of $V(\phi)$, while dashed curves show the small-field expansion of
$V(\phi)$ (the lower line of Eq. (4)). The empty, colored dotes indicate the
initial value of the $\phi$ field, for which one gets 50 e-folds of inflation.

The Lagrangian density for the Abelian vector DM $X_{\mu}$ is

$\displaystyle{\cal L}_{\textsc{\rm DM}}$
$\displaystyle=-\frac{1}{4}X_{\mu\nu}X^{\mu\nu}+\frac{1}{2}m_{X}^{2}X_{\mu}X^{\mu},$
(5)

where $X_{\mu\nu}\equiv\partial_{\mu}X_{\nu}-\partial_{\nu}X_{\mu}$ denotes
the field strength tensor. The mass for the DM vector boson, $m_{X}$, is
generated via an Abelian Higgs mechanism with a large expectation value of a
dark Higgs field $\Phi$ so that the radial Higgs mode is heavy and therefore
is integrated out.

We consider the following form of interaction Lagrangian,

$\displaystyle{\cal L}_{\text{int}}$
$\displaystyle=-\bigg{\\{}\frac{h^{\mu\nu}}{M_{\rm
Pl}}\Big{[}T_{\mu\nu}^{\phi}+T_{\mu\nu}^{\textsc{\rm
SM}}+T_{\mu\nu}^{\textsc{\rm DM}}\Big{]}+g_{h\phi}M_{\rm
Pl}\,\phi|{\bm{h}}|^{2}$ $\displaystyle\qquad+\frac{{\cal
C}_{X}^{\phi}\,m_{X}^{2}}{2M_{\rm Pl}}\phi X_{\mu}X^{\mu}+\frac{{\cal
C}_{X}^{\bm{h}}\,m_{X}^{2}}{2M_{\rm
Pl}^{2}}X_{\mu}X^{\mu}|{\bm{h}}|^{2}\bigg{\\}},$ (6)

where $T_{\mu\nu}^{\phi}$, $T_{\mu\nu}^{\textsc{\rm SM}}$, and
$T_{\mu\nu}^{\textsc{\rm DM}}$ are the energy-momentum tensors for the
inflaton, SM, and DM sectors, respectively. Above $h_{\mu\nu}$ denotes the
graviton field and ${\bm{h}}$ is the SM Higgs doublet, which in the linear
parametrization can be written as

${\bm{h}}=\frac{1}{\sqrt{2}}\begin{pmatrix}h_{2}+ih_{3}\\\
h_{0}+ih_{1}\end{pmatrix},$ (7)

where $h_{i}$ $(i\\!=\\!0,1,2,3)$ are the four real scalar components. The
dimensionless constant $g_{h\phi}$ parametrizes interactions between the
inflaton and the SM Higgs, while ${\cal C}_{X}^{\phi}$ and ${\cal
C}_{X}^{\bm{h}}$ are the dimensionless Wilson coefficients for the DM–inflaton
and DM–Higgs effective interactions, respectively. In what follows, we assume
that all couplings i.e., $g_{h\phi}$, $\mathcal{C}_{X}^{\phi}$ and
$\mathcal{C}_{X}^{{\bm{h}}}$ are real and positive. Note that higher powers of
$\phi$ could also appear in the interaction Lagrangian (6); however, for
simplicity, we will consider only the lowest-dimensional operators that allow
the inflaton to communicate with the SM and DM sectors. In the case of the
vector DM discussed here, the lowest order direct interaction between the SM
and DM appears through the dim-6 operator ${\cal
C}_{X}^{\bm{h}}m_{X}^{2}/(2M_{\rm Pl}^{2})X_{\mu}X^{\mu}|{\bm{h}}|^{2}$ [38]
and therefore must be suppressed by the UV cut-off scale assumed here to be
$M_{\rm Pl}$. Moreover, the interaction Lagrangian contains only linear terms
in the graviton field $h_{\mu\nu}$ which induces tree-level interactions
between the SM and DM via a single graviton s-channel exchange [32, 34, 33].

Note that the last two terms in Eq.(6) can be written in a proper gauge
invariant form by coupling the inflaton field and the SM Higgs doublet to the
kinetic term of the dark Higgs field $\Phi$, i.e.,

$\displaystyle{\cal L}_{\text{int}}\supset-\frac{1}{2}\frac{{\cal
C}_{X}^{\phi}}{M_{\rm Pl}}\phi\,|D_{\mu}\Phi|^{2}-\frac{1}{2}\frac{{\cal
C}_{X}^{\bm{h}}}{M_{\rm Pl}^{2}}|{\bm{h}}|^{2}|D_{\mu}\Phi|^{2},$ (8)

where $D_{\mu}\Phi=(\partial_{\mu}-ig_{X}X_{\mu})\Phi$ is the covariant
derivative of the $\Phi$ field and $g_{X}$ is the $U(1)_{X}$ gauge coupling.
When the dark Higgs field $\Phi$ acquires vacuum-expectation-value (vev)
$\langle\Phi\rangle$, the dark $U(1)_{X}$ symmetry is spontaneously broken,
and a mass term for the vector DM field is generated with mass $m_{X}\equiv
g_{X}\langle\Phi\rangle$. Assuming the radial mode of the dark Higgs is much
heavier than the dark gauge boson, we can integrate out so that effective
interaction terms for the vector DM with the inflaton and the SM Higgs doublet
are generated as in the interaction Lagrangian (6). Furthermore, one can also
write a dim-5 effective interaction between the inflaton and vector DM of the
form, $(\phi/M_{\rm Pl})X_{\mu\nu}X^{\mu\nu}$. However, note that the
corresponding vertex would involve the momentum of the gauge bosons, which is
of the same order as the DM mass since the DM is not in thermal equilibrium.
Therefore, effectively the above operator is similar to the one considered in
Eq. (6).

It is also important to emphasize that the graviton coupling to the SM and DM
energy-momentum tensors (6) leads to an indirect interaction between these two
sectors proportional to $1/M_{\rm Pl}^{2}$, which is of the same order as the
interaction via the effective DM–SM Higgs portal operator for ${\cal
C}_{X}^{\bm{h}}\\!\sim\\!{\cal O}(1)$. This was one of our primary motivations
to consider such effective operators in the model. On the other hand, spin-0
and spin-1/2 DM particles would interact with the SM sector through effective
operators of dim-4 and dim-5, respectively. Then, DM would couple to the SM
sector much stronger than its coupling induced by a single graviton exchange.
In such a scenario, one needs tuning of the DM-SM effective couplings to be
suppressed to get a similar strength to the gravitational interaction.
Therefore, we have found these options less interesting than the spin-1 case.

Note also that a possible dim-4 mixing of the SM hypercharge and Abelian dark
gauge bosons, of the form $\epsilon B_{\mu\nu}X^{\mu\nu}$, is absent in our
model due to dark charge conjugation symmetry, which ensures the stability of
DM. However, for mixing parameter $\epsilon$ small enough, such mixing would
be allowed if the DM lifetime was longer than the age of the Universe.

It is crucial to notice that dynamics of the $\phi$ field during the
inflationary epoch could be modified by its couplings to the SM Higgs and DM
specified in ${\cal L}_{\rm int}$ (6). Therefore, let us first discuss
perturbativity limits imposed upon the strength of the $g_{h\phi}$ coupling
and the Willson coefficients ${\cal C}_{X}^{\phi}$ and ${\cal
C}_{X}^{\bm{h}}$. Requiring, for instance, that amplitude for ${\bm{h}}$
scattering in a background of an inflaton classical field with its double
insertion is smaller than the corresponding amplitude with a single insertion
implies

$g_{h\phi}\lesssim\left(\frac{\Lambda^{2}}{\phi M_{\rm Pl}}\right),$ (9)

where $\phi$ denotes the strength of the external classical inflaton field
determined by its equation of motion (EOM). Similar reasoning implies

${\cal C}_{X}^{\phi}\lesssim\left(\frac{M_{\rm
Pl}}{\phi}\right)\left(\frac{\Lambda}{m_{X}}\right)^{2},$ (10)

where it has been assumed that $\Lambda\gtrsim m_{X}$. On the other hand,
perturbativity of the $X_{\mu}X^{\mu}|{\bm{h}}|^{2}$ operator together with
the requirement that the maximum temperature $T_{\rm max}$ of the thermal bath
during reheating remains below the cut-off scale, $\Lambda$, for the effective
operator $(D_{\mu}\Phi)^{\dagger}(D^{\mu}\Phi)|{\bm{h}}|^{2}$ implies

${\cal C}_{X}^{\bm{h}}\lesssim\min\left\\{\left(\frac{M_{\rm
Pl}}{m_{X}}\right)^{2},\left(\frac{M_{\rm Pl}}{T_{\rm
max}}\right)^{2}\right\\}.$ (11)

Before moving forward, let us discuss here the role of the SM Higgs field
during the inflation and reheating phases in more detail. Note that typical
energy scales for inflation and reheating are much larger than the scale of
electroweak interactions. Therefore, the Higgs potential during inflation and
reheating periods can be approximated as

$\displaystyle V({\bm{h}})$
$\displaystyle=\lambda_{h}\bigg{(}|{\bm{h}}|^{2}-\frac{v_{\rm
EW}^{2}}{2}\bigg{)}^{2}\simeq\lambda_{h}|{\bm{h}}|^{4},$ (12)

where $\lambda_{h}$ is the Higgs-boson quartic coupling. At the electroweak
scale $v_{\rm EW}\simeq 246\;\hbox{GeV}$, the value of Higgs quartic coupling
is $\lambda_{h}\simeq 0.13$ corresponding to the Higgs mass
$m_{h_{0}}^{\tiny{{\rm EW}}}=125\;\hbox{GeV}$. Due to quantum corrections,
within the SM, the Higgs quartic coupling $\lambda_{h}$ runs down to negative
values at energy scales $\sim 10^{10}\;\hbox{GeV}$ making the Higgs potential
unstable [41]. However, in the presence of the new physics interactions as
considered in this work in Eq. (6), the Higgs stability can be achieved for
larger Higgs field values [41]. In particular, the inflaton–Higgs interaction
term, $g_{h\phi}M_{\rm Pl}\,\phi\,|{\bm{h}}|^{2}$, generates a
$\phi$-dependent Higgs-boson mass $m_{h_{0}}^{2}\\!=\\!g_{h\phi}M_{\rm
Pl}|\phi|$, with $|\phi|\\!\sim\\!{\cal O}({\rm few})M$ during inflation,
which leads to the stability/positivity of the Higgs potential. The maximal
Higgs field strength up to which the potential is stable could be estimated as
$|{\bm{h}}_{\rm max}|^{2}\sim g_{h\phi}M_{\rm Pl}\phi/\lambda_{h}$.

During the inflationary phase, besides the inflaton potential $V(\phi)$, one
could consider contributions to inflaton dynamics that originate from the SM
Higgs quartic interactions and the inflaton–Higgs coupling. However, in this
work, we consider a parameter space where the inflationary dynamics is
dominated by the cosmological constant term $\Lambda^{4}$. Thus, we require
that both $\lambda_{h}|{\bm{h}}|^{4}$ and $g_{h\phi}M_{\rm
Pl}\,\phi\,|{\bm{h}}|^{2}$ terms are smaller than $\Lambda^{4}$ in the whole
region of stability, i.e., up to the largest allowed Higgs-field strength
${\bm{h}}={\bm{h}}_{\rm max}$. It is easy to see that the required condition
reads

$g_{h\phi}\lesssim\sqrt{\lambda_{h}}\left(\frac{\Lambda^{2}}{\phi M_{\rm
Pl}}\right).$ (13)

Which gives a stronger bound than Eq. (9) for $\lambda_{h}<1$, and hence we
adopt the (13) constraint in the following analysis.

Furthermore, as shown in Ref. [41], for the Higgs mass
$m_{h_{0}}\\!>\\!3H_{I}/2$ the Higgs field fluctuations during inflation are
strongly suppressed ensuring stability, where $H_{I}$ denotes the Hubble
parameter during inflation. This condition could be written in terms of a
lower limit for $g_{h\phi}$, i.e.,

$g_{h\phi}\gtrsim\frac{3}{4}\sqrt{6\alpha}\left(\frac{\Lambda^{2}}{\phi M_{\rm
Pl}}\right)^{2}\left(\frac{\phi}{M}\right).$ (14)

For the benchmark values of $\Lambda\\!=\\!3\times 10^{-3}M_{\rm Pl}$ and
$\alpha\\!=\\!1/6$, and assuming that $\phi\\!\sim\\!M\\!=\\!M_{\rm Pl}$, we
obtain the following consistency region for $g_{h\phi}$:

$6\cdot 10^{-11}\lesssim g_{h\phi}\lesssim 3\cdot 10^{-6}$ (15)

Hereinafter, for numerical calculations, we adopt the EW value of the Higgs
quartic coupling $\lambda_{h}=0.13$. Note that the above limits are a subject
of $\lambda_{h}$; increasing $\lambda_{h}$ implies larger $g_{h\phi}$ allowed
values.

##  III Inflaton dynamics

In this section, we discuss the dynamics of the inflaton field. Based on the
arguments above, we ignore Higgs boson contributions to the inflaton dynamics
so that the classical equation of motion for $\phi$ in the FLRW background is
given by

$\ddot{\phi}+3H\dot{\phi}+V_{,\phi}(\phi)=0,\,$ (16)

where $H\equiv\dot{a}/a$ is the Hubble rate, the $\it{overdot}$ denotes
derivative w.r.t. the cosmic time, $t$, whereas $V_{,\phi}(\phi)$ is
derivative of the potential w.r.t. the $\phi$ field. In the above EoM, we have
neglected spatial derivatives of the inflaton, regarding $\phi$ as a
spatially-homogeneous scalar field. Later on, for collision terms in the
Boltzmann equations, we will calculate $S$-matrix elements in the presence of
classical inflaton fields that are solutions of (16). In other words, we will
treat the SM and DM fields as small perturbations as compared to solutions of
(16) and their possible back-reaction on the inflaton field will be neglected
333This is similar, e.g., to the Rutherford scattering which is an elastic
scattering of a charged particle by a static Coulomb point-particle potential
that is a solution of Maxwell equations in empty space undisturbed by the
presence of any other particles.. As discussed above, the inflaton-Higgs
interaction term, i.e., $g_{h\phi}M_{\rm Pl}\phi|\bm{h}|^{2}$, is suppressed
during inflation; therefore, Eq. (16) is valid during this period as well as
during the early stages of reheating. However, once the reheating resumes, the
Higgs field is produced through the coherent oscillations of the inflaton
field. As we will see below, the Higgs dynamics will be taken into account in
the Boltzmann equation. Assuming that the SM radiation and DM energy densities
are negligible in the primordial Universe, one can write the first Friedmann
equation as

$H^{2}\simeq\frac{\rho_{\phi}}{3M_{\mathrm{Pl}}^{2}}.\,$ (17)

The energy density, $\rho_{\phi}$, and the pressure, $p_{\phi}$, for the
homogenous, scalar field are

$\displaystyle\rho_{\phi}$ $\displaystyle=\frac{1}{2}\dot{\phi}^{2}+V(\phi),$
$\displaystyle p_{\phi}$ $\displaystyle=\frac{1}{2}\dot{\phi}^{2}-V(\phi).$
(18)

We can now discuss solutions to Eq. (16) during and after inflation. As we
have already pointed out above, during the phase of the accelerated expansion,
$V(\phi)$ is approximately flat. In this period $\ddot{\phi}$ is negligible
and $V(\phi)\\!\gg\\!\dot{\phi}^{2}/2$, which in turn implies $H^{2}\simeq
V(\phi)/(3M_{\rm Pl}^{2})$. Using these two simplifying assumptions during the
inflationary phase one finds constant solutions for the inflaton field
$\phi_{I}(t)\simeq\phi(a_{\star})$ and the Hubble rate
$H_{I}\simeq\Lambda^{2}/(\sqrt{3}M_{\rm Pl})$. The slow-roll evolution of the
inflaton field can be parameterized by the so-called potential slow-roll
parameters $\epsilon_{V}$ and $\eta_{V}$, defined as

$\displaystyle\epsilon_{V}$
$\displaystyle=\frac{M_{\mathrm{Pl}}^{2}}{2}\left(\frac{V_{,\phi}(\phi)}{V(\phi)}\right)^{2},$
$\displaystyle\eta_{V}$
$\displaystyle=M_{\mathrm{Pl}}^{2}\frac{V_{,\phi\phi}(\phi)}{V(\phi)}.$ (19)

Note that the two cosmological observables, the spectral index $n_{s}$ and the
tensor-to-scalar ratio $r$, can be related to the potential slow-roll
parameters as

$\displaystyle r$ $\displaystyle=16\epsilon_{V},$ $\displaystyle n_{s}$
$\displaystyle=1-6\epsilon_{V}+2\eta_{V}.$ (20)

For the $\alpha$-attractor inflaton potential (4), the slow-roll potential
parameters are

$\displaystyle\epsilon_{V}$ $\displaystyle=8n^{2}\left(\frac{M_{\rm
Pl}}{M}\right)^{2}\operatorname{csch}^{2}\left(\frac{2|\phi|}{M}\right),$ (21)
$\displaystyle\eta_{V}$ $\displaystyle=8n\left(\frac{M_{\rm
Pl}}{M}\right)^{2}\left[2n-\cosh\left(\frac{2|\phi|}{M}\right)\right]\operatorname{csch}^{2}\left(\frac{2|\phi|}{M}\right).$
(22)

To achieve sufficiently many e-folds ($\sim 50-60$) of expansion, the slow-
roll parameters are required to be small, i.e., $\epsilon_{V},|\eta_{V}|\ll
1$, for a long enough period. The accelerated expansion of the Universe ends
when $\ddot{a}=0$, which corresponds to $\epsilon_{V}\simeq 1$. This condition
determines the value of $\phi(a_{e})$ at the end of inflation, i.e., at
$a\\!=\\!a_{e}$.

After the end of inflation, the inflaton field $\phi$ starts to oscillate
around the minimum of its potential with decreasing amplitude. Moreover, the
character of the oscillations strongly depends on the shape of $V(\phi)$ in
the vicinity of the minimum. Note that the time scale of the inflaton
oscillations is typically much shorter than the variation of its amplitude.
Thus, in this phase, a generic solution to Eq. (16) can be written as a
product of two functions [9, 42, 13, 26],

$\displaystyle\phi(t)=\varphi(t)\cdot\mathcal{P}(t),$ (23)

where $\mathcal{P}(t)$ is a quasi-periodic, fast-oscillating function, while
$\varphi(t)$ denotes a slowly-varying (w.r.t. the time scale of the
oscillations) envelope function defined by the following condition:

$\displaystyle\rho_{\phi}=V(\varphi)=\Lambda^{4}\Big{(}\frac{\varphi}{M_{\rm
Pl}}\Big{)}^{2n}\,,$ (24)

where the inflaton potential has been expanded in a form applicable during the
reheating phase with $\alpha\\!=\\!1/6$. In Fig. 2 we show the numerical
solution for $\phi(a)$, obtained after solving two coupled equations (16) and
(17). The rapid oscillations of $\mathcal{P}(a)$ are damped by the decreasing
envelope $\varphi(a)$ due to the redshift of the Universe. The frequency of
oscillations as well as their amplitude depends on the slope of the inflaton
potential $n$. The explicit form of $\varphi$ and $\mathcal{P}$ will be
determined below.

Figure 2: Numerical solutions for the inflaton EoM (16) and (17) as a function
of the scale factor $a$ for various choices on $n$. The vertical dashed line
indicates the end of inflation, defined by the condition
$\epsilon_{V}(\phi(a_{e}))=1$. Colored dots present the value of $\phi(a_{e})$
for different values of $n$.

Our goal now is to obtain approximate analytical solutions of (16) during the
oscillatory phase. Before we do that, let us first find a relation between the
inflaton energy density $\rho_{\phi}$ and pressure $p_{\phi}$. To that end, we
differentiate $\rho_{\phi}$ (18) w.r.t. time and using (16) we obtain the
continuity equation

$\displaystyle\dot{\rho}_{\phi}+3H(1+w)\rho_{\phi}$ $\displaystyle=0,$
$\displaystyle w$ $\displaystyle\equiv p_{\phi}/\rho_{\phi},$ (25)

where the barotropic parameter $w$ is, in general, time-dependent. Note that
at this stage, the effects of quantum particle production are ignored, so that
(25) holds during inflation and the beginning of reheating as long as the
inflaton energy density dominates, i.e., the inflaton decay rate is smaller
than the Hubble rate. Ignoring expansion, assuming periodicity, and averaging
over one period of oscillations, we can express the inflaton energy density
and pressure (18) as [9, 42, 13, 26],

$\displaystyle\langle\rho_{\phi}\rangle$
$\displaystyle=V(\varphi)=\rho_{\phi}$ (26) $\displaystyle\langle
p_{\phi}\rangle$
$\displaystyle=\frac{n-1}{n+1}V(\varphi)=\frac{n-1}{n+1}\rho_{\phi},$ (27)

where we have used Eq. (24) along with the relation
$\langle\lvert\mathcal{P}^{2n}\rvert\rangle\\!=\\!1/(n+1)$. Hereinafter, the
following definition for the time-average of a quantity $f(t)$ over one period
of inflaton oscillation is employed,

$\displaystyle\langle f(t_{i})\rangle=\frac{1}{{\cal
T}(t_{i})}\int_{t_{i}}^{t_{i}+{\cal T}(t_{i})}\\!\\!\\!dt\;f(t),$ (28)

where $t_{i}$ is some reference time at a particular instant during the
oscillatory phase. Above, ${\cal T}(t_{i})$ denotes the period of the inflaton
oscillations, which in general can be time-dependent. We can now define the
averaged equation-of-state parameter during the reheating phase as

$\displaystyle\bar{w}\equiv\frac{\langle
p_{\phi}\rangle}{\langle\rho_{\phi}\rangle}=\frac{n-1}{n+1}.$ (29)

In particular, we have $\bar{w}\\!\in\\!\\{-1/5,0,1/5,1/3\\}$ for
$n\\!\in\\!\\{2/3,1,3/2,2\\},$ respectively. The evolution of the equation of
state $w\equiv\rho_{\phi}/p_{\phi}$ and the its time-averaged value $\bar{w}$
is shown in Fig. 3 as a function of the scale factor $a$. During inflation,
i.e., for $a\\!<\\!a_{e}$, the potential term $V(\phi)$ dominates over the
kinetic term, $\dot{\phi}^{2}/2$, and the $\rho_{\phi}/p_{\phi}$ ratio is
constant, i.e., $w\simeq-1$. After the end of inflation, $w$ starts to
oscillate between $-1$ and $1$, while $\bar{w}$ quickly approaches a constant
limit, consistent with the prediction (29).

Figure 3: The equation-of-state parameter, $w$, as a function of the scale
factor $a$. Thin colored curves show the full numerical solutions, i.e.,
$p_{\phi}/\rho_{\phi}$, while thick lines present the time-averaged barotropic
parameter, i.e., $\bar{w}\\!=\\!\langle
p_{\phi}\rangle/\langle\rho_{\phi}\rangle$ (29). The horizontal axes are
normalized to the scale factor at the end of inflation $a_{e}$.

Using the continuity equation (25) together with the definition (24) and (29)
one can derive the equation of motion for the envelope

$\displaystyle\dot{\varphi}(t)=-\frac{3}{n+1}H\varphi(t).$ (30)

Form of the above equation proves that indeed $\varphi$ is a slowly varying
function of time and fast oscillations are relegated to the $\mathcal{P}$
function 444Note that the definition of the envelope function (24) implies
that $d\langle\rho_{\phi}\rangle/dt=\langle\dot{\rho}_{\phi}\rangle$ assuming
periodicity and neglecting cosmological time evolution while averaging over a
period.. The solution to this equation as a function of the scale factor $a$
reads

$\displaystyle\varphi(a)=\varphi_{e}\,\Big{(}\frac{a_{e}}{a}\Big{)}^{\frac{3}{n+1}},$
(31)

where $\varphi_{e}\equiv\varphi(a_{e})$ denotes the initial value of the
envelope determined by the condition $\epsilon_{V}(\phi(a_{e}))=1$ at the end
of inflation, cf. Eq.(21), and for $\alpha=1/6$ it is given by

$\displaystyle\varphi_{e}\equiv\varphi(a_{e})=\frac{M_{\rm
Pl}}{2}\operatorname{ArcSinh}\big{(}2\sqrt{2}\,n\big{)}.$ (32)

Note that for $n\\!>\\!0$, $\varphi(a)$ decreases with time during the
oscillatory phase. Using the Friedmann equation (17) together with Eq. (24),
we obtain the following equation that governs the evolution of the scale
factor $a$ with $t$,

$\displaystyle\frac{da}{a}=\frac{1}{\sqrt{3}M_{\rm
Pl}}\Lambda^{2}\left(\frac{\varphi}{M_{\rm Pl}}\right)^{n}dt.$ (33)

Using Eq. (31) we get

$\displaystyle
a(t)=a_{e}\bigg{[}1+\frac{3n}{n+1}\frac{\Lambda^{2}}{\sqrt{3}M_{\rm
Pl}}\Big{(}\frac{\varphi_{e}}{M_{\rm
Pl}}\Big{)}^{n}(t-t_{e})\bigg{]}^{\frac{n+1}{3n}},$ (34)

where $t_{e}$ denotes the cosmic time at end of inflation. Note that the above
relation implies $\varphi(t)\\!\sim\\!t^{-1}$ for $n\\!=\\!1$, whereas, for
the quartic potential, i.e., $n\\!=\\!2$, we get
$\varphi(t)\\!\sim\\!t^{-1/2}$.

In order to find a dynamical equation for the quasi-periodic function
$\mathcal{P}(t)$ we differentiate (23) w.r.t. time, obtaining

$\displaystyle\dot{\phi}=\dot{\varphi}\mathcal{P}+\varphi\dot{\mathcal{P}},$
(35)

which, in turn, implies

$\displaystyle\dot{\mathcal{P}}(t)=\pm\frac{\sqrt{2(\rho_{\phi}-V(\phi))}}{\varphi}+\frac{3n}{n+1}H\mathcal{P}(t),$
(36)

where the $\pm$ sign corresponds to $\dot{\phi}\\!>\\!0$ and
$\dot{\phi}\\!<\\!0$, respectively. Note that the second term in the r.h.s. of
(36) is negligible as long as the time scale of the inflaton oscillations is
much shorter than the time scale of expansion. Consequently, one can simplify
Eq.(36) as

$\displaystyle\dot{\mathcal{P}}$
$\displaystyle\simeq\pm\frac{m_{\phi}}{\sqrt{n(2n-1)}}\sqrt{1-\lvert\mathcal{P}\rvert^{2n}},$
(37)

where we have introduced an effective mass of the inflaton field, see also
[43, 13, 26],

$\displaystyle
m_{\phi}^{2}\\!\equiv\\!\frac{\partial^{2}V(\phi)}{\partial\phi^{2}}\bigg{|}_{\phi=\varphi}$
$\displaystyle\\!=\\!2n(2n-1)\frac{\Lambda^{4}}{M_{\rm
Pl}^{2}}\Big{(}\frac{\rho_{\phi}}{\Lambda^{4}}\Big{)}^{\\!\\!\frac{n-1}{n}},$
(38)

which is time-dependent for $n\\!\neq\\!1$. In particular, since $\rho_{\phi}$
decreases due to the expansion, $m_{\phi}$ increases for $n\\!<\\!1$, while
for $n\\!>\\!1$ it decreases with time in the oscillatory phase. Using Eqs.
(24), (31) and (38) we can obtain an explicit formula for $m_{\phi}$ as a
function of the scale factor $a$

$\displaystyle m_{\phi}(a)=\sqrt{2n(2n-1)}\frac{\Lambda^{2}}{M_{\rm
Pl}}\Big{(}\frac{\varphi_{e}}{M_{\rm
Pl}}\Big{)}^{n-1}\Big{(}\frac{a_{e}}{a}\Big{)}^{\frac{3(n-1)}{n+1}}\,.$ (39)

Since $m_{\phi}$ varies on the time scale much larger than the oscillation
time scale, we can solve Eq. (37) assuming that $m_{\phi}$ is constant during
one period of oscillations. Then the generic solution for the quasi-periodic
function ${\cal P}(t)$ can be written in terms of the inverse of the
regularized incomplete beta function ${\cal I}^{-1}_{z}(i,j)$ as

$\displaystyle\mathcal{P}(a)=\bigg{[}{\cal
I}^{-1}_{z}\\!\Big{(}\frac{1}{2n},\frac{1}{2}\Big{)}\bigg{]}^{\frac{1}{2n}},$
(40)

with

$\displaystyle z$ $\displaystyle\equiv 1-\frac{4}{\cal T}(t-t_{0})$
$\displaystyle\\!=\\!1-\frac{4}{\mathcal{T}}\frac{n+1}{\sqrt{3}n}\frac{M_{\rm
Pl}}{\Lambda^{2}}\Big{(}\frac{M_{\rm
Pl}}{\varphi_{e}}\Big{)}^{\\!n}\bigg{[}\Big{(}\frac{a}{a_{e}}\Big{)}^{\frac{3n}{n+1}}-1\bigg{]},$
(41)

where the period of the oscillations, $\mathcal{T}$, is given by

$\mathcal{T}\equiv\frac{2\pi}{\omega}=\frac{\sqrt{4\pi}}{m_{\phi}}\sqrt{\frac{2n-1}{n}}\frac{\Gamma\left(\frac{1}{2n}\right)}{\Gamma\left(\frac{n+1}{2n}\right)}\,,$
(42)

with frequency $\omega$. For the quadratic potential i.e., $n\\!=\\!1$ we
recover the well-known results: $\mathcal{T}\\!=\\!2\pi/m_{\phi}$ and
$\mathcal{P}(t)\\!=\\!\cos(m_{\phi}t)$ [10]. Note that the slow variation of
$m_{\phi}$ w.r.t. time is also a source of the variation of the period.
Because $m_{\phi}$ increases (decreases) for $n\\!<\\!1$ ($n\\!>\\!1$),
$\mathcal{T}$ decreases (increases) with time. In the next section, we employ
these results obtained for the inflaton profile as well as its time-dependent
mass to study the reheating dynamics, including the inflaton interactions with
the Higgs boson.

##  IV Reheating dynamics

The reheating dynamics involve the coherently oscillating inflaton field
$\phi$, the SM radiation produced through decays of the inflaton field to the
SM Higgs boson, and the vector DM $X_{\mu}$. To track the evolution of this
system, one has to consider three coupled Boltzmann equations averaged over
the inflaton oscillations. The energy or number density of each considered
species changes due to the expansion of the Universe, controlled by terms
proportional to the Hubble rate $H$, and interactions between those three
sectors. In what follows, we assume that $\phi$ can directly “decay” 555The
inflaton “decay” to the SM Higgs boson and vector DM can be understood as the
production of the pairs of SM Higgs bosons and DM vector bosons in a quantum
process out of the vacuum in the presence of the classical inflaton field
(23). to the SM Higgs pairs, with an averaged (over each oscillation period)
decay rate $\langle\Gamma_{\phi\rightarrow hh}\rangle$, and to DM vectors,
with an averaged rate $\langle\Gamma_{\phi\rightarrow XX}\rangle$. Due to the
fact that the latter process arises from the dim-5 operator,
$\phi(D_{\mu}\Phi)^{\dagger}(D^{\mu}\Phi)/M_{\rm Pl}$, the
$\phi\\!\rightarrow\\!XX$ process is subdominant, and the reheating dynamics
is dominantly governed by perturbative decays of the inflaton field to the SM
Higgs bosons.

During the reheating period, the Higgs potential can be written as

$\displaystyle V({\bm{h}})$
$\displaystyle=\mu_{h}^{2}\,|{\bm{h}}|^{2}+\lambda_{h}\,|{\bm{h}}|^{4}\,,$
(43)

where the Higgs mass parameter is

$\mu_{h}^{2}(\phi)\equiv g_{h\phi}M_{\rm Pl}\,\phi(t).$ (44)

In the above expression we have neglected the Higgs mass contributions
proportional to the electroweak vev $v_{\rm EW}\\!\simeq\\!246\;\hbox{GeV}$.
In the presence of an oscillating inflating field, the Higgs field undergoes
rapid phase transitions acquiring a $\phi(t)$-dependent vev during the
reheating phase, i.e.,

$\displaystyle v_{h}$ $\displaystyle=\begin{dcases}~{}0\,,&\phi(t)>\\!0\,,\\\
\sqrt{\frac{|\mu_{h}^{2}(\phi)|}{\lambda_{h}}}\,,&\phi(t)<\\!0.\end{dcases}$
(45)

Note that $v_{h}$ is a function of the inflaton field $\phi$, therefore at the
end of the reheating phase, this non-trivial inflaton-induced vacuum-
expectation-value (vev) vanishes, and the Higgs field remains in the symmetric
phase until the electroweak phase transition at the temperature scale ${\cal
O}(100)\;\hbox{GeV}$. Furthermore during reheating all SM massive particles
receive a non-zero mass due to their coupling to the Higgs field, i.e.,

$\displaystyle m_{\rm{SM}}=\frac{m_{\rm{SM}}^{\rm{EW}}}{v_{\rm{EW}}}v_{h},$
(46)

where $m_{\rm{SM}}^{\rm{EW}}$ denotes the EW mass of SM particles.

Before moving ahead, we would like to comment on the gravitational production
of the SM radiation from the inflaton background through a graviton exchange
[26, 24, 28, 44]. This production mechanism is universal and independent of
any direct coupling of the inflaton to the SM sector. However, as we show in
Appendix B, in the presence of the inflaton-Higgs interaction (6) with
$g_{h\phi}$ satisfying the bound (227) the gravitational production of SM
radiation through the inflaton background field is subdominant. Therefore, we
do not discuss the gravitational production of the SM radiation in the
following and only consider the production of the SM radiation through
inflaton decays to the SM Higgs.

During the reheating epoch the DM particles can be produced through the
following types of processes:

  * •

Inflaton decays: In this case, the DM $X_{\mu}$ particles are produced from
the inflaton decays through the effective interaction (6).

  * •

Gravitational production from the inflaton field: This is an irreducible DM
production mechanism where the inflaton background field transfers its energy
density to the DM through the graviton exchange. This production mechanism is
a result of graviton universal coupling to energy-momentum tensors of the
inflaton and DM fields (6).

  * •

Gravitational freeze-in: Annihilation of SM particles to the DM through
$s$-channel graviton $h_{\mu\nu}$ exchange can be considered as an irreducible
mechanism of DM production, which does not require any additional coupling
between these two sectors. An amplitude for this process is suppressed by
$1/M_{\rm Pl}^{2}$ as a result of graviton universal coupling to energy-
momentum tensors of the SM and DM sectors (6).

  * •

Higgs portal freeze-in and Higgs decays: In this case, the dark sector is
produced through annihilations of Higgs bosons due to the dim-6 operator (6).
In the massive reheating scenario, Higgs boson acquires the inflaton-induced
vev (45) and therefore the contact operator ${\cal
C}_{X}^{\bm{h}}\,m_{X}^{2}/(2M_{\rm
Pl}^{2})X_{\mu}X^{\mu}\lvert{\bm{h}}\rvert^{2}$, after expanding around the
vev $v_{h}$, generates an additional Higgs-DM interaction responsible for
$h_{0}\to XX$ decays, i.e.,

$\displaystyle\frac{{\cal C}_{X}^{\bm{h}}\,m_{X}^{2}}{2M_{\rm
Pl}^{2}}v_{h}h_{0}X_{\mu}X^{\mu}+\cdots\,.$

Furthermore, the above term opens a new DM production channel, in which $X$
particles are created from the annihilation of SM species with $s$-channel
Higgs exchange, i.e., ${\rm SM\,SM}\rightarrow h_{0}\rightarrow XX$. Here we
should point out that both processes $h_{0}\\!\rightarrow\\!XX$ and ${\rm
SM\,SM}\rightarrow h_{0}\rightarrow XX$ are possible only in the broken phase,
i.e., half-period of the inflaton oscillations when $\phi\\!<\\!0$ and the
Higgs boson vev is non-zero. Note that in the massless reheating scenario DM
particles are produced only through the Higgs annihilations via the contact
diagram $hh\\!\rightarrow XX$.

In Fig. 4 we collect all possible interactions among the inflaton, SM and DM
fields. Where “cross” denotes the inflaton background interaction with the
graviton field.

Figure 4: Feynman diagrams for all possible interactions between the inflaton,
SM, and DM, where “cross” indicate the inflaton background interaction with
the graviton field.

We write the time-averaged Boltzmann equations (BEq) for the inflaton energy
density $\rho_{\phi}$, the SM radiation energy density $\rho_{\rm SM}$, and
the DM number density $n_{X}$ during the reheating phase as

$\displaystyle\dot{\rho}_{\phi}\\!+\\!\frac{6n}{n\\!+\\!1}H\rho_{\phi}$
$\displaystyle=-\langle\Gamma_{\phi}\rangle\rho_{\phi},$ (47)
$\displaystyle\dot{\rho}_{\rm SM}+4H\rho_{\rm SM}$
$\displaystyle\\!=\\!\langle\Gamma_{\phi\to{\rm SM}\,{\rm
SM}}\rangle\rho_{\phi}\\!-\\!2\langle
E_{X}\rangle\mathcal{S}_{\rm{SM}}\\!-\\!\langle
E_{h_{0}}\rangle\mathcal{D}_{h_{0}},$ (48) $\displaystyle\dot{n}_{X}+3Hn_{X}$
$\displaystyle=\mathcal{D}_{\phi}+\mathcal{S}_{\phi}+\mathcal{S}_{\rm{SM}}+\mathcal{D}_{h_{0}},$
(49)

where the Hubble rate is

$\displaystyle H^{2}=\frac{1}{3M_{\rm Pl}^{2}}\left(\rho_{\phi}+\rho_{\rm
SM}+\rho_{X}\right).$ (50)

Above $\langle\cdots\rangle$ stands for thermal or (and) time-averaged
quantity. The total inflaton decay rate, $\langle\Gamma_{\phi}\rangle$,
contains two terms induced by the contact inflaton interactions with the Higgs
field and DM vectors, as well as terms generated by indirect interactions
through gravity. For details see Appendix B. The first two terms in the third
Boltzmann equation, i.e., $\mathcal{D}_{\phi}$ and $\mathcal{S}_{\phi}$,
account for the DM pair production from the vacuum in the presence of the
inflaton background via the effective dim-5 operator (6) and through gravity,
respectively. The explicit form of the these two terms has been derived in
Appendix B where we have used the notation ${\cal D}_{\phi}\equiv{\cal
D}^{(1)}_{\phi\to XX}$, and ${\cal S}_{\phi}\equiv{\cal D}^{(2)}_{\phi\to
XX}$. Next two terms in the r.h.s. of Eq.(49) are defined as

$\displaystyle\mathcal{S}_{\rm{SM}}$
$\displaystyle\\!\equiv\\!\bar{n}_{X}^{2}\,\langle\sigma_{XX\rightarrow{\rm
SM\,SM}}\lvert v\rvert\rangle,$ $\displaystyle\mathcal{D}_{h_{0}}$
$\displaystyle\\!\equiv\\!\bar{n}_{h_{0}}\langle\Gamma_{h_{0}\to XX}\rangle,$
(51)

with $\sigma_{XX\rightarrow{\rm SM\,SM}}$ being the DM annihilation cross-
section to the SM, whereas $\Gamma_{h_{0}\to XX}$ denotes the Higgs decay rate
to DM particles. Above $\langle E_{X}\rangle$ and $\langle E_{h_{0}}\rangle$
denote averaged energies of the $X$ and $h_{0}$ particles, respectively. Note
that in the above relations (51) we have assumed that the number density of
the Higgs field follows the thermal equilibrium value, i.e.,
$n_{h_{0}}\\!\simeq\\!\bar{n}_{h_{0}}$, while the dark sector is assumed to be
out-of-thermal equilibrium, i.e., $n_{X}\\!\ll\\!\bar{n}_{X}$. Ignoring
quantum statistics and adopting the Maxwell-Boltzmann distribution function,
the equilibrium number density of species $i$ with spin $J_{i}$ is given by

$\displaystyle\bar{n}_{i}$
$\displaystyle=\frac{(2J_{i}+1)}{2\pi^{2}}m_{i}^{2}\,T\,K_{2}\\!\left(\frac{m_{i}}{T}\right),$
(52)

where $K_{2}(x)$ denotes the modified Bessel function of the second kind.
Finally, let us emphasize that in the above Boltzmann equations we have
assumed that during one oscillation period variations of the Hubble rate $H$
and the effective mass of the inflaton field $m_{\phi}$ could be neglected.

In this section, we assume that the presence of DM interactions does not
modify substantially the reheating dynamics, which is mainly driven by the
first two Boltzmann equations (47)-(48). In particular, we assume that
inflaton predominantly transfers its energy into the radiation sector through
the Higgs portal, since the inflaton’s higher dimensional interactions through
gravity to the SM and DM as well as direct coupling to the vector DM are
negligible as compared to a relevant (dim-3) Higgs portal operator. Such that
$\langle\Gamma_{\phi}\rangle\simeq\langle\Gamma_{\phi\rightarrow{\rm{SM}}\,{\rm{SM}}}\rangle\simeq\langle\Gamma_{\phi\rightarrow
hh}\rangle$. Furthermore, in the evolution of inflation and SM radiation
energy densities, we ignore the DM scatting and Higgs decays to DM particles.
This is a reasonable approximation for the freeze-in DM scenario with a
negligible DM abundance at the onset of the reheating period.

There are several other comments here in order. Firstly, notice that we are
assuming that the whole SM is created from Higgs boson pairs emerging in the
process of inflaton “decays”, i.e., the Higgs bosons decay and scatter to
produce the rest of the SM. Secondly, we assume the SM particles produced are
thermalized instantaneously. This is a reasonable assumption given that once
the Higgs bosons are produced, they would immediately decay/scatter to the
rest of the SM particles, leading to instantaneous thermalization [45].
Thirdly, note that the r.h.s. of (47) was not present when we have discussed
the inflaton dynamics in Sec. III, cf. (25). Here it describes a quantum
process of the particle pair production (in other words, a transition from the
quantum vacuum to the two-particles final state) in the presence of the
oscillating classical inflaton field. Fourthly, let us define the temperature
of the thermal bath, $T$, through the following well-known relation,

$\displaystyle\rho_{\rm SM}=\frac{\pi^{2}}{30}g_{\star}(T)\,T^{4},$ (53)

where $g_{\star}(T)$ counts an effective number of relativistic degrees of
freedom at the temperature $T$. Here we approximate $g_{\star}(T)$ by the
constant value of $106.75$.

Finally, we would like to comment on possible thermal effects on the Higgs
potential (43) and in particular thermal corrections to the Higgs mass
parameter $\mu_{h}^{2}$ during the reheating phase. In the large temperature
limit, i.e., $T\\!\gg\\!m_{\rm{SM}}$, the thermal corrections to the Higgs
mass $\delta\mu_{h}^{2}(T)$ are [46],

$\displaystyle\delta\mu_{h}^{2}(T)$
$\displaystyle\simeq\bigg{[}\frac{3g^{2}}{16}+\frac{g^{\prime
2}}{16}+\frac{y_{t}^{2}}{4}+\frac{\lambda_{h}}{2}\bigg{]}\,T^{2}\simeq
0.4\,T^{2},$ (54)

where $g/g^{\prime}$ and $y_{t}$ are the SM electroweak gauge and top Yukawa
couplings, respectively. The temperature-dependent contribution to the Higgs
mass arises from the Higgs interactions with the high-temperature bath of
relativistic particles. In the following analysis, we neglect such thermal
corrections which become relevant only when
$T\\!\gtrsim\\!|\mu_{h}(\varphi)|$.

###  IV.1 Evolution of the inflaton and SM energy densities

In the following, our goal is to find approximate analytical solutions for the
inflaton and the SM radiation energy densities by solving the first two
Boltzmann Eqs. (47)-(48). After employing the simplifications discussed above,
the Boltzmann Eqs. (47)-(48) take the simple form,

$\displaystyle\dot{\rho}_{\phi}\\!+\\!\frac{6n}{n\\!+\\!1}H\rho_{\phi}$
$\displaystyle=0,$ (55) $\displaystyle\dot{\rho}_{\rm SM}+4H\rho_{\rm SM}$
$\displaystyle=\langle\Gamma_{\phi\to hh}\rangle\rho_{\phi}.$ (56)

Above in the first equation, we have neglected the inflaton decay rate
$\langle\Gamma_{\phi}\rangle$ compared to the Hubble rate $H$ during the
reheating phase. We use the above set of simplified Boltzmann equations to
analytically determine the reheating dynamics, however, note that for the
numerical analysis we consider solutions of the exact Boltzmann Eqs.
(47)-(48). The initial conditions for the inflaton energy density
$\rho_{\phi}$ and the SM energy density $\rho_{\rm SM}$ at the onset of the
reheating phase are

$\displaystyle\rho_{\phi_{e}}$
$\displaystyle\equiv\rho_{\phi}(a_{e})=\Lambda^{4}\Big{(}\frac{\varphi_{e}}{M_{\rm
Pl}}\Big{)}^{2n},$ $\displaystyle\rho_{\rm SM}(a_{e})$ $\displaystyle=0,$ (57)

where we have employed Eq. (24). Hereinafter, we use the scale factor $a$ as
an independent time variable rather than physical time $t$, where the two are
related during the reheating phase through Eq. (34). Therefore, the inflaton
energy density during reheating evolves as

$\displaystyle\rho_{\phi}(a)$
$\displaystyle=\rho_{\phi_{e}}\Big{(}\frac{a_{e}}{a}\Big{)}^{\frac{6n}{n+1}},$
for $\displaystyle a\in[a_{e},a_{\tiny\text{rh}}],$ (58)

where $a_{\tiny\text{rh}}$ denotes the scale factor at the end of reheating
which we define by the equality of the inflaton and radiation energy
densities, i.e., $\rho_{\phi}(a_{\tiny\text{rh}})=\rho_{\rm
SM}(a_{\tiny\text{rh}})$. As we see in the following, the end of reheating
roughly coincides with the condition $\langle\Gamma_{\phi}\rangle\\!\sim\\!H$.

Since the inflaton energy density dominates the total energy density during
reheating, we can solve the Friedmann equation (50), obtaining

$\displaystyle H(a)$
$\displaystyle=H_{e}\left(\frac{a_{e}}{a}\right)^{\frac{3n}{n+1}},$ for
$\displaystyle a\in[a_{e},a_{\tiny\text{rh}}],$ (59)

where $H_{e}\equiv H(a_{e})=\sqrt{\rho_{\phi_{e}}}/(\sqrt{3}M_{\rm Pl})$. It
is instructive to rewrite the Boltzmann equation for the SM energy density
(56) as

$\displaystyle\frac{d(a^{4}\rho_{\rm SM})}{da}$
$\displaystyle\simeq\frac{\langle\Gamma_{\phi\to
hh}\rangle}{H}a^{3}\rho_{\phi}.$ (60)

After using the solutions for the inflaton energy density (58) and the Hubble
rate (59) during the reheating phase, we can rewrite Eq. (60) as

$\displaystyle\frac{d(a^{4}\rho_{\rm SM})}{da}\simeq 3M_{\rm
Pl}^{2}H_{e}\langle\Gamma_{\phi\rightarrow hh}\rangle
a_{e}^{3}\,\Big{(}\frac{a}{a_{e}}\Big{)}^{\frac{3}{n+1}},$ (61)

which is straightforward to solve when the inflaton decay rate to SM Higgs
boson pairs is known. As we will show in the following subsections, in
general, the width is a function of time [12]. In the next subsection, we
calculate the inflaton “decay” to a pair of Higgs bosons, taking into account
the non-trivial $\phi$-depenent Higgs mass.

###  IV.2 Higgs boson production from oscillating inflaton

In this subsection, we calculate the energy gain $E$ of the SM radiation due
to the Higgs pair production at the expense of the inflaton field. The process
that is considered is a quantum production from a vacuum state, $\ket{0}$, to
the two-Higgs-boson final state, $\ket{\bm{h}\bm{h}}$, in the presence of the
classical inflaton field, see also [47, 13, 26, 30, 48]. Due to the energy
conservation, the energy density gained by the Higgs pair production from the
vacuum must be equal to the energy loss of the inflaton field. To put it
another way, the energy density, initially accumulated in the coherent
oscillations of the inflaton field, is transferred to the SM radiation sector
during reheating due to the cubic inflaton-Higgs coupling
$\phi|{\bm{h}}|^{2}$. Note that in this term, $\phi=\varphi\,\mathcal{P}$
should be interpreted as a given time-dependent coefficient, so effectively it
makes the Higgs boson mass vary with time. For interactions linear in the
$\phi$ field the energy gain per unit volume $V$ and per unit time due to the
pair production of final state particles $f$ with mass $m_{f}$ can be
calculated as

$\displaystyle\frac{1}{V}\frac{dE}{dt}$
$\displaystyle\\!=\\!\frac{\varphi^{2}(t)}{8\pi}\sum_{k=1}^{\infty}k\omega\lvert\mathcal{P}_{k}\rvert^{2}\Big{\lvert}\mathcal{M}_{\phi\rightarrow
ff}(k)\Big{\rvert}^{2}{\rm
Re}\\!\bigg{[}\\!\sqrt{1-\frac{4m_{f}^{2}}{k^{2}\omega^{2}}}\bigg{]},$ (62)

where $V$ is the volume factor. The matrix element squared
$\lvert\mathcal{M}_{\phi\rightarrow ff}\rvert^{2}$, summed over
spin/polarization of final states, describes the quantum process of production
of pair of particles out of the vacuum in the presence of the classical
inflaton background, while $\mathcal{P}_{k}$ denotes Fourier coefficients in
the expansion of $\mathcal{P}(t)$:

$\displaystyle\mathcal{P}(t)=\sum_{k=-\infty}^{\infty}\mathcal{P}_{k}e^{-ik\omega
t}\,,$ (63)

with the oscillation frequency/energy $\omega=2\pi/{\cal T}$ (42) and Fourier
mode number $k$,

$\displaystyle\mathcal{P}_{k}=\frac{1}{\mathcal{T}(t_{0})}\int_{t_{0}}^{t_{0}+\mathcal{T}(t_{0})}\,dt\,\mathcal{P}(t)e^{ik\omega
t}.$ (64)

Note that for $n\\!=\\!1$ the only non-zero Fourier coefficient is
$\mathcal{P}_{1}=1/2$, while for $n\\!=\\!2/3$, $n\\!=\\!3/2$ and $n\\!=\\!2$
all even coefficients are zero. Moreover, the value of $\mathcal{P}_{k}$
quickly decreases with $k$ such that the sum $\sum_{k}|\mathcal{P}_{k}|^{2}$
converges around $k$ of the order ten. The numerical values of
$\sum_{k}|\mathcal{P}_{k}|^{2}$ are collected in Table 1.

$\begin{array}[]{|c||c|c|}\hline\cr
n&\sum_{k}|\mathcal{P}_{k}|^{2}&\sum_{k}|\mathcal{P}_{k}^{2n}|^{2}\\\
\hline\cr 2/3&0.2605&0.1883\\\ 1&1/4&1/16\\\ 3/2&0.2377&0.1455\\\
2&0.2286&0.0635\\\ \hline\cr\end{array}$

Table 1: Numerical values of the summation factors appearing in Eq. (62)
(first column) and Eq. (102) (second column) for different values of $n$.

It is understood that the above process is possible only if it is
kinematically allowed, i.e., when $kw\\!>\\!2m_{f}$. Then, the collision term
in Eq. (47) could be written by defining the inflaton “decay width” as

$\frac{1}{V}\frac{dE}{dt}\equiv\rho_{\phi}\,\Gamma_{\phi}.$ (65)

Adopting the r.h.s. of the above expression, one can mimic the standard form
of the collision term for decaying particles with the “width” calculated using
the classical solution of (16). It would be an acceptable iterative procedure
that starts with the classical solution. Nevertheless, we will adopt an
alternative approach. It is important to realize that the “width” defined by
(65) non-trivially depends on $\rho_{\phi}$, i.e., on the function we are
seeking by solving the respective Boltzmann equation. Therefore we find it
more appropriate to express the whole r.h.s. of (47) in terms of the inflaton
density $\rho_{\phi}$ adopting (24). Solutions of the Boltzmann equations in
both cases approximately coincide during reheating when the “width” is smaller
than the Hubble rate. Formally the difference is of higher order in powers of
coupling constant that enters the “width”.

Our next step is to calculate the time-averaged inflaton decay width,
$\langle\Gamma_{\phi\to hh}\rangle$, using the strategy described above. For
the cubic inflaton-Higgs interaction, we obtain the following expression by
averaging over the fast oscillations of the inflaton field,

$\displaystyle\langle\Gamma_{\phi\rightarrow hh}\rangle$
$\displaystyle=\frac{g_{h\phi}^{2}}{32\pi}\frac{M_{\rm
Pl}^{2}}{m_{\phi}}\gamma_{h}\,,$ (66)

with

$\displaystyle\gamma_{h}$ $\displaystyle\\!\equiv\\!\frac{8n\sqrt{\\!\pi
n(2n\\!-\\!1)}\Gamma\\!\left(\\!\frac{n+1}{2n}\\!\right)}{\Gamma\\!\left(\frac{1}{2n}\right)}\\!\\!\sum_{i=0}^{3}\sum_{k=1}^{\infty}k\lvert\mathcal{P}_{k}\rvert^{2}\bigg{\langle}\\!\\!{\rm
Re}\\!\bigg{[}\\!\sqrt{1-\frac{4m_{h_{i}}^{2}}{k^{2}\omega^{2}}}\bigg{]}\\!\bigg{\rangle},$
(67)

where we have employed Eq. (62) and (65). Above the inflaton decay is summed
over four real components, $h_{i}$, of the Higgs doublet. Note that
$\langle\Gamma_{\phi\rightarrow hh}\rangle$ has two sources of time-
dependence:

  * (i)

the effective mass of the inflaton field $m_{\phi}$ (38),

  * (ii)

the inflaton-induced mass of the SM Higgs $m_{h_{i}}$ (44).

The time-averaged inflaton decay rate is time-independent only in the massless
reheating scenario, i.e., $m_{h_{i}}\\!=\\!0$ with the quadratic inflaton
potential, i.e., $n\\!=\\!1$. In this special case $\gamma_{h}\\!=\\!4$ due to
the four real massless scalar components of the Higgs field.

Since the energy scale of the electroweak symmetry breaking is much smaller
than the energy scale of reheating, the EW mass of the Higgs boson can be
neglected during this period. However, all Higgs doublet components, $h_{i}$,
acquire non-zero masses via their coupling to the inflaton field,

$\displaystyle m_{h_{i}}^{2}=g_{h\phi}M_{\rm
Pl}\,\varphi\begin{cases}\lvert\mathcal{P}\rvert,&\mathcal{P}>0,\,i=0,1,2,3,\\\
\\!\\!2\lvert\mathcal{P}\rvert,&\mathcal{P}<0,\,i=0,\\\
\infty,&\mathcal{P}<0,\,i=1,2,3.\end{cases}$ (68)

We see that $m_{h_{i}}^{2}$ is subjected to the short-scale oscillations of
the inflaton field through $\mathcal{P}$ and that is why we have to average
the square root in (66). Though, the summation and averaging performed in (66)
deserves an explanation. During one half of the inflaton oscillation period,
when $\mathcal{P}\\!>\\!0$, the electroweak symmetry remains unbroken, and
each Higgs doublet component receives a mass
$m_{h_{i}}\\!=\\!\sqrt{g_{{\bm{h}}}M_{\rm Pl}\varphi|\mathcal{P}|}$. In the
second half of the period, when $\mathcal{P}$ becomes negative, the physical
Higgs mass is $m_{h_{0}}\\!=\\!\sqrt{2g_{{\bm{h}}}M_{\rm
Pl}\varphi|\mathcal{P}|}$, while the three Goldstone bosons become the
longitudinal components of the SM gauge bosons and decouple from the inflaton.
In the unitary gauge that corresponds to an infinite mass of $h_{i}$ for
$i=1,2,3$, these are non-dynamical modes. The averaging performed in (66)
takes the above properties of $m_{h_{i}}^{2}$ into account. It is instructive
to write the averaged masses of the Higgs boson components as

$\langle m_{h_{i}}\\!\rangle\\!=\\!\sqrt{\\!g_{h\phi}M_{\rm
Pl}\,\varphi}\,\frac{\Gamma\\!\left(\\!\frac{3}{4n}\\!\right)\Gamma\\!\left(\\!\frac{n+1}{2n}\\!\right)}{\Gamma\\!\left(\\!\frac{2n+3}{4n}\\!\right)\Gamma\\!\left(\\!\frac{1}{2n}\\!\right)}\\!\\!\begin{cases}\\!\frac{1+\sqrt{2}}{2},&i\\!=\\!0,\\\
\\!\frac{1}{2},&i\\!=\\!1,2,3,\end{cases}$ (69)

which depends on time only through the slow varying envelope function
$\varphi(a)$ given in Eq. (31).

Note that the non-trivial Higgs mass effects appear only through the phase-
space factor $\langle\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\rangle$ in $\gamma_{h}$
(67). In particular, the phase-space kinematic effects are only relevant when
the energy $k\omega$ of the inflaton mode $k$ is larger than twice Higgs boson
mass $m_{h_{i}}$. We also note that for a general value of $n$ the dominant
contribution originates from the first Fourier mode, i.e., $k\\!=\\!1$, as the
amplitude ${\cal P}_{k}$ drops quickly for larger mode numbers. Therefore it
is instructive to write the ratio

$\displaystyle\frac{\langle m_{h}\rangle}{\omega}$
$\displaystyle\\!\simeq\\!\frac{\sqrt{\\!g_{h\phi}M_{\rm
Pl}\,\varphi}}{m_{\phi}}\sqrt{\frac{2n-1}{\pi\,n}}\frac{\Gamma\\!\left(\\!\frac{3}{4n}\\!\right)}{\Gamma\\!\left(\\!\frac{2n+3}{4n}\\!\right)},$
$\displaystyle\\!=\\!\frac{\sqrt{\\!g_{h\phi}}}{\sqrt{2\pi}\,n}\frac{M_{\rm
Pl}^{2}}{\Lambda^{2}}\frac{\Gamma\\!\left(\\!\frac{3}{4n}\\!\right)}{\Gamma\\!\left(\\!\frac{2n+3}{4n}\\!\right)}\Big{(}\frac{\varphi_{e}}{M_{\rm
Pl}}\Big{)}^{\\!\frac{3-2n}{2}}\Big{(}\frac{a}{a_{e}}\Big{)}^{\\!\frac{3(2n-3)}{2(n+1)}},$
(70)

where for brevity we have neglected the Higgs field $h_{i}$ component-
dependent contribution $(1\\!+\\!\sqrt{2})/2$ for $i=0$ and $1/2$ for
$i=1,2,3$. In the second line of Eq. (70) we have employed the solution for
$\varphi(a)$ (31), which is valid during reheating, i.e.,
$\langle\Gamma_{\phi\rightarrow hh}\rangle\\!\ll\\!H$, and the expression for
the inflaton mass $m_{\phi}$ (39). Note that the time-dependence of $\langle
m_{h}\rangle/\omega$ through the scale factor $a$ disappears for
$n\\!=\\!3/2$, while for $n\\!<\\!3/2$ ($n\\!>\\!3/2$) it decreases
(increases) with time during reheating. This implies that the kinematic
suppression amplifies over time for $n\\!>\\!3/2$ and stays constant for
$n\\!=\\!3/2$. However, for $n\\!<\\!3/2$, the role of the kinematic
suppression becomes less and less important as the Universe expands.

The crucial factor in the calculation of the inflaton width (66) is
$\langle\\!\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\rangle$, which captures the
kinematic phase-space suppression due to non-trivial Higgs mass. We find the
phase-space factor after time-averaging over fast inflaton oscillations as

$\displaystyle\bigg{\langle}\\!\\!\sqrt{1\\!-\\!\frac{4m_{h_{i}}^{2}}{k^{2}\omega^{2}}}\\!\bigg{\rangle}$
$\displaystyle\approx\bigg{\lceil}\\!\frac{9\,\kappa_{i}\,n\Gamma\\!\left(\frac{n+1}{2n}\right)}{16\sqrt{\pi}\,\Gamma\\!\left(\frac{1}{2n}\right)}\frac{k^{2}\omega^{2}}{g_{h\phi}M_{\rm
Pl}\varphi}\\!\bigg{\rceil},$
$\displaystyle=\bigg{\lceil}\frac{9\sqrt{\pi}\kappa_{i}n^{2}\Gamma\\!\left(\frac{n+1}{2n}\right)^{3}}{16(2n-1)\Gamma\\!\left(\frac{1}{2n}\right)^{3}}\frac{k^{2}m_{\phi}^{2}}{g_{h\phi}M_{\rm
Pl}\varphi}\bigg{\rceil},$ (71)

where for $i\\!=\\!0$ component of the Higgs field $\kappa_{0}\\!=\\!1$,
whereas for $i\\!=\\!1,2,3$ components $\kappa_{1,2,3}\\!=\\!1/2$ due to the
fact that these modes become Goldstone modes in the broken phase. Above
‘ceiling’ function $\lceil\cdots\rceil$ is defined as follows:

$\lceil x\rceil=\begin{cases}1,&x\geq 1,\\\ x,&x<1,\end{cases}$ (72)

where $\lceil\cdots\rceil=1$ implies massless final states, i.e., no kinematic
suppression. Note that in the above expression, for a generic $n$, the time
dependence enters through the inflaton envelope function $\varphi(a)$ (31) and
the inflaton mass $m_{\phi}$ (39). We can also rewrite Eq. (71) as a function
of the scale factor:

$\displaystyle\bigg{\langle}\\!\\!\sqrt{1\\!-\\!\frac{4m_{h_{i}}^{2}}{k^{2}\omega^{2}}}\\!\bigg{\rangle}\\!\approx\\!\bigg{\lceil}\frac{g_{k}}{g_{h\phi}}\Big{(}\\!\frac{a_{e}}{a}\\!\Big{)}^{\\!\frac{3(2n-3)}{(n+1)}}\bigg{\rceil},$
(73)

where we have used $M\\!=\\!M_{\rm Pl}$ ($\alpha\\!=\\!1/6$),
$\varphi_{e}\\!\sim\\!M_{\rm Pl}$ and defined

$\displaystyle
g_{k}\\!\equiv\\!\frac{9\sqrt{\pi}\kappa_{i}k^{2}n^{3}\Lambda^{4}\,\Gamma\\!\left(\\!\frac{n+1}{2n}\\!\right)^{3}}{8\,M_{\rm
Pl}^{4}\,\Gamma\\!\left(\\!\frac{1}{2n}\\!\right)^{3}}.$ (74)

From Eq. (73) we see that the phase-space factor
$\langle\\!\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\rangle$ increases (decreases)
with time for $n\\!<\\!3/2$ ($n\\!>\\!3/2$), whereas it stays constant for
$n\\!=\\!3/2$ as shown in Fig. 5.

Figure 5: Evolution of the phase-space factor
$\big{\langle}\\!\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\big{\rangle}$ for the
inflaton decay to a pair Higgs bosons as a function of the scale factor $a$
for $k\\!=\\!1$ and $g_{h\phi}\\!=\\!10^{-5}$ and $10^{-10}$ in the left- and
right-panels, respectively. Solid (dashed) lines are for the $i\\!=\\!0$
($1,2,3$) component of the SM Higgs doublet.

Moreover, one can estimate the scale factor $a_{k}$ at which the kinematic
suppression becomes irrelevant (for $n\\!<\\!3/2$) or relevant (for
$n\\!>\\!3/2$) by the condition $\lceil\cdots\rceil=1$ in the above expression
as

$a_{k}\equiv
a_{e}\bigg{(}\frac{g_{k}}{g_{h\phi}}\bigg{)}^{\frac{(n+1)}{3(2n-3)}},\qquad
n\neq 3/2.$ (75)

Whereas for the special case $n\\!=\\!3/2$ the phase-space (73) suppression is
a time-independent constant factor when $g_{h\phi}>g_{k}$ for $k=1$. Note that
depending on the strength of $g_{h\phi}$, the kinematic suppression may end
during the reheating phase for $n\\!<\\!3/2$. In this case, for scale factor
$a_{k}\\!<\\!a\\!\leq\\!a_{\rm rh}$ the reheating dynamics would be the same
as that of massless Higgs final states. On the other hand, the kinematic
suppression becomes more and more significant with time for $n>3/2$. This is
because even if the averaged Higgs mass was much smaller than the inflaton
energy at the onset of reheating, it would soon become of the order of
$k\omega$, and hence the kinematic effects would be important until the end of
the reheating phase. With the above definition of $a_{k}$, we can recast the
phase-phase factor (73) as

$\displaystyle\bigg{\langle}\\!\\!\sqrt{1\\!-\\!\frac{4m_{h_{i}}^{2}}{k^{2}\omega^{2}}}\\!\bigg{\rangle}\\!\approx\\!\begin{dcases}\\!\\!\bigg{\lceil}\Big{(}\\!\frac{a_{k}}{a}\\!\Big{)}^{\\!\frac{3(2n-3)}{(n+1)}}\bigg{\rceil},&n\\!\neq\\!3/2,\\\
\\!\\!\bigg{\lceil}\frac{g_{k}}{g_{h\phi}}\bigg{\rceil},&n\\!=\\!3/2.\end{dcases}$
(76)

In Fig. 5 we have plotted exact numerical results for the phase-space factor
$\langle\\!\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\rangle$ with $k\\!=\\!1$ as a
function of the scale factor $a$ for two values of the inflaton-Higgs coupling
$g_{h\phi}=10^{-5}$ (left panel) and $10^{-10}$ (rightpanel). First of all,
let us note that if the mass of the SM Higgs is small as compared to the
inflaton mode energy $k\omega$, the phase-space factor
$\langle\\!\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\rangle\\!=\\!1$. For
$n\\!<\\!3/2$ and sufficient strong inflaton-Higgs coupling, e.g.,
$g_{h\phi}\\!=\\!10^{-5}$, the
$\langle\\!\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\rangle$ factor initially grows
with the scale factor $a$, until the inflaton-induced mass of the Higgs boson
becomes negligible and the kinematical suppression is no longer relevant. From
the left-panel of Fig. 5 for $g_{h\phi}\\!=\\!10^{-5}$ the number of e-folds
for kinematic suppression is $N_{k}\\!\sim\\!4.5$ and $9$ for $n\\!=\\!2/3$
and $n\\!=\\!1$, respectively. As noted above for $n\\!=\\!3/2$, the phase-
space factor is time-independent and is given by
$\langle\\!\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\rangle\\!\sim\\!\lceil\kappa_{i}\Lambda^{4}/(4g_{h\phi}M_{\rm
Pl}^{4})\rceil$. This means that the kinematical suppression is present during
the whole period of reheating unless the $g_{h\phi}$ coupling is very weak and
the ceiling function approaches 1. For $n\\!>\\!3/2$, the role of the non-zero
mass of the Higgs boson increases with time so that
$\langle\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\rangle$ decreases as the scale
factor increases. This means that even if the effects of the non-zero Higgs
mass are negligible at the onset of reheating period with increasing time, the
kinematical suppression becomes more and more relevant. The numerical results
for the phase-space factor agree with our approximate analytic results (73) to
very good precision.

In Fig. 6 we show the evolution of the mass of the 0-component of the Higgs
doublet, $m_{h_{0}}$ (68), its average $\langle m_{h_{0}}\rangle$ (69), and
the effective inflaton energy $\omega$ (42) as a function of the scale factor
during the first few e-folds of reheating. In these plots we have fix
$n\\!=\\!2/3$ (left-panel) and $3/2$ (right-panel) for two values of
$g_{h\phi}\\!=\\!10^{-5},10^{-7}$. Note that kinematic phase-space suppression
discussed above is relevant for regions where $\langle
m_{h_{0}}\rangle\\!\gtrsim\\!m_{\phi}$. Furthermore, in Fig. 6 we also have
compared $m_{h_{0}}(a)$ with the evolution of the thermal bath temperature $T$
(see below for details). When the Higgs mass is larger than the SM bath
temperature it behaves as a non-relativistic matter. We note that the
$\phi$-dependent mass of the SM Higgs drops below the temperature, roughly
corresponding to the $a\simeq a_{k}$.

Figure 6: The evolution of the Higgs mass $m_{h_{0}}(a)$ for two values of the
inflaton-Higgs coupling: $g_{h\phi}=10^{-5}$ (blue) and $g_{h\phi}=10^{-7}$
(red) for $n\\!=\\!2/3$ and $3/2$ in the left and right-panel, respectively.
The colored (blue, red) dashed curves show the evolution of the averaged Higgs
mass, i.e., $\langle m_{h_{0}}(a)\rangle$. The thick solid purple and light
blue lines denote $T(a)$ for $g_{h\phi}=10^{-5}$ and $10^{-7}$, respectively,
whereas the loosely dashed green line illustrates the evolution of the
inflaton mass $m_{\phi}$.

###  IV.3 Reheating through massless Higgs bosons

In this subsection, we consider the reheating dynamics by neglecting the
$\phi$-dependent Higgs mass, which we refer to as the massless reheating
scenario. In other words, we neglect the contribution to the Higgs boson mass
originating from the cubic interaction term $\phi|{\bm{h}}|^{2}$ in the
Lagrangian (6). This is, of course, merely a reference point introduced to
illustrate kinematical mass effects in the realistic case, i.e., with
$m_{h_{i}}\neq 0$. As shown in the following subsection, there are significant
differences between the massless scenario and the case where we consider the
inflaton-induced Higgs mass. Our main goal is to find an approximate
analytical solution to the Boltzmann equation for the SM radiation (61), which
then allows us to determine the temperature evolution during the reheating
phase. In the massless reheating scenario, there is no kinematical suppression
of the Higgs boson production. Thus, the only source of time-dependence of the
inflaton decay rate (66) is through the effective inflaton mass $m_{\phi}$
(39). In this case, the inflaton decay can be written as

$\displaystyle\langle\Gamma_{\phi\rightarrow
hh}^{(0)}(a)\rangle=\widetilde{\Gamma}_{\phi\rightarrow
hh}^{(0)}\Big{(}\frac{a}{a_{e}}\Big{)}^{\frac{3(n-1)}{n+1}},$ (77)

with the constant width

$\displaystyle\widetilde{\Gamma}_{\phi\rightarrow hh}^{(0)}$
$\displaystyle\equiv\frac{g_{h\phi}^{2}}{\sqrt{2\pi}}\frac{n\Gamma\left(\frac{n+1}{2n}\right)}{\Gamma\left(\frac{1}{2n}\right)}\frac{M_{\rm
Pl}^{3}}{\Lambda^{2}}\Big{(}\frac{M_{\rm
Pl}}{\varphi_{e}}\Big{)}^{\\!n-1}\sum_{k=1}^{\infty}k|{\cal P}_{k}|^{2},$ (78)

where superscript “$(0)$” denotes the massless Higgs boson case. Note that for
a generic value of $n\\!\neq\\!1$ the inflaton decay width is a function of
time which can lead to non-trivial consequences for reheating dynamics [12].
In particular, for $n\\!<\\!1$, $\Gamma_{\phi\rightarrow hh}^{(0)}$ decreases
with the growth of the scale factor, while for $n\\!>\\!1$ it increases.
Therefore, one can expect that for fixed values of the model parameters
($g_{h\phi},\Lambda$), the reheating is the most efficient for $n\\!>\\!1$,
whereas for $n\\!<\\!1$, the duration of reheating is prolonged. Moreover, for
the quadratic inflaton potential, i.e., $n\\!=\\!1$ we reproduce the well-
known expression for the time-independent decay width, i.e.,
$\Gamma_{\phi\rightarrow hh}^{(0)}\\!=\\!(g_{h\phi}M_{\rm Pl})^{2}/(8\pi
m_{\phi})$, with $m_{\phi}\\!=\\!\sqrt{2}\Lambda^{2}/M$ and ${\cal
P}_{k}\\!=\\!\delta_{1k}/2$.

Inserting formula (77) into Eq. (61) one obtains the following solution for
the SM radiation energy density,

$\displaystyle\rho_{\rm SM}^{(0)}(a)\\!=\\!3M_{\rm
Pl}^{2}H_{e}\,\widetilde{\Gamma}_{\phi\rightarrow
hh}^{(0)}\,\frac{n+1}{4n+1}\\!\bigg{[}\\!\Big{(}\frac{a_{e}}{a}\Big{)}^{\\!\frac{3}{n+1}}-\Big{(}\frac{a_{e}}{a}\Big{)}^{\\!4}\bigg{]},$
(79)

for $a\\!\in\\![a_{e},a_{\tiny\text{rh}}]$. Note that for positive $n$, the
first term in the square brackets above dominates over the second one during
the reheating phase, i.e., $a_{e}\\!<\\!a\leq a_{\rm rh}$. Consequently, the
temperature of the thermal bath (53) during reheating evolves as

$\displaystyle T^{(0)}(a)$ $\displaystyle\simeq\bigg{[}\frac{90(n+1)M_{\rm
Pl}^{2}H_{e}\,\widetilde{\Gamma}_{\phi\rightarrow
hh}^{(0)}}{(4n+1)\pi^{2}\,g_{\star}(T)}\bigg{]}^{\\!1/4}\Big{(}\frac{a_{e}}{a}\Big{)}^{\frac{3}{4(n+1)}}.$
(80)

For $n\\!=\\!1$ we obtain standard scaling of temperature w.r.t. the scale
factor, i.e., $T\\!\propto\\!a^{-3/8}$, however for $n\\!<\\!1$ ($n\\!>\\!1$)
the dependence of the temperature on the scale factor is more (less) steeper.
The explicit dependence on the scale factor of the SM energy density
$\rho_{\rm SM}(a)$ and temperature $T(a)$ for different $n$ is listed in Table
2.

Another quantity of physical importance is the maximum temperature of the SM
bath $T_{{\rm max}}^{(0)}$. In the massless reheating scenario, $T_{{\rm
max}}^{(0)}$ is typically reached just after the end of inflation, at

$\displaystyle a_{{\rm
max}}^{(0)}=a_{e}\bigg{(}\\!\frac{4(n+1)}{3}\\!\bigg{)}^{\\!\frac{n+1}{4n+1}},$
(81)

Thus, the maximum temperature in the massless reheating scenario is

$\displaystyle T_{{\rm max}}^{(0)}$
$\displaystyle\\!=\\!\bigg{[}\\!\frac{90(n+1)M_{\rm
Pl}^{2}H_{e}\,\widetilde{\Gamma}_{\phi\rightarrow
hh}^{(0)}}{(4n+1)\pi^{2}\,g_{\star}(T)}\\!\bigg{]}^{\\!1/4}\\!\\!\bigg{(}\\!\frac{3}{4(n+1)}\\!\bigg{)}^{\\!\frac{3}{4(4n+1)}}.$
(82)

###  IV.4 Reheating through massive Higgs bosons

In this subsection, we discuss the reheating through inflaton decays to Higgs
boson pairs, taking into account the effects of the inflaton-induced Higgs
boson mass on the reheating dynamics. The important difference of this
scenario as compared to the reheating with the massless Higgs, discussed
above, is the non-trivial kinematical phase-space suppression of the inflaton
decay width, $\langle\Gamma_{\phi\rightarrow hh}\rangle$ (66) through
$\gamma_{h}$ factor (67). In particular, the main difference is the time-
averaged phase-space factor
$\langle\\!\sqrt{1-(2m_{h_{i}}/k\omega)^{2}}\rangle$ that generates a
kinematical suppression due to the non-zero Higgs mass, as discussed above in
Sec. IV.2. The explicit form of the time-dependence of the phase-space factor
has been calculated in Eq. (73) for the four components of the SM Higgs
doublet. Moreover, as it has been pointed out above, the $\gamma_{h}$ factor
(67) is time-independent during the reheating phase for $n\\!=\\!3/2$, whereas
for $n\\!<\\!3/2\,(n\\!>\\!3/2)$, $\gamma_{h}$ increases (decreases) during
reheating. For $n\\!<\\!3/2$ the kinematic suppression gradually disappears as
the averaged Higgs mass becomes smaller compared to the inflaton energy, i.e.,
$\langle m_{h}\rangle\lesssim\omega$. The scale factor at which the kinematic
suppression vanishes is denoted by $a_{k}$. Note that, in principle, the
kinematic suppression can be present throughout the reheating phase, or it may
vanish during this phase, i.e., $a_{k}\\!<\\!a_{\rm rh}$. In the latter case,
the inflaton-induced mass of the Higgs field suppresses the radiation
production for $a_{e}\\!<\\!a\\!\leq\\!a_{k}$, while for $a_{k}<a\leq a_{\rm
rh}$, the kinematic suppression is irrelevant, and the reheating dynamics
follows the same evolution as that of the massless reheating scenario
discussed above. On the other hand, for $n\\!>\\!3/2$ it can happen that at
the onset of reheating, the kinematical suppression has a negligible impact on
the reheating dynamics, but at some point, $a_{k}$, when the averaged Higgs
mass approaches the inflaton energy, i.e., $\langle m_{h}\rangle\sim\omega$,
it starts to affect the reheating dynamics until the end of reheating phase.

The averaged inflaton decay width Eq. (66) with inflaton-induced Higgs mass is
given by,

$\displaystyle\langle\Gamma_{\phi\rightarrow hh}(a)\rangle$
$\displaystyle\simeq\widetilde{\Gamma}_{\phi\rightarrow
hh}\Big{(}\frac{a}{a_{e}}\Big{)}^{\frac{3(n-1)}{n+1}}$
$\displaystyle\times\sum_{k=1}^{\infty}k|{\cal
P}_{k}|^{2}\begin{dcases}\\!\\!\bigg{\lceil}\\!\Big{(}\\!\frac{a_{k}}{a}\\!\Big{)}^{\\!\frac{3(2n-3)}{(n+1)}}\bigg{\rceil},&n\\!\neq\\!3/2,\\\
\\!\\!\bigg{\lceil}\frac{g_{k}}{g_{h\phi}}\bigg{\rceil},&n\\!=\\!3/2,\end{dcases}$
(83)

where constant width $\widetilde{\Gamma}_{\phi\rightarrow hh}$ is defined as

$\displaystyle\widetilde{\Gamma}_{\phi\rightarrow hh}$
$\displaystyle\equiv\frac{g_{h\phi}^{2}}{\sqrt{2\pi}}\frac{n\Gamma\left(\frac{n+1}{2n}\right)}{\Gamma\left(\frac{1}{2n}\right)}\frac{M_{\rm
Pl}^{3}}{\Lambda^{2}}\Big{(}\frac{M_{\rm Pl}}{\varphi_{e}}\Big{)}^{\\!n-1}.$
(84)

Consequently, in the massive reheating scenario, the Boltzmann equation for
the SM radiation energy density (61) can be analytically solved with the
inflaton decay rate (83). In the following, we calculate the SM radiation
energy density and temperature for three different cases,
$n\\!\lessgtr\\!3/2$, and $n=3/2$.

In this case when $n\\!\lessgtr\\!3/2$ for $a\\!\gtrless\\!a_{k}$ the inflaton
decay rate to the SM Higgs boson (83) has no phase-space suppression, and it
is equal to that of the massless Higgs boson case Eq. (77). Therefore for
$n\\!\lessgtr\\!3/2$ and $a\\!\gtrless\\!a_{k}$ the SM radiation energy
density is the same as the massless case Eq. (79). However, for
$a\\!\lessgtr\\!a_{k}$ with $n\\!\lessgtr\\!3/2$ we have a phase-space
kinematic suppression in the inflaton decay rate (83) and, therefore, the SM
energy density is modified in comparison to the massless case Eq. (79).
Similarly, for the case when $n=3/2$, there is no time-dependent kinematical
suppression for the inflaton decay rate; however, there is a constant
suppression proportional to $\lceil g_{k}/g_{h\phi}\rceil$. Therefore, for
$n=3/2$, the time dependence of the SM energy density and the temperature
would be the same as that of the massless case up to constant kinematical
suppression. Hence we can write the SM energy density for different values of
$n$ as

$\displaystyle\rho_{\rm SM}(a)$ $\displaystyle\simeq 3M_{\rm
Pl}^{2}H_{e}\widetilde{\Gamma}_{\phi\rightarrow
hh}\begin{dcases}\\!\frac{n+1}{2(5-n)}\\!\bigg{[}\\!\Big{(}\frac{a_{e}}{a}\Big{)}^{\\!\frac{6(n-1)}{n+1}}-\Big{(}\frac{a_{e}}{a}\Big{)}^{\\!4}\bigg{]}\\!\sum_{k=1}^{\infty}k|{\cal
P}_{k}|^{2}\Big{(}\frac{a_{k}}{a_{e}}\Big{)}^{\\!\frac{3(2n-3)}{(n+1)}},\quad&a\\!\lessgtr\\!a_{k},\,n\\!\lessgtr\\!3/2,\\\
\frac{n+1}{4n+1}\bigg{[}\\!\Big{(}\frac{a_{e}}{a}\Big{)}^{\\!\frac{3}{n+1}}-\Big{(}\frac{a_{e}}{a}\Big{)}^{\\!4}\bigg{]}\\!\sum_{k=1}^{\infty}k|{\cal
P}_{k}|^{2},&a\\!\gtrless\\!a_{k},\,n\\!\lessgtr\\!3/2,\\\
\frac{n+1}{4n+1}\bigg{[}\\!\Big{(}\frac{a_{e}}{a}\Big{)}^{\\!\frac{3}{n+1}}-\Big{(}\frac{a_{e}}{a}\Big{)}^{\\!4}\bigg{]}\\!\sum_{k=1}^{\infty}k|{\cal
P}_{k}|^{2}\Big{\lceil}\frac{g_{k}}{g_{h\phi}}\Big{\rceil},&\qquad
n\\!=\\!3/2,\end{dcases}$ (85)

where the time-dependence through the scale factor is dominated by the first
terms in the square brackets above. Thus, during the kinematical suppression
phase, i.e., $a\\!\lessgtr\\!a_{k}$ for $n\\!\lessgtr\\!3/2$, the SM radiation
energy density and the temperature scale as

$\displaystyle\rho_{\rm SM}$ $\displaystyle\propto a^{\frac{6(1-n)}{n+1}}\,,$
$\displaystyle T$ $\displaystyle\propto a^{\frac{3(1-n)}{2(n+1)}},$
$\displaystyle a$ $\displaystyle\\!\lessgtr\\!a_{k},\,n\\!\lessgtr\\!3/2.$
(86)

Whereas, during the region of parameters where the time-dependent kinematical
suppression is absent, we get the scaling as that of the massless Higgs case,
i.e.,

$\displaystyle\rho_{\rm SM}\\!\sim\\!\rho_{\rm SM}^{(0)}$
$\displaystyle\propto a^{\frac{-3}{n+1}}\,,$ $\displaystyle T$
$\displaystyle\\!\sim\\!T^{(0)}\propto a^{\frac{-3}{4(n+1)}},$ $\displaystyle
a$ $\displaystyle\\!\gtrless\\!a_{k},\,n\\!\lesseqgtr\\!3/2.$ (87)

It is important to note that the above scaling behavior is very unusual for
the SM radiation energy density. Thus, the temperature of the thermal bath
during reheating when the kinematical suppression is present. It turns out
that in the massive scenario for $n\\!=\\!2/3$, both $\rho_{\rm SM}$ and $T$
increases with time, whereas for $n\\!=\\!1$ they are nearly constant. As we
have discussed earlier on, for $n\\!=\\!3/2$ the inflaton decay rate depends
on time only though $m_{\phi}$ and thus, in both massless and massive cases,
we obtain the same scaling. Finally, for $n\\!>\\!3/2$, the $\rho_{\rm SM}$
dependence on the scale factor is steeper in the massive reheating scenario in
comparison to the massless Higgs case. The approximate scalings of averaged
inflaton decay width, $\langle\Gamma_{\phi\rightarrow hh}\rangle$, the SM
radiation energy density, $\rho_{\rm SM}$, and the temperature, $T$, with the
scale factor, $a$, are collected in Table 2 for the massless (when the
kinematic suppression is irrelevant) and massive reheating scenarios (when the
kinematic suppression is relevant).

$\begin{array}[]{|c||c|c|c||c|c|c|}\hline\cr n&\langle\Gamma_{\phi\rightarrow
hh}^{(0)}\rangle&\rho_{\rm SM}^{(0)}&T^{(0)}&\langle\Gamma_{\phi\rightarrow
hh}\rangle&\rho_{\rm SM}&T\\\ \hline\cr
2/3&a^{-3/5}&a^{-9/5}&a^{-9/20}&a^{12/5}&a^{6/5}&a^{3/10}\\\
1&a^{0}&a^{-3/2}&a^{-3/8}&a^{3/2}&a^{0}&a^{0}\\\
3/2&a^{3/5}&a^{-6/5}&a^{-3/10}&a^{3/5}&a^{-6/5}&a^{-3/10}\\\
2&a^{1}&a^{-1}&a^{-1/4}&a^{0}&a^{-2}&a^{-1/2}\\\ \hline\cr\end{array}$

Table 2: Approximate scaling of the averaged inflaton decay width,
$\langle\Gamma_{\phi\rightarrow hh}\rangle$, the SM radiation energy density,
$\rho_{\rm SM}$, and the temperature, $T$, with the scale factor, $a$, in the
massless (with superscript $(0)$) and massive reheating scenario when the
kinematic suppression is relevant.

###  IV.5 Tachyonic resonant production of Higgs bosons

Before we go onto the numerical analysis of the reheating dynamics we would
like to make some comments here regarding the possible tachyonic resonant
production of Higgs boson during the early phases of reheating. In our model
Higgs boson Lagrangian has the form,

${\cal
L}\supset\sqrt{-g}\Big{[}g^{\mu\nu}\left(\partial_{\nu}\bm{h}\right)^{\dagger}\left(\partial_{\mu}\bm{h}\right)-\mu_{h}^{2}|\bm{h}|^{2}-\lambda_{h}|\bm{h}|^{4}\Big{]}$
(88)

where the inflaton-induced Higgs mass $\mu_{h}^{2}\equiv
g_{h\phi}M_{\mathrm{Pl}}\phi$ (43) could be larger than the inflaton effective
mass $m_{\phi}^{2}$ (38) at the onset of reheating phase, i.e.
$\phi\\!=\\!\phi_{e}\\!\sim\\!M_{\rm Pl}$ for values of $g_{h\phi}$ satisfying
(13) and (14), i.e.

$\sqrt{\lambda_{h}}\left(\frac{\Lambda^{2}}{\phi M_{\rm Pl}}\right)\gtrsim
g_{h\phi}\gtrsim\frac{3}{4}\left(\frac{\Lambda^{2}}{\phi M_{\rm
Pl}}\right)^{2}\left(\frac{\phi}{M_{\rm Pl}}\right).$ (89)

for $\alpha=1/6$. Therefore tachyonic resonant production of Higgs boson can
constitute an efficient source of (p)reheating [49, 50]. This can be easily
seen from the equation of motion of the Higgs field in the unitary gauge as,

$\left(\frac{d^{2}}{dt^{2}}-\frac{\nabla^{2}}{a^{2}}+3H\frac{d}{dt}+g_{h\phi}M_{\mathrm{Pl}}\phi\right)h_{0}=-\lambda_{h}h_{0}^{3}$
(90)

where the non-linear term is present due to Higgs self-interactions. General
solution of the above non-linear equation with the expanding background is
difficult to handle and require lattice simulations. However, in the following
we employ the Hartree approximation to replace the non-linear term by a linear
term as $\lambda_{h}h_{0}^{3}\to 3\lambda_{h}\langle h_{0}^{2}\rangle h_{0}$,
where $\langle h_{0}^{2}\rangle$ is variance of the Higgs field which can be
calculated using the linear solution (or iteratively in orders of non-
linearity parameter $\lambda_{h}$) for the Higgs mode function $h_{0,k}$ with
mode momentum $k$ as,

$\left\langle
h_{0}^{2}\right\rangle(t)=\frac{1}{2\pi^{2}}\int_{0}^{\infty}dk\,k^{2}\left|h_{0,k}(t)\right|^{2}.$
(91)

The Higgs mode function $h_{0,k}$ satisfies the mode equation

$\ddot{h}_{0,k}+3H\dot{h}_{0,k}+\left(\frac{k^{2}}{a^{2}}+g_{h\phi}M_{\mathrm{Pl}}\phi+3\lambda_{h}\left\langle
h_{0}^{2}\right\rangle\right)h_{0,k}=0,$ (92)

where we have adopted the Hartree approximation. Applying the field
redefinition, $\tilde{h}_{0,k}\equiv a^{3/2}h_{0,k}$, we can recast the above
mode equation as

$\ddot{\tilde{h}}_{0,k}+\omega_{k}^{2}\tilde{h}_{0,k}=0,$ (93)

where the dispersion relation is

$\omega_{k}^{2}\equiv\frac{k^{2}}{a^{2}}+\frac{9}{4}H^{2}w+g_{h\phi}M_{\rm
Pl}\phi+3\lambda_{h}\left\langle h_{0}^{2}\right\rangle,$ (94)

with $w=p/\rho$ which can be in the range $w\in[-1,1]$ for our choice of
inflaton potential (4). In the above dispersion relation $\omega_{k}^{2}$
second and third terms can be negative and hence source the non-perturbative
production of Higgs modes, i.e. preheating dynamics. However, note that second
term $H^{2}w$ is much smaller compared to the inflaton-induced Higgs mass when
$\phi$ is oscillating with an amplitude $\sim\\!M_{\rm Pl}$ and $g_{h\phi}$
satisfying (89). Therefore tachyonic resonant Higgs production is mainly
driven by large negative Higgs mass term $\mu_{h}^{2}\sim-g_{h\phi}M_{\rm
Pl}|\phi|$. However, the positive large self-interaction contribution term
$\lambda_{h}\left\langle h_{0}^{2}\right\rangle$ to mode frequency shuts off
the tachyonic production for any $k$ mode if

$3\lambda_{h}\left\langle h_{0}^{2}\right\rangle\gtrsim|\mu_{h}^{2}|.$ (95)

This gives a lower-bound on the Higgs quartic coupling for which the tachyonic
resonant production is irrelevant for $\mu_{h}^{2}<0$. We calculate the
variance of the Higgs field $\langle h_{0}^{2}\rangle$ (91) for
$\mu_{h}^{2}<0$ by using the leading linear solution for the mode function
$h_{0,k}(t)$, i.e. neglecting $\lambda_{h}$ term in the mode equation (92),
see also [49], as

$\langle h_{0}^{2}\rangle\sim\frac{|\mu_{h}|\Lambda^{2}}{16\pi^{2}M_{\rm
Pl}}\exp\\!\bigg{(}\frac{4|\mu_{h}|M_{\rm Pl}}{\Lambda^{2}}\bigg{)}.$ (96)

Hence the lower-bound on the Higgs quartic coupling (95) reads,

$\displaystyle\lambda_{h}$
$\displaystyle\gtrsim\frac{16\pi^{2}}{3}\frac{\sqrt{g_{h\phi}|\phi|}M_{\rm
Pl}^{3/2}}{\Lambda^{2}}\exp\\!\bigg{(}\\!\\!-\\!\frac{4\sqrt{g_{h\phi}|\phi|}M_{\rm
Pl}^{3/2}}{\Lambda^{2}}\bigg{)}.$ (97)

For the values of Higgs quartic coupling employed throughout this work,
$\lambda_{h}\sim 0.1$, the lower limit for $g_{h\phi}$ to block the tachyonic
resonant production of Higgs boson reads,

$g_{h\phi}\gtrsim\frac{4\Lambda^{4}}{M_{\rm Pl}^{3}|\phi|}.$ (98)

Note that the above lower limit for $g_{h\phi}$ is slightly stronger than the
one in (89). On the other hand we find that for

$\frac{3}{4}\left(\frac{\Lambda^{4}}{M_{\rm Pl}^{3}|\phi|}\right)\lesssim
g_{h\phi}\lesssim 4\left(\frac{\Lambda^{4}}{M_{\rm Pl}^{3}|\phi|}\right).$
(99)

and $\lambda_{h}\sim 0.1$ tachyonic resonant production of Higgs modes is
possible. A detailed analysis of this possibility is an interesting
possibility, however it is beyond the scope of present work.

To summarize, our analysis above (with an order of magnitude approximations)
concludes that in the parameter space considered in this work (with
$g_{h\phi}$ respecting lower bound (98)) the tachyonic resonant Higgs
production is inefficient due to Higgs self-interactions. Therefore in the
following numerical analysis we only focus on the perturbative production of
Higgs bosons as discussed in the subsections above.

###  IV.6 Numerical analysis

In this subsection, we present a numerical analysis of the reheating dynamics
due to the inflaton decays to the SM Higgs boson through cubic interaction of
the form $g_{h\phi}M_{\rm Pl}\,\phi|{\bm{h}}|^{2}$ (6). In particular, we
provide solutions of the first two Boltzmann equations (47) and (48) for the
inflaton $\rho_{\phi}(a)$ and SM radiation $\rho_{\rm SM}(a)$ energy
densities. The reheating dynamics depend only on the form of the inflaton
potential and inflaton-Higgs interaction term 666We assume that the inflaton
interactions with other SM fields and DM are much smaller than the leading
inflaton-Higgs coupling through the dim-3 operator.. For the numerical
analysis, we consider the $\alpha$-attractor T-model of inflaton with
potential (4), where we have fixed $\alpha=1/6$, such that $M=M_{\rm Pl}$, and
the scale of inflation $\Lambda=3\times 10^{-3}M_{\rm Pl}$. Whereas, the
benchmark values for $n$ have been chosen to be $2/3,1,3/2,2$ which correspond
to the equation of state during the reheating phase $\bar{w}=-1/5,0,1/5,1/3$,
respectively. Furthermore, we consider two values for the inflaton-Higgs
interaction $g_{h\phi}\\!=\\!10^{-5}$ and $10^{-10}$, which are close to its
upper and lower bounds, see Eq. (15).

Figure 7: Solutions of the first two Boltzmann equations (47) and (48) for the
inflaton $\rho_{\phi}$ and radiation $\rho_{\rm SM}$ energy densities are
plotted in the upper panel as a function of the scale factor $a$. The empty
and filled dots indicate the end of reheating if produced Higgs boson pairs
are massless and massive, respectively. The lower panel shows thermal bath
temperature $T$ (blue), the Hubble parameter $H$ (green) and the averaged
inflaton width $\langle\Gamma_{\phi\rightarrow hh}\rangle$ (red) for the
massless and massive cases.The empty and filled stars show where the maximal
temperature during reheating has been obtained for the massless and massive
cases, respectively. Adopted parameters are specified above.

Figure 8: Same as in Fig. 7 except for $g_{h\phi}=10^{-10}$.

In the upper panel of Fig. 7 and Fig. 8 we have plotted the numerical
solutions for the inflaton (orange) and radiation (purple) energy densities
for fixed $g_{h\phi}=10^{-5}$ and $10^{-10}$ in Fig. 7 and Fig. 8,
respectively. The lower panel of each figure presents the evolution of the
thermal bath temperature (blue), the Hubble rate (green), and the averaged
inflaton decay width (red) as functions of the scale factor. The solid and
dashed curves correspond to the massive and massless reheating scenario,
respectively. The empty (filled) dots indicate the end of reheating phase,
i.e., the moment at which the inflaton energy density becomes equal to the SM
radiation energy density in the massless (massive) scenario. Note that the
$\rho_{\phi}(a_{\text{rh}})\\!=\\!\rho_{\rm SM}(a_{\text{rh}})$ equality
roughly coincides with $\langle\Gamma_{\phi\rightarrow hh}\rangle\\!\sim\\!H$
in both considered scenarios.

First of all, let us recall that in the both considered reheating scenarios
the inflaton energy density scales as $a^{-6n/(n+1)}$ (58) for
$a\\!<\\!a_{\text{rh}}$, while after the end of reheating $\rho_{\phi}$
sharply drops as it is seen in upper panels of Fig. 7 and Fig. 8. However, the
evolution of the SM radiation energy density during the reheating period is
drastically different in the massive reheating scenario as compared to the
case where the inflaton-induced Higgs mass has been neglected. This non-
trivial behavior not only changes the SM radiation energy density but also
affects the duration of reheating in the two cases. Note that when inflaton-
induced Higgs mass effects are taken into account, the SM radiation energy
density is suppressed compared to the massless case, which results in
elongation of the reheating period for the massive reheating case. Such
effects are more significant for relatively large inflaton–Higgs coupling
$g_{h\phi}$, which is manifestly shown in Fig. 7 and Fig. 8. These numerical
results for the SM radiation energy density agree well with our analytic
results Eq. (79) for massless and Eq. (85) for massive reheating cases.
Similarly, the SM bath temperature $T$, shown in the lower panel of these
plots, also agrees with our analytic results obtained above. Finally, we
should also emphasize that after the end of reheating, i.e., $a>a_{\rm rh}$,
the energy density of the Universe is mainly dominated by the SM radiation
energy density which scales as $a^{-4}$ in both considered scenarios.

As it is shown in the lower panel of these plots, the Hubble rate $H(a)$
scales as $a^{-3n/(n+1)}$ during the reheating phase and as $a^{-2}$
afterward, independent of the Higgs mass effects. However, as discussed above,
the most significant effect of the non-trivial Higgs mass during reheating is
on the inflaton decay rate to the Higgs boson due to phase space suppression.
In the scenario with the massless Higgs boson, the inflaton decay rate scales
as $a^{-3(n-1)/(n+1)}$ (77), which is constant for $n=1$. Whereas, in the
presence of the kinematical suppression the evolution of the inflaton decay
rate presented in Fig. 7 and Fig. 8 agrees well with analytic result (83).
Note, however, that the above analytical results are only valid during the
reheating period ($a_{e}\leq a\leq a_{\rm rh}$) when the inflaton energy
density (58) is analytically calculated. Since after the end of reheating
$a>a_{\rm rh}$ the inflaton energy sharply vanishes, therefore the inflaton
decay rate scale accordingly. In particular, as it has been discussed below
Eq. (65), the inflaton width $\langle\Gamma_{\phi\rightarrow hh}\rangle$
adopted in the RHS of Eq. (47) depends on $\rho_{\phi}$ through this inflaton
mass (38). Thus, for $n\\!>\\!1$, according to Eq. (66) the vanishing inflaton
energy density $\rho_{\phi}$ after the end of reheating implies divergent
averaged “width”, which is indeed seen in Fig. 7 and Fig. 8 for $n=3/2,2$. The
Higgs mass effects are clearly important for large inflaton-Higgs coupling
e.g., $g_{h\phi}=10^{-5}$, however as discussed in Sec. IV.2 even for
relatively small coupling e.g., $g_{h\phi}=10^{-10}$ for $n>3/2$ the
kinematical suppression can be significant. This effect is can be seen in Fig.
8 with $n=2$.

It should also be pointed out that when the kinematic suppression of the
inflaton decay width disappears at $a\\!\gtrless\\!a_{k}$ for
$n\\!\lessgtr\\!3/2$, the averaged inflaton decay rate
$\langle\Gamma_{\phi\rightarrow hh}\rangle$ for massless and massive cases
nearly merge. As mentioned above the convergence of the two cases is not
perfect as for the massless reheating case the inflaton “decays” to four
massless Higgs components. On the other hand, for the massive case (even if
the mass contribution is tiny), during one half of the oscillation period
(${\cal P}\\!<\\!0$ the electroweak symmetry broken phase), the decay final
state is just one real Higgs particle $h_{0}$, while during the other half
(${\cal P}\\!>\\!0$ unbroken phase) the final state is made of four massive
Higgs components $h_{i}$ degenerate in mass. Therefore, eventually, the
averaged inflaton decay rate $\langle\Gamma_{\phi\to hh}\rangle$ for the
massive case, in the limit when mass effects are negligible, is a factor 5/8
smaller than that of the massless case, see Fig. 7 and Fig. 8.

Figure 9: Left panel: Relation between reheating numbers of e-folds $N_{\rm
rh}$ and the value of the inflaton-Higgs coupling $g_{h\phi}$. Right panel:
Relation between the maximal temperature, $T_{\rm max}$, obtained during
reheating and the value of the inflaton-Higgs coupling $g_{h\phi}$. Filled
(empty) dotes present pairs of $(g_{h\phi},N_{\rm rh})$ obtained assuming the
massive (massless) reheating scenario, whereas solid (dashed) lines are
corresponding linear interpolations to the data. In the upper panel, colored
regions show constraints on $N_{\rm rh}$ for different values of $n$, see
Table (3).

We have seen that the consequences of the non-zero SM Higgs mass during the
reheating phase can be dramatic. Not only does it change the evolution of the
radiation energy density and the temperature, but it also affects the duration
of the reheating period, cf. the left panel of Fig. 9, where we have compared
the number of e-folds for reheating phase, $N_{\rm rh}$, as a function of the
inflaton-Higgs coupling, $g_{h\phi}$, in both the massless and massive cases.
In the massive case, reheating is not only delayed but also strongly
suppressed. In the scenario with the massive Higgs field, the production of
the SM radiation is always less efficient than in the massless case.
Therefore, the thermal bath temperature, measured by $\rho_{\rm SM}$, in the
massive scenario is reduced compared to the massless case during the reheating
phase, shown in the right-panel of Fig. 9 as a function of $g_{h\phi}$. It is
worth noting that in the standard scenario with massless Higgs, the maximal
temperature, $T_{\rm max}$, is typically attained shortly after the end of
inflation. However, in the massless case and for $n\\!=\\!2/3$ it is reached
after a few e-folds of reheating when kinematic suppression becomes irrelevant
and $\langle\Gamma_{\phi\rightarrow hh}\rangle$ approaches
$\langle\Gamma_{\phi\rightarrow hh}^{(0)}\rangle$. The relation between the
maximal temperature of the thermal bath and the value of the inflaton-Higgs
coupling in both reheating scenarios is shown in the lower panel in Fig. 9.
The largest discrepancy between the massless and massive cases is again
observed for the $n\\!=\\!2$ curve. As the strength of the inflaton-Higgs
interactions decreases, the role of the non-zero mass of the Higgs boson
becomes less and less relevant for the dynamics of reheating, and as we can
see in the left panel of Fig. 9, all solid lines slowly approach the
corresponding dashed lines. For the small value of the inflaton-Higgs
coupling, e.g., $g_{h\phi}=10^{-10}$, the distinction between solid and the
corresponding dashed lines is relatively small. However, the $n\\!=\\!2$ line
is an exception to this rule, and even for $g_{h\phi}=10^{-10}$ we observe a
non-negligible deviation from the massless case. This is caused by the fact
that for the $n\\!=\\!2$ case, the slope of $\rho_{\rm SM}$ as a function of
$a$ is steeper in the massive scenario, and thus it takes more e-folds of
reheating to drop $\rho_{\phi}$ below $\rho_{\rm SM}$.

##  V DM production during reheating

In this section, we study implications of Higgs-boson-induced reheating
discussed in the above section for the production of vector DM. As outlined in
Sec. II and specified in (6), the vector DM $X_{\mu}$ interacts directly with
the inflaton field and the SM Higgs bosons through dim-5
($\mathcal{C}_{X}^{\phi}m_{X}^{2}/M_{\rm Pl}\phi X_{\mu}X^{\mu}$) and dim-6
($\mathcal{C}_{X}^{{\bm{h}}}m_{X}^{2}/M_{\rm
Pl}^{2}|{\bm{h}}|^{2}X_{\mu}X^{\mu}$) operators, respectively. Moreover, the
vector DM couples indirectly to the energy-momentum tensor of the inflaton and
the SM through $s$-channel graviton exchange. Furthermore, in the massive
reheating scenario, the Higgs portal operator
$\mathcal{C}_{X}^{{\bm{h}}}m_{X}^{2}/M_{\rm
Pl}^{2}|{\bm{h}}|^{2}X_{\mu}X^{\mu}$ expanded around $v_{h}$ provides a term
$\propto v_{h}h_{0}X_{\mu}X^{\mu}$, which accounts not only for the Higgs
boson decays to DM pairs, but it also generates indirect DM-SM interactions
mediated by the Higgs exchange. As already noted in Sec. II, such processes
can only occur in the one-half of the inflaton oscillations period, when the
symmetry is broken and $v_{h}$ is non-zero. To sum up, ignoring the inflaton-
induced Higgs mass, i.e., in the massless reheating scenario, the dark sector
can be populated either due to gravitational interactions with the inflaton
and SM particles or a result of direct higher-dimensional interactions with
the $\phi$ and Higgs field. On the other hand, in the massive scenario, there
exist two additional DM production channels, i.e., direct decays of the SM
Higgs boson $h_{0}$ and freeze-in from SM particles via s-channel $h_{0}$
exchange.

The DM dynamics is governed by the Boltzmann equation (49), which can be
recast in the following form in terms of the comoving number density
$N_{X}\equiv n_{X}a^{3}$:

$\displaystyle\frac{dN_{X}}{da}$
$\displaystyle=\frac{a^{2}}{H}\Big{[}\mathcal{D}_{\phi}+\mathcal{S}_{\phi}+\mathcal{S}_{\rm{SM}}+\mathcal{D}_{h_{0}}\Big{]}.$
(100)

The first term on the r.h.s. of the above equation takes into account DM
production through inflaton decays. The source term $\mathcal{S}_{\phi}$
describes the gain of DM number density due to gravitational interactions with
the inflaton. In the scattering term, $\mathcal{S}_{\rm{SM}}$ (51), the
annihilation cross-section $\sigma_{XX\rightarrow{\rm SM\,SM}}$ includes
contributions from the graviton exchange and the Higgs portal operator.
Finally, the decay term $\mathcal{D}_{h_{0}}$, defined in Eq. (51), accounts
for Higgs boson decays into pairs of DM vectors.

All the terms on the r.h.s of (100) can be equally important for DM
production. However, their origin is very different, and thus it is convenient
to discuss each channel separately. In particular, the gravitational DM
production, through the graviton exchange with the inflaton background field
and the SM radiation bath, can be treated as an irreducible production
mechanism that is always present, regardless of other DM interactions. Let us
also emphasize that these two sources of gravitational production of DM are
very different despite their deceptive resemblance. In this work, we treat
$\phi$ as a homogenous, classical field, not a quantum particle. Within this
framework, DM vectors are produced from the vacuum in the oscillating
background of the inflaton field. To put it another way, the energy density of
the $\phi$ field is transferred to the dark sector indirectly through gravity
which couples to $T_{\phi}^{\mu\nu}$ and $T_{X}^{\mu\nu}$. Note that in this
case, DM particles are produced non-thermally since the inflaton field is not
in thermal equilibrium. Thus, this kind of DM production mechanism is
insensitive to the thermal history of the Universe and depends mainly on the
initial value of the inflaton energy density and its evolution during the
reheating period. Contrarily, the gravitational production from the SM
particles realizes a standard freeze-in DM scenario, in which SM quantum
states couple to $s$-channel virtual graviton which subsequently couples to
pairs of DM species. In this case, SM species are assumed to be in thermal
equilibrium, which implies a non-trivial dependence of the DM relic abundance
on the evolution of the SM temperature $T$.

In what follows, we assume that DM production does not have a significant
impact on the evolution of the first two Boltzmann equations, i.e., we adopt
two assumptions: (i) the inflaton decays mainly to the SM, and (ii) DM
particles are not in equilibrium with the SM thermal bath. Therefore, DM
production does not substantially affect the evolution of the SM bath
temperature and the Hubble rate. This requires a small DM branching ratio,
i.e., $\langle\Gamma_{\phi\rightarrow
XX}\rangle/\langle\Gamma_{\phi}\rangle\ll 1$, and sufficiently weak
interactions between the SM and DM. It is worth emphasizing that such
suppression emerges naturally in our model with vector DM since the
inflaton–DM as well as SM–DM interactions are sourced by the higher-
dimensional terms (6). Thus the above two assumptions seem to be well
justified and robust within our model. Therefore, in the following, we solve
the Boltzmann equation (100) for DM evolution with the inflaton and SM energy
densities as well as other related parameters obtained in the previous
section.

###  V.1 Inflaton induced gravitational DM production

Let us start our discussion with the most generic production mechanism, namely
the gravitational DM production in the background of the oscillating inflaton
field. This scenario does not require any additional (besides gravitational)
coupling of $X$ particles to the inflaton and thus can be treated as a kind of
irreducible “background” to any DM production mechanism, which is present
regardless of other DM interactions. In particular, the gravitational DM
production from the $\phi$ field provides an unavoidable contribution to the
relic DM abundance, which should be taken into account in any DM scenario.
Recently, it was shown [28] that in the case of scalar and fermionic DM, pure
gravitational production from the inflaton field can account for the observed
abundance of DM particles,

$\displaystyle\Omega_{X}^{\rm obs}h^{2}=0.1198\pm 0.0012\,,$ (101)

measured by the Planck Collaboration [40]. In this work, we focus on spin-1 DM
field, for which the source term $\mathcal{S}_{\phi}$ takes the following
form,

$\displaystyle\mathcal{S}_{\phi}$
$\displaystyle\\!=\\!\frac{\rho_{\phi}^{2}}{8\pi M_{\rm
Pl}^{4}}\sum_{k=1}^{\infty}|\mathcal{P}_{k}^{2n}|^{2}\bigg{(}\\!1-\frac{4m_{X}^{2}}{(k\omega)^{2}}+\frac{12m_{X}^{4}}{(k\omega)^{4}}\\!\bigg{)}\sqrt{1\\!-\\!\frac{4m_{X}^{2}}{(k\omega)^{2}}},$
(102)

where we have expressed $V(\phi)$ using the envelope $\varphi$ and the quasi-
periodic function $\mathcal{P}$ as
$V(\phi)\simeq\Lambda^{4}(\varphi/M)^{2n}\mathcal{P}^{2n}(t)$ and decomposed
$\mathcal{P}^{2n}(t)$ into the Fourier modes as

$\displaystyle\mathcal{P}^{2n}(t)=\sum_{k}\mathcal{P}_{k}^{2n}\,e^{-ik\omega
t},$ (103)

with

$\displaystyle\mathcal{P}_{k}^{2n}=\frac{1}{\mathcal{T}(t_{0})}\int_{t_{0}}^{t_{0}+\mathcal{T}(t_{0})}dt\mathcal{P}^{2n}(t)\,e^{ik\omega
t},$ (104)

and the inflaton mode energy/frequency $\omega=2\pi/{\cal T}$. For $n\neq 1$,
the inflaton mode energy/frequency $\omega$ is time-dependant, where this
dependency is the same as the inflaton effective mass time-dependance (39).
Therefore, it is instructive to write the inflaton mode energy/frequency
$\omega$ as

$\displaystyle\omega(a)=\omega_{e}\Big{(}\frac{a_{e}}{a}\Big{)}^{\frac{3(n-1)}{n+1}}\,,$
(105)

where $\omega_{e}\\!\equiv\\!\omega(a_{e})$ is the mode frequency at the onset
of reheating phase. The values of the $\mathcal{P}_{k}^{2n}$ coefficients
decrease with $k$ for all considered values of $n$. Consequently, the
$\sum_{k}|\mathcal{P}_{k}^{2n}|^{2}$ sum quickly converges. Moreover, for
$n\\!=\\!1$ the only non-zero Fourier coefficient is
$\mathcal{P}_{2}^{2}=1/4$, whereas for $n\\!=\\!3/2$ ($n\\!=\\!2$) all even
(odd) coefficients are zero. The numerical values of the
$\sum_{k}|\mathcal{P}_{k}^{2n}|^{2}$ sum for different values of $n$ are
collected in Table 1. For details see Appendix B where ${\cal
S}_{\phi}\equiv{\cal D}^{(2)}_{\phi\to XX}$.

For the quadratic inflaton potential, i.e., $n\\!=\\!1$, the frequency
$\omega=m_{\phi}=\sqrt{2}\Lambda^{2}/M_{\rm Pl}$ is time-independent during
the reheating phase and the only non-zero Fourier coefficient for $|{\cal
P}_{k}^{2n}|$ is for $k=2$. In this case, the source term Eq. (102) simplifies
as

$\displaystyle{\cal
S}_{\phi}\overset{n=1}{=}\frac{1}{128\pi}\frac{\rho_{\phi}^{2}}{M_{\rm
Pl}^{4}}\bigg{(}\\!1-\frac{m_{X}^{2}}{m_{\phi}^{2}}+\frac{3m_{X}^{4}}{m_{\phi}^{4}}\\!\bigg{)}\sqrt{1\\!-\\!\frac{m_{X}^{2}}{m_{\phi}^{2}}}.$
(106)

Note that the above expression receives time-dependence only through
$\rho_{\phi}$. Thus, the rate of DM production is expected to be the largest
at the onset of reheating. Similar behavior is observed for other values of
$n$ for light DM (LDM) such that $m_{X}\\!<\\!\omega_{e}$. However, in generic
case, i.e., $n\\!\neq\\!1$, the $S_{\phi}$ term has another source of time-
dependence — the frequency $\omega$, which increases (decreases) with time for
$n\\!<\\!1$ ($n\\!>\\!1$). For instance, for the $n\\!=\\!1$ case, the time-
independence of $\omega$ implies that there exists a constant mass threshold
during the whole period of reheating. Due to the fact that in this case, the
only non-zero Fourier coefficient is $\mathcal{P}_{2}^{2}$, the phase space
factor $\sqrt{1-(m_{X}/m_{\phi})^{2}}$ requires $m_{X}\\!<\\!m_{\phi}$. This
means that DM particles with masses exceeding the effective mass of the
inflaton field cannot be gravitationally produced from the inflaton
background. For the remaining values of $n$, we have infinitely many non-zero
coefficients $\mathcal{P}_{k}^{2n}$ that contribute to the source term
$\mathcal{S}_{\phi}$. In this case, we can compensate for the smallness of the
frequency $\omega$ considering higher harmonics. By increasing the value of
the denominator in the $\sqrt{1-(m_{X}/k\omega)^{2}}$ factor, we circumvent
the standard kinematical suppression. However, in this case,
$\mathcal{S}_{\phi}$ receives another suppression that comes from the Fourier
coefficients, since $\mathcal{P}_{k}^{2n}$ rapidly decreases with increasing
$k$. Consequently, the main contribution to the $\mathcal{S}_{\phi}$ term
comes from the first (minimal) non-zero harmonic $k_{\min}$, which for $n=2/3$
and $n=3/2$ is $k_{\min}=1$, whereas for $n=1$ and $n=2$ is $k_{\min}=2$.

We should also stress that for $n\\!>\\!1$ and for a given DM mass, the
$m_{X}/\omega$ ratio is the smallest at the onset of reheating and increases
with time. Thus, the DM production rate is the largest just after the end of
inflation. Contrarily, for the $n<1$ case, e.g., $n\\!=\\!2/3$, the
$m_{X}/\omega$ ratio decreases with time. Thus, in this case for heavy DM
(HDM), i.e., $\omega_{e}<2m_{X}$, the source term $\mathcal{S}_{\phi}$ is the
largest at some later moment $a_{X}\\!>\\!a_{e}$ when $\omega(a_{X})\sim
2m_{X}$. This behavior is manifestly shown in Fig. 10, where we have plotted
the evolution of the $a^{2}\mathcal{S}_{\phi}/H^{2}$ terms as a function of
the scale factor for LDM (orange curve) and HDM (purple and green curves)
cases. For $n\\!<\\!1$, the inflaton mode energy/frequency $\omega$ increases
with time, see Eq. (105). Thus, even if the DM mass is larger than the mode
frequency, i.e., $\omega<2m_{X}$ at the onset of reheating phase, the gradual
increase in the inflaton mode energy $\omega$ can provide the phase for the
production of DM before the end of reheating. This, in particular, means that
the production of HDM is kinematically disfavoured at the onset of reheating,
but if $2m_{X}\\!<\\!\omega(a_{\rm{rh}})$ such particles can be nevertheless
produced abundantly during the reheating period. It is useful to write
explicitly the form of $a_{X}$, i.e., $k_{\min}\omega(a_{X})=2m_{X}$, using
Eq. (105) as

$\displaystyle a_{X}$
$\displaystyle=a_{e}\begin{dcases}\Big{\lfloor}\frac{2m_{X}}{\omega_{e}}\Big{\rfloor}^{\frac{1+n}{3(1-n)}},&n<1,\\\
1,&n\geq 1.\end{dcases}$ (107)

Above ‘floor’ function $\lfloor\cdots\rfloor$ is defined as follows:

$\lfloor x\rfloor=\begin{cases}1,&x\leq 1,\\\ x,&x>1.\end{cases}$ (108)

In other words, $a_{X}$ indicates the lowest value of the scale factor for
which the lowest non-zero harmonics $k_{\min}$ with a given $n$ can produce DM
species with mass $m_{X}$.

Figure 10: The evolution of the comoving number density
$dN_{X}/da=a^{2}\mathcal{S}_{\phi}/H$ as a function of the scale factor $a$
for $n\\!=\\!2/3$. For LDM, i.e., $2m_{X}\\!<\\!\omega(a_{e})$ (orange curve),
the source term $\mathcal{S}_{\phi}$ reaches the maximal value at the
beginning of reheating, i.e., $a=a_{e}$. Whereas, for HDM, i.e.,
$2m_{X}\\!>\\!\omega(a_{e})$ (purple and green curves), the source term
$\mathcal{S}_{\phi}$ is not maximized at $a=a_{e}$ due to the kinematical
suppression, but instead grows with time until $\omega(a)\sim m_{X}$.

The number density of DM species produced gravitationally from the inflaton
background can be obtained by solving the Boltzmann equation (100). Keeping
only the gravitationally produced DM from the inflaton field, i.e.,
$\mathcal{S}_{\phi}$ term, we get,

$\displaystyle
n_{X}(a)=\frac{1}{a^{3}}\int_{a_{e}}^{a}d\tilde{a}\,\frac{\tilde{a}^{2}}{H(\tilde{a})}\mathcal{S}_{\phi}(\tilde{a}).$
(109)

The comoving number density, $N_{X}$, approaches a constant limit around
$a_{f}\simeq a_{\rm{rh}}$, which results from the fact that after the end of
reheating, the inflaton is depleted and the contribution of $\rho_{\phi}$ to
the total energy density quickly becomes negligible. Thus, the present-day
number density of DM particles, $n_{X}(a_{0})$, can be well approximated by
the solution of (109) with the upper limit of the integral taken as
$a=a_{\rm{rh}}$. Inserting solutions for $H(a)$ and $\rho_{\phi}(a)$, obtained
in the previous section, we find

$\displaystyle n_{X}(a_{0})$
$\displaystyle\simeq\frac{n+1}{6n-3}\frac{\sqrt{3}}{8\pi}\bigg{(}\frac{\sqrt{\rho_{\phi_{e}}}}{M_{\rm
Pl}}\bigg{)}^{\\!3}\,\sum_{k=1}^{\infty}|\mathcal{P}_{k}^{2n}|^{2}$
$\displaystyle\qquad\times\Big{(}\frac{a_{e}}{a_{X}}\Big{)}^{\\!\frac{6n-3}{n+1}}\Big{(}\frac{a_{e}}{a_{0}}\Big{)}^{\\!3}\bigg{[}1-\Big{(}\frac{a_{X}}{a_{\rm{rh}}}\Big{)}^{\frac{6n-3}{n+1}}\bigg{]}.$
(110)

In the above equation, we have neglected the phase space factor proportional
to $m_{X}^{2}/(k\omega)^{2}$. Moreover, we have also assumed that
$N_{X}(a_{\rm rh})\simeq N_{X}(a_{0})$, which reflects the fact that the
efficient DM gravitational production in the inflaton background is possible
only during reheating.

Now let us discuss the solution (110) in more detail. Firstly, we should
emphasize that for $n\\!>\\!1/2$, as employed in this work, the second factor
in the squared bracket is always small, i.e.,
$(a_{X}/a_{\rm{rh}})^{6n-3/(n+1)}\\!\ll\\!1$ since $a_{X}\leq a_{\rm rh}$,
therefore this factor can neglected. Secondly, for $n\\!\geq\\!1$ we should
always take $a_{X}\\!=\\!a_{e}$, since in this case the inflaton mode energy
$\omega$ is the largest at $a=a_{e}$, which implies that for $a>a_{e}$ the
phase space suppression increases. Finally, for $n\\!<\\!1$, $a_{X}$ can be
larger or equal $a_{e}$ depending on the relation between $2m_{X}$ and
$\omega_{e}$. In particular, the for heavy DM, i.e., $2m_{X}>\omega$, the
$a_{X}$ is related to the DM mass through the following relation

$\displaystyle a_{X}$
$\displaystyle=a_{e}\Big{(}\frac{2m_{X}}{\omega_{e}}\Big{)}^{\frac{1+n}{3(1-n)}},$
$\displaystyle n$ $\displaystyle<1.$ (111)

Taking into account all simplifications we can recast Eq. (110) as

$\displaystyle n_{X}(a_{0})$
$\displaystyle\simeq\frac{\sqrt{3}}{8\pi}\bigg{(}\frac{\sqrt{\rho_{\phi_{e}}}}{M_{\rm
Pl}}\bigg{)}^{3}\Big{(}\frac{a_{e}}{a_{0}}\Big{)}^{3}\frac{n+1}{6n-3}\sum_{k}|\mathcal{P}_{k}^{2n}|^{2}$
$\displaystyle\qquad\times\begin{dcases}\Big{\lceil}\frac{\omega_{e}}{2m_{X}}\Big{\rceil}^{\\!\frac{2n-1}{1-n}},&n\\!<\\!1,\\\
1,&n\\!\geq\\!1.\end{dcases}$ (112)

In any successful model, the predicted value of a DM abundance today,

$\displaystyle\Omega_{X}=\frac{m_{X}n_{X}(a_{0})}{\rho_{c}},$ (113)

should match with the observed value (101). Above, $\rho_{c}=1.054\times
10^{-5}h^{2}\;{\rm GeV}\,{\rm cm}^{-3}$ denotes the critical density. Using
Eq. (112) we can estimate the present-day amount of DM produced purely
gravitationally from the inflaton field as

$\displaystyle\Omega_{X}^{\rm{gr},\phi}=\frac{\sqrt{3}}{8\pi}\bigg{(}\frac{\sqrt{\rho_{\phi_{e}}}}{M_{\rm
Pl}}\bigg{)}^{3}\Big{(}\frac{a_{e}}{a_{\rm{rh}}}\Big{)}^{3}\frac{s_{0}}{s_{\rm
rh}}\frac{1}{\rho_{c}}$
$\displaystyle\quad\times\frac{n+1}{6n-3}\sum_{k}|\mathcal{P}_{k}^{2n}|^{2}\begin{dcases}m_{X}\Big{\lceil}\frac{\omega_{e}}{2m_{X}}\Big{\rceil}^{\\!\frac{2n-1}{1-n}},&n\\!<\\!1,\\\
m_{X},&n\\!\geq\\!1.,\end{dcases}$ (114)

where $s$ is the entropy-density, and $s_{0}\\!=\\!2970\,{\rm cm}^{-3}$
denotes its present-day value and $s_{\rm rh}\\!\equiv\\!s(a_{\rm rh})$ is the
entropy density at the end of reheating. In the above relation we have used
the conservation of the entropy, i.e., $s_{\rm
rh}/s(a_{0})\simeq\frac{a_{0}^{3}}{a_{\rm{rh}}^{3}}$, neglecting the change of
$g_{\star,S}$ in the period between $a_{\rm{rh}}$ and $a_{0}$.

For $n\\!=\\!2/3$ and heavy DM, i.e., $2m_{X}>\omega_{e}$ we note a very
peculiar behaviour of $\Omega_{X}^{\rm{gr},\phi}h^{2}$. It turns out that in
this case, the predicted abundance of DM particles does not depend on $m_{X}$,
as the ‘ceiling’ function in Eq. (114) results in $\omega_{e}/(2m_{X})$ and
the $m_{X}$ factors cancel out. This mass-independence of
$\Omega_{X}^{\rm{gr},\phi}h^{2}$ for heavy DM is observed only if reheating
lasts long enough such that $\omega(a)$ equals $2m_{X}$ well before the end of
reheating. Thus, the obtained approximation works only for DM particles with
mass $\omega_{e}\\!\lesssim\\!2m_{X}\\!\lesssim\\!\omega(a_{\rm{rh}})$. Note
also that in other cases, for fixed value of $\Lambda$, $\alpha$, the present-
day abundance of DM species depends only on the value of $n$, $m_{X}$ and thus
on $a_{i}$ and $a_{\rm{rh}}$. Therefore, it is rather not sensitive to the
evolution of the thermal bath but depends on the duration of the reheating
period.

###  V.2 Gravitational freeze-in production

Next, we discuss the second irreducible DM production mechanism, namely the
gravitational freeze-in from the SM particles. Gravitational freeze-in is a
universal phenomenon in which the SM particles annihilate to produce DM
through the s-channel graviton exchange. This mechanism only requires
gravitational interaction between the SM and DM. Similar to the mechanism
discussed in the previous subsection, the gravitational freeze-in from SM is
most efficient during the reheating period. However, in this case, DM species
are produced in a quantum process from two SM particles in thermal
equilibrium. Thus, the abundance of DM particles produced from SM quanta is
sensitive to the thermal history of the Universe. Moreover, since the
amplitude of such a process is proportional to $s/M_{\rm Pl}^{2}$, where $s$
is the Mandelstam variable, efficient DM production requires high energy or
temperature of the SM annihilating particles. In principle, it was shown,
e.g., [34, 33], that pure gravitational DM production from the SM sector can
account for all the observed DM relic density (101) if the reheating
efficiency and the thermal bath temperature are sufficiently high. In contrast
to the production in the inflaton background, the freeze-in DM production
depends on the evolution of the radiation sector, or more precisely, on the
evolution of the thermal bath temperature during the reheating phase. Thus,
the efficiency of DM gravitational freeze-in production depends on the
dynamics of reheating, which is determined by the scale of inflaton $\Lambda$,
the inflaton–Higgs coupling $g_{h\phi}$, $n$ and the inflaton-induced mass of
the Higgs boson.

For the gravitational freeze-in DM production, the Boltzmann equation
simplifies as

$\displaystyle\frac{dN_{X}}{da}=\frac{a^{2}}{H}\mathcal{S}_{\rm{SM}}=\frac{a^{2}}{H}\bar{n}_{X}^{2}\langle\sigma\lvert
v\rvert\rangle_{XX\rightarrow h_{\mu\nu}\rightarrow\rm{SM\,SM}},$ (115)

where we have kept only the $\mathcal{S}_{\rm{SM}}$ term on the r.h.s. of
(100). The thermally averaged cross-section $\langle\sigma\lvert
v\rvert\rangle$ for any $2\to 2$ annihilation process can be calculated using
the standard formula [51],

$\displaystyle\langle\sigma\lvert v\rvert\rangle_{ii\to ff}$
$\displaystyle=\frac{(2J_{i}+1)^{2}}{\bar{n}_{i}^{2}}\frac{T}{32\pi^{4}}\int_{4m_{\rm
max}^{2}}^{\infty}ds\;\sqrt{s}(s-4m_{i}^{2})\,\sigma_{ii\to
ff}(s)K_{1}\left(\frac{\sqrt{s}}{T}\right),$ (116)

where $m_{\rm max}\equiv{\rm max}(m_{i},m_{f})$, and $\sigma_{ii\to ff}$
denotes the cross-section averaged over possible spin states of the DM
particles $i$ in the initial state and SM particles $f$ in the final state.
For the gravitational freeze-in production of DM, there are four relevant
annihilation processes,

$\displaystyle\sigma_{XX\rightarrow h_{\mu\nu}\rightarrow h_{i}h_{i}}$
$\displaystyle=\frac{1}{9}\frac{1}{2}\frac{1}{480\pi s^{3}M_{\rm
Pl}^{4}}\frac{\sqrt{s-4\bar{m}_{h_{i}}^{2}}}{\sqrt{s-4m_{X}^{2}}}\bigg{[}2\bar{m}_{h_{i}}^{4}\left(84m_{X}^{4}+68m_{X}^{2}s+19s^{2}\right)$
$\displaystyle\hskip
142.26378pt-2\bar{m}_{h_{i}}^{2}s(6m_{X}^{2}+s)(2m_{X}^{2}+7s)+3s^{2}(6m_{X}^{4}+2m_{X}^{2}s+s^{2})\bigg{]},$
(117) $\displaystyle\sigma_{XX\rightarrow
h_{\mu\nu}\rightarrow\bar{\psi}\psi}$
$\displaystyle=\frac{1}{9}\frac{1}{960\pi s^{2}M_{\rm
Pl}^{4}}\frac{\sqrt{s-4\bar{m}_{\psi}^{2}}}{\sqrt{s-4m_{X}^{2}}}\\!\bigg{(}\\!1\\!-\\!\frac{4\bar{m}_{\psi}^{2}}{s}\bigg{)}\\!\bigg{[}\\!2\bar{m}_{\psi}^{2}\left(84m_{X}^{4}+68m_{X}^{2}s+19s^{2}\right)\\!+\\!s(48m_{X}^{4}+56m_{X}^{2}s+13s^{2})\\!\bigg{]},$
(118) $\displaystyle\sigma_{XX\rightarrow h_{\mu\nu}\rightarrow VV}$
$\displaystyle=\frac{1}{9}\frac{1}{\delta_{V}}\frac{1}{480\pi s^{3}M_{\rm
pl}^{4}}\frac{\sqrt{s-4\bar{m}_{V}^{2}}}{\sqrt{s-4m_{X}^{2}}}\bigg{[}\bar{m}_{V}^{4}(504m_{X}^{4}+408m_{X}^{2}s+114s^{2})+\bar{m}_{V}^{2}(408m_{X}^{4}s+536m_{X}^{2}s^{2}+118s^{3})$
$\displaystyle\hskip
142.26378pt+114m_{X}^{4}s^{2}+118m_{X}^{2}s^{3}+29s^{4}\bigg{]},$ (119)
$\displaystyle\sigma_{XX\rightarrow h_{\mu\nu}\rightarrow\gamma\gamma}$
$\displaystyle=\frac{1}{9}\frac{1}{2}\frac{1}{240\pi M_{\rm
Pl}^{4}}\frac{48m_{X}^{4}+56m_{X}^{2}s+13s^{2}}{\sqrt{s(s-4m_{X}^{2})}},$
(120)

which account for gravitational annihilation of the $X$ species into SM Higgs
bosons, fermions, and massive and massless vectors, respectively. Above
$\delta_{V}\\!=\\!2$ for the identical gauge bosons in the final states,
otherwise $\delta_{V}=1$. In the above annihilation cross sections, we have
taken into account the non-trivial SM particle masses due to inflaton-induced
electroweak vev $v_{h}$, where $\bar{m}$ denotes the time-averaged (over fast
inflaton oscillations) mass. In the limit $m_{\rm{SM}}\rightarrow 0$, the
corresponding thermally-averaged cross-sections are [34],

$\displaystyle\langle\sigma|v|\rangle_{XX\rightarrow h_{\mu\nu}\rightarrow
h_{i}h_{i}}$ $\displaystyle\\!=\\!\frac{9}{32\pi^{4}}\frac{1}{120\pi M_{\rm
Pl}^{4}}\frac{m_{X}^{6}T^{2}}{\bar{n}_{X}^{2}}\bigg{\\{}3K_{1}^{2}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}+2K_{2}^{2}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}\bigg{[}1+4\Big{(}\frac{T}{m_{X}}\Big{)}^{2}\bigg{]}+4\frac{T}{m_{X}}K_{1}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}K_{2}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}\bigg{\\}},$
(121) $\displaystyle\langle\sigma|v|\rangle_{XX\rightarrow
h_{\mu\nu}\rightarrow\bar{\psi}\psi}$
$\displaystyle\\!=\\!\frac{9}{32\pi^{4}}\frac{1}{90\pi M_{\rm
Pl}^{4}}\frac{m_{X}^{6}T^{2}}{\bar{n}_{X}^{2}}\bigg{\\{}11K_{1}^{2}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}+K_{2}^{2}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}\bigg{[}9+26\Big{(}\frac{T}{m_{X}}\Big{)}^{2}\bigg{]}+13\frac{T}{m_{X}}K_{1}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}K_{2}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}\bigg{\\}},$
(122) $\displaystyle\langle\sigma|v|\rangle_{XX\rightarrow
h_{\mu\nu}\rightarrow\gamma\gamma}$
$\displaystyle\\!=\\!\frac{9}{32\pi^{4}}\frac{3}{135\pi M_{\rm
Pl}^{4}}\frac{m_{X}^{6}T^{2}}{\bar{n}_{X}^{2}}\bigg{\\{}11K_{1}^{2}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}+K_{2}^{2}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}\bigg{[}9+26\Big{(}\frac{T}{m_{X}}\Big{)}^{2}\bigg{]}+13\frac{T}{m_{X}}K_{1}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}K_{2}\\!\Big{(}\\!\frac{m_{X}}{T}\\!\Big{)}\bigg{\\}}.$
(123)

On the other hand, in the massive reheating scenario, i.e., for
$m_{\rm{SM}}\\!\neq\\!0$, the integral over $s$ in Eq. (116) cannot be
obtained analytically because of the non-zero mass of the SM species,
generated by their coupling to the inflaton or the Higgs doublet. In this case
we integrate (116) numerically using Eqs. (117)–(119) and adopting the
solution for $T(a)$ which we have found in the previous section. It is
important to notice that even in the massive reheating scenario, SM photons
and gluons are still massless during reheating. Thus, one expects that in this
case the freeze-in DM production is dominated by the annihilation of the
massless SM gauge bosons, since other annihilation channels are Boltzmann
suppressed if $m_{\rm{SM}}/T\\!>\\!1$. In the limit of asymptotically small or
large DM mass the source term, $\mathcal{S}_{XX\rightarrow
h_{\mu\nu}\rightarrow\rm{SM}\,\rm{SM}}=\bar{n}_{X}^{2}\langle\sigma|v|\rangle_{XX\rightarrow
h_{\mu\nu}\rightarrow\rm{SM}\,\rm{SM}}$, for massless SM particles behaves as

$\displaystyle\mathcal{S}_{XX\rightarrow
h_{\mu\nu}\rightarrow\rm{SM}\,\rm{SM}}\\!\simeq\\!\begin{dcases}\\!\frac{A_{i}}{40\pi^{5}}\frac{T^{8}}{M_{\rm
Pl}^{4}}\,,&m_{X}\\!\ll\\!T,\\\
\\!\frac{B_{i}}{32\pi^{4}}\frac{m_{X}^{5}T^{3}}{M_{\rm
Pl}^{4}}e^{-2m_{X}/T}\,,&m_{X}\\!\gg\\!T,\end{dcases}$ (124)

where $A_{i}\\!=\\!3,13,26$ and $B_{i}\\!=\\!3/16,1,2$, for SM scalars,
fermions and vectors, i.e., $i=0,1/2,1$, respectively. The asymptotic behavior