# Gamma-ray Signal from $Z_{N\geq 3}$ Dark Matter-Companion Models

Jun Guo<EMAIL_ADDRESS>College of Physics and Communication Electronics,
Jiangxi Normal University, Nanchang 330022, China  Zhaofeng Kang
<EMAIL_ADDRESS>School of physics, Huazhong University of Science and
Technology, Wuhan 430074, China  Ji-Gang Zhao<EMAIL_ADDRESS>College of
Physics and Communication Electronics, Jiangxi Normal University, Nanchang
330022, China

###### Abstract

In Ref. Guo:2021rre , we proposed to replace the final dark matter (DM)
particle in the semi-annihilation mode $\rm DM+DM\to antiDM+Higgs~{}boson$
with its $Z_{N\geq 3}$ companion, thus reducing DM number density without DM-
nucleon scattering. In this work, we study the indirect detection signals from
DM annihilation, the Higgs boson pair with one of them from the companion
decay being on- or off- shell, depending on the DM-companion mass splitting.
We generate the photon spectrum by using PYTHIA8 and study the properties of
the spectrum, to find that the hard part of the spectrum in our model is
mainly shaped by the direct Higgs boson and thus does not differ much from
that of the conventional semi-annihilation mode. Using the Fermi-LAT data of
white dwarfs, we derive the current limit of the DM annihilation cross section
for ${\rm DM+DM\to companion^{*}+Higgs~{}boson}$, and for the relatively light
DM, it reaches the typical thermal cross section. However, for the TeV scale
DM, we have to rely on the Cherenkov Telescope Array, which is able to rule
out the whole parameter space except for the coannihilation region.

##  I Introduction

The nature of the putative particle dark matter (DM) remains an puzzle in the
field of particle physics, and the weakly interacting massive particles (WIMP)
have been the most promising candidates due to the WIMP miracle to understand
the DM relic density $\Omega h^{2}\simeq 0.12$ today. In general, the WIMP DM
is supposed to scatter with nucleons. The accordingly devised DM direct
detection experiments such as XENON1T XENON:2017lvq  and PandaX PandaX:2014mem
have pushed the detection target to the multi-ton-scale level to hunt for the
rare scattering signals between DM and the target atoms, however, no
significant signals have been observed yet. Therefore, they merely impose
stringent bounds on the parameter space for the typical WIMP DM, whose
annihilation cross section is closely tied with the DM-nucleon scattering
rate, rendering the WIMP miracle questionable and even the WIMP crisis.

It might be a crisis for the benchmark WIMP DM candidate such as neutralino in
supersymmetry and the singlet scalar DM through the Higgs portal. But there
are still many WIMPs living a carefree life that can break the general
relations between DM-DM annihilation and DM-nucleon scattering. An example is
the $Z_{N\geq 3}$ symmetric DM-companion model we proposed before Guo:2021rre
, where DM has a companion; moreover, the quantum number assignment admits the
semi-annihilation mode: ${\rm DM+DM\to companion^{*}+Higgs~{}boson}$, where
the anti-companion will decay back to DM plus some other SM particles, for
instance an on or off shell Higgs boson depending on the mass splitting
between DM and the companion, $\Delta m$. Remarkably, such kind of semi-
annihilation does not give rise to scattering between WIMP DM and nucleon and
thus is able to maintain the WIMP miracle and easily evade the direct
detection bounds 111Due to the characterized DM and DM instead of DM anti-DM
annihilation mode, such models provide a good opportunity to seed matter
asymmetry Chen:2024arl .. Then, we need to search for this type of WIMP DM by
indirectly detect its Higgs boson signals from DM annihilation today.

The best way to look for these products in the sky is to collect their gamma
rays injected into the cosmic ray. In this work, based on the simplest model
with two identical (except for mass) complex scalar fields with respect to
$Z_{3}$, we simulated the resulting gamma spectra for different DM-companion
mass splitting $\Delta m$, to find that for most interesting cases, the
contribution from the direct Higgs boson plays the dominant role while the
secondary Higgs boson is relatively soft. As a consequence, our signals, in
particular for a sufficiently small $\Delta m$, do not differ much from the
conventional $Z_{3}$ symmetric DM model DEramo:2010keq ; Belanger:2012vp ;
Ko:2014nha ; Guo:2015lxa ; Guo:2023kqt ; DiazSaez:2022nhp ; Liu:2023kil ,
which produces the single Higgs boson signal ${\rm DM+DM\to
DM^{*}+Higgs~{}boson}$. Therefor, our study also apply to these models; see
some variants with more complicated model extensions Queiroz:2019acr ;
Toma:2021vlw .

Then, we derive the indirect detection bound on this model using the Fermi-LAT
observation, which covers the energy range between 20 MeV and 300 GeV and has
been observing the gamma-ray sky for 15 years, drawing the most detailed map
to date. Many works have made a deep analysis of the data, focusing on a wide
variety of targets. For instance, Ref. Li:2013qya ; Karwin:2020tjw  tried to
search for the DM annihilation signals in the direction of the Andromeda
galaxy (M31), while Ref. Fermi-LAT:2015sau ; Karwin:2016tsw ; Fermi-
LAT:2017opo ; DiMauro:2021raz  studied the data from the Galactic Center and
founded hints for gamma-ray excess, but it is not a significant DM signal. In
particular, Ref. Fermi-LAT:2010cni ; Fermi-LAT:2011vow ; Geringer-
Sameth:2014qqa ; Fermi-LAT:2015att ; Fermi-LAT:2015ycq ; Fermi-LAT:2016uux ;
Calore:2018sdx ; Hoof:2018hyn ; DiMauro:2021qcf ; McDaniel:2023bju  analyzed
the Milky Way dwarf spheroidal satellite galaxies (dSphs) using the LAT data,
and placed the most stringent constraint on signals from DM activities. The
dSphs of the MilkyWay are the most DM dominated celestial bodies and lack of
non-thermal astrophysical processes, so they are the most prominent targets
for searching for DM decay or annihilation. In this work, we will use the
dSphs of the MilkyWay to obtain the upper bound on the cross section for ${\rm
DM+DM\to companion^{*}+Higgs~{}boson}$. We find that the current bound reaches
the typical thermal WIMP DM cross section for DM below the TeV scale. Whereas
for the even heavier DM above the TeV scale, sensitivity improvement of order
of magnitude is required, and the Cherenkov Telescope Array (CTA)
CTAConsortium:2017dvg  offers an opportunity to cover this parameter space.

The paper is organized as follows, in Section II, we will briefly introduce
our $Z_{3}$ symmetric DM-companion model. In Section III, we will study the
spectral properties of dark matter annihilation in our model. In Section IV.2,
we will first present the data that we will use, and then derive the upper
limit of the DM annihilation cross section. Finally, we will give our
conclusions in Section V.

##  II $Z_{3}$ symmetric DM models: from Higgs to companion portal

The $Z_{3}$ symmetric model is characterized by the presence of semi-
annihilation mode. The simplest model involves a complex scalar $S_{1}$
transforms under $Z_{3}$ as $S_{1}\rightarrow e^{ik2\pi/3}S_{1}$ with $k=1$ or
$2$, and the resulting Lagrangian is given by

$\displaystyle-{\cal L}_{Z_{3}}\supset$ $\displaystyle
m_{1}^{2}S_{1}S_{1}^{*}+\lambda_{1h}|S_{1}|^{2}|\Phi|^{2}+\left(\frac{A_{1}S_{1}^{3}}{3}+c.c\right),$
(1)

the real scalar part component of $\Phi$ is the Higgs boson $H$, and the cubic
term is new compared to the usual $Z_{2}$ symmetric Higgs portal dark matter
model and gives rise to the characterized semi-annihilation mode $SS\to
S^{*}H$, which, through a properly large $A_{1}$, provides a way to attenuate
the stringent exclusion to the usual $Z_{2}$ Higgs portal dark matter from
direct matter detection. But this attenuation is not remarkable Athron:2018ipf
since the effective dark matter annihilation rate still requires a sizable
Higgs portal coupling $\lambda_{1h}$.

A complete attenuation is available via the mechanism proposed by us in Ref.
Guo:2021rre . It considers at least one more dark species charged under
$Z_{N}$ (dubbed as companion), and then dark matter can gain correct relic
density via the semi-annihilation into the companion; the Higgs portal that
generates DM-nucleon interaction can be removed by hand or by symmetry. The
mechanism can be realized simply by extending the model specified by Eq. (1)
with one more scalar $S_{2}$, a copy of $S_{1}$ in quantum number 222 The
model has an updated version if dark matter is instead a Dirac fermion $\Psi$,
which automatically does not have Higgs portal interaction at all. However,
such models will have very different DM indirect detection signals than the
photon signals discussed in the current work, and we leave more details in the
discussion.. The resulting most general Lagrangian is tedious, and it is
unnecessary to list them all. So, we just present the relevant terms in the
mass basis of $S_{1,2}$,

$\displaystyle-{\cal L}_{Z_{3}}\supset
m_{S_{1}}^{2}S_{1}S_{1}^{*}+m_{S_{2}}^{2}S_{2}S_{2}^{*}+\left(\frac{A_{1}}{3}S_{1}^{3}+\mu_{12}S_{1}S_{2}^{*}H+c.c.\right).$
(2)

Other terms, in particular $|S_{1}|^{2}|\Phi|^{2}$, are assumed to be small
and irrelevant. The key feature is that the dark matter candidate $S_{1}$ does
not have the conventional Higgs portal term as in Eq. (1). However, as long as
the masses lying in the window

$\displaystyle m_{S_{2}}>m_{S_{1}}>(m_{S_{2}}+m_{H})/2>m_{S_{2}}/2,$ (3)

terms in the bracket will give rise to two-body semi-annihilation mode for a
pair of dark matter particles $S_{1}$:

$\displaystyle S_{1}S_{1}\to S^{*}_{2}(\to S^{*}_{1}+H~{}{\rm
or}~{}H^{off})+H,$ (4)

with $H^{off}$ denoting for the off-shell Higgs boson, which is the situation
if the masses of dark matter and companion are sufficiently close to each
other, $\Delta m\equiv m_{S_{2}}-m_{S_{1}}<m_{H}=125.6$ GeV; it is a typical
case in the relatively light dark matter region.

In particular, in the highly degenerate region ($\Delta m/m_{S_{1}}\lesssim
20\%$), the model actually resembles the previous model in the sense of
indirect detection. This case corresponds to the coannihilation region
analyzed in Ref. Guo:2021rre , which may lead to the reduction of
$\langle\sigma v\rangle_{S_{1}S_{1}\to S_{2}^{*}H}$ given a large
$\langle\sigma v\rangle_{S_{2}S_{2}^{*}\to HH}$.

##  III DM annihilation and photon spectrum

In this section we analyze the photon spectra from the Higgs bosons produced
in different annihilation modes: ${\rm DM+DM\to DM^{*}+Higgs~{}boson}$ (from
the conventional $Z_{3}$ Higgs portal model), ${\rm DM+DM\to
companion^{*}+Higgs~{}boson}$ (from the $Z_{3}$ DM-companion model) and ${\rm
DM+DM\to Higgs~{}boson~{}pair}$ (from the conventional $Z_{2}$ Higgs portal
model). The photon spectra of the first and third modes are determined solely
by DM mass and their features are clear, while the second case involves one
more parameter and is more complicated. Hence, our subsequent analysis will
focus on the second mode.

For the model we studied, the DM annihilation mode is $S_{1}S_{1}\rightarrow
HS_{2}^{*}$, and the followed decay of $S^{*}_{2}$ depends on the mass
splitting $\Delta m$. When $\Delta m>m_{H}$, $S_{2}^{*}$ decays into DM plus
an on-shell Higgs boson (the secondary Higgs boson), while in the case of
$\Delta m\leq m_{H}$, the on-shell Higgs boson becomes off shell, denoted as
$H^{off}$. For both cases, DM annihilation today yields a characteristic
photon spectrum. To get these spectra, we take the following procedure:

  1. 1.

First, we generate the model file using FeynRules Alloul:2013bka .

  2. 2.

Then, we generate the process $S_{1}S_{1}\rightarrow HS_{2}^{*}$ using
MadGraph Alwall:2014hca  without hardronization and detector simulation.

  3. 3.

Next, we pass the events to PYTHIA8 Sjostrand:2014zea  to implement
hardronization, and finally obtain the DM annihilation photon spectrum.

In the rest part of this section, we will discuss the two cases separately,
and we will also compare them with those from other two cases.

###  III.1 Photon spectrum with an off-shell Higgs boson from companion decay

Both DM mass $m_{S_{1}}$ and DM-companion mass splitting $\Delta m$ can affect
the shape of the photon spectrum from DM annihilation. When $m_{S_{1}}$
becomes larger, both the annihilation products, the DM companion $S_{2}^{*}$
and the Higgs boson $H$, would be more boosted and then give rise to more high
energy photons. This is clearly seen in the left panel of Fig. 1, the photon
spectrum for 300 GeV DM is much softer than that for the 800 GeV DM.

The photon spectrum is affected by the DM-companion mass splitting from two
aspects, the primary (accessible) decay channels of $H^{off}$ and as well the
boost of its decay products. In this work, we consider two levels of splitting
and the corresponding dominant decay modes of $S_{2}^{*}$: 1) A relatively
large splitting with $m_{H}>\Delta m\geq 2m_{b}$, and then the dominant decay
channel of $S^{*}_{2}$ is $S^{*}_{2}\rightarrow S^{*}_{1}b\bar{b}$; 2) A
relatively small splitting with $2m_{\tau}<\Delta m<2m_{b}$, and then
$S^{*}_{1}\tau\bar{\tau}$ dominates the decay products of $S^{*}_{2}$, as show
on the Fig. 1 left panel, we compare the cases of $\Delta m=5$ GeV and $\Delta
m=10$ GeV for the DM mass of 800 GeV, and find that the difference between
them is tiny. Even smaller splitting is not considered since it means a
significant fine-tuning of the DM-companion mass spectrum. The impact on
boosting the fermion pair can be seen from the energy of $S_{2}^{*}$ (in the
center mass frame) produced in dark matter annihilation:

$\displaystyle E_{S^{*}_{2}}\simeq m_{S_{1}}+\frac{(m_{S_{1}}+\Delta
m)^{2}-m_{H}^{2}}{4m_{S_{1}}}.$ (5)

Therefore, for a fixed DM mass, the resulting fermion pair becomes more
energetic as $\Delta m$ increases. Its impact on the final photon spectrum is
complicated, since the fermion pair comes from $S^{*}_{2}$ three-body decay,
and thus in the following we give numerical samples, shown on the left panel
of Fig. 1.

Our samples are for a given DM mass of 800 GeV, with two cases of mass
splitting: the smaller $\Delta m=10$ GeV and the larger $\Delta m=100$ GeV. It
is found that the latter generates a slightly softer photon spectrum, which is
mainly ascribed to the soften of the direct Higgs boson $H$. To demonstrate
this, we show the components of the photon spectrum in the right panel of Fig.
1. From it one can see that, in each case, the component from the direct Higgs
boson (dashed line) always exceeds that from $S_{2}^{*}$ or the off-shell
Higgs boson (dotted line), over the whole spectrum; in particular in the
energy region above around 10 GeV, the latter is almost negligible. Moreover,
compared to $\Delta m=10$ GeV, in the case of $\Delta m=100$ GeV, the $H$
produces fewer photons starting from $E_{\gamma}\gtrsim 5$ GeV. Although in
the latter case, $H^{off}$ generates more photons starting from
$E_{\gamma}\gtrsim 1$ GeV, this secondary component only gives a mild
modification to the spectrum generated by $H$ 333One can understand this along
this line: a larger $\Delta m$ means a larger mass of DM companion
$S_{2}^{*}$, which then takes away more kinetic energy according to Eq. (6) ,
but at the same time, the direct Higgs boson, the main source of photons, is
less boosted, so the harder photons from its cascade decay is reduced. .

Therefore, in this sense, the indirect detection signals from the ${\rm
DM}+HH^{off}$ mode resembles that from the ${\rm DM}+H$ mode contaning a
single direct Higgs boson. This is in accordance with the expectation that the
$Z_{3}$ symmetric DM-companion model reduces to the conventional $Z_{3}$ Higgs
portal model when $\Delta m\ll m_{S_{1}}$; explicitly, for the given DM mass
800 GeV, we compare the photon spectra from the former with that from the
latter (dotted pink) in the left panel of Fig. 1. We find that, for a
sufficiently small mass splitting, they coincide in the hard part and only
show a mild difference in the soft part below the GeV scale. Hence, the
constraint on $\langle\sigma v\rangle_{S_{1}S_{1}\rightarrow HS_{2}^{*}}$ at a
given dark matter mass made in this limit can be well applied to
$\langle\sigma v\rangle_{S_{1}S_{1}\rightarrow HS_{1}^{*}}$ from the
conventional $Z_{3}$ Higgs portal model.

The photon spectrum from direct Higgs boson pair mode is also shown in the
left panel of Fig. 1, and one can see that it is obviously more energetic than
other cases: roughly speaking, the hardness of its spectrum is roughly twice
than others in the relatively high energy region. The reason is not difficult
to understand. For this DM annihilation mode, all the energy of the
annihilating DM pair is transferred to the Higgs boson pair, whereas in the
semi-annihilation modes, both in the conventional $Z_{3}$ Higgs portal model
or its DM-companion variant, the DM particle contained in the final state
always takes away a substantial portion of the initial DM energy. This
argument means that the photon spectrum corresponding to the combination of a
direct Higgs boson and secondary on-shell Higgs boson discussed in the
following will not give a much significant difference.

Figure 1: Left: Photon spectra for different dark matter mass $m_{S_{1}}$ and
DM-companion mass splittings $\Delta m$ for the $HS_{2}^{*}(\to H^{off}+\rm
DM)$ mode. To make a comparison, the photon spectra from the $HS$ mode (dotted
pink) and $HH$ mode (dotted black) are shown. Right: two samples of the
components of the photon spectrum for the $HS_{2}^{*}(\to H^{off}+\rm DM)$
mode with $\Delta m=10$ GeV (green) and $\Delta m=100$ GeV (blue), taking
$m_{S_{1}}=800$ GeV, where the total spectrum (solid) is a sum of the spectra
from the direct Higgs boson (dashed) and secondary Higgs boson (dotted).

###  III.2 Photon spectrum with an on-shell Higg boson from companion decay

In the case of $\Delta m>m_{H}$, in addition to the direct Higgs boson, DM
direct annihilation produces another on-shell Higgs boson, the secondary Higgs
boson from $S_{2}^{*}$ two-body decay. Unlike the previous case, here the mass
splitting $\Delta m$ can only affect the decay spectrum in one way, the
kinetic boost of the Higgs bosons.

The shape of the photon spectrum is again controlled by $m_{S_{1}}$ and
$\Delta m$. From the left panel of Fig. 2, it is ready to find that the shape
mainly depends on the DM mass. As for the mass splitting, its impact is not so
significant compared to the previous case containing $H^{off}$. To understand
the reason, in the right panel of Fig. 2, we show the photon spectra for
$\Delta m=130$ GeV and $\Delta m=500$ GeV, respectively, both having DM mass
$m_{S_{1}}=800$ GeV. We can see that, the photon spectrum again is almost
dominated by the component from the direct Higgs boson, however, for the large
mass splitting, when $E_{\gamma}\gtrsim{\cal O}(40)$ GeV, the component from
the secondary Higgs boson begins to be comparable and eventually dominated at
the high-energy tail region.

Once again, we compare the photon spectra in the ${\rm DM}+2H$ mode with other
spectra in Fig. 2, taking a 800 GeV dark matter. It is seen that, as
previously commented, they are indeed similar to the shape of the DM$+H$ mode
in the relatively higher energy range, and it is even harder for
$E_{\gamma}\gtrsim{\cal O}(10)$ GeV, by virtue of its more boosted direct
Higgs boson. Therefore, the real reason is traced back to why the secondary
Higgs boson again plays a minor role. As mentioned before, it is due to the
reason that most of the energy of $S_{2}^{*}$ is taken away by DM, so the
secondary Higgs boson is relatively soft. Only when the mass splitting is
large enough, $\Delta m\gg m_{H}$, can it be significantly accelerated, thus
hardening its component in the photon spectrum. At the same time, the direct
Higgs boson becomes softer by increasing $\Delta m$, which is seen from the
expression of its energy,

$\displaystyle E_{H}\simeq m_{S_{1}}-\frac{(m_{S_{1}}+\Delta
m)^{2}-m_{H}^{2}}{4m_{S_{1}}}.$ (6)

It is substantially reduced when $\Delta m$ becomes comparable to $m_{S_{1}}$.

Figure 2: Left: Photon spectra for different DM mass and DM-companion mass
splittings for the $HS_{2}^{*}(\to H+\rm DM)$ mode. Right: two samples of the
components of the spectrum with $\Delta m=130$ GeV (green) and $\Delta m=500$
GeV (blue) with $m_{S_{1}}=800$ GeV, where the total spectrum (solid) is a sum
of the spectra from the direct (dashed) and secondary (dotted) Higgs bosons,
respectively.

##  IV Indirect Detection

After analyzing the characteristic photon spectrum from the Higgs boson
“pair”, injected by DM semi-annihilation into the cosmic ray, the main subject
of this section is to apply the FERMI-LAT data to make constraint on such kind
of DM models. But the prospect at CTA will also be briefly studied.

###  IV.1 Injection, propogation and observation

The indirect detection experiments such as FERMI-LAT observe the gamma-ray
flux signal in a particular region of the sky. The expected signal from the DM
hypothesis of the observed sky map region in the energy range $[E_{\rm
min},E_{\rm max}]$ is factorizied as:

$\displaystyle\phi(\Delta\Omega;E_{\rm min},E_{\rm
max})=\frac{1}{4\pi}\frac{\langle\sigma v\rangle}{2m_{\rm DM}^{2}}\int_{E_{\rm
min}}^{E_{\rm max}}\frac{dN_{\gamma}}{dE_{\gamma}}\times\rm J,$ (7)

where $\langle\sigma v\rangle$ is the thermally averaged annihilation cross
section of DM times the relative velocity, and the
$\frac{dN_{\gamma}}{dE_{\gamma}}$ is the DM annihilation photon spectrum of a
particular DM model. The $\rm J$-factor is the integration of DM density
square along the line of sight $s$ and solid angle $\Delta\Omega$, given as:

$\displaystyle{\rm J}=\int_{\Delta\Omega}\int_{l.o.s}d\Omega
ds\rho^{2}(r(s,l,b)),$ (8)

where $r$ is the radial distance from the Galactic Center,

$\displaystyle r=\sqrt{s^{2}+R^{2}-2sR\cos{l}\cos{b}},$ (9)

where $R$ is the distance from the Sun to the Galactic Center, s is the
distance to the Earth and $(l,b)$ is the Galactic coordinates. $\rho(r)$ is
the DM density profile, in this study, we are setting the limitation bound by
using the dSphs data, and it is found that the $\rm J$-factor is fairly
insensitive to it Martinez:2009jh ; Strigari:2012acq . Therefore, we ignore
the effect of DM density profile to the indirect detection signal, and the
J-factor we use could be found in Fermi-LAT:2015att ; McDaniel:2023bju .

###  IV.2 Fermi-LAT dSphs Data Analysis and Limitations

We will adopt the available Fermi-LAT dSphs gamma-ray data for the analysis in
this work. In Fermi-LAT:2015att , the authors analyzed the 6 years of Fermi-
LAT dSphs data processed with $\rm pass$ 8\. They analyzed 15 dSphs with
measured J-factor, and provided the spectrum energy distribution (SED) for
each dSph. For each SED, the energy is divided into 24 energy bins, and a
likelihood table of the photon energy flux values is provided in each bin.
Recently, authors of Ref. McDaniel:2023bju  have generously provided the
latest 14.3 years of LAT data. They analyzed 50 dSphs, and divided them into
three sets: 1) a measured set which contains 30 dSphs with measured J-factors;
2) the benchmark set which contains 42 dSphs, some of which only have
estimated J-factors; 3) the inclusive set which contains all the 50 dSphs. In
their analysis, all the dSphs gamma-ray data is selected from the $\rm
P8R3\\_SOURCE\\_V3$ class, with the SEDs ranging from 500 MeV to 1 TeV, also
divided into 24 bins.

We calculate the upper limit of Fermi-LAT dSphs data to our model by
performing a joint likelihood analysis, which is $\Pi_{i}\mathcal{L}_{i}$, the
product of the likelihood of dSphs of interest. For each such dSph, the
likelihood function for DM hypothesis is a product of the combined likelihood
of dSph photon flux and the likelihood of dSph J-factor,

$\displaystyle\mathcal{L}_{i}(m_{\rm DM},\langle\sigma
v\rangle|\mathcal{D}_{i})=\mathcal{L}_{i}(m_{\rm DM},\langle\sigma
v\rangle|\mathcal{D}_{i})\mathcal{L}_{J}(J_{i}|J_{obs,i},\sigma_{i}).$ (10)

The former is the production over the likelihood of each energy bin $E_{n}$,
$\prod_{n}\mathcal{L}_{E_{n}}(\frac{d\phi}{dE},E_{n})$, where
$\frac{d\phi}{dE}$ is the DM photon-flux, a function of $m_{\rm DM}$,
$\langle\sigma v\rangle$ and the J-factor. The likelihood function of dSph
J-factor takes the form of

$\displaystyle\mathcal{L}_{J}(J_{i}|J_{obs,i},\sigma_{i})=\frac{1}{\ln{10}J_{obs,i}\sqrt{2\pi}\sigma_{i}}\exp(-\left(\log_{10}(J_{i})-\log_{10}(J_{obs,i})\right)^{2}/2\sigma_{i}^{2}).$
(11)

In our analysis, we will use the benchmark set to analyze the 14-years data
provided in McDaniel:2023bju , with the J-factor value $J_{obs,i}$ and the
uncertainty $\sigma_{i}$ shown in Table. 2. For the 6-years data analysis, we
use the 15 dSphs provided in Fermi-LAT:2015att , with the J-factor value
$J_{obs,i}$ and the uncertainty $\sigma_{i}$ given in Table. 1.

To determine the best-fit model under the maximum likelihood estimation, we
use the test statistic (TS) which is defined as Cowan:2010js

$\displaystyle\rm{TS}=-2\ln{\frac{\mathcal{L}_{max}}{\mathcal{L}(m_{DM},\langle\sigma
v\rangle|\mathcal{D})}}$ (12)

where $\mathcal{L}_{\rm max}$ is the maximum likelihood.

Table 1:  Milky way dSphs’ properties for Fermi-LAT:2015att  Name | $\ell$ | $b$ | Dist. | $\log_{10}({J_{\textrm{obs}}})$
---|---|---|---|---
| (deg) | (deg) | (kpc) | ($\log_{10}(\text{GeV}^{2}\text{cm}^{-5})$)
Bootes I | 358.1 | 69.6 | 66 | $18.8\pm 0.22$
Canes Ven. II | 113.6 | 82.7 | 160 | $17.9\pm 0.25$
Carina | 260.1 |  $-$22.2 | 105 | $18.1\pm 0.23$
Coma Ber. | 241.9 | 83.6 | 44 | $19.0\pm 0.25$
Draco | 86.4 | 34.7 | 76 | $18.8\pm 0.16$
Fornax | 237.1 |  $-$65.7 | 147 | $18.2\pm 0.21$
Hercules | 28.7 | 36.9 | 132 | $18.1\pm 0.25$
Leo II | 220.2 | 67.2 | 233 | $17.6\pm 0.18$
Leo IV | 265.4 | 56.5 | 154 | $17.9\pm 0.28$
Sculptor | 287.5 |  $-$83.2 | 86 | $18.6\pm 0.18$
Segue 1 | 220.5 | 50.4 | 23 | $19.5\pm 0.29$
Sextans | 243.5 | 42.3 | 86 | $18.4\pm 0.27$
Ursa Maj. II | 152.5 | 37.4 | 32 | $19.3\pm 0.28$
Ursa Minor | 105.0 | 44.8 | 76 | $18.8\pm 0.19$
Willman 1 | 158.6 | 56.8 | 38 | $19.1\pm 0.31$
Table 2:  Milky way dSphs’ properties for McDaniel:2023bju  Name | $\ell$ | $b$ | Dist. | $\log_{10}({J_{\textrm{obs}}})$ | Name | $\ell$ | $b$ | Dist. | $\log_{10}({J_{\textrm{obs}}})$
---|---|---|---|---|---|---|---|---|---
| (deg) | (deg) | (kpc) | ($\log_{10}(\text{GeV}^{2}\text{cm}^{-5})$) |  | (deg) | (deg) | (kpc) | ($\log_{10}(\text{GeV}^{2}\text{cm}^{-5})$)
Aquarius II | 338.48 | -9.33 | 108.0 | $17.80\pm 0.55$ | Boötes II | 209.51 | 12.86 | 42.0 | $18.30\pm 0.95$
Canes Venatici I | 202.01 | 33.55 | 218.0 | $17.42\pm 0.16$ | Pictor I | 70.95 | -50.29 | 114.0 | $18.00\pm 0.60$
Carina | 100.41 | -50.96 | 105.0 | $17.83\pm 0.10$ | Carina II | 114.11 | -58.0 | 36.0 | $18.25\pm 0.55$
Coma Berenices | 186.75 | 23.91 | 44.0 | $19.00\pm 0.35$ | Draco | 260.07 | 57.92 | 76.0 | $18.83\pm 0.12$
Draco II | 238.17 | 64.58 | 22.0 | $18.93\pm 1.54$ | Eridanus II | 56.09 | -43.53 | 380.0 | $16.60\pm 0.90$
Fornax | 39.96 | -34.5 | 147.0 | $18.09\pm 0.10$ | Grus I | 344.18 | -50.18 | 120.0 | $16.50\pm 0.80$
Hercules | 247.77 | 12.79 | 132.0 | $17.37\pm 0.53$ | Horologium I | 43.88 | -54.12 | 79.0 | $19.00\pm 0.81$
Hydrus I | 37.39 | -79.31 | 28.0 | $18.33\pm 0.36^{b}$ | Leo I | 152.11 | 12.31 | 254.0 | $17.64\pm 0.13$
Leo II | 168.36 | 22.15 | 233.0 | $17.76\pm 0.20$ | Leo IV | 173.24 | -0.55 | 154.0 | $16.40\pm 1.08$
Leo V | 172.79 | 2.22 | 178.0 | $17.65\pm 0.97$ | Pegasus III | 336.1 | 5.41 | 215.0 | $18.30\pm 0.93$
Pisces II | 344.63 | 5.95 | 182.0 | $17.30\pm 1.04$ | Reticulum II | 53.92 | -54.05 | 30.0 | $18.90\pm 0.38$
Sagittarius II | 298.16 | -22.07 | 69.0 | $17.35\pm 1.36^{d}$ | Segue 1 | 151.75 | 16.08 | 23.0 | $19.12\pm 0.53$
Sextans | 153.26 | -1.61 | 86.0 | $17.73\pm 0.12$ | Tucana II | 342.98 | -58.57 | 58.0 | $18.97\pm 0.54$
Tucana IV | 0.73 | -60.85 | 48.0 | $18.40\pm 0.55^{e}$ | Ursa Major I | 158.77 | 51.95 | 97.0 | $18.26\pm 0.28$
Ursa Major II | 132.87 | 63.13 | 32.0 | $19.44\pm 0.40$ | Ursa Minor | 227.24 | 67.22 | 76.0 | $18.75\pm 0.12$
Boötes IV | 233.69 | 43.73 | 209.0 | $17.25\pm 0.60$ | Carina III | 114.63 | -57.9 | 28.0 | $19.70\pm 0.60$
Centaurus I | 189.59 | -40.9 | 116.0 | $18.14\pm 0.60$ | Cetus II | 19.47 | -17.42 | 30.0 | $19.10\pm 0.60$
Cetus III | 31.33 | -4.27 | 251.0 | $17.30\pm 0.60$ | Columba I | 82.86 | -28.01 | 183.0 | $17.60\pm 0.60$
Grus II | 331.02 | -46.44 | 53.0 | $18.40\pm 0.60$ | Phoenix II | 355.0 | -54.41 | 83.0 | $18.30\pm 0.60$
Canes Venatici II | 194.29 | 34.32 | 160.0 | $17.82\pm 0.47$ | Pictor II | 101.18 | -59.9 | 46.0 | $18.83\pm 0.60$
Reticulum III | 56.36 | -60.45 | 92.0 | $18.20\pm 0.60$ | Tucana V | 354.35 | -63.27 | 55.0 | $18.90\pm 0.60$

Figure 3: Left: Excluding lines with different DM degree of Mass splittings
for the $HS_{2}^{*}(\to H^{off}+\rm DM)$ mode. Right: Excluding lines with
different DM degree of Mass splittings for the $HS_{2}^{*}(\to H+\rm DM)$
mode. The 14-year (6-year) data limitations is shown in solid (dotted) lines,
and the yellow band is the 1-$\sigma$ uncertainty of CTA expected exclusion
sensitivity of $W^{+}W^{-}$ annihilation channel.

In the Fig. 3, using both the 6-year and 14-year Fermi-LAT dSphs data, we
demonstrate the upper bounds on the $m_{S_{1}}-\langle\sigma v\rangle$ plane
at a 95% confidence level (C.L.) with $\rm TS=-2.71$. First of all, the
current exclusion is already close to the typical thermal DM cross section
with $\langle\sigma v\rangle\simeq 3\times 10^{-26}\rm cm^{-3}s^{-1}$ for the
weak scale DM. But for the relatively heavy DM of sub-TeV or even heavier, we
still need to increase sensitivity orders of magnitude. Next, it is evident
that in most of the DM mass region, the 6-year data sets an even tighter
constraint than the 14-year data; only in the relatively heavy region of a few
TeVs does the latter begin to become more stringent. This is because in the
14-year data, some dSphs, such as Reticulum II, show inconspicuous excesses at
the $2-3\sigma$ level 444It was also noticed in the previous studies Fermi-
LAT:2015ycq ; DiMauro:2021qcf ; McDaniel:2023bju . But there is no single one
that meets the detection requirement $\rm TS>25$ adopted by Fermi-LAT, and
thus can not be regarded as evidence for DM., and therefore a larger
$\langle\sigma v\rangle$ is allowed in the 95% C.L. exclusion. We have to
stress that, for simplicity, in this work we did not perform uncertainty
analysis associated with the background cosmic rays. If one takes the blank-
field analysis, which takes into account both the unknown gamma-ray sources
and the uncertainty of backgrounds, then the 14-year data could instead
provide the most stringent restrictions McDaniel:2023bju .

The upper bounds in the parameter space with an off-shell Higgs boson signal
are shown in the left panel of Fig. 3, and we fix the degree of mass splitting
$k\equiv\Delta m/m_{S_{1}}$, taking three typical $k=$0.05 (blue), $0.1$ (red)
and $0.3$ (green); the corresponding region of dark matter mass is limited by
Eq. (3), and a larger $k$ allows for a narrower interval, with a smaller upper
limit. One can see that, for a given DM mass, a larger mass splitting results
in a little more stringent bound in most DM mass region, as according to the
spectral analysis in the previous section, it produces more photons in the
$E_{\gamma}\lesssim 5$ GeV region; moreover, such photons tend to yield the
relatively stronger bound in the likelihood analysis. But overall, the off-
shell case behaves quite like the single Higgs boson mode from the semi-
annihilation model (gray), which is almost indistinguishable with the $k=0.1$
case, especially for the 6-years data.

The upper bounds on the case with two on-shell Higgs bosons are show on the
right panel of Fig. 3. Similarly, a larger mass splitting yields a slightly
stronger constraint, except for the DM mass region above $\sim 3$ TeV , where
the $k=0.05$ case instead yields a stronger constraint than the $k=0.1$ case
(but the $k=0.3$ case still gives the strongest one). Note that it occurs only
for the 14-year data, so this is a combined result between spectrum and data.
The reason is probably due to that the $S_{2}^{*}$ on-shell channel is about
to open, and the direct product Higgs is quite soft, and as discussed in the
previous section, the $S_{2}^{*}$ decay produced Higgs is even softer, so the
gamma-ray spectrum produced by these two Higgs boson is soft.

We end up this section with the comment on comparison with the ${\rm DM}+H$
and $HH$ modes. As expected, the $HS$ mode confronts with very close exclusion
lines to the off-shell mode with a small $k$. But there are still some tiny
difference. For instance, a larger $\Delta m\lesssim m_{H}$ will change the
situation, as manifest for the TeV scale DM. For the on-shell mode, the bounds
on the ${\rm DM}+H$ mode is weaker with exception in the multi-TeV region. The
$HH$ mode is always subject to stricter constraints, and its upper limit of
cross section is approximately half of the others.

###  IV.3  The heavy dark matter prospect at CTA

For the TeV-scale heavy dark matter, the limitation of Fermi-LAT to these
Higgs bosons in the sky remains relatively weak, because it just observes the
gamma-ray energy between 20 MeV and 300 GeV.

In recent years, new indirect detection experiments aiming at probing heavy
dark matter are on the table. For instance, the Cherenkov Telescope Array
(CTA) CTAConsortium:2017dvg , which has high sensitivity to the photon in the
energy range from a few tens of GeV to above 300 TeV, is currently in trial
operation. So, it is able to give a strict bound on very heavy DM, and for the
$W^{+}W^{-}$ channel, the expected exclusion line (with 1$\sigma$ uncertainty)
by analysis the Galactic Center data can reach the typical thermal cross
section, $\langle\sigma v\rangle\simeq 1\times 10^{-26}\rm cm^{3}s^{-1}$. The
specific sensitivity curvtures for the modes we studied need future study, and
in this work we take that of the $W^{+}W^{-}$ channel as an approximation
CTA:2020qlo , and demonstrate it in the plots. It is seen that the complete
parameter space of the models studied in this work can be covered, except for
the region with substantial coannihilation effect in the DM-companion model.

##  V Conclusion

The $Z_{N}$ symmetric DM-companion model offers WIMP DM candidates that may be
free of DM direct detection signals, but they still leave indirect detection
signals. In this work we study the $Z_{3}$ symmetric DM-companion model, where
the scalar dark matter annihilates into a scalar anti-DM companion and the SM
Higgs boson, and the subsequent decay production of the companion could be a
DM plus an on-shell (off-shell) Higgs boson. Similar to the conventional semi-
annihilation model, our case could evades the stringent constrains from DM
direct detection experiments, so we need to search for the photon signal
produced by the DM annihilation according to the indirect DM searches.

We study the spectrum produced by our model, and analysis the relationship
between the mass splitting $\Delta m$ and the shape of the spectrum, we also
make a comparison between our case, Higgs-portal mode and conventional semi-
annihilation mode. According to the decay of the DM companion $S_{2}$, our
study is divided into two cases according to the decay of anti-DM companion,
the off-shell case and the on-shell case, which is determined by the mass
splitting $\Delta m$. We find the photon spectrum in the off-shell case will
affected by both mass splitting and DM mass, and it resembles that from the
conventional semi-annihilation case. In the on-shell case, however, the
spectrum is dominantly affected by the DM mass, and the second Higgs boson
from the anti-DM companion decay plays only a subsidiary role in the photon
generation.

Finally, we set the upper bound on the $m_{S_{1}}-\langle\sigma v\rangle$
plane at a $95\%$ C.L level by using both 6-year 14-year Fermi-LAT dSphs data.
We find that the 6 year constraint is even tighter than the 14 year one, due
to the fact that some inconspicuous excesses at the $2-3\sigma$ level are
present in the latter. The $95\%$ CL upper limitation of the $\langle\sigma
v\rangle$ in our model reaches about $\simeq 10^{-25}$ $\rm cm^{3}/s$ at
$m_{S_{1}}\simeq 300$ GeV, which is close to the typical thermal DM cross
section, but for the heavy mass region of a few TeV, Fermi-LAT is still not
sensitive enough, we expect that the future CTA could cover this region and we
will have a detailed study in the future.

We end up the work with a comment on the case that the WIMP DM candidate of
the $Z_{N}$ symmetric DM-companion model is instead a fermion. Then, the
resulting indirect detection signals depend on model. A good case in point is
provided in Ref. Chen:2024arl , where the companion field further couples to
the SM leptonic sector with the help of a doublet $Z_{3}$ companion, carrying
identical SM quantum numbers with the SM Higgs doublet, and then these
couplings open decay channel of companion into DM plus neutrino. We can search
for such neutrino signal, other than the gamma-ray signal from the direct
Higgs boson discussed in this work.

## acknowledgments

Thanks for the useful discussion with Qiang Yuan. This work was supported by
the National Natural Science Foundation of China under Grant No. 12305111 and
11775086.

## References

  * (1) E. Aprile et al. [XENON], Eur. Phys. J. C 77 (2017) no.12, 881 doi:10.1140/epjc/s10052-017-5326-3 [arXiv:1708.07051 [astro-ph.IM]].
  * (2) X. Cao et al. [PandaX], Sci. China Phys. Mech. Astron. 57 (2014), 1476-1494 doi:10.1007/s11433-014-5521-2 [arXiv:1405.2882 [physics.ins-det]].
  * (3) F. D’Eramo and J. Thaler, JHEP 06, 109 (2010).
  * (4) G. Belanger, K. Kannike, A. Pukhov and M. Raidal, JCAP 04, 010 (2012).
  * (5) P. Ko and Y. Tang, JCAP 05, 047 (2014).
  * (6) J. Guo, Z. Kang, P. Ko and Y. Orikasa, Phys. Rev. D 91, no.11, 115017 (2015).
  * (7) J. Guo, L. Wu and B. Zhu, Phys. Lett. B 840 (2023), 137853 doi:10.1016/j.physletb.2023.137853 [arXiv:2302.06159 [hep-ph]].
  * (8) B. Díaz Sáez and K. Ghorbani, JCAP 02, 002 (2023).
  * (9) A. Liu, Z. L. Han, Y. Jin and H. Li, Phys. Rev. D 108, no.7, 075021 (2023).
  * (10) J. Guo, Z. Kang and P. Zhang, Phys. Lett. B 826, 136885 (2022).
  * (11) F. S. Queiroz and C. Siqueira, JCAP 04 (2019), 048 doi:10.1088/1475-7516/2019/04/048 [arXiv:1901.10494 [hep-ph]].
  * (12) T. Toma, Phys. Rev. D 105 (2022) no.4, 043007 doi:10.1103/PhysRevD.105.043007 [arXiv:2109.05911 [hep-ph]].
  * (13) B. Abi et al. [DUNE], [arXiv:2002.03005 [hep-ex]].
  * (14) E. Richard et al. [Super-Kamiokande], Phys. Rev. D 94 (2016) no.5, 052001 doi:10.1103/PhysRevD.94.052001 [arXiv:1510.08127 [hep-ex]].
  * (15) K. Abe et al. [Hyper-Kamiokande], [arXiv:1805.04163 [physics.ins-det]].
  * (16) Z. Li, X. Huang, Q. Yuan and Y. Xu, JCAP 12 (2016), 028 doi:10.1088/1475-7516/2016/12/028 [arXiv:1312.7609 [astro-ph.CO]].
  * (17) C. M. Karwin, S. Murgia, I. Moskalenko, S. Fillingham, A. K. Burns and M. Fieg, Phys. Rev. D 103 (2021) no.2, 023027 doi:10.1103/PhysRevD.103.023027 [arXiv:2010.08563 [astro-ph.HE]].
  * (18) M. Ajello et al. [Fermi-LAT], Astrophys. J. 819 (2016) no.1, 44 doi:10.3847/0004-637X/819/1/44 [arXiv:1511.02938 [astro-ph.HE]].
  * (19) C. Karwin, S. Murgia, T. M. P. Tait, T. A. Porter and P. Tanedo, Phys. Rev. D 95 (2017) no.10, 103005 doi:10.1103/PhysRevD.95.103005 [arXiv:1612.05687 [hep-ph]].
  * (20) M. Ackermann et al. [Fermi-LAT], Astrophys. J. 840 (2017) no.1, 43 doi:10.3847/1538-4357/aa6cab [arXiv:1704.03910 [astro-ph.HE]].
  * (21) M. Di Mauro, Phys. Rev. D 103 (2021) no.6, 063029 doi:10.1103/PhysRevD.103.063029 [arXiv:2101.04694 [astro-ph.HE]].
  * (22) A. A. Abdo et al. [Fermi-LAT], Astrophys. J. 712 (2010), 147-158 doi:10.1088/0004-637X/712/1/147 [arXiv:1001.4531 [astro-ph.CO]].
  * (23) M. Ackermann et al. [Fermi-LAT], Phys. Rev. Lett. 107 (2011), 241302 doi:10.1103/PhysRevLett.107.241302 [arXiv:1108.3546 [astro-ph.HE]].
  * (24) A. Geringer-Sameth, S. M. Koushiappas and M. G. Walker, Phys. Rev. D 91 (2015) no.8, 083535 doi:10.1103/PhysRevD.91.083535 [arXiv:1410.2242 [astro-ph.CO]].
  * (25) M. Ackermann et al. [Fermi-LAT], Phys. Rev. Lett. 115 (2015) no.23, 231301 doi:10.1103/PhysRevLett.115.231301 [arXiv:1503.02641 [astro-ph.HE]].
  * (26) A. Drlica-Wagner et al. [Fermi-LAT and DES], Astrophys. J. Lett. 809 (2015) no.1, L4 doi:10.1088/2041-8205/809/1/L4 [arXiv:1503.02632 [astro-ph.HE]].
  * (27) A. Albert et al. [Fermi-LAT and DES], Astrophys. J. 834 (2017) no.2, 110 doi:10.3847/1538-4357/834/2/110 [arXiv:1611.03184 [astro-ph.HE]].
  * (28) F. Calore, P. D. Serpico and B. Zaldivar, JCAP 10 (2018), 029 doi:10.1088/1475-7516/2018/10/029 [arXiv:1803.05508 [astro-ph.HE]].
  * (29) S. Hoof, A. Geringer-Sameth and R. Trotta, JCAP 02 (2020), 012 doi:10.1088/1475-7516/2020/02/012 [arXiv:1812.06986 [astro-ph.CO]].
  * (30) M. Di Mauro and M. W. Winkler, Phys. Rev. D 103 (2021) no.12, 123005 doi:10.1103/PhysRevD.103.123005 [arXiv:2101.11027 [astro-ph.HE]].
  * (31) A. McDaniel, M. Ajello, C. M. Karwin, M. Di Mauro, A. Drlica-Wagner and M. A. Sánchez-Conde, Phys. Rev. D 109 (2024) no.6, 063024 doi:10.1103/PhysRevD.109.063024 [arXiv:2311.04982 [astro-ph.HE]].
  * (32) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, Comput. Phys. Commun. 185 (2014), 2250-2300 doi:10.1016/j.cpc.2014.04.012 [arXiv:1310.1921 [hep-ph]].
  * (33) S. L. Chen, Z. Kang, Z. K. Liu and P. Zhang, [arXiv:2405.05694 [hep-ph]].
  * (34) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, JHEP 07 (2014), 079 doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * (35) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen and P. Z. Skands, Comput. Phys. Commun. 191 (2015), 159-177 doi:10.1016/j.cpc.2015.01.024 [arXiv:1410.3012 [hep-ph]].
  * (36) G. D. Martinez, J. S. Bullock, M. Kaplinghat, L. E. Strigari and R. Trotta, JCAP 06 (2009), 014 doi:10.1088/1475-7516/2009/06/014 [arXiv:0902.4715 [astro-ph.HE]].
  * (37) L. E. Strigari, Phys. Rept. 531 (2013), 1-88 doi:10.1016/j.physrep.2013.05.004 [arXiv:1211.7090 [astro-ph.CO]].
  * (38) G. Cowan, K. Cranmer, E. Gross and O. Vitells, Eur. Phys. J. C 71 (2011), 1554 [erratum: Eur. Phys. J. C 73 (2013), 2501] doi:10.1140/epjc/s10052-011-1554-0 [arXiv:1007.1727 [physics.data-an]].
  * (39) B. S. Acharya et al. [CTA Consortium], WSP, 2018, ISBN 978-981-327-008-4 doi:10.1142/10986 [arXiv:1709.07997 [astro-ph.IM]].
  * (40) A. Acharyya et al. [CTA], JCAP 01 (2021), 057 doi:10.1088/1475-7516/2021/01/057 [arXiv:2007.16129 [astro-ph.HE]].
  * (41) P. Athron, J. M. Cornell, F. Kahlhoefer, J. Mckay, P. Scott and S. Wild, Eur. Phys. J. C 78, no.10, 830 (2018).