Dark Matter Constraints from a Joint Analysis of Dwarf Spheroidal Galaxy Observations with VERITAS

TL;DR

This paper uses joint analysis of gamma-ray observations from four dwarf galaxies with VERITAS to set new upper limits on WIMP dark matter annihilation cross sections, improving detection sensitivity through spectral and spatial photon analysis.

Contribution

It introduces an optimized photon weighting technique for IACT data, enhancing dark matter search sensitivity by combining multiple dwarf galaxy observations.

Findings

No gamma-ray emission detected from individual or joint analysis.

Upper limits on dark matter annihilation cross section are established at various final states.

Method improves constraints compared to previous analyses.

Abstract

We present constraints on the annihilation cross section of WIMP dark matter based on the joint statistical analysis of four dwarf galaxies with VERITAS. These results are derived from an optimized photon weighting statistical technique that improves on standard imaging atmospheric Cherenkov telescope (IACT) analyses by utilizing the spectral and spatial properties of individual photon events. We report on the results of $\sim$230 hours of observations of five dwarf galaxies and the joint statistical analysis of four of the dwarf galaxies. We find no evidence of gamma-ray emission from any individual dwarf nor in the joint analysis. The derived upper limit on the dark matter annihilation cross section from the joint analysis is $1.35\times 10^{-23} {\mathrm{ cm^3s^{-1}}}$ at 1 TeV for the bottom quark ($b\bar{b}$) final state, $2.85\times 10^{-24}{\mathrm{ cm^3s^{-1}}}$ at 1 TeV for the tau lepton ($\tau^{+}\tau^{-}$) final state and $1.32\times 10^{-25}{\mathrm{ cm^3s^{-1}}}$ at 1 TeV for the gauge boson ($\gamma\gamma$) final state.

Figures

Tables

Table 1. Table 1: Dwarf galaxy zenith and azimuth range, total accumulated exposure and energy range after cuts are applied. Azimuth is measured east from north. Upper energy range is defined as energies where the uncertainty in effective area is less than 10%.

Dwarf	Zenith	Azimuth	Exposure	Energy Range
	[deg]	[deg]	[hours]	[GeV]
Segue 1	15-35	100-260	92.0	80 - 50000
Draco	25-40	320-40	49.8	120 - 70000
Ursa Minor	35-45	340-30	60.4	160 - 93000
Boötes 1	15-30	120-249	14.0	100 - 41000
Willman 1	20-30	340-40	13.6	100 - 43000

Table 2. Table 2: Dwarf galaxy detection significance (generalized Li & Ma method) and integral flux upper limit with 95% confidence level above 300 GeV, assuming a spectral index of -2.4. The last two columns are the heliocentric distance to each object and the inferred value of J 𝐽 J -profile integrated within a cone with half-angle of 0.17 ∘ superscript 0.17 0.17^{\circ} (i.e. over the ON region), errors denote the 16th and 84th percentiles on the posterior Geringer-Sameth et al. ( 2015a ) . Note that this analysis uses the J 𝐽 J -profile convolved with the VERITAS instrument response as discussed in Section IV.

Dwarf	$N_{ON}$	$N_{OFF}$	$\overline{\alpha }$	Significance	$N^{95\%}$	${\Phi }^{95\%}$	$D$	${\log }_{10}J\left({0.17}^{\circ }\right)$
	[counts]	[counts]		[$\sigma$]	[counts]	[10-12cm2s-1]	[kpc]	[GeV2 cm-5]
Segue 1	15895	120826	0.131	0.7	235.8	0.34	23	19.2${}_{-0.3}{}^{}{}_{}{}^{+0.3}$
Draco	4297	39472	0.111	-1.0	33.5	0.15	76	18.3${}_{-0.1}{}^{}{}_{}{}^{+0.1}$
Ursa Minor	4181	35790	0.119	-0.1	91.6	0.37	76	18.9${}_{-0.3}{}^{}{}_{}{}^{+0.3}$
Boötes 1	1206	10836	0.116	-1.0	34.5	0.40	66	18.3${}_{-0.4}{}^{}{}_{}{}^{+0.3}$
Willman 1	1926	18187	0.108	-0.6	23.5	0.39	38	N/A

Equations

$$\frac{dJ(\hat{\mathbf{n}})}{d\Omega }=\int {\rho }^2(\ell \hat{\mathbf{n}})d\ell ,$$

$$\rho (r)={\rho }_s[r/r_s{]}^{-\gamma }[1+(r/r_s{)}^{\alpha }{]}^{(\gamma -\beta )/\alpha },$$

$$T=\sum _i{w}_i,$$

$$w=\log \left[1+\frac{s}{b}\right],$$

$$\frac{dF(E,\hat{\mathbf{n}})}{dEd\Omega }=\frac{⟨\sigma v⟩}{8\pi M^2}\frac{d{N}_{\gamma }(E)}{dE}\frac{dJ(\hat{\mathbf{n}})}{d\Omega },$$

$$\frac{d{N}_{\gamma }(E)}{dE}=\sum _i{B}_i\frac{d{N}_{\gamma ,i}(E)}{dE},$$

$$\frac{dN(E,\hat{\mathbf{n}})}{dEd\Omega }=\underset{E_t}{\int }\underset{{\Omega }_t}{\int }d{E}_{t}d{\Omega }_t\frac{dF(E_t,{\hat{\mathbf{n}}}_t)}{d{E}_{t}d{\Omega }_t}R(E,\hat{\mathbf{n}}|E_t,{\hat{\mathbf{n}}}_t),$$

$$R(E,\hat{\mathbf{n}}|E_t,{\hat{\mathbf{n}}}_t)=\sum _{\text{runs }}\tau A_{\mathrm{eff}}(E_t)\mathrm{PSF}(\hat{\mathbf{n}}|E_t,{\hat{\mathbf{n}}}_t)D(E|E_t).$$

$$s(\nu ,E,\theta )=\frac{dN(\nu ,E,\theta )}{dEd\Omega }dE2\pi \sin (\theta )d\theta ,$$

$${\Phi }_{\gamma }(>E_{\min })$$

Full text

I Introduction

The search for standard model particles resulting from the annihilation of dark matter particles provides an important complement to the efforts of direct searches for dark matter interactions and searches for dark matter production at particle accelerators. Among the theoretical candidates for the dark matter particle above a few GeV, Weakly Interacting Massive Particles (WIMPs) are well motivated Jungman et al. (1996); Servant & Tait (2003) as they naturally provide the measured present day cold dark matter density Planck Collaboration et al. (2014); Steigman (1979); Zeldovic et al. (1965); Zel’dovich (1965); Chiu (1966). In such models, the WIMPs either decay or annihilate into standard model particles that produce mono-energetic gamma-ray lines and/or a continuum of gamma rays with energies up to the dark matter particle mass.

Attractive targets for indirect dark matter searches are nearby massive objects with high inferred dark matter content and that are not expected to be sources of very-high-energy gamma rays. Dwarf spheroidal galaxies (dSphs) are relatively close ($\sim$20 to 200 kpc) to Earth and lack conventional astrophysical high-energy sources of gamma rays Mateo (1998). Five dwarf galaxies have been observed with the Very Energetic Radiation Imaging Telescope Array System (VERITAS) between 2007 and 2013, for a total of 230 hours of high quality data.

In this paper we perform a joint statistical analysis of dwarf galaxies observed with VERITAS. We find no evidence of dark matter annihilation in any of the dwarf galaxies individually observed with VERITAS or in a joint analysis of four of the dwarfs. We place upper limits on the emitted flux and derive upper limits on the annihilation cross section.

II Observations

VERITAS is an array of four imaging atmospheric Cherenkov telescopes (IACTs), each 12 m in diameter, located at the Fred Lawrence Whipple Observatory in southern Arizona, USA (31.68*∘* N, 110.95*∘* W, 1.3 km above sea level). Each VERITAS camera contains 499 pixels (0.15*∘* diameter) and has a field of view of 3.5*∘. VERITAS began full array operations in the spring of 2007. The instrument has gone through a number of upgrades since then to improve performance. In the summer of 2009, the first telescope (“T1”) was moved to its current location in the array to provide a more uniform distance between telescopes, improving the sensitivity of the system Perkins et al. (2009). The telescope-level trigger was replaced with a faster system in the fall of 2011 Zitzer & for the VERITAS Collaboration (2013b), allowing for greater night-sky background (NSB) reduction during all operating modes of the experiment. The VERITAS camera pixels were replaced in summer 2012 with higher quantum efficiency photomultiplier tubes (PMTs), allowing for a lowered energy threshold D. B. Kieda for the VERITAS Collaboration (2013). VERITAS is sensitive to gamma rays from approximately 85 GeV (after camera upgrade) to greater than 30 TeV with a typical energy resolution of $15-25\%$ and an angular resolution (68% containment) of $<$0.1 degrees per event. The flux sensitivity of the standard analysis is such that a source with a flux of order of 1% of the Crab Nebula flux can be detected in approximately 25 hours of observation. The looser event selection criteria (commonly referred to as “cuts”) used in this work described later in this section resulted in a slightly larger energy resolution (25%-30% at 1 TeV) and angular resolution ($\sim$0.12∘* at 1 TeV).

From the beginning of four-telescope operations in 2007 to the summer of 2013, five dwarf galaxies in the northern hemisphere have been observed by VERITAS: Segue 1, Ursa Minor, Draco, Boötes and Willman 1. Quality data for this analysis requires moonless and clear atmospheric (based on infrared temperature measurements) conditions and operation of all four telescopes. Dwarf galaxy data used in this work were taken during three different epochs of VERITAS operations: data taken before the move of T1, data taken after the T1 move, and data taken after the camera upgrade. In all three epochs, data were obtained with the wobble pointing strategy, where the camera center is offset by 0.5 degrees from the target position Fomin et al. (1994). The wobble mode allows for simultaneous background estimation and source observation, reducing the systematic uncertainties in the background estimation as opposed to using separate pointings for background estimation.

The data reduction mostly follows the standard techniques employed by VERITAS Holder et al. (2006), with the notable exceptions being the methodology of the cosmic-ray background estimate, the adopted statistical approach based on individual photon weighting, and the method of image characterization for shower reconstruction. Images recorded by the VERITAS cameras are calibrated by the photomultiplier tube (PMT) gains. Traditionally the showers are characterized by their second moments Hillas (1985). In this work each Cherenkov shower image is fit with a two-dimensional elliptical Gaussian function to get the parameter characterization of the shower Christiansen & VERITAS Collaboration (2012). This fitting method for Cherenkov images is advantageous because the two-dimensional elliptical Gaussian fit allows for better point-spread function (PSF) characterization at high energies, and is less biased to images that are truncated at the edge of the camera or by dead pixels or suppressed pixels due to bright stars. This method of fitting has also been shown to reduce the time for a weak point source to reach 5$\sigma$ by 20% Christiansen & VERITAS Collaboration (2012). The stereo reconstruction of the event’s arrival direction and energy is accomplished by combining parameters from multiple telescopes Krawczynski et al. (2006). The hadronic cosmic-ray background is reduced by applying mean scaled width and mean scaled length cuts Krawczynski et al. (2006). The cuts were optimized a priori using data from known weak and soft-spectral very-high-energy sources. These “soft” cuts were selected to give the lowest possible energy threshold, which increases sensitivity to dark matter searches by allowing more low energy events to be used for the analysis. An additional cut is applied on the angle between the target position and the reconstructed arrival position, $\theta<0.17$ degrees, thus defining the signal search region or “ON region”.

Many IACT analyses select background events from one or more OFF regions in the camera field of view Berge et al. (2007). Two methods for forming an OFF region are commonly used. In the reflected region method (also called a wobble analysis), the source is offset from the telescope tracking position, and OFF regions consist of regions with the same size as the ON region with the same offset. In the ring background method the OFF region is an annulus surrounding the ON region.

This analysis requires a larger sample of the measured background and to determine its energy spectrum, therefore a third method is introduced. We name this new method the “crescent” background method (CBM) Zitzer & for the VERITAS Collaboration (2013a). This method was previously described in Berge et al. (2007) but this is the first time it has been applied to IACT data. Background events are selected from an annulus similar to the ring background. However, the annulus is centered on the tracking position as opposed to the source position (see Figure 1). This gives roughly a factor of two more background events than from standard reflected regions (depending on the field of view of the array pointing). The ring background method typically used is not suitable for this analysis, due to the energy dependence in IACT acceptances. Those acceptances are symmetrical around the tracking position to first order Berge et al. (2007). By selecting events only from a region at approximately the same angular distance from the tracking position, we reduce the energy dependence of the background scaling factor, $\alpha$.

Visible starlight may bias the background estimate and is removed by defining circular background exclusion regions centered around stars with apparent magnitudes of $m_{\mathrm{B}}<8$. The size of the exclusion region used varies with the brightness of the star; for example an exclusion region of 0.4 degrees is set around the 3.5-apparent magnitude star $\eta$ Leonis in the field of Segue 1. The central region of radius 0.3 degrees around each dwarf is also excluded.

The scaling factor of each background event, $\alpha$, used to calculate the gamma-ray excess and significance Li & Ma (1983) is determined by the ratio of the integral of the cosmic-ray acceptance within the ON region to the integral of the acceptance within the crescent-shaped OFF region. To better account for background systematics associated with deep exposures, an acceptance function was derived using the zenith angle of observation as well as the angular distance from the tracking direction. The procedure is similiar to the one described in the appendix of Rowell (2003) and is described in more detail in B. Zitzer for the VERITAS Collaboration (2015). An acceptance gradient in the VERITAS cameras was determined by utilizing a smoothed map of the ratio of counts using the total data set for each dSph in each skymap bin to the azimuthally-symmetric acceptance in that map bin, a parameter we refer to as flatness. If the radial-only acceptance adequately describes the cosmic-ray background, then the flatness map should be uniform within statistical errors across the field of view, i.e. it should not correlate with zenith angle or any other external parameters. A second map was produced with the mean difference of the zenith angle of the reconstructed photon direction from the zenith angle of the array tracking direction at the time the event was recorded. We will refer to this as the mean zenith map for simplicity. A scatter plot of the contents of each bin for the mean zenith map and the flatness map was made, showing a strong correlation for each field of view. That correlation was fit with a fourth-degree polynomial which was used to re-weight each bin in the spatial acceptance map and re-calculate $\alpha$. The difference between $\alpha$ with and without the zenith correction is $\lesssim$ 1%.

III Dark Matter Distribution within the dwarfs

The strength of the predicted gamma-ray signal is proportional to the dark matter distribution within dwarf galaxies. In general, this is characterized by the $J$-profile, defined as

[TABLE]

where $\ell$ is the line-of-sight distance along the $\hat{\mathbf{n}}$ direction, $d\Omega$ is the solid angle, and $\rho$ is the mass density profile of the dwarf galaxy.

The distribution of dark matter in dwarf galaxies is obtained using line-of-sight velocity and position measurements of stars that are gravitationally bound within the dwarf galaxy potential well Simon & Geha (2007); Walker, Mateo & Olszewski (2009). Distributions of stellar velocities and positions are functions of the gravitational potential as described by the Jeans equation Strigari et al. (2007, 2008); Binney & Tremaine (2008); Walker (2013); Battaglia et al. (2013); Strigari (2013).

We adopt the observational constraints on $J$-profiles as derived by  Geringer-Sameth et al. (2015a). The density profile of each dwarf is modeled as a “generalized” NFW (Navarro-Frenk-White) profile Zhao (1996),

[TABLE]

with five free parameters. A likelihood function relates the five parameters (and a sixth nuisance parameter specifying the stellar velocity anisotropy) to the observables through the Jeans equation. The parameter space is explored, giving rise to a chain of posterior sample halos.

This analysis generates many realizations of halos which reasonably fit the stellar kinematic data. This produces a systematic uncertainty for the dark matter search. When we present the results of the search and limits on the annihilation cross section we will separate this systematic uncertainty from the statistical uncertainty induced by our finite event statistics. This is done by repeating the analysis separately for different realizations of halo parameters. The systematic uncertainty “band” that results from this repetition should be thought of as reflecting our imperfect knowledge of the dwarf density profiles. See Section IX.C of  Geringer-Sameth et al. (2015b) for details.

Use of the Jeans equation requires the assumption that stellar tracers are in dynamical equilibrium and the analysis of Geringer-Sameth et al. (2015a) further assumes spherical symmetry, Plummer light profiles, and velocity anisotropy that is constant with radius. These are approximations, and all real systems will violate them at some level. Bonnivard et. al. Bonnivard et al. (2015) have studied the biases introduced by these effects. While the statistical uncertainty due to finite kinematic sample sizes dominates the errors in $J$ for ultrafaint dwarfs (e.g. Segue 1, Boötes 1, Willman 1), the assumption of spherical symmetry may cause a moderate bias (comparable to the statistical error bar) for the classical dwarfs (e.g. Draco, Ursa Minor). In the combined analysis, the uncertainties for Segue 1 dominate the error budget and our results will be insensitive to the other systematic effects mentioned above.

The stellar population of Willman 1 shows irregular kinematics, which may be due to ongoing tidal disruption of the satellite Willman et al. (2011). Regardless of the cause, the observations strongly suggest that Willman 1 is not in dynamical equilibrium, violating a core assumption of the Jeans equation. This object was excluded from the analysis of Geringer-Sameth et al. (2015a), who considered the inferred $J$-profile to be unreliable with no handle on the magnitude of the error. In the present work, we therefore exclude Willman 1 from results which require an estimate of its $J$-profile.

Additionally, Bonnivard et. al. Bonnivard et al. (2016) have pointed out the possibility of contamination of the stellar samples used to perform the Jeans analysis. Milky Way interlopers mistakenly included in the spectroscopic sample of dwarf member stars will inflate the inferred velocity dispersion and may bias $J$-profiles toward large expected annihilation signals. In particular, there are indications that Segue 1 may suffer from such contamination: the removal of several ambiguous stars from Segue 1 sample can have drastic (i.e. orders of magnitude) effects on $J$. Compared with classical dwarfs, this issue will be most severe for ultrafaint dwarfs, which have much smaller spectroscopic samples. While several groups have begun extending the Jeans analysis framework to encompass foreground contaminationBonnivard et al. (2016)Ichikawa et al. (2016)Zhu et al. (2016), no uniform analysis of the dwarf population has been performed, though several groups have begun extending the analysis framework to encompass this effect Ichikawa et al. (2016)Zhu et al. (2016). Notably, the issue of contamination has not been observationally checked for any ultrafaint dwarfs apart from Segue 1 and the recently discovered Reticulum II. Ichikawa et. al. Ichikawa et al. (2016), simulating future spectroscopic observations, find that contamination may bias $J$ high by factors of $\sim 3$ for the classical dwarfs Draco and Ursa Minor. Therefore, we caution that the uncertainties in our particle physics limits may be underestimated due to this additional astrophysical systematic uncertainty.

IV Event weighting

We employ a newly-developed event weighting method Geringer-Sameth et al. (2015b) to simultaneously analyze the data from all five dwarf fields. This technique improves on standard IACT analyses by utilizing the spectral and spatial properties of the individual events. It also takes into account the expected properties of the annihilation signal and the instrumental and astrophysical backgrounds, to perform an “optimal” analysis (see Geringer-Sameth et al. (2015b) for further details and a theoretical development of the technique).

Given the reconstructed events in an ON region we seek an optimal way to extract a possible dark matter signal. Each reconstructed event is assigned a weight based on three parameters: the dwarf field $\nu$ it came from, its reconstructed energy $E$, and its reconstructed angular separation from the dwarf galaxy $\theta$. The test statistic $T$ is defined as

[TABLE]

where the index $i$ runs over all ON events from all dwarf fields and $w_{i}=w(\nu_{i},E_{i},\theta_{i})$ is the weight of the $i$th event.

The weight function $w(\nu,E,\theta)$ can be an arbitrary function of the event properties. For example, a conventional ON/OFF analysis (see e.g. Aliu et al. (2012)) is recovered if $w=1$ for all events within the ON region of a particular dwarf and $w=0$ for all other events. In this case the test statistic is just the number of observed events in the ON region.

The weight function can be designed to distinguish, as efficiently as possible, the difference between background and background plus a dark matter signal. An intuitive solution is to weight different events according to how likely they are to be due to dark matter compared to background.

It has been shown Geringer-Sameth et al. (2015b) that when testing a simple null hypothesis (background only) against a simple alternative (signal plus background) the optimal form of the weight function $w(\nu,E,\theta)$ is

[TABLE]

where $s(\nu,E,\theta)$ is the expected number of signal events with properties $(\nu,E,\theta)$, and $b(\nu,E,\theta)$ is the expected number of background events due to all other processes besides dark matter annihilation (e.g. hadronic air showers, leptonic air showers and diffuse astrophysical gamma rays). The test statistic derived from this weighting is optimal in the sense that it maximizes the statistical power of the hypothesis test; if a dark matter signal is hidden in the data this test statistic is most likely to turn up a detection (see Geringer-Sameth et al. (2015b) for details).

The functions $s(\nu,E,\theta)$ and $b(\nu,E,\theta)$ are differential quantities, namely the expected number of events from dwarf $\nu$ with energies between $E$ and $E+dE$ and angular separations between $\theta$ and $\theta+d\theta$. We use the events in the OFF region of each dwarf to estimate the function $b$. The energy spectrum of these background events is modeled as a piecewise function. For energies below 1 TeV we replace each event with a Gaussian of width 3% of the measured energy, giving a kernel density estimate. This is a requirement of the kernel estimator and is unrelated to the VERITAS energy dispersion. Above 1 TeV we splice on a power law with exponential cutoff. The form is $f(E)=f_{0}(E/E_{0})^{\gamma}\exp[(E-E_{0})/E_{\mathrm{cut}}]$, where $E_{0}=\text{1 TeV}$ and $f_{0}$ is the kernel density estimate of the spectrum at 1 TeV. The choice of 3% of the measured energy as well as 1 TeV for the energy cutoff are arbitrary and do not affect the statistical significances of the search or the coverage of the limits. The parameters $\gamma$ and $E_{\mathrm{cut}}$ are obtained using the unbinned maximum likelihood. We choose this smooth fitting function to avoid noise in the kernel density estimator due to the relatively low number of observed events with high energies. The corrected solid angle ratios $\alpha$ between OFF and ON regions are used to predict the expected number of background events in the ON region for each dwarf. The background is assumed to be isotropic within the ON region so the $\theta$ dependence of $b(\nu,E,\theta)$ is proportional to $\sin(\theta)d\theta$.

The expected signal $s(\nu,E,\theta)$ is determined by convolving the dark matter annihilation flux with the VERITAS instrument response. The gamma-ray flux from annihilation, i.e. flux of photons from direction $\hat{\mathbf{n}}$ per energy per solid angle, is given by

[TABLE]

where $M$ is the dark matter particle mass, $\langle\sigma v\rangle$ is the velocity-averaged annihilation cross section, and $dN_{\gamma}/dE$ is the spectrum of gamma rays from a single annihilation event. This last spectrum is determined by the branching ratios $B_{i}$ into the various Standard Model final states:

[TABLE]

where $dN_{\gamma,i}/dE$ is the number of gamma rays produced per annihilation per gamma-ray energy by the products of channel $i$. We adopt the annihilation spectra given in Cirelli et al. (2011), including electroweak corrections. For annihilation into a two-photon final state we model the energy spectrum as a gaussian of width 10% of the dark matter mass and an amplitude of two photons. This width is always less than the VERITAS energy resolution.

The number of events reconstructed with energy $E$ and angular separation $\theta$ is given by the convolution

[TABLE]

where the subscript $t$ denotes true energies and directions and the function $R$ is the response of VERITAS. For clarity we have omitted a subscript $\nu$ from the quantities in Eq. 7, but the predicted dark matter flux and VERITAS response depend on which dwarf is being considered.

The response $R(E,\hat{\mathbf{n}}|E_{t},\hat{\mathbf{n}}_{t})dEd\Omega$ is the probability (per incident flux) that a gamma ray with true energy $E_{t}$ and direction $\hat{\mathbf{n}}_{t}$ will be reconstructed with an energy in the interval $dE$ around $E$ and in the solid angle $d\Omega$ around direction $\hat{\mathbf{n}}$. It is the product (summed over VERITAS observation runs) of the effective area $A_{\mathrm{eff}}$, live time per observation run $\tau$, instrument PSF, and energy dispersion $D$:

[TABLE]

These four factors are computed for each observation run. Because the considered $J$-profiles and PSFs are azimuthally symmetric in $\hat{\mathbf{n}}$ (i.e. $dJ/d\Omega$ only depends on the angle between $\hat{\mathbf{n}}$ and the dwarf and the PSF only depends on the angle between $\hat{\mathbf{n}}$ and $\hat{\mathbf{n}}_{t}$), the expected number of events is also azimuthally symmetric and depends only on $\theta$, the angle between the reconstructed direction $\hat{\mathbf{n}}$ and the direction of the dwarf.

The VERITAS point spread function, ${\rm PSF}(\theta|E_{t})$ (probability per solid angle of detecting a photon of true energy $E_{t}$ an angular distance $\theta$ away from its true direction) is derived from gamma-ray simulations. The reason that simulations were used instead of data from a bright source (for example, the Crab Nebula) is that simulations provide much larger statistics, and therefore better characterization at all energies. The simulated PSF agrees well with Crab Nebula data, to within $\lesssim$10% in the energy range where VERITAS is most sensitive. The same quality and background rejection cuts are applied to the simulated events, which are then binned in $\theta$ from 0*∘* to 2*∘* and in $E$ in the range from 0.01 TeV to 100 TeV, covering the entire VERITAS energy range. At each energy, the binned histogram is normalized over $\theta$, forming the probability distribution function, ${\rm PSF}(\theta|E_{t})$. The VERITAS epoch, the energy and the zenith angle are the only simulated parameters that have an impact on the shape of the PSF in this work, although others were investigated. Azimuthal angle and background noise dependencies have a negligible effect for this analysis. Examples of the energy dependence are shown in the left panel of Figure  2. The differences in the curves are due to differences in zenith angle and the epochs the dSphs were observed in.

The effective collection area, $A_{\mathrm{eff}}(E_{t})$ is a function of the true gamma-ray energy $E_{t}$, and it depends on the zenith and azimuth angles of observations, the amount of background noise present, VERITAS configuration epoch, offset of the source from the target position, and the gamma-ray cuts Mohanty et al. (1998). The right panel of Figure  2 depicts the average effective area curves of the observing conditions (zenith, azimuth, NSB and epochs) for all dwarf galaxies included in this study.

The line spread function, or energy dispersion ($D(E|E_{t})$) quantifies the energy resolution and bias of VERITAS. It is constructed by generating Monte Carlo gamma-ray showers at a true energy and putting the simulated showers through a simulated detector and the same reduction and cuts as the data. The shower reconstruction algorithm of the data analysis assigns the event a reconstructed energy $E$ Mohanty et al. (1998). Simulated showers that survive the “soft” cuts described above are put into a two dimensional histogram of reconstructed and true energy. Each bin of $E_{t}$ is normalized to unity to produce a probability density function.

Finally, the expected number of dark matter events from a dwarf with reconstructed energy between $E$ and $E+dE$ and separation between $\theta$ and $\theta+d\theta$ is simply

[TABLE]

with $dN/dEd\Omega$ given by Eq. (7).

To conduct a search for annihilation or set limits on the cross section we compute the probability distribution for measuring the test statistic under various hypotheses. For example, to conduct a search for dark matter annihilation, the observed value of the test statistic $T_{\mathrm{obs}}$ is compared with the probability distribution for $T$ due to background processes only $\mathrm{P}(T|\text{bg-only})$. The significance of the detection is defined as the probability that $T$ is less than $T_{\mathrm{obs}}$ under the background-only hypothesis. It is convenient to convert this probability into a “sigma value” using percentiles of a standard Gaussian distribution.

Alternatively, to construct upper limits on the annihilation cross section we compute the distribution for $T$ given a particular dark matter model, which includes specifying values for the particle mass $M$, cross section $\langle\sigma v\rangle$, and the branching fractions $B_{i}$ (see Eqs. (5) and (6)).

The method for computing the probability distribution for $T$ under any dark matter hypothesis (i.e. $\langle\sigma v\rangle\neq 0$), is detailed in Geringer-Sameth et al. (2015b). An abbreviated description follows. The test statistic is the sum of two independent quantities $T_{s}$ and $T_{b}$: the sum of the weights of events due to dark matter (signal) and all other sources (background). The weights of individual signal events are statistically independent and they are independent of the weights of background events. Further, in this study we assume that background events are all independent of each other.

Under these conditions, the variables $T_{s}$ and $T_{b}$ are described by compound Poisson distributions: the sum of independent random variables (the weights) where the number of terms in the sum is a Poisson distributed variable. All that is required to construct the distribution is the expected number of events that will be detected with each weight. This is found by discretizing the $(\nu,E,\theta)$ space in a finite number of bins and computing the expected number of events in each bin (Eq. 9) and the weight assigned to events in each bin (Eq. 4). Then a histogram is formed over the weight variable.

For the background events we consider the same discretized $(\nu,E,\theta)$ space. The weight of events in each bin is computed as above. The expected number of background events in each bin is computed using the empirical energy distribution of the OFF events and assuming the background events will be isotropic within the ON region. Specifically, each OFF event from dwarf $\nu$ with reconstructed energy in bin $E$ contributes $\alpha d\Omega_{j}/\Omega$ expected events to the $(\nu,E,\theta_{j})$ bin, where $\alpha$ is the ON/OFF ratio for the run, $d\Omega_{j}$ is the solid angle of the $j$-th $\theta$-bin, and $\Omega$ is the total solid angle of the ON region. This procedure is equivalent to a background model where events are sampled from OFF regions (with replacement) and distributed isotropically within the ON region; the probability of selecting an OFF event is proportional to its $\alpha$ value.

The probability distribution for $T$ is the convolution of the probability distributions for $T_{s}$ and $T_{b}$ (since $T=T_{s}+T_{b}$). The compound Poisson distributions and the convolutions are efficiently calculated using standard Fast Fourier Transform techniques.

In principle, the statistical power of the analysis can be increased by having an eventfls weight depend on the run in which it was detected (in addition to it’s energy, angular separation, and which dwarf field it was detected in). This generalization would automatically and optimally “downgrade” runs which had poor observing conditions (smaller effective area, larger background flux). However, this requires having accurate background models and response functions on a run by run basis and current datasets are not large enough to allow this. In general, the search becomes more sensitive as the event weights are allowed to depend on more observables.

V Results

V.1 Search for annihilation in individual dwarfs

The search for dark matter annihilation is performed by measuring $T_{\mathrm{obs}}$ and comparing this with the probability distribution for $T$ due to background. A search in an individual dwarf field is performed by setting the weights of events from all other dwarfs to zero. The weight function Eq. (4) requires a signal hypothesis $s(\nu,E,\theta)$ which depends on the dark matter parameters $M$, $\langle\sigma v\rangle$, and $B_{i}$. We perform a search for dark matter of each mass and annihilation channel (assuming $B_{i}=1$) in heavy quarks ($b\bar{b}$) and leptons ($\tau^{+}\tau^{-}$) as well as a two photon final state. The cross section $\langle\sigma v\rangle$ is a measure of the expected signal amplitude and must be specified in order to assign weights. A specific value $\langle\sigma v\rangle_{90}$ is used: it is the value of the cross section for which there is a 90% chance of making a 3$\sigma$ detection, where $\sigma$ is defined as number of standard deviations above the background. In VHE astronomy, 5$\sigma$ is typically required for a discovery. In practice, the search is essentially independent of the specific value of $\langle\sigma v\rangle$ used in the weighting, but $\langle\sigma v\rangle_{90}$ is chosen to make the search as sensitive as possible to cross sections that are on the verge of being detectable by the instrument.

Figure 3 shows the results for the search in the individual dwarfs. No evidence of dark matter annihilation at any mass has been observed in any one of the dwarfs. Note that annihilation into a two photon final state terminates at the highest energy of the event sample as shown in the last column of Table 1. These run from the lowest reconstructed energy for an off source event to an upper energy where the uncertainty in the effective area is 10%. The limits given here are insensitive to these energy thresholds.

V.2 Flux upper limits

Due to the lack of any detectable signal and in order to compare with complementary experiments we derive a flux upper limit $\Phi_{\gamma}(E>E_{\mathrm{min}})$, as

[TABLE]

where $N_{\gamma,{\mathrm{obs}}}$ is the total observed number of events along the direction of a dwarf, $\tau_{j}$ and ${A_{\mathrm{eff}}}_{,j}(E)$ are the observation time and effective area of each $j$ run, respectively, and $dN_{\gamma}/dE$ is the assumed source differential energy spectrum. The energy threshold $E_{\mathrm{min}}$ is defined here as the maximum of the efficiency curve which is defined as the effective area curve multiplied by the assumed source differential spectrum. In this case, the assumed differential spectrum is a power law of index -2.4. The bounded profile likelihood ratio statistical method of Rolke et al. Rolke et al. (2005) is used in this analysis to determine the upper limit on the number of gamma rays from the direction of each dwarf. The last column in Table 2 shows the resulting upper limits.

V.3 Combined search

Compared with examining individual dwarfs, pooling the data from all of them yields a search sensitive to weaker annihilation cross sections. The ON events from Boötes 1, Draco, Segue 1, and Ursa Minor are weighted according to Eq. (4) and summed according to Eq. (3). We do not include Willman 1 in the joint analyses because its irregular kinematics preclude a reliable determination of its $J$-profile via the Jeans equation (see discussion in Section III and Geringer-Sameth et al. (2015a)).

In this approach, the $J$-profiles must be taken into account since they are no longer degenerate with the cross section. We incorporate the systematic uncertainties in the dark matter distributions in the dwarfs by performing an ensemble of searches. For each, we assign each dwarf a $J$-profile from the posterior distribution of halo parameters Geringer-Sameth et al. (2015b). The scatter of the search resulting from many such realizations gives a measure of the systematic uncertainty due to our incomplete understanding of the density profiles in the dwarfs.

The results of the combined search are shown in Figure 4. The dashed lines bound 68% of the halo profile realizations and the solid line is the median significance. The combined observation shows no sign of dark matter annihilation in any channel.

V.4 Upper limits on the cross section

We slightly modify the procedure of Geringer-Sameth et al. (2015b) to compute cross section upper limits. In that work 95% confidence limits were generated using the Neyman construction of confidence belts. There, a hypothesis test is performed at every value of the cross section. The $\langle\sigma v\rangle$-space is divided into two regions where the hypothesis can and cannot be rejected at 95% confidence, with high enough values of $\langle\sigma v\rangle$ always being rejected. The boundary between the regions constitutes a 95% upper limit on the cross section. The hypothesis test is performed by asking, for a given value of $\langle\sigma v\rangle$, whether the probability that $T<T_{\mathrm{obs}}$ is less than 5%. If it is, then this value of the cross section is rejected.

In this work we adopt the ${\rm CL}_{s}$ technique Junk (1999); Read (2002) (sometimes called modified frequentist analysis) to produce upper limits. This method is strictly more conservative than the Neyman construction described above, i.e. always gives a larger upper limit, but has the benefit of being immune to downward fluctuations of background causing the upper limits to be much lower than the experimental sensitivity. That is, in the scheme described above, if there is a strong enough negative fluctuation of background so that $\mathrm{P}(T<T_{\mathrm{obs}}|\langle\sigma v\rangle=0)<5\%$ even the $\langle\sigma v\rangle=0$ hypothesis will be rejected causing the $\langle\sigma v\rangle$ upper limit to be zero.

The 95% confidence level upper limits on the annihilation cross section are presented in Figures  5 and  7. Each panel constrains dark matter with a 100% branching fraction into various Standard Model final states. The shaded band represents the $1\sigma$ systematic uncertainty induced by our imperfect knowledge of the dwarfs’ density profiles. They are produced by repeating the limit calculation over an ensemble of realizations of the dwarf halos from the distribution described in Section III. The lower, upper, and center of the band correspond to the 16th, 84th, and 50th percentiles of the distribution of limits over halo realizations. All other systematic uncertainties are negligible in this work in comparison and have been ignored.

As discussed in Section III, recent work has questioned the reliability of the $J$-profile of Segue 1 because of possible foreground contamination of its spectroscopic sample. By excluding Segue 1 from the combined analysis (i.e. setting its dark matter density to zero) we can bracket the effect that this unmodeled systematic uncertainty has on the particle physics constraints. Cross section limits are substantially weakened below a particle mass of about 400 GeV due to the lower energy threshold for the Segue 1 observations as compared to Draco and Ursa Minor (see Figure 2). Depending on the annihilation channel, excluding Segue 1 increases the $\langle\sigma v\rangle$ limit by a factor between 9-14 at 100 GeV, 4-7 at 200 GeV, 2-5 at 400 GeV, 2-3.3 at 1 TeV, and 1.2-2 above 10 TeV. Combined limits with and without Segue 1 included in Figures  5 and  7.

V.5 Statistical fluctuations

Hypothetically, if we were to repeat the measurement many times while holding the $J$-profiles of the dwarfs fixed, we would still obtain a distribution of limits due to statistical fluctuations intrinsic to a finite data set. We quantify the impact of the statistical uncertainty by looking at the distribution of the test statistic under the background-only hypothesis. That is, without using the events in the ON region, we take $T_{\mathrm{obs}}$ to be a given quantile of $\mathrm{P}(T\mid\langle\sigma v\rangle=0)$ and find the upper limit that would be obtained if this value had actually been measured. By taking the $0,\pm 1\sigma,\pm 2\sigma$ quantiles we find ranges where the observed limit is likely to lie. These are plotted in Figures  6 and  7. Specifically, due to random fluctuations of the background in the ON region, there is a 68% chance that the observed limit lies in the green band and a 95% chance that it lies in the yellow band. The dashed line is the median expected limit: there is a 50% chance that the observed limit is stronger than this. The solid black curve is the observed limit using the data from the ON region. This plot contains similar information to Figures  3 and 4. It shows how consistent the observations are with the background-only hypothesis. These plots were made using a particular set of $J$-profiles for the dwarfs, chosen to align well with Figures 5 and  7, and are meant to illustrate the experimental sensitivity of VERITAS and show the effect of background fluctuations on the cross section limits. The median limits for all channels are shown in Figure  8.

VI Conclusions and Discussion

The VERITAS limits in comparison with other concurrent gamma-ray instruments as well as older VERITAS results are shown in Figure 9. For the first time in an IACT DM search, this work uses the individual direction in addition to energy information of each event in the construction of the test statistic. The VERITAS results shown in this work are a substantial improvement over the entire WIMP mass range over the previous result with 48 hours on Segue 1 Aliu et al. (2015). VERITAS has a diverse dark matter program: observing time is divided between both the classical and ultrafain dSphs since we still have an imperfect knowledge of dwarf spheroidals and their J-profiles and their systematic uncertainties. This is especially important in light of the considerable uncertainty in the reconstruction of dwarf dark matter density profiles (see Section III and Figure 5). The strategy taken here of combining multiple targets in a single dark matter search mitigates sensitivity to future findings about particular galaxies. Pointed telescopes that rely heavily on a single target (e.g. Segue 1) may find their results susceptible to large, unaccounted systematic uncertainties. The Fermi-LAT, with a large duty cycle on all dSphs and low backgrounds, sets more stringent limits in the low mass range; however, the IACTs (VERITAS, MAGIC and HESS) put more stringent limits at the high mass range ($M\gtrsim$ 1 TeV), where Fermi-LAT has very low statistics.

Although no future hardware upgrades are currently planned for VERITAS, several advanced analysis techniques are starting to be deployed for VERITAS data. These techniques (e.g. boosted decision trees for $\gamma$/Hadron separationKrause et al. (2017)) could boost dark matter sensitivity by 30-50%. Additionally, the cuts used for this analysis were “point-like”, optimized for the detection of point sources. Nearly all the dark matter profiles for dwarf galaxies extend larger than the ON source region used in this work. An extended source analysis using a larger signal region could boost dark matter sensitivity by as much as a factor of two, dependent on the J-profile for each dSph. Dwarfs and other dark matter targets remain high-priority targets for the remainder of the lifetime of VERITAS.

The current upper limits on the annihilation cross section are about two orders of magnitude away from the relic abundance value ($\langle\sigma v\rangle\approx 10^{-26}{\mathrm{cm}}^{2}{\mathrm{s}}^{-1}$). This highlights the importance of improving both the instrumental sensitivity and the particle physics analysis. It is vital to extract all information present in the data to push experiments to the limit of their capability. The event weighting method, applied to IACT analysis for the first time, is a powerful and efficient way to combine multiple data sets and use our knowledge of the dark matter distribution and particle properties to perform optimal searches. For the first time, the event angular direction is used in addition to the energy of individual events for an IACT dark matter search.

It should be noted that the dark matter annilihation limits in this work were independently cross checked with a variation of the Full Likelihood utilized by the MAGIC collaborationAleksić et al. (2014) for a single halo realization for each dSph. The only major difference is that DM profiles were convolved with the VERITAS PSF described in this work, giving an integrated $J$-factor that is a function of energy. The combined dwarf limits of the two methods agreed within both the expected limits and $J$-factor systematic limits for the entire DM mass range used in this work.

To reach the thermal relic cross section, it may be necessary to combine all data taken from several gamma-ray telescopes into a single, deep search, expanding on the example that has been demonstrated by the MAGIC and Fermi-LAT collaborations MAGIC collaboration (2016). The methods we employed here may help prepare the experimental astroparticle physics community to accomplish this with upcoming experiments such as the Cherenkov Telescope Array (CTA) Acharya et al. (2013).

VII Acknowledgments

This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, and by NSERC in Canada. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument. SMK acknowledges support from the Department of Energy through Grant DE-SC0010010, and thanks the Aspen Center for Physics and the Center for Experimental and Theoretical Underground Physics (CETUP∗) for hospitality where part of this work was completed. This research used computational resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The VERITAS Collaboration is grateful to Trevor Weekes for his seminal contributions and leadership in the field of very high energy gamma-ray astrophysics, which made this study possible.