Boosting (In)direct Detection of Dark Matter

In this thesis, I study the expected direct and indirect detection signals of dark matter. More precisely, I study three aspects of dark matter; I use hydrodynamic simulations to extract properties of weakly interacting dark matter that are relevant for both direct and indirect detection signals, and construct viable dark matter models with interesting experimental signatures. First, I analyze the full scale Illustris simulation, and find that Galactic indirect detection signals are expected to be largely symmetric, while extragalactic signals are not, due to recent mergers and the presence of substructure. Second, through the study of the high resolution Milky Way simulation Eris, I find that metal-poor halo stars can be used as tracers for the dark matter velocity distribution. I use the Sloan Digital Sky Survey to obtain the first empirical velocity distribution of dark matter, which weakens the expected direct detection limits by up to an order of magnitude at masses GeV. Finally, I expand the weakly interacting dark matter paradigm by proposing a new dark matter model called boosted dark matter. This novel scenario contains a relativistic component with interesting hybrid direct and indirect detection signatures at neutrino experiments. I propose two search strategies for boosted dark matter, at Cherenkov-based experiments and future liquid-argon neutrino detectors.

Boosting (In)direct Detection of Dark Matter


Lina Necib

Submitted to the Department of Physics
in partial fulfillment of the requirements for the degree of

Doctor of Philosophy

at the


June 2017

© Massachusetts Institute of Technology 2017. All rights reserved.

Department of Physics
May 19, 2017

Certified by
Jesse Thaler
Associate Professor
Thesis Supervisor

Accepted by
Nergis Mavalvala
Associate Department Head of Physics

Boosting (In)direct Detection of Dark Matter


Lina Necib



Submitted to the Department of Physics

on May 19, 2017, in partial fulfillment of the

requirements for the degree of

Doctor of Philosophy


Thesis Supervisor: Jesse Thaler
Title: Associate Professor


I am very grateful to my amazing family: Sonia Bouzgarrou and Mohamed, Nada, Abdelkarim and Yasmine Necib for their unconditional support, even when they told their friends I had a "desk job at a university," and could not understand why I had not discovered dark matter yet. I am so thankful to my partner and best friend George Brova for his love and kindness, no matter how far apart we are, for teaching me how to code, and for just making every day better, easier, and happier.

I would like to thank my incredible experimentalist friends Sylvia Lewin, Alex Leder, Gabriel Collins, Nancy Aggrawal, and Amy Chell, for always being there for me, no matter what time of day or night it was, and for watching over me when life got too stressful. I am really thankful to my amazing fourth floor CTP friends Sarah Geller, Jasmine Brewer, Hongwan Liu, and Andrew Turner for making me look forward to coming to work every single day this past year, for being so kind, and making me feel like I belong here. I am also thankful to my CTP friends Nick Rodd, Gherardo Vita, Patrick Fitzpatrick, Kiaran Dave, Nikhil Raghuram, Dan Kolodrubetz, and Zachary Thomas, for making me feel part of a group, teaching me physics, and simply adding joy to the halls of the CTP. MIT is a wonderful place because of them, and I am saddened by the thought of leaving. I would also like to thank my amazing friend and housemate Alex Ji for being my unofficial astrophysics advisor, and helping me delve into this fascinating new area. I am more than lucky to have been surrounded by this amazing group of people.

While at MIT, I had the privilege of working with great collaborators. In particular, I would like to thank Kaustubh Agashe, Nicolas Bernal, Yanou Cui, Frederic Dreyer, Nayara Fonseca, Jonah Herzog-Arbeitman, Gordan Krnjaic, Piero Madau, Pedro Machado, Siddharth Mishra-Sharma, Jarrett Moon, Ben Safdi, Paul Torrey, Mark Vogelsberger, Taritree Wongjirad, and Jesus Zavala. I am also grateful to my undergraduate advisor Ken Lane for his endless support and kindness over the years.

I would also like to thank in particular the brilliant female physicists, Tracy Slatyer, Mariangela Lisanti, and Janet Conrad, with whom I had the great opportunity to work, and who have been exceptional role models in the early stages of my career.

Also, I will be always grateful to my officemate and great friend Ian Moult, for spending day and night in our office teaching me physics, getting me cookies when I needed them, keeping me afloat when everything was beyond stressful, and making our office the great place that it was.

Finally, I am certain that I would not have become the physicist I am today if it was not for my astounding advisor Jesse Thaler, whom I thank for his ubiquitous presence, the preponderance of hours he spent teaching me how to do research, how to give a presentation, how to write papers, and how to stay sane, for which I am eternally grateful.

List of Figures

Chapter 1 Introduction

1.1 Dark Matter

1.1.1 Introduction

One of the most important questions in physics today is Dark Matter (DM), a non-baryonic matter that represents a quarter of the energy budget of the universe, and is five times more abundant than the matter that makes up our visible universe. Understanding the origin of DM is therefore crucial; it is a question at the intersection of the smallest and largest scales, of particle physics and astrophysics. Making a consistent theory requires studying Nature across both the largest and smallest scales, which is the theme of this thesis.

Gravitational evidence for DM has been discovered when the movement of the luminous matter (gas, stars, galaxies) was inconsistent with the motion extrapolated from their brightness. Fritz Zwicky first measured the velocities of the stars of the Coma Cluster in 1933, and suggested that Dark Matter, a term coined by Henri Poincare in the early 20th century, might be more abundant than baryonic matter [Zwicky:1933gu]. DM became the spotlight of modern physics with the seminal work of Vera Rubin, in which she calculated the circular velocity of several galaxies, and found that the stellar velocity is independent of distance to the center of the host galaxies, a result that is unexpected if all the matter in a given galaxy only originated from the observable gas and stars [Rubin:1970zza, Rubin:1980zd, Rubin:1985ze]. The evolution of cosmology only helped solidify the existence of DM, especially with the extremely accurate measurements of the cosmic microwave background [Spergel:2003cb, Planck:2015xua, deSwart:2017heh].

Two plausible solutions to the discrepancy between the amount of luminous matter and the kinematics of the stars emerged: either a modification to gravity is required, or there exists an additional substance, which has yet to be detected, that makes up for the missing mass. What has driven the field to investigate further the theory of this new particle is the study of the Bullet Cluster [Clowe:2006xq]. A merger of the main Bullet Cluster with a subcluster has been observed, and the gas, stars, and the mass distribution have each been traced separately. Gas, representing 90 of the baryonic mass, has been detected in the middle of the merger, in collision between the two clusters. The total mass distribution of the mergers, traced through weak lensing, has been found to follow the collisionless stellar component, which should only account for of the total mass. It has thus been concluded that there exists a collisionless unseen matter that dominates the mass of the clusters, which in this case passed essentially unperturbed through the merging event.

Different search methods have been set up to detect DM: direct detection experiments, in which a DM particle scatters off a heavy nucleus causing it to recoil and produce a detectable signal, indirect detection experiments, where standard model (SM) particles that are the product of DM annihilation/decay are detected, and collider experiments that attempt to produce DM particles. In this thesis, I will focus on predicting the behavior of DM from simulations in both direct and indirect detection experiments, and then build a novel DM model with interesting hybrid (in)direct detection signals.

1.1.2 Simulations

body simulations are realizations of our known universe, or parts of it, tracing the evolution of DM “particles" from early times to the present day. These DM “particles" are generally solar masses, and they are not to be confused with the fundamental DM particle. Fundamental DM particles could include Weakly Interacting Dark Matter (WIMP) particles, or the DM candidate I introduced with collaborators, called Boosted Dark Matter (BDM) [Agashe:2014yua]; in both of these scenarios, the fundamental DM particle could have mass of several GeV (the mass of a few protons). These simulations have recently evolved into hydrodynamic simulations, which are body simulations that also include star and gas particles [1992ApJ...391..502K, 2005MNRAS.361..776S, Springel:2011yw, Dubois:2014lxa, 2015MNRAS.450.1349K]. Including baryonic physics furthers our understanding of the physics of DM, as the latter might be influenced by the kinematics of the baryons.

Figure 1-1: Evolution of the abundance of a WIMP DM of mass GeV, and an annihilation cross section to SM of cm/s.

1.1.3 WIMP Freeze-Out

One of the most commonly discussed DM models the WIMP scenario, which we therefore reference throughout this thesis. In this class of models, the weakly interacting DM particles are in thermal equilibrium with the SM at early times. For simplicity, let’s assume that the process that maintains this equilibrium is the DM annihilation , where denotes the DM particle. The evolution of the DM abundance as a function of inverse temperature of the universe is shown in Fig. 1-1. Initially, this process occurs along with SM SM, and therefore equilibrium is maintained between the DM and SM sectors. As the universe expands and cools down, the process (SM SM) is less and less efficient, and therefore the number density of DM is slowly depleted, as shown in Fig. 1-1 for , where is the DM mass. This would lead to the annihilation of all DM particles were it not for the process of freeze-out; as the universe expands, it is less and less likely for DM particles to find one another and annihilate, and therefore the abundance of DM (per co-moving volume) freezes out around into the density that we have today.

1.2 Spherical Cows of Dark Matter Indirect Detection

In hydrodynamic simulations, DM particles interact only gravitationally with baryons. Although these simulations have large particles, which instead of being thought of as one object like a star, can be thought of as a stellar/DM population, they can consistently reproduce large scale features of the universe. Some zoom-in simulations, which are realizations of the Milky-Way, can successfully simulate the spiral arms of the Milky Way as well as the independence of the velocity from the distance to the center of the galaxy [Guedes:2011ux]. This leads us to believe that some of the features of DM are well reproduced in these simulations.

In Chapter 2, I present work performed with Nicolás Bernal and Tracy Slatyer, where I performed the analyses on the simulations [Bernal:2016guq]. I use the cosmological simulation Illustris [Vogelsberger:2014dza, Genel:2014lma, Nelson:2015dga] to predict the morphology of DM indirect detection signals. Galaxies reside in the centers of large DM clumps referred to as DM halos. Previous studies of DM halos have shown that they tend to be triaxially distributed [1991ApJ...377..365K, Dubinski:1993df, Debattista:2007yz, Pedrosa:2009rw, Tissera:2009cm, 2011ApJ...740...25B, Zemp:2011nk, Bryan:2012mw]. Indirect detection signals, however, are sensitive to the DM morphology in projection, perpendicular to the line of sight. In Chapter 2, I perform two analyses of the expected projected morphology of DM from simulations: a Galactic analysis, in which I place an observer at 8.5 kpc from the halo center for halos of comparable size to the Milky Way, and an extragalactic analysis, in which I place the observer outside the halo. I found that the expected Galactic indirect detection signals are highly symmetric, while the extragalactic signals, especially those of annihilating DM are almost uniformly distributed across axis ratios.

In these analyses, I built a new metric for the moment of inertia tensor that is weighed in DM “luminosity". This metric is easily applicable to data. Indeed, I compared the simulated axis ratios with those of the Milky Way in Gamma rays detected by the Fermi-LAT telescope [Atwood:2009ez]. This gamma ray data is of particular interest since, a few years ago, an excess of gamma rays has been detected at the center of the Milky Way [Goodenough:2009gk, Daylan:2014rsa, Calore:2014xka, Zhou:2014lva, Linden:2016rcf], and many have theorized about the possibility of a DM origin for this excess [Modak:2013jya, Abdullah:2014lla, Berlin:2014pya, Boehm:2014bia, Martin:2014sxa, Fonseca:2015rwa, TheFermi-LAT:2015kwa]. I therefore compared the full gamma ray Fermi sky map in energies GeV, as well as the residual gamma rays attributed to the excess with the predicted Illustris result. I found that the excess is in agreement with the expected morphology of a DM signal.

A stacked analysis of clusters in X-rays has shown the presence of a line at 3.5 keV, unaccounted for by atomic emission lines [Bulbul:2014sua, Bulbul:2014ala, Franse:2016dln]. A few studies attributed these X-rays to DM being a sterile neutrino [Ruchayskiy:2015onc, Iakubovskyi:2015dna]. Although it is difficult to isolate the photons contributing to the lines, I studied the full X-ray data from clusters and found that they tend to be more symmetric than we expected the background to be from the active merger history of Illustris. This should be followed closely, and will be of particular interest with the improvement in spatial resolutions of future telescopes.

1.3 Empirical Determination of Dark Matter Velocity Using Metal Poor Stars

Another method of detecting DM is direct detection. The detection rate of DM scattering off heavy nuclei depends on astrophysical quantities, and in particular the velocity distribution of DM. Most DM direct detection limits use the Standard Halo Model (SHM) which assumes a Maxwellian DM velocity distribution, obtained from a collisionless isothermal density distribution with a flat rotation curve. Using a zoom-in simulation called Eris [Guedes:2011ux, Pillepich:2014jfa], I found that metal poor stars, as they originate in satellite galaxies that merge into the Milky Way halo similarly to DM, are excellent tracers for the kinematics of DM.

I used the velocity distribution of the stellar halo which is formed by old metal-poor stars kinematically unassociated with the stellar disk, obtained from the Sloan Digital Survey (SDSS) [Juric:2005zr, Ivezic:2008wk, Bond:2009mh], to infer that of the DM, as shown in Chapter 3, which is a published work performed in collaboration with Jonah Herzog-Arbeitman, Mariangela Lisanti, and Piero Madau [eris_paper]. This is important in setting direct detection limits. I found that at low masses, the limits set using the SHM are an order of magnitude stronger than those found using the SDSS distribution. This work can be extended to more local stellar catalogs such as RAVE-TGAS [2006AJ....132.1645S, 2017AJ....153...75K] to study the local DM velocity distribution, as well as DM substructure [ravepaper].

1.4 Boosted Dark Matter

All these previous analyses target collisionless DM (or the WIMP paradigm in the limit where the DM interacts weakly), but it is important to investigate other DM scenarios. One such scenario is that of BDM, developed in Chapter 4, which presents published work performed in collaboration with Kaustubh Agashe, Yanou Cui, and Jesse Thaler [Agashe:2014yua]. It consists of two DM components and , where is heavier and dominates the DM abundance. annihilates to which, due to the mass hierarchy, gets a Lorentz boost today. Such boosted particles can scatter off electrons in neutrino experiments, and cause an excess of electron events over muons, which can be studied in experiments such as Super-Kamiokande (Super-K) and its future upgrades [Fukuda:2002uc, Galkin:2008qe, Kearns:2013lea]. We are therefore directly detecting the component, while indirectly detecting the DM particle.

The boosted particles generally emit electrons in the forward direction, which therefore point back to the DM origin —regions dense in DM, such as the Galactic Center, or dwarf galaxies [Necib:2016aez]. This is useful in Cherenkov-based experiments like Super-K in limiting the backgrounds, but can be used to set more constraining limits on BDM from DM point sources such as dwarf galaxies at liquid-argon detectors like the Deep Underground Neutrino Experiment (DUNE) [Acciarri:2015uup]. Argon-based detectors have excellent particle identification and spatial resolution, and can therefore extend the reach of observable BDM signals, as shown in Chapter 5, which is published work performed in collaboration with Jarrett Moon, Taritree Wongjirad, and Janet Conrad [Necib:2016aez]. I present the limits on BDM in a simplified constant amplitude model across multiple experimental technologies, and propose a search strategy that the limits atmospheric backgrounds for each technology separately.

1.5 Additional Topics

During my time at MIT, I worked on multiple aspects of dark matter research, beyond the ones I described above. In particular, I studied an asymmetric DM model with interesting indirect detection signals in gamma rays, analyzed Fermi-LAT gamma ray data at high latitudes to extract the fraction of the extragalactic gamma ray background (EGB) dominated by point sources, studied the effect of an ultralight scalar field on neutrino oscillation parameters, and finally developed new jet substructure observables to discriminate boosted and bosons as well as top quarks from backgrounds.

In Ref. [Fonseca:2015rwa], my collaborators and I built an asymmetric DM model with a hybrid thermal history. This particular DM model is both asymmetric and WIMP-like, and thus the final DM abundance is set by a hybrid mechanism at the intersection of thermal freeze-out and asymmetric DM. Such model has interesting potential indirect detection signals, which is generally lacking in asymmetric DM models. In particular, asymmetric DM candidates studied produced gamma rays consistent with the GeV excess at the center of the galaxy [Daylan:2014rsa].

I have also studied gamma rays from the Fermi-LAT telescope in order to estimate the point source contribution to the extragalactic gamma ray background (EGB) [Lisanti:2016jub]. This is of particular interest since an excess of high energy neutrinos has recently been detected, the origin of which is still unknown [Aartsen:2013jdh]. If these high energy neutrinos originated from star forming galaxies, these same galaxies would have also emitted gamma rays that would contribute to the EGB. In our work, we used a novel analysis method called Non Poissonian Template Fit [Malyshev:2011zi, Mishra-Sharma:2016gis] to estimate the number of point sources in the EGB. This sets a limit on the fraction of diffuse gamma rays emitted by star forming galaxies. We alleviated some of the tension found in previous work [Zechlin:2015wdz], which set stronger bounds on the diffuse emission from the star forming galaxies and were therefore inconsistent with the interpretation that these high energy neutrinos originated in star forming galaxies.

With collaborators, I considered sub-eV scalar DM coupling to neutrinos as , inducing temporal variations on neutrino parameters [Krnjaic:2017zlz]. The scalar has a local field value at the spacetime coordinate (, ), and that can be written as


where is a quantity set by the coupling of to neutrinos, the mass of , and the local DM density. denotes the virial velocity of DM. The oscillations of the scalar field, although faster than , can introduce distorted neutrino oscillations by shifting the square mass differences as well as the neutrino angles by a quantity proportional to . From precise measurements of neutrino parameters, we can set bounds on the mass and coupling of the scalar , over decades in masses, from eV to eV, expecting better improvement by the future experiments DUNE [Acciarri:2015uup] and JUNO [Li:2014qca].

Finally, another way of detecting DM is to produce it at colliders. Colliders are overwhelmingly dominated by Quantum Chromodynamics (QCD) backgrounds. Sprays of QCD particles are collimated into jets, some of which might have an inner substructure that helps identify its origin. Multiple jet substructure techniques have been developed to study multi-prong objects [Abdesselam:2010pt, Altheimer:2012mn, Altheimer:2013yza, Adams:2015hiv, Larkoski:2015uaa]. Some of these are based on energy correlation functions, which are functions of energy fractions and angles of the different components of the jet [Larkoski:2013eya]. In order to disentangle new physics from QCD backgrounds, my collaborators and I introduced new jet substructure observables that discriminate 2 and 3 prong jets from QCD backgrounds [Moult:2016cvt]. We used power counting techniques [Larkoski:2014gra, Larkoski:2015kga] to develop a new set of observables, , , and that show excellent discrimination power, but also stability under changes in mass and cuts. Such a property is desirable in highly desirable in experimental searches [Dolen:2016kst]. These observables are currently being investigated by the CMS collaboration of the Large Hadron Collider, with a recent analysis published just days ago [CMS:2017dhi].

Chapter 2 Spherical Cows of Dark Matter

2.1 Introduction

Gravitational evidence for dark matter (DM) is well established [Zwicky:1933gu, Rubin:1970zza, Markevitch:2003at], yet DM still evades all other means of detection [Akerib:2015rjg, Agnese:2015nto, Ackermann:2015zua, TheFermi-LAT:2015kwa, Khachatryan:2014rra]. A current focus of the search for DM is indirect detection: DM annihilation or decay could produce observable Standard Model (SM) pchapters, including photons. If such a signal was detected, the direction of the incoming photons could be used to map out the morphology of DM halos.

N-body simulations constitute a valuable tool for studying the expected DM distribution [1961AJ.....66..590V, 1963MNRAS.126..223A, 1992ApJ...391..502K, 1983ApJ...270..390M, Bryan:2012mw], and can be used to predict the properties of indirect signals from DM [Stoehr:2003hf, Bernal:2014mmt, Calore:2015oya]. Hydrodynamic simulations include baryonic matter as well as DM, and thus can probe the impact of baryonic feedback on the DM distribution [Bryan:2012mw]. With recent hydrodynamic simulations that generate large ensembles of DM halos, we can make statistical statements about the general properties of DM halos with and without baryons [Dubinski:1991bm, Warren:1992tr, Dubinski:1993df, Jing:1994sg, Jing:2002np, Kazantzidis:2004vu, Allgood:2005eu, Debattista:2007yz, Schneider:2011ed, Bruderer:2015rnh]. In particular, as we demonstrate in this work, we can map out the full distribution of properties relevant to indirect DM searches, rather than relying on a small number of example halos.

In this chapter, we focus on studying the morphology of indirect detection signals using N-body simulations. More specifically, we study sphericity/asymmetry of signals after projection along the line of sight. We perform a statistical analysis of the annihilation/decay signatures of a large number of halos in two simulations: the hydrodynamic simulation Illustris-1, which includes DM and baryons, and its DM-only equivalent Illustris-1-Dark [Vogelsberger:2014dza, Vogelsberger:2014kha]. We predict the shape of annihilation/decay DM signals from Galactic and extragalactic (EG) sources, and diagnose the effect of baryons on the asymmetry and sphericity of these signals. For the remainder of the text, we will refer to signals as “spherical” if they could be produced by the line-of-sight projection of a spherical 3D source of photons; i.e. they are symmetric under rotation of the sky around the line-of-sight pointing toward their center.

Several potential signals have appeared in indirect DM searches over the past few years. An anomalous emission line at 3.5 keV has been found in a stacked analysis of 73 galaxy clusters [Bulbul:2014sua] and in other regions [Ruchayskiy:2015onc, Iakubovskyi:2015dna, Bulbul:2014ala, Franse:2016dln]. Analysis of data from the Fermi Gamma-Ray Telescope (hereafter Fermi) [Atwood:2009ez] has shown an unexplained spherically symmetric excess of (GeV) gamma rays at the center of the Galaxy [Goodenough:2009gk, Abazajian:2012pn, Daylan:2014rsa, Zhou:2014lva, Calore:2014xka, TheFermi-LAT:2015kwa, Linden:2016rcf]. Studying expected properties could help discriminate DM against astrophysical backgrounds.

This chapter is organized as follows. First we introduce our methodology in Sec. 2.2; we describe the Illustris simulation, and the related computations of DM density, and define the metrics used for the determination of halo shapes. We then perform two analyses of annihilation and decay signals, one where the observer is situated at a location 8.5 kpc from the center of the halo (Sec. 2.3), and one where the observer is outside the halo (Sec. 2.4). In each of these sections, we present the overall distributions for asymmetry and axis ratio. For the former (Galactic) analysis, we focus on the subcategory of Milky-Way (MW) type halos. For the latter (extragalactic) analysis, we focus on cluster-sized halos. In both cases, we also study possible correlations between halo axes and the baryonic disk. In Sec. 2.5, we discuss two case studies of the morphology of astrophysical backgrounds for DM searches, first considering the gamma-ray background and signal for the Fermi inner galaxy excess, and then the clusters in which the keV line is detected. We summarize and conclude in Sec. 2.6.

2.2 Methodology

In the context of indirect DM searches for photons or neutrinos, the quantity of interest for decay (annihilation) is the integrated DM density (density squared) of DM particles along the line of sight. This is referred to as the -factor; to compute it within the Illustris simulation, we must define it in the context of the discrete representation of the underlying matter distribution.

2.2.1 Illustris Simulation

The Illustris simulation is a publicly available111 hydrodynamic simulation that traces the evolution of DM particles, as well as gas, stars and black holes across redshifts from to today [Vogelsberger:2013eka, Torrey:2013pwa, Vogelsberger:2014dza, Vogelsberger:2014kha, Genel:2014lma]. The Illustris simulation employs a comprehensive suite of baryon physics including stellar evolution and feedback, gas recycling, supermassive black hole growth, and feedback from active galactic nuclei [Vogelsberger:2014dza]. In this work, we focus on the last snapshot at , which reflects the simulated state of today’s Universe [Nelson:2015dga]. The simulation is conducted at 3 different resolution levels, Illustris-1, Illustris-2, and Illustris-3. It also includes the same set of simulations for DM only particles, at the same resolution levels, Illustris-1-Dark, Illustris-2-Dark, and Illustris-3-Dark. The simulations cover a total volume of . In this work, we focus on the highest resolution simulations Illustris-1 and Illustris-1-Dark. Parameters of the simulations are shown in Table 2.1, including the mass of the DM and baryon particles, and the spatial resolutions of the simulations. The mass of a DM particle is fixed throughout the simulation, but that of a baryonic particle (which sums the mass of the gas, stars and black holes) is not conserved, but kept within a factor of 2 of the quoted baryonic mass . Gravity is included with softening of the potential at small scales to avoid numerical two-body particle scattering [Springel:2009aa]. The softening lengths for both DM and baryons, and , are also shown in Table 2.1.

Illustris-1 1.4 0.7
Illustris-1-Dark 1.4
Table 2.1: The particle masses and softening lengths for the Illustris-1 and Illustris-1-Dark simulations [illustris]; “DM” subscripts label DM, while “b” subscripts label baryons.

The Illustris simulation used the friends-of-friends (FOF) algorithm to identify DM halos [Davis:1985rj]. The Illustris-1 simulation has 7713601 halos, and 4366546 identified as subhalos. The Illustris-1-Dark simulation includes 4263625 halos and 4872374 subhalos. To limit the impact of poorly resolved objects on our results, we only examine halos with at least 1000 DM particles. This cut leaves and halos for Illustris-1 and Illustris-1-Dark respectively. The halo mass function of the Illustris simulation described in Ref. [Vogelsberger:2014dza] is in good agreement with the empirical data. Deviations from observations are present at the low and high end of the resolved masses, where the details of the implementation of the stellar and AGN feedback are important. The mass range of halos that pass the 1000 DM particle cut is to .

2.2.2 Computing -factors

The quantity of interest in this analysis is the -factor, defined by


where is the density of DM, and the integral is along the line of sight and over solid angles. In order to integrate the local DM density (or density squared) given the discrete particle distribution, we use a kernel summation interpolant to reconstruct a continuous DM density field [1985A&A...149..135M, Hernquist:1988zk, Springel:2011yw]. More explicitly, for a field , one can define a smoothed interpolated version , related to through a kernel function


where defines the length of the smoothing. The kernel function approaches a delta function as . We use a cubic spline to compute the local density :


The smoothing kernel in this case is , where is chosen to be the distance to the 33 nearest neighbor of the point considered. In order to compute the density along a particular line of sight defined by the galactic coordinates , we sum over the density of the 32 nearest neighbors to a particular point, and adjust the next step in the integral to be the newly found . We have checked that doubling the number of neighbor particles employed in this procedure from 32 to 64 does not alter our results. This is due to the fact that the contribution of further particles is proportional to .

We proceed to construct sky maps of annihilation/decay -factors for each halo, by placing an observer at kpc from the center of the halo along the -axis of the simulation, and compute the -factor for different values of galactic coordinates . The center of the halo is defined as the location of the gravitational potential minimum. We use the package HEALPix222 to divide the sky into equal area pixels [2005ApJ...622..759G]. The total number of pixels in a map is defined by


where nside is an input parameter that defines the pixelation.

Figure 2-1: Logarithmic map for DM annihilation for halo labeled “0” in the Illustris-1 simulation. Lighter colors signify higher DM luminosity. We used HEALPix with nside (see Eq. 2.4). This map is taken by positioning the observer at 8.5 kpc from the center of the halo. The mass of this halo is .

We show an example of such constructed maps in Fig. 2-1, for the annihilation signal of a particular halo. For the example chosen, the halo radius, defined as the radius such that the average density interior to that radius is 200 times the critical density of the universe,333The distances in the simulation are presented in units of kpc/, where is the reduced Hubble constant so that km/sec/Mpc. The cosmological parameters used are , and [Vogelsberger:2014dza]. is 1659 kpc kpc, much larger than the observer distance from the center of the halo kpc, and therefore, there is still sizable signal from high latitudes.

2.2.3 Asymmetry Parameterization

Axis Ratio

One measure of sphericity is the axis ratio. First, we summarize previous methods of finding the axis ratio and the major axis. Unlike previous analyses that analyzed the halo shapes through the 3D moment of inertia tensor [1986ApJ...304...15B, 1987ApJ...319..575B, Dubinski:1991bm, Jing:2002np, Allgood:2005eu, Schneider:2011ed], we compute the 2D projection of the inertia tensor along the plane perpendicular to the line between the observer and the halo center, as all indirect detection signals are found in projection. That is defined as [Tenneti:2014poa]


where the sum is taken over the DM particles of the halo, and correspond to the coordinates of the particle projected on the plane perpendicular to the observer.444See Ref. [2012MNRAS.420.3303B] for a discussion of the different definitions of the inertia tensor that occur in the literature. For example, if the observer is located along the axis (where the center of the Cartesian coordinate system is at the center of the halo, defined in the simulation as the location of the most-bound particle), run over the four combinations of the and coordinates for each of the DM particles of the halo. The axis ratio is defined as the ratio of the square root of the eigenvalues of the inertia tensor, where in this work we use the convention where the axis ratio is always less than 1. The major axis of the halo is the eigenvector corresponding to the largest eigenvalue. In this notation, the axis ratio of a spherical halo is 1.

We introduce a variation on the inertia tensor defined in Eq. 2.5 that is adaptable to indirect detection signals. This new inertia tensor uses the same information that we would have looking at a DM annihilation/decay signal. In this case, the DM particle coordinates are weighed by luminosity in DM signal, which is the -factor at that location. The new inertia tensor that we call the -tensor is therefore


where the coordinates are obtained from scanning through the pixels in the sky and inferring the Cartesian coordinates of that particular pixel (assuming we live in a sphere), and is given by Eq. 2.1 at a point in the sky given by the coordinates . With this approach, all particles within the same line of sight contribute only once but their contribution is weighed with the observed intensity of the signal. As above, the observed axis ratio is defined as the ratio of the square roots of the eigenvalues of the -tensor, and the halo’s major axis is the eigenvector with the largest eigenvalue.

Quadrant Analysis

As a second parameterization of the observed asymmetry of DM signals, we divide the observed sky into four equal quadrants, with the origin of the coordinate system lying along the line of sight to the center of the halo. The halos are oriented randomly relative to the quadrant boundaries unless otherwise stated.

We then determine the -factor associated with each quadrant as discussed in Sec. 2.2.2


where labels the different quadrants, and is the list of pixels in quadrant . is the value of the -factor found at pixel . The quadrants are labeled such that quadrant 1 is adjacent to 2 and 4, and opposite to 3. We define the following ratios, describing the relative predicted emission in pairs of opposite or adjacent quadrants


For signals that appear spherical from the point of view of the observer, the ratios defined in Eqs. 2.8 and 2.9 . For signals that appear strongly elongated or asymmetric, or depending on which quadrants dominate the DM signal.

2.3 Galactic Analysis

Figure 2-2: Histogram of the conventional 2-dimensional axis ratio (left) and the newly defined observed axis ratio (right) for annihilation and decay, comparing both Illustris-1 and Illustris-1-Dark (see Eq. 2.6).

In this section, we perform a statistical analysis of the sphericity of annihilation/decay signals as observed from a location kpc from the halo center, similar to the Earth’s separation from the center of the Milky Way [2009A&A...498...95V, Gillessen:2008qv]. We first study the distribution of the observed axis ratio in annihilation and decay as defined in Sec. 2.2.3, as well as the distribution of the ratio in intensity of opposite and adjacent quadrants as introduced in Sec. 2.2.3. We plot the histograms of probability densities in each distribution, counting each halo just once unless otherwise stated. We then focus our analysis on Milky-Way-like halos, orienting the observer to be on the halo disk, and study the axis ratio distribution. We finally examine the possibility of correlations between the halo minor axis and the baryonic disk.

2.3.1 Observed Axis Ratio

For comparison, we first illustrate the distribution of the axis ratio as obtained from the two-dimensional inertia tensor defined in Eq. 2.5. As shown in Fig. 2-2 (left), we find that the distribution peaks at axis ratio , which is consistent with results found from the projected shapes of DM halos inferred from the position of galaxies in the Sloan Digital Sky Survey Data Release 4 [Wang:2007ud]. This suggests that halos are mostly symmetric in Cartesian projection, which is independent of the observer’s distance. In the same figure, we show the distributions of axis ratio in the DM-only and the DM+baryon simulations. We only find a minor tendency for halos in the DM+baryons simulation to be more symmetric in projection.

Figure 2-3: The distribution of the asymmetry parameters (left) and (right) for DM annihilation as observed from a point 8.5 kpc from the halo center, in halos taken from the DM-only and DM+baryon simulations of Illustris-1. -factors are computed over all halos (blue) as well as when omitting the inner disk (red), and through the inner only (green). All regions are centered on the halo center. The inset shows a zoom-in of the region of small /.

Moving to the -tensor, we plot the distributions of the observed axis ratio for indirect detection signals in Fig. 2-2 (right), for observers at distances from the center of the halo comparable to the solar circle radius.555Since in this figure we include all halos, not only MW-sized halos, the solar circle radius may be much smaller or larger as a fraction of the virial radius than it is in the Milky Way. However, we have checked that for all but a handful of halos, the observer is still within the halo virial radius at this distance; furthermore, we show results specifically for MW-sized halos in Sec. 2.3.4. These distributions are generally more peaked towards higher axis ratios, with the peak now at axis ratio . The distributions for annihilation signals are broader than those for decay; annihilation signals accentuate anisotropies due to their dependence on the square of the DM density. In the case of annihilating DM, the effect of including baryons is more pronounced; DM-only halos are less spherical/more elongated than those including baryons. This effect has been previously studied and our results are consistent with previous analyses [1991ApJ...377..365K, Dubinski:1993df, Debattista:2007yz, Pedrosa:2009rw, Tissera:2009cm, Zemp:2011nk, 2011ApJ...740...25B, Bryan:2012mw, Kelso:2016qqj].

2.3.2 Quadrant Analysis

In this section, we compute the distributions of and as defined in Eqs. 2.8 and 2.9, for the halos of Illustris-1 and Illustris-1-Dark. This measure is of particular interest as it is easy to test in template analyses (See for example Refs. [Finkbeiner2004, Daylan:2014rsa]). For every halo, we build a map of the -factors at different points in the sky in annihilation and decay following the procedure in Sec. 2.2.2. We picked nside (see Eq. 2.4), but tested that the results are stable under a change of nside. In order to characterize the spatial distribution of indirect detection signals, we consider three different regions of each halo. First, we take the halos as a whole, then we omit the inner cone of half angle and finally we look at the inner cone of half angle . 10 degrees from the perspective of an observer at a distance of 8.5 kpc covers one softening length kpc . We show the results for DM annihilation in Fig. 2-3 as a probability distribution of and . Similar distributions can be found for decay as shown in App. A.1. We have included Gaussian error bars, shown as the shaded regions of Fig. 2-3.

For the case of annihilation, shown in Fig. 2-3, we find that there is a notable difference in the sphericity of the DM-only simulation compared to the DM+baryons simulation; there are more asymmetric halos ()666By we mean and . in the DM-only case. This result confirms the finding in the previous section, that inclusion of baryons tend to make DM halos more spherical, although in both cases, most halos have mostly spherical decay/annihilation signals.

In order to understand the effect of the inner 10 degrees, we compare the histogram of the whole halo minus the inner 10 degrees with that of the distribution that includes the entire halo. We find that deviations only happen at larger values of and which occur with probabilities less than a few percent. This leads us to believe that our distributions are largely unaffected by the inner few softening lengths where resolution artefacts might play a larger role. The innermost region of the halo tends to also be more spherical than outer regions, as shown in Fig. 2-3 when comparing the distribution in which we omitted the inner with the distribution across the whole halo. The distributions of are slightly more peaked towards zero in the inner cone of half angle 30 degrees. The differences between the three regions intensify in the tail of the distribution. We note that due to the resolution of the simulation, we cannot make a statement on the sphericity of the inner few degrees around the Galactic Center.

2.3.3 Correlation with Baryon Disk

In this section, we examine the distribution of the angle between the angular momentum vector of the baryonic disk and halo’s minor axis found in projection. We first outline how to compute each of these axes.

In order to find the orientation of the baryon disk, we compute the three dimensional angular momentum vector of the star forming gas (SG). We first determine the location of the gas particles with a positive star forming rate, since the gas must have cooled to form stars and contribute to the angular momentum of the disk, and we then compute the 3D angular momentum vector for a particular halo as


where , with () the coordinates of the particle (the center of the halo), and () is the mass (3D velocity) of the gas cell. We then project the angular momentum vector on the plane perpendicular to the line between the observer and the halo, and label the new 2D angular momentum vector .

We now turn to computing the halo’s minor axis. As shown in Sec. 2.2.3, there are multiple ways to compute the inertia tensor in projection: (1) projecting the particle coordinates onto the plane perpendicular to the line between the observer and the center of the halo, (2) computing the -tensor in annihilation and (3) computing the -tensor in decay. In each of these cases, we compute the eigenvalues and eigenvectors of the tensor. The eigenvector that corresponds to the smallest eigenvalue is taken to be the halo’s minor axis .

In order to compute the angle between the minor halo axis and the angular momentum vector, we consider the normalized inner product . In Fig. 2-4, we compare the newly found distribution of the angle to that of a flat distribution in using the three definitions of the minor axis. We find that these three measures are consistent with a slight correlation between the halo’s minor axis (found in projection) and the angular momentum vector; there is a slight preference for the two axes to be aligned to each other, i.e. . This implies a slight preference for the halo’s major axis to be aligned with the baryonic disk.

Figure 2-4: Histogram of the angle between the minor axis of the DM halo and the angular momentum vector of the star forming gas. The minor axis is the eigenvector corresponding to the smallest eigenvalue of the inertia tensor found from regular projection (Eq. 2.5) or the weighed inertia tensor (Eq. 2.6) in the case of annihilation and decay. is the mean of the distribution (a flat distribution in ).

Our results found in projection are consistent with previous 3D analyses [Bailin:2005xq, AragonCalvo:2006ay, Hahn:2010ma, Debattista:2013rm, Velliscig:2015ffa, 2012MNRAS.420.3303B, Codis:2014awa, Dubois:2014lxa, Chisari:2015qga, Velliscig:2015ixa, Kiessling:2015sma, Debattista:2015hia, 2011MNRAS.413.1973W] (Ref. [Tenneti:2014poa] shows a 2D projected misalignment angle of order ).

2.3.4 Milky-Way-like Halos

We now focus on the subset of MW-like halos, to see if they share consistent sphericity properties with the overall sample. This is crucial, as were we to discover DM through its annihilation/decay to SM particles in the MW, the signal/background could be analyzed exactly in the same way we analyze the Illustris data. To that aim, we require the following:

  • Total mass: The total mass of the halo lies in the range (see for example Ref. [Schaller:2015mua])


    where is the mass of the halo enclosed in a sphere with a mean density 200 times the critical density of the Universe today. is the solar mass. The number of halos in Illustris-1 within this mass range is 1652.

  • Stellar mass: The total stellar mass lies within the range [McMillan:2011wd, Calore:2015oya]


    This further drops the number of MW-like halos in the Illustris-1 simulation to 650.

We then perform the analysis of Sec. 2.3.1 on this restricted sample of halos. We find that indeed the distributions shown in Fig. 2-5 are consistent with the more general results shown in Fig. 2-2, though with lower statistics. The DM signal is expected to be spherical, and peaks at values , although with a more peaked distribution.

Figure 2-5: Histogram of the observed axis ratio for annihilation and decay of MW-like halos as defined by the requirements in Eqs. 2.11 and 2.12. The distribution of the MW-like halos is shown in solid lines while the distribution of rotated halos to increase statistics in shown in dotted lines.

In order to increase statistics, we study the observed axis ratio for MW-like halos from 12 different projections by placing the observer at different locations along the sphere centered at the halo center with radius  kpc. We plot the new distribution in Fig. 2-5 in dotted lines. We find that the distributions are preserved but more smoothed out.

We will discuss in Sec. 2.5.1 the axis ratio of the gamma-ray sky as observed by Fermi [Atwood:2009ez], as compared to the distribution of observed axis ratios in MW-like DM halos.

2.4 Extragalactic Analysis

In this section, we perform a similar analysis to Sec. 2.3, but we now situate the observer outside the halo in consideration. As an example, we set the observer at a distance


where is the distance from the halo center at which the overdensity of the halo is 200 times the critical density of the Universe. We check that our results are independent of the distance between the observer and the center of the halo as long as . We increase nside to 512 in this analysis to be able to resolve smaller structures of the halos (See Eq. 2.4), then downgrade the maps to nside = 32 for computational efficiency in the analysis.777If the maps are generated originally at nside=32, the lines of sight through the center of each pixel do not adequately describe the average emission from that pixel, as large variations in the brightness can occur on scales smaller than a pixel. Consequently, changes in the pixelation can markedly change the results. To resolve this problem, we generate the maps at higher resolution, and use these higher-resolution maps to determine the total emission in each (nside=32) pixel. Once this is done, our results are stable with respect to the choice of pixel size. With this choice of , the halos cover of the map, which is higher than most extragalactic signals, but we do so in order to resolve the inner structure.

Figure 2-6: Histogram of the observed axis ratio for annihilation and decay for extragalactic sources, comparing both Illustris-1 and Illustris-1-Dark (see Eq. 2.6).
Figure 2-7: The distribution of the asymmetry parameters (left) and (right) for DM annihilation as observed from a point outside the halo, taken from DM-only and DM+baryons simulations. -factors are computed over all halos (blue) and through the inner only (green).

2.4.1 Axis Ratio

As with our previous analysis (Sec. 2.3), we study the distribution of both the relevant axis ratios and the quadrant parameters. In Fig. 2-6, we plot the distribution of observed axis ratio in the case of annihilation and decay of DM, for the DM-only simulation as well as DM+baryons simulation. In both decay and annihilation, the distributions of axis ratio are flatter than in the Galactic analysis; while decay signals still generally have fairly spherical profiles, the distribution of axis ratio for annihilation signals is nearly flat, although slightly peaked around 0.9. An interesting feature is the non-negligible fraction of halos with axis ratio . As we will explore in Sec. 2.4.4, this behavior is due to halo mergers.

Serving both as a consistency check and as a study of the baryonic effects, the DM-only simulation exhibits similar features to the baryonic simulation, with the distributions shifted slightly towards lower values of the axis ratio.

2.4.2 Quadrant Analysis

As shown in Fig. 2-7, the ratios of opposite and adjacent quadrants show that DM signals are less spherical when observed at a larger distance. This is reasonable as all features of the halo are at an equivalent distance from the observer, while in the Galactic analysis, it is harder to resolve small anisotropies that are at a larger distance from the observer. These results suggest that especially for extragalactic annihilation signals, observation of an elongated morphology could not be used to disfavor a DM hypothesis, and there is no reason to expect highly spherical signals that could easily be distinguished from astrophysical sources with complex and non-spherical distributions. (However, if the primary astrophysical backgrounds were near-spherical, a highly elongated profile might provide a hint for a DM origin.)

In order to omit possible signals from secondary subhalos which are off the center of the halo, defined by the most bound particle, we analyze the ratios of the quadrants within a cone of half angle . We find that the distributions within the cone do indeed appear more spherical, but the effect generally dominates at the tail of the distribution, where the asphericity is more extreme.

Figure 2-8: Histogram of the angle between the minor axis of the DM halo for the case of extragalactic signals and the angular momentum vector of the star forming gas. The minor axis is the eigenvector corresponding to the lowest eigenvalue of the inertia tensor found from regular projection (Eq. 2.5) or the weighed inertia tensor (Eq. 2.6) in the case of annihilation and decay. is the mean of the distribution (a flat distribution in ).

2.4.3 Correlation with Baryon Disk

Similarly to Sec. 2.3.3, we plot in Fig. 2-8 the distribution of the angle between the halo’s minor axis (found in annihilation and decay of the DM particles) with the angular momentum vector of the star forming gas, this time analyzing the minor axis from the extragalactic maps. We find a strong correlation between the DM minor axis and the angular momentum vector, as the two tend to be aligned. Therefore, the halo’s major axis is tangent to the baryonic disk. This is more obvious in this analysis compared to the Galactic analysis of Sec. 2.3.3 since Galactic DM signals are more spherical and therefore harder to orient in a particular direction; correlating the direction of a mostly spherical signal is done at random (See Sec. 2.3.3 for a comparison with previous work.).

2.4.4 Halo Mergers / Subhalos

Figure 2-9: Logarithmic map for DM annihilation for halo labeled “499” in the Illustris-1 simulation from the point of view of an observer external to the halo. We used HEALPix with nside (see Eq. 2.4). The observer is located at a distance kpc of the center. The halo mass is .

Many of the halos in the simulation have experienced a recent merger or interaction; an example is shown in Fig. 2-9. To test the effect of these mergers/large subhalos on our sphericity distributions, we study the observed axis ratio in two different sets of subsamples of the data.

First, we omit from the analysis halos where the second-largest subhalo (the first one being the host halo) has a mass fraction higher than () of the total mass of the halo. As an example, halo “499”, shown in Fig. 2-9, encompasses the main host halo of mass fraction 0.49, and a second subhalo of mass fraction 0.44. Removing these halos leads to a more steeply falling axis ratio distribution for small axis ratios, , as shown in Fig. 2-10; compared with the original distribution in Fig. 2-6, the low-axis-ratio longer tail of the distribution is diminished. When the cut is strengthened to remove all subhalos with more than of the total mass of the halo, this tail is removed almost completely.

Second, we perform the quadrant analysis on the inner of the halo, shown in Fig. 2-7, which should only pick out the subhalo with the deepest potential well, as the location of halos in the Illustris simulation is set by the most bound particle. The distributions of and are peaked closer to zero (sphericity) when considering only pixels within the inner of the halo. The distribution is still fairly flat and not especially peaked at near-sphericity.

We see that a non-negligible fraction of the halos are expected to have elongated DM distributions due to recent mergers and/or massive subhalos. In many cases, the presence of such mergers should be apparent from the baryonic matter, but in cases where the merging halo was a low-mass system, the peak of the annihilation/decay signal might be substantially displaced from the center of the potential well inferred from the baryonic matter. This is consistent with previous work (see for example Ref. [Moore:2003eq]). Alternatively, one can also try to understand the virialization of the halos through a virialization parameter such as the one given in Ref. [Wise:2007jc], though we do not do so in this work as it is computationally intensive.

Figure 2-10: Histogram of the newly defined observed axis ratio for annihilation and decay for extragalactic sources, comparing both Illustris-1 and Illustris-1-Dark (see Eq. 2.6) after having omitted the “merger” halos (see text).

2.5 Comparison to Photon Data

In this section, we compare the anisotropy/sphericity distributions for DM halos, from the Illustris simulation, with the astrophysical backgrounds for potential DM signals.

2.5.1 Fermi data

For this analysis, we use Pass 8 data from Fermi collected between August 4, 2008 and June 3, 2015 [Atwood:2009ez, Atwood:2013rka]. We employ the recommended data quality cuts: zenith angle , instrumental rocking angle , DATA_QUAL 0, LAT_CONFIG=1. We use the Ultraclean event class and select the top quartile of events by point spread function [Portillo:2014ena]. We divide these photons into thirty equally logarithmically-spaced energy bins between 0.3 and 300 GeV; we restrict our analysis to the eight energy bins covering the range from GeV, as in this energy range the point spread function of the telescope is small and stable, and the gamma-ray excess has been clearly detected [Goodenough:2009gk, Abazajian:2012pn, Daylan:2014rsa, Calore:2014xka, TheFermi-LAT:2015kwa].

First, we analyze the full Fermi data with no additional cuts, as it is dominated by background. We pixelize the sky using HEALPix with nside , and adopt the same strategy outlined in Sec. 2.3 for the analysis of signals from within a halo, where the center of the halo is located at 8.5 kpc from the observer. Computing the -tensor defined in Eq. 2.6, we find an average axis ratio of , with a few percent spread across the different energy bins. In our default orientation, none of the 650 MW-like halos in the sample, shown in Fig. 2-5, had an axis ratio this small or smaller, in either annihilation or decay. When we tested the effect of viewing the halos from different directions, we still found no halos with this level of elongation in decay signals, but for annihilation, two halos (out of 650) attained this level of elongation for specific orientations, corresponding to 10 samples out of tests.

Second, we isolate the residual Fermi signal in the energy bin that dominates the signal GeV, and study its morphology. The residual signal map888We thank Nicholas Rodd for providing us with the residual maps. is obtained through a similar analysis strategy as that used in Ref. [Linden:2016rcf]. The region of interest in this analysis is and , as the diffuse background templates are optimized to this region. We utilize standard template fitting methods (as in [Linden:2016rcf] for example) to determine the contribution of the following templates: a uniform isotropic template, a diffuse background model by Fermi’s diffuse model p6v11, a bubbles template map and an NFW template for the DM contribution. The residual map is a HEALPix map with nside=256, obtained after subtraction of the non-DM contributions with a coefficient of their best fit. We find that the axis ratio in this region is , confirming previous results [Daylan:2014rsa, Calore:2014xka] that the signal is indeed spherical.

For a proper understanding of the origin of the Fermi signal and background, we perform the same analysis of Sec. 2.3.1 but with the distributions of gas and stars of the simulation instead of DM. 999More precisely, we compute these distributions using the formalism of DM decay. We also place the observer on the baryonic disk, defined by the plane that passes through the center of the halo and perpendicular to the angular momentum vector found in Sec. 2.3.3. We show the histograms of the axis ratio of the star and gas in Fig. 2-11. We find consistent results in which the DM is more spherical/less elongated that the gas and the stars. We note that the Fermi gamma-ray emission, which largely traces the gas distribution of the Milky Way, is still quite non-spherical compared to the gas distribution of most Illustris halos. It would be interesting to understand if this reflects a general tendency for the baryonic component of Illustris halos to be more spherical and less disk-like than in reality, at least for spiral galaxies (which are known to be difficult to reproduce in cosmological simulations [Vogelsberger:2014kha]). To do so one could refine the criteria imposed to select the Milky-Way-like halos defined in Sec. 3.4 (total mass and stellar mass) and even impose further constraints, e.g. the local dark matter surface density or the rotation curves. However that would have decreased even more the number of halos, limiting the present statistical analysis. Furthermore, any disk could appear ellipsoidal if observed at an angle, and it is worth noting that the angular momentum vectors computed in the analysis of Illustris have significant errors, and therefore the observer could be placed slightly off the disk.

Figure 2-11: Histogram of the observed axis ratio for annihilation and decay of DM in MW-like halos as defined by the requirements in Eqs. 2.11 and 2.12. We also show the histograms of the distribution of gas and stars, computed in similar manner as DM decay. We finally show the axis ratio of the Fermi background data, as well as the residual Fermi signal as discussed in Sec. 2.5.1.

2.5.2 Cluster Data

Figure 2-12: Histogram of the observed axis ratio for annihilation and decay of DM in cluster-like halos with masses larger than . We also show the histograms of the distribution of stars, computed in similar manner as DM decay signals. We finally show the axis ratio of the X-ray data. In order to make a fair comparison, we show in blue (red) arrows the location of the annihilation (decay) axis ratio of the clusters that match in mass those observed, with a cut of . We find 22 our of 352 halos that satisfy that criteria.

As an example of a potential extragalactic DM signal, we use X-ray images of 78 clusters taken by the telescope XMM-Newton [Struder:2001bh, Turner:2000jy]. In a recent analysis [Bulbul:2014sua, Bulbul:2014ala], the stacked spectrum of 73 of these 78 galaxy clusters has shown a line at keV. This sample includes clusters with a high number of counts, but to avoid the closer clusters from dominating the stacked analysis, the sample includes clusters with at least counts if the redshift and counts if the redshift is . This sample finally yielded clusters with low redshift (less than 0.35) and masses larger than and therefore we compare them to Illustris maps computed at .

The X-ray images are obtained from XMM-Newton data101010We thank Esra Bulbul for providing us with the X-ray images. She should be contacted for any image requests., with a field-of-view with radius 14 arcminutes, and an angular resolution of 6 arcseconds; typically the clusters in this sample have a radius of a few arcminutes. We set the center of the cluster to be the center of mass, where pixel brightness is the mass equivalent. We then compute the -tensor given by Eq. 2.6 across a rectangle of pixels centered around .

In Fig. 2-12, we show the normalized distribution of the observed axis ratios in the cluster data, alongside the DM annihilation and decay signals expected for halos of masses larger than in order to increase statistics. We show the halos with the mass cut that matches that of the cluster data as arrows in blue (red) for annihilation (decay) in Fig. 2-12. Although with lower statistics, the sample with the same mass cut as the observed data as well as the extended sample show a tendency for axes ratios to extend to lower values that the observed data.

As we explain in App. A.2, there is a trend for smaller halos to be more spherical, so we impose a mass cut on the Illustris halos to compare to the cluster data. However, note that the cluster sample may not be a representative sample of all clusters of similar masses. The observed clusters are quite symmetric about their centers of mass, and so it would be difficult to distinguish the astrophysical X-ray emission from a DM signal based on gross morphology alone, although the most spherical clusters (in the 0.9-1 bin) appear more symmetric than () of the halos studied in annihilation (decay) signals. In the same figure, we show the normalized distribution of the stars. The gas distribution was not included since it does not reproduce observational constraints in the case of clusters [Genel:2014lma]. The X-ray data appears to be even more spherical than the star population, so indirect detection studies should not assume a spherical morphology for DM signals. More specific studies, such as analyzing possible signals using gravitational lensing, are required to understand extragalactic DM signals [Graham:2015yga].

2.6 Conclusions

In this chapter, we studied morphological properties of DM Galactic and extragalactic indirect detection annihilation/decay signals, using the high-statistics Illustris simulation to map out the expected distribution of those properties. To understand the morphology of DM signals, we introduced two parametrizations for the asymmetry/elongation of an observed signal. The first is an analog of the inertia tensor called the -tensor; it weighs every pixel’s contribution to the inertia tensor with the observed (DM) brightness. We also divided the sky into quadrants and studied the ratios of predicted signal brightness across opposite and adjacent quadrants. The advantage of these two methods is they are compatible with indirect detection observations.chapter

We explored the DM signal morphology in two cases. In the first scenario, the observer is situated inside the halo, at a distance of 8.5 kpc from its center. In this analysis we showed both results for the full halo sample, and for a subsample with DM mass and stellar mass comparable to the Milky Way. In this case, annihilation and decay signals are expected to be fairly symmetric, with the distribution of observed axis ratio peaking at . Halo substructure is more prominent in DM annihilation and the predicted signals appear slightly less symmetric compared to the case of decay, but these effects are minor. Baryons also play a role in making decay and annihilation distributions appear more spherical, but the effects are generally quite small, as the fraction of halos at any given axis ratio changes by a few percent. We find that our results are fairly robust when the center or outer regions of the halos are excluded, and are only slightly affected by the presence of baryons.

In the second scenario we studied, the observer is external to the halo; this is the relevant analysis for searches for extragalactic DM signals. Both decay and annihilation signals are more frequently non-spherical than in the Galactic case; this is especially true for annihilation, where the distribution of axis ratio is very flat, and a sizable fraction of halos have a small axis ratio in the range . We believe that this tail can be largely attributed to halos possessing large subhalos, possibly due to recent halo mergers. Once halos with a substantial second subhalo are removed, the peak of the distribution shifts towards values of the axis ratio closer to 1.

We examined the possible correlation between the baryonic disk and the principal axis of the decay/annihilation signal. We found that in the Galactic analysis, the signal’s minor axis tended to be aligned with the angular momentum vector of the baryons, i.e. the direction perpendicular to the baryonic plane, although this correlation was quite mild (depending on the method of calculation, there were roughly more halos with than expected from the uncorrelated case, where is the angle between the DM signal’s minor axis and the angular momentum vector). In the extragalactic analysis, we find a stronger correlation between the direction of the minor axis and the angular momentum vector of the baryons, as there is an excess of over the flat distribution for an angle between the halo’s minor axis and the angular momentum vector. We think that the correlation is more pronounced in the extragalactic case first because the halos are largely non-spherical and therefore do have a preferred direction that does correlate with the baryons. This correlation is more pronounced for halos with a massive second subhalo. We think that this is due to the process of virialization; after the merger has occurred, the new subhalo is slowly getting virialized with the rest of the halo, and that process is sensitive to the presence of the baryons.

Finally, we used two sets of observational data to study the degree to which DM signals might be distinguishable from astrophysical backgrounds: gamma-ray data from the Fermi Gamma-Ray Space Telescope as an example for the Galactic analysis, and X-ray cluster data as a case study of potential extragalactic signals. The Fermi all-sky data in the 2-12 GeV band have an axis ratio , which is smaller than the axis ratio for Galactic decay signals from all tested MW-like halos (650 halos in 12 different orientations), and smaller than the axis ratio for Galactic annihilation signals more than of the time. When we remove estimates of the astrophysical backgrounds and examine the “GeV excess”, focusing on the region around the Galactic center, we find that the residual is almost perfectly spherical, consistent with expectations for possible Galactic DM signals. Compared to the distributions of gas and stars, the background is closer to being part of the gas distribution, while the signal is more likely a DM signal. It is however difficult to exactly reproduce the MW morphology with the Illustris simulation. In contrast, the cluster X-ray maps are quite spherical, suggesting that it would be difficult to reliably exclude a DM origin for signals distributed similarly to the background, based on this approach alone.

To summarize, this study quantifies the degree of asymmetry to be expected in Galactic and extragalactic signals of DM annihilation or decay, putting the use of morphological data to separate potential signals from astrophysical background on a firmer footing.

Chapter 3 Empirical Determination of Dark Matter Velocities using Metal-Poor Stars

3.1 Introduction

The velocity distribution of dark matter (DM) in the Milky Way provides a fossil record of the galaxy’s evolutionary history. In the CDM paradigm, the Milky Way’s DM halo forms from the hierarchical merger of smaller subhalos [1978ApJ...225..357S]. As a subhalo falls into, and then orbits, its host galaxy, it is tidally disrupted and continues to shed mass until it completely dissolves. With time, this tidal debris virializes and becomes smoothly distributed in phase space. Debris from more recent mergers that has not equilibrated can exhibit spatial or kinematic substructure [2005ApJ...635..931B, 2010MNRAS.406..744C, 2011ApJ...733L...7H, 2007AJ....134.1579K, 2009ApJ...694..130M, 2009MNRAS.399.1223S, 2010A&ARv..18..567K, Lisanti:2011as, Kuhlen:2012fz, Lisanti:2010qx].

Knowledge of the DM velocity distribution is required to interpret results from direct detection experiments [PhysRevD.31.3059, Drukier:1986tm], which search for DM particles that scatter off terrestrial targets. The scattering rate in these experiments depends on both the local number density and velocity of the DM [1996PhR...267..195J, Freese:2012xd]. In the Standard Halo Model (SHM), the velocity distribution is modeled as a Maxwell-Boltzmann, which assumes that the DM distribution is isotropic and in equilibrium [Drukier:1986tm]. Deviations from these assumptions can be important for certain classes of DM models (see [Freese:2012xd] for a review).

-body simulations, which trace the build-up of Milky Way–like halos in a cosmological context, do find differences with the SHM. In DM-only simulations, this is most commonly manifested as an excess of high-velocity particles as compared to a Maxwellian distribution with the same peak velocity [Vogelsberger:2008qb, MarchRussell:2008dy, Kuhlen:2009vh]. However, full hydrodynamic simulations, which include gas and stars, find that the presence of baryons makes the DM halos more spherical and the velocities more isotropic, consistent with the SHM [Ling:2009eh, Kuhlen:2013tra, Bozorgnia:2016ogo, Kelso:2016qqj, Sloane:2016kyi].

In this chapter, we demonstrate that the DM velocity distribution can be empirically determined using populations of metal-poor stars in the Solar neighborhood. This proposal relies on the fact that these old stars share a merger history with DM in the CDM framework, and should therefore exhibit similar kinematics. The hierarchical formation of DM halos implies that the Milky Way’s stellar halo also formed from the accretion, and eventual disruption, of dwarf galaxies [Searle:1978gc, Johnston:1996sb, Helmi:1999uj, Helmi:1999ks, Bullock:2000qf, Bullock:2005pi]. For example, the chemical abundance patterns of the stellar halo can be explained by the accretion—nearly 10 Gyr ago—of a few  M DM halos hosting dwarf-irregular galaxies [Robertson:2005gv, Font:2005qs, Font:2005rm]. The stars from these accreted galaxies would have characteristic chemical abundances.

A star’s abundance of iron, Fe, and -elements (O, Ca, Mg, Si, Ti) depends on its host galaxy’s evolution. Core-collapse supernova (SN), like Type II, result in greater -enrichment relative to Fe over the order of a few Myr. Thermonuclear SN, such as Type Ia, however, act on longer time scales and produce large amounts of Fe relative to elements. For a galaxy that experiences only a brief star-formation period, the enrichment of its interstellar medium is dominated by explosions of core-collapse SN, suppressing Fe abundances. Observations indicate that the Milky Way’s inner stellar halo, which extends out to 20 kpc, is metal-poor, with an iron abundance of and -enhancement of  [1991AJ....101.1835R, 1991AJ....101.1865R, 1995AJ....109.2757M, 2006ApJ...636..804A, 2004AJ....128.1177V, Ivezic:2008wk].111The stellar abundance of element relative to is defined as: where is the number density of the element.

To demonstrate the correlation between the stellar and DM velocity distributions, we use the Eris simulation, one of the highest resolution hydrodynamic simulations of a Milky Way–like galaxy [Guedes:2011ux]. We show that the velocity distribution of metal-poor halo stars in Eris successfully traces that of the virialized DM component in the Solar neighborhood. Using results from the Sloan Digital Sky Survey (SDSS), we then infer the local velocity distribution for the smooth DM component in our Galaxy. The result differs from the SHM in important ways, and suggests that current limits on spin-independent DM may be too strong for masses below 10 GeV.

3.2 The Eris Simulation

Eris is a cosmological zoom-in simulation that employs smoothed particle hydrodynamics to model the DM, gas, and stellar distributions in a Milky Way–like galaxy from to today [Guedes:2011ux, Guedes:2012gy]. It employs the TreeSPH code Gasoline [Wadsley:2003vm] to simulate the evolution of the galaxy in a WMAP cosmology [Spergel:2006hy]. The mass resolution is and  M for each DM and gas ‘particle,’ respectively. An overview of the simulation is provided in Refs. [Guedes:2011ux, Guedes:2012gy, Pillepich:2014jfa, 2016arXiv161202832S, 2015ApJ...807..115S], and we summarize the relevant aspects for our study here.

The Eris DM halo has a virial mass of  M and radius  kpc, and experienced no major mergers after . Within , there are , , and DM, gas, and star particles, respectively. At , the DM halo hosts a late-type spiral galaxy. The disk has a scale length of 2.5 kpc and exponential scale height of 490 pc at 8 kpc from the galactic center. The properties of the Eris disk and halo are comparable to their Milky Way values [Guedes:2011ux, Pillepich:2014jfa]

A star ‘particle’ of mass  M is produced if the local gas density exceeds atoms/cm. The star formation rate depends on the gas density, , as , where is the stellar density and is the dynamical time. Metals are redistributed by stellar winds and Type Ia and Type II SNe [2016arXiv161202832S, 2015ApJ...807..115S]. The abundances of Fe and O are tracked as the simulation evolves, while the abundances of all other elements are extrapolated assuming their measured solar values [Asplund:2009fu].

Stars may either be bound to the main host halo or to its satellites when they form. We are primarily interested in the latter, as these stars share a common origin with the DM. The vast majority of halo stars in Eris originated in satellites and are older than those born in the host [Pillepich:2014jfa]. They are more metal-poor than disk stars, on average, and we take advantage of this difference to distinguish the two components in the Eris galaxy.

Figure 3-1: The density distribution as a function of Galactocentric radius for the dark matter (black) and all stars (blue) in Eris. The distributions for subsamples of stars with and are also shown (dotted brown, dashed red, and solid orange, respectively). The density of the most metal-poor stellar population exhibits the same dependence on radius as the dark matter near the Sun’s position,  kpc.

3.3 Stellar Tracers for Dark Matter

Figure 3-1 shows the density distribution of the DM and stars in Eris as a function of Galactocentric radius. The distribution for all stars is steeper than that for DM. However, this includes contributions from thin and thick disk, as well as halo stars. To select the stars that are most likely to be members of the halo, we place cuts on both the Fe and -element abundances. Figure 3-1 illustrates what happens when progressively stronger cuts are placed on [Fe/H], while keeping . As the cut on iron abundance varies from to , the density fall-off becomes noticeably more shallow.

Figure 3-2: Distributions of the three separate velocity components of the DM (solid black) and stars in Eris. The velocities are in the galactocentric frame, where the -axis is oriented along the stellar angular momentum vector. The stellar distributions are shown separately for different metallicities, with and iron abundance varying from (dotted brown) to (solid orange). The distribution for all stars—dominated primarily by the disk—is also shown (solid blue). All distributions are shown for  kpc; the DM is additionally required to lie within 2 kpc of the plane. To guide the eye, the orange shading highlights the differences between the DM and distributions. The discrepancy in the distributions is due to the preferential disruption of subhalos on prograde orbits in Eris; observations of the Milky Way halo do not see such pronounced prograde rotation [Carollo:2007xh, Bond:2009mh].

Because the focus of this work is the DM distribution in the Solar neighborhood, we consider galactocentric radii in the range  kpc, where  kpc is the Sun’s position. In this range, the DM distribution falls off as , which is essentially consistent with the best-fit power-law for the most metal-poor subsample, which falls off as . This illustrates that the stars with lower iron abundance are adequate tracers for the underlying DM density distribution (see also Ref. [Tissera:2014uwa]). The correspondence between the density distributions breaks down above  kpc, indicating a transition from the inner to the outer halo that is consistent with observations [Carollo:2007xh].

Figure 3-2 compares the velocity distribution of candidate halo stars in Eris with that of the DM.222Throughout, we define the -axis to be oriented along the angular momentum vector of the stars. For comparison, we also show the stellar distribution with no metallicity cuts; it is dominated by disk stars with a characteristic peak at  km/s and narrow dispersions in the radial and vertical directions. All distributions are shown for  kpc. Because direct detection experiments are only sensitive to DM within the Solar neighborhood, we restrict its vertical displacement from the disk to be  kpc. The stellar distributions are shown with no cut on the vertical displacement; we find that the results do not change if we restrict the metal-poor population to vertical displacements greater than 2 kpc. Unfortunately, there are too few metal-poor star particles within 2 kpc of the disk in Eris to restrict to this region.

The and distributions show an excellent correspondence between the halo stars and the DM. Indeed, as increasingly more metal-poor stars are selected, their velocity distribution approaches that of the DM exactly. We apply the two-sided Kolmogorov-Smirnov test to establish whether the DM and halo stars share the same and probability distributions. The null hypothesis that the DM and stars share the same parent distribution is rejected at 95% confidence if the -value is less than . The -values for the distributions are for , suggesting that its velocity distribution is indistinguishable from that of the DM in the radial and vertical directions.

Interpreting the distribution of azimuthal velocities requires more care. As illustrated in Fig. 3-2, the azimuthal velocities are skewed to positive values for both the DM and halo stars. The prograde rotation in the DM distribution is attributable to the ‘dark disk,’ which comprises 9% of all the DM in the Solar neighborhood in Eris [Pillepich:2014jfa]. Dark disks form from the disruption of subhalos as they pass through the galactic disk. Subhalos on prograde orbits are preferentially disrupted due to dynamical friction, leading to a co-rotating DM disk [Read:2008fh]. The effect on the stars is similar, and—indeed—more pronounced due to dissipative interactions between halo stars and the disk [Pillepich:2014jfa]. The end result is that the halo stars systematically under-predict the DM distribution at negative azimuthal velocities.

Current observations suggest that our own Milky Way has an inner halo with either modest or vanishing prograde rotation [Carollo:2007xh, Bond:2009mh], and constrain the possible contributions from a dark disk [Read:2014qva]. This suggests that the mergers that resulted in Eris’ prograde halo might not have occurred in our own Galaxy, making the comparison of the DM and halo azimuthal motions more straightforward in realization. In the absence of such mergers, we assume that the DM and metal-poor stars have distributions that match just as well as those in the and cases.

We have verified that the results presented in Fig. 3-2 are robust even as the spatial and [/Fe] cuts are varied. We consider , remove the [/Fe] cut altogether, and study the region where  kpc. In all these cases, the conclusions remain the same.

3.4 Empirical Velocity Distribution.

We now look to the kinematic properties of the Milky Way’s inner halo to infer the local DM velocities by extrapolating the correspondence argued above to our Galaxy. Spatial, chemical, and kinematic properties of the smooth inner halo have been characterized by SDSS [Juric:2005zr, Ivezic:2008wk, Bond:2009mh]. The sample includes stars with -band magnitude and heliocentric distances of 100 pc to 10 kpc that cover 6500 deg of sky at latitudes  [Bond:2009mh]. The metallicity of the halo stars is well-modeled by a Gaussian with mean and standard deviation 0.30 dex [Juric:2005zr]. The Galactic velocity distribution is provided for candidate halo stars with :


where  km/s in spherical coordinates. Over the volume probed, the velocity ellipsoid does not exhibit a tilt in the spherical coordinate system and the dispersions are constant. Additionally, the azimuthal velocities (in cylindrical coordinates) exhibit no prograde motion, in contrast to Eris.

Figure 3-3 shows the speed distribution for a mock catalog of halo stars generated using Eq. 3.1, with a spread that corresponds to varying the dispersions within their errors. The peak velocity is located at 130 km/s. For comparison, the Eris DM speed distribution is shown for the region  kpc and  kpc. ErisDark is a DM-only simulation generated with the same initial conditions as Eris and described in Ref. [Kuhlen:2013tra]; its DM speed distribution, plotted for  kpc, is included as an example of a DM-only simulation result, which typically yields lower peak speeds. For comparison, the SHM is also included in Fig. 3-3. For an isotropic dispersion (), Eq. 3.1 simplifies to the Maxwell-Boltzmann distribution . This corresponds to a collisionless isothermal distribution with density , and yields a flat rotation curve with circular velocity , where  km/s.

Figure 3-3: Galactocentric speed distribution for SDSS inner-halo stars (solid blue), generated from Eq. 3.1. For comparison, we show the Standard Halo Model (dashed red), and the dark matter speed distributions in the Eris (dot-dashed black) and ErisDark halos (dot-dashed gray). The inset shows the expected background-free 95% C.L. limit on the DM spin-independent scattering cross section, assuming the exposure and energy threshold of the LUX experiment [Akerib:2016vxi] for the SDSS and SHM velocity distributions.

If the SDSS halo stars are adequate tracers for the local DM, then Fig. 3-3 suggests that the DM speeds may be slower, on average, than what is expected in the SHM. This can lead to noticeable differences in the predicted signal rate for direct detection experiments. If a DM particle of mass scatters off a nucleus with momentum transfer and effective cross section , the scattering rate is


where is the recoil energy of the nucleus, is the local DM density, is the DM-nucleus reduced mass, is the exponential nuclear form factor [1996PhR...267..195J], is the minimum velocity needed to scatter, and is the velocity of the lab frame relative to the Galactic frame. Assuming the exposure of the LUX experiment, with kg days and a minimum energy threshold of keV [Akerib:2016vxi], we derive the 95% one-sided Poisson C.L bound (3.0 events) on the scattering cross section as a function of the DM mass. The result is shown in the inset of Fig. 3-3 for the SHM and SDSS distributions. The bounds on the lightest DM are significantly weakened when the empirical distribution is used rather than the SHM.

There are several important caveats to keep in mind. First, the SDSS distribution is obtained for candidate halo stars with , and we have yet to demonstrate that these truly probe the kinematics of the primordial population of the halo. To achieve this, we must understand how the distribution evolves as progressively tighter cuts are placed on the iron abundance. If the distribution remains stable, then we can feel confident in extrapolating the results to DM based on the intuition garnered from Eris. Second, Eq. 3.1 only describes the smooth component of the inner halo in the SDSS volume, and does not account for any spatial or kinematic substructure.

3.5 Conclusions

In this chapter, we propose that DM velocities can be determined empirically using metal-poor stars in the Solar neighborhood. Low metallicity stars are typically born in galaxies outside our own. Like DM, they are dragged into the Milky Way through mergers, and predominantly populate the halo surrounding the disk. We demonstrate the close correlation between the distributions of DM and metal-poor stars using the Eris simulation, and conclude that the kinematics of the primordial stellar population tracks that of the virialized DM. To verify the generality of these findings and understand their dependence on the merger history, this study should be repeated with other hydrodynamic simulations of Milky Way–like halos. It would also be pertinent to understand whether the correspondence holds in generalizations of CDM, such as self-interacting DM.

The velocity distribution of the smooth inner halo has been characterized by SDSS and can be used to infer the local DM velocities. The corresponding speed distribution has a lower peak velocity and smaller dispersion than what is typically assumed in the SHM. This affects predictions for the DM scattering rate in direct detection experiments. Specifically, the empirical DM distribution weakens published limits on the spin-independent cross section by nearly an order of magnitude at masses below 10 GeV. The wealth of data from Gaia [Perryman2001] will allow us to better understand whether the SDSS distribution is an accurate descriptor of the most metal-poor stars in the Solar neighborhood, and whether any additional substructure exists from recent mergers. We explore this subject in greater detail in a follow-up study [ravepaper].

Chapter 4 (In)Direct Detection of Boosted Dark Matter

4.1 Introduction

A preponderance of gravitational evidence points to the existence of dark matter (DM) [Zwicky:1933gu, Begeman:1991iy, Bertone:2010zza]. Under the compelling assumption that DM is composed of one or more species of massive particles, DM particles in our Milky Way halo today are expected to be non-relativistic, with velocities . Because of this small expected velocity, DM indirect detection experiments are designed to look for nearly-at-rest annihilation or decay of DM, and DM direct detection experiments are designed to probe small nuclear recoil energies on the order of ( is the reduced mass of the DM-nucleus system, is the nucleus mass). In addition, these conventional detection strategies are based on the popular (and well-motivated) assumption that DM is a weakly-interacting massive particle (WIMP) whose thermal relic abundance is set by its direct couplings to the standard model (SM).

In this chapter, we explore a novel possibility that a small population of DM (produced non-thermally by late-time processes) is in fact relativistic, which we call “boosted DM”. As a concrete example, consider two species of DM, and (which need not be fermions), with masses . Species constitutes the dominant DM component, with no direct couplings to the SM. Instead, its thermal relic abundance is set by the annihilation process111To our knowledge, the first use of to set the relic abundance of appears in the assisted freeze-out scenario [Belanger:2011ww]. As an interesting side note, we will find that assisted freeze-out of can lead to a novel “balanced freeze-out” behavior for . In App. B.1, we show that the relic abundance can scale like (unlike for standard freeze-out). In this chapter, of course, we are more interested in the boosted population, not the thermal relic population.


At the present day, non-relativistic particles undergo the same annihilation process in the galactic halo today, producing relativistic final state particles, with Lorentz factor . These boosted DM particles can then be detected via their interactions with SM matter at large volume terrestrial experiments that are designed for detecting neutrinos and/or proton decay, such as Super-K/Hyper-K [Fukuda:2002uc, Abe:2011ts], IceCube/PINGU/MICA [Ahrens:2002dv, Aartsen:2014oha, MICA], KM3NeT [Katz:2006wv], and ANTARES [Collaboration:2011nsa], as well as recent proposals based on liquid Argon such as LAr TPC and GLACIER [Bueno:2007um, Badertscher:2010sy], and liquid scintillator experiments like JUNO [PhysRevD.88.013008, Li:2014qca]. In such experiments, boosted DM can scatter via the neutral-current-like process


similar to high energy neutrinos. This boosted DM phenomenon is generic in multi-component DM scenarios and in single-component DM models with non-minimal stabilization symmetries), where boosted DM can be produced in DM conversion [DEramo:2010ep, SungCheon:2008ts, Belanger:2011ww], semi-annihilation (where is a non-DM state) [DEramo:2010ep, Hambye:2008bq, Hambye:2009fg, Arina:2009uq, Belanger:2012vp], self-annihilation [Carlson:1992fn, deLaix:1995vi, Hochberg:2014dra], or decay transition .

In order to be detectable, of course, boosted DM must have an appreciable cross section to scatter off SM targets. Based on Eq. (4.1) alone and given our assumption that is isolated from the SM, one might think that could also have negligible SM interactions. In that case, however, the dark sector would generally have a very different temperature from the SM sector, with the temperature difference depending on details related to reheating, couplings to the inflaton, and entropy releases in the early universe [Berezhiani:1995am, Berezhiani:2000gw, Ciarcelluti:2010zz, Feng:2008mu]. So if we want to preserve the most attractive feature of the WIMP paradigm—namely, that the thermal relic abundance of is determined by its annihilation cross section, insensitive to other details—then must have efficient enough interactions with the SM to keep in thermal equilibrium at least until freezes out. Such -SM couplings then offer a hope for detecting the dark sector even if the major DM component has no direct SM couplings.

As a simple proof of concept, we present a two-component DM model of the above type, with / now being specified as fermions. The dominant DM component has no (tree-level) interactions with the SM, such that traditional DM searches are largely insensitive to it. In contrast, the subdominant DM component has significant interactions with the SM via a dark photon that is kinetically-mixed with the SM photon. The two processes related to the (in)direct detection of the / dark sector are illustrated in Fig. 4-1. In the early universe, the process on the left, due to a contact interaction between and , sets both the thermal relic abundance of as well as the production rate of boosted in the galactic halo today. The resulting boosted population has large scattering cross sections off nuclei and electrons via dark photon exchange, shown on the right of Fig. 4-1. Assuming that itself has a small thermal relic abundance (which is expected given a large SM scattering cross section), and is light enough to evade standard DM detection bounds, then (direct) detection of boosted via (indirect) detection of annihilation would offer the best non-gravitational probe of the dark sector.222Because has no direct coupling to the SM, the solar capture rate is suppressed. By including a finite -SM coupling, one could also imagine boosted DM coming from annihilation in the sun. The possibility of detecting fast-moving DM emerging from the sun has been studied previously in the context of induced nucleon decay [Huang:2013xfa], though not with the large boost factors we envision here which enable detection via Cherenkov radiation. Note, however, that particles are likely to become trapped in the sun due to energy loss effects (see Sec. 4.4.4), limiting solar capture as a viable signal channel.

Figure 4-1: (Left) Production of boosted particles through annihilation in the galactic center: . This process would be considered “indirect detection” of . (Right) Scattering of off terrestrial electron targets: . This process would be considered “direct detection” of .

Beyond just the intrinsic novelty of the boosted DM signal, there are other reasons to take this kind of DM scenario seriously. First, having the dominant DM component annihilate into light stable particles (i.e. assisted freeze-out [Belanger:2011ww]) is a novel way to “seclude” DM from the SM while still maintaining the successes of the thermal freeze-out paradigm of WIMP-type DM.333For variations such as annihilating to dark radiation or to dark states that decay back to the SM, see for instance Refs. [Pospelov:2007mp, ArkaniHamed:2008qn, Ackerman:mha, Nomura:2008ru, Mardon:2009gw, AB_CMB]. Such a feature enables this model to satisfy the increasingly severe constraints from DM detection experiments. A key lesson from secluded DM scenarios [Pospelov:2007mp] is that it is often easier to detect the “friends” of DM (in this case ) rather than the dominant DM component itself [Bjorken:2009mm]. Second, our study here can be seen as exploring the diversity of phenomenological possibilities present (in general) in multi-component DM scenarios. Non-minimal dark sectors are quite reasonable, especially considering the non-minimality of the SM (with protons and electrons stabilized by separate - and -number symmetries). Earlier work along these lines includes, for instance, the possibility of a mirror DM sector [Hodges:1993yb, Mohapatra:2000qx, Berezhiani:2000gw, Foot:2001hc]. Recently, multi-component DM scenarios have drawn rising interest motivated by anomalies in DM detection experiments [Fairbairn:2008fb, Zurek:2008qg, Profumo:2009tb] and possible new astrophysical phenomena such as a “dark disk” [Fan:2013yva]. Boosted DM provides yet another example of how the expected kinematics, phenomenology, and search strategies for multi-component DM can be very different from single-component DM.

The outline of the rest of this chapter is as follows. In Sec. 4.2, we present the above model in more detail. In Sec. 4.3, we describe the annihilation processes of both and , which sets their thermal relic abundances and the rate of boosted DM production today, and we discuss the detection mechanisms for boosted DM in Sec. 4.4. We assess the discovery prospects at present and future experiments in Sec. 4.5, where we find that Super-K should already be sensitive to boosted DM by looking for single-ring electron events from the galactic center (GC). We summarize the relevant constraints on this particular model in Sec. 4.6, and we conclude in Sec. 4.7 with a discussion of other DM scenarios with similar phenomenology. More details are relegated to the appendices.

4.2 Two Component Dark Matter

Consider two species of fermion DM and with Dirac masses , which interact via a contact operator444Via a Fierz rearrangement, we can rewrite this operator as where .


This operator choice ensures an -wave annihilation channel [Cui:2010ud], as in Fig. 4-1, which is important for having a sizable production rate of boosted today. A UV completion for such operator is shown in Fig. (a)a in App. B.2. Other Lorentz structures are equally plausible (as long as they lead to -wave annihilation).

As an extreme limit, we assume that Eq. (4.3) is the sole (tree-level) interaction for at low energies and that is the dominant DM component in the universe today. We assume that both and are exactly stable because of separate stabilizing symmetries (e.g. a ).

The subdominant species is charged under a dark gauge group, with charge for definiteness. This group is spontaneously broken, giving rise to a massive dark photon with the assumed mass hierarchy


We will take the gauge coupling of the dark to be sufficiently large (yet perturbative) such that the process efficiently depletes and gives rise to a small thermal relic abundance (see Eq. (4.12) below).

Via kinetic mixing with the SM photon [Holdom:1985ag, Okun:1982xi, Galison:1983pa] (strictly speaking, the hypercharge gauge boson),


acquires -suppressed couplings to SM fields. In this way, we can get a potentially large cross section for to scatter off terrestrial SM targets, in particular from exchange (with large and suitable ) as in Fig. 4-1. In principle, we would need to account for the possibility of a dark Higgs boson in the spectrum, but for simplicity, we assume that such a state is irrelevant to the physics we consider here, perhaps due to a Stuckelberg mechanism for the [Stueckelberg:1900zz, Kors:2005uz] or negligible couplings of to matter fields.

The parameter space of this model is defined by six parameters


Throughout this chapter, we will adjust to yield the desired DM relic abundance of , assuming that any DM asymmetry is negligible. Because the process has homogeneous scaling with and , the dominant phenomenology depends on just the three mass parameters: , , and . To achieve a sufficiently large flux of boosted particles, we need a large number density of particles in the galactic halo. For this reason, we will focus on somewhat low mass thermal DM, with typical scales:


Constraints on this scenario from standard DM detection methods are summarized later in Sec. 4.6. This includes direct detection and CMB constraints on the thermal relic population. In addition, can acquire couplings to through a -loop, thus yielding constraints from direct detection of , and we introduce a simple UV completion for Eq. (4.3) in App. B.2 which allows us to compute this effect without having to worry about UV divergences.

There are a variety of possible extensions and modifications to this simple scenario. One worth mentioning explicitly is that and/or could have small Majorana masses which lead to mass splittings within each multiplet (for this would appear after breaking) [TuckerSmith:2001hy, Cui:2009xq]. As discussed in Refs. [Finkbeiner:2009mi, Graham:2010ca, Pospelov:2013nea], both components in an inelastic DM multiplet can be cosmologically stable, such that the current day annihilation is not suppressed. These splittings, however, would typically soften the bounds on the non-relativistic component of from conventional direct detection experiments, since the scattering would be inelastic (either endothermic or exothermic). This is one way to avoid the direct detection of bounds discussed in Sec. 4.6.

4.3 Thermal Relic Abundances and Present-Day Annihilation

To find the relic density of /, we need to write down their coupled Boltzmann equations. In App. B.1, we provide details about this Boltzmann system (see also Refs. [Belanger:2011ww, Bhattacharya:2013hva, Modak:2013jya]), as well as analytic estimates for the freeze-out temperature and relic abundance in certain limits. Here, we briefly summarize the essential results.

The annihilation channel not only determines the thermal freeze-out of the dominant DM component but also sets the present-day production rate for boosted particles in Milky Way. Considering just the operator from Eq. (4.3), the thermally-averaged cross section in the -wave limit is:


As discussed in App. B.1, the Boltzmann equation for effectively decouples from when . In this limit, the relic density takes the standard form expected of WIMP DM (assuming -wave annihilation):