# Two Real Scalar WIMP Model in the Assisted Freeze-Out Scenario

Bastián Díaz Sáez<EMAIL_ADDRESS>Physik-Department, Technische
Universität München, James-Franck-Straße, 85748 Garching, Germany  Kilian
Möhling<EMAIL_ADDRESS>Institut für Kern- und Teilchenphysik,
TU Dresden, Zellescher Weg 19, 01069 Dresden, Germany  Dominik Stöckinger
<EMAIL_ADDRESS>Institut für Kern- und Teilchenphysik, TU
Dresden, Zellescher Weg 19, 01069 Dresden, Germany

###### Abstract

We study a simple dark matter model given by two interacting real singlet
scalars, with only one of them coupled to the Higgs. The model therefore
presents a minimal assisted-freeze-out framework: both scalars contribute to
the dark matter relic density, but only one of them takes part in the elastic
scattering with nuclei. This reduces the expected interaction rate in direct
detection experiments, in such a way that the model in some regions of the
parameter space may evades XENON1T constraints. We explore the model under
theoretical (perturbativity, stability potential and unitarity) and
experimental constraints (Higgs to invisible, relic density, direct
detection), with precise approximations near the Higgs resonance to calculate
the averaged annihilation cross section. We show that the model is viable in
the Higgs resonance region and for scalar singlet masses of hundreds of GeV.

###### Contents

  1. 1 Introduction
  2. 2 Model and Theoretical Constraints
  3. 3 Relic Abundance
    1. 3.1 Boltzmann equations
    2. 3.2 Assisted freeze-out
      1. 3.2.1 Higgs resonance
      2. 3.2.2 Outside the Higgs resonance ($m_{1}>m_{h}/2$)
  4. 4 Dark Matter Phenomenology
    1. 4.1 Experimental Constraints
    2. 4.2 Scan Results
  5. 5 Discussion and Summary
  6. 6 Acknowledgment
  7. A Average Annihilation Cross Sections

##  1 Introduction

While cosmological observations over the past years have made it more and more
apparent that Dark Matter (DM) is an important part of our universe [1],
direct detection experiments have put strong constraints on the parameter
space of many WIMP models, in particular one of the most popular DM models:
the real singlet scalar or also called Singlet Higgs Portal (SHP) [2, 3]. Even
when considering the singlet scalar as just a partial DM source together with
additional non-thermal contributions (e.g. axions [4]), most of the parameter
space around a few hundred GeV up to the TeV scale is ruled-out. Following
this minimal approach, models with two real singlet scalars have been studied
in the last years, reopening the possibility of having scalar DM near the weak
scale. Two classes of models have been explored: one scalar stabilized by a
$Z_{2}$ symmetry with the other scalar being unstable [5, 6, 7, 8, 9, 10, 11]
or both scalars stabilized by two symmetries $Z_{2}\times Z_{2}^{\prime}$ [10,
12, 13, 14, 15, 16] 111Fermions and vectors have also been considered in
similar frameworks, e.g. [17, 18, 19].. In the latter case, when the two
scalars are WIMPS with sizeable interaction between them [15], the model may
give the correct relic density on large portions of the parameter space, but
most of it is excluded by XENON1T direct detection results. Only the Higgs
resonance region, where a more careful treatment of the thermally averaged
annihilation cross section is necessary, seems to remain viable. In view of
this, it may be of interest to look for regions of parameter space in this
framework which evade the strong direct detection constraints, possibly
allowing DM with masses of several hundred GeV up to a few TeV. For instance,
in [10] such a two-component scalar model is explored under the assumption
that one of the scalars couples very weakly to the Higgs ($\lambda\sim
10^{-5}$), giving rise to a chronology in the decoupling process (exchange
driven freeze-out), with DM conversion being a pivotal process in the
determination of the relic abundance. They show that this mechanism
successfully evades XENON1T bounds, allowing DM with masses of hundreds of
GeV, provided certain quartic coupling take values above the unity.

In this work we study the two real scalar model, both stable gauge singlets,
with one of the Higgs portals set to zero. In this regime the model becomes a
minimal realization of an appealing idea, the so-called assisted freeze-out
mechanism [12]. The motivation of this approach is that only one of the
scalars interacts directly with baryonic matter through its Higgs portal,
whereas the other scalar may contribute most to the DM budget. In this way,
the Higgs coupled scalar may pass the direct detection constraints more easily
due to the rescaling of the expected DM direct detection rate by its partial
abundance: $\Omega_{coupled}/\Omega_{total}$. Qualitatively, this framework
can be viewed as an approximation of the exchange driven freeze-out [10].

We study the parameter space in detail to identify viable coupling and mass
ranges of the assisted freeze-out framework. For masses above the Higgs
resonance, we find that large couplings $\lambda\gtrsim 1$ are required in
order to pass both direct detection and relic density constraints, in
agreement with [10]. We then further impose theoretical constraints on the
resulting parameter space in order to see to what extent the resulting high
couplings are consistent with perturbative, scalar potential stability and
unitarity constraints.

We also specifically study the resonance effects for DM masses around 60 GeV,
half of the Higgs boson mass, where the real singlet scalar in SHP is not yet
ruled out. Here, the averaged annihilation cross section is strongly enhanced,
and we will show that a precise treatment involving the Narrow Width
Approximation (NWA) is necessary in order to get the correct relic abundance
of both scalars in this region, even when one of them is not directly coupled
to the Higgs.

In the following section, we present the model along with the theoretical
constraints. In section 3, we present the corresponding Boltzmann equation,
approximations to the averaged annihilation cross section and we study the
predicted relic abundance in detail, specifically in comparison with the SHP.
In section 4, we show some scan results inside and outside the resonance
region, taking into account the various constraints, and finally give a brief
summary and the conclusions in section 5.

##  2 Model and Theoretical Constraints

The model introduces two real gauge singlet scalars $S_{1}$ and $S_{2}$, each
one stabilized by $\mathbb{Z}_{2}\times\mathbb{Z^{\prime}}_{2}$ symmetry, with
the fields transforming as:

$\displaystyle\mathbb{Z}_{2}$ $\displaystyle\colon S_{1}\to-S_{1},\ S_{2}\to
S_{2},\ \text{SM}\to\text{SM},$ (1a) $\displaystyle\mathbb{Z}^{\prime}_{2}$
$\displaystyle\colon S_{1}\to S_{1},\ S_{2}\to-S_{2},\ \text{SM}\to\text{SM}.$
(1b)

Consequently, the scalar potential for this model is

$V(H,S_{1},S_{2})=\mu^{2}_{H}|H|^{2}+\lambda_{H}|H|^{4}+\sum_{i=1,2}\left(\frac{\mu^{2}_{i}}{2}S^{2}_{i}+\frac{\lambda_{i}}{4!}S^{4}_{i}+\frac{\lambda_{Hi}}{2}|H|^{2}S^{2}_{i}\right)+\frac{\lambda_{12}}{4}S^{2}_{1}S^{2}_{2}.$
(2)

The tree-level scalar masses after EWSB are

$\displaystyle m_{h}^{2}$ $\displaystyle=2v_{H}^{2}\lambda_{H},$ (3)
$\displaystyle m_{i}^{2}$
$\displaystyle=\mu_{i}^{2}+\dfrac{1}{2}v_{H}^{2}\lambda_{Hi},$ (4)

with $v_{H}$ the Higgs vacuum expectation value of 246 GeV. In order for the
singlet scalars to not acquire a vacuum expectation value, we set
$\mu_{i}^{2}>0$, then the physical scalar states are: $\\{h,S_{1},S_{2}\\}$.

According to the normalization of our scalar potential eq. (2), the stability
of the potential at tree-level requires [20, 21, 22]

$\displaystyle\lambda_{H}>0,\quad\lambda_{1}>0,\quad\lambda_{2}>0,\quad\lambda_{a}>0,\quad\lambda_{b}>0,\quad\lambda_{c}>0,$
(5)
$\displaystyle\sqrt{2\lambda_{H}\lambda_{1}\lambda_{2}/9}+\lambda_{H1}\sqrt{\lambda_{2}/3}+\lambda_{H2}\sqrt{\lambda_{1}/3}+\sqrt{2\lambda_{a}\lambda_{b}\lambda_{c}}>0,$
(6)

where $\lambda_{a}\equiv\lambda_{H1}+\sqrt{2\lambda_{H}\lambda_{1}/3}$,
$\lambda_{b}\equiv\lambda_{H2}+\sqrt{2\lambda_{H}\lambda_{2}/3}$ and
$\lambda_{c}\equiv\lambda_{12}+\sqrt{\lambda_{1}\lambda_{2}/9}$222From these
last constraints we obtain, for instance, that if $\lambda_{H1}$ is negative,
it must to satisfy $\lambda_{H1}^{2}<2\lambda_{H}\lambda_{2}/3$, and if it is
positive, then $\lambda_{H1}>0$.. Taking $m_{h}=$ 125 GeV in eq. (3) gives
$\lambda_{H}=0.13$.

On the other hand, tree-level perturbative unitarity puts constraints on the
zeroth partial wave of the system [23, 24], which for a $2\rightarrow 2$
elastic scattering is given by

$\displaystyle
a_{0}=\frac{1}{32\pi}\sqrt{\frac{4|\vec{p}_{f}||\vec{p}_{f}|}{s}}\int_{-1}^{+1}\mathcal{M}_{2\rightarrow
2}d\cos\theta,$ (7)

with $s$ the center-of-mass frame (CM) energy and $p_{i,f}^{\text{CM}}$ the
tree-momenta of initial and final particles in the CM. Unitarity then requires
$|\text{Re}(a_{0})|\leq\frac{1}{2}$. Considering the six-channel system

$\displaystyle\\{hh,S_{1}S_{1},S_{2}S_{2},S_{1}S_{2},hS_{1},hS_{2}\\},$ (8)

the zero-partial wave matrix at high energies $s\gg m_{h},m_{1},m_{2}$ is

$\displaystyle
a_{0}\rightarrow\frac{1}{16\pi}\begin{pmatrix}3\lambda_{H}&\lambda_{H1}&\lambda_{H2}&0&0&0\\\
\lambda_{H1}&\lambda_{1}&\lambda_{12}&0&0&0\\\
\lambda_{H2}&\lambda_{12}&\lambda_{2}&0&0&0\\\ 0&0&0&\lambda_{12}&0&0\\\
0&0&0&0&\lambda_{H1}&0\\\ 0&0&0&0&0&\lambda_{H2}\\\ \end{pmatrix}.$ (9)

Perturbative unitarity constraints on the matrix $a_{0}$ result in bounds on
its eigenvalues $e_{i}$. In this work we concentrate on the case in which
$\lambda_{H2}=0$, which then results in the following constraint on the
eigenvalues of the $3\times 3$ upper block

$\displaystyle|e_{i}(\lambda_{1},\lambda_{2},\lambda_{H1},\lambda_{12})|\leq
8\pi,\quad i=1,2,3.$ (10)

For simplicity we omit the analytical expressions for each eigenvalue.

Another theoretical constraint is given by perturbativity, which implies
approximate upper limits for all coupling parameters in order for the
perturbative expansion to be reliable. At this point we do not choose any
specific value for these upper limit. In the following discussions we take
into account all three theoretical constraints (stability, unitarity and
perturbativity) and analyse their impact on the viable parameter space.

Figure 1: Unitarity constraints on the pair $(\lambda_{H1},\lambda_{12})$ for
$\lambda_{1,2}=1,4\pi,7\pi$, respectively, with the blue region excluded by
this constraint. The region inside of each red square may be interpreted as
the allowed parameter space by perturbativity when the perturbative limit on
the couplings is $1,4\pi,7\pi$, respectively.

Considering only positive values for $\lambda_{H1}$ and $\lambda_{12}$
(negative values do not change the analysis of our work) then the stability
potential constraints eq. 5 do not give any new information. In Fig. 1 we show
the unitarity bounds on $\lambda_{H1}$ and $\lambda_{12}$ for the choices
$\lambda_{1,2}=1,4\pi$ and $7\pi$, respectively, with the blue region excluded
by this constraint. The bigger the values of the quartic couplings
$\lambda_{1,2}$ the stronger are the constraints on the remaining pair
$\lambda_{H1},\lambda_{12}$. In particular, as the middle panel shows, it is
not possible to have the quartic couplings $\lambda_{1,2}$ and
$\lambda_{H1},\lambda_{12}$ all simultaneously with values above $4\pi$.
Additionally, the region inside the red squares indicates the regions allowed
by different choices of the perturbativity limit. The limits are set to $1$,
$4\pi$, and $7\pi$, respectively, i.e. to the same values as the quartic
couplings in each panel. Clearly, perturbativity places limits on each
individual coupling, while unitarity places limits on combinations of
couplings. If very strong perturbativity limits are imposed as in the left
panel, unitarity does not imply additional constraints, while very loose
perturbativity limits, as in the right panel, are much less stringent than
unitarity limits.

Therefore, the analysis shows that the theoretical constraints depend on both
the assumed values for the quartic couplings $\lambda_{1,2}$ as well as the
choice of the perturbativity limit. In the last part of this work, when we
contrast theoretical and experimental constraints on the model, we discuss the
viability of parameter space depending on the chosen criteria.

##  3 Relic Abundance

In this section we discuss the main formulas describing the evolutions of the
DM number density throughout the thermal history. We introduce precise
approximations for the averaged annihilation cross section and study the
behaviour of the DM relic density both in and outside the Higgs resonance
region. We also point out that a careful treatment of the averaged
annihilation cross section near the Higgs resonance is necessary in order to
obtain the correct results.

###  3.1 Boltzmann equations

The coupled Boltzmann equations for the DM number densities of $S_{1}$ and
$S_{2}$ are given by

$\displaystyle\begin{split}\dfrac{dn_{1}}{dt}+3Hn_{1}=&-\braket{\sigma_{11XX}v}\Bigl{[}n_{1}^{2}-n_{1e}^{2}\Bigr{]}-\braket{\sigma_{1122}v}\Bigl{[}n_{1}^{2}-\dfrac{n_{1e}^{2}}{n_{2e}^{2}}n_{2}^{2}\Bigr{]}\Theta(m_{1}-m_{2})\\\
&+\braket{\sigma_{2211}v}\Bigl{[}n_{2}^{2}-\dfrac{n_{2e}^{2}}{n_{1e}^{2}}n_{1}^{2}\Bigr{]}\Theta(m_{2}-m_{1}),\\\
\dfrac{dn_{2}}{dt}+3Hn_{2}=&-\braket{\sigma_{22XX}v}\Bigl{[}n_{2}^{2}-n_{2e}^{2}\Bigr{]}-\braket{\sigma_{2211}v}\Bigl{[}n_{2}^{2}-\dfrac{n_{2e}^{2}}{n_{1e}^{2}}n_{1}^{2}\Bigr{]}\Theta(m_{2}-m_{1})\\\
&+\braket{\sigma_{1122}v}\Bigl{[}n_{1}^{2}-\dfrac{n_{1e}^{2}}{n_{2e}^{2}}n_{2}^{2}\Bigr{]}\Theta(m_{1}-m_{2}),\end{split}$
(11)

with $X$ being some SM particle, $\braket{\sigma v}$ the thermal average of
the cross section times the relative velocity of a given process
$ii\rightarrow jj$, with $i,j=1,2,X$, and, assuming a Maxwell-Boltzmann
thermal distribution, the equilibrium number densities $n_{ie}$ as a function
of the thermal bath temperature $T$ are given by

$\displaystyle
n_{ie}(T)=g_{i}\dfrac{m_{i}^{2}}{2\pi^{2}}TK_{2}(\tfrac{m_{i}}{T}),$ (12)

with the internal degrees of freedom $g_{i}$, and the modified Bessel function
of the second kind $K_{2}$. Furthermore, the Hubble parameter $H$ and entropy
density $s$ in the early universe can be written as

$\displaystyle H(T)=\sqrt{\dfrac{4\pi^{3}G}{45}g_{*}(T)}\cdot T^{2},\quad
s(T)=\dfrac{2\pi^{2}}{45}g_{*s}(T)\cdot T^{3},$ (13)

with $G$ is the Newton gravitational constant, and $g_{*}(T)$ and $g_{*s}(T)$
are the effective degrees of freedom contributing respectively to the energy
and the entropy density at temperature $T$. For our numerical evaluations we
interpolate $g_{*}$ and $g_{*s}$ from the values given in [25]. We have
introduced Heaviside functions in eq. (11) in order to select the formulation
of the interaction term in which $n_{ie}^{2}/n_{je}^{2}\rightarrow 0$ at small
$T$, making the numerical evaluation more convenient.

Lastly, introducing $x:=\mu/T$, where $\mu=m_{1}m_{2}/(m_{1}+m_{2})$ and
$Y_{i}:=n_{i}/s$, we rewrite eq. (11) in the usual dimensionless form

$\displaystyle\begin{split}\dfrac{dY_{1}}{dx}=&-\lambda_{11XX}\cdot[Y_{1}^{2}-Y_{1e}^{2}]-\lambda_{1122}\cdot[Y_{1}^{2}-r_{12}^{2}Y_{2}^{2}]\cdot\Theta(m_{1}-m_{2})\\\
&+\lambda_{2211}\cdot[Y_{2}^{2}-r_{21}^{2}Y_{1}^{2}]\cdot\Theta(m_{2}-m_{1}),\\\
\dfrac{dY_{2}}{dx}=&-\lambda_{22XX}\cdot[Y_{2}^{2}-Y_{2e}^{2}]-\lambda_{2211}\cdot[Y_{2}^{2}-r_{21}^{2}Y_{1}^{2}]\cdot\Theta(m_{2}-m_{1})\\\
&+\lambda_{1122}\cdot[Y_{1}^{2}-r_{12}^{2}Y_{2}^{2}]\cdot\Theta(m_{1}-m_{2}),\end{split}$
(14)

with

$\displaystyle
r_{ij}(x):=\dfrac{Y_{ie}(x)}{Y_{je}(x)},\quad\lambda_{iijj}(x):=\dfrac{\langle\sigma_{iijj}v\rangle(x)\cdot
s(T)}{x\cdot H(T)}\qquad\text{for}\leavevmode\nobreak\ i,j=1,2,X.$ (15)

The thermal average annihilation cross section is given by

$\displaystyle\langle\sigma_{iijj}v\rangle(x)=\int_{4m_{i}^{2}}^{\infty}\dfrac{s\sqrt{s-4m_{i}^{2}}K_{1}(\sqrt{s}\tfrac{x}{\mu})(\sigma
v)_{iijj}}{16\frac{\mu}{x}m_{i}^{4}K_{2}^{2}(m_{i}\tfrac{x}{\mu})}ds,$ (16)

where $K_{1}$ and $K_{2}$ are modified Bessel functions of the second kind.
Properly evaluating this expression plays a significant role in determining
the DM relic density, especially near resonances and thresholds. However, the
numerical evaluation is computationally quite expensive and does not give
immediate insight into the analytic behavior. We will therefore briefly
discuss some approximations which speed up the calculations significantly and
provide a nice basis for analytical understanding. In particular, we consider
the following two scenarios: (i) the cross section has a resonance at some
$s\gtrsim 4m_{i}^{2}$ or (ii) the cross section is sufficiently smooth near
$s\gtrsim 4m_{i}^{2}$.

  1. (i)

The relevant resonance is the Higgs resonance, corresponding to the decay
$h\to S_{i}S_{i}$. A well known approximation for a cross section near a
resonance is the replacement of the propagator by a Dirac delta function
(Narrow Width Approximation):

$\displaystyle\dfrac{1}{(s-m_{h}^{2})^{2}+m_{h}^{2}\Gamma^{2}}\approx\dfrac{\pi}{m_{h}\Gamma}\delta(s-m_{h}^{2}),$
(17)

where $\Gamma$ is the total Higgs decay width, including the invisible decays
into DM particles. Using this approximation, the integral in eq. (16) becomes
an algebraic expression.

  2. (ii)

Since the Bessel functions under the integral in eq. (16) suppresses all
contributions of the cross section for $s\gg 4m_{i}^{2}$, as long as the cross
section is sufficiently smooth near $s\gtrsim 4m_{i}^{2}$, one can write to
good approximation:

$\displaystyle\sigma v(s)\approx\sigma v(s_{0})+(\sigma
v)^{\prime}(s_{0})\cdot(s-s_{0}).$ (18)

Choosing $s_{0}=\langle s\rangle$ we find $\langle\sigma v\rangle\approx\sigma
v(\langle s\rangle)$ which we will call the naive approximation333 $\langle
s\rangle=4m_{i}^{2}+12m_{i}T\cdot
K_{3}(\tfrac{m_{i}}{T})/K_{2}(\tfrac{m_{i}}{T})\approx
4m_{i}^{2}(1+\tfrac{3\mu}{m_{i}x})$. Often times the small $x$-dependent
correction term in $\langle s\rangle$ is ignored in which case we refer to it
as the constant approximation $\langle\sigma v\rangle\approx\sigma v(4m^{2})$
.

Figure 2: Comparison of $\braket{\sigma_{11XX}v}$444The annihilation cross
section for the Higgs portal can be found e.g. in the appendix of [2]. The
annihilation cross section for the scalar interaction is given by
$\displaystyle\sigma_{iijj}v(s)=\frac{\sqrt{s-4m_{j}^{2}}}{8\pi
s\sqrt{s}}\left[\frac{\lambda_{H1}^{2}\lambda_{H2}^{2}v_{H}^{4}}{(s-m_{h}^{2})^{2}+m_{h}^{2}\Gamma^{2}}+\frac{2(s-m_{h}^{2})\lambda_{H1}\lambda_{H2}\lambda_{12}v_{H}^{2}}{(s-m_{h}^{2})^{2}+m_{h}^{2}\Gamma^{2}}+\lambda_{12}^{2}\right].$
(19) as a function of $m_{1}$, calculated with NWA (red), the naive
approximation (green solid) and the constant approximation (green dashed). In
the left plot was set $x=5$, whereas in the right one $x=20$. In both plots we
consider $\lambda_{H1}=0.1$ and a solid blue line showing the micrOMEGAs
prediction. The meaning of (i) and (ii) are explained in the text.

As was previously highlighted, for scalar masses $m_{i}$ close to the Higgs
resonance (i.e. $m_{i}\lesssim 60$ GeV) the peak in the annihilation cross
section contributes significantly to the integrand in eq. (16) and, especially
at high temperatures where the Bessel suppression is significant only at very
large $s$, the thermal average is therefore up to several orders of magnitude
larger than the cross section itself.

In order to quantify the effects of this enhancement, in Fig. 4 we show a
comparison of the three mentioned approximations at different temperatures (or
$x$), considering also the result obtained with the public tool
micrOMEGAs[26], which solves 16 with a precision of about 1%. In the left
plot, setting $x=5$, the NWA (solid red curve) provides an accurate
approximation that matches very well the micrOMEGAs result (solid blue) and
shows the enhancement of $\braket{\sigma_{11XX}v}$ over the approximations
based on eq. (18) (green lines). In this way, the naive and constant
approximations are unreliable and underestimate the true result in this
region. At higher $x$, as it is shown in the right plot for $x=20$, the
resonance affects the integrand less due to the Bessel suppression, resulting
in the NWA approaching the constant approximation. Away from resonance, the
continuous green line reproduces the micrOMEGAs result very well even for
small $x$, while the dashed green line shows significant deviations.

Therefore, the value of $\braket{\sigma_{11XX}v}$ in region (i) (i.e. for
$m_{i}<m_{H}/2$) can be approximated by the maximum value between the NWA and
constant approximation, and in region (ii) by the naive approximation. This
allows for efficient computation of large scans of the parameter space due to
the fact that we do not have to deal with the integral eq. (16), but only with
simple algebraic expressions for $\braket{\sigma_{11XX}v}$. As an outlook, we
remark that improvements on the annihilation cross section considering higher
order corrections may be done for the threshold effects [2, 27]. Here, we do
not develop a treatment of this, but the corrections coming from such effects
are considered in our full calculation given in section 4. We give a more
detailed discussion on the DM relic density in the assisted freeze-out and its
dependence on the resonance effects as well as the free parameters in the
following section.

Finally, after solving the Boltzmann equations the DM density parameter is
given by

$\displaystyle\Omega
h^{2}=\Omega_{1}h^{2}+\Omega_{2}h^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \qquad\Omega_{i}h^{2}=\dfrac{2.9713\cdot
10^{9}}{10.5115\leavevmode\nobreak\ \text{GeV}}\cdot Y_{i,0}\cdot m_{i},$ (20)

where $Y_{i,0}$ is the yield of $S_{i}$ today, i.e. after freeze-out.

###  3.2 Assisted freeze-out

We now focus on the assisted freeze-out scenario, i.e. we set one Higgs portal
coupling to zero, $\lambda_{H2}=0$. In this way, from now on only $S_{1}$ will
be coupled to the SM, whereas $S_{2}$ becomes decoupled. Consequently,
$\lambda_{22XX}=0$ in eq. (14) and the relevant parameters of this scenario
are $m_{1},m_{2},\lambda_{H1}$ and $\lambda_{12}$. We allow sizable couplings
$\lambda_{H1}$ and $\lambda_{12}$, such that the coupled scalar starts in
thermal equilibrium with the SM, and at the same time brings the decoupled one
into the equilibrium too555The case in which $\lambda_{H2}$ takes
$\mathcal{O}(1)$ values was studied in major detail in [15]. If $\lambda_{H2}$
takes much lower values, $\lesssim 10^{-5}$ [10], the interactions of $S_{2}$
with the SM become negligible..

In order to study the evolution of the DM yield $Y_{i}$, we separate the
analysis in two regions: $m_{1}\approx m_{h}/2$ (coupled scalar is at Higgs
resonance), where resonance effects have to be taken into account, and
$m_{1}>m_{h}/2$. Due to the obvious correspondence, in some parts of our
analysis we compare the coupled scalar $S_{1}$ to the single real singlet
scalar $S$ of the Singlet Higgs Portal (SHP) [2]. In the following comparisons
we will therefore always assume the same mass and Higgs portal coupling for
$S_{1}$ and $S$.

####  3.2.1 Higgs resonance

Figure 3: Comparison of relic abundance calculations for evolution of $Y_{1}$
(blue curves) and $Y_{2}$ (red curves) and for the real singlet scalar $S$
(grey). The parameters are $(m_{1},\lambda_{H1},\lambda_{12})=(60,0.01,0.1)$
and $m_{2}=50$ GeV (left) and $m_{2}=70$ GeV (right). The continuous lines
correspond to the predictions using the NWA, the dashed lines using the
constant approximation, and the dotted points are the respective equilibrium
densities. In the legend we also give the corresponding relic abundance
$\Omega_{i}$ for each plot.

As was shown in last section, the averaged annihilation cross section is
drastically enhanced near the Higgs resonance. The subsequently large
interaction rates result in the DM particles dropping out of equilibrium at
later times, leading to a lower relic density.

In Fig. 3, we exemplify this fact taking $m_{1}=60$ GeV, a value well close to
$m_{h}/2$. The left plot shows a typical evolution of the DM densities in the
inverted mass hierarchy ($m_{1}>m_{2}$, i.e. heavy coupled scalar), while the
right plot shows the evolution in the normal mass hierarchy ($m_{1}<m_{2}$,
i.e. light coupled scalar). The relic density is significantly smaller when
$\braket{\sigma v}$ is calculated using the NWA instead of using the constant
approximation. As shown in Fig. 3, this affects both scalars in the inverted
hierarchy (left plot), whereas in the normal hierarchy (right plot) the
decoupled scalar $S_{2}$ is mostly unaffected. This is explained by the fact
that in the inverted hierarchy ($m_{2}<m_{1}$) the evolution of $Y_{2}$
depends on $Y_{1}$ ($dY_{2}/dx\approx\lambda_{1122}Y_{1}^{2}$) much more
strongly than in the normal hierarchy
($dY_{2}/dx\approx-\lambda_{2211}Y_{2}^{2}$). Since the NWA gives much larger
values for $\braket{\sigma_{11SM}v}$ than the constant approximation, we
expect $S_{1}$ to freeze out later in both hierarchies, which means that
$S_{2}$ also stays in equilibrium for longer. Hence we note, that a precise
treatment of the averaged annihilation cross section near the Higgs resonance
leads to results significantly differing from those using just the constant
approximation and is therefore necessary for obtaining accurate predictions.

Lastly, in comparison with the singlet scalar $S$, we notice an increase in
the relic density of the coupled scalar $S_{1}$ in both hierarchies, making
evident the effect of the scalar interaction on the thermal evolution of
$Y_{1}$. This suggests, that in general for a given mass $m_{1}$ and coupling
$\lambda_{H1}$ the assisted freeze-out produces a much larger DM relic than
the SHP. The remaining parameters $\lambda_{12}$ and $m_{2}$ must therefore be
chosen carefully in order to obtain the observed relic density. In particular
we expect $\lambda_{12}$ to take relatively large values in order to keep the
particles in equilibrium for longer. In the following we discuss the
dependence of the DM relic density on these free parameters in more detail.

Figure 4: DM density parameter $\Omega h^{2}$ (green), $\Omega_{1}h^{2}$
(blue) and $\Omega_{2}h^{2}$ (red) as a function of $m_{2}$ (left) and
$\lambda_{12}$ (right), both plots with $m_{1}$ = 60 GeV and $\lambda_{H1}$ =
0.01. In the left plot we have set $\lambda_{12}=1$, and in the right plot
$m_{2}=50$ GeV (dashed) and $m_{2}=70$ GeV (solid).

As was already pointed out, the dynamics of the abundance evolution is highly
dependent on the mass hierarchy of the two singlet scalars. The behavior of
the resulting relic densities can be seen in more detail in Fig. 4 (left),
which fixes $m_{1}=60$ GeV and varies $m_{2}$. For $m_{2}<m_{1}$ (inverted
hierarchy) the overwhelmingly dominant contribution to the total relic comes
from $\Omega_{2}$, while in the normal mass hierarchy ($m_{2}>m_{1}$)
$\Omega\approx\Omega_{2}$ only at very high mass differences with a much
smaller relative difference between $\Omega_{1}$ and $\Omega_{2}$. Again, this
behavior can be understood by looking at the corresponding Boltzmann Equation.
In the normal hierarchy we already mentioned that if $m_{2}\gg m_{1}$, we
expect $Y_{2}$ to evolve independently of $Y_{1}$. In fact, since
$\langle\sigma_{2211}v\rangle\sim m_{2}^{-2}$, eq. (14) effectively decouples,
resulting in $Y_{1}\sim const.$ (i.e. $\Omega_{1}\sim const.$) and $Y_{2}\sim
m_{2}^{2}$ ($\Omega_{2}\sim m_{2}^{3}$), as it is shown in Fig. 4 (left) with
the blue and red lines, respectively. Although in the normal hierarchy at high
mass differences we have $\Omega_{2}\gg\Omega_{1}$, we typically obtain a too
small overall DM relic abundance. This suggests that we would have to go to
smaller couplings $\lambda_{H1}$, which in turn will produce larger relative
contribution of $\Omega_{1}$ to $\Omega$. In the inverted hierarchy,
constraints on the DM relic abundance also put very strong constraints on the
possible masses $m_{2}$ for given $m_{1},\lambda_{H1}$ and $\lambda_{12}$.
This will be an important point to keep in mind when choosing the mass
resolution in parameter space scans.

Lastly, considering the $\lambda_{12}$ dependence of the DM relic abundance,
we would expect a larger coupling to lead to a later freeze out and therefore
lower relic densities. For small couplings $\lambda_{12}\sim\lambda_{H1}$ this
behavior can indeed be seen in Fig. 4 (right), however, at large couplings
$\lambda_{12}\gg\lambda_{H1}$, $\Omega$ becomes asymptotically constant. In
the normal mass hierarchy this is again explained by the fact that
$dY_{2}/dx\sim-\lambda_{2211}Y_{2}^{2}$, and therefore $Y_{2}\sim
1/\lambda_{12}^{2}$. Since the interaction term is subsequently roughly
proportional to $\lambda_{12}^{2}Y_{2}^{2}\sim 1/\lambda_{12}^{2}$, we find
$Y_{1}\sim const.$ at large couplings. Interestingly, in the inverted mass
hierarchy the total relic abundance behaves very much in the same way,
although the independence of $\Omega$ on $\lambda_{12}$ at large couplings now
comes from the fact that $Y_{2}\sim const$. The explanation is very similar to
that in the normal hierarchy: we have $dY_{1}/dx\sim-\lambda_{1122}Y_{1}^{2}$,
i.e. $Y_{1}\sim 1/\lambda_{12}^{2}$ and therefore $Y_{2}\sim const$. Overall,
the inverted mass hierarchy with large couplings $\lambda_{12}$ therefore seem
to be the most promising region, giving both a large enough overall abundance
$\Omega$ as well as $\Omega_{1}/\Omega\ll 1$.

####  3.2.2 Outside the Higgs resonance ($m_{1}>m_{h}/2$)

Figure 5: The same than the previous figure, but this time with $m_{1}$ = 100
GeV and $\lambda_{H1}$ = 0.1. In the left plot we have set
$\lambda_{12}\leavevmode\nobreak\ =\leavevmode\nobreak\ 1$, and in the right
plot $m_{2}=90$ GeV (dashed) and $m_{2}=110$ GeV (solid).

Outside of the Higgs resonance, the DM into SM annihilation rate is much
smaller, which leads to an earlier freeze-out and a much larger relic
abundance (at equal couplings $\lambda_{H1}$ and $\lambda_{12}$). The
dependence of $\Omega_{i}$ on the coupling $\lambda_{12}$ and on the mass
$m_{2}$ still qualitatively follows our discussion in the last section, except
that the effect of the mass hierarchy is strongly enhanced, as it is shown in
Fig. 5 (left). This behavior is expected, since as we discussed in the last
section, the mass dependence is mainly dictated by the scalar interaction
term, which now makes up a larger relative contribution to the eq. (14), due
to the fact that $\langle\sigma_{11XX}v\rangle$ is much smaller outside the
Higgs resonance region. This fact can also be seen in Fig. 5 (right) where,
despite larger coupling $\lambda_{H1}$, $\Omega$ begins to flatten at smaller
$\lambda_{12}\sim 0.1$ than previously, again is explained by the scalar
interaction term having a much stronger relative effect on the evolution. We
have now developed a sufficient understanding of the relic abundance behavior
and can move on to some actual predictions of the model, taking into account
the various experimental constraints on its parameter space.

##  4 Dark Matter Phenomenology

In this section, we start reviewing relevant experimental constraints, such as
Higgs decay to invisible at the LHC, relic density and direct detection, to
subsequently applied them onto the promising parameter space of the model.
Similarly to the last section, we focus on the Higgs resonance region
separately and contrast the results with the SHP predictions. Finally, we
present some scans on the parameter space outside the Higgs resonance,
contrasting those points that fulfill the experimental constraints with
theoretical bounds.

###  4.1 Experimental Constraints

The strongest constraint on the parameter space is given by the observed DM
relic abundance. We use the value given by the most recent Planck data [1]:

$\displaystyle\Omega_{c}h^{2}=0.11933\pm 0.00091,\qquad\text{at }2\sigma.$
(21)

We will denote this measured value of the DM relic by $\Omega_{c}$ and the DM
relic predicted by the model by $\Omega$. This gives a parameter space
constraint $\Omega=\Omega_{c}$, which we applied with a tolerance of $\pm 5\%$
accounting for uncertainties in our computation.

The second constraint we need to pass is given by direct detection experiments
such as XENON1T. These give an upper bound on the spin-independent (SI) DM
nucleon scattering cross section $\sigma_{SI}$,

$\displaystyle\sigma_{SI,i}=\dfrac{\lambda_{Hi}^{2}f_{N}^{2}}{4\pi}\dfrac{\mu_{i}^{2}m_{n}^{2}}{m_{h}^{4}m_{i}^{2}},\quad
i=1,2,$ (22)

where $f_{N}$ is the effective Higgs-nucleon coupling, $m_{n}$ is the nucleon
mass and $\mu_{i}=m_{n}m_{i}/(m_{n}+m_{i})$ is the DM-nucleon reduced mass [2,
3]. Since in our case only one of the scalars interacts with the SM, i.e.
$S_{1}$, the differential detection rate is
$dR/dE\propto\sigma_{SI}n_{1\odot}=\sigma_{SI}(\Omega_{1}/\Omega)n_{\odot}$,
where $n_{\odot}$ denotes the local DM density. This effectively rescales the
direct detection constraint on our model to
$(\Omega_{1}/\Omega)\sigma_{SI,1}\leq\sigma_{Xe}$, where $\sigma_{Xe}$ is the
upper bound given by XENON1T [28]. In the SHP, it has been shown that indirect
bounds set much weaker constraints on the parameter space than the direct
detection limits [2]. In our scenario, this still holds, since the only
parameter space passing the stringent XENON1T constraints is in the inverted
mass hierarchy where $\Omega_{1}/\Omega\ll 1$, therefore primary fluxes
produced from the annihilation of a pair of $S_{1}$ are highly suppressed.

Lastly, for $m_{1}\leq m_{h}/2$, the Higgs boson can decay into two $S_{1}$,
with an invisible decay width given by

$\displaystyle\Gamma_{inv}(h\to
S_{1}S_{1})=\frac{\lambda^{2}_{H1}v_{H}^{2}}{32\pi
m_{h}}\sqrt{1-\frac{4m^{2}_{1}}{m^{2}_{h}}}.$ (23)

This contributes to the invisible branching rate
$\text{Br}(h\to\text{inv})=\Gamma_{\text{inv}}/(\Gamma_{\text{SM}}+\Gamma_{\text{inv}})$,
with $\Gamma_{\text{SM}}=4.07$ MeV [29]. Experimental searches put strong
constraints on this quantity, with the most stringent value being given by
$\text{Br}(h\to\text{inv})<0.19$ at $95\%$ C.L. [29].

###  4.2 Scan Results

Figure 6: Contours of fixed relic abundance in the assisted freeze-out model
(solid) and in the SHP (dashed) in the parameter space $(m_{1},\lambda_{H1})$
around the Higgs resonance. All contours of the assisted freeze-out model
fulfil $\Omega=\Omega_{c}$ with $\lambda_{12}=1$, for different mass
splittings $m_{2}-m_{1}$. The light red region is ruled out due to the Higgs
to invisible constraint, and the region above the gold contour is ruled out in
the SHP by the XENON1T data. The red solid lines denote the regions ruled out
by the direct detection constraints on the assisted freeze out model.

In Fig. 6 we show some numerical values in the resonance region to exemplify
the general behaviour here. In this plot, we have considered contours with
fixed mass differences $m_{2}=m_{1}+1$ GeV (blue), $m_{1}-1$ GeV (green) and
$m_{1}-5$ GeV (grey), which give the observed DM relic abundance at
$\lambda_{12}=1$. As expected, in the normal mass hierarchy $S_{1}$ behaves
very similar to $S$, however in the inverted hierarchy, $\Omega_{1}$ drops
enough such that direct detection constraint becomes effectively negligible,
only ruling out parameter space outside the Higgs resonance region. For
comparison, we have also plotted contours where the SHP could account for the
100%, 10% and 1% of the correct relic abundance with a black, blue and red
dashed lines, respectively. As expected, due to the additional degrees of
freedom coming from $m_{2}$ and $\lambda_{12}$, the two scalar model can
explain the entire DM relic abundance in that region but also even at slightly
smaller masses $m_{1}$. Thus, for masses near the Higgs resonance the SHP
parameter space opens up in the two scalar model and it remains to be seen,
whether WIMPs in the region can be fully ruled out as DM candidates.

Figure 7: Random scan on the inverted mass hierarchy parameter space for fixed
$m_{1}=100$ GeV and $\lambda_{H2}=0$. The red points fulfill the direct
detection constraints, while the blue-green points fulfill
$\Omega\leq\Omega_{c}$ (with $\Omega$ indicated by the color bar). The left
plot shows points projected on the 3d parameter space
$(m_{1},\lambda_{H1},\lambda_{12})$, the middle plot shows the projection on
the plane $m_{1}-\lambda_{H1}$ plane, and the right plot shows just the
overlap of points fulfilling the direct detection constraints and the relaxed
abundance constraint $\Omega\leq\Omega_{c}$.

Outside the Higgs resonance region, the most promising parameter space is in
the inverted mass hierarchy, as by our analysis from section 3.2, both the
abundance and direct detection constraints are most easily fulfilled there. We
thus ran larger scans of points $(m_{2},\lambda_{H1},\lambda_{12})$ for fixed
$m_{1}=100$ GeV (Fig. 7) and $m_{2}<m_{1}$. In this region, points fulfilling
the direct detection constraints (red) are either at very small couplings
$\lambda_{H1}$ and arbitrary $\lambda_{12}$, or, for larger couplings
$\lambda_{H1}$ and $\lambda_{12}$ at masses $m_{2}$ mostly below that
fulfilling the abundance constraint for the given couplings. Since in the
inverted mass hierarchy a lower relic abundance at equal couplings can only be
achieved by a smaller mass difference, points fulfilling $\Omega<\Omega_{c}$
lie opposite from the points fulfilling the direct detection constraints, with
the points fulfilling $\Omega=\Omega_{c}$ forming a boundary. Hence the
overlap of these two regions does not increase significantly when relaxing the
abundance constraint. This fact is shown in Fig. 7 (middle) where the small
overlap of points fulfilling the direct detection and relaxed abundance
constraint is shown. In particular from Fig. 7 (right) we see no overlap for
$\lambda_{12}\lesssim 6$, since smaller $\lambda_{12}$ also requires small
$\lambda_{H1}$ to pass the direct detection constraint which in turn produces
a DM over-abundance.

Figure 8: Allowed points by experimental constraints considering $m_{1}=100$
GeV (left) and 500 GeV (right) projected in the plane
$(\lambda_{H1},\lambda_{12})$. The color bar indicates the mass of $S_{2}$.
The dashed lines correspond to the maximum value allowed for $\lambda_{12}$
considering that the theoretical constraints allowed quartic couplings up to
the corresponding value indicated on the legend.

To get a more precise picture of the interplay between the parameter space
allowed by the experimental and the theoretical constraints, in Fig. 8 we show
a plot of the allowed points with $m_{1}=100$ and 500 GeV projected in the
plane ($\lambda_{H1},\lambda_{12}$), with a color bar for $m_{2}$. The left
plot shows the same parameter space as Fig. 7, and confirms that sizable
values for $\lambda_{12}$ are required in order to fulfill the experimental DM
constraints. Further, as shown in the right plot, for a higher $m_{1}$, the
smallest possible values for $\lambda_{12}$ increase as well. Interestingly,
the smallest achievable $\lambda_{12}$ couplings in each case are around 6 and
8 (left/right plots), and they are achieved for very particular mass
splittings, 1 GeV in the left case and around 5 GeV in the right plot.

As was argued in Section 2, it is necessary that the coupling values be in
agreement with perturbativity and unitarity, which depend on the values of the
quartic couplings. In the Fig. 8 we show the contours of the unitarity
constraints for different $\lambda_{1,2}$. The points below each dashed line
are allowed by unitarity for a certain range of values of $\lambda_{1,2}$. For
instance, all the points below the yellow dashed line are allowed for
$\lambda_{1,2}$ up to $5\pi$. If the values of $\lambda_{1,2}$ are below
$4\pi$, then the corresponding unitarity dashed line rises above the red
dashed line, and the allowed parameter space will depend only on the
perturbative limit assumed on $\lambda_{H1},\lambda_{12}$.

##  5 Discussion and Summary

In this paper we investigated to what extent the parameter space of the
Singlet Higgs Portal can be reopened by introducing a second DM candidate
decoupled from the SM: Assisted freeze-out framework. The outcome is positive.
In studying the two real stable scalar model, we have implemented our own
numerical solution of the Boltzmann equation using efficient but precise
approximations for the averaged annihilation cross section, which allowed us
to compute large scans of the parameter space of our model.

We conclude that, in order to evade direct detection and relic abundance
constraints outside the Higgs resonance region it is required that one of the
key new parameters, $\lambda_{12}$, take high values, which under some
theoretical conditions based on stability potential, perturbativity and
unitarity, the parameter space reopen around scales of hundreds of GeV.
Additionally, our results are in concordance with what was found [10], where
they studied the model in an equivalent framework with results exemplified for
a 1% of split difference between the scalar masses. In the present work, we
explored in major detail the parameter space, concluding that outside the
resonance region, it is allowed to have a split mass difference up to 10% in
the most optimistic cases, i.e. if the coupling $\lambda_{12}$ is very large,
above 10. This is allowed by tree-level unitarity provided other Higgs
potential parameters are not above $4\pi$.

Similarly to the SHP, the two real scalar model in the assisted freeze-out
regime shows that parts of the Higgs resonance region remains not (yet) ruled
out by experimental data. First, even when one of the scalars is not coupled
directly to the Higgs, we have shown that the relic abundances of both scalars
near the Higgs resonance region may depend strongly on resonance effects,
changing their values in more than one order of magnitude in some cases using
the NWA with respect to other approaches for the average annihilation cross
section (e.g. constant approximation). Furthermore, due to the additional
degrees of freedom $m_{2}$ and $\lambda_{12}$, the two scalar model in the
assisted freeze-out framework could completely explain DM for a broader range
of values $m_{1}$ and $\lambda_{H1}$, however, the explanation would still
require $m_{1}$ in the region of $m_{h}/2$ and simultaneously a very small
mass splitting between $m_{1}$ and $m_{2}$.

##  6 Acknowledgment

B.D.S. would like to thank Roger Galindo, Felipe Rojas and Sebastian Norero
for useful discussions, and to CONICYT Grant No. 74200120.

##  Appendix A Average Annihilation Cross Sections

For a scattering process $ij\rightarrow ab$, the average cross section times
velocity following a Maxwell-Boltzmann thermal distribution is given by

$\displaystyle\langle\sigma_{ijab}v\rangle=\dfrac{1}{n_{ie}n_{ie}}\int\dfrac{d^{3}p_{i}}{(2\pi)^{3}}\int\dfrac{d^{3}p_{j}}{(2\pi)^{3}}(\sigma_{ijab}v)\cdot
e^{-(E_{i}+E_{j})/T},$ (24)

where $v$ denotes the Møller velocity, and $n_{ie}$ the equilibrium densities
given in 12. The cross sections for the processes we are going to consider can
be written as functions of
$s:=(E_{i}+E_{j})^{2}-(\vec{p_{i}}+\vec{p_{j}})^{2}$. The integrals in Eq. 24
then break down to just a single integral over s

$\displaystyle\langle\sigma
v\rangle=\int_{(m_{i}+m_{j})^{2}}^{\infty}ds\dfrac{(\sigma v)^{(cm)}\cdot
s\cdot
K_{1}(\tfrac{\sqrt{s}}{T})}{16Tm_{i}^{2}m_{j}^{2}K_{2}(\tfrac{m_{i}}{T})K_{2}(\tfrac{m_{j}}{T})}\sqrt{\dfrac{(m_{i}^{2}-m_{j}^{2})^{2}}{s}+s-2(m_{i}^{2}+m_{b}^{2})},$
(25)

where $(\sigma v)^{(cm)}$ denotes the calculation in the center of mass frame
and $K_{\alpha}$ is the modified Bessel function of the second kind.

Similarly, $\langle s\rangle$ can now be calculated in reverse. Since
$s=4m_{i}^{2}+4p_{(cm)}^{2}$, then

$\displaystyle\langle
s\rangle=\dfrac{1}{n_{ie}}\int\dfrac{d^{3}p}{(2\pi)^{3}}(4m_{i}^{2}+4p^{2})e^{-E/T}=4m_{i}^{2}+12m_{i}TK_{3}(\tfrac{m_{i}}{T})/K_{2}(\tfrac{m_{i}}{T}).$
(26)

## References

  * [1] Planck collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” (2018) , arXiv:1807.06209 [astro-ph.CO].
  * [2] J. M. Cline, K. Kainulainen, P. Scott, and C. Weniger, “Update on scalar singlet dark matter,” Phys. Rev. D 88 (2013) 055025, arXiv:1306.4710 [hep-ph]. [Erratum: Phys.Rev.D 92, 039906 (2015)].
  * [3] GAMBIT Collaboration, P. Athron et al., “Status of the scalar singlet dark matter model,” Eur. Phys. J. C 77 no. 8, (2017) 568, arXiv:1705.07931 [hep-ph].
  * [4] B. Dasgupta, E. Ma, and K. Tsumura, “Weakly interacting massive particle dark matter and radiative neutrino mass from Peccei-Quinn symmetry,” Phys. Rev. D 89 no. 4, (2014) 041702, arXiv:1308.4138 [hep-ph].
  * [5] A. Abada, D. Ghaffor, and S. Nasri, “A Two-Singlet Model for Light Cold Dark Matter,” Phys. Rev. D83 (2011) 095021, arXiv:1101.0365 [hep-ph].
  * [6] A. Ahriche, A. Arhrib, and S. Nasri, “Higgs Phenomenology in the Two-Singlet Model,” JHEP 02 (2014) 042, arXiv:1309.5615 [hep-ph].
  * [7] J. A. Casas, D. G. Cerdeño, J. M. Moreno, and J. Quilis, “Reopening the Higgs portal for single scalar dark matter,” JHEP 05 (2017) 036, arXiv:1701.08134 [hep-ph].
  * [8] V. R. Shajiee and A. Tofighi, “Electroweak Phase Transition, Gravitational Waves and Dark Matter in Two Scalar Singlet Extension of The Standard Model,” Eur. Phys. J. C79 no. 4, (2019) 360, arXiv:1811.09807 [hep-ph].
  * [9] A. Arhrib and M. Maniatis, “The two-real-singlet Dark Matter model,” Physics Letters B 796 (2019) 15 – 19, arXiv:1807.03554 [hep-ph].
  * [10] T. N. Maity and T. S. Ray, “Exchange driven freeze out of dark matter,” Phys. Rev. D 101 no. 10, (2020) 103013, arXiv:1908.10343 [hep-ph].
  * [11] M. Maniatis, “Light dark matter in a minimal extension with two additional real singlets,” Phys. Rev. D 103 no. 1, (2021) 015010, arXiv:2005.13443 [hep-ph].
  * [12] G. Belanger and J.-C. Park, “Assisted freeze-out,” JCAP 1203 (2012) 038, arXiv:1112.4491 [hep-ph].
  * [13] F. Bazzocchi and M. Fabbrichesi, “A simple inert model solves the little hierarchy problem and provides a dark matter candidate,” Eur. Phys. J. C 73 no. 2, (2013) 2303, arXiv:1207.0951 [hep-ph].
  * [14] K. P. Modak, D. Majumdar, and S. Rakshit, “A Possible Explanation of Low Energy $\gamma$-ray Excess from Galactic Centre and Fermi Bubble by a Dark Matter Model with Two Real Scalars,” JCAP 03 (2015) 011, arXiv:1312.7488 [hep-ph].
  * [15] S. Bhattacharya, P. Poulose, and P. Ghosh, “Multipartite Interacting Scalar Dark Matter in the light of updated LUX data,” JCAP 1704 no. 04, (2017) 043, arXiv:1607.08461 [hep-ph].
  * [16] M. Pandey, D. Majumdar, and K. P. Modak, “Two Component Feebly Interacting Massive Particle (FIMP) Dark Matter,” JCAP 1806 no. 06, (2018) 023, arXiv:1709.05955 [hep-ph].
  * [17] L. Bian, R. Ding, and B. Zhu, “Two Component Higgs-Portal Dark Matter,” Phys. Lett. B728 (2014) 105–113, arXiv:1308.3851 [hep-ph].
  * [18] S. Bhattacharya, A. Drozd, B. Grzadkowski, and J. Wudka, “Two-Component Dark Matter,” JHEP 10 (2013) 158, arXiv:1309.2986 [hep-ph].
  * [19] K. Agashe, Y. Cui, L. Necib, and J. Thaler, “(In)direct Detection of Boosted Dark Matter,” JCAP 10 (2014) 062, arXiv:1405.7370 [hep-ph].
  * [20] K. Kannike, “Vacuum Stability Conditions From Copositivity Criteria,” Eur. Phys. J. C 72 (2012) 2093, arXiv:1205.3781 [hep-ph].
  * [21] K. Kannike, “Vacuum Stability of a General Scalar Potential of a Few Fields,” Eur. Phys. J. C 76 no. 6, (2016) 324, arXiv:1603.02680 [hep-ph]. [Erratum: Eur.Phys.J.C 78, 355 (2018)].
  * [22] P. M. Ferreira, “The vacuum structure of the Higgs complex singlet-doublet model,” Phys. Rev. D 94 no. 9, (2016) 096011, arXiv:1607.06101 [hep-ph].
  * [23] B. W. Lee, C. Quigg, and H. B. Thacker, “Weak Interactions at Very High-Energies: The Role of the Higgs Boson Mass,” Phys. Rev. D 16 (1977) 1519.
  * [24] G. Cynolter, E. Lendvai, and G. Pocsik, “Note on unitarity constraints in a model for a singlet scalar dark matter candidate,” Acta Phys. Polon. B 36 (2005) 827–832, arXiv:hep-ph/0410102.
  * [25] L. Husdal, “On Effective Degrees of Freedom in the Early Universe,” Galaxies 4 no. 4, (2016) 78, arXiv:1609.04979 [astro-ph.CO].
  * [26] G. Belanger, F. Boudjema, and A. Pukhov, “micrOMEGAs : a code for the calculation of Dark Matter properties in generic models of particle interaction,” in Theoretical Advanced Study Institute in Elementary Particle Physics: The Dark Secrets of the Terascale, pp. 739–790.  2013\.  arXiv:1402.0787 [hep-ph].
  * [27] J. M. Cline and K. Kainulainen, “Electroweak baryogenesis and dark matter from a singlet Higgs,” JCAP 01 (2013) 012, arXiv:1210.4196 [hep-ph].
  * [28] J. e. a. Aprile, E. Aalbers, “Dark matter search results from a one ton-year exposure of xenon1t,” Physical Review Letters 121 no. 11, (Sep, 2018) . http://dx.doi.org/10.1103/PhysRevLett.121.111302.
  * [29] CMS Collaboration, A. M. Sirunyan et al., “Search for invisible decays of a Higgs boson produced through vector boson fusion in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Lett. B 793 (2019) 520–551, arXiv:1809.05937 [hep-ex].


that orbit-finitely spanned spaces are closed under duals. The space
$V\Rightarrow W$ is defined to be the space of finitely supported linear maps
from $V$ to $W$, with the natural vector space structure. This can be seen as
a subspace of a larger space, call it $V\Rrightarrow W$, which contains all
linear maps, not necessarily finitely supported. A standard result in linear
algebra is that $V\Rrightarrow W$ is linearly isomorphic to $V\otimes
W\Rrightarrow 1$, where $1$ is the field. This isomorphism is equivariant, and
it preserves the property of being finitely supported. Therefore, this gives
us an equivariant linear isomorphism between $V\Rightarrow W$ is isomorphic
and $V\otimes W\Rightarrow 1$. The latter space is the dual of the orbit-
finitely spanned space $V\otimes W$, and therefore it is orbit-finitely
spanned by [6, Corollary VI.5].

We would like to remark that a similar result as Theorem 6.23 was observed in
[przybyłek2024note, Theorem 3.8], but using the smaller category of vector
spaces that admit an orbit-finite basis.