# WIMPs during Reheating

Nicolás Bernal  and Yong Xu

###### Abstract

Weakly Interacting Massive Particles (WIMPs) are among the best-motivated dark
matter candidates. In the standard scenario where the freeze-out occurs well
after the end of inflationary reheating, they are in tension with the severe
experimental constraints. Here, we investigate the thermal freeze-out of WIMPs
occurring during reheating, while the inflaton $\phi$ coherently oscillates in
a generic potential $\propto\phi^{n}$. Depending on the value of $n$ and the
spin of the inflaton decaying products, the evolution of the radiation and
inflaton energy densities can show distinct features, therefore, having a
considerable impact on the freeze-out behavior of WIMPs. As a result of the
injection of entropy during reheating, the parameter space compatible with the
observed DM relic abundance is enlarged. In particular, the WIMP thermally
averaged annihilation cross-section can be several magnitudes lower than that
in the standard case. Finally, we discuss the current bounds from dark matter
indirect detection experiments, and explore future challenges and
opportunities.

PI/UAN-2022-722FT

##  1 Introduction

The existence of non-baryonic dark matter (DM) in the Universe is very
compelling, as suggested by both astrophysical and cosmological observations
[1]. A viable DM candidate has to satisfy several properties: It has to be
electromagnetically neutral, stable at the cosmological scales, and non-
relativistic at the matter-radiation equality in order to allow structure
formation. Finally, it must feature a relic density
$\Omega_{\text{dm}}h^{2}\simeq 0.12$, accounting for 27% of the total energy
budget of the Universe [2, 3].

The most prominent DM production mechanism in the early universe corresponds
to the weakly interacting massive particle (WIMP) paradigm [4, 5]. In the
standard WIMP scenario, DM has a mass at the electroweak scale and couples to
the standard model (SM) thermal plasma with a sizable strength, typical of
electroweak interactions. WIMPs reach thermal equilibrium with the SM thermal
plasma and eventually freeze out, giving rise to the observed DM relic
abundance. This DM freeze-out is commonly assumed to occur well after the end
of reheating, when the Universe energy density is dominated by SM radiation.
To match observations, a thermally averaged annihilation cross-section
$\langle\sigma v\rangle=\text{few}\times 10^{-26}$ cm3/s is typically required
[6]. The WIMP mechanism is particularly interesting because it could be tested
in a number of complementary ways, including direct, indirect, and collider
probes. However, current null experimental results and severe constraints on
the natural parameter space motivate some quests beyond the standard WIMP
paradigm [4].

There are several possibilities to evade the strong experimental constraints.
For example, DM might have feeble interactions with the SM thermal plasma so
that it is produced out of equilibrium by the so-called freeze-in mechanism
(FIMP) [7, 8, 9, 10, 11, 12, 13]. Alternatively, WIMPs might have decoupled
during a non-standard cosmological epoch [14], with an annihilation cross
section much smaller than the one enforced in the standard case [15]. In
general, the production of DM in scenarios with a non-standard expansion phase
has recently gained increasing interest; see, e.g., Refs. [16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
40, 41, 15, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52]. For earlier work, see
also Refs. [53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64]. Moreover,
scenarios where freeze-out occurs during reheating are also viable, and in
such cases, the freeze-out temperature $T_{\text{fo}}$ is higher than the
reheating temperature $T_{\text{rh}}$; note that the latter characterizes the
onset of the radiation-dominated era. We note that WIMP scenarios with low
reheating temperature have been widely investigated in the literature [65, 66,
60, 67, 68, 28], usually triggered by the decay of a long-lived scalar field
with an equation-of-state (EoS) parameter $\omega=0$ and a constant decay rate
$\Gamma_{\phi}$.111For studies on baryogenesis with a low reheating
temperature or during an early matter-dominated phase, see Refs. [69, 65, 70,
71, 72] and [73, 74, 75], respectively. Furthermore, primordial gravitational
wave production in scenarios with an early matter era has recently received
particular attention [76, 77, 78, 79, 80, 81, 82]. Finally, we note that this
kind of scenario in which DM reaches chemical equilibrium with the SM bath
eliminates the possible overproduction of DM during inflation [83, 84].

In this work, different from the aforementioned literature, we investigate the
phenomenology of WIMP freezing out during reheating, where the inflaton field
$\phi$ oscillates around the minimum of a generic potential $\propto\phi^{n}$.
Such potentials naturally arise in a number of inflationary scenarios like the
$\alpha$-attractor $T$-model [85, 86], or the Starobinsky model [87, 88, 89,
90]. In such a case, both the inflaton EoS and its decay rate are modified
compared to the usually assumed case (i.e. $n=2$). In particular, the EoS
depends on the shape of the potential $\omega=(n-2)/(n+2)$ [91], and
$\Gamma_{\phi}$ develops a time dependency [92, 93]. As a consequence, the
inflaton and SM energy densities exhibit a non-standard evolution during
reheating [94], which has a strong impact on WIMP decoupling.222We note that
the phenomenology of FIMPs [92, 93, 95, 96, 97] and the QCD axion DM [98] with
a time-dependent decay rate has recently been analyzed. Depending on the value
of $n$, the inflaton decay products, as well as the reheating temperature, we
find that due to the large injection of entropy during reheating,
$\langle\sigma v\rangle$ can be several orders of magnitude smaller than in
the standard case, thus evading the experimental constraints. We highlight
that some regions of the new parameter space are currently being tested by
several experiments. Furthermore, we notice that the future CTA experiment
might increase our ability to probe different scenarios.

The remainder of the paper is organized as follows. In Sect. 2 we revisit
inflationary reheating, in particular, the evolution of the inflaton and the
radiation energy densities, with an inflaton oscillating in a potential
$\propto\phi^{n}$. In Sect. 3, we investigate WIMP DM production focusing on
the case where freeze-out occurs during reheating. Our findings are summarized
in Sect. 4.

##  2 Cosmology during Reheating

We consider an inflaton $\phi$ that, after the end of inflation, oscillates at
the bottom of a monomial potential $V$ of the form

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

where $\lambda$ is a dimensionless coupling and $\Lambda$ an energy scale.
Potentials with these types of minima naturally arise in a number of
inflationary scenarios, such as the $\alpha$-attractor $T$-model [85, 86], or
the Starobinsky model [87, 88, 89, 90]. The equation of motion for the
oscillating inflaton field is given by [91]

$\ddot{\phi}+(3\,H+\Gamma_{\phi})\,\dot{\phi}+V^{\prime}(\phi)=0\,,$ (2.2)

where $H$ denotes the Hubble expansion rate and $\Gamma_{\phi}$ the inflaton
decay rate. The dots and primes correspond to derivatives with respect to time
$t$ and $\phi$, respectively. By multiplying the above expression by
$\dot{\phi}$, Eq. (2.2) can be rewritten as

$\frac{d}{dt}\left[\frac{1}{2}\dot{\phi}^{2}+V(\phi)\right]=-(3\,H+\Gamma_{\phi})\,\dot{\phi}^{2}\,.$
(2.3)

By defining the energy density and pressure of $\phi$ as
$\rho_{\phi}\equiv\frac{1}{2}\,\dot{\phi}^{2}+V(\phi)$ and
$p_{\phi}\equiv\frac{1}{2}\,\dot{\phi}^{2}-V(\phi)$, together with the
equation-of-state parameter $\omega\equiv p_{\phi}/\rho_{\phi}=(n-2)/(n+2)$
[91], Eq. (2.3) gives the evolution of $\rho_{\phi}$ as

$\frac{d\rho_{\phi}}{dt}+\frac{6\,n}{2+n}\,H\,\rho_{\phi}=-\frac{2\,n}{2+n}\,\Gamma_{\phi}\,\rho_{\phi}\,.$
(2.4)

During reheating, that is, in the range $a_{I}\ll a\ll a_{\text{rh}}$, where
$a$ is the scale factor, and $a_{I}$ and $a_{\text{rh}}$ correspond to the
scale factors at the beginning and end of the reheating, respectively, Eq.
(2.4) admits the analytical solution

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

As during reheating the Hubble expansion rate is dominated by the inflaton
energy density, it follows that

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

At the end of reheating (i.e., at $a=a_{\text{rh}}$), the inflaton and
radiation energy densities are equal,
$\rho_{R}(a_{\text{rh}})=\rho_{\phi}(a_{\text{rh}})=3\,M_{P}^{2}\,H_{\text{rh}}^{2}$.
Note that in order not to spoil BBN, the reheating temperature must satisfy
$T_{\text{rh}}\geq T_{\text{BBN}}\simeq 4$ MeV [99, 100, 101, 102, 103].

During reheating, the inflaton transfers its energy density to SM radiation,
the end of reheating corresponding to the onset of the SM radiation-dominated
era. The SM radiation energy density $\rho_{R}$ is governed by the Boltzmann
equation [93]

$\frac{d\rho_{R}}{dt}+4\,H\,\rho_{R}=+\frac{2\,n}{2+n}\,\Gamma_{\phi}\,\rho_{\phi}\,.$
(2.7)

Using Eq. (2.5), one can solve Eq. (2.7) and further obtain the solution for
the evolution of the radiation energy density during reheating

$\rho_{R}(a)\simeq\frac{2\,\sqrt{3}\,n}{2+n}\,\frac{M_{P}}{a^{4}}\int_{a_{I}}^{a}\Gamma_{\phi}(a^{\prime})\,\sqrt{\rho_{\phi}(a^{\prime})}\,a^{\prime
3}\,da^{\prime},$ (2.8)

where a general scale factor dependence of $\Gamma_{\phi}$ has been assumed.
Such a dependence can come, for example, from the inflaton mass parameter. The
effective mass $m_{\phi}$ for the inflaton field, understood as the second
derivative of its potential, is given by

$m_{\phi}^{2}\equiv\frac{d^{2}V}{d\phi^{2}}=n\,(n-1)\,\lambda\,\frac{\phi^{n-2}}{\Lambda^{n-4}}\simeq
n\,(n-1)\,\lambda^{\frac{2}{n}}\,\Lambda^{\frac{2\,(4-n)}{n}}\rho_{\phi}(a)^{\frac{n-2}{n}}.$
(2.9)

It is interesting to note that for $n\neq 2$, $m_{\phi}$ features a field
dependence that, in turn, would lead to an inflaton decay rate with a scale
factor (or time) dependence.

In the following, two different reheating scenarios will be investigated in
which the inflaton either decays to a pair of fermions or bosons via a
trilinear coupling.

###  2.1 Fermionic Reheating

In the case where the inflaton only decays into a pair of fermions $\Psi$ and
$\bar{\Psi}$ via a trilinear interaction $y\,\phi\,\bar{\Psi}\,\Psi$, with $y$
being the corresponding Yukawa coupling, the decay rate is given by

$\Gamma_{\phi}=\frac{y_{\text{eff}}^{2}}{8\pi}\,m_{\phi}\,,$ (2.10)

where the effective coupling $y_{\text{eff}}\neq y$ (for $n\neq 2$) is
obtained after averaging over oscillations [104, 93, 105]. We will return to
our treatment for effective coupling shortly. The evolution of the SM energy
density in Eq. (2.8) becomes

$\rho_{R}(a)\simeq\frac{\sqrt{3}}{8\pi}\,\frac{n\sqrt{n(n-1)}}{7-n}\,\frac{y_{\text{eff}}^{2}\,\lambda^{\frac{1}{n}}\,M_{P}\,\rho_{\phi}(a_{\text{rh}})^{\frac{n-1}{n}}}{\Lambda^{\frac{n-4}{n}}}\left(\frac{a_{\text{rh}}}{a}\right)^{\frac{6(n-1)}{2+n}}\left[1-\left(\frac{a_{I}}{a}\right)^{\frac{2(7-n)}{2+n}}\right],$
(2.11)

and, therefore, the temperature of the SM bath evolves as

$T(a)\simeq
T_{\text{rh}}\times\begin{dcases}\left(\frac{a_{\text{rh}}}{a}\right)^{\frac{3}{2}\frac{n-1}{n+2}}&\text{
for }n<7\,,\\\ \frac{a_{\text{rh}}}{a}&\text{ for }n>7\,.\end{dcases}$ (2.12)

The Hubble expansion rate (cf. Eq. (2.6)) can be rewritten as a function of
$T$ as

$H(T)\simeq
H_{\text{rh}}\begin{dcases}\left(\frac{T}{T_{\text{rh}}}\right)^{\frac{2\,n}{n-1}}&\text{
for }n<7\,,\\\ \left(\frac{T}{T_{\text{rh}}}\right)^{\frac{3\,n}{n+2}}&\text{
for }n>7\,.\end{dcases}$ (2.13)

Before closing this subsection, we note that for the case with $n=2$ where the
inflaton oscillates in a quadratic potential with an EoS parameter $\omega=0$,
the standard dependences with the scale factor $\rho_{R}(a)\propto a^{-3/2}$
and $T(a)\propto a^{-3/8}$ are reproduced.

###  2.2 Bosonic Reheating

Alternatively, if the inflaton only decays into a pair of bosons $S$ through
the interaction $\mu\,\phi\,S\,S$, with $\mu$ being a coupling with mass
dimension, the decay rate is instead

$\Gamma_{\phi}=\frac{\mu_{\text{eff}}^{2}}{8\pi\,m_{\phi}}\,,$ (2.14)

where again the effective coupling $\mu_{\text{eff}}\neq\mu$ (if $n\neq 2$)
can be obtained after averaging over oscillations. Using a procedure similar
to the previous fermionic case, one sees that the SM energy density scales as

$\rho_{R}(a)\simeq\frac{\sqrt{3}}{8\pi}\frac{1}{1+2\,n}\sqrt{\frac{n}{n-1}}\,\frac{\mu_{\text{eff}}^{2}}{\lambda^{1/n}}\,M_{P}\,\Lambda^{\frac{n-4}{n}}\,\rho_{\phi}(a_{\text{rh}})^{\frac{1}{n}}\left(\frac{a_{\text{rh}}}{a}\right)^{\frac{6}{2+n}}\left[1-\left(\frac{a_{I}}{a}\right)^{\frac{2\,(1+2n)}{2+n}}\right],$
(2.15)

with which the SM temperature evolves as

$T(a)\simeq
T_{\text{rh}}\left(\frac{a_{\text{rh}}}{a}\right)^{\frac{3}{2}\frac{1}{2+n}},$
(2.16)

and Hubble as

$H(T)\simeq H_{\text{rh}}\left(\frac{T}{T_{\text{rh}}}\right)^{2\,n},$ (2.17)

during reheating.

Figure 1: Left: Numerical evolution of the energy densities of the SM (solid
black) and $\phi$ (dashed blue), as a function of the scale factor, for
$T_{\text{rh}}=1$ GeV. Right: Corresponding evolution of the SM temperature.

In Fig. 1, the left panels depict the evolution of both the SM radiation
(solid black) and the inflaton (dashed blue) energy densities, as a function
of the scale factor, for $T_{\text{rh}}=1$ GeV. The corresponding evolution of
the SM bath temperature is shown in the right panels. We note that the bump at
$T\sim 0.1$ GeV is due to the QCD phase transition, where the latter leads to
a sudden decrease of effective relativistic degrees of freedom in the thermal
plasma. The upper panels correspond to the case $n=2$, where the effects of
bosonic and fermionic decays are equivalent. During reheating,
$\rho_{\phi}\propto a^{-3}$, $\rho_{R}\propto a^{-3/2}$, and $T\propto
a^{-3/8}$. The middle and lower panels correspond to $n=4$ with fermionic or
bosonic decay, respectively. For the first case, $\rho_{R}\propto a^{-3}$ (and
$T\propto a^{-3/4}$) while for the latter case, $\rho_{R}\propto a^{-1}$ (and
$T\propto a^{-1/4}$). Note that to numerically solve the Boltzmann equations
for $\rho_{\phi}$ and $\rho_{R}$, the values of the effective couplings
appearing in the decay rates are needed, which shall be obtained numerically
[93]. Instead, here we fix $T_{\text{rh}}$ and then calculate the effective
coupling values so that the reheating temperature matches the desired value.
Let us remind the reader that the reheating temperature has been defined as
$\rho_{R}(T_{\text{rh}})=\rho_{\phi}(T_{\text{rh}})=3\,M_{P}^{2}\,H_{\text{rh}}^{2}$.

Before proceeding, a comment on the preheating effect is necessary. With the
help of lattice simulations, it has been shown that, during reheating, the EoS
parameter $\omega\to 1/3$ for $n\gtrsim 3$ due to inflaton fragmentation and
parametric resonance [106, 107, 108, 109]. This implies that preheating
efficiently drives the Universe to a radiation-dominated phase. However, to
fully deplete the inflaton energy, perturbative decay through trilinear
couplings between the inflaton and the daughter particles is still required,
which is expected to occur in the last stage of the heating process after
inflation [107]. Note that for efficient inflaton energy transfer via
preheating, sizable couplings are usually required, which could spoil the
inflaton potential and inflationary predictions due to loop corrections.
Furthermore, in the literature, coherent oscillations of the inflaton
condensate have been found to break down [110], which could delay preheating.
It has also been shown that depending on the spin of the daughter particles,
preheating might be shut down due to the large effective vacuum expectation
value of the Higgs field [111]. In this work, we are interested in low-
reheating-temperature scenarios, so we assume small inflaton couplings to
daughter particles, in which case preheating effects are not expected to be
efficient. Thereafter, throughout this work, we will focus on the perturbative
reheating scenario.

Having understood possible cosmological histories for the background, in the
next section the evolution of the DM number density will be studied, taking
particular care to the case where DM freezes out during reheating.

##  3 Dark Matter Production

The evolution of the DM number density $n_{\text{dm}}$ is governed by the
Boltzmann equation

$\frac{dn_{\text{dm}}}{dt}+3\,H\,n_{\text{dm}}=-\langle\sigma
v\rangle\left(n_{\text{dm}}^{2}-n_{\mathrm{eq}}^{2}\right),$ (3.1)

with the Hubble expansion rate

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

$n_{\text{eq}}(T)$ corresponds to the equilibrium number density given by

$n_{\text{eq}}(T)\simeq
g\left(\frac{m\,T}{2\pi}\right)^{3/2}e^{-\frac{m}{T}}\,$ (3.3)

for non-relativistic particles, where $m$ and $g$ denote the mass and the
number of degrees of freedom of the DM field, respectively. In this work, we
focus on a temperature-independent cross section, which corresponds to an
$S$-wave dominated interaction. We note that direct production of DM out of
the inflaton decay may also be relevant. In that case, Eq. (3.1) would have an
additional term $\propto 2\,\text{Br}\,\Gamma_{\phi}\,n_{\phi}$, where Br
denotes the branching ratio of the inflaton decay into a couple of DM
particles, and $n_{\phi}$ is the inflaton number density. Here we focus on the
effect of the time dependence of $\Gamma_{\phi}$ on the evolution of the WIMP
thermal production, and therefore, in the following analysis this contribution
will be omitted.333This is typically a good assumption as long as
$\text{Br}\lesssim 10^{-4}\,m/(100~{}\text{GeV})$ [28, 15]. Finally, for a
discussion regarding DM production during the SM thermalization process, we
refer the reader to Refs. [112, 113, 114, 115, 116].

###  3.1 Freeze out after Reheating

We first revisit the standard scenario, where DM freeze-out occurs well after
reheating, in the standard radiation-dominated era. In this case, since the SM
entropy is conserved during and after the DM freeze-out, Eq. (3.1) can
conveniently be rewritten as

$\frac{dY}{d\,x}=-\frac{\langle\sigma
v\rangle\,s}{H\,x}\left(Y^{2}-Y^{2}_{\text{eq}}\right),$ (3.4)

where $x\equiv m/T$, the DM yield $Y\equiv n_{\text{dm}}/s$, and the SM
entropy density $s(T)=\frac{2\pi^{2}}{45}\,g_{\star s}\,T^{3}$, with $g_{\star
s}(T)$ the number of relativistic degrees of freedom that contribute to the SM
entropy. Taking into account that in a radiation-dominated era
$H(T)=\sqrt{\frac{\rho_{R}}{3\,M_{P}^{2}}}=\frac{\pi}{3}\sqrt{\frac{g_{\star}(T)}{10}}\,\frac{T^{2}}{M_{P}}$,
with $g_{\star}(T)$ being the number of relativistic degrees of freedom
contributing to the SM energy density, Eq. (3.4) admits the standard
approximate solution

$Y_{0}\simeq\frac{15}{2\pi\,g_{\star
s}}\sqrt{\frac{g_{\star}}{10}}\,\frac{1}{M_{P}\,T_{\text{fo}}\,\langle\sigma
v\rangle}\,,$ (3.5)

where $Y_{0}$ corresponds to the DM yield at present, long after freeze-out.
To match the whole observed DM relic density, it is required that

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

with $\rho_{c}\simeq 1.05\times 10^{-5}\,h^{2}$ GeV/cm3 being the critical
energy density, $s_{0}\simeq 2.69\times 10^{3}$ cm-3 the present entropy
density [117], and $\Omega_{\text{dm}}h^{2}\simeq 0.12$ the observed DM relic
abundance [2].

The temperature $T_{\text{fo}}$ at which the DM freeze-out occurs, or
equivalently $x_{\text{fo}}\equiv\frac{m}{T_{\text{fo}}}$, is defined by the
equality $\left.\frac{n_{\text{eq}}\,\langle\sigma
v\rangle}{H}\right|_{x_{\text{fo}}}=1$, and given by

$x_{\text{fo}}=-\frac{1}{2}\,\mathcal{W}_{-1}\left[-\frac{8\pi^{5}}{45}\frac{g_{\star}}{g^{2}}\frac{1}{\left(M_{P}\,m\,\langle\sigma
v\rangle\right)^{2}}\right],$ (3.7)

where $\mathcal{W}_{-1}$ is the $-1$ branch of the Lambert function. The
observed DM relic abundance is typically matched for cross-sections at the
electroweak scale $\langle\sigma v\rangle=\mathcal{O}\left(10^{-9}\right)$
GeV-2 which corresponds to $\mathcal{O}\left(10^{-26}\right)$ cm3/s and
$x_{\text{fo}}\sim 25$, with a small logarithmic dependence on the DM mass
[6].

###  3.2 Freeze out during Reheating

Even if the DM freeze-out is typically assumed to occur well after the end of
reheating, when the Universe is radiation dominated, this is not guaranteed.
Alternatively, DM freeze-out could occur during reheating. In this case, the
SM entropy is not conserved because of the decay of the inflaton into SM
particles. Therefore, instead of Eq. (3.4), it is more convenient to rewrite
Eq. (3.1) as

$\frac{dN}{da}=-\frac{\langle\sigma
v\rangle}{H\,a^{4}}\left(N^{2}-N_{\text{eq}}^{2}\right),$ (3.8)

where $N\equiv a^{3}\,n_{\text{dm}}$ and $N_{\text{eq}}\equiv
a^{3}\,n_{\text{eq}}$. Using Eq. (2.6) for the Hubble parameter during
reheating, Eq. (3.8) can be analytically solved as

$N(a_{\text{rh}})\simeq\frac{6}{2+n}\,\frac{H_{\text{rh}}}{\langle\sigma
v\rangle}\,a_{\text{rh}}^{3}\left(\frac{a_{\text{rh}}}{a_{\text{fo}}}\right)^{\frac{-6}{2+n}},$
(3.9)

where $a_{\text{fo}}$ and $a_{\text{rh}}$ correspond to the scale factor when
$T=T_{\text{fo}}$ and $T=T_{\text{rh}}$, respectively. After reheating,
namely, for $a>a_{\text{rh}}$, the SM entropy is conserved, and therefore
$Y_{0}\simeq Y(a_{\text{rh}})$, implying

$Y_{0}=\frac{N(a_{\text{rh}})}{s(a_{\text{rh}})\,a_{\text{rh}}^{3}}\simeq\frac{6}{2+n}\,\frac{45}{2\pi^{2}\,g_{\star
s}}\,\frac{H_{\text{rh}}}{\langle\sigma
v\rangle\,T_{\text{rh}}^{3}}\times\begin{dcases}\left(\frac{T_{\text{rh}}}{T_{\text{fo}}}\right)^{\frac{4}{n-1}}&\text{
for fermionic decay,}\\\
\left(\frac{T_{\text{rh}}}{T_{\text{fo}}}\right)^{4}&\text{ for bosonic
decay,}\end{dcases}$ (3.10)

for $n<7$ in the case of fermionic decay. This assumption will be followed
from now on. We note that, in the era after the freeze-out and before the end
of reheating, the DM yield evolves as $Y(T)\propto(a\,T)^{-3}$, which
corresponds to

$Y(T)\propto\begin{dcases}T^{\frac{7-n}{n-1}}&\text{for fermionic decay,}\\\
T^{1+2\,n}&\text{for bosonic decay.}\end{dcases}$ (3.11)

Figure 2: Examples of the evolution of the DM abundance, for different
benchmark points described in the text that fit the total observed abundance.
The solid black lines show the DM yield $Y$, whereas the dashed black lines
show the equilibrium density. Furthermore, the dotted red vertical lines
correspond to $T=T_{\text{fo}}$ and $T=T_{\text{rh}}$, and the gray horizontal
bands depict the observed DM abundance $\Omega_{\text{dm}}h^{2}\simeq 0.12$.

Examples of the evolution of DM abundance are shown in Fig. 2, for reference
points that fit the total observed abundance: the upper panel corresponds to
$n=2$, $T_{\text{rh}}=1$ GeV, $m=100$ GeV and $\langle\sigma
v\rangle=10^{-11}$ GeV-2, the lower left panel corresponds to a fermionic
decay with $n=4$, $T_{\text{rh}}=0.1$ GeV, $m=100$ GeV and $\langle\sigma
v\rangle=6\times 10^{-10}$ GeV-2, and the lower right panel corresponds to a
bosonic decay with $n=4$, $T_{\text{rh}}=1$ GeV, $m=40$ GeV and $\langle\sigma
v\rangle=5\times 10^{-11}$ GeV-2. The solid black lines depict the evolution
of the DM yield $Y$, while the dashed black lines show the equilibrium yield
$Y_{\text{eq}}$. Furthermore, the red dotted vertical lines show
$T=T_{\text{fo}}$ and $T=T_{\text{rh}}$, and the gray horizontal bands
correspond to the observed DM abundance $\Omega_{\text{dm}}h^{2}\simeq 0.12$
[2].

The evolution of the DM yield depicted in Fig. 2 is obtained by numerically
solving the Boltzmann equations (2.4), (2.7), (3.1) and the Friedmann equation
(3.2). We have checked that the analytical estimate of the asymptotic value of
$Y_{0}$ shown in Eq. (3.10) agrees well with the numerical result.
Furthermore, the evolution of the DM yield for $T_{\text{fo}}\gg T\gg
T_{\text{rh}}$ (that is, $Y\propto x^{-5}$ for $n=2$, and $Y\propto x^{-1}$ or
$Y\propto x^{-9}$ for fermionic or bosonic decays, respectively, and $n=4$)
also matches the analytical expressions of Eq. (3.11). It is worth noticing
that in the lower right panel, the DM only reaches chemical equilibrium with
the SM bath at $x\simeq 3$.

Figure 3: Parameter space that matches the whole observed DM abundance, for
$T_{\text{rh}}=T_{\text{BBN}}$ (solid lines), $10^{-1}$ GeV (dashed), $10^{0}$
GeV (dash-dotted), and $10^{1}$ GeV (dotted). The upper panel corresponds to
$n=2$, while the lower to $n=4$ and fermionic (left) or bosonic (right) decays
of the inflaton. The thick black lines correspond to the standard scenario,
where DM freeze-out in the radiation-dominated era. The red areas are
disregarded as they generate a DM underabundance ($\Omega_{\text{dm}}<0.12$),
the reheating occurs after BBN ($T_{\text{rh}}<T_{\text{BBN}}$), or DM does
not reach chemical equilibrium with the SM (‘No freeze-out’).

The contours for the reheating temperature required to match the whole
measured DM abundance are shown in Fig. 3, in the parameter space
$[m,\,\langle\sigma v\rangle]$, for bosonic and fermionic decays of the
inflaton, and different values of $n$. The thick black lines correspond to the
standard case, where DM freezes out in a radiation-dominated era, whereas the
thin black lines to $T_{\text{rh}}=T_{\text{BBN}}$, $10^{-1}$ GeV, $10^{0}$
GeV and $10^{1}$ GeV. Some constraints apply and are shown as red regions:
$i)$ Above the thick black line, corresponding to higher cross sections, DM
decouples very late and is therefore underproduced
($\Omega_{\text{dm}}h^{2}<0.12$). $ii)$ Reheating temperatures below
$T_{\text{BBN}}$ are in conflict with cosmological observations, which
corresponds to $m\lesssim\mathcal{O}(10^{-1})$ GeV. And $iii)$, very small
thermally averaged cross sections $\langle\sigma v\rangle$ could not be enough
to guarantee chemical equilibrium between the dark and visible sectors. It is
worth mentioning that in this case, labeled ‘No freeze-out’, the whole DM
relic abundance could be matched, however, the production would not correspond
to the WIMP mechanism, but rather to a FIMP scenario, with the usual strong
dependence on initial conditions. This case will not be considered in the
present analysis.

In Fig. 3, it can be seen that a DM freeze-out during reheating allows
exploring smaller cross-sections $\langle\sigma v\rangle$ compared to the
usual case where the freeze-out occurs in the standard radiation-dominated
case. For example, with $n=2$ and $T_{\text{rh}}=10$ GeV, $\langle\sigma
v\rangle$ can be as small as $\mathcal{O}(10^{-15})$ GeV-2 for a DM mass of
$\sim 7$ TeV. In a scenario with larger $T_{\text{rh}}$, even smaller
$\langle\sigma v\rangle$ are allowed for heavier WIMPs. We also note that for
a fixed $T_{\text{rh}}$, $\langle\sigma v\rangle$ decreases with increasing
$m$. This is because for heavier WIMPs, freeze-out occurs earlier (i.e., at a
larger $T_{\text{fo}}$), which implies a smaller $\langle\sigma v\rangle$.
Finally, it should be noted that with increasing $n$, the parameter space
shrinks, as depicted in Fig. 3. The main reason is that the dilution effect
becomes less prominent with larger $n$ (that is, the factor
$\left(a_{\text{rh}}/a_{\text{fo}}\right)^{\frac{-6}{2+n}}$ in Eq. (3.9)), and
hence a larger $\langle\sigma v\rangle$ is needed to compensate for this. A
larger $\langle\sigma v\rangle$ corresponds to a later freeze-out, where the
WIMP densities are smaller, and therefore a smaller dilution is required.

Figure 4: Freeze-out temperature required to match the whole observed DM
abundance, as a function of the DM mass, for $T_{\text{rh}}=T_{\text{BBN}}$
(solid lines), $10^{-1}$ GeV (dashed lines), $10^{0}$ GeV (dotted lines in
dashes) and $10^{1}$ GeV (dotted lines). The upper panel corresponds to $n=2$,
while the lower left to a fermionic decay with $n=4$, and the lower right to a
bosonic decay with $n=4$. The thick black lines represent the standard freeze-
out in radiation domination.

In the present case, DM freezes out during reheating, at a temperature given
by

$x_{\text{fo}}\simeq\begin{dcases}\frac{3+n}{2\,(1-n)}\,\mathcal{W}_{-1}\left[\frac{2\,(1-n)}{3+n}\left(\frac{g^{2}}{(2\pi)^{3}}\,\frac{m^{2}\,T_{\text{rh}}^{4}\,\langle\sigma
v\rangle^{2}}{H_{\text{rh}}^{2}}\right)^{\frac{1-n}{3+n}}\left(\frac{m}{T_{\text{rh}}}\right)^{\frac{4}{3+n}}\right]&\text{for
fermions,}\\\
\frac{3-4\,n}{2}\,\mathcal{W}_{-1}\left[\frac{2}{3-4\,n}\left(\frac{g^{2}}{(2\pi)^{3}}\,\frac{m^{2}\,T_{\text{rh}}^{4}\langle\sigma
v\rangle^{2}}{H_{\text{rh}}^{2}}\right)^{\frac{1}{3-4n}}\left(\frac{m}{T_{\text{rh}}}\right)^{\frac{4\,(1-n)}{3-4n}}\right]&\text{for
bosons.}\end{dcases}$ (3.12)

Figure 4 shows the freeze-out temperature required to match the entire
observed DM abundance, as a function of the DM mass, for
$T_{\text{rh}}=T_{\text{BBN}}$ (solid lines), $10^{-1}$ GeV (dashed lines),
$10^{0}$ GeV (dash-dotted lines) and $10^{1}$ GeV (dotted lines). The upper
panel corresponds to $n=2$, the lower left to a fermionic decay with $n=4$,
and the lower right to a bosonic decay with $n=4$. The thick black lines
represent the standard freeze-out in radiation domination; cf. Eq. (3.7).
Additionally, the lower red bands depict a relativistic freeze-out, i.e.
$x_{\text{fo}}\leq 3$. Note that for $n=2$, the freeze-out temperature can be
as low as $x_{\text{fo}}\sim 6$ due to the entropy dilution effect. For $n=2$
and $n=4$ (bosonic), the lines are cut to ensure that WIMPs thermalize; lower
values of $x_{\text{fo}}$ would correspond to the ‘No freeze-out’ regime
depicted in Fig. 3. For a fermionic decay with $n=4$, $x_{\text{fo}}$ tends to
be independent of $m$, with only a small dependence on the relativistic
degrees of freedom. Note that the temperature during reheating (for the
fermionic case) tends to feature the same scaling as the free radiation for
large values of $n$, cf. Eq. (2.12); this implies that the corresponding
freeze-out behavior becomes similar to the standard case. Actually, such a
tendency can also be seen in Fig. 4. However, for the bosonic case, due to a
different temperature scaling (cf. Eq. (2.16)), the allowed parameter space
becomes distinct from the fermionic case with an increase of $n$. Note that,
as argued earlier, in both fermonic and bosonic cases, with increasing $n$, a
larger $\langle\sigma v\rangle$ (corresponding to a later freeze-out with
lower $T_{\text{fo}}$) is needed, implying that $x_{\text{fo}}$ increases with
$n$, as depicted in Fig. 4.

Figure 5: Same as Fig. 3 with current experimental bounds (blue) and future
sensitivities (green).

Before closing this section, we comment on some present and future
possibilities of testing the current scenario using the different channels
offered by DM indirect detection experiments. In addition to the information
given in Fig. 3, Fig. 5 also shows in blue present constraints coming from:
$i)$ CMB spectral distortions from WIMP annihilation into charged particles
[118] (labeled ‘CMB’), $ii)$ the combined analysis (labeled ‘combined’) of
$\gamma$-ray data using Fermi-LAT, HAWC, H.E.S.S., MAGIC, and VERITAS by
considering WIMP annihilation to $\tau^{+}\tau^{-}$ (upper line) and
$b\bar{b}$ (lower line) [119], and $iii)$ AMS-02 measurements of antimatter in
cosmic rays, in the positron [120] and the antiproton [121] channels.444It is
interesting to recall that the bounds presented depends on the DM annihilation
channels, and suffer from large uncertainties arising from e.g. assumptions on
the DM density profile, the local DM density, and the propagation of charged
particles in the interstellar medium [122, 123, 124, 125, 126, 127, 128].
Additionally, Fig. 5 shows in green the projected sensitivity of the ground-
based CTA experiment, which could be capable of testing WIMP DM at the TeV
scale and above. The upper (lower) green dotted line shows its expected
sensitivity for WIMP annihilation to $b\bar{b}$ ($W^{+}W^{-}$) [129]. Note
that CTA might be able to probe a large fraction of the TeV-scale WIMPs in the
standard scenario, and even a chunk of the parameter space where the freeze-
out occurs during reheating, for both fermionic and bosonic decays.

##  4 Conclusions

Despite the huge experimental efforts over the past decades, the nature of
dark matter (DM) remains one of the most challenging and fundamental questions
in particle physics. In particular, scenarios where DM is a weakly interacting
massive particle (WIMP) reaching thermal equilibrium with the standard model
(SM) bath have received immense attention both theoretically and
experimentally, but unfortunately no overwhelming evidence for WIMP DM has
been found. This motivates a quest beyond the standard WIMP paradigm. An
alternative is to renounce the thermal nature of DM, as in the FIMP paradigm.
However, these kinds of scenario suffer from a strong dependence on the
largely unknown initial conditions after cosmological inflation.

In this work, we focus on another avenue where DM remains a thermal relic, but
its freeze-out (that is, its departure from chemical equilibrium) occurs
during (and not after) inflationary reheating. In particular, we investigated
a scenario where the inflaton field $\phi$ coherently oscillates at the bottom
of a generic potential $V(\phi)\propto\phi^{n}$, while decaying into SM
particles to reheat the Universe. Depending on the details of the reheating
(i.e. the shape of the potential and the coupling between the inflaton and the
SM particles), the inflaton and the SM energy density could feature a distinct
evolution. Consequently, in this case the behavior of WIMP freeze-out differs
from the case where it occurs well after reheating. Due to the dilution effect
during reheating, smaller thermally averaged cross-sections $\langle\sigma
v\rangle$ becomes compatible with the measured DM relic abundance.

In particular, if DM freezes-out after reheating, one requires $\langle\sigma
v\rangle\sim\mathcal{O}(10^{-9})$ GeV${}^{-2}\sim\mathcal{O}(10^{-26})$ cm3/s
to match the observed DM abundance, and a DM freeze-out temperature
$x_{\text{fo}}\equiv m/T_{\text{fo}}\sim 25$. Instead, if the freeze-out
occurs during reheating while the inflaton oscillates in a quadratic potential
($n=2$), if $T_{\text{rh}}\sim 1$ GeV, $\langle\sigma v\rangle$ could be of
the order $\mathcal{O}(10^{-15})$ GeV-2 for a WIMP at the TeV scale.
Interestingly, a sizable part of the parameter space is in the range of
current experiments and future proposals. For instance, for TeV-scale WIMPs,
the upcoming CTA could probe the expected cross-sections for DM decoupling
after reheating, and a portion of the plane $[m,\,\langle\sigma v\rangle]$
relevant for the case where DM decouples during reheating.

## Acknowledgments

We would like to thank Manuel Drees once again for enlightening discussions.
NB received funding from the Patrimonio Autónomo - Fondo Nacional de
Financiamiento para la Ciencia, la Tecnología y la Innovación Francisco José
de Caldas (MinCiencias - Colombia) grants 80740-465-2020 and 80740-492-2021.
NB is also funded by the Spanish FEDER/MCIU-AEI under grant FPA2017-84543-P.
This project has received funding/support from the European Union’s Horizon
2020 research and innovation program under the Marie Skłodowska-Curie grant
agreement No 860881-HIDDeN.

## References

  * [1] G. Bertone and D. Hooper, _History of dark matter_ , _Rev. Mod. Phys._ 90 (2018) 045002 [1605.04909].
  * [2] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [3] M. Drees, _Dark Matter Theory_ , _PoS_ ICHEP2018 (2019) 730 [1811.06406].
  * [4] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre et al., _The waning of the WIMP? A review of models, searches, and constraints_ , _Eur. Phys. J. C_ 78 (2018) 203 [1703.07364].
  * [5] L. Roszkowski, E.M. Sessolo and S. Trojanowski, _WIMP dark matter candidates and searches—current status and future prospects_ , _Rept. Prog. Phys._ 81 (2018) 066201 [1707.06277].
  * [6] G. Steigman, B. Dasgupta and J.F. Beacom, _Precise Relic WIMP Abundance and its Impact on Searches for Dark Matter Annihilation_ , _Phys. Rev. D_ 86 (2012) 023506 [1204.3622].
  * [7] J. McDonald, _Thermally generated gauge singlet scalars as selfinteracting dark matter_ , _Phys. Rev. Lett._ 88 (2002) 091304 [hep-ph/0106249].
  * [8] K.-Y. Choi and L. Roszkowski, _E-WIMPs_ , _AIP Conf. Proc._ 805 (2005) 30 [hep-ph/0511003].
  * [9] A. Kusenko, _Sterile neutrinos, dark matter, and the pulsar velocities in models with a Higgs singlet_ , _Phys. Rev. Lett._ 97 (2006) 241301 [hep-ph/0609081].
  * [10] K. Petraki and A. Kusenko, _Dark-matter sterile neutrinos in models with a gauge singlet in the Higgs sector_ , _Phys. Rev. D_ 77 (2008) 065014 [0711.4646].
  * [11] L.J. Hall, K. Jedamzik, J. March-Russell and S.M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080 [0911.1120].
  * [12] F. Elahi, C. Kolda and J. Unwin, _UltraViolet Freeze-in_ , _JHEP_ 03 (2015) 048 [1410.6157].
  * [13] N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen, _The Dawn of FIMP Dark Matter: A Review of Models and Constraints_ , _Int. J. Mod. Phys. A_ 32 (2017) 1730023 [1706.07442].
  * [14] R. Allahverdi et al., _The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe_ , _Open J.Astrophys._ 04 (2021)  [2006.16182].
  * [15] P. Arias, N. Bernal, A. Herrera and C. Maldonado, _Reconstructing Non-standard Cosmologies with Dark Matter_ , _JCAP_ 10 (2019) 047 [1906.04183].
  * [16] G. Kane, K. Sinha and S. Watson, _Cosmological Moduli and the Post-Inflationary Universe: A Critical Review_ , _Int. J. Mod. Phys. D_ 24 (2015) 1530022 [1502.07746].
  * [17] R.T. Co, F. D’Eramo, L.J. Hall and D. Pappadopulo, _Freeze-In Dark Matter with Displaced Signatures at Colliders_ , _JCAP_ 12 (2015) 024 [1506.07532].
  * [18] H. Davoudiasl, D. Hooper and S.D. McDermott, _Inflatable Dark Matter_ , _Phys. Rev. Lett._ 116 (2016) 031303 [1507.08660].
  * [19] L. Randall, J. Scholtz and J. Unwin, _Flooded Dark Matter and S Level Rise_ , _JHEP_ 03 (2016) 011 [1509.08477].
  * [20] A. Berlin, D. Hooper and G. Krnjaic, _PeV-Scale Dark Matter as a Thermal Relic of a Decoupled Sector_ , _Phys. Lett. B_ 760 (2016) 106 [1602.08490].
  * [21] T. Tenkanen and V. Vaskonen, _Reheating the Standard Model from a hidden sector_ , _Phys. Rev. D_ 94 (2016) 083516 [1606.00192].
  * [22] J.A. Dror, E. Kuflik and W.H. Ng, _Codecaying Dark Matter_ , _Phys. Rev. Lett._ 117 (2016) 211801 [1607.03110].
  * [23] A. Berlin, D. Hooper and G. Krnjaic, _Thermal Dark Matter From A Highly Decoupled Sector_ , _Phys. Rev. D_ 94 (2016) 095019 [1609.02555].
  * [24] F. D’Eramo, N. Fernandez and S. Profumo, _When the Universe Expands Too Fast: Relentless Dark Matter_ , _JCAP_ 05 (2017) 012 [1703.04793].
  * [25] S. Hamdan and J. Unwin, _Dark Matter Freeze-out During Matter Domination_ , _Mod. Phys. Lett. A_ 33 (2018) 1850181 [1710.03758].
  * [26] L. Visinelli, _(Non-)thermal production of WIMPs during kination_ , _Symmetry_ 10 (2018) 546 [1710.11006].
  * [27] J.A. Dror, E. Kuflik, B. Melcher and S. Watson, _Concentrated dark matter: Enhanced small-scale structure from codecaying dark matter_ , _Phys. Rev. D_ 97 (2018) 063524 [1711.04773].
  * [28] M. Drees and F. Hajkarim, _Dark Matter Production in an Early Matter Dominated Era_ , _JCAP_ 02 (2018) 057 [1711.05007].
  * [29] F. D’Eramo, N. Fernandez and S. Profumo, _Dark Matter Freeze-in Production in Fast-Expanding Universes_ , _JCAP_ 02 (2018) 046 [1712.07453].
  * [30] D. Maity and P. Saha, _Connecting CMB anisotropy and cold dark matter phenomenology via reheating_ , _Phys. Rev. D_ 98 (2018) 103525 [1801.03059].
  * [31] N. Bernal, C. Cosme and T. Tenkanen, _Phenomenology of Self-Interacting Dark Matter in a Matter-Dominated Universe_ , _Eur. Phys. J. C_ 79 (2019) 99 [1803.08064].
  * [32] E. Hardy, _Higgs portal dark matter in non-standard cosmological histories_ , _JHEP_ 06 (2018) 043 [1804.06783].
  * [33] D. Maity and P. Saha, _CMB constraints on dark matter phenomenology via reheating in Minimal plateau inflation_ , _Phys. Dark Univ._ 25 (2019) 100317 [1804.10115].
  * [34] T. Hambye, A. Strumia and D. Teresi, _Super-cool Dark Matter_ , _JHEP_ 08 (2018) 188 [1805.01473].
  * [35] N. Bernal, C. Cosme, T. Tenkanen and V. Vaskonen, _Scalar singlet dark matter in non-standard cosmologies_ , _Eur. Phys. J. C_ 79 (2019) 30 [1806.11122].
  * [36] A. Arbey, J. Ellis, F. Mahmoudi and G. Robbins, _Dark Matter Casts Light on the Early Universe_ , _JHEP_ 10 (2018) 132 [1807.00554].
  * [37] M. Drees and F. Hajkarim, _Neutralino Dark Matter in Scenarios with Early Matter Domination_ , _JHEP_ 12 (2018) 042 [1808.05706].
  * [38] A. Betancur and Ó. Zapata, _Phenomenology of doublet-triplet fermionic dark matter in nonstandard cosmology and multicomponent dark sectors_ , _Phys. Rev. D_ 98 (2018) 095003 [1809.04990].
  * [39] C. Maldonado and J. Unwin, _Establishing the Dark Matter Relic Density in an Era of Particle Decays_ , _JCAP_ 06 (2019) 037 [1902.10746].
  * [40] A. Poulin, _Dark matter freeze-out in modified cosmological scenarios_ , _Phys. Rev. D_ 100 (2019) 043022 [1905.03126].
  * [41] T. Tenkanen, _Standard model Higgs field and hidden sector cosmology_ , _Phys. Rev. D_ 100 (2019) 083515 [1905.11737].
  * [42] C. Han, _Higgsino Dark Matter in a Non-Standard History of the Universe_ , _Phys. Lett. B_ 798 (2019) 134997 [1907.09235].
  * [43] N. Bernal, F. Elahi, C. Maldonado and J. Unwin, _Ultraviolet Freeze-in and Non-Standard Cosmologies_ , _JCAP_ 11 (2019) 026 [1909.07992].
  * [44] P. Chanda, S. Hamdan and J. Unwin, _Reviving $Z$ and Higgs Mediated Dark Matter Models in Matter Dominated Freeze-out_, _JCAP_ 01 (2020) 034 [1911.02616].
  * [45] D. Mahanta and D. Borah, _TeV Scale Leptogenesis with Dark Matter in Non-standard Cosmology_ , _JCAP_ 04 (2020) 032 [1912.09726].
  * [46] G. Arcadi, S. Profumo, F.S. Queiroz and C. Siqueira, _Right-handed Neutrino Dark Matter, Neutrino Masses, and non-Standard Cosmology in a 2HDM_ , _JCAP_ 12 (2020) 030 [2007.07920].
  * [47] P. Konar, A. Mukherjee, A.K. Saha and S. Show, _A dark clue to seesaw and leptogenesis in a pseudo-Dirac singlet doublet scenario with (non)standard cosmology_ , _JHEP_ 03 (2021) 044 [2007.15608].
  * [48] D. Bhatia and S. Mukhopadhyay, _Unitarity limits on thermal dark matter in (non-)standard cosmologies_ , _JHEP_ 03 (2021) 133 [2010.09762].
  * [49] P. Arias, D. Karamitros and L. Roszkowski, _Frozen-in fermionic singlet dark matter in non-standard cosmology with a decaying fluid_ , _JCAP_ 05 (2021) 041 [2012.07202].
  * [50] B. Barman, P. Ghosh, F.S. Queiroz and A.K. Saha, _Scalar multiplet dark matter in a fast expanding Universe: Resurrection of the desert region_ , _Phys. Rev. D_ 104 (2021) 015040 [2101.10175].
  * [51] G. Arcadi, J.P. Neto, F.S. Queiroz and C. Siqueira, _Roads for right-handed neutrino dark matter: Fast expansion, standard freeze-out, and early matter domination_ , _Phys. Rev. D_ 105 (2022) 035016 [2108.11398].
  * [52] D. Mahanta and D. Borah, _WIMPy Leptogenesis in Non-Standard Cosmologies_ , 2208.11295.
  * [53] J.D. Barrow, _Massive Particles as a Probe of the Early Universe_ , _Nucl. Phys. B_ 208 (1982) 501.
  * [54] M. Kamionkowski and M.S. Turner, _Thermal Relics: Do we Know their Abundances?_ , _Phys. Rev. D_ 42 (1990) 3310.
  * [55] J. McDonald, _WIMP Densities in Decaying Particle Dominated Cosmology_ , _Phys. Rev. D_ 43 (1991) 1063.
  * [56] P. Salati, _Quintessence and the relic density of neutralinos_ , _Phys. Lett. B_ 571 (2003) 121 [astro-ph/0207396].
  * [57] D. Comelli, M. Pietroni and A. Riotto, _Dark energy and dark matter_ , _Phys. Lett. B_ 571 (2003) 115 [hep-ph/0302080].
  * [58] F. Rosati, _Quintessential enhancement of dark matter abundance_ , _Phys. Lett. B_ 570 (2003) 5 [hep-ph/0302159].
  * [59] C. Pallis, _Massive particle decay and cold dark matter abundance_ , _Astropart. Phys._ 21 (2004) 689 [hep-ph/0402033].
  * [60] G.B. Gelmini and P. Gondolo, _Neutralino with the right cold dark matter abundance in (almost) any supersymmetric model_ , _Phys. Rev. D_ 74 (2006) 023510 [hep-ph/0602230].
  * [61] G. Gelmini, P. Gondolo, A. Soldatenko and C.E. Yaguna, _The Effect of a late decaying scalar on the neutralino relic density_ , _Phys. Rev. D_ 74 (2006) 083514 [hep-ph/0605016].
  * [62] A. Arbey and F. Mahmoudi, _SUSY constraints from relic density: High sensitivity to pre-BBN expansion rate_ , _Phys. Lett. B_ 669 (2008) 46 [0803.0741].
  * [63] T. Cohen, D.E. Morrissey and A. Pierce, _Changes in Dark Matter Properties After Freeze-Out_ , _Phys. Rev. D_ 78 (2008) 111701 [0808.3994].
  * [64] A. Arbey and F. Mahmoudi, _SUSY Constraints, Relic Density, and Very Early Universe_ , _JHEP_ 05 (2010) 051 [0906.0368].
  * [65] G.F. Giudice, E.W. Kolb and A. Riotto, _Largest temperature of the radiation era and its cosmological implications_ , _Phys. Rev. D_ 64 (2001) 023508 [hep-ph/0005123].
  * [66] N. Fornengo, A. Riotto and S. Scopel, _Supersymmetric dark matter and the reheating temperature of the universe_ , _Phys. Rev. D_ 67 (2003) 023514 [hep-ph/0208072].
  * [67] M. Drees, H. Iminniyaz and M. Kakizaki, _Abundance of cosmological relics in low-temperature scenarios_ , _Phys. Rev. D_ 73 (2006) 123502 [hep-ph/0603165].
  * [68] L. Roszkowski, S. Trojanowski and K. Turzyński, _Neutralino and gravitino dark matter with low reheating temperature_ , _JHEP_ 11 (2014) 146 [1406.0012].
  * [69] S. Davidson, M. Losada and A. Riotto, _A New perspective on baryogenesis_ , _Phys. Rev. Lett._ 84 (2000) 4284 [hep-ph/0001301].
  * [70] R. Allahverdi, B. Dutta and K. Sinha, _Baryogenesis and Late-Decaying Moduli_ , _Phys. Rev. D_ 82 (2010) 035004 [1005.2804].
  * [71] A. Beniwal, M. Lewicki, J.D. Wells, M. White and A.G. Williams, _Gravitational wave, collider and dark matter signals from a scalar singlet electroweak baryogenesis_ , _JHEP_ 08 (2017) 108 [1702.06124].
  * [72] R. Allahverdi, P.S.B. Dev and B. Dutta, _A simple testable model of baryon number violation: Baryogenesis, dark matter, neutron–antineutron oscillation and collider signals_ , _Phys. Lett. B_ 779 (2018) 262 [1712.02713].
  * [73] N. Bernal and C.S. Fong, _Hot Leptogenesis from Thermal Dark Matter_ , _JCAP_ 10 (2017) 042 [1707.02988].
  * [74] S.-L. Chen, A. Dutta Banik and Z.-K. Liu, _Leptogenesis in fast expanding Universe_ , _JCAP_ 03 (2020) 009 [1912.07185].
  * [75] M. Chakraborty and S. Roy, _Baryon asymmetry and lower bound on right handed neutrino mass in fast expanding Universe: an analytical approach_ , 2208.04046.
  * [76] H. Assadullahi and D. Wands, _Gravitational waves from an early matter era_ , _Phys. Rev. D_ 79 (2009) 083511 [0901.0989].
  * [77] R. Durrer and J. Hasenkamp, _Testing Superstring Theories with Gravitational Waves_ , _Phys. Rev. D_ 84 (2011) 064027 [1105.5283].
  * [78] L. Alabidi, K. Kohri, M. Sasaki and Y. Sendouda, _Observable induced gravitational waves from an early matter phase_ , _JCAP_ 05 (2013) 033 [1303.4519].
  * [79] F. D’Eramo and K. Schmitz, _Imprint of a scalar era on the primordial spectrum of gravitational waves_ , _Phys. Rev. Research._ 1 (2019) 013010 [1904.07870].
  * [80] N. Bernal and F. Hajkarim, _Primordial Gravitational Waves in Nonstandard Cosmologies_ , _Phys. Rev. D_ 100 (2019) 063502 [1905.10410].
  * [81] D.G. Figueroa and E.H. Tanin, _Ability of LIGO and LISA to probe the equation of state of the early Universe_ , _JCAP_ 08 (2019) 011 [1905.11960].
  * [82] N. Bernal, A. Ghoshal, F. Hajkarim and G. Lambiase, _Primordial Gravitational Wave Signals in Modified Cosmologies_ , _JCAP_ 11 (2020) 051 [2008.04959].
  * [83] M.A.G. Garcia, M. Pierre and S. Verner, _Scalar Dark Matter Production from Preheating and Structure Formation Constraints_ , 2206.08940.
  * [84] K. Kaneta, S.M. Lee and K.-y. Oda, _Boltzmann or Bogoliubov? Approaches compared in gravitational particle production_ , _JCAP_ 09 (2022) 018 [2206.10929].
  * [85] R. Kallosh and A. Linde, _Universality Class in Conformal Inflation_ , _JCAP_ 07 (2013) 002 [1306.5220].
  * [86] R. Kallosh, A. Linde and D. Roest, _Superconformal Inflationary $\alpha$-Attractors_, _JHEP_ 11 (2013) 198 [1311.0472].
  * [87] A.A. Starobinsky, _A New Type of Isotropic Cosmological Models Without Singularity_ , _Phys. Lett._ 91B (1980) 99.
  * [88] A.A. Starobinsky, _Nonsingular Model of the Universe with the Quantum Gravitational de Sitter Stage and its Observational Consequences_ , in _Second Seminar on Quantum Gravity_ , pp. 103–128, 1, 1981.
  * [89] A.A. Starobinsky, _The Perturbation Spectrum Evolving from a Nonsingular Initially De-Sitter Cosmology and the Microwave Background Anisotropy_ , _Sov. Astron. Lett._ 9 (1983) 302.
  * [90] L. Kofman, A.D. Linde and A.A. Starobinsky, _Inflationary Universe Generated by the Combined Action of a Scalar Field and Gravitational Vacuum Polarization_ , _Phys. Lett. B_ 157 (1985) 361.
  * [91] M.S. Turner, _Coherent Scalar Field Oscillations in an Expanding Universe_ , _Phys. Rev. D_ 28 (1983) 1243.
  * [92] M.A.G. Garcia, K. Kaneta, Y. Mambrini and K.A. Olive, _Reheating and Post-inflationary Production of Dark Matter_ , _Phys. Rev. D_ 101 (2020) 123507 [2004.08404].
  * [93] M.A.G. Garcia, K. Kaneta, Y. Mambrini and K.A. Olive, _Inflaton Oscillations and Post-Inflationary Reheating_ , _JCAP_ 04 (2021) 012 [2012.10756].
  * [94] R.T. Co, E. Gonzalez and K. Harigaya, _Increasing Temperature toward the Completion of Reheating_ , _JCAP_ 11 (2020) 038 [2007.04328].
  * [95] A. Ahmed, B. Grzadkowski and A. Socha, _Implications of time-dependent inflaton decay on reheating and dark matter production_ , _Phys. Lett. B_ 831 (2022) 137201 [2111.06065].
  * [96] B. Barman, N. Bernal, Y. Xu and Ó. Zapata, _Ultraviolet freeze-in with a time-dependent inflaton decay_ , _JCAP_ 07 (2022) 019 [2202.12906].
  * [97] A. Ahmed, B. Grzadkowski and A. Socha, _Higgs boson induced reheating and ultraviolet frozen-in dark matter_ , 2207.11218.
  * [98] P. Arias, N. Bernal, J.K. Osiński and L. Roszkowski, _Dark Matter Axions in the Early Universe with a Period of Increasing Temperature_ , 2207.07677.
  * [99] S. Sarkar, _Big bang nucleosynthesis and physics beyond the standard model_ , _Rept. Prog. Phys._ 59 (1996) 1493 [hep-ph/9602260].
  * [100] M. Kawasaki, K. Kohri and N. Sugiyama, _MeV scale reheating temperature and thermalization of neutrino background_ , _Phys. Rev. D_ 62 (2000) 023506 [astro-ph/0002127].
  * [101] S. Hannestad, _What is the lowest possible reheating temperature?_ , _Phys. Rev. D_ 70 (2004) 043506 [astro-ph/0403291].
  * [102] F. De Bernardis, L. Pagano and A. Melchiorri, _New constraints on the reheating temperature of the universe after WMAP-5_ , _Astropart. Phys._ 30 (2008) 192.
  * [103] P.F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti, _Bounds on very low reheating scenarios after Planck_ , _Phys. Rev. D_ 92 (2015) 123534 [1511.00672].
  * [104] Y. Shtanov, J.H. Traschen and R.H. Brandenberger, _Universe reheating after inflation_ , _Phys. Rev. D_ 51 (1995) 5438 [hep-ph/9407247].
  * [105] K. Ichikawa, T. Suyama, T. Takahashi and M. Yamaguchi, _Primordial Curvature Fluctuation and Its Non-Gaussianity in Models with Modulated Reheating_ , _Phys. Rev. D_ 78 (2008) 063545 [0807.3988].
  * [106] K.D. Lozanov and M.A. Amin, _Equation of State and Duration to Radiation Domination after Inflation_ , _Phys. Rev. Lett._ 119 (2017) 061301 [1608.01213].
  * [107] D. Maity and P. Saha, _(P)reheating after minimal Plateau Inflation and constraints from CMB_ , _JCAP_ 07 (2019) 018 [1811.11173].
  * [108] P. Saha, S. Anand and L. Sriramkumar, _Accounting for the time evolution of the equation of state parameter during reheating_ , _Phys. Rev. D_ 102 (2020) 103511 [2005.01874].
  * [109] S. Antusch, D.G. Figueroa, K. Marschall and F. Torrenti, _Energy distribution and equation of state of the early Universe: matching the end of inflation and the onset of radiation domination_ , _Phys. Lett. B_ 811 (2020) 135888 [2005.07563].
  * [110] R. Easther, R. Flauger and J.B. Gilmore, _Delayed Reheating and the Breakdown of Coherent Oscillations_ , _JCAP_ 04 (2011) 027 [1003.3011].
  * [111] K. Freese, E.I. Sfakianakis, P. Stengel and L. Visinelli, _The Higgs Boson can delay Reheating after Inflation_ , _JCAP_ 05 (2018) 067 [1712.03791].
  * [112] K. Harigaya, M. Kawasaki, K. Mukaida and M. Yamada, _Dark Matter Production in Late Time Reheating_ , _Phys. Rev. D_ 89 (2014) 083532 [1402.2846].
  * [113] K. Harigaya, K. Mukaida and M. Yamada, _Dark Matter Production during the Thermalization Era_ , _JHEP_ 07 (2019) 059 [1901.11027].
  * [114] M. Drees and B. Najjari, _Energy spectrum of thermalizing high energy decay products in the early universe_ , _JCAP_ 10 (2021) 009 [2105.01935].
  * [115] M. Drees and B. Najjari, _Multi-Species Thermalization Cascade of Energetic Particles in the Early Universe_ , 2205.07741.
  * [116] K. Mukaida and M. Yamada, _Cascades of high-energy SM particles in the primordial thermal plasma_ , 2208.11708.
  * [117] Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * [118] R.K. Leane, T.R. Slatyer, J.F. Beacom and K.C.Y. Ng, _GeV-scale thermal WIMPs: Not even slightly ruled out_ , _Phys. Rev. D_ 98 (2018) 023016 [1805.10305].
  * [119] Hess, HAWC, VERITAS, MAGIC, H.E.S.S., Fermi-LAT collaboration, _Combined dark matter searches towards dwarf spheroidal galaxies with Fermi-LAT, HAWC, H.E.S.S., MAGIC, and VERITAS_ , _PoS_ ICRC2021 (2021) 528 [2108.13646].
  * [120] L. Bergstrom, T. Bringmann, I. Cholis, D. Hooper and C. Weniger, _New Limits on Dark Matter Annihilation from AMS Cosmic Ray Positron Data_ , _Phys. Rev. Lett._ 111 (2013) 171101 [1306.3983].
  * [121] F. Calore, M. Cirelli, L. Derome, Y. Genolini, D. Maurin, P. Salati et al., _AMS-02 antiprotons and dark matter: Trimmed hints and robust bounds_ , _SciPost Phys._ 12 (2022) 163 [2202.03076].
  * [122] M. Zemp, J. Diemand, M. Kuhlen, P. Madau, B. Moore, D. Potter et al., _The Graininess of Dark Matter Haloes_ , _Mon. Not. Roy. Astron. Soc._ 394 (2009) 641 [0812.2033].
  * [123] M. Pato, O. Agertz, G. Bertone, B. Moore and R. Teyssier, _Systematic uncertainties in the determination of the local dark matter density_ , _Phys. Rev. D_ 82 (2010) 023531 [1006.1322].
  * [124] N. Bernal, J.E. Forero-Romero, R. Garani and S. Palomares-Ruiz, _Systematic uncertainties from halo asphericity in dark matter searches_ , _JCAP_ 09 (2014) 004 [1405.6240].
  * [125] M. Boudaud et al., _A new look at the cosmic ray positron fraction_ , _Astron. Astrophys._ 575 (2015) A67 [1410.3799].
  * [126] N. Bernal, J.E. Forero-Romero, R. Garani and S. Palomares-Ruiz, _Systematic uncertainties from halo asphericity in dark matter searches_ , _Nucl. Part. Phys. Proc._ 267-269 (2015) 345.
  * [127] N. Bernal, L. Necib and T.R. Slatyer, _Spherical Cows in Dark Matter Indirect Detection_ , _JCAP_ 12 (2016) 030 [1606.00433].
  * [128] M. Benito, N. Bernal, N. Bozorgnia, F. Calore and F. Iocco, _Particle Dark Matter Constraints: the Effect of Galactic Uncertainties_ , _JCAP_ 02 (2017) 007 [1612.02010].
  * [129] CTA collaboration, _Sensitivity of the Cherenkov Telescope Array to a dark matter signal from the Galactic centre_ , _JCAP_ 01 (2021) 057 [2007.16129].


11institutetext: Université Paris-Saclay, CentraleSupélec, Mathématiques et
Informatique pour la Complexité et les Systèmes, Inria Saclay, 91190, Gif-sur-
Yvette, France. 11email<EMAIL_ADDRESS>22institutetext:
Université Paris-Saclay, Institut Gustave Roussy, Inserm, Radiothérapie
Moléculaire et Innovation Thérapeutique, 94800, Villejuif,
France.33institutetext: Departement de Radiology, Hôpital Cochin, AP-HP
Centre, Université de Paris, 27 Rue du Faubourg Saint-Jacques, 75014 Paris,
France44institutetext: TheraPanacea, Pépiniere Santé Cochin, Paris, France.

# Exploring Deep Registration Latent Spaces

Théo Estienne 1122 Maria Vakalopoulou 11 Stergios Christodoulidis 11 Enzo
Battistella 1122 Théophraste Henry 22 Marvin Lerousseau 1122 Amaury Leroy 2244
Guillaume Chassagnon 33 Marie-Pierre Revel 33 Nikos Paragios 44 Eric Deutsch
22

###### Abstract

Explainability of deep neural networks is one of the most challenging and
interesting problems in the field. In this study, we investigate the topic
focusing on the interpretability of deep learning-based registration methods.
In particular, with the appropriate model architecture and using a simple
linear projection, we decompose the encoding space, generating a new basis,
and we empirically show that this basis captures various decomposed
anatomically aware geometrical transformations. We perform experiments using
two different datasets focusing on lungs and hippocampus MRI. We show that
such an approach can decompose the highly convoluted latent spaces of
registration pipelines in an orthogonal space with several interesting
properties. We hope that this work could shed some light on a better
understanding of deep learning-based registration methods.

###### Keywords:

Deep Learning-based Medical Image Registration Deformable Registration
Explainability

##  1 Introduction

Deep learning methods provide the state of the art performance for various
applications currently. This is due to their inherent property to generate
highly abstract representations hierarchically. These representations are
building on top of each other, making it possible to encode highly non-linear
manifolds. Even though such hierarchies can outperform traditional methods,
they lack explainability, making their translation difficult to solve real-
life problems. This drawback is of great significance in the medical field and
especially for the algorithms that are intended to be adapted to clinical
practice, addressing problems of precision medicine [12, 2]. For these
reasons, it is essential to identify ways to understand better the high
throughput operations that are applied.

Recently, with the introduction of the differentiable spatial transformer
[13], trainable deep learning registration methods are becoming more and more
popular, reducing computational times while reporting similar to traditional
methods performance [1, 25, 17]. Meanwhile, the deformation field, which is
one of the products of deformable registration methods, has been shown to
encode not only the spatial correspondences but also clinical relevant
information that could add valuable aspects to a variety of problems related
to survival assessment or anomaly detection [18]. Indeed, encoding information
between subjects can be very informative for various medical tasks such as
medical image segmentation [9]. However, according to our knowledge, there are
not many efforts focusing on understanding and analysing this encoding
information which could initiate the explainability of deep learning-based
registration methods.

In this study, we propose a framework for interpreting the encoded
representation of deep learning-based registration methods. In particular,
with the appropriate model architecture and by using a simple linear
projection, we decompose the encoding space, generating a new basis that
captures various geometrical operations. This decomposed encoding space is
then driving the generation of the deformation field. The contributions of
this work are twofold: (i) to the best of our knowledge, this study is one of
the first to explore the explainability of deep learning-based registration
methods through their encodings using linear projections, (ii) we show
empirically, using two different datasets, one focusing on lungs and the other
on the brain hippocampus that our projections are associated with different
types of deformations and in particular rigid transformations. We hope that
this work can highlight the very challenging topic of explainability of deep
neural networks.

##  2 Related Work

Explaining how deep neural networks function is a matter of extensive research
the recent years. GradCam [21] is one of the most popular methods that can
provide some insights on deep neural networks for many applications, including
medical imaging. GradCam highlights the region of the original input that
contributes the most to the final prediction, producing coarse heatmaps based
on the gradients. Similar to GradCam, there are many additional methods based
on the gradient [27, 3, 24] that are commonly used for the explainability of
the models. Moreover, in [10] the authors proposed a general framework of
explanations as meta-predictors while they also reinterpret the network’s
saliency providing a natural generalisation of the gradient-based saliency
techniques. Even though such approaches can provide information on where the
models attend, they can be mostly utilised in classification or detection
schemes.

Representation disentangling methodologies is a concurrent field of research
also investigating explainability topics. Such approaches are mainly focusing
on generating interpretable latent representations by enforcing several
constraints. This can be achieved either using architecture tricks [22, 20] or
with appropriate loss functions [4, 14]. In medical image computing, several
studies focus on approaches for generating disentangled or decomposed
representations. In [19] for example, the authors proposed a multimodal image
registration method by decomposing the volumes into a common latent shape
space and separate latent appearance spaces via image-to-image translation
approach and generative models. Our method shares many common points with the
approaches mentioned above, yet it focuses on exploring the registration
latent space decomposition.

##  3 Methodology

Deep learning-based registration methods have received much attention in the
last few years [1, 25]. Formally, let us consider two volumes, the moving $M$
and the fixed $F$. The goal of deep learning-based registration methods is to
obtain the best parameters $\theta^{*}$ for the network $g_{\theta}$ that will
map most accurately $M$ to $F$ using the predicted deformation grid $\Phi$.
The network $g_{\theta}$ usually is composed of an encoding $E_{\psi}$ and a
decoding $D_{\omega}$ part.

[width=]tikz_schema

Figure 1: Overall overview of our proposed framework. The different subjects
($X$,$Y$) are projected on the latent representation by the encoder $E_{\psi}$
and then a linear decomposition of this latent space is calculated to identify
a new vector space $(\vec{u}_{i})$.

There are multiple ways to fuse the input volumes in deep learning-based
registration approaches. Most of the methods use an early fusion strategy on
which the two volumes are concatenated before they pass through the
$g_{\theta}$. However, some methods investigate late fusion strategies [11, 7]
where the two volumes pass independently through the encoder, and their
merging operation is achieved in the encoding representation using various
operations such as concatenation or subtraction. Thanks to this formulation,
each volume has a unique encoding representation. In this study, we adopt the
second strategy using the subtraction operation to encode each volume
independently and calculate its latent space’s linear decomposition. In Figure
1, the overall scheme is presented.

###  3.1 Deep learning-based registration scheme

To perform our experiments and obtain our embeddings, we defined a network
based on a 3D UNet architecture [5]. The encoder and the decoder are composed
of a fixed number of blocks with 3D convolution layers (stride $3$, padding
$1$), instance normalisation layer and leaky ReLU activation function. The
down and up-sampling operations are performed with a 3D convolution layer with
stride and padding of $2$. One of the main differences in our architecture was
the absence of skip connections. Indeed, we want to enforce that all
information passes through the last encoding layer without any leak due to the
skip connections. This modification led us to reduce the downsampling
operations from four to three for the lung dataset, to maintain the spatial
resolution of the bottleneck.

Different formulations have been proposed to generate the deformation from
deep learning schemes, such as displacement field formulation [1],
diffeomorphic formulations [6, 15] and formulations based on the spatial
gradients [25]. In this work, we focused on the last one, with our network
regressing the spatial gradients $\nabla_{x}\Phi_{x}$, $\nabla_{y}\Phi_{y}$
and $\nabla_{z}\Phi_{z}$, while the final deformation field is obtained
through a cumulative sum operation. We also followed the symmetric formulation
proposed in [8], predicting both the forward and backward deformations:
$\nabla\Phi_{M\rightarrow F}=D_{\omega}(E_{\psi}(M)-E_{\psi}(F))$ and
$\nabla\Phi_{F\rightarrow M}=D_{\omega}(E_{\psi}(M)-E_{\psi}(F))$.

The network was trained with a combination of four losses, one focusing on the
intensity similarity using normalised cross-correlation ($\mathcal{L}_{sim}$),
one focusing on anatomical structures using dice loss ($\mathcal{L}_{seg}$)
and two losses for regularisation of the displacements. The first one was the
Jacobian loss which is exploited on different works such as [17, 16, 26]
($\mathcal{L}_{jac}$) and the second one enforcing smooth gradients similar to
[8] ($\mathcal{L}_{smooth}$). As such our final loss is:
$\mathcal{L}=(\mathcal{L}_{sim}+\mathcal{L}_{seg}+\alpha\mathcal{L}_{smooth}+\beta\mathcal{L}_{jac})_{\small{M\rightarrow
F}}+(\mathcal{L}_{sim}+\mathcal{L}_{seg}+\alpha\mathcal{L}_{smooth}+\beta\mathcal{L}_{jac})_{F\rightarrow
M}$ with $\alpha$ and $\beta$ being the weights of the regularisation losses.

###  3.2 Decomposition of latent space

Let $A_{train}=\\{X_{i}|i\in\left[0,n\right]\\}$ be the set of our $n$
training samples. The proposed formulation apply the encoder independently to
each volume, and thus we can obtain the set of latent vectors:
$E_{\psi}(A_{train})=\\{E_{\psi}(X_{i})|i\in\left[0,n\right]\\}$. Then, we
decompose this space using principal components analysis (PCA). That way, we
obtain a set of principal vectors
$\mathcal{U}_{K}=(\overrightarrow{u_{1}},\cdots,\overrightarrow{u_{K}})$ with
$K$ being a hyperparameter fixing the number of principal components. It worth
noting that each vector $\vec{u}_{i}$ has the same size as the activation map
of the encoder’s last layer. This size depends on the number of channels, the
size of the input images and the number of downsampling operations. We flatten
each encoding representation from its four dimensions representation (channel
dimension and the three spatial dimensions) to a one-dimensional array to
perform the PCA. Thus, the PCA is not calculated channel-wise, but all the
channels are considered together. Each principal vector $\vec{u}_{i}$ can be
converted to a deformation grid $\phi_{i}$ using the corresponding decoder
$D_{\omega}$: $\phi_{i}=D_{\omega}(\vec{u}_{i})$. Therefore, we obtained a set
of elementary transformations $\\{\phi_{i}\\}_{i=1\cdots K}$. These elementary
transformations generate a basis that can be used to approximate and decompose
every new deformation. Using such a decomposition, we can obtain a
representation in small dimensions of every training volume $X_{i}$. These
representations are obtained by the projection of $E_{\psi}(X_{i})$ to each
principal vector: $a^{j}_{i}=E_{\psi}(X_{i})\cdot\overrightarrow{u_{j}}$. For
every volume of our training set we have the approximation:
$E_{\psi}(X_{i})\approx\sum_{j=1}^{K}a^{j}_{i}\overrightarrow{u_{j}}$. After
calculating the vector of the principal components $\mathcal{U}_{K}$ with the
training set, we projected each image of the validation set to obtain its PCA
representation.

###  3.3 Implementation and Training Details

The Adam optimiser was used for our training, with a constant learning rate
set to $1e^{-4}$, a batch size equal to $4$ and $8$ for lung and hippocampus,
respectively. Our models were trained for $600$ epochs, and it last
approximately $4$ and $9$ hours for the lung and hippocampus dataset.
Concerning data augmentation, we applied random flip, rotation, translation
and zoom. Moreover, the weights of the different loss components were set to
$1$ except the loss for smoothness set to $\alpha=0.1$ for both datasets and
the weight for the jacobian loss $\beta$ that was discarded for the
hippocampus dataset. During the training process, we registered random pairs
of different patients. Our training has been performed using the framework
PyTorch and one GPU card Nvidia Tesla V100 with 32G memory. The PCA
decomposition was calculated using the library scikit-learn, and the number of
principal components $K$ was set to $32$. Using $32$ components, our
decomposition covered $95\%$ and $93\%$ of the variance ratio for the lung and
hippocampus dataset, respectively, while $42\%$ and $62\%$ are covered by the
first four components for each dataset, respectively.

##  4 Experiments and Results

We performed our experiments on two different datasets, one public and one
private. Starting with the public dataset, we conduct experiments with the
hippocampus111http://medicaldecathlon.com/ [23]. This dataset comprises $394$
MRI with the segmentations of two small structures, the head and the body
hippocampus. The images have been cropped around the hippocampus into small
patches of $64\times 64\times 64$ voxels. The second dataset is composed of
$41$ lung MRI patients ($12$ healthy and $29$ diseased with pulmonary
fibrosis) together with their lung segmentations. Each patient had been
acquired in two states, the inspiration and the expiration. Each volume has
been resampled to $1.39$mm on the x and z-axis and $1.69$ on the y-axis and
cropped to $128\times 64\times 128$ volumes. The same normalisation strategy
has been applied for the two datasets: $\mathcal{N}(0,1)$ standardisation,
clip to $\left[-5,5\right]$ to remove outliers values and min-max
normalisation to $(0,1)$. Both datasets were split into training and
validation, resulting in $200$ and $60$ patients for hippocampus and $28$ and
$13$ patients for the lung dataset.

As the first step of our evaluation, we benchmarked the performance of the
registration network $g_{\theta}$, on which our decomposition is based on.
More specifically, we obtained a Dice coefficient of $0.90\pm 0.04$ for the
lungs and $0.76\pm 0.05$ for the hippocampus, while the initial unregistered
cases reported a Dice of $0.74\pm 0.14$ and $0.59\pm 0.15$ respectively.
Moreover, we calculated the registration for $g_{\theta}$ with the skip
connections to measure their impact on the registration. The Dice is then
equal to $0.92\pm 0.02$ and $0.85\pm 0.03$ respectively. Thus, by removing the
skip-connections, we decrease the performance of the registration, slightly on
the lungs, more importantly, on the hippocampus. However, both strategies
register the pair of volumes properly.

###  4.1 Qualitative Evaluation

To understand and evaluate the calculated components of $\mathcal{U}_{K}$ per
dataset, we perform a qualitative analysis. In particular, for each principal
vector $\overrightarrow{u_{i}}$, we calculated the corresponding deformation
$\\{\phi_{i}\\}_{i=1\cdots K}$ and we applied to the moving image $M$ together
with its corresponding segmentation map $M_{seg}$. More formally, the deformed
contour correspond to
$\mathcal{W}(D_{\omega}(\lambda\overrightarrow{u_{i}}),M_{seg})$ with
$\mathcal{W}$ being the warping operation and $\lambda$ the parameter to
control the strength of the displacements for better visualisation.

Figure 2: Visualisation of the displacements following the first four
principal components. For each component, we depicted coronal and sagittal
views. In red, the contours of the lungs of the $M$, and in gold, the
$\mathcal{W}(D_{\omega}(\lambda\overrightarrow{u_{i}}),M)$ lung’s contours.
The deformation field is represented with arrows. The arrows’ norm is
represented with a colour map, red being the smallest and white the largest.
Other patients and components are displayed on supplementary materials.

Figure 3: Representation of the deformed MR together with its lung contours
following the first and fourth principal components $u_{1}$, $u_{4}$. The red
contour represents the position of the lung segmentation of the input image
while the gold contour the position of the deformed lung. The values of lambda
range from: $-200$, $-100$, $0$, $100$ and $200$ (left to right). Negative
values of lambda correspond to an upward translation, while positive values to
a downward translation.

In Figure 2, we show the principal components obtained for one validation
subject for the lung dataset. Interestingly, one can observe that each
$\phi_{i}$ corresponds to a different elementary transformation. More
precisely, the $1^{st}$ component is associated with translation, the $2^{nd}$
with a deformation focusing on the bottom of the lungs, the $3^{rd}$ with a
deformation on the right lung focusing also on the heart region and lastly the
$4^{th}$ with a deformation focusing on the top region of the lung and
shoulders.

Figure 4: Visualisation of the displacements following the first four
principal components. We depicted a sagittal view of one patient of the
validation set of the hippocampus dataset. We represented the ground truth
hippocampus contours (red) and the deformed one (gold), following the
principal components.

In Figure 3, we show the effect of the values of $\lambda$. In the figure, we
present the lung contours of the scaled component (in red) and the
corresponding component of the warped of the first ad third components. As we
have indicated, the $1st$ component is associated with translation, which we
can also be observed in this visualisation. In particular, for this experiment
we sample $\lambda$ from the values $\\{-200,-100,0,100,200\\}$. One can
observe that we retrieve a near identity deformation for a value of $0$, while
for negative and positive values, the lung moves up and down, respectively. On
the other hand, the fourth component is responsible for deforming the
shoulders and the top of the lungs. Indeed, one can observe that through the
different $\lambda$ values, the top lungs region is the one that reports the
most changes. In Figure 7, similarly, the $4$ deformations produced by the
first $4$ principal components of the hippocampus dataset are presented. In
this case, the $1^{st}$ component seems to capture rotation on the sagittal
plane, the $2^{nd}$ translation and shrinking towards the bottom right, while
the $3^{rd}$ seems to be the same operation towards the top left corner.
Finally, the $4^{th}$ seems to be related to scaling, inflating both the
hippocampus’s head and tail. We observed that the decomposition of the two
datasets created different elementary transformations $\phi_{i}$, with
transformations closer to affine for the hippocampus and more complex for the
lung.

Figure 5: Visualisation of the differences between the components of a
reference image and the same image to which we applied a predefined
transformation. The first ten components have been displayed. From right to
left: rotation, translation along Z-axis and scaling.

Finally, to verify the obtained decomposition, we performed a case study for
all the validation subjects of the hippocampus dataset. More specifically, we
applied some predefined translation using $10$ pixels on the $z$ axis,
rotation using $20$ degrees on the $z$ axis and scaling using a factor of
$0.2$, transforming each subject $X$ to $X^{\prime}$. Then we calculated the
difference between the projection of $E_{\psi}(X)$ and $E_{\psi}(X^{\prime})$
on the PCA decomposition. In Figure 5, a box plot for all the validation
subjects of the absolute difference is presented. Specifically, the amount
$||a^{j}_{E_{\psi}(X)}-a^{j}_{E_{\psi}(X^{\prime})}||$ is shown for each
principal component $j$, with $a^{j}_{E_{\psi}(X)}$ being the projection of
$E_{\psi}(X)$ on the principal vectors $\mathcal{U}_{K}$, for the three
different applied deformations. One can observe that for rotation and
translation, only one component is significantly different from the rest. In
the case of scaling, however, two components seem to be more activated.
Moreover, these findings are in accordance with Figure 7 for the rotation and
translation. In supplementary materials, we upgraded the Figure 5 by comparing
the network with and without skip-connections. Contrary to our proposed
formulation, many components are activated with the skip-connections,
demonstrating the necessity of removing them to have a good decomposition.

##  5 Discussion & Conclusion

In this work, we proposed an approach to decompose and explain the
representations of deep learning-based registration methods. The proposed
method utilises a linear decomposition on the latent space projecting it to
principal components closely associated with anatomically aware deformations.
Our method’s dynamics are demonstrated in two different MRI datasets, focusing
on lung and hippocampus anatomies. We hope that these results will take some
steps towards a better understanding of latent representations learned by the
deep learning registration architectures. We also explored a direct
application of the PCA on the deformation’s grid instead of the latent
representation. However, we did not observe any qualitative correlations with
types of deformations, which is the case for our proposed formulation. One of
the main limitations of our approach is the difficulty of quantitative
evaluation. Our future steps include the more extensive evaluation of our
method, including new anatomies such as abdominal volumes and its clinical
significance. More specifically, we want to apply our approach to multi-
temporal follow-up of patients, monitoring diseases’ progression.

## Funding

This work has been partially funding by the ARC: Grant SIGNIT201801286, the
Fondation pour la Recherche Médicale: Grant DIC20161236437, SIRIC-SOCRATE 2.0,
ITMO Cancer, Institut National du Cancer (INCa) and Amazon Web Services (AWS).

## References

  * [1] Balakrishnan, G., Zhao, A., Sabuncu, M.R., Guttag, J., Dalca, A.V.: Voxelmorph: a learning framework for deformable medical image registration. IEEE transactions on medical imaging 38(8), 1788–1800 (2019)
  * [2] Castro, D.C., Walker, I., Glocker, B.: Causality matters in medical imaging. Nature Communications 11(1), 1–10 (2020)
  * [3] Chattopadhay, A., Sarkar, A., Howlader, P., Balasubramanian, V.N.: Grad-cam++: Generalized gradient-based visual explanations for deep convolutional networks. In: 2018 IEEE Winter Conference on Applications of Computer Vision (WACV). pp. 839–847. IEEE (2018)
  * [4] Chen, X., Duan, Y., Houthooft, R., Schulman, J., Sutskever, I., Abbeel, P.: Infogan: Interpretable representation learning by information maximizing generative adversarial nets. arXiv preprint arXiv:1606.03657 (2016)
  * [5] Çiçek, Ö., Abdulkadir, A., Lienkamp, S.S., Brox, T., Ronneberger, O.: 3d u-net: learning dense volumetric segmentation from sparse annotation. In: International conference on medical image computing and computer-assisted intervention. pp. 424–432. Springer (2016)
  * [6] Dalca, A.V., Balakrishnan, G., Guttag, J., Sabuncu, M.R.: Unsupervised learning for fast probabilistic diffeomorphic registration. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 729–738. Springer (2018)
  * [7] Estienne, T., Lerousseau, M., Vakalopoulou, M., Alvarez Andres, E., Battistella, E., Carré, A., Chandra, S., Christodoulidis, S., Sahasrabudhe, M., Sun, R., et al.: Deep learning-based concurrent brain registration and tumor segmentation. Frontiers in computational neuroscience 14, 17 (2020)
  * [8] Estienne, T., Vakalopoulou, M., Battistella, E., Carré, A., Henry, T., Lerousseau, M., Robert, C., Paragios, N., Deutsch, E.: Deep Learning Based Registration Using Spatial Gradients and Noisy Segmentation Labels. In: Segmentation, Classification, and Registration of Multi-Modality Medical Imaging Data. pp. 87–93. Lecture Notes in Computer Science, Cham (2021)
  * [9] Estienne, T., Vakalopoulou, M., Christodoulidis, S., Battistela, E., Lerousseau, M., Carre, A., Klausner, G., Sun, R., Robert, C., Mougiakakou, S., et al.: U-resnet: Ultimate coupling of registration and segmentation with deep nets. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 310–319. Springer (2019)
  * [10] Fong, R.C., Vedaldi, A.: Interpretable explanations of black boxes by meaningful perturbation. In: Proceedings of the IEEE International Conference on Computer Vision. pp. 3429–3437 (2017)
  * [11] Heinrich, M.P.: Closing the gap between deep and conventional image registration using probabilistic dense displacement networks. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 50–58. Springer (2019)
  * [12] Holzinger, A., Langs, G., Denk, H., Zatloukal, K., Müller, H.: Causability and explainability of artificial intelligence in medicine. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery 9(4), e1312 (2019)
  * [13] Jaderberg, M., Simonyan, K., Zisserman, A., Kavukcuoglu, K.: Spatial Transformer Networks. arXiv:1506.02025 [cs] (Feb 2016)
  * [14] Kingma, D.P., Welling, M.: Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114 (2013)
  * [15] Krebs, J., Delingette, H., Mailhé, B., Ayache, N., Mansi, T.: Learning a probabilistic model for diffeomorphic registration. IEEE transactions on medical imaging 38(9), 2165–2176 (2019)
  * [16] Kuang, D., Schmah, T.: FAIM – A ConvNet Method for Unsupervised 3D Medical Image Registration. arXiv:1811.09243 [cs] (Jun 2019)
  * [17] Mok, T.C., Chung, A.: Fast symmetric diffeomorphic image registration with convolutional neural networks. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 4644–4653 (2020)
  * [18] Ou, Y., Weinstein, S.P., Conant, E.F., Englander, S., Da, X., Gaonkar, B., Hsieh, M.K., Rosen, M., DeMichele, A., Davatzikos, C., et al.: Deformable registration for quantifying longitudinal tumor changes during neoadjuvant chemotherapy. Magnetic resonance in medicine 73(6), 2343–2356 (2015)
  * [19] Qin, C., Shi, B., Liao, R., Mansi, T., Rueckert, D., Kamen, A.: Unsupervised Deformable Registration for Multi-modal Images via Disentangled Representations. In: Chung, A.C.S., Gee, J.C., Yushkevich, P.A., Bao, S. (eds.) Information Processing in Medical Imaging. pp. 249–261 (2019)
  * [20] Sahasrabudhe, M., Shu, Z., Bartrum, E., Alp Guler, R., Samaras, D., Kokkinos, I.: Lifting autoencoders: Unsupervised learning of a fully-disentangled 3d morphable model using deep non-rigid structure from motion. In: Proceedings of the IEEE/CVF International Conference on Computer Vision Workshops. pp. 0–0 (2019)
  * [21] Selvaraju, R.R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., Batra, D.: Grad-cam: Visual explanations from deep networks via gradient-based localization. In: Proceedings of the IEEE international conference on computer vision. pp. 618–626 (2017)
  * [22] Shu, Z., Sahasrabudhe, M., Guler, R.A., Samaras, D., Paragios, N., Kokkinos, I.: Deforming autoencoders: Unsupervised disentangling of shape and appearance. In: Proceedings of the European conference on computer vision (ECCV). pp. 650–665 (2018)
  * [23] Simpson, A.L., Antonelli, M., Bakas, S., Bilello, M., Farahani, K., Van Ginneken, B., Kopp-Schneider, A., Landman, B.A., Litjens, G., Menze, B., et al.: A large annotated medical image dataset for the development and evaluation of segmentation algorithms. arXiv preprint arXiv:1902.09063 (2019)
  * [24] Springenberg, J.T., Dosovitskiy, A., Brox, T., Riedmiller, M.: Striving for simplicity: The all convolutional net. arXiv preprint arXiv:1412.6806 (2014)
  * [25] Stergios, C., Mihir, S., Maria, V., Guillaume, C., Marie-Pierre, R., Stavroula, M., Nikos, P.: Linear and Deformable Image Registration with 3D Convolutional Neural Networks. In: Image Analysis for Moving Organ, Breast, and Thoracic Images. pp. 13–22. Lecture Notes in Computer Science, Springer International Publishing, Cham (2018)
  * [26] Zhang, S., Liu, P.X., Zheng, M., Shi, W.: A diffeomorphic unsupervised method for deformable soft tissue image registration. Computers in Biology and Medicine 120, 103708 (May 2020). https://doi.org/10.1016/j.compbiomed.2020.103708
  * [27] Zhou, B., Khosla, A., Lapedriza, A., Oliva, A., Torralba, A.: Learning deep features for discriminative localization. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 2921–2929 (2016)

## Supplementary Material

Figure 6: Extension of the figure 2 with two others patients and the component
1 to 4.

Figure 7: Extension of the Figure 4 with two other patients and the components
1 to 8 (first row 1-4, second row 5-8).

Figure 8: Comparison between the networks with and without skip connections.
We displayed the difference between the components of an image and the same
image to which we applied a predefined transformation. Only one or two
components are activated without the skip-connections, while many of them are
with the skip. The first ten components have been displayed. The results are
in blue and pink for respectively without and with the skip-connections. From
right to left: rotation, translation along Z-axis and scaling.