UMN-TH-4129/22, FTPI-MINN-22/20

August 2022

# Post-Inflationary Dark Matter Bremsstrahlung

Yann Mambrini  Keith A. Olive  Jiaming Zheng

###### Abstract

Dark matter may only interact with the visible sector efficiently at energy
scales above the inflaton mass, such as the Planck scale or the grand
unification scale. In such a scenario, the dark matter is mainly produced out
of equilibrium during the period of reheating, often referred to as UV freeze-
in. We evaluate the abundance of the dark matter generated from bremsstrahlung
off the inflaton decay products assuming no direct coupling between the
inflaton and the dark matter. This process generally dominates the production
of dark matter for low reheating temperatures where the production through the
annihilations of particle in the thermal plasma becomes inefficient. We find
that the bremsstrahlung process dominates for reheating temperatures
$T_{\mathrm{RH}}\lesssim 10^{10}$ GeV, and produces the requisite density of
dark matter for a UV scale $\simeq 10^{16}$ GeV. As examples, we calculate
numerically the yield of the dark matter bremsstrahlung through gravitation
and dimension-6 vector portal effective interactions.

##  1 Introduction

Dark Matter (DM) is only indirectly observed by its gravitational effects and
its nature remains a mystery. A wide variety of well-motivated extensions to
the Standard Model contain stable fundamental particles as candidates for DM
and may link it to other profound physics questions such as the strong CP
problem, the hierarchy problem, charge quantization and gauge coupling
unification. Direct searches [1, 2, 3, 4, 5, 6] for particle DM have excluded
a large portion of the parameter space of the weakly interactive massive
particle (WIMP) scenario [7] where the DM was in chemical equilibrium with the
Standard Model (SM) sector in the early universe (see [8, 9] for reviews on
the subject).

Dark matter can also be produced out of chemical equilibrium [10, 11, 12, 13,
14, 15, 16, 17, 18], dubbed "freeze-in" production, if it only directly
interacts with particles in the thermal plasma through a very weak
interaction, even gravity [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
32, 33]. This scenario can also appear in grand unified theories (GUTs) where
the DM particles only interact with the visible sector through a heavy
mediator at the GUT scale [15, 34, 35, 36] or an intermediate scale associated
with moduli fields [37], the seesaw mechanism [38, 39] or even the inflaton
sector [40, 41]. In these scenarios, the dark matter is most efficiently
generated at higher temperatures because these interactions are effectively
irrelevant operators at low energy. This kind of DM production is usually
called “UV freeze-in” in the literature.

The DM - plasma particle interactions directly induce the annihilation of
plasma particles $\cal P$ into DM particles $\cal X$ as illustrated in Fig.
1(a) and has been widely considered in the literature. In a generic
inflationary cosmology, the plasma particle itself is produced by the inflaton
during reheating. The ${\cal P}{\cal P}-{\cal X}{\cal X}$ interaction
unavoidably induces radiative dark matter production processes during inflaton
decay or scattering, even in the absence of a direct coupling of the inflaton
to the dark matter sector as we show in Fig. 1. For example, the authors of
[42] calculated the decay of the inflaton to dark matter through loop
processes as in Fig. 1(c) and found that it can be sufficient to generate the
right relic abundance. However, the loop process is not always operative for
every form of DM-plasma particle interaction. In particular, the loop process
vanishes if the DM-plasma interaction is mediated by a heavy vector boson or
by a spin-$2$ particle such as the graviton. However, in addition to the loop
process, there is the possibility of a four-body final state contribution to
inflaton decay as shown in Fig. 1(b), which is a bremsstrahlung-like process.
As long as the interaction in (a) is present, the four-body decay will have a
non-zero branching ratio and should be taken into account in the total
production of dark matter.

In this paper, we study the bremsstrahlung production of dark matter during
inflaton decay as in Fig. 1(b). As we will show, this process tends to
dominate the production of DM in a vast parameter space, especially with
relatively low reheating temperature. We will then focus on the comparison
between Fig. 1(a) and Fig. 1(b), and provide numerical results for the DM
bremsstrahlung rate in two UV scenarios where the loop process in Fig. 1(c)
vanishes, i.e., 1) when the DM only interacts gravitationally, and 2) when the
DM interacts through a heavy vector boson.

Figure 1: DM production processes. The interaction (black bulb) that induces
the annihilation of a pair of plasma particle ${\cal P}$ to a pair of dark
matter ${\cal X}$ in panel (a) unavoidably induces DM bremsstrahlung from
inflaton decay in panel (b). Depend on the form of the interaction, the loop
process in panel (c) may also be induced.

The organization of this work is as follows. In Section 2, we provide general
qualitative estimates of the DM bremsstrahlung. In Section 3 we consider the
minimal (unavoidable) extension, where the only interaction between the
visible and dark sector is through gravity. Finally, we analyze in Section 4
the DM bremsstrahlung production when the DM-plasma interaction is induced by
heavy mediators, and comment briefly on possible UV completions. We conclude
in Section 5.

##  2 Bremsstrahlung production versus plasma annihilation: general
considerations

We start with a general estimate of the DM bremsstrahlung during the reheating
process and compare it with the production of DM through the scattering of
particles in the plasma. We consider a minimal scenario where the reheating is
induced by the decay of the inflaton $\phi$ into two particles ${\cal P}$ that
rapidly thermalizes in a hot plasma111We discuss the production of dark matter
during the thermalization phase in the Appendix.. The plasma particles ${\cal
P}$ can interact with the dark matter particles ${\cal X}$ through an
effective interaction with energy scale $\Lambda$. We restrict our discussion
to $\Lambda\gg m_{\phi}$, so that the production rate is not sensitive to the
production threshold of the mediator. Note that $\cal P$ and $\cal X$ are
symbols that represent general plasma (SM) and dark matter particles
respectively disregarding their detailed properties. The ${\cal P}-{\cal X}$
interaction induces an out-of-equilibrium annihilation process ${\cal
PP}\rightarrow{\cal XX}$, as shown in Fig. 1(a). The same effective
interaction also induces bremsstrahlung production of an ${\cal X}$ pair along
with ${\cal P}$ final states during the reheating process through the inflaton
decay $\phi\rightarrow{\cal P}{\cal P}\rightarrow{\cal PPXX}$, as shown in
Fig. 1(b). We will refer to this process as the bremsstrahlung production
throughout our work. It turns out that the bremsstrahlung production of dark
matter generally dominates when the inflaton decay rate is slow relative to
the expansion rate and results in a low reheating temperature, as we will now
show.

For concreteness, we consider a simplified scenario of reheating where the
real scalar inflaton $\phi$ mainly decays to a pair of complex scalars $s$ or
a pair of chiral fermions $f$ during reheating,222The chirality of $f$
suppresses $\phi\rightarrow f\bar{f}$.

$\phi\rightarrow s^{*}s,ff~{}({\rm or}~{}\bar{f}\bar{f})\,.$ (2.1)

From CPT, the inflaton decay rate to the $ff$ final state should be equal to
the $\bar{f}\bar{f}$ final state. For the brevity of notation, we will simply
write $ff$ to include both channels in the rest of this paper. Examples of
inflaton decay channels are a pair of SM Higgs particles or a pair of right-
handed neutrinos which can be easily realized by relevant operators $\phi
H^{\dagger}H$ or $(\phi\nu_{R}^{T}C\nu_{R}+h.c.)$ for a singlet inflaton
$\phi$, where $C$ is the charge conjugate operator. We assume $m_{s},m_{f}\ll
m_{\phi}$ so that the efficiency of the reheating is not blocked by phase
space factors, and the inflaton decay products thermalize the Standard Model
plasma instantaneously [43, 44, 45, 46] (see however [47, 48, 49, 50]).

The interactions responsible for the inflaton decays are described by the
operators,

$\displaystyle{\cal L}_{\phi s}$ $\displaystyle=\mu_{\phi s}|s|{}^{2}\phi\,,$
(2.2a) $\displaystyle{\cal L}_{\phi f}$ $\displaystyle=\frac{c_{\phi
f}}{2}f^{T}Cf\phi+h.c.~{}\,.~{}$ (2.2b)

The corresponding inflaton decay rate to a final state $i$ is

$\Gamma_{\phi}^{i}=\frac{\lambda_{i}^{2}}{16\pi}m_{\phi}\,,$ (2.3)

with the effective coupling $\lambda_{i}$ defined as

$\lambda_{s}=\frac{\mu_{\phi s}}{m_{\phi}}\,,\qquad\lambda_{f}=c_{\phi f}\,,$
(2.4)

for $\phi\rightarrow s^{*}s$ and $\phi\rightarrow ff$, respectively. We would
expect the natural range of $\lambda_{s}$ from a UV model to be comparable to
that of $\lambda_{f}$.

The evolution of the inflaton, the radiation(thermal plasma) and the dark
matter after the end of inflation (defined when $\ddot{a}=0$, where $a$ is the
cosmological scale factor) is governed by the system of the Friedmann and
Boltzmann equations [51, 52, 53],

$\displaystyle\frac{{\rm d}\rho_{\phi}}{{\rm d}t}+3H\rho_{\phi}$
$\displaystyle=-\Gamma_{\phi}\rho_{\phi}\,,$ (2.5a) $\displaystyle\frac{{\rm
d}\rho_{R}}{{\rm d}t}+4H\rho_{R}$ $\displaystyle=\Gamma_{\phi}\rho_{\phi}\,,$
(2.5b) $\displaystyle\frac{{\rm d}n_{{\cal X}}}{{\rm d}t}+3Hn_{{\cal X}}$
$\displaystyle=R_{{\cal X}}\,,$ (2.5c)

where $R_{{\cal X}}=R_{\phi{\cal X}}+R_{{\cal PX}}$ is the total production
rate of the dark matter that includes the contributions from the annihilation
of plasma particles $R_{{\cal PX}}$ and from the DM bremsstrahlung
$R_{\phi{\cal X}}$. Note that we have assumed a quadratic form of the inflaton
potential during the oscillation phase in Eq.(2.5a). The equation of state and
thus (2.5a) should be modified for more complex potentials with a minimum
approximated by $V\sim\phi^{k}$ [54, 55].

After inflation, the universe is dominated by an oscillating inflaton field
whose energy density dilutes with the scale factor approximately as $a^{-3}$.
During this epoch, the radiation temperature from the inflaton decay scales as
$T\propto a^{-3/8}$ [51, 52, 53, 54, 55, 56, 57, 58]. If we define ${\cal
Y}_{{\cal X}}\equiv n_{{\cal X}}/T^{8}$, the dark matter density during
inflaton domination evolves as (see, e.g. [28]),

$\frac{{\rm d}{\cal Y}_{\cal X}}{dT}=-\frac{8}{3}\frac{R_{\cal X}}{HT^{9}}\,,$
(2.6)

and

$H^{2}(T)=\frac{\rho_{\phi}(T)}{3M_{P}^{2}}=\frac{\rho_{\phi}({T_{\mathrm{RH}}})}{3M_{P}^{2}}\left(\frac{T}{T_{\mathrm{RH}}}\right)^{8}=\frac{\alpha
T_{\mathrm{RH}}^{4}}{3M_{P}^{2}}\left(\frac{T}{T_{\mathrm{RH}}}\right)^{8},$
(2.7)

where $\alpha\equiv\pi^{2}g_{*}/30$ with $g_{*}$ the effective degrees of
freedom for the thermal plasma density at $T_{\mathrm{RH}}$ and ${M}_{P}\equiv
1/\sqrt{8\pi G_{N}}\simeq 2.4\times 10^{18}$ GeV is the reduced Planck mass.
Note that we have defined the reheating temperature $T_{\rm RH}$ as the plasma
temperature when
$\rho_{\phi}(T_{\mathrm{RH}})=\rho_{R}(T_{\mathrm{RH}})=\rho_{\rm RH}$. With
the approximate matter domination down to $T_{\mathrm{RH}}$, the reheating
temperature can be analytically expressed as [53, 54],

$T_{\mathrm{RH}}\approx\left(\frac{12}{25\alpha}\right)^{1/4}\left(\Gamma_{\phi}M_{P}\right)^{1/2}\,.$
(2.8)

Therefore, it is sufficient to know the temperature dependence of the
production rate, $R_{\cal X}(T)$ to solve equation (2.6) which determines the
dark matter abundance. The temperature dependence is different for production
initiated from the thermal bath compared to the bremsstrahlung production
process introduced here. The latter will dominate in different regimes
depending on the value of $T_{\mathrm{RH}}$.

For the annihilation process ${\cal PP}\rightarrow{\cal XX}$ as in Fig. 1(a),
we express the dark matter production rate for $T\gg m_{{\cal X}}$ as,

$R_{\cal PX}=\frac{T^{6+n}}{\Lambda^{2+n}}\,,$ (2.9)

where $\Lambda$ represents the UV scale. In fact, if one wants to be more
precise, the interactions between the dark matter ${\cal X}$ and the plasma
particles ${\cal P}$ can be parametrized by

$\frac{1}{\Lambda^{2+n}}=\frac{{\cal O}(1)}{\pi^{5}}\frac{1}{M^{2+n}}\,,$
(2.10)

where the factors of $\pi$ are being counted to monitor the phase space
factors of both the initial particle densities and the final states. These can
make up a few orders of magnitude and become important when comparing with DM
bremsstrahlung that is further suppressed by phase space. This
parameterization is applicable to a wide variety of scenarios. Indeed, $M$ can
then directly be identified with the mass of a heavy mediator in $SO(10)$ [15,
34, 36], or if a graviton is exchanged [28, 29, 30], $M\simeq M_{P}$, or
related to the SUSY breaking scale in supergravity constructions [53, 59, 60]
where $M\simeq\sqrt{m_{3/2}M_{P}}$, $m_{3/2}$ being the gravitino mass.

The same scale-suppressed effective interaction responsible for the
annihilation production of DM inevitably induces the bremsstrahlung
$\phi\rightarrow{\cal PPXX}$ as shown in Fig. 1(b). We can estimate the
branching ratio of the bremsstrahlung process by a mere dimensional analysis
with phase space counting. For example, for $n=2$

$\displaystyle{\rm Br}_{{\cal X}}^{{\rm brem}}$
$\displaystyle=\frac{\Gamma_{\phi\rightarrow{\cal
PPXX}}}{\Gamma_{\phi}}\equiv\tilde{c}\left(\frac{1}{4\pi^{2}}\right)^{2}\frac{\left(m_{\phi}/4\right)^{4}}{M^{4}}\,.$
(2.11)

For other values of $n$, the branching ratio will scale as
$(m_{\phi}^{2+n}/M^{2+n})$. The interaction determines the numerical factor
$\tilde{c}$ and may contain an implicit dependence on $n$. We will show some
typical examples of $\tilde{c}\sim{\cal O}(1)$ for $n=2$ in Sections 3 and 4.
The expression above has adopted the parametrization in Eq. (2.10). This
expression can be compared with the loop process given in Fig. 1c (assuming
the loop process does not vanish). The loop branching fraction was calculated
in [42] leading to

$\displaystyle{\rm Br}_{{\cal X}}^{{\rm loop}}$
$\displaystyle=\frac{\Gamma_{\phi\rightarrow{\cal
XX}}}{\Gamma_{\phi}}=\left(\frac{1}{4\pi^{2}}\right)^{2}\left(1+\frac{\pi^{4}}{4}\right)\frac{m_{\phi}^{4}}{M^{4}}\,.$
(2.12)

for $n=2$. In the case of bremsstrahlung we counted the phase space of each
extra particle in addition to the two body phase space as $1/(4\pi^{2})$
(these are the same factors of $\pi$ which appear in the loop calculation),
the factor $m_{\phi}/4$ reflects the typical energy of the decay product. As
one can see the loop process is typically two orders of magnitude larger than
bremsstrahlung (for ${\tilde{c}}\sim 1$). However as we have said repeatedly,
not all couplings lead to a non-zero rate for the loop in Fig. 1c.

The production rate (the number of dark matter particles produced per unit
volume and unit of time) for bremsstrahlung can then be written as

$R_{\phi{\cal X}}=2\Gamma_{\phi}{\rm Br}_{\cal X}^{\rm
brem}\rho_{\phi}/m_{\phi}\,.$ (2.13)

Dark matter can also be produced from plasma annihilation as long as it is
lighter than the maximum temperature $T_{{\rm max}}$ during the epoch of
inflaton oscillations [51, 52, 53]. If we combine Eq.(2.6), (2.7), and (2.9),
we see that $d{\cal Y}_{\cal X}/dT\propto T^{n-7}$. Thus, for $n<6$, DM is
mainly produced at lower temperature near the threshold at $T_{c}\simeq{\rm
Max}[T_{\mathrm{RH}},m_{{\cal X}}]\ll T_{{\rm max}}$, while for $n\geq 6$, the
dark matter is produced more efficiently at higher temperature up to $T_{{\rm
max}}$. The estimated dark matter density at the reheating temperature
$T_{\mathrm{RH}}$ is obtained by integrating Eq.(2.6) between $T_{\rm max}$
and $T_{c}$, with the rate given by Eq.(2.9). We obtain [42, 55],

$\displaystyle n_{{\cal X}}^{{\rm ann}}(T_{\mathrm{RH}})\simeq$
$\displaystyle\frac{8}{\sqrt{3\alpha}(6-n)}\frac{{M}_{P}T_{\mathrm{RH}}^{10}}{\Lambda^{n+2}T_{c}^{6-n}}\,,\quad
n<6\,,$ (2.14a) $\displaystyle n_{{\cal X}}^{{\rm
ann}}(T_{\mathrm{RH}})\simeq$
$\displaystyle\frac{8}{\sqrt{3\alpha}}\frac{{M}_{P}T_{\rm
RH}^{10}}{\Lambda^{8}}\ln\frac{T_{{\rm max}}}{T_{c}}\,,\quad n=6\,,$ (2.14b)
$\displaystyle n_{{\cal X}}^{{\rm ann}}(T_{\mathrm{RH}})\simeq$
$\displaystyle\frac{8}{\sqrt{3\alpha}(n-6)}\frac{{M}_{P}T_{\mathrm{RH}}^{10}T_{{\rm
max}}^{n-6}}{\Lambda^{n+2}}\,,\quad n>6\,.$ (2.14c)

It should be noted that there is an ${\cal O}(1)$ correction in the exact dark
matter relic density for $n\leq 6$ since the inflaton domination approximation
breaks down at $T_{\mathrm{RH}}$ [61, 53]. There is also a correction for the
production of dark matter produced before the radiation reaches the
temperature $T_{\rm max}$ [28]. Nevertheless, the important feature of these
solutions, is that for $n\geq 6$, the production is enhanced by the maximum
temperature during inflaton oscillations, i.e., at temperatures which can be
much larger than $T_{\mathrm{RH}}$.

For DM bremsstrahlung, integrating Eq. (2.6) with the rate given by Eq.(2.13)
between $T_{\mathrm{RH}}$ and $T_{{\rm max}}$ gives

$\displaystyle n_{{\cal X}}^{{\rm brem}}(T_{\mathrm{RH}})$
$\displaystyle\simeq\frac{10\alpha}{3}{\rm Br}_{{\cal X}}^{{\rm
brem}}\frac{T_{\mathrm{RH}}^{4}}{m_{\phi}}\,,$ (2.15)

and

$\displaystyle n_{{\cal X}}^{{\rm brem}}/n_{{\cal X}}^{{\rm ann}}\simeq$
$\displaystyle\frac{10\alpha^{3/2}(6-n)\tilde{c}\pi}{512\sqrt{3}}\frac{\left(m_{\phi}/4\right)^{n+1}T_{c}^{6-n}}{{M}_{P}T_{\mathrm{RH}}^{6}}\,,\quad
n<6\,,$ (2.16a) $\displaystyle n_{{\cal X}}^{{\rm brem}}/n_{{\cal X}}^{{\rm
ann}}\simeq$
$\displaystyle\frac{10\alpha^{3/2}\tilde{c}\pi}{512\sqrt{3}}\frac{\left(m_{\phi}/4\right)^{7}}{{M}_{P}T_{\mathrm{RH}}^{6}}\left(\ln\frac{T_{{\rm
max}}}{T_{c}}\right)^{-1}\,,\quad n=6\,,$ (2.16b) $\displaystyle n_{{\cal
X}}^{{\rm brem}}/n_{{\cal X}}^{{\rm ann}}\simeq$
$\displaystyle\frac{10\alpha^{3/2}(n-6)\tilde{c}\pi}{512\sqrt{3}}\frac{\left(m_{\phi}/4\right)^{n+1}}{{M}_{P}T_{\mathrm{RH}}^{6}T_{{\rm
max}}^{n-6}}\,,\quad n>6\,,$ (2.16c)

assuming 2 $\chi$’s are produced in each process.

In Fig. 2, we show the parameter space which achieves the observed dark matter
relic density inferred from CMB observations [62] for $n=0$, $n=2$, $n=4$ and
$n=6$ with $m_{\phi}=10^{13}$GeV, $g_{*}=100$ and $\tilde{c}=1$. For $n=6$ we
set $\ln(T_{\rm max}/T_{c})=1$ for illustration. The black contours show the
required value of $\Lambda$ in GeV for each choice of $m_{\chi}$ and
$T_{\mathrm{RH}}$. The red-dashed lines show the relative importance of the
bremsstrahlung contribution relative to production via plasma annihilation,
$n_{{\cal X}}^{{\rm brem}}/n_{{\cal X}}^{{\rm ann}}$, as given in Eq. (2.16).
For $n=0$, the ratio $n_{{\cal X}}^{{\rm brem}}/n_{{\cal X}}^{{\rm ann}}\simeq
4\times 10^{-5}(T_{c}/T_{\mathrm{RH}})^{6}$ and thus the bremsstrahlung
process only dominates when $m_{\chi}>T_{\mathrm{RH}}$. 333The tiny wiggles
near the first turning point in the black contours in Fig. 2a are artifacts of
the limited numerical resolution of the plot.

For larger $n$, it is interesting to note that the bremsstrahlung process
dominates when $m_{\chi}<T_{\mathrm{RH}}$, for sufficiently low
$T_{\mathrm{RH}}$. For example, for $n=2$, for reheating temperatures below
$T_{\mathrm{RH}}\lesssim 10^{10}$ GeV (with $m_{\chi}<10^{10}$ GeV), the
bremsstrahlung production process, which to date has never been taken into
account in the literature, is the dominant process which populates the dark
matter in the Universe. Indeed, taking a dimension-six effective interaction
($n=2$) between DM and the plasma (SM) and $m_{{\cal X}}\ll T_{\mathrm{RH}}$,
one finds for example

$\frac{n_{{\cal X}}^{{\rm brem}}}{n_{{\cal X}}^{{\rm
ann}}}\simeq{1.7}\,\tilde{c}\times\left(\frac{g_{*}}{100}\right)^{3/2}\left(\frac{10^{10}\,{\rm
GeV}}{T_{\mathrm{RH}}}\right)^{2}\left(\frac{m_{\phi}}{10^{13}{\rm
GeV}}\right)^{3}\,.~{}$ (2.17)

For $m_{\phi}=10^{13}\,$GeV, the bremsstrahlung processes dominate for low
reheating temperature $T_{\mathrm{RH}}\lesssim 10^{10}$ GeV corresponding to
$\lambda_{i}\lesssim 4\times 10^{-5}$, and agrees with our numerical result
showed in Fig. 2. This is one of the main results of our work. Moreover, we
show in the Appendix that this result is quite robust, even when taking into
account the possibility for production through annihilation during the
thermalization phase of reheating [47, 46].

We can also observe a change in the slope of $T_{\mathrm{RH}}$ vs. $m_{\chi}$
for fixed $\Lambda$, when the black curves cross $n_{{\cal X}}^{{\rm
brem}}/n_{{\cal X}}^{{\rm ann}}\approx 1$ as shown by the red lines. This is
due to the strong difference in the dependence on $T_{\mathrm{RH}}$ for the
two processes. Whereas the bremsstrahlung process can be seen as a mere decay
process (rescaled by a constant branching fraction), which implies $n_{\cal
X}^{\rm brem}\propto T_{\mathrm{RH}}^{4}$, the annihilation channels has a
stronger dependence on the temperature due to its origin, proportional to the
density squared of the particles in the bath, thus $n_{\cal X}^{\rm
ann}\propto T_{\mathrm{RH}}^{10}/T_{c}^{6-n}$, see Eqs.(2.14). For $m_{\cal
X}<T_{\mathrm{RH}}$, $T_{c}=T_{\mathrm{RH}}$, and $n_{\cal X}^{\rm ann}\propto
T_{\mathrm{RH}}^{4+n}$. Otherwise, $n_{\cal X}^{\rm ann}\propto
T_{\mathrm{RH}}^{10}/m_{\cal X}^{6-n}$. In the specific case of $n=0$, the
slope for $T_{\mathrm{RH}}$ vs. $m_{\cal X}$ for the bremsstrahlung and
annihilation processes are the same. This is because annihilations dominate
when $m_{\cal X}<T_{\mathrm{RH}}$, and $n_{\cal X}^{\rm ann}\propto
T_{\mathrm{RH}}^{4}$ as it does for the bremsstahlund process which dominates
at higher masses. This is clearly seen in the upper left panel of Fig. 2. For
$n\neq 0$, the density from annihilations has a steeper dependence on the
temperature. The relic density today as a fraction of critical density,
$\Omega_{\cal X}h^{2}\propto n_{\cal X}m_{\cal X}/T_{\mathrm{RH}}^{3}$ and we
obtain $\Omega^{\rm brem}_{\cal X}h^{2}\propto m_{\cal X}T_{\mathrm{RH}}$ for
the bremsstrahlung production, and $\Omega^{\rm ann}_{\cal X}h^{2}\propto
m_{\cal X}T_{\mathrm{RH}}^{1+n}$ for production from the thermal plasma (when
$m_{\cal X}<T_{\mathrm{RH}}$). Curves of constant relic density will then
follow a steeper slope for the bremsstrahlung process than for the
annihilation, and this corresponds to the result we obtained numerically in
Fig. 2. Note also the change of sign in the slope for $n=0$ when one reaches
the point $m_{\cal X}=T_{\mathrm{RH}}$. This comes from the fact that, for a
given mass, the time needed to produce dark matter is reduced (from
$T_{\mathrm{max}}$ to $m_{\cal X}$ rather than from $T_{\mathrm{max}}$ to
$T_{\mathrm{RH}}$). In other words, to produce the same amount of dark matter,
one needs larger $T_{\mathrm{max}}$, and as a consequence, larger
$T_{\mathrm{RH}}$. The slope recovers the $m_{\cal X}\propto
1/{T_{\mathrm{RH}}}$ behavior characteristic of inflaton decay, when
$R_{\phi\cal X}>R_{{\cal P}\cal X}$.

Figure 2:  Parameter values that achieve the observed dark matter relic
density [62] for $n=0$, $n=2$, $n=4$, and $n=6$ with $m_{\phi}=10^{13}$GeV,
$g_{*}=100$ and $\tilde{c}=1$. The black contours show the required parameter
$\Lambda$ in GeV; The red line presents the relative importance $n_{{\cal
X}}^{{\rm brem}}/n_{{\cal X}}^{{\rm ann}}$ of bremsstrahlung vs. plasma
annihilation production for the choice of $\Lambda$ that achieves a realistic
DM relic density.

It is also interesting to note the change in behavior of the ratio $n_{{\cal
X}}^{{\rm brem}}/n_{{\cal X}}^{{\rm ann}}$ as function of the dark matter
mass. Indeed, whereas the bremsstrahlung production depends only on the
reheating temperature (assuming we are well below the kinematic threshold),
thermal production for $n<6$ depends strongly on the lower bound of the
integration in Eq. (2.6), $T_{c}$. As long as $m_{\cal X}<T_{\mathrm{RH}}$,
the lower bound is $T_{\mathrm{RH}}$ and the thermal production is independent
of $m_{\cal X}$, as we can see in Fig. 2 with $n=2$, $4$ and $6$, where this
regime corresponds to the horizontal red dashed lines. However, if $m_{\cal
X}>T_{\mathrm{RH}}$, the thermal production is less efficient because it takes
place over a shorter period of time, and for larger values of $m_{\cal X}$,
the period is shorter. This explains why, in the case of large dark matter
masses, the bremsstrahlung production largely dominates the thermal one. To be
more precise, For $T_{\rm max}\gtrsim m_{{\cal X}}\gtrsim T_{\mathrm{RH}}$,
the production from plasma annihilation is suppressed by a factor of
$T_{\mathrm{RH}}^{4}/m_{\cal X}^{4}$ due to the production threshold at
$T\simeq m_{\cal X}$, whilst for bremsstrahlung production, the only threshold
is $2m_{{\cal X}}<m_{\phi}$. This suppression factor stems from the factor of
$T_{c}^{6-n}$ in Eq. (2.16a).

This domination of bremsstrahlung production is more pronounced for higher $n$
as can be seen in Fig. 2 for $n=4$ and $6$. For $m_{\chi}<T_{\mathrm{RH}}$ and
$n\leq 6$ the ratio in Eq. (2.16) is proportional to
$\left({m_{\phi}}/{T_{\mathrm{RH}}}\right)^{n}$. For $n=6$, bremsstrahlung
production nearly dominates all of the parameter space with
$T_{\mathrm{RH}}<m_{\phi}$. For even larger values of $n>6$, ${n_{{\cal
X}}^{{\rm brem}}}/{n_{{\cal X}}^{{\rm ann}}}$ is independent of $m_{\chi}$ (so
long as it is kinematically allowed) and receives additional enhancement
compared to $n\leq 6$ as is evident from Eq. (2.16c). This ratio can easily
reach $10^{7}$ for typical values of
$\frac{T_{\mathrm{max}}}{T_{\mathrm{RH}}}=100$. We conclude therefore that not
only is bremsstrahlung production dominant for reheating temperatures
$T_{\mathrm{RH}}\lesssim 10^{10}$ GeV in the case $n=2$ (corresponding to the
exchange of a massive particle) but the dominance is even stronger for larger
values of $n$. This is another important result of our work. We recall as
well, that for $n=0$, bremsstrahlung dominates when
$m_{\chi}<T_{\mathrm{RH}}$.

Thus, bremsstrahlung production dominates in a wide range of reheating
temperatures. Combining Eqs. (2.11) and Eq. (2.15) while writing the relic
abundance as [9]

$\Omega_{\cal X}h^{2}=1.6\times
10^{8}\,\frac{g_{0}}{g_{*}}\frac{n(T_{\mathrm{RH}})}{T_{\mathrm{RH}}^{3}}\frac{m_{\cal
X}}{1~{}{\rm GeV}}\,,$ (2.18)

with $g_{0}=3.91$, we obtain

$\frac{\Omega^{\rm brem}_{\cal X}h^{2}}{0.12}\simeq\tilde{c}\,\frac{m_{\cal
X}}{{2.3}\times 10^{8}\,{\rm GeV}}\left(\frac{m_{\phi}}{10^{13}\,{\rm
GeV}}\right)^{3}\frac{T_{\mathrm{RH}}}{10^{10}\,{\rm
GeV}}\left(\frac{10^{16}{\rm GeV}}{\Lambda}\right)^{4}\,,$ (2.19)

or

$\frac{\Omega^{\rm brem}_{\cal X}h^{2}}{0.12}\simeq\tilde{c}\,\frac{m_{\cal
X}}{{7.0}\times 10^{10}\,{\rm GeV}}\left(\frac{m_{\phi}}{10^{13}\,{\rm
GeV}}\right)^{3}\frac{T_{\mathrm{RH}}}{10^{10}\,{\rm
GeV}}\left(\frac{10^{16}{\rm GeV}}{M}\right)^{4}\,,$ (2.20)

for $n=2$ where we have normalized to the observed cold dark matter density
[62]. The threshold condition $m_{\cal X}<m_{\phi}\simeq 10^{13}\,{\rm GeV}$
sets $\Lambda\lesssim 4\times 10^{17}$ GeV for $T_{\mathrm{RH}}<10^{12}$ GeV.
Therefore, to achieve the observed DM abundance, the parameter $\Lambda$ of
the interaction can lie between $m_{\phi}$ and $4\times 10^{17}$ GeV. Note
that this range is largely independent of the reheating temperature, since for
a given $m_{\cal X}$, $\Lambda\propto T_{\mathrm{RH}}^{\frac{1}{4}}$. Within
the viable parameter space according to this estimate, probably the most
interesting choice for the interaction scale is the GUT scale, $\Lambda\simeq
10^{16}\,$GeV. In this case, the dark sector and the visible sector is only
connected by GUT scale interactions such as through heavy gauge bosons of a
unification group. This induces vector-type effective interaction between the
plasma particle and the DM. As one can see from Eq. (2.20), the relic density
(for $n=2$) depends on four parameters,
$m_{\cal\chi},m_{\phi},T_{\mathrm{RH}}$, and $M$. Lower dark matter masses can
be obtained raising $T_{\mathrm{RH}}$ and lowering $M$ for fixed inflaton
mass, $m_{\phi}=10^{13}$ GeV. However if we raise $T_{\mathrm{RH}}$, the
bremsstrahlung process will no longer dominate. Thus the lowest dark matter
mass determined by bremsstrahlung is obtained for $T_{\mathrm{RH}}\simeq
10^{10}$ GeV and $M\gtrsim m_{\phi}$ giving $m_{\cal X}\gtrsim 70$ MeV, for
$\tilde{c}\sim 1$.

For a more general inflaton mass, the bremsstrahlung production may still
dominate for sufficiently small $T_{\mathrm{RH}}$ and/or $T_{\mathrm{RH}}\ll
m_{\cal X}$. However, the DM relic density poses a lower bound to $m_{\phi}$
because of the production threshold $m_{\phi}>2m_{\cal X}$. Imposing
$m_{\phi}\gg m_{\cal X}/2,\,T_{\mathrm{RH}}$, Eq.(2.20) sets

$m_{\phi}\gtrsim 1.1\times 10^{12}\,{\rm GeV}\left(\frac{M}{10^{16}\,{\rm
GeV}}\right)^{4/5}$ (2.21)

for bremsstrahlung to dominate the production of the observed DM relic.

The estimates we have made up to now are mostly qualitative in that the model
dependent factor $\tilde{c}$ is not guaranteed to be ${\cal O}(1)$. In the
remainder of the paper, we perform numerical calculations for the rate of DM
bremsstrahlung and plasma annihilation for dark matter produced through a
gravity portal or a vector portal. We will see that the result is in close
agreement with the estimate of $\tilde{c}\sim{\cal O}(1)$ in this section.

##  3 Gravitational bremsstrahlung

A well motivated choice for the scale $\Lambda$ is the Planck scale. Indeed,
DM will always interact with the inflaton or plasma particle through the
exchange of a graviton, $h_{\mu\nu}$. In this section, we evaluate the dark
matter abundance produced by gravitational bremsstrahlung that is unavoidable
in inflationary cosmology. As the inflaton decays, DM particles can be emitted
from the decay products or from the inflaton itself through a virtual
graviton, leading to the 4-body final state,

$\phi\rightarrow ss^{*}(\text{or}\,ff)+XX(\text{or}\,\bar{\chi}\chi)\,,$ (3.1)

where we denote by $X$ a real scalar, or $\chi$, a Dirac fermion as the dark
matter candidate. The relevant diagrams are shown in Fig. 3. The analysis can
be easily generalized to complex scalars or Majorana fermions as DM by
counting degrees of freedoms. We note that contracting the two plasma particle
lines in Fig. 3 to close a loop generating a process $\phi\rightarrow{\cal
P}_{\rm loop}\rightarrow h_{\mu\nu}\rightarrow{\cal X}~{}\cal{X}$ (depicted in
Fig. 1(c) replacing the effective coupling by a graviton) gives zero since the
spin-$2$ massless graviton cannot mix with the spin-$0$ inflaton. This is in
contrast to the cases treated in [42]. Therefore, the gravitational
bremsstrahlung is the leading radiative production process of dark matter. We
will show that even when bremsstrahlung is the dominant gravitational
production mechanism in some parts of the parameter space, it cannot saturate
the observed dark matter abundance alone. In general, the dominant
gravitational source of dark matter comes from the gravitational annihilation
of inflaton [28, 29]. In a sense, our result confirms the approximate result
obtained in the previous section, since the exchange of a graviton can be seen
as an effective model with $\Lambda=M_{P}>10^{17}$ GeV.

Figure 3: .The diagrams for gravitationally induced DM bremsstrahlung from
inflaton decay.

To evaluate the bremsstrahlung rate, we use the Monte Carlo event generator
MadGraph [63] to calculate the matrix element and simulate the final 4-body
phase space. MadGraph supports all of the Lorentz structure of vertices [64]
for a massless graviton444The massive graviton propagator does not reduce to
the massless one in the zero mass limit because of the extra degrees of
freedom. [65, 66] that appear in Fig. 3 which is crucial in the cancellation
of collinear divergences. We use FeynRules [67] to generate the UFO model
files for MadGraph.

Some comments on the infrared behavior of the process is in order to justify
the validity of the numerical simulation. A massless mediator may lead to soft
or collinear divergences like in QED or QCD which needs to be regulated with a
cut in the particle momenta. Fortunately, gravity has no collinear divergences
[68] and thus no subtle treatment is required. Indeed, we have checked
analytically that all of the soft and collinear divergences in Fig. 3 cancel
but only after summing over all of the four contributing diagrams. This
introduces subtleties in numerical evaluation with MadGraph. By default,
MadGraph handles the particle widths with Breit-Wigner propagators. 555Many
thanks to Olivier Mattelaer for pointing this out in the MadGraph5_aMC@NLO
question and answer launchpad:
https://answers.launchpad.net/mg5amcnlo/+question/698109 The width in
denominators of these propagators mix contributions of different orders of the
perturbative calculation and spoils the delicate cancellation of the IR
divergences. To maintain the cancellation at a fixed order, we set the widths
of all particles to zero when simulating with MadGraph. The cancellation of
large numbers from almost on-shell propagators between diagrams could pose
another challenge to the numerical stability. To evade this issue, we applied
a cut on the DM invariant mass, $(\sum p_{{\rm DM}})^{2}>10^{-8}m_{\phi}^{2}$,
to get rid of very low energy DM in the numerical simulation. Due to the
absence of soft divergences, the DM particles are rarely produced with soft or
near collinear momenta. Thus, this cut guarantees the numerical stability with
only mild effects on the resultant inflaton partial width. 666Practically,
such a cut is not even needed in a wide range of inflaton/DM mass ratios after
scaling down $M_{P}$ thanks to the numerical precision of MadGraph. This can
be used as a stringent consistency check of the implementation of the model.
The numerical result for the inflaton partial width due to the DM
bremsstrahlung, in the limit of $m_{{\cal X}}\ll m_{\phi}$ is

$\displaystyle\Gamma_{\phi,{\rm brem}}^{G}$ $\displaystyle={\cal C}_{{\cal
PX}}^{G}\lambda_{i}^{2}\frac{m_{\phi}^{5}}{{M}_{P}^{4}}\,.~{}$ (3.2)

where the coupleings $\lambda_{i=s,f}$ are given in Eq. (2.4). We can now
evaluate the coefficients ${\tilde{c}}$ for $M=M_{P}$ and $n=2$,

${\tilde{c}}=2^{16}\pi^{5}{\cal C}_{{\cal PX}}^{G}\,,$ (3.3)

and the numerical factor ${\cal C}_{{\cal PX}}^{G}$ is obtained by simulation,

$\displaystyle{\cal C}_{sX}^{G}$ $\displaystyle\simeq 2.4\times 10^{-9}\,,$
(3.4a) $\displaystyle{\cal C}_{s\chi}^{G}$ $\displaystyle\simeq 4.1\times
10^{-9}\,,$ (3.4b) $\displaystyle{\cal C}_{fX}^{G}$ $\displaystyle\simeq
2.6\times 10^{-9}\,,$ (3.4c) $\displaystyle{\cal C}_{f\chi}^{G}$
$\displaystyle\simeq 4.3\times 10^{-9}\,,$ (3.4d)

for $\phi\rightarrow ss^{*}XX$, $ss^{*}\bar{\chi}\chi$, $ffXX$, and
$ff\bar{\chi}\chi$, respectively. From Eq. (2.20), we see that the
gravitational bremsstrahlung process alone cannot produce enough DM particles
to saturate the observed relic abundance when $M=M_{P}$ and $m_{\cal
X},T_{\mathrm{RH}}<m_{\phi}$.

The dark matter particle can also be produced gravitationally by the
annihilation of the inflaton during its oscillationary epoch or by the
annihilation of a pair of plasma particle during or after reheating. These two
scenarios have been considered widely in previous works such as [19, 20, 22,
23, 28, 29]. The production from inflaton annihilation was found to be
dominant whenever the gravitational production saturates the observed
abundance777see also [69, 70] for an alternative treatment of the
gravitational production.. The gravitational annihilation of the inflaton may
also leave a significant impact on the maximal or the reheating temperature of
the thermal plasma [29, 71, 72, 30]. We briefly summarize the result for
inflaton annihilation [28, 29] as follows. For inflaton annihilation
$\phi\phi\rightarrow XX$ or $\bar{\chi}\chi$ during oscillations about a
quadratic potential, the DM production rate that enters the Boltzmann equation
Eq. (2.5c) is

$R_{\phi\phi}^{G}=\frac{\rho_{\phi}^{2}}{128\pi{M}_{P}^{4}}\beta_{\phi\phi}\,,$
(3.5)

where $\beta_{\phi\phi}=1$ for $\phi\phi\rightarrow XX$ and
$\beta_{\phi\phi}=m_{\chi}^{2}/m_{\phi}^{2}$ for
$\phi\phi\rightarrow\bar{\chi}\chi$, assuming $m_{\chi,X}\ll m_{\phi}$. All
possible DM degrees of freedom have been summed in this inclusive production
rate. The production of fermionic DM is suppressed by its mass compared to the
case of scalar DM since the DM particle, anti-particle pair must have the same
helicity or the opposite chiralities. In contrast, the production rate from
bremsstrahlung does not suffer from this suppression because the pair of
plasma particles ($ss^{*}$ or $ff$ ) in the decay product carry angular
momentum so fermionic DM does not have to be in the same helicity state. The
produced DM density at $T_{\mathrm{RH}}$ is [28, 29]

$n_{{\cal X}}^{\phi\phi}(T_{\mathrm{RH}})\simeq\frac{\sqrt{3}\rho_{\rm
RH}^{\frac{3}{2}}}{192\pi{M}_{P}^{3}}\left(\frac{\rho_{{\rm end}}}{\rho_{\rm
RH}}\right)^{\frac{1}{2}}\beta_{\phi\phi}\,,$ (3.6)

where $\rho_{{\rm end}}$ is the energy density stored in the inflaton field at
the end of the slow roll inflation. At the start of inflaton oscillations, the
high energy density of the inflaton condensate significantly enhances the
production rate of dark matter. Using Eq. (2.18) and Eq. (3.6), we obtain [28,
29]

$\frac{\Omega_{\chi}^{\phi\phi}h^{2}}{0.12}\simeq\left(\frac{10^{13}~{}{\rm
GeV}}{m_{\phi}}\right)^{2}\left(\frac{T_{\mathrm{RH}}}{10^{10}~{}{\rm
GeV}}\right)\left(\frac{\rho_{\rm end}}{10^{64}~{}{\rm
GeV}^{4}}\right)^{\frac{1}{2}}\left(\frac{m_{\chi}}{6.4\times 10^{10}~{}{\rm
GeV}}\right)^{3}$ (3.7)

for the relic abundance of a fermionic dark matter candidate, $\chi$. The
range of DM masses for which the bremsstrahlung process dominates the
production of dark matter relative to inflaton scattering is found by
comparing Eqs. (3.7) and (2.20). For ${\tilde{c}}\simeq 0.08$, this leads to

$m_{\chi}\lesssim
0.3~{}\left(\frac{M_{P}}{M}\right)^{2}\left(\frac{m_{\phi}}{10^{13}~{}{\rm
GeV}}\right)^{5/2}\left(\frac{10^{16}~{}{\rm GeV}}{\rho_{\rm
end}^{\frac{1}{4}}}\right)~{}{\rm PeV}\,,$ (3.8)

However, for these masses, the production is much too low to match the CMB-
inferred relic density [62], as we can see from Eq. (2.20). The production of
scalars, $X$, is even less efficient since the gravitational inflaton
scattering does not depend on the mass of the dark matter, and hence always
dominates over the bremsstrahlung process. We can then safely neglect the
gravitational bremsstrahlung process for the rest of our work.

##  4 Vectorial bremsstrahlung

In GUT-inspired models, there are typically gauge bosons with mass $\gtrsim
10^{15}$ GeV, i.e., much above the mass of the inflaton or $T_{\mathrm{max}}$.
This type of heavy mediator also appears in see-saw constructions. In this
section, we proceed to discuss DM bremsstrahlung during inflaton decay in the
general cases where the DM ${\cal X}$ and the plasma particle ${\cal P}$
interact through a heavy vector mediator. The dark matter candidate can be a
_complex_ scalar ${\bf X}$ or a Dirac spinor $\chi$. We will stay within the
limit of $m_{{\rm{\cal X}}}\ll m_{\phi}$, i.e., the decay products of the
inflaton are effectively massless in the amplitude.

We consider the following dimension-six vector portal effective interactions
suppressed by a cut-off scale $M\gg m_{\phi}$,

$\displaystyle{\cal L}_{s{\bf X}}^{V}$ $\displaystyle=\frac{g_{s{\bf
X}}^{2}}{M^{2}}\left(s^{*}\partial^{\mu}s-s\partial^{\mu}s^{*}\right)\left({\bf
X}^{*}\partial_{\mu}{\bf X}-{\bf X}\partial_{\mu}{\bf X}^{*}\right)\,,$ (4.1a)
$\displaystyle{\cal L}_{s\chi}^{V}$
$\displaystyle=\frac{g_{s\chi}^{2}}{M^{2}}{\rm
i}\left(s^{*}\partial_{\mu}s-s\partial_{\mu}s^{*}\right)\bar{\chi}\gamma_{\mu}\chi\,,$
(4.1b) $\displaystyle{\cal L}_{f{\bf X}}^{V}$ $\displaystyle=\frac{g_{f{\bf
X}}^{2}}{M^{2}}{\rm i}\bar{f}\gamma_{\mu}f\left({\bf X}^{*}\partial_{\mu}{\bf
X}-{\bf X}\partial_{\mu}{\bf X}^{*}\right)\,,$ (4.1c) $\displaystyle{\cal
L}_{f\chi}^{V}$
$\displaystyle=\frac{g_{f\chi}^{2}}{M^{2}}\bar{f}\gamma_{\mu}f\bar{\chi}\gamma_{\mu}\chi\,.$
(4.1d)

The fermionic bilinears with $\gamma_{5}$ are equivalent to the ones without
it in the massless limit since the phases of chiral spinors can be redefined
independently. These interactions can be obtained by integrating out a heavy
gauge boson that interact with both the DM particle and the plasma particle as
we will demonstrate briefly at the end of this section. In such a scenario,
the real scalar inflaton condensate cannot annihilate efficiently through the
vector portal because it requires the presence of another “partner scalar”
field to form a complex scalar to interact with a gauge boson at the tree
level. Although the inflaton is often the real part of a complex field, the
imaginary part is far less abundant during inflaton oscillations. Thus, the
annihilation of the inflaton to DM can only occur through one-loop diagrams
with two virtual gauge bosons and one virtual “partner scalar”, and is thus
suppressed by at least an additional factor of $m_{\phi}^{4}/M^{4}$. Moreover,
the radiative decay of the inflaton to dark matter $\phi\rightarrow{\bf
X^{*}X}$ or $\phi\rightarrow\bar{\chi}\chi$ induced by a plasma particle loop
as in Fig. 1(c) vanishes for vector type effective interactions because a
massive scalar cannot mix with a vector mediator. Technically, since the
plasma particle loop only has one Lorentz index, it must be proportional to
the inflaton 4-momentum $p_{\phi}^{\mu}$ and vanishes after contracting with
the current of the final state DM particles. Thus, as in the case of
gravitational portal, bremsstrahlung is the leading radiative production
mechanism.

Following the same procedure as in the previous section, the DM bremsstrahlung
rate of $\phi\rightarrow{\cal PPXX}$ induced by the above interactions for
$M\gg m_{\phi}\gg m_{s},m_{f}$ is parameterized as,

$\Gamma_{\phi\rightarrow{\cal PPXX}}^{V}={\cal C}_{{\cal PX}}^{V}g_{{\cal
PX}}^{4}\lambda_{i}^{2}\frac{m_{\phi}^{5}}{M^{4}}\,,$ (4.2)

with ${\cal C}_{{\cal PX}}^{V}$ a numerical factor and $g_{{\cal PX}}$ the
dimensionless coefficient of the effective operators as defined in Eq. (4.1).
In this case,

${\tilde{c}}=2^{16}\pi^{5}{\cal C}_{{\cal PX}}^{V}g_{{\cal PX}}^{4}\,.$ (4.3)

As we did for the gravitational coupling, we use MadGraph for numerical
simulation of the 4-body decay. We obtained

$\displaystyle{\cal C}_{s{\bf X}}^{V}$ $\displaystyle=1.0\times 10^{-8}\,,$
(4.4a) $\displaystyle{\cal C}_{s\chi}^{V}$ $\displaystyle=4.1\times
10^{-8}\,,$ (4.4b) $\displaystyle{\cal C}_{f{\bf X}}^{V}$
$\displaystyle=6.9\times 10^{-9}\,,$ (4.4c) $\displaystyle{\cal
C}_{f\chi}^{V}$ $\displaystyle=2.8\times 10^{-8}\,,$ (4.4d)

for $\phi\rightarrow ss^{*}{\bf X}{\bf X}^{*}$, $\phi\rightarrow
ss^{*}\chi\bar{\chi}$, $\phi\rightarrow f\bar{f}{\bf X}{\bf X}^{*}$ and
$\phi\rightarrow f\bar{f}\chi\bar{\chi}$, respectively.

On the other hand, the rate of DM production by annihilation of particles from
the plasma ${\cal P}$ is evaluated as [29]

$R^{\cal P}_{\rm ann}=\sum_{i_{1},i_{2}}\frac{1}{1024\pi^{6}}\int{\rm
d}E_{1}{\rm d}E_{2}E_{1}E_{2}f_{i_{1}}(E_{1})f_{i_{2}}(E_{2}){\rm
d}\cos\theta_{12}{\rm d}\Omega_{13}\left|{\cal M}^{i_{1}i_{2}}\right|^{2}\,,$
(4.5)

where $i_{1}$ and $i_{2}$ run over all annihilating particle species. The
angles and the energies are related to the Mandelstam variables by
$t=s(\cos\theta_{13}-1)/2$ and $s=2E_{1}E_{2}(1-\cos\theta_{12})$. For the
vector portal interactions parametrized by Eqs. (4.1), the production rates
are given by

$\displaystyle R_{\rm ann}^{s{\bf X}}$ $\displaystyle=\frac{g_{s{\bf
X}}^{4}\pi^{3}T^{8}}{8100M^{4}}\,,$ (4.6a) $\displaystyle R_{\rm
ann}^{s{\chi}}$
$\displaystyle=\frac{g_{s{\chi}}^{4}\pi^{3}T^{8}}{2025M^{4}}\,,$ (4.6b)
$\displaystyle R_{\rm ann}^{f{\bf X}}$ $\displaystyle=\frac{49g_{f{\bf
X}}^{4}\pi^{3}T^{8}}{259200M^{4}}\,,$ (4.6c) $\displaystyle R_{\rm
ann}^{f{\chi}}$
$\displaystyle=\frac{49g_{f{\chi}}^{4}\pi^{3}T^{8}}{64800M^{4}}\,.$ (4.6d)

We are now able to compare the three modes of dark matter production : the
gravitational production through inflaton scattering, the bremsstrahlung
production, and the scattering from particles in the plasma. In Fig. 4, we
present the parameter space for vector-type interactions that achieves the
observed dark matter relic density considering all production modes. The upper
and lower panels are for ${\cal L}_{f\chi}^{V}$ and ${\cal L}_{f{\bf X}}^{V}$,
respectively. The left panels and the right panels corresponds to $\rho_{{\rm
end}}=(10^{16}{\rm GeV})^{4}$ and $(10^{14}{\rm GeV})^{4}$. In all cases, we
assumed $m_{\phi}=10^{13}$ GeV and $g_{*}=100$. The black contours show the
required effective interaction scale $M/|g_{f{\bf X}}|$ or $M/|g_{f\chi}|$.
The gray shaded region is where inflaton annihilation through gravity
overproduces the DM. As this contribution can not be removed, the shaded
region is excluded. Along the boundary of the shaded region, the relic density
produced by inflaton annihilations mediated by gravity matches the observed
cold dark matter density without any additional contributions. The red lines
are the contours when the relative importance of bremsstrahlung vs other
production mechanisms combined $n_{{\cal X}}^{{\rm brem}}/(n_{{\cal X}}^{{\rm
other}})=0.01$, $1$, and $100$ for the choice of $M/g_{f{\cal X}}$ that
attains a realistic DM relic density equal to that inferred from CMB
observations [62]. The red lines may have sharp turns at $m_{\cal
X}=T_{\mathrm{RH}}$ due to the suppression of plasma annihilation at $m_{\cal
X}>T_{\mathrm{RH}}$. As we have seen previously, the DM bremsstrahlung
dominates at smaller reheating temperatures of $T_{\mathrm{RH}}\lesssim
10^{10}$GeV, whereas to achieve the observed DM relic density, the scale $M$
of the effective interaction can lie in a wide range above $m_{\phi}$ up to
$M\sim 10^{16}$GeV, in accordance with our estimate in Section 2. As expected,
the allowed parameter space is larger for lower values of $\rho_{\rm end}$ due
to the smaller amount of dark matter produced gravitationally.

Figure 4: Parameter space that realizes the observed dark matter relic density
for dimension-6 operators with $m_{\phi}=10^{13}$GeV. The upper and lower
panels are for ${\cal L}_{f\chi}^{V}$ and ${\cal L}_{f{\bf X}}^{V}$,
respectively. The left panels and right panels corresponds to $\rho_{{\rm
end}}=(10^{16}{\rm GeV})^{4}$ and $(10^{14}{\rm GeV})^{4}$. The black contours
show the required effective interaction scale, $M/g_{f{\cal X}}$ in GeV. The
red line presents the relative importance $n_{{\cal X}}^{{\rm brem}}/(n_{{\cal
X}}^{{\rm other}})$ of bremsstrahlung vs other production mechanisms combined
for the choice of $M/g_{f{\cal X}}$ that achieves a realistic DM relic. The
gray shaded region is where inflaton annihilation through gravity overproduces
the DM.

As a prove of concept, we outline a construction of a gauged $U(1)_{B-L}$
extension of the Standard Model that contains all of the ingredients discussed
in this section. The model can be realized as an intermediate subgroup of the
$SO(10)$ grand unification [73, 74]. In this model, the $U(1)_{B-L}$ is
spontaneously broken by the vacuum expectation value (vev) $v_{s}$ of a heavy
Higgs field $S$ with $(B-L)_{S}=2$, which can be the Majoron [75, 38, 39]. The
right-handed neutrinos $N_{i}$ have $(B-L)_{N}=-1$ and obtain a large majorana
mass through its Yukawa coupling to $S$ as in the seesaw mechanism [76, 77,
78, 79, 80, 81]. In addition, we introduce a Dirac spinor $\chi$ with
$(B-L)_{N}=-1$ as the dark matter and impose a $\mathbb{Z}{}_{2}$ symmetry for
its stability. 888The $\mathbb{Z}{}_{2}$ symmetry can be realized as a remnant
discrete symmetry of the broken $SO(10)$ unification group [82, 83, 84, 85].
The inflaton is taken to be a real singlet $\phi$. The interaction terms
relevant to the discussion are

${\cal L}\supset\bar{\chi}{\rm i}D_{\mu}\gamma^{\mu}\chi+{\rm
i}\bar{N}D_{\mu}\gamma^{\mu}N+\frac{y_{N}}{2}N^{T}CNS+\mu_{s}\left|S\right|^{2}\phi\,,$
(4.7)

where $D_{\mu}$ is the gauge covariant derivative and we have suppressed the
generation index of $N_{i}$. Expanding $S=v_{s}+s$, the $U(1)_{B-L}$ gauge
boson $V$ obtains a mass $m_{V}=2\sqrt{2}g_{B-L}v_{s}$. For $m_{V}$ and
$m_{s}$ much heavier than $m_{\phi}$, we can integrate out the heavy gauge
boson and the Higgs mode. The effective interactions become,

${\cal L}_{{\rm
eff}}\supset\frac{g_{B-L}^{2}}{m_{V}^{2}}\bar{\chi}\gamma^{\mu}\chi\bar{N}\gamma_{\mu}N+\frac{m_{N}\mu_{s}}{m_{s}^{2}}N^{T}CN\phi\,,$
(4.8)

where $m_{N}\equiv y_{N}v_{s}$ is the Majorana neutrino mass. Here, we have
obtained the Yukawa operator responsible for reheating and the dimension-6
operator for dark matter production. Comparing with Eq. (4.1) and Eq. (2.2b),
we can identify

$\displaystyle\lambda_{s}$ $\displaystyle=2\frac{m_{N}\mu_{s}}{m_{s}^{2}}\,,$
(4.9a) $\displaystyle\frac{g_{f\chi}^{2}}{M^{2}}$
$\displaystyle=\frac{g_{B-L}^{2}}{m_{V}^{2}}\,.$ (4.9b)

In this case, we have checked that if one considers a Standard Mode-like Higgs
mechanism, with $\mu_{s}\simeq m_{s}$, for right handed neutrino masses below
$\lesssim 10^{-4}m_{s}$ (leading to $T_{\mathrm{RH}}\lesssim 10^{10}$ GeV),
bremsstrahlung production radiated from a neutrino leg dominates all other
type of dark matter production. Note that in Eq. (4.7), we only present terms
that respect the chiral symmetry of $\chi$. It is technically natural to set
the chiral symmetry breaking parameters $m_{\chi}$ and the coupling of
$\bar{\chi}\chi\phi$ to small values to suppress the direct decay of $\phi$ to
a DM pair. In this case, the dominant DM production process is either the
plasma annihilation, the DM bremsstrahlung or the gravitationally annihilation
of the inflaton, as we have presented in Fig. 4.

##  5 Conclusion and Discussion

In this paper, we have computed the DM density produced by bremsstrahlung off
an inflaton decay product. We evaluated the rate of such production in general
scenarios of DM - plasma particle interaction, including gravity and vector
portal dimension-6 effective operators. The DM bremsstrahlung dominates when
the reheating temperature is low so that the plasma annihilation to DM becomes
less efficient. This significantly reduces the required dark matter mass to
achieve the observed DM abundance for low reheating temperatures.

We stress that such dark bremsstrahlung process generally occurs in the UV
freeze-in DM scenario widely considered in the literature. It requires only
two ingredients: 1) the decay of the inflaton to particles in the thermal
plasma which is a general and natural assumption in inflationary cosmology and
necessary for reheating; and 2) a feeble interaction between the DM and the
plasma particle mediated by a heavy mediator, commonly found in unification
theories and UV freeze-in DM models. Therefore, this effect needs to be
considered in almost every UV freeze-in DM scenarios. Besides the dark matter
abundance, the dark bremsstrahlung can also produce very energetic dark
radiation that could leave imprints in cosmological observations[86].

Finally, we note that the DM production can also be enhanced by the
thermalization process [44, 45, 47, 46]. A full treatment of the
thermalization process relies on the detail of the thermal plasma is beyond
the scope of this paper. We include a qualitative comparison between the DM
bremsstrahlung and the thermalization production to the appendix and conclude
that while thermalization may compete with the production from thermal plasma,
the DM bremsstrahlung always dominates the regime of low reheating
temperature.

## Acknowledgements

We would like to thank S. Clery, Y. Xu, H. He and N. Bernal for useful
discussions. This project has received support from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skodowska-Curie
Grant Agreement No 860881-HIDDeN and the IN2P3 Master Projet UCMN. The work of
K.A.O. was supported in part by DOE grant DE-SC0011842 at the University of
Minnesota. J.Z. was supported in part by the NSF of China (No.11675086 and
11835005).

##  Appendix A Bremsstrahlung vs. Thermalization

In the thermalization process, the decay products ${\cal P}$ of the inflaton
condensate are still in a very energetic state and may significantly increase
the annihilation cross section that scales as $E^{n}/\Lambda^{n+2}$. In this
appendix, we briefly summarize the result for this process from the literature
and compare it to DM bremsstrahlung. Ref. [47] considered the annihilation of
two energetic(hard) particles from the inflaton decay into a pair of DM before
the thermalization time scale $t_{{\rm
th}}{\rm\approx\alpha^{-16/5}}m_{\phi}^{4/5}\Gamma_{\phi}^{-3/5}{M}_{P}^{-6/5}$.
The DM density generated by such a process by the end of reheating can be
recast as

$n_{{\rm th},0}=c_{{\rm
th}}\frac{T_{\mathrm{RH}}^{34/5}m_{\phi}^{n-6/5}}{{M}_{P}^{3/5}\Lambda^{n+2}}\,.$
(A.1)

For a rough order of magnitude estimate, the numerical factor $c_{{\rm
th}}\lesssim{\cal O}(10^{9-n})$ for $g_{*}\approx 100$. The relative
importance between the DM bremsstrahlung and this pre-thermalization
production scales as $m_{\phi}^{11/5}{M}_{P}^{3/5}T_{\mathrm{RH}}^{-14/5}$.
Thus, the DM bremsstrahlung always dominates at reheating temperature below
$m_{\phi}$.

Throughout the reheating process, the hard(energetic) decay products from
inflaton decay can also annihilate with particles in the plasma into a pair of
DM. Ref. [46] found that for $m_{\cal X}^{2}\lesssim T_{\mathrm{RH}}m_{\phi}$,
the relative importance of such a hard-plasma annihilation process to the pre-
thermalization hard-hard annihilation is

$\frac{n_{\rm th}^{\prime}}{n_{{\rm
th},0}}\sim\begin{cases}\left(\frac{T_{{\rm
th}}}{T_{\mathrm{RH}}}\right)^{3/2}\left(\frac{m_{\phi}}{T_{\mathrm{RH}}}\right)^{1/2}\,,&\text{for
}n=0\,,\\\ \left(\frac{T_{{\rm
th}}}{T_{\mathrm{RH}}}\right)^{3/2}\left(\frac{m_{\phi}}{T_{\mathrm{RH}}}\right)^{(2-n)/2}\,,&\text{for
}0<n\leq 5\,,\\\ \left(\frac{T_{\rm
th}}{m_{\phi}}\right)^{(n-2)/2}\,,&\text{for }n>5\,,\end{cases}$ (A.2)

where

$T_{{\rm
th}}\equiv\alpha_{0}^{4/5}\left(\frac{\Gamma_{\phi}{M}_{P}^{2}}{m_{\phi}^{3}}\right)^{2/5}m_{\phi}$
(A.3)

is the thermalization temperature and $\alpha_{0}\sim{\cal O}(10^{-2})$ the
fine structure constant of the interaction that thermalizes the inflaton decay
product. Note that for $m_{\cal X}\gtrsim T_{\mathrm{RH}}$, the production is
further suppressed by a factor of $(T_{\mathrm{RH}}m_{\phi}/m_{\cal
X}^{2})^{3/2}$ because of the production threshold. We stress that the
discussion of thermalization here is an order of magnitude estimation at best.

In Fig. 5, we show the parameter space similar to Fig. 4, but now including
the contribution from DM produced during thermalization estimated by Eq.
(A.2). We see that the bremsstrahlung dominant regions are shifted towards
lower $T_{\mathrm{RH}}$ compare to Fig. 4. The effect is most prominent at
$m_{\cal X}>T_{\mathrm{RH}}$ where the production is suppressed in plasma
annihilation but not during the thermalization era due to the energy stored in
the pre-thermalized particles. The interpretation of these figure needs some
caution. The enhancement factor $(T_{\rm th}/T_{\mathrm{RH}})^{\frac{3}{2}}$
for dimension-6 interactions in Eq. (A.2) taken from Ref. [46] is only precise
up to power counting. It ignored probable numerical factors from the phase
space and the number of degrees of freedom. Actually, in the blue shaded
region with $T_{\mathrm{RH}}>m_{\cal X}$ in Fig. 5, we get opposite result
from Ref. [46] regarding the relative importance between plasma annihilation
and thermalization production since we have considered various numerical
factors in both $n_{\rm th,0}$ and the plasma production according to [47]
while [46] did not. These factors are enough to modify the result by a few
orders of magnitudes. Therefore, at this stage we are unclear whether the
$T_{\mathrm{RH}}>m_{\cal X}$ parts of the blue shaded regions in Fig. 5 is
indeed dominated by thermalization or plasma annihilation. What is certain is
that deep in the white region the dominant production mechanism must be DM
bremsstrahlung. This is because the relative importance between bremsstrahlung
and thermalization at $m_{\cal X}<T_{\mathrm{RH}}$ scales as
$T_{\mathrm{RH}}^{-\frac{13}{10}}$. Bremsstrahlung always dominates by orders
of magnitudes at low reheating temperature.

Figure 5:  Parameter space that realizes the observed dark matter relic
density for dimension-6 operators with contribution from thermalization
production of dark matter. The upper and lower panels are for ${\cal
L}_{f\chi}^{V}$ and ${\cal L}_{f{\bf X}}^{V}$, respectively. The left panels
and right panels corresponds to $\rho_{{\rm end}}=(10^{16}{\rm GeV})^{4}$ and
$(10^{14}{\rm GeV})^{4}$. The black contours show the required effective
interaction scale, $M/g_{f{\cal X}}$ in GeV; The red line presents the
relative importance $n_{{\cal X}}^{{\rm brem}}/(n_{{\cal X}}^{{\rm other}})$
of bremsstrahlung vs other production mechanisms combined for the choice of
$M/g_{f{\cal X}}$ that achieves a realistic DM relic. The gray shaded region
is where inflaton annihilation through gravity overproduces the DM. The blue
shaded areas are dominated by thermalization production.

## References

  * [1] LUX-ZEPLIN Collaboration, J. Aalbers et al., “First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment,” arXiv:2207.03764 [hep-ex].
  * [2] PandaX-4T Collaboration, Y. Meng et al., “Dark Matter Search Results from the PandaX-4T Commissioning Run,” Phys. Rev. Lett. 127 no. 26, (2021) 261802, arXiv:2107.13438 [hep-ex].
  * [3] XENON Collaboration, E. Aprile et al., “Dark Matter Search Results from a One Ton-Year Exposure of XENON1T,” Phys. Rev. Lett. 121 no. 11, (2018) 111302, arXiv:1805.12562 [astro-ph.CO].
  * [4] LUX Collaboration, D. S. Akerib et al., “Results from a search for dark matter in the complete LUX exposure,” Phys. Rev. Lett. 118 no. 2, (2017) 021303, arXiv:1608.07648 [astro-ph.CO].
  * [5] DarkSide Collaboration, P. Agnes et al., “Low-Mass Dark Matter Search with the DarkSide-50 Experiment,” Phys. Rev. Lett. 121 no. 8, (2018) 081307, arXiv:1802.06994 [astro-ph.HE].
  * [6] CRESST Collaboration, A. H. Abdelhameed et al., “First results from the CRESST-III low-mass dark matter program,” Phys. Rev. D 100 no. 10, (2019) 102002, arXiv:1904.00498 [astro-ph.CO].
  * [7] M. Srednicki, R. Watkins, and K. A. Olive, “Calculations of Relic Densities in the Early Universe,” Nucl. Phys. B 310 (1988) 693.
  * [8] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, “The waning of the WIMP? A review of models, searches, and constraints,” Eur. Phys. J. C 78 no. 3, (2018) 203, arXiv:1703.07364 [hep-ph].
  * [9] Y. Mambrini, “Particles in the dark Universe,” Springer Ed., ISBN 978-3-030-78139-2 (2021) .
  * [10] J. R. Ellis, J. S. Hagelin, D. V. Nanopoulos, K. A. Olive, and M. Srednicki, “Supersymmetric Relics from the Big Bang,” Nucl. Phys. B 238 (1984) 453–476.
  * [11] D. V. Nanopoulos, K. A. Olive, and M. Srednicki, “After Primordial Inflation,” Phys. Lett. B 127 (1983) 30–34.
  * [12] M. Y. Khlopov and A. D. Linde, “Is It Easy to Save the Gravitino?,” Phys. Lett. B 138 (1984) 265–268.
  * [13] K. A. Olive, D. N. Schramm, and M. Srednicki, “Gravitinos as the Cold Dark Matter in an $\Omega$ = 1 Universe,” Nucl. Phys. B 255 (1985) 495–504.
  * [14] L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, “Freeze-In Production of FIMP Dark Matter,” JHEP 03 (2010) 080, arXiv:0911.1120 [hep-ph].
  * [15] Y. Mambrini, K. A. Olive, J. Quevillon, and B. Zaldivar, “Gauge Coupling Unification and Nonequilibrium Thermal Dark Matter,” Phys. Rev. Lett. 110 no. 24, (2013) 241306, arXiv:1302.4438 [hep-ph].
  * [16] N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, “The Dawn of FIMP Dark Matter: A Review of Models and Constraints,” Int. J. Mod. Phys. A 32 no. 27, (2017) 1730023, arXiv:1706.07442 [hep-ph].
  * [17] B. Barman, D. Borah, and R. Roshan, “Effective Theory of Freeze-in Dark Matter,” JCAP 11 (2020) 021, arXiv:2007.08768 [hep-ph].
  * [18] B. Barman, N. Bernal, A. Das, and R. Roshan, “Non-minimally coupled vector boson dark matter,” JCAP 01 no. 01, (2022) 047, arXiv:2108.13447 [hep-ph].
  * [19] Y. Tang and Y.-L. Wu, “On Thermal Gravitational Contribution to Particle Production and Dark Matter,” Phys. Lett. B 774 (2017) 676–681, arXiv:1708.05138 [hep-ph].
  * [20] Y. Ema, R. Jinno, K. Mukaida, and K. Nakayama, “Gravitational Effects on Inflaton Decay,” JCAP 05 (2015) 038, arXiv:1502.02475 [hep-ph].
  * [21] Y. Ema, R. Jinno, K. Mukaida, and K. Nakayama, “Gravitational particle production in oscillating backgrounds and its cosmological implications,” Phys. Rev. D 94 no. 6, (2016) 063517, arXiv:1604.08898 [hep-ph].
  * [22] Y. Ema, K. Nakayama, and Y. Tang, “Production of Purely Gravitational Dark Matter,” JHEP 09 (2018) 135, arXiv:1804.07471 [hep-ph].
  * [23] Y. Ema, K. Nakayama, and Y. Tang, “Production of purely gravitational dark matter: the case of fermion and vector boson,” JHEP 07 (2019) 060, arXiv:1903.10973 [hep-ph].
  * [24] M. Garny, M. Sandora, and M. S. Sloth, “Planckian Interacting Massive Particles as Dark Matter,” Phys. Rev. Lett. 116 no. 10, (2016) 101302, arXiv:1511.03278 [hep-ph].
  * [25] M. Garny, A. Palessandro, M. Sandora, and M. S. Sloth, “Theory and Phenomenology of Planckian Interacting Massive Particles as Dark Matter,” JCAP 02 (2018) 027, arXiv:1709.09688 [hep-ph].
  * [26] M. Chianese, B. Fu, and S. F. King, “Impact of Higgs portal on gravity-mediated production of superheavy dark matter,” JCAP 06 (2020) 019, arXiv:2003.07366 [hep-ph].
  * [27] M. Chianese, B. Fu, and S. F. King, “Interplay between neutrino and gravity portals for FIMP dark matter,” JCAP 01 (2021) 034, arXiv:2009.01847 [hep-ph].
  * [28] Y. Mambrini and K. A. Olive, “Gravitational Production of Dark Matter during Reheating,” Phys. Rev. D 103 no. 11, (2021) 115009, arXiv:2102.06214 [hep-ph].
  * [29] S. Clery, Y. Mambrini, K. A. Olive, and S. Verner, “Gravitational portals in the early Universe,” Phys. Rev. D 105 no. 7, (2022) 075005, arXiv:2112.15214 [hep-ph].
  * [30] S. Clery, Y. Mambrini, K. A. Olive, A. Shkerin, and S. Verner, “Gravitational portals with nonminimal couplings,” Phys. Rev. D 105 no. 9, (2022) 095042, arXiv:2203.02004 [hep-ph].
  * [31] B. Barman and N. Bernal, “Gravitational SIMPs,” JCAP 06 (2021) 011, arXiv:2104.10699 [hep-ph].
  * [32] M. R. Haque and D. Maity, “Gravitational dark matter: Free streaming and phase space distribution,” Phys. Rev. D 106 no. 2, (2022) 023506, arXiv:2112.14668 [hep-ph].
  * [33] S. Aoki, H. M. Lee, A. G. Menkara, and K. Yamashita, “Reheating and dark matter freeze-in in the Higgs-R2 inflation model,” JHEP 05 (2022) 121, arXiv:2202.13063 [hep-ph].
  * [34] Y. Mambrini, N. Nagata, K. A. Olive, J. Quevillon, and J. Zheng, “Dark matter and gauge coupling unification in nonsupersymmetric SO(10) grand unified models,” Phys. Rev. D 91 no. 9, (2015) 095010, arXiv:1502.06929 [hep-ph].
  * [35] N. Nagata, K. A. Olive, and J. Zheng, “Weakly-Interacting Massive Particles in Non-supersymmetric SO(10) Grand Unified Models,” JHEP 10 (2015) 193, arXiv:1509.00809 [hep-ph].
  * [36] G. Bhattacharyya, M. Dutra, Y. Mambrini, and M. Pierre, “Freezing-in dark matter through a heavy invisible Z’,” Phys. Rev. D 98 no. 3, (2018) 035038, arXiv:1806.00016 [hep-ph].
  * [37] D. Chowdhury, E. Dudas, M. Dutra, and Y. Mambrini, “Moduli Portal Dark Matter,” Phys. Rev. D 99 no. 9, (2019) 095028, arXiv:1811.01947 [hep-ph].
  * [38] E. Dudas, Y. Mambrini, and K. A. Olive, “Monochromatic neutrinos generated by dark matter and the seesaw mechanism,” Phys. Rev. D 91 (2015) 075001, arXiv:1412.3459 [hep-ph].
  * [39] E. Dudas, L. Heurtier, Y. Mambrini, K. A. Olive, and M. Pierre, “Model of metastable EeV dark matter,” Phys. Rev. D 101 no. 11, (2020) 115029, arXiv:2003.02846 [hep-ph].
  * [40] L. Heurtier and F. Huang, “Inflaton portal to a highly decoupled EeV dark matter particle,” Phys. Rev. D 100 no. 4, (2019) 043507, arXiv:1905.05191 [hep-ph].
  * [41] L. Heurtier, “The Inflaton Portal to Dark Matter,” JHEP 12 (2017) 072, arXiv:1707.08999 [hep-ph].
  * [42] K. Kaneta, Y. Mambrini, and K. A. Olive, “Radiative production of nonthermal dark matter,” Phys. Rev. D 99 no. 6, (2019) 063508, arXiv:1901.04449 [hep-ph].
  * [43] S. Davidson and S. Sarkar, “Thermalization after inflation,” JHEP 11 (2000) 012, arXiv:hep-ph/0009078.
  * [44] K. Harigaya and K. Mukaida, “Thermalization after/during Reheating,” JHEP 05 (2014) 006, arXiv:1312.3097 [hep-ph].
  * [45] K. Harigaya, M. Kawasaki, K. Mukaida, and M. Yamada, “Dark Matter Production in Late Time Reheating,” Phys. Rev. D 89 no. 8, (2014) 083532, arXiv:1402.2846 [hep-ph].
  * [46] K. Harigaya, K. Mukaida, and M. Yamada, “Dark Matter Production during the Thermalization Era,” JHEP 07 (2019) 059, arXiv:1901.11027 [hep-ph].
  * [47] M. A. G. Garcia and M. A. Amin, “Prethermalization production of dark matter,” Phys. Rev. D 98 no. 10, (2018) 103504, arXiv:1806.01865 [hep-ph].
  * [48] M. Drees and B. Najjari, “Energy spectrum of thermalizing high energy decay products in the early universe,” JCAP 10 (2021) 009, arXiv:2105.01935 [hep-ph].
  * [49] M. Drees and B. Najjari, “Multi-Species Thermalization Cascade of Energetic Particles in the Early Universe,” arXiv:2205.07741 [hep-ph].
  * [50] M. A. G. Garcia, K. Kaneta, Y. Mambrini, K. A. Olive, and S. Verner, “Freeze-in from preheating,” JCAP 03 no. 03, (2022) 016, arXiv:2109.13280 [hep-ph].
  * [51] D. J. H. Chung, E. W. Kolb, and A. Riotto, “Production of massive particles during reheating,” Phys. Rev. D 60 (1999) 063504, arXiv:hep-ph/9809453.
  * [52] G. F. Giudice, E. W. Kolb, and A. Riotto, “Largest temperature of the radiation era and its cosmological implications,” Phys. Rev. D 64 (2001) 023508, arXiv:hep-ph/0005123.
  * [53] M. A. G. Garcia, Y. Mambrini, K. A. Olive, and M. Peloso, “Enhancement of the Dark Matter Abundance Before Reheating: Applications to Gravitino Dark Matter,” Phys. Rev. D 96 no. 10, (2017) 103510, arXiv:1709.01549 [hep-ph].
  * [54] M. A. G. Garcia, K. Kaneta, Y. Mambrini, and K. A. Olive, “Inflaton Oscillations and Post-Inflationary Reheating,” JCAP 04 (2021) 012, arXiv:2012.10756 [hep-ph].
  * [55] M. A. G. Garcia, K. Kaneta, Y. Mambrini, and K. A. Olive, “Reheating and Post-inflationary Production of Dark Matter,” Phys. Rev. D 101 no. 12, (2020) 123507, arXiv:2004.08404 [hep-ph].
  * [56] B. Barman, N. Bernal, Y. Xu, and O. Zapata, “Ultraviolet freeze-in with a time-dependent inflaton decay,” JCAP 07 no. 07, (2022) 019, arXiv:2202.12906 [hep-ph].
  * [57] A. Ahmed, B. Grzadkowski, and A. Socha, “Higgs boson induced reheating and ultraviolet frozen-in dark matter,” arXiv:2207.11218 [hep-ph].
  * [58] A. Ahmed, B. Grzadkowski, and A. Socha, “Implications of time-dependent inflaton decay on reheating and dark matter production,” Phys. Lett. B 831 (2022) 137201, arXiv:2111.06065 [hep-ph].
  * [59] E. Dudas, T. Gherghetta, Y. Mambrini, and K. A. Olive, “Inflation and High-Scale Supersymmetry with an EeV Gravitino,” Phys. Rev. D 96 no. 11, (2017) 115032, arXiv:1710.07341 [hep-ph].
  * [60] E. Dudas, T. Gherghetta, K. Kaneta, Y. Mambrini, and K. A. Olive, “Gravitino decay in high scale supersymmetry with R -parity violation,” Phys. Rev. D 98 no. 1, (2018) 015030, arXiv:1805.07342 [hep-ph].
  * [61] J. Ellis, M. A. G. Garcia, D. V. Nanopoulos, K. A. Olive, and M. Peloso, “Post-Inflationary Gravitino Production Revisited,” JCAP 03 (2016) 008, arXiv:1512.05701 [astro-ph.CO].
  * [62] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [63] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer, and T. Stelzer, “MadGraph 5 : Going Beyond,” JHEP 06 (2011) 128, arXiv:1106.0522 [hep-ph].
  * [64] P. de Aquino, W. Link, F. Maltoni, O. Mattelaer, and T. Stelzer, “ALOHA: Automatic Libraries Of Helicity Amplitudes for Feynman Diagram Computations,” Comput. Phys. Commun. 183 (2012) 2254–2263, arXiv:1108.2041 [hep-ph].
  * [65] H. van Dam and M. J. G. Veltman, “Massive and massless Yang-Mills and gravitational fields,” Nucl. Phys. B 22 (1970) 397–411.
  * [66] V. I. Zakharov, “Linearized gravitation theory and the graviton mass,” JETP Lett. 12 (1970) 312.
  * [67] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, “FeynRules 2.0 - A complete toolbox for tree-level phenomenology,” Comput. Phys. Commun. 185 (2014) 2250–2300, arXiv:1310.1921 [hep-ph].
  * [68] S. Weinberg, “Infrared photons and gravitons,” Phys. Rev. 140 (1965) B516–B524.
  * [69] M. A. G. Garcia, M. Pierre, and S. Verner, “Scalar Dark Matter Production from Preheating and Structure Formation Constraints,” arXiv:2206.08940 [hep-ph].
  * [70] K. Kaneta, S. M. Lee, and K.-y. Oda, “Boltzmann or Bogoliubov? Approaches Compared in Gravitational Particle Production,” arXiv:2206.10929 [astro-ph.CO].
  * [71] M. R. Haque and D. Maity, “Gravitational Reheating,” arXiv:2201.02348 [hep-ph].
  * [72] R. T. Co, Y. Mambrini, and K. A. Olive, “Inflationary Gravitational Leptogenesis,” arXiv:2205.01689 [hep-ph].
  * [73] H. Fritzsch and P. Minkowski, “Unified Interactions of Leptons and Hadrons,” Annals Phys. 93 (1975) 193–266.
  * [74] H. Georgi, “The State of the Art—Gauge Theories,” AIP Conf. Proc. 23 (1975) 575–582.
  * [75] I. Z. Rothstein, K. S. Babu, and D. Seckel, “Planck scale symmetry breaking and majoron physics,” Nucl. Phys. B 403 (1993) 725–748, arXiv:hep-ph/9301213.
  * [76] T. Yanagida, “Horizontal Symmetry and Masses of Neutrinos,” Prog. Theor. Phys. 64 (1980) 1103.
  * [77] P. Minkowski, “$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?,” Phys. Lett. B 67 (1977) 421–428.
  * [78] M. Gell-Mann, P. Ramond, and R. Slansky, “Complex Spinors and Unified Theories,” Conf. Proc. C 790927 (1979) 315–321, arXiv:1306.4669 [hep-th].
  * [79] R. N. Mohapatra and G. Senjanovic, “Neutrino Mass and Spontaneous Parity Nonconservation,” Phys. Rev. Lett. 44 (1980) 912.
  * [80] J. Schechter and J. W. F. Valle, “Neutrino Masses in SU(2) x U(1) Theories,” Phys. Rev. D 22 (1980) 2227.
  * [81] J. Schechter and J. W. F. Valle, “Neutrino Decay and Spontaneous Violation of Lepton Number,” Phys. Rev. D 25 (1982) 774.
  * [82] L. M. Krauss and F. Wilczek, “Discrete Gauge Symmetry in Continuum Theories,” Phys. Rev. Lett. 62 (1989) 1221.
  * [83] L. E. Ibanez and G. G. Ross, “Discrete gauge symmetry anomalies,” Phys. Lett. B 260 (1991) 291–295.
  * [84] L. E. Ibanez and G. G. Ross, “Discrete gauge symmetries and the origin of baryon and lepton number conservation in supersymmetric versions of the standard model,” Nucl. Phys. B 368 (1992) 3–37.
  * [85] S. P. Martin, “Some simple criteria for gauged R-parity,” Phys. Rev. D 46 (1992) R2769–R2772, arXiv:hep-ph/9207218.
  * [86] G. Ballesteros, M. A. G. Garcia, and M. Pierre, “How warm are non-thermal relics? Lyman-$\alpha$ bounds on out-of-equilibrium dark matter,” JCAP 03 (2021) 101, arXiv:2011.13458 [hep-ph].