

\section{Introduction}

The nature of dark matter (DM) represents a fundamental question for the understanding of our Universe.
Observational hints at cosmological and galaxy scales such as the discrepancy between the measured rotation curves of galaxies and their theoretical predictions, the formation of large structures, and the anisotropies of the Comic Microwave Background show that DM makes up about 85\% of the total matter.

The search for DM has therefore become a priority in the scientific community where a collective effort has been made in indirect, direct, and collider searches in order to unravel its mystery. Its detection would also be a milestone in searches for Physics beyond the Standard Model (SM). 
In this talk, we focus on the indirect detection using the observations made by five gamma-ray experiments towards twenty dwarf spheroidal galaxies (dSphs).
These dSphs represent one of the most promising targets for DM indirect searches due to their high DM content and their negligible astrophysical background~\cite{PhysRevD.69.123501}.
They are all located at high Galactic latitude and no sign of very high energy emission has been detected so far in the dSphs' directions.
Gamma rays have the advantage of being neutral and do not get deflected by magnetic fields. Thus, the regions of $\gamma$-ray production can be traced back from the incident direction. The observations are therefore performed based on this property where the telescopes are directly pointing to the sources.

This work represents a collective effort between three imaging atmospheric Cherenkov telescope (IACT) arrays  H.E.S.S., MAGIC, and VERITAS, the water Cherenkov array HAWC, and the space-borne telescope {\it Fermi}-LAT, which agreed on sharing their data previously published individually. 
The goal of our study is to combine the individual upper limits published by each collaboration in order to optimize the statistics and increase the sensitivity to potential DM signals. The combination brings the novelty of extending the upper and lower boundaries of the energy range and the derivation of the upper limits on the DM annihilation cross-section over the widest DM particle mass range ever.
In this work, each of the five collaborations performed the analysis of their own data sets using a common DM model to optimize the data handling at different energy, angular resolutions, and sensitive energy ranges of the various instruments.
By following this procedure, we also avoid the need for sharing raw data and instrument response functions (IRFs) outside the collaborations.
As no significant excess was detected from the selected sources, nor in their combination, we derive upper limits on the DM annihilation cross-section as a function of the DM particle mass by combining the likelihood functions of all dSphs and all experiments.



\section{Experiments} \label{Sec:Experiments}

\subsection{Fermi-LAT}

The {\it Fermi}-Large Area Telescope ({\it Fermi-LAT}) is the collaboration which operates the pair conversion Large Area Telescope (LAT) carried by the Fermi satellite orbiting the Earth at an altitude of 565 km since 2008. The telescope has a wide field of view covering about 20\% of the sky and scans the whole sky every 3 hours in the energy range between $\sim 20$ MeV and 1~TeV. {\it Fermi}-LAT thus covers the lowest energy region of this study. Detailed descriptions of the detector and its performance can be found in \cite{2012ApJS..203....4A}.

\subsection{HAWC}

The High-Altitude Water Cherenkov (HAWC) detector is a high-energy $\gamma$-ray telescope located at Sierra Negra, Mexico at 4100~m altitude and consists of an array consisting of 300 water Cherenkov detectors (WCD) covering an area of 22,000~$\rm{m}^2$. The WCD are sensitive to $\gamma$-ray events of energies ranging from 300 GeV to a couple hundred TeV \cite{abeysekara}. The experiment covers a field of view of 15\% of the sky at all times.

\subsection{H.E.S.S.}

The High Energy Stereosocpic System (H.E.S.S.) experiment is an array consisting of five IACTs designed to detect brief and faint flashes of Cherenkov radiation generated by very high energy $\gamma$ rays between $\sim$30 GeV and $\sim 100$ TeV. The telescope array is located in central Namibia in the Khomas Highland region at 1,800 m above sea level~\cite{Aharonian:2006pe} at 110 km south west of Windhoek. The four small telescopes are equipped with a 12~m reflector while the central one is 28~m. The array collects the $\gamma$ rays within a field of view of 5$^\circ$.

\subsection{MAGIC}
The Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescope array consists of two telescopes of 17 m diameter reflector situated at the Roque de los Muchachos Observatory on the Canary Island of La Palma, Spain, 2,200~m above sea level. MAGIC is sensitive to very high energy $\gamma$-ray events above $\sim 50$ GeV~\cite{Sitarek} and is equipped with fast imaging cameras with a field of view of 3.5\degree.

\subsection{VERITAS}
The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is an array of four telescopes of 12~m reflector located at the Fred Lawrence Whipple Observatory in Southern Arizona. The telescope array is sensitive to a very high energetic band from $\sim$ 85~GeV up to $\sim 30$~TeV whose events are recorded within a $3.5\degree$ field of view~\cite{Park}.



\section{DM signal} \label{sec:DM}

The differential flux of $\gamma$ rays from the self-annihilation of Majorana DM particles is given by:
\begin{equation}
    \label{eq:DM_flux}
    \frac{\text{d}^2\Phi \left(\langle \sigma v \rangle, J \right)}{\text{d}E \text{d}\Omega} = \frac{1}{4\pi} \frac{\langle \sigma v \rangle}{2m^{2}_{\chi}} \sum_{f} \text{BR}_{f}\frac{\text{d}N_{f}}{\text{d}E} \times \frac{ \text{d}J}{\text{d}\Omega}.
\end{equation}

The first term contains the mass $m_{\chi}$ of the DM particles in GeV and their annihilation cross-section averaged over the velocity distribution $\langle \sigma v \rangle$ in $\rm{cm}^3\rm{s}^{-1}$. It also carries the differential spectrum $dN_f /dE$ for a given annihilation channel $f$.
Since we do not assume any specific particle physics model, each channel is treated individually where the branching ratio $\text{BR}_{f} = 100\%$.
The second term, known as the astrophysical $J$ factor, describes the amount of DM annihilations within a source or a region of the sky. 

The differential $J$-factor is defined as the integral of the square of the DM density distribution $\rho_{\rm{DM}}$ along the line-of-sight (l.o.s.):
\begin{equation}
  \frac{\text{d}J}{\text{d}\Omega} = \displaystyle{\int_\mathrm{l.o.s.} \: \rho_{\textrm{DM}}^2(r(s,\theta)) \:ds },
  \label{dwarf_J}
\end{equation}
where $\rho_{\textrm{DM}}$ is assumed to be spherically symmetric for all considered dSphs and depends on the distance to the center of the source $r$. This distance can also be expressed in terms of the distance $s$ from Earth along the line of sight, and the angular distance $\theta$ with respect to the center of the source, as $r(s, \theta) = (s^2 +d^2 -2sd\cos\theta)^{1/2}$, where $d$ is the distance between the Earth and the source. 
The $J$ factor computation is usually performed through a Jeans analysis based on the spherical Jeans equations~\cite{Bonnivard:2014kza, Geringer-Sameth:2014yza, binney2011galactic, Bonnivard:2015xpq}.
This technique relies on the spectroscopic data and assumes that the galaxies are in steady-state hydrodynamic equilibrium, have a spherical symmetry, and are non-rotating systems to reconstruct the galactic dynamics.
In this work, we use the $J$ factors produced by Geringer-Sameth \emph{et al.}~\cite{Geringer-Sameth:2014yza}.








\section{Joint likelihood analysis}\label{Sec:Likelihood}

\subsection{Dataset}

Twenty classical and ultrafaint dwarf spheroidal galaxies are selected for the combination. All were observed by one or more instruments and previously published by individual collaborations. Table~\ref{tab:dsphs_list} presents the list of dwarf galaxies used this project and the experiments with which they were observed.

\begin{table}[h!]
\small{
\centering{
\begin{tabular}{cccc}
\hline
\hline
Source name          & Experiments  & Distance & $\log_{10}J$           \\
 & & \scriptsize{(kpc)} &  \scriptsize{$\log_{10}(\rm{GeV}^2 \rm{cm}^{-5}\rm{sr})$}   \\
\hline
\CellTopTwo{}
Bootes I             & {\it Fermi}-LAT, HAWC, VERITAS &$66$ & $18.24^{+0.40}_{-0.37}$  \\
\CellTopTwo{}
Canes Venatici I     & {\it Fermi}-LAT   &$218$ & $17.44^{+0.37}_{-0.28}$           \\
\CellTopTwo{}
Canes Venatici II     & {\it Fermi}-LAT, HAWC & $160$& $17.65^{+0.45}_{-0.43}$    \\
\CellTopTwo{}
Carina     & {\it Fermi}-LAT, H.E.S.S. & $105$&  $17.92^{+0.19}_{-0.11}$  \\
\CellTopTwo{}
Coma Berenices     & {\it Fermi}-LAT, HAWC, H.E.S.S., MAGIC & $44$&  $19.02^{+0.37}_{-0.41}$   \\
\CellTopTwo{}
Draco     & {\it Fermi}-LAT, HAWC, MAGIC, VERITAS & $76$ &  $19.05^{+0.22}_{-0.21}$  \\
\CellTopTwo{}
Fornax    & {\it Fermi}-LAT, H.E.S.S.   & $147$ & $17.84^{+0.11}_{-0.06}$  \\
\CellTopTwo{}
Hercules   & {\it Fermi}-LAT, HAWC  & $132$ & $16.86^{+0.74}_{-0.68}$  \\
\CellTopTwo{}
Leo I     & {\it Fermi}-LAT, HAWC  & $254$ & $17.84^{+0.20}_{-0.16}$  \\
\CellTopTwo{}
Leo II     & {\it Fermi}-LAT, HAWC  & $233$ & $17.97^{+0.20}_{-0.18}$   \\
\CellTopTwo{}
Leo IV   & {\it Fermi}-LAT, HAWC  & $154$&  $16.32^{+1.06}_{-1.70}$  \\
\CellTopTwo{}
Leo T     & {\it Fermi}-LAT  &$417$& $17.11^{+0.44}_{-0.39}$    \\
\CellTopTwo{}
Leo V     & {\it Fermi}-LAT & $178$&  $16.37^{+0.94}_{-0.87}$  \\
\CellTopTwo{}
Sculptor     & {\it Fermi}-LAT, H.E.S.S.  & $86$& $18.57^{+0.07}_{-0.05}$  \\
\CellTopTwo{}
Segue I     & {\it Fermi}-LAT, HAWC, MAGIC, VERITAS   & $23$ &  $19.36^{+0.32}_{-0.35}$  \\
\CellTopTwo{}
Segue II     & {\it Fermi}-LAT   & $35$ & $16.21^{+1.06}_{-0.98}$ \\
\CellTopTwo{}
Sextans    & {\it Fermi}-LAT, HAWC &$86$ & $17.92^{+0.35}_{-0.29}$  \\
\CellTopTwo{}
Ursa Major I     & {\it Fermi}-LAT, HAWC  & $97$ & $17.87^{+0.56}_{-0.33}$  \\
\CellTopTwo{}
Ursa Major II     & {\it Fermi}-LAT, HAWC, MAGIC  &$32$& $19.42^{+0.44}_{-0.42}$   \\
\CellTopTwo{}
Ursa Minor     & {\it Fermi}-LAT, VERITAS  & $76$ & $18.95^{+0.26}_{-0.18}$  \\
\hline
\hline
\end{tabular}
}
\caption{Summary of the relevant properties of the dSphs included in the combination of {\it Fermi}-LAT, HAWC, H.E.S.S., MAGIC, and VERITAS likelihood functions. The list of the observed dwarf galaxies is presented in column 1 with the instruments that performed the observations in column 2. Their heliocentric distance and  $J$ factor with their estimated $\pm 1\sigma$ uncertainties are listed in columns 3 and 4 respectively. The $J$ factors are given for a source extension truncated at the outermost observed star with their estimated $\pm 1\sigma$ uncertainties. }
\label{tab:dsphs_list}}
\end{table}


\subsection{Combination principle}

Our search for DM is carried out using a technique of maximum likelihood in which the profile likelihood ratio $\lambda$ is a function of the annihilation cross-section $\langle \sigma v \rangle$, \emph{i.e.} the parameter of interest, and reads as:

\begin{equation}
    \label{Eq:likelihood_profile}
\lambda \left( \langle \sigma \mathit{v} \rangle \mid \boldsymbol{\mathcal{D}_{\text{dSphs}}} \right) = \frac{\mathcal{L} \left( \langle \sigma \mathit{v} \rangle ; \boldsymbol{\hat{\hat{\nu}}} \mid \boldsymbol{\mathcal{D}_{\text{dSphs}}} \right)}{\mathcal{L} \left( \widehat{\langle\sigma v \rangle} ; \boldsymbol{\hat{\nu}} \mid \boldsymbol{\mathcal{D}_{\text{dSphs}}} \right)} ,
\end{equation}
where $\boldsymbol{\mathcal{D}_{\text{dSphs}}}$ is the dataset, $\boldsymbol{\nu}$ represents the nuisance parameters, $\widehat{\langle\sigma v \rangle}$ and $\boldsymbol{\hat{\nu}}$ are the values that maximize $\mathcal{L}$ globally, and $\boldsymbol{\hat{\hat{\nu}}}$ the values that maximize $\mathcal{L}$ for a given value of $\langle \sigma \mathit{v} \rangle$. 
The joint likelihood function $\mathcal{L}$ describing all measurements is the product of the individual likelihood functions of all instruments and all dSphs and is given by:
\begin{equation}
    \label{Eq:dSph_combination}
    \mathcal{L} \left( \langle \sigma \mathit{v} \rangle ; \boldsymbol{\nu} \mid \boldsymbol{\mathcal{D}_{\text{dSphs}}} \right) = \prod_{l=1}^{N_{\text{dSphs}}} \mathcal{L}_{\text{dSph},l} \left( \langle \sigma \mathit{v} \rangle ; J_{l}, \boldsymbol{\nu_{l}} \mid \boldsymbol{\mathcal{D}_{l,\text{measured}}} \right) \times \mathcal{J}_l \left( J_{l} \mid J_{l,\text{obs}}, \sigma_{\log{J_l}} \right).
\end{equation}
The quantity $ N_{\text{dSphs}} = 20 $ is the total number of dSphs; $ \boldsymbol{\mathcal{D}_{l,\text{measured}}} $ is the dataset from gamma-ray observations for the $l$-th dSph; $ \boldsymbol{\nu_{l}} $ is the set of nuisance parameters associated to the $l$-th dSph, excluding $ J_{l} $; and $ J_{l} $ is the total $J$ factor of the $l$-th dSph, whose value can be found in Tab.~\ref{tab:dsphs_list}
; $\log_{10} J_{l,\text{obs}}$ and $ \sigma_{\log{J_l}} $ are obtained from the fit (see Jeans analysis in Sec.~\ref{sec:DM}) of a log-normal function of $J_{l,\text{obs}}$ to the posterior distribution of $ J_{l} $ \cite{2015PhRvL.115w1301A}. 

\subsection{Shared data format}
In order to perform the combination of the observations, a table of test statistic (TS) values is provided by each experiment for the annihilation channels, $b\bar{b}$ and $\tau^+\tau^-$, for each set of $m_{\chi}$ and $\langle \sigma v \rangle$. All collaborations agreed on 63 DM masses ranging from 10~GeV to 100~TeV for all continuum channels following the mass spacing of~\cite{Cirelli} to avoid an interpolation. The $\langle \sigma v \rangle$ range is defined between $10^{-28}~\rm{cm}^3\rm{s}^{-1}$ and $10^{-18}~\rm{cm}^3\rm{s}^{-1}$ and is logarithmically spaced in 1001 values. 








\subsection{Statistical uncertainty bands}
The $68\%$ (1$\sigma$) and $95\%$ (2$\sigma$) containment bands are derived by individual experiments by performing 300 Poisson realizations of the background events. Each collaboration provides the results of their statistical uncertainties in the same format as for the nominal values which are then combined following the same procedure as the combination of the nominal upper limits.




\section{Results and discussion}\label{Sec:Results}
No significant DM signal has been observed by any of the five instruments.
We therefore present the results of the combined upper limits at 95\% C.L. on the DM annihilation cross-section $\langle \sigma v \rangle$ in the case of two annihilation channels, $b\bar{b}$ and $\tau^+ \tau^-$, using all the data collected towards the twenty dSphs. We note that we selected these hadronic and leptonic channels as the follow up of our previous results presented at ICRC 2019~\cite{Oakes:2019ywx}.
We set our upper limits by solving $\mathrm{TS} = -2 \ln{ \lambda(\langle \sigma \mathit{v} \rangle)}$ for $\langle \sigma \mathit{v} \rangle$, with $\mathrm{TS}= 2.71 $. The value 2.71 represents the 95\% confidence level of a one-sided distribution assuming the test statistics behaves like a $\chi^2$ distribution with one degree of freedom. The combination is performed using two independent public analysis software packages, \texttt{gLike} \cite{javier_rico_2021_4601451} and \texttt{LklCombiner} \cite{tjark_miener_2021_4597500}, that provide compatible results.
The combined upper limits are presented in Fig.~\ref{fig:bb_tautau_UL} and are given with their 68\% ($1\sigma$) and 95\% ($2\sigma$) containment bands. These limits (solid black lines) are expected to be close to the median limit (dashed black lines) as no signal is present. 
We obtain upper limits within the 2 $\sigma$ expected bands for the two annihilation channels $b\bar{b}$ and $\tau^+ \tau^-$.
The individual limits produced by each experiment are also indicated in the figures as a comparison to our new combined results.
Below \textasciitilde500~GeV, the DM limits are largely dominated by the {\it Fermi}-LAT experiment.
Between \textasciitilde500 GeV to \textasciitilde10 TeV, {\it Fermi}-LAT continues to dominate for the hadronic DM channel then above \textasciitilde10 TeV, the IACTs (H.E.S.S., MAGIC, and VERITAS) and HAWC take over.
In the case of the leptonic channel, both the IACTs and HAWC contribute significantly to the DM limit from \textasciitilde1 TeV to \textasciitilde100 TeV.

\vspace{-0.2cm}

\begin{figure}[h!]
\includegraphics[width=0.93\textwidth]{gd_ICRC2021_v3_crop.png}
\vspace{-0.1cm}

\caption{Upper limits at 95\% confidence level on $\langle \sigma \mathit{v} \rangle$ as a function of the DM mass for the annihilation channels $b\bar{b}$ (left) and $\tau^+ \tau^-$ (right), using the set of $J$ factors from Ref.~\cite{Geringer-Sameth:2014yza}. The black solid line represents the observed combined limit, the black dashed line is the median of the null hypothesis corresponding to the expected limit, while the green and yellow bands show the 68\% and 95\% containment bands. Combined upper limits for each individual detector are also indicated as solid, colored lines. }

\label{fig:bb_tautau_UL}
\end{figure}

We observe that the combined DM constraints from all five telescopes are 2 to 3 times stronger than any individual telescope for multi-TeV DM.
The selection of multiple targets increases statistics used to probe these sources and allows us to derive upper limits spanning the largest mass range of any WIMP DM search. We note that these limits depend on the choice of the annihilation channels and are driven by the objects with the highest $J$ factors that can be observed. The ultrafaint dSphs, containing a few tens of bright stars only, can be subject to large systematic uncertainties for the determination of their $J$-factors such as Segue I.
The derivation of upper limits through 6 additional annihilation channels is currently in progress, with 5 other continuum channels and the monoenergetic channel $\gamma \gamma$. A further analysis using a second $J$ factor set derived by~\cite{Bonnivard:2014kza, Bonnivard:2015xpq} is also yet to come in order to study the systematics induced by the choice of $J$ factor.







{\small
\bibliographystyle{ieeetr}
