compat=1.0.0

# Fermion Singlet Dark Matter in a Pseudoscalar Dark Matter Portal

Bastián Díaz Sáez<EMAIL_ADDRESS>Physik-Department, Technische
Universität München, James-Franck-Straße, 85748 Garching, Germany  Patricio
Escalona<EMAIL_ADDRESS>Departamento de Física, Universidad
Técnica Federico Santa María, Avenida España 1680, Valparaíso, Chile
Sebastián Norero<EMAIL_ADDRESS>Instituto de Física, Pontificia
Universidad Católica de Chile, Avenida Vicuña Mackenna 4860, Santiago, Chile
Alfonso Zerwekh<EMAIL_ADDRESS>Departamento de Física y Centro
Científico-Tecnológico de Valparaíso, Universidad Técnica Federico Santa
María, Avenida España 1680, Valparaíso, Chile

###### Abstract

We explore a simple extension to the Standard Model containing two gauge
singlets: a Dirac fermion and a real pseudoscalar. In some regions of the
parameter space both singlets are stable without the necessity of additional
symmetries, then becoming a possible two-component dark matter model. We study
the relic abundance production via freeze-out, with the latter determined by
annihilations, conversions and semi-annihilations. Experimental constraints
from invisible Higgs decay, dark matter relic abundance and direct/indirect
detection are studied. We found three viable regions of the parameter space,
and the model is sensitive to indirect searches.

###### Contents

  1. 1 Introduction
  2. 2 The model
  3. 3 Relic Density
    1. 3.1 Parameter Space
    2. 3.2 Boltzmann equation
      1. 3.2.1 $m_{s}<m_{h}$
      2. 3.2.2 $m_{s}>m_{h}$
  4. 4 Phenomenology
    1. 4.1 Experimental Constraints
    2. 4.2 Scan
    3. 4.3 Indirect detection
      1. 4.3.1 Box-shape gamma ray
      2. 4.3.2 Upper bounds
  5. 5 Discussion and Conclusions
  6. A Annihilation Cross Sections

##  1 Introduction

Astrophysical evidence of dark matter (DM) has been accumulating for more than
forty years now, but its fundamental nature remains unknown. From the particle
physics points of view, different approaches have been carried out over the
years to account for the elusive DM (for a review see [1]), and in particular,
simple extensions to the SM containing gauge singlets look appealing for their
simplicity, DM predictions and testable phenomenology [2, 3, 4, 5, 6, 7, 8].
Nowadays, those WIMP minimal extensions have been very constrained, especially
by direct detection, motivating other possible alternatives. For instance, the
interplay of a (pseudo)scalar and a fermion, both gauge singlets, open up the
possibilities in many aspects: multi-component DM, new interaction channels,
novel experimental signatures, small-scale structures, among others [9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29].

Keeping minimality, in this work we study a gauge singlet sector composed of a
real pseudoscalar and a Dirac fermion. Depending on the coupling values and
mass hierarchy between the singlets, the model admits a variety of DM
productions with either one or the two singlets being stable. Two new
interactions are present in the model: a Higgs portal and a dark sector
coupling. The Higgs interaction is key because it regulates to what extent the
dark sector is coupled to the SM, whereas the internal dark sector coupling
only regulates the coupling between the two singlets. In our knowledge, we
study for the first time the WIMP regime of this framework in which both
couplings take sizable values such that both singlets were in thermal
equilibrium with the SM bath in the early universe. Interestingly, in certain
mass hierarchy between the two singlets, the stability of both fields is
guaranteed by a parity symmetry without the necessity of introducing new ad
hoc discrete symmetries. Models with this last feature or accidental
symmetries have been studied in different DM context, such as Minimal Dark
Matter [30], spontaneous symmetry breaking [31], two DM components [32],
vector DM [33, 34] and rank-two fields [35].

In the model under consideration, the DM relic abundance is triggered by
annihilations, DM conversions [36] and semi-annihilations [37, 38], showing
remarkable features in some regions of the parameter space. Further, we
constraint the model considering the measured relic abundance in the universe,
Higgs invisible decay and direct/indirect detection. For the latter, we
explore box-shaped gamma ray spectra [39, 40], and we confront the available
parameter space with Fermi-LAT data, CTA projections and AMS-02 bounds.

The paper is organized as follows. In section 2 we present the model and its
theoretical constraints. In section 3 we explore the possible DM relic
abundance mechanisms presents in the model, with a precise analysis of the
two-component freeze-out scenario. In section 4 we review experimental
constraints and the available parameter space along with indirect detection
signals. In the last section we discuss and state our conclusions.

##  2 The model

The model adds to the SM two gauge singlets: one Dirac fermion $\psi$ and a
pseudo-scalar $s$. Under a parity transformation, the fields transform as
$\psi\rightarrow\gamma^{0}\psi$ and $s\rightarrow-s$, giving rise to the
following Lagrangian:

$\displaystyle\mathcal{L}=\mathcal{L}_{SM}+\bar{\psi}(i\not{\partial}-m_{\psi})\psi+\frac{1}{2}(\partial_{\mu}s)^{2}+ig_{\psi}s\bar{\psi}\gamma_{5}\psi-V(H,s),$
(1)

where the scalar potential is given by

$\displaystyle V(H,s)$ $\displaystyle=$
$\displaystyle\mu^{2}|H|^{2}+\lambda_{H}|H|^{4}+\frac{\mu_{s}^{2}}{2}s^{2}+\frac{\lambda_{s}}{4!}s^{4}+\frac{\lambda_{hs}}{2}|H|^{2}s^{2},$
(2)

with $H$ being the Higgs doublet111One may consider $\psi$ as a sterile
neutrino that mixes with the active ones. If the mixing is small enough, the
sterile neutrino may be stable on cosmological scales and can be produced
through active-sterile oscillations. In this work we assume that the mixing
angle is sufficiently small to avoid these effects. For a discussion in this
direction see [41, 22, 42].. We make $g_{\psi}$ real under a fermion field
redefinition such that the theory is CP-conserving. We assume that the singlet
scalar does not acquire vacuum expectation value (vev), and after EWSB in the
unitary gauge, $H=(0,(v_{H}+h)/\sqrt{2})^{T}$ with $v_{H}=246$ GeV the Higgs
vev, the scalar potential may be rewritten as

$\displaystyle V(h,s)$ $\displaystyle=$
$\displaystyle\frac{\mu_{s}^{2}}{2}s^{2}+\frac{\lambda_{s}}{4!}s^{4}+\frac{\lambda_{hs}}{4}(h^{2}+2v_{H}h)s^{2},$
(3)

with the mass of the scalars given by

$\displaystyle m_{h}^{2}$ $\displaystyle=$ $\displaystyle
2v_{H}^{2}\lambda_{H},$ (4) $\displaystyle m_{s}^{2}$ $\displaystyle=$
$\displaystyle\mu_{s}^{2}+\lambda_{hs}v_{H}^{2}/2.$ (5)

Here we consider $m_{h}=125$ GeV. The stability of the fermion is easily
recognized due to the fact that it appears in pairs in the Lagrangian. In the
scalar sector $s$ appears in pairs, forbidding its decay. The linear term in
$s$ in 1, imply that as $m_{s}\geq 2m_{\psi}$ the pseudoscalar may decay into
a pair of $\psi$, whereas as $m_{s}<2m_{\psi}$ the scalar singlet becomes
stable at tree level and at all orders in perturbation theory. In the
following we argument the latter fact.

Figure 1: General radiative $s$ decay. The values for $n^{\prime}$ refers to
internal $s$ lines in the closed fermionic loop, $n$ for outgoing $s$ lines
and $N$ for Higgs lines. The black globe takes into account arbitrary
interactions among $s$ and $h$ fields.

The possible decay of $s$ at an arbitrary number of loops is represented in
Fig. 1, with the decay of $s$ followed by a singlet fermion closed loop and
the black circle representing possible interactions between an arbitrary
number of $s$ and $h$. In the figure $n^{\prime},n$ and $N$ simply represent
the number of scalar lines. For simplicity, let us start assuming
$n^{\prime}=0$. If $n$ is even, the resulting fermion trace will at most
contain terms of the form
$\epsilon^{\mu\nu\rho\dots}p_{\mu}p_{\nu}p_{\rho}\dots$, which vanishes
exactly. If $n$ is odd, the trace is different than zero, but there is no way
to connect an odd number of outgoing $s$ with an arbitrary number of $N$ Higgs
boson in the black bloop, due to the presence of the CP symmetry of the scalar
potential. Now, if $n^{\prime}\neq 0$, the previous arguments still remains,
because internal lines of $s$ in the fermion loop only add an even number of
both $\gamma^{5}$ and fermion propagators to the trace, only adding vanishing
contributions without changing the final result. In consequence, from a
perturbative point of view, the stability is guaranteed for the pseudoscalar
singlet 222One can introduce an axion-like anomalous 5-dimensional effective
operator that conserves CP [43],
$\mathcal{L}\supset\frac{\lambda}{\Lambda}sG_{\mu\nu}\tilde{G}^{\mu\nu}$, with
$\lambda$ a dimensionless coupling, $\Lambda$ some high energy scale, $G$ the
field strength of any gauge field and $\tilde{G}$ its dual. This operator
induce the decay of the pseudoscalar singlet into gauge bosons. For instance,
considering $m_{s}=10^{3}\>\mathrm{GeV}$, we require $\lambda\lesssim 10^{-7}$
for a Planck scale induced operator in order to have a cosmologically stable
particle. For a GUT scale induced operator we require $\lambda\lesssim
10^{-10}$. Other possible low energy origin of such operator requires the
introduction of heavy vector-like fermions [44], which are not part of our
model construction..

Finally, theoretical constraints put bounds on the free parameters of the
model. The stability of the electroweak vacuum for $s$ and $h$ imposes that
[45]

$\displaystyle\lambda_{s}>0,\quad\lambda_{hs}>-\sqrt{\frac{2}{3}\lambda_{H}\lambda_{s}},$
(6)

with $\lambda_{H}=m_{h}^{2}/(2v_{H}^{2})\simeq 0.1$. From perturbativity we
set that $|g_{\psi}|,|\lambda_{s}|<4\pi$ to ensure that loop corrections are
smaller than tree-level processes [45]. Unitarity constraints are less
stringent than the upper limits based on our perturbativity criteria [46, 5].
The signs of $g_{\psi}$ and $\lambda_{hs}$ are not relevant for the analysis
in this work due to the fact that the relevant processes depend quadratically
on them, therefore the theoretical constraints set $0<\lambda_{hs}<4\pi$ and
$0<g_{\psi}<4\pi$.

Figure 2: Mass plane, with the green region indicating $s$ unstable, whereas
in the orange and blue regions both singlets are stable. The blue one
indicates where $s$-channel semi-annihilations are present, with $\Delta\equiv
2m_{\psi}-m_{s}-m_{h}\geq 0$ (for details see Sec. 3.2).

##  3 Relic Density

###  3.1 Parameter Space

Depending on the intensity of the couplings and the mass hierarchy between
$m_{\psi}$ and $m_{s}$, the model may present different DM scenarios with one
or two stable particles. When $\lambda_{hs}$ takes very small values, i.e.
$\lambda_{hs}\sim 10^{-12}-10^{-6}$, the singlet sector never thermalize with
the SM, then the DM production may occurs via freeze-in and/or dark freeze-out
[42, 47]. On the contrary, as the Higgs portal coupling takes sizable values,
the singlet scalar enters into thermal equilibrium with the SM (this work).
Based on the latter fact, two possible coupling regimes concerning $g_{\psi}$
may be present:

  * •

$g_{\psi}\lesssim 10^{-6}$: The two singlets in the dark sector will interact
feebly. If $m_{s}>2m_{\psi}$, $s$ becomes unstable, and $\psi$ may be produced
via freeze-in through the decay of $s$ and $2\rightarrow 2$ scattering
processes. The green region in Fig. 2 shows this parameter space in terms of
the mass hierarchy. For $m_{s}\leq 2m_{\psi}$, $s$ becomes stable and define a
DM candidate identical to that of the Singlet Higgs Portal model (SHP) [5, 8],
since $\psi$ does not interfere in the dynamic of the former due to their
feebly interactions. The SHP model has been exhaustively studied previously,
and displays a highly constrained parameter space around the EW scale.

  * •

$g_{\psi}\gtrsim 10^{-6}$: $s$ will bring $\psi$ fast into the thermal
equilibrium. From the relic abundance point of view, the only relevant case
here is when both singlets are stable, i.e. $m_{s}<2m_{\psi}$ (orange and blue
region in Fig. 2 ), otherwise $\psi$ would not have any channel to annihilate,
giving rise to an overabundance (see Feynman diagrams in Fig. 3). Based on the
latter point, different type of interactions appear that determine the relic
abundance of each singlet.

In this work, we focus on this last scenario, in which both singlets are
stable. It is worth to mention that even when the first case with one DM
candidate via freeze-in is a perfectly viable DM model, it presents a
challenge phenomenology [48] (see however [49]). Recently, novel scenarios
with inverse semi-production have been proposed with interesting dark matter
production and phenomenology [50, 51].

###  3.2 Boltzmann equation

$\psi$$\bar{\psi}$$s$$s$$h$(a)$\psi$$\psi$$s$$s$$h$(b)$\psi$$\psi$$\psi$$s$$s$(c)$\psi$$\psi$$\psi$$s$$s$(d)$s$$s$$h$SMSM(e)$s$$s$$s$$h$$h$(f)$s$$s$$s$$h$$h$(g)$s$$s$$h$$h$(h)
Figure 3: Relevant processes at the freeze-out time in the two DM component
model. Diagrams (a) and (b) are the $s$\- and $t$-channel semi-annihilations,
respectively, (c) and (d) correspond to DM conversions, and (e) to (h) are the
annihilations of the pseudoscalar into SM and Higgs particles.

As we pointed out before, the two-component scenario of our interest occurs
for $m_{s}<2m_{\psi}$ and when the couplings $(g_{\psi},\lambda_{hs})$ are
sizable. In this case, at temperatures higher than the individual masses of
the singlet, both DM components are in thermal equilibrium with the SM. The
departure of the equilibrium occurs once the temperature goes below the masses
of the singlets, and three types of scattering participate in this process:
annihilations, semi-annihilations, and dark matter conversions (Fig. 3). Based
on [52], the evolution of the individual singlet abundances $Y_{i}\equiv
n_{i}/s$, with $i=\psi,s$, as a function of the temperature $x\equiv\mu/T$,
with $\mu=m_{\psi}m_{s}/(m_{\psi}+m_{s})$, is given by

$\displaystyle\frac{dY_{\psi}}{dx}$ $\displaystyle=$
$\displaystyle-\lambda_{\psi\bar{\psi}ss}\left(Y_{\psi}^{2}-Y_{s}^{2}\frac{Y_{\psi,e}^{2}}{Y_{s,e}^{2}}\right)-\lambda_{\psi\bar{\psi}sh}\left(Y_{\psi}^{2}-Y_{s}\frac{Y_{\psi,e}^{2}}{Y_{s,e}}\right),$
$\displaystyle\frac{dY_{s}}{dx}$ $\displaystyle=$
$\displaystyle-\lambda_{ssXX}\left(Y_{s}^{2}-Y_{s,e}^{2}\right)+\lambda_{\psi\bar{\psi}ss}\left(Y_{\psi}^{2}-Y_{s}^{2}\frac{Y_{\psi,e}^{2}}{Y_{s,e}^{2}}\right)$
(7) $\displaystyle+$
$\displaystyle\frac{1}{2}\lambda_{\psi\bar{\psi}sh}\left(Y_{\psi}^{2}-Y_{s}\frac{Y_{\psi,e}^{2}}{Y_{s,e}}\right)-\frac{1}{2}\lambda_{s\psi\psi
h}\left(Y_{s}Y_{\psi}-Y_{\psi}Y_{s,e}\right),$

where we have defined

$\displaystyle\lambda_{ijkl}(x):=\dfrac{\langle\sigma_{ijkl}v\rangle(x)\cdot
s(T)}{x\cdot H(T)},\qquad\text{for}\leavevmode\nobreak\ i,j,k,l=\psi,s,h,X.$
(8)

with $X$ referring to a SM particle, $\braket{\sigma v}$ the thermally
averaged cross section, and the entropy density $s$ and Hubble rate $H$ in a
radiation dominated universe given by

$\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},$ (9)

where $G$ is the Newton gravitational constant, and $g_{*}$ and $g_{s*}$ are
the effective degrees of freedom contributing respectively to the energy and
the entropy density at temperature $T$ [53]. The equilibrium densities,
$Y_{i,e}\equiv n_{i,e}/s$, are calculated using the Maxwell-Boltzmann
distribution, whose number density is given by

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

with $g_{i}$ the internal spin degrees of freedom, and $K_{2}$ is the modified
Bessel function of the second kind. The DM density parameter is given by

$\displaystyle\Omega h^{2}=\sum_{i=\psi,s}\Omega_{i}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},$ (11)

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

In the following we analyze the behavior of the relic abundance obtained by
solving Eq. (3.2) making use of micrOMEGAs 5.2.7a [52], which calculates
automatically all $\braket{\sigma v}$ in Eq. 3.2. We pay special attention to
those regions in which $\Omega_{s}h^{2}$ acquires smaller values than in the
SHP [5] for a fixed $(m_{s},\lambda_{hs})$. This will be a necessary
requirement in order to evade stringent direct detection bounds on the
pseudoscalar DM from a few GeV to the TeV scale. Due to the subtle behavior of
semi-annihilations in the $t$-channel, we divide the analysis in two parts,
depending on the relative hierarchy between $m_{s}$ and $m_{h}$. Finally, for
convenience, we introduce a new parameter $\Delta\equiv
2m_{\psi}-m_{s}-m_{h}$, such that as $\Delta\geq 0$ $s$-channel semi-
annihilations are present in the relic abundance calculation. This is
represented in blue in Fig. 2 .

Figure 4: (above) Relic abundance for $m_{s}=40$ and $\lambda_{hs}=1$, for
different values of $g_{\psi}$ specified at the top of each plot. The green
and red region are forbidden by the instability of $s$ and relic
overabundance. The vertical dashed grey lines remark the points at which
$m_{\psi}=m_{s}$ and $\Delta=0$. (below) Relic abundance as a function of
$g_{\psi}$ for $m_{s}=40$ and $m_{\psi}=30,40$ and $100$ GeV, respectively,
with $\lambda_{hs}=1$.

####  3.2.1 $m_{s}<m_{h}$

In this region all the processes in Fig. 3 may be present but the semi-
annihilation in the $t$-channel. Now, when $m_{\psi}$ is small enough, i.e.
$\Delta<0$, then the $s$-channel semi-annihilation is not present either, and
only pseudoscalar annihilations and DM conversions are present. Normally in
this regime, an overabundance occurs as $m_{s}>m_{\psi}$, since $\psi$ does
not have effective annihilation channels. As $m_{\psi}>m_{s}$, DM conversions
of the type $\psi\bar{\psi}\rightarrow ss$ become effective, decreasing the
overabundance of $\Omega_{\psi}$. If $m_{\psi}$ is big enough such that
$\Delta\geq 0$, semi-annihilations in the $s$-channel become efficient,
decreasing $\Omega_{\psi}$ even more. These characteristics can be seen in
Fig. 4 (above), where the abundances are depicted as a function of $m_{\psi}$,
for $m_{s}=40$ GeV, $g_{\psi}=(0.1,1,10)$ and $\lambda_{hs}=1$, with the green
and red regions representing $s$ unstable and DM overabundance, respectively.
The effects of conversions at $m_{\psi}=m_{s}$ becomes sharper as $g_{\psi}$
increases, since $\braket{\sigma_{\psi\psi ss}v}\sim g_{\psi}^{4}$ (see Appx.
A). Equivalently, $s$-channel semi-annihilations near $\Delta=0$ pushes down
$\Omega_{\psi}$, and since both conversions and s-channel semi-annihilations
depends inversely on $m_{\psi}^{2}$, their effectiveness decreases with the
growing of $m_{\psi}$, then making $\Omega_{\psi}$ to increase for higher
values of $m_{\psi}$. The changes of $\Omega_{s}$ are stronger only for
$m_{s}>m_{\psi}$, where $ss\rightarrow\psi\bar{\psi}$ conversions are
effective, and due to $\braket{\sigma_{ss\psi\psi}v}\propto g_{\psi}^{4}$,
$\Omega_{s}$ decreases notoriously for such high couplings $g_{\psi}$ shown in
the Fig. 4 (above) right plot. We do not show the dependence of the relic
abundances on $\lambda_{hs}$, but the changes are minor due to the fact that
this parameter affects only the semi-annihilation intensity and modulates
$\Omega_{s}$.

In Fig. 4 (below) we show the abundances as a function of $g_{\psi}$, keeping
$\lambda_{hs}=1$. The relic abundance behavior depends strongly on the mass
hierarchy and the magnitude of $g_{\psi}$. As it was previously discussed, for
$m_{\psi}<m_{s}$ (left plot) it is $\Omega_{\psi}$ that dominates the relic,
independently of $g_{\psi}$. In this case, $\lambda_{hs}$ would only change
the relative abundance of the scalar singlet, but keeping very high values for
$\Omega_{\psi}$. Oppositely, in the cases $m_{\psi}>m_{s}$ (middle and right
plot), the relic hierarchy do depends on $g_{\psi}$, showing the effectiveness
of conversions and $s$-channel semi-annihilation, respectively, with a
notorious fall of $\Omega_{\psi}$ in each case.

####  3.2.2 $m_{s}>m_{h}$

In this case, the $t$-channel semi-annihilation $\psi+s\rightarrow\psi+h$ is
present, and it may participates strongly in the determination of $Y_{s}$ for
$\Delta<0$. The effectiveness of this semi-annihilation on $Y_{s}$ becomes
highly sharp, due to the fact that once both DM components decouple from the
SM plasma, $Y_{s}$ follows in a good approximation

$\displaystyle\frac{dY_{s}}{dx}\propto-\frac{1}{2}x^{-2}\lambda_{s\psi\psi
h}Y_{\psi}Y_{s},\quad\quad x\gtrsim 10,$ (12)

Figure 5: DM yields as a function of $x$ in the region $m_{s}>m_{h}$ and
$\Delta<0$, considering (continuous lines) and not (dashed lines) the
$t$-channel semi-annihilation. In both cases we have set $m_{s}=130$ GeV and
$g=\lambda_{hs}=1$, with $m_{\psi}=120$ GeV (left) and $m_{\psi}=110$ GeV
(right). In the legends are specified the parameter densities for each case.
The dotted lines correspond to the equilibrium densities of each singlet.

assuming $Y_{\psi}\approx$ constant and $Y_{s,e}\approx 0$. The solution of
Eq. (12) gives an exponential suppression for $Y_{s}$, highly sensitive to the
mass difference between the singlets. As an example of this behavior, in Fig.
5 we show the evolution of the densities $Y_{\psi}$ (blue lines) and $Y_{s}$
(red curves) as a function of $x=\mu/T$, for $m_{s}=130$ GeV, and
$m_{\psi}=120$ GeV (left plot) and 110 GeV (right plot). As Fig. 5(left)
suggests, $Y_{s}$ depends strongly on the $t$-channel semi-annihilation, as it
can be compared the dashed and solid red curves, with the former not
considering the process in the Boltzmann equation and the latter containing
it. As the mass difference between the two singlets increases, the effects
becomes sharper, as shown in Fig. 5(right). $Y_{\psi}$ is almost independent
on this process, as it can be seen through the overlapping of the dashed and
solid blue lines. This effect is particularly interesting due to the fact that
as $\Omega_{s}$ becomes negligible, no direct detection bounds will apply on
the pseudoscalar DM. This behavior was also seen in a two-component DM model
consisting of two complex scalars stabilized by a $Z_{5}$ symmetry [54].

For $\Delta>0$, the two semi-annihilations enter into the coupled Boltzmann
equations, and it is not longer possible to assume $Y_{\psi}$ constant,
therefore Eq. (12) is not a good approximation to determine $Y_{\psi}$. In
Fig. 6 we show the dependence of the relic abundance as a function of
$m_{\psi}$ (left) and $m_{s}$ (right). Contrary to the case in the low mass
regime of Fig. 4, in the left plot of Fig. 6 we observe that conversions occur
at higher values of $m_{\psi}$ than $s$-channel semi-annihilations. When the
latter opens up, the total relic abundance decreases in various order of
magnitude. There is an interesting effect we want to point out, namely the
fall of $\Omega_{s}$ as $m_{s}\gtrsim m_{\psi}$. In both plots of Fig. 6 is
possible to observe this effect, but in the right plot is more clear the role
that conversions are playing, with the red dashed line making a tiny well from
the soft growing as $m_{s}$ increases. This is, conversions of the type
$ss\rightarrow\psi\psi$ start to be effective making a deviation from the well
known behavior of the SHP [5] (for works with related behaviors see [36, 55,
56]). This type of wells, along with the Higgs resonance at $m_{s}\approx
m_{h}/2$, will be of particular importance to evade direct detection
constraints.

Figure 6: (left) Relic abundance as a function of $m_{\psi}$, for $m_{s}=$ 500
GeV, $g_{\psi}=10$ and $\lambda_{hs}=0.25$. (right) Relic abundance as a
function of $m_{s}$, for $m_{\psi}=500$ GeV, and the same couplings than in
the left plot. The blue region in the right corresponds to $\Delta<0$.

In conclusion, we have analyzed the effectiveness of annihilations,
conversions and semi-annihilations on the relic abundance of $\psi$ and $s$ in
the parameter space. A random scan under theoretical and experimental
constraints becomes useful in order to explore in major detail the viable
regions of the model, and this is what we analyze in Sec. 4.

##  4 Phenomenology

In this section we present the relevant experimental constraints on the model
along with a full scan of the parameter space. Once obtained the allowed
parameter space, we explore indirect signals and we set upper limits from it.

###  4.1 Experimental Constraints

When $m_{s}\leq m_{h}/2$, the Higgs boson can decay into two $s$ DM particles,
with a decay width given by

$\Gamma_{inv}(h\to ss)=\frac{\lambda^{2}_{hs}v_{H}^{2}}{32\pi
m_{h}}\sqrt{1-\frac{4m^{2}_{s}}{m^{2}_{h}}}.$ (13)

This contributes to the invisible branching rate
$\text{Br}(h\to\text{inv})=\Gamma_{\text{inv}}/(\Gamma_{\text{SM}}+\Gamma_{\text{inv}})$,
where the total decay width of the Higgs into SM particles is given by
$\Gamma_{\text{SM}}\leavevmode\nobreak\ =\leavevmode\nobreak\ 4.07$ MeV[57].
Experimental searches put stringent constraints on this quantity, and the most
strict value is given by $\text{Br}(h\to\text{inv})<0.19$ at $95\%$ C.L. [57].

From the DM sector we have the following constraints. First, the measurement
of the DM relic density today, given by $\Omega_{c}h^{2}=0.120\pm 0.001$ [58].
Based on the uncertainties in our computation, we apply this constraint with a
tolerance of $\sim\leavevmode\nobreak\ $10%, i.e.
$\Omega_{c}h^{2}\in[0.11,0.12]$. Secondly, direct detection, which in
multicomponent DM scenarios the interaction of DM with nucleus matter through
the spin-independent (SI) cross section comes from weighting it by a factor
which takes into account the relative abundance $\Omega_{i}$ over the Planck
measured $\Omega_{c}$, i.e.,

$\displaystyle\hat{\sigma}_{\text{SI,i}}=\left(\frac{\Omega_{i}}{\Omega_{\text{c}}}\right)\sigma_{\text{SI,i}},\quad
i=\psi,s,$ (14)

where $\sigma_{\text{SI,i}}$ is the spin-independent cross section. At tree
level, it is only the pseudoscalar DM which interacts with nuclei via the
Higgs portal, with its SI cross section given by [5]

$\sigma_{\text{SI,s}}=\frac{\lambda^{2}_{hs}f^{2}_{N}}{4\pi}\frac{\mu_{i}^{2}m^{2}_{n}}{m^{4}_{h}m_{s}^{2}},$
(15)

with $m_{n}$ as the nucleon mass, $f_{N}$ is a factor proportional to the
nucleon matrix elements, and $\mu\leavevmode\nobreak\ =\leavevmode\nobreak\
m_{n}m_{s}/(m_{n}+m_{s})$ is the DM-nucleon reduced mass. Upper bounds on
$\sigma_{\text{SI}}$ are given by XENON1T data [59] and the projections of
XENONnT [60].

###  4.2 Scan

In this section we show a scan on the parameter space in the range
$m_{\psi,s}\in[10,1000]$ GeV (subject to $m_{s}<2m_{\psi}$) and
$g_{\psi},\lambda_{hs}\in[0.001,4\pi]$. The full scan highlight interesting
features that were already anticipated in Sec. 3.2. In Fig. 7 we show the scan
of points projected on the plane $(m_{\psi},m_{s})$, with the density color
representing the total abundance $\Omega h^{2}$ (left), $\Omega_{\psi}h^{2}$
(middle) and $\Omega_{s}h^{2}$, respectively. The green region (top left)
corresponds to those points with $m_{s}>2m_{\psi}$, then making $s$ unstable.
We pay attention to four regions based on Fig. 7 (left):

Figure 7: Random scan points projected in the mass plane with the bar color
indicating the total relic abundance $\Omega h^{2}$ (left),
$\Omega_{\psi}h^{2}$ (middle) and $\Omega_{s}h^{2}$ (right). The homogeneous
green region in the top left of each plot corresponds to $s$ unstable, i.e.
$m_{s}>2m_{\psi}$. The red dashed line indicates $m_{\psi}=m_{s}$ and the
yellow one the frontier $\Delta=0$.

  1. (a)

Dark red band. In general, $\Omega_{\psi}$ get too big values for
$m_{\psi}<m_{s}$, since no effective annihilation channels for the singlet
fermion are present (see Sec. 3.2.1). However, there is a smooth relic density
transition from the deep red to lower relic densities, due to the thermal tail
distribution of conversions and $s$-channel semi-annihilation as
$\Delta\lesssim 0$. With respect to $\Omega_{s}$, it is always sub-abundant,
specially for $m_{s}>m_{h}$ where the $t$-channel semi-annihilations are
present, and as we have described in Sec. 3.2.2, tiny mass shift between
$\psi$ and $s$ makes $\Omega_{s}$ to decrease strongly (deep blue region in
Fig. 7 (right)).

  2. (b)

Green area in below. In this region the relic is mostly dominated by the
pseudoscalar, as it can be seen from the plots. As it was shown in Fig. 4,
$\Omega_{\psi}<\Omega_{s}$ for conversion-driven processes of the type
$\psi\bar{\psi}\rightarrow ss$ and moderate couplings ($g_{\psi}\sim 1$). For
very high $m_{\psi}$ is possible to observe that $\Omega_{\psi}$ tend to
increase, as conversions and semi-annihilations are less effective.

  3. (c)

Higgs resonance and thresholds. In the region 50 GeV $\lesssim m_{s}\leq
m_{h}$ the Higgs resonance and SM thresholds are present, then pseudoscalar
annihilations are enhanced. These effects are more clear in Fig. 7 (right) for
$\Omega_{s}$. For $m_{s}\leavevmode\nobreak\ \lesssim\leavevmode\nobreak\
m_{h}/2$, $s$-channel annihilations of the singlet pseudoscalar via Higgs
resonance tend to be very effective, decreasing $\Omega_{s}$ substantially, as
it can be seen with the deep blue horizontal line. For $m_{h}/2<m_{s}<m_{W}$
the abundance of $s$ lift up, and for $m_{s}\lesssim m_{W}$ the relic
decreases again. The thresholds at $m_{Z}$ and $m_{h}$ are also present but
they are less notorious.

  4. (d)

Blue-green upper region: For the singlet masses $\gtrsim$ 100 GeV, the relic
of both components tends to be low (blue points), but as the masses increase
the effectiveness of conversions and semi-annihilations become less strong,
i.e. $\braket{\sigma v}\sim m_{\psi,s}^{-2}$ (Appx. A), then raising the relic
of both components. In this mass region the relic abundance seems to acquire
moderate values in order to fulfill the Planck measurement.

The conjunction of the three plots suggests that most of the points that could
give the correct relic abundance are just in very specific sectors, and those
regions must to shrink even more when other constraints are taken into
account. In Fig. 8 (left) we show the resulting parameter space points after
imposing the constraints from relic abundance, Higgs to invisible and direct
detection bounds. Specifically, the regions are: (i) Higgs resonance, (ii)
region with $m_{s}>m_{h}$ and $\Delta<0$, and (iii) in the high mass with
$\Delta>0$ and $m_{\psi}<m_{s}$, indicating the corresponding $\Omega_{\psi}$
with a color bar. Region (i) is the only one that allows $m_{\psi}>m_{s}$ and
where the Higgs to invisible constraint applies 333In the presence of a
resonance, the smallness of the Higgs portal coupling may breakdown the
assumption of local thermal equilibrium when chemical decoupling is taking
place. In our case we have not considered those effects, where a more involved
treatment must be carried [61].. The points in region (ii) easily evade
XENON1T constraints due to the fact that $\Omega_{s}\rightarrow 0$, then the
condition $\Omega_{\psi}=\Omega_{c}$ is enough to fulfill all the experimental
constraints. Finally, in the high mass region we found some points anticipated
by the analysis around Fig. 6, where the power of conversions and $s$-channel
semi-annihilations makes $\Omega_{s}$ to decrease enough evading XENON1T. The
cost of the latter implies values for $g_{\psi}$ near the perturbativity limit
criteria (this point is discussed in the last section). In Fig. 8 (right) we
project the scan points on the plane $(m_{s},\hat{\sigma}_{SI,s})$, contrasted
with the limits given by XENON1T and the XENONnT projection [60]. Even when
XENON experiments rule out most of the orange points (i.e. regions (i) and
(iii), a fraction of red points ($\Delta<0$ and $m_{s}>m_{h}$) easily evade
the strongest direct detection constraints. It is worth to say that this
inverted peak at the Higgs resonance is a typical characteristic of Higgs
portals, hardly to be ruled out by experiments.

Figure 8: (left) Points fulfilling direct detection and relic abundance in all
the scanned parameter space. (right) The same selected random points shown in
the left plot, now projected on $(m_{s},\hat{\sigma}_{SI,s})$, where we
distinguish points with $\Delta>0$ (orange) and $\Delta<0$ (red). The
continuous blue line corresponds to the upper limit given by XENON1T (1
t$\cdot$y)[59], whereas the green dashed line a projection for XENONnT (20
t$\cdot$y) [60].

To end this subsection, we would like to point out about the limitations of
our random scan, which was based on overlaying exclusion limits (mainly from
relic abundance and direct detection). It is well known that this type of
exclusions present some drawbacks, such as the blindness to the influence of
the choice of some parameters (e.g. DM halo distribution or the Higgs mass
pole) on the experimental limits, or the fact that the allowed regions do not
present additional information about which points are more favored than
others. Considering those limitations, statistical analysis based on combined
likelihoods functions may give a more accurate and realistic way to cure those
disadvantages, although beyond the scope of the present analysis [8, 62].

###  4.3 Indirect detection

In the previous subsection we found that three regions of the parameter space
fulfill invisible Higgs decay, relic abundance and direct detection
constraints, with two of them being practically in the region $\Delta>0$. In
that case, $s$-channel semi-annihilations are expected to give sizable fluxes
of particles today through its $s$-wave nature. The $t$-channel semi-
annihilation also goes in the $s$-wave, but we have checked that in general
its cross section is lower than the corresponding $s$-channel, and indirect
signals with a pair of $s$ in the initial state have shown to be not sizable
outside the Higgs resonance [5]. Additionally, $\Omega_{\psi}h^{2}$ reaches
$\sim 50\%$ to $99\%$ of the total DM budget in the high mass regime (see Fig.
8(left)), consequently the scaling factor $(\Omega_{\psi}/\Omega_{c})^{2}$
does not considerable suppress the flux produced by the $s$-channel semi-
annihilation. In the following, we show the box-shaped differential spectra
[63] for this channel with its subsequent decay $h\rightarrow\gamma\gamma$,
and secondly, restrictions on the parameter space coming from bounds based on
searches of gamma rays (Fermi-LAT), anti-protons (AMS-02), and projections
from the Cherenkov Telescope Array (CTA) are presented.

####  4.3.1 Box-shape gamma ray

The differential flux of photons produced in fermionic DM annihilations and
received at earth from a given solid angle in the sky $\Delta\Omega$ with a
detector of area $A$ is given by

$\displaystyle\frac{d\Phi_{\gamma}}{dE_{\gamma}}=\frac{1}{A}\frac{dN_{\gamma}}{dE_{\gamma}dt}=\frac{\braket{\sigma_{sh}v}}{8\pi
m_{\psi}^{2}}\left(\text{Br}_{h}\frac{dN}{dE_{\gamma}}\right)\left(\frac{1}{\Delta\Omega}\int_{\Delta\Omega}Jd\Omega\right),$
(16)

where $\braket{\sigma_{sh}v}$ is the corresponding average annihilation cross
section times velocity of the process ${\psi+\bar{\psi}\rightarrow s+h}$,
$\text{Br}_{h}$ is the branching ratio of the Higgs into two photons, and
$dN/dE_{\gamma}$ the corresponding normalized spectra. The $J$ factor is the
integral of the squared DM density $\rho_{DM}$ along the line of sight
$J=\int_{\text{l.o.s}}ds\rho^{2}_{DM}$. We consider as our main region of
interest the galactic center, which features $\Delta\Omega=1.30$ sr,
$\int_{\Delta\Omega}Jd\Omega=9.2\times 10^{22}$ GeV2cm-5, assuming a NFW
profile normalized to a local DM density of 0.4 GeV/cm3 [64]. To determine the
normalized spectra of emitted photons, we first note that their energies, in
the fermion DM collision frame444Today, fermion DM $\psi$ moves non-
relativistically, therefore colliding practically in the earth rest frame.,
are given by

$\displaystyle
E_{\gamma,1}=\frac{m_{h}^{2}/2}{E_{h}-\sqrt{E_{h}^{2}-m_{h}^{2}}\cos\theta},\quad\quad
E_{\gamma,2}=E_{h}-E_{\gamma,1},$ (17)

where $\theta$ corresponds to the angle sustained by one of the photons and
the in-flight Higgs. The energy of the emitted scalar particles in the rest
frame of the collision of $\psi$ are given by

$\displaystyle
E_{h}=m_{\psi}\left(1-\frac{m_{s}^{2}-m_{h}^{2}}{4m_{\psi}^{2}}\right),\quad
E_{s}=m_{\psi}\left(1+\frac{m_{s}^{2}-m_{h}^{2}}{4m_{\psi}^{2}}\right).$ (18)

For a fixed $E_{h}$, the energy of the photon received at earth depends only
on $\theta$, with a maximum (minimum) energy for $\theta=0(\pi/2)$ 555It is
possible to receive the two photons at earth as $\theta\rightarrow 0$, or
equivalently, when both photons are emitted transverse to the Higgs direction.
This scenario requires a highly boosted Higgs and a detailed analysis of the
emitted photons. then displaying a box-shaped spectrum centered at
$E_{c}\equiv(E(0)+E(\pi))/2=E_{h}/2$ and with a width $\Delta E\equiv
E(0)-E(\pi)=\sqrt{E_{h}^{2}-m_{h}^{2}}$. Therefore, the normalized spectra can
be written as 666Note that in the $\psi$ annihilation center-of-mass frame,
due to the conservation of angular momentum, the two photons must be emitted
back to back if they have the same polarization, and co-linearly if they have
opposite polarization; the conservation of linear momentum requires $s$ to be
emitted along with one of the photons in the former case, and in the direction
opposite to the photons in the latter. See also [65]..

$\displaystyle\frac{dN}{dE_{\gamma}}=\frac{2}{\Delta
E}\Theta\left(E_{\gamma}-E_{c}+\frac{1}{2}\Delta
E\right)\Theta\left(E_{c}+\frac{1}{2}\Delta E-E_{\gamma}\right),$ (19)

where the factor multiplicative factor 2 accounts for the two emitted photons
(then $E_{\gamma}$ any of the two photons in 17), and $\Theta$ are Heaviside
functions.

In Fig. 9 (left) we show the box-shaped spectra for two random points allowed
by the relic density abundance and XENON1T, contrasted with the data given by
Fermi-LAT (red dashed line), where we have taken the background fitting
function in order to estimate the data signal: $d\Phi_{\gamma}/dE=2.4\times
10^{-5}(E_{\gamma}/$GeV$)^{-2.55}$ GeV-1cm-2s-1sr-1 [64]. As it was
anticipated above, the differential fluxes for the selected points are well
below the red dashed line, in principle not showing any possible tension. The
blue and orange dot curves show the spectra when the two DM components become
degenerated, widening the boxes and therefore peaking at higher energy
photons. Only points with masses (and high couplings) as small as
$m_{\psi}\sim 100$ GeV may approach significantly to the red dashed line, as
it is depicted by the grey dashed line in Fig. 9 (left). It is known that the
power of box-shaped gamma rays goes in their potential deviation with respect
to the power-law Fermi-LAT background, then the previous result is premature
to discard a possible constraint on the parameter space. In order to do this,
in the following we use the recent upper bounds on gamma rays produced in
semi-annihilations based on Fermi-LAT data, along with anti-protons flux
measurements.

Figure 9: Differential gamma ray flux originated from the process
$\psi+\bar{\psi}\rightarrow s+h$, followed by the decay
$h\rightarrow\gamma\gamma$, for two random points fulfilling relic abundance
and direct detection (continuous lines), whereas the dotted curves represent
hypothetical points in which the two DM components are completely degenerated.
The dashed red line corresponds to the flux given by Fermi-LAT [63], and the
gray an hypothetical point with low $(m_{\psi},m_{s})$.

####  4.3.2 Upper bounds

Upper limits on $\braket{\sigma_{hs}v}$ have not been constructed yet in order
to constraint the parameter space via indirect searches. However, they can be
inferred from existing upper bounds on the average annihilation cross section
times velocity for $DM+DM\rightarrow DM+h$ process based on gamma-ray signals
[66], and from $DM+DM\rightarrow b\bar{b}(W^{+}W^{-})$ based on anti-proton
flux measurements [67]. In the following we derive the useful algebraic
relations to translate the existing upper bounds to our average semi-
annihilation cross section times velocity $\braket{\sigma_{sh}v}$.

First, let us consider the process $\psi\bar{\psi}\rightarrow\psi^{\prime}h$,
with $\psi$ and $\psi^{\prime}$ arbitrary DM particles, and $h$ the Higgs
boson. This process presents an average cross section times relative velocity
given by

$\displaystyle\braket{\sigma_{\psi^{\prime}h}v}=\frac{16\pi}{J}m_{\psi}^{2}\Phi_{h},\quad\quad\text{with}\quad\quad
E_{h}=m_{\psi}\left(1-\frac{m_{\psi^{\prime}}^{2}-m_{h}^{2}}{4m_{\psi}^{2}}\right)$
(20)

with $J$ an arbitrary $J$-factor, $m_{\psi}$ the mass of the initial states,
$\Phi_{h}$ is the flux, and the outgoing Higgs having an energy $E_{h}$
(equivalently to 18). Now, it is possible to have the same flux $\Phi_{h}$
with an outgoing Higgs having the same energy $E_{h}$ in another process given
by $\psi\bar{\psi}\rightarrow hh$:

$\displaystyle\Phi_{h}=\frac{1}{2}\frac{\braket{\sigma_{hh}v}}{8\pi
E_{h}^{2}}J,$ (21)

where the 1/2 factor comes from the fact that we are considering one outgoing
Higgs, and we have set $m_{\psi}=E_{h}$. Combining 20 and 21 we obtain

Figure 10: Points with $\Delta>0$ fulfilling direct detection constraints as a
function of the re-scaled average annihilation cross section times velocity
given by expressions 23 and 26 multiplied by the relative abundance
$(\Omega_{\psi}/\Omega_{c})^{2}$ of the annihilating fermion DM. Upper bounds
for Fermi-LAT [66], CTA [66] and AMS-02 [67] are given by the red continuous
lines. For details see the main text.

$\displaystyle\braket{\sigma_{\psi^{\prime}h}v}=\frac{1}{(1-\xi_{\psi^{\prime}})^{2}}\braket{\sigma_{hh}v},$
(22)

with $\xi_{x}\equiv\left(m_{x}^{2}-m_{h}^{2}\right)/(4m_{\psi}^{2})$, with $x$
being some particle in the final state. From the last relation $\psi^{\prime}$
is an arbitrary state, then in eq. 22 taking $\psi^{\prime}=\psi$ and
$\psi^{\prime}=s$ and combining the resulting expressions, it follows that

$\displaystyle\left(\frac{1-\xi_{s}}{1-\xi_{\psi}}\right)^{2}\braket{\sigma_{sh}v}=\braket{\sigma_{\psi
h}v},$ (23)

Upper limits on $\braket{\sigma_{\psi h}v}$ [66], now can be compared to our
$\braket{\sigma_{sh}v}$ for certain values of $m_{\psi}$ and $m_{s}$ through
eq. 23. In Fig. 10 (top row) we show the projection of under-abundant points
(green) and those that give approximately the correct relic abundance (blue
points) for the resulting average annihilation cross section times velocity
(left side of eq. 23) times the relative abundance of the annihilating DM
particles. The red line in both plots corresponds to the upper bound found in
[66] (right side of eq. 23). All these points are in the region $\Delta>0$ and
fulfill direct detection constraints. None of the bounds touch the parameter
space points, then not showing any possible exclusion. Note that the blue
points well below the upper bounds (red line) correspond to the Higgs
resonance region, since lower couplings $(g_{\psi},\lambda_{hs})$ are
necessary to compensate the resonance effect, then reducing the ID signals.

On the other hand, AMS-02 anti-protons measurement set upper bounds on
$\braket{\sigma_{\text{DM DM}\rightarrow
X\bar{X}}v}\equiv\braket{\sigma_{X\bar{X}}v}$ with
$X\bar{X}=b\bar{b},W^{+}W^{-}$ [67]. From the semi-annihilation with the Higgs
decaying into $b\bar{b}$ we have that

$\displaystyle\text{Br}(h\rightarrow
b\bar{b})\braket{\sigma_{sh}v}=\frac{8\pi}{J}m_{\psi}^{2}\Phi_{b\bar{b}}$ (24)

Additionally, the flux $\Phi_{b\bar{b}}$ in 24 can be produced by an arbitrary
interaction $\psi\bar{\psi}\rightarrow b\bar{b}$:

$\displaystyle\Phi_{b\bar{b}}=\frac{\braket{\sigma_{b\bar{b}}v}}{8\pi
m_{\psi}^{2}}J$ (25)

Combining 24 and 25 we obtain

$\displaystyle\text{Br}(h\rightarrow
b\bar{b})\braket{\sigma_{sh}v}=\braket{\sigma_{b\bar{b}}v}$ (26)

Equivalently for $W^{+}W^{-}$. In Fig. 9 (bottom row) we show the projection
of points considering the semi-annihilation with Higgs decay into $b\bar{b}$
and $W^{+}W^{-}$, contrasted with the upper bounds given by AMS-02. From the
resulting plots, $b\bar{b}$ search set the stringent bounds on the parameter
space of the model, discarding points with the correct relic abundance with
masses below $\sim 500-700$ GeV, approximately.

##  5 Discussion and Conclusions

In this paper we have explored a simple extension to the SM containing two
gauge-singlet fields: a Dirac fermion and a real pseudoscalar. From the DM
point of view, the model present multiple scenarios with one or two DM
components. Previous works have focused on the scenario when the Higgs portal
coupling takes very small values, in such a way that the dark sector evolves
decoupled from the SM, with the DM relic abundance produced via dark freeze-
out/freeze-in. In the present work we found that as the Higgs portal take
sizable values it is possible to produce a one-component DM candidate via
freeze-in or a two-component DM via freeze-out. We have focused on the latter
scenario, which only presents four free parameters: two masses and two
couplings.

The stability of the two singlets is guaranteed by parity arguments without
the necessity of invoking additional symmetries. Furthermore, the relic
abundance of both singlets is determined via freeze-out through annihilations,
DM conversions and semi-annihilations. The appearance of the latter is
contrary to the standard belief that this type of processes appear only in the
presence of symmetries larger than $Z_{2}$. We have explored the relic density
abundance of both DM components in the parameter space, founding interesting
behaviors depending on the mass hierarchy and couplings values. Semi-
annihilations and DM conversions play an important role in two of the three
available regions, making the pseudoscalar relic abundance low enough to evade
direct detection XENON1T bounds, and then allowing DM with masses of hundreds
of GeV up to the TeV scale. As it is usual with Higgs portals, the resonance
region will not be discarded completely, even with the powerful projection of
XENONnT.

We have complemented our analysis with indirect detection signatures and
bounds. Firstly, we have explored the box-shape gamma rays signals appearing
from a fermion DM semi-annihilation, showing small boxes signals for the
fermion DM $s$-channel semi-annihilation. Furthermore, we have translated
bounds from gamma-ray searches from Fermi-LAT and CTA projections onto the
semi-annihilation, not showing any possible tension, although the sensitivity
of both experiments rely in the ballpark of part of the parameter space of the
model. In fact, considering that our numerical methods are not very precise, a
more detailed and sophisticated analysis such as a global-fit could be
perfectly sensitive to these gamma-rays upper bounds in view of the expected
quantitative differences between the two methods. As a final analysis, we have
tested the model with anti-protons upper bounds from AMS-02 experiment,
showing an exclusion of fermion DM masses below $\sim 500-700$ GeV.

Finally, considering that the viability of the model requires high values for
the dark sector coupling $g_{\psi}$ in some regions of the parameter space,
there are two important points to be taken into account, although the precise
calculation of them are beyond the scope of the present work. The first is
related to the possible appearance of Landau poles at low energy scales,
implying an urgent UV completion of the model (e.g. embedding the pseudoscalar
into a complex singlet which acquires vev). The second point is related to the
spin-independent one-loop amplitude in direct detection for the singlet
fermion, which due to the high coupling values it may give a sizable
contribution, possibly excluding even more parameter space of the model, in
particular the special region in which the relic abundance of the pseudoscalar
drops to zero.

## Acknowledgments

We would like to thank Alejandro Ibarra, Camilo García, Alexander Belyaev and
Claudio Dib for useful discussions, and to Alexander Pukhov for helping with
MicrOMEGAS. B.D.S would like also to ANID (ex CONICYT) Grant No. 74200120.
A.Z. has been partially founded by ANID (Chile) PIA/APOYO AFB 180002, and P.
E. has been partially founded by project FONDECYT N° 1170171.

##  Appendix A Annihilation Cross Sections

The exact values of $\braket{\sigma v}$ for the different $2\rightarrow 2$
processes were calculated with micrOMEGAs 5.2.7a. In order to show the
dependence of the thermally average cross section on the parameters of the
model, in this appendix we show some expressions for $\braket{\sigma v}$ in
the limit in which $\braket{\sigma v}\approx(\sigma
v)|_{s=(m_{1}+m_{2})(1+v^{2}/4)}$, with $m_{1}$ and $m_{2}$ the masses of the
annihilating particles, keeping the $s$-wave only when it is present:

$\displaystyle\braket{\sigma v}_{\psi+\bar{\psi}\rightarrow s+h}$
$\displaystyle=$
$\displaystyle\frac{g_{\psi}^{2}\lambda_{hs}^{2}v_{H}^{2}\sqrt{-2m_{h}^{2}(m_{s}^{2}+4m_{\psi}^{2})+m_{h}^{4}+(m_{s}^{2}-4m_{\psi}^{2})^{2}}}{64\pi
m_{\psi}^{2}(4m_{\psi}^{2}-m_{s}^{2})^{2}},$ (27) $\displaystyle\braket{\sigma
v}_{\psi+\bar{\psi}\rightarrow s+s}$ $\displaystyle=$
$\displaystyle\frac{g_{\psi}^{4}m_{\psi}(m_{\psi}^{2}-m_{s}^{2})^{5/2}}{24\pi(m_{s}^{2}-2m_{\psi}^{2})^{4}}v^{2}$
(28) $\displaystyle\braket{\sigma v}_{s+s\rightarrow\psi+\bar{\psi}}$
$\displaystyle=$
$\displaystyle\frac{g_{\psi}^{4}(m_{s}^{2}-m_{\psi}^{2})^{3/2}(2m_{s}^{2}+3m_{\psi}^{2})}{60\pi
m_{s}^{7}}v^{4}$ (29) $\displaystyle\braket{\sigma
v}_{\psi+s\rightarrow\psi+h}$ $\displaystyle=$
$\displaystyle-\frac{g_{\psi}^{2}\lambda_{hs}^{2}v_{H}^{2}\sqrt{(m_{s}^{2}-m_{h}^{2})((m_{s}+2m_{\psi})^{2}-m_{h}^{2})}\left(2m_{\psi}-\frac{-m_{h}^{2}+m_{s}^{2}+2m_{s}m_{\psi}+2m_{\psi}^{2}}{m_{s}+m_{\psi}}\right)}{32\pi
m_{s}(m_{s}+m_{\psi})^{2}\left(m_{\psi}\frac{-m_{h}^{2}+m_{s}^{2}+2m_{s}m_{\psi}+2m_{\psi}^{2}}{m_{s}+m_{\psi}}+m_{s}^{2}-2m_{\psi}^{2}\right)^{2}}$
(30)

The corresponding expressions for the annihilations of a pair of $s$ into SM
particles can be found in the appendix of [5].

## References

  * [1] G. Bertone and D. Hooper, “History of dark matter,” Rev. Mod. Phys. 90 (Oct, 2018) 045002. https://link.aps.org/doi/10.1103/RevModPhys.90.045002.
  * [2] V. Silveira and A. Zee, “SCALAR PHANTOMS,” Phys. Lett. B 161 (1985) 136–140.
  * [3] J. McDonald, “Gauge singlet scalars as cold dark matter,” Phys. Rev. D 50 (1994) 3637–3649, arXiv:hep-ph/0702143.
  * [4] C. P. Burgess, M. Pospelov, and T. ter Veldhuis, “The Minimal model of nonbaryonic dark matter: A Singlet scalar,” Nucl. Phys. B 619 (2001) 709–728, arXiv:hep-ph/0011335.
  * [5] 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)].
  * [6] Y. G. Kim, K. Y. Lee, and S. Shin, “Singlet fermionic dark matter,” JHEP 05 (2008) 100, arXiv:0803.2932 [hep-ph].
  * [7] M. Escudero, A. Berlin, D. Hooper, and M.-X. Lin, “Toward (Finally!) Ruling Out Z and Higgs Mediated Dark Matter Models,” JCAP 12 (2016) 029, arXiv:1609.09079 [hep-ph].
  * [8] 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].
  * [9] M. Pospelov, A. Ritz, and M. B. Voloshin, “Secluded WIMP Dark Matter,” Phys. Lett. B 662 (2008) 53–61, arXiv:0711.4866 [hep-ph].
  * [10] Q.-H. Cao, E. Ma, J. Wudka, and C. P. Yuan, “Multipartite dark matter,” arXiv:0711.3881 [hep-ph].
  * [11] K. M. Zurek, “Multi-Component Dark Matter,” Phys. Rev. D 79 (2009) 115002, arXiv:0811.4429 [hep-ph].
  * [12] S. Baek, P. Ko, and W.-I. Park, “Search for the Higgs portal to a singlet fermionic dark matter at the LHC,” JHEP 02 (2012) 047, arXiv:1112.1847 [hep-ph].
  * [13] L. Lopez-Honorez, T. Schwetz, and J. Zupan, “Higgs portal, fermionic dark matter, and a Standard Model like Higgs at 125 GeV,” Phys. Lett. B 716 (2012) 179–185, arXiv:1203.2064 [hep-ph].
  * [14] M. Heikinheimo, A. Racioppi, M. Raidal, C. Spethmann, and K. Tuominen, “Dark Supersymmetry,” Nucl. Phys. B 876 (2013) 201–214, arXiv:1305.4182 [hep-ph].
  * [15] K. Ghorbani, “Fermionic dark matter with pseudo-scalar Yukawa interaction,” JCAP 01 (2015) 015, arXiv:1408.4929 [hep-ph].
  * [16] Y. G. Kim, K. Y. Lee, C. B. Park, and S. Shin, “Secluded singlet fermionic dark matter driven by the Fermi gamma-ray excess,” Phys. Rev. D 93 no. 7, (2016) 075023, arXiv:1601.05089 [hep-ph].
  * [17] S. Bhattacharya, A. Drozd, B. Grzadkowski, and J. Wudka, “Two-Component Dark Matter,” JHEP 10 (2013) 158, arXiv:1309.2986 [hep-ph].
  * [18] S. Esch, M. Klasen, and C. E. Yaguna, “Detection prospects of singlet fermionic dark matter,” Phys. Rev. D 88 (2013) 075017, arXiv:1308.0951 [hep-ph].
  * [19] S. Esch, M. Klasen, and C. E. Yaguna, “A minimal model for two-component dark matter,” JHEP 09 (2014) 108, arXiv:1406.0617 [hep-ph].
  * [20] S. Ipek, D. McKeen, and A. E. Nelson, “A Renormalizable Model for the Galactic Center Gamma Ray Excess from Dark Matter Annihilation,” Phys. Rev. D 90 no. 5, (2014) 055021, arXiv:1404.3716 [hep-ph].
  * [21] Y. Cai and A. P. Spray, “Fermionic Semi-Annihilating Dark Matter,” JHEP 01 (2016) 087, arXiv:1509.08481 [hep-ph].
  * [22] J. König, A. Merle, and M. Totzauer, “keV Sterile Neutrino Dark Matter from Singlet Scalar Decays: The Most General Case,” JCAP 11 (2016) 038, arXiv:1609.01289 [hep-ph].
  * [23] A. Ahmed, M. Duch, B. Grzadkowski, and M. Iglicki, “Multi-Component Dark Matter: the vector and fermion case,” Eur. Phys. J. C 78 no. 11, (2018) 905, arXiv:1710.01853 [hep-ph].
  * [24] F. Kahlhoefer, K. Schmidt-Hoberg, and S. Wild, “Dark matter self-interactions from a general spin-0 mediator,” JCAP 08 (2017) 003, arXiv:1704.02149 [hep-ph].
  * [25] S. Baek, P. Ko, and J. Li, “Minimal renormalizable simplified dark matter model with a pseudoscalar mediator,” Phys. Rev. D 95 no. 7, (2017) 075011, arXiv:1701.04131 [hep-ph].
  * [26] G. Arcadi, M. Lindner, F. S. Queiroz, W. Rodejohann, and S. Vogl, “Pseudoscalar Mediators: A WIMP model at the Neutrino Floor,” JCAP 03 (2018) 042, arXiv:1711.02110 [hep-ph].
  * [27] K. Ghorbani and P. H. Ghorbani, “Leading Loop Effects in Pseudoscalar-Higgs Portal Dark Matter,” JHEP 05 (2019) 096, arXiv:1812.04092 [hep-ph].
  * [28] S. Bhattacharya, P. Ghosh, and N. Sahu, “Multipartite Dark Matter with Scalars, Fermions and signatures at LHC,” JHEP 02 (2019) 059, arXiv:1809.07474 [hep-ph].
  * [29] M. Duch, B. Grzadkowski, and D. Huang, “Strong Dark Matter Self-Interaction from a Stable Scalar Mediator,” JHEP 03 (2020) 096, arXiv:1910.01238 [hep-ph].
  * [30] M. Cirelli, N. Fornengo, and A. Strumia, “Minimal dark matter,” Nucl. Phys. B 753 (2006) 178–194, arXiv:hep-ph/0512090.
  * [31] D. G. E. Walker, “Dark Matter Stabilization Symmetries from Spontaneous Symmetry Breaking,” arXiv:0907.3146 [hep-ph].
  * [32] N. Bernal, D. Restrepo, C. Yaguna, and O. Zapata, “Two-component dark matter and a massless neutrino in a new $B-L$ model,” Phys. Rev. D 99 no. 1, (2019) 015038, arXiv:1808.03352 [hep-ph].
  * [33] B. D. Sáez, F. Rojas-Abatte, and A. R. Zerwekh, “Dark Matter from a Vector Field in the Fundamental Representation of $SU(2)_{L}$,” Phys. Rev. D 99 no. 7, (2019) 075026, arXiv:1810.06375 [hep-ph].
  * [34] A. Belyaev, G. Cacciapaglia, J. Mckay, D. Marin, and A. R. Zerwekh, “Minimal Spin-one Isotriplet Dark Matter,” Phys. Rev. D 99 no. 11, (2019) 115003, arXiv:1808.10464 [hep-ph].
  * [35] O. Cata and A. Ibarra, “Dark Matter Stability without New Symmetries,” Phys. Rev. D 90 no. 6, (2014) 063509, arXiv:1404.0432 [hep-ph].
  * [36] G. Belanger and J.-C. Park, “Assisted freeze-out,” JCAP 03 (2012) 038, arXiv:1112.4491 [hep-ph].
  * [37] F. D’Eramo and J. Thaler, “Semi-annihilation of Dark Matter,” JHEP 06 (2010) 109, arXiv:1003.5912 [hep-ph].
  * [38] G. Belanger, K. Kannike, A. Pukhov, and M. Raidal, “Impact of semi-annihilations on dark matter phenomenology - an example of $Z_{N}$ symmetric scalar dark matter,” JCAP 04 (2012) 010, arXiv:1202.2962 [hep-ph].
  * [39] Y. Nomura and J. Thaler, “Dark Matter through the Axion Portal,” Phys. Rev. D 79 (2009) 075008, arXiv:0810.5397 [hep-ph].
  * [40] A. Ibarra, A. S. Lamperstorfer, S. López-Gehler, M. Pato, and G. Bertone, “On the sensitivity of CTA to gamma-ray boxes from multi-TeV dark matter,” JCAP 09 (2015) 048, arXiv:1503.06797 [astro-ph.HE]. [Erratum: JCAP 06, E02 (2016)].
  * [41] A. Merle, A. Schneider, and M. Totzauer, “Dodelson-Widrow Production of Sterile Neutrino Dark Matter with Non-Trivial Initial Abundance,” JCAP 04 (2016) 003, arXiv:1512.05369 [hep-ph].
  * [42] M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, “Observational Constraints on Decoupled Hidden Sectors,” Phys. Rev. D 94 no. 6, (2016) 063506, arXiv:1604.02401 [astro-ph.CO]. [Erratum: Phys.Rev.D 96, 109902 (2017)].
  * [43] Y. Mambrini, S. Profumo, and F. S. Queiroz, “Dark Matter and Global Symmetries,” Phys. Lett. B 760 (2016) 807–815, arXiv:1508.06635 [hep-ph].
  * [44] H. M. Lee, M. Park, and V. Sanz, “Interplay between Fermi gamma-ray lines and collider searches,” JHEP 03 (2013) 052, arXiv:1212.5647 [hep-ph].
  * [45] R. N. Lerner and J. McDonald, “Gauge singlet scalar as inflaton and thermal relic dark matter,” Phys. Rev. D 80 (2009) 123507, arXiv:0909.0520 [hep-ph].
  * [46] 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.
  * [47] K. Kainulainen, S. Nurmi, T. Tenkanen, K. Tuominen, and V. Vaskonen, “Isocurvature Constraints on Portal Couplings,” JCAP 06 (2016) 022, arXiv:1601.07733 [astro-ph.CO].
  * [48] M. Blennow, E. Fernandez-Martinez, and B. Zaldivar, “Freeze-in through portals,” JCAP 01 (2014) 003, arXiv:1309.7348 [hep-ph].
  * [49] G. Bélanger et al., “LHC-friendly minimal freeze-in models,” JHEP 02 (2019) 186, arXiv:1811.05478 [hep-ph].
  * [50] T. Bringmann, P. F. Depta, M. Hufnagel, J. T. Ruderman, and K. Schmidt-Hoberg, “Pandemic Dark Matter,” arXiv:2103.16572 [hep-ph].
  * [51] A. Hryczuk and M. Laletin, “Dark matter freeze-in from semi-production,” arXiv:2104.05684 [hep-ph].
  * [52] G. Bélanger, F. Boudjema, A. Pukhov, and A. Semenov, “micrOMEGAs4.1: two dark matter candidates,” Comput. Phys. Commun. 192 (2015) 322–329, arXiv:1407.6129 [hep-ph].
  * [53] L. Husdal, “On Effective Degrees of Freedom in the Early Universe,” Galaxies 4 no. 4, (2016) 78, arXiv:1609.04979 [astro-ph.CO].
  * [54] G. Bélanger, A. Pukhov, C. E. Yaguna, and O. Zapata, “The Z5 model of two-component dark matter,” JHEP 09 (2020) 030, arXiv:2006.14922 [hep-ph].
  * [55] 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].
  * [56] B. D. Sáez, K. Möhling, and D. Stöckinger, “Two Real Scalar WIMP Model in the Assisted Freeze-Out Scenario,” arXiv:2103.17064 [hep-ph].
  * [57] 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].
  * [58] Planck collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” (2018) , arXiv:1807.06209 [astro-ph.CO].
  * [59] 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.
  * [60] XENON Collaboration, E. Aprile et al., “Physics reach of the XENON1T dark matter experiment,” JCAP 04 (2016) 027, arXiv:1512.07501 [physics.ins-det].
  * [61] T. Binder, T. Bringmann, M. Gustafsson, and A. Hryczuk, “Early kinetic decoupling of dark matter: when the standard way of calculating the thermal relic density fails,” Phys. Rev. D 96 no. 11, (2017) 115010, arXiv:1706.07433 [astro-ph.CO]. [Erratum: Phys.Rev.D 101, 099901 (2020)].
  * [62] S. S. AbdusSalam et al., “Simple and statistically sound strategies for analysing physical theories,” arXiv:2012.09874 [hep-ph].
  * [63] A. Ibarra, S. Lopez Gehler, and M. Pato, “Dark matter constraints from box-shaped gamma-ray features,” JCAP 07 (2012) 043, arXiv:1205.0007 [hep-ph].
  * [64] G. Vertongen and C. Weniger, “Hunting Dark Matter Gamma-Ray Lines with the Fermi LAT,” JCAP 05 (2011) 027, arXiv:1101.2610 [hep-ph].
  * [65] A. Ghosh, A. Ibarra, T. Mondal, and B. Mukhopadhyaya, “Gamma-ray signals from multicomponent scalar dark matter decays,” JCAP 01 (2020) 011, arXiv:1909.13292 [hep-ph].
  * [66] F. S. Queiroz and C. Siqueira, “Search for Semi-Annihilating Dark Matter with Fermi-LAT, H.E.S.S., Planck, and the Cherenkov Telescope Array,” JCAP 04 (2019) 048, arXiv:1901.10494 [hep-ph].
  * [67] A. Reinert and M. W. Winkler, “A Precision Search for WIMPs with Charged Cosmic Rays,” JCAP 01 (2018) 055, arXiv:1712.00002 [astro-ph.HE].