Boosting (In)direct Detection of Dark Matter

TL;DR

This thesis investigates dark matter detection signals using simulations and proposes a new boosted dark matter model with unique experimental signatures, improving understanding of detection prospects and constraints.

Contribution

It introduces a novel boosted dark matter model and provides empirical velocity distributions, enhancing detection strategies and theoretical understanding.

Findings

Galactic signals are symmetric; extragalactic signals are asymmetric.

Dark matter velocity distribution derived from SDSS data weakens detection limits.

Proposed boosted dark matter model predicts hybrid detection signatures.

Abstract

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 $\lesssim 10$ 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.

Tables

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

Simulation	$m_{\mathrm{DM}}\left(M_{☉}\right)$	$m_b\left(M_{☉}\right)$	${ϵ}_{\mathrm{DM}}\left(\mathrm{kpc}\right)$	${ϵ}_b\left(\mathrm{kpc}\right)$
Illustris-1	$6.3\times {10}^6$	$1.3\times {10}^6$	1.4	0.7
Illustris-1-Dark	$7.5\times {10}^6$	$-$	1.4	$-$

Table 2. Table 4.1: List of experiments studied in this chapter, their angular resolutions θ e res superscript subscript 𝜃 𝑒 res \theta_{e}^{\rm res} on the Cherenkov-emitted electron direction, and the typical minimum energy threshold E e thresh superscript subscript 𝐸 𝑒 thresh E_{e}^{\rm thresh} of the detected electron. We note here that IceCube has too high of an energy threshold for our analysis, but we are interested in its future low-energy extensions such as PINGU and MICA. For PINGU, we have scaled the nominal volume (1 MTon) down by a factor of 2 to estimate particle identification efficiency. The MICA values are speculative at present, since there is not yet a technical design report.

Experiment	Volume (MTon)	$E_e^{\mathrm{thresh}}$ (GeV)	${\theta }_e^{\mathrm{res}}$ (degree)	Refs.
Super-K	$2.24\times {10}^{-2}$	0.01	$3^{\circ }$	[4]
Hyper-K	$0.56$	0.01	$3^{\circ }$	[63]
IceCube	${10}^3$	100	${30}^{\circ }$	[261, 262]
PINGU	$0.5$	1	${23}^{\circ }\left(\text{at GeV scale}\right)$	[191]
MICA	$5$	0.01	${30}^{\circ }\left(\text{at 10 MeV scale}\right)$	[192, 263]

Table 3. Table 5.1: Detectors included in this analysis. We use the exposure time of Super-K as a reference for comparison with the rest of the experiments.

Name	Number target $e^{-}$	Energy Threshold	Angular Resolution	Exposure Time	Refs.
		(MeV)	(deg)	(years)
Super-K	$7.45\times {10}^{33}$	100	5	13.6	[61]
Hyper-K	$1.86\times {10}^{35}$	100	5	13.6	[189, 63]
DUNE-10 kton	$2.70\times {10}^{33}$	30	1	13.6	[65]
DUNE-40 kton	$1.08\times {10}^{34}$	30	1	13.6	[65]

Table 4. Table 5.2: A summary of the production frequency of free hadrons in collisions between atmospheric electron (anti)neutrinos and argon-40.

Final State Hadron	0 Produced ($\%$)	1 Produced ($\%$)	$>$ 1 Produced ($\%$)
p	17.7	50.4	31.8
n	36.6	33.8	29.6
${\pi }^{\pm ,0}$	73.0	21.2	5.8
$K^{\pm ,0}$	99.4	0.5	0.1
Heavier Hadrons	98.9	1.1	0.00

Table 5. Table 5.3: Kinetic energy thresholds for DUNE to be able to detect various hadrons [ 7 ] .

Hadron	Detection Threshold (MeV)
p	21
${\pi }^{\pm ,0}$	10
$K^{\pm ,0}$	17

Table 6. Table 5.4: Table of dSphs’s locations, distances and J-factors, compiled in Refs. [ 8 , 9 ]

Name	l	$b$	Distance (kpc)	${\log }_{10}\left(J_{\text{ann}}\right)$	${\log }_{10}\left(J_{\text{dec}}\right)$	Refs.
	(deg)	(deg)	(kpc)	(${\log }_{\text{10}}$ [GeV2 cm-5])	(${\log }_{\text{10}}$ [GeV cm-2])
Bootes I	358.1	69.6	66	$18.8\pm 0.22$	$17.9\pm 0.26$	[331]
Carina	260.1	-22.2	105	$18.1\pm 0.23$	$17.9\pm 0.17$	[332]
Coma Berenices	241.9	83.6	44	$19.0\pm 0.25$	$18.0\pm 0.25$	[333]
Draco	86.4	34.7	76	$18.8\pm 0.16$	$18.5\pm 0.12$	[334]
Fornax	237.1	-65.7	147	$18.2\pm 0.21$	$17.9\pm 0.05$	[332]
Hercules	28.7	36.9	132	$18.1\pm 0.25$	$16.7\pm 0.42$	[333]
Reticulum II	265.9	-49.6	32	$19.6\pm 1.0$	$18.8\pm 0.7$	[335, 330]
Sculptor	287.5	-83.2	86	$18.6\pm 0.18$	$18.2\pm 0.07$	[332]
Segue 1	220.5	50.4	23	$19.5\pm 0.29$	$18.0\pm 0.31$	[336]
Sextans	243.5	42.3	86	$18.4\pm 0.27$	$17.9\pm 0.23$	[332]
Ursa Major I	159.4	54.4	97	$18.3\pm 0.24$	$17.6\pm 0.38$	[333, 9]
Ursa Major II	152.5	37.4	32	$19.3\pm 0.28$	$18.4\pm 0.27$	[333]
Ursa Minor	105.0	44.8	76	$18.8\pm 0.19$	$18.0\pm 0.16$	[334]
Willman 1	158.6	56.8	38	$19.1\pm 0.31$	$17.5\pm 0.84$	[337, 338]

Equations

$$\varphi (x,t)={\eta }_{\varphi }\cos (m_{\varphi }(t-\vec{v}.\vec{x})),$$

$$J=\{\begin{array}{ll}\int {\rho }^{2}dsd\Omega & \text{for annihilation},\\ \int \rho dsd\Omega & \text{for decay},\end{array}$$

$$F_s(\vec{r})=\int F(\overrightarrow{r^{\prime }})W(\vec{r}-\overrightarrow{r^{\prime }},d)d{r}^{\prime },$$

$$w(q)=\frac{8}{\pi d^3}\{\begin{array}{ll}1-6q^2+6q^3;& 0\le q\le \frac{1}{2}\\ 2(1-q{)}^3;& \frac{1}{2}<q\le 1\\ 0;& q>1.\end{array}$$

$$n_{\mathrm{pix}}=12\times {\text{nside}}^2,$$

$$I_{i,j}=\sum _n{x}_{n,i}{x}_{n,j},$$

$${\mathcal{J}}_{i,j}=\sum _{n}J(z_{n,i},z_{n,j})z_{n,i}{z}_{n,j},$$

$$J_k=\sum _{i\in R_k}{J}_i,$$

$$R_{\mathrm{opp}}$$

$$R_{\mathrm{adj}}$$

$${\vec{L}}_3=\sum _{i\in \text{SG}}{m}_{\mathrm{gas}}{\vec{r}}_i\times \vec{v},$$

$$5\times 10^{11}{M}_{\odot }<M_{200}<2.5\times 10^{12}{M}_{\odot },$$

$$4.5\times 10^{10}{M}_{\odot }<M_{\text{Stars}}<8.3\times 10^{10}{M}_{\odot }.$$

$$r=2R_{200},$$

$$[X/Y]={\log }_{10}\left(N_X/N_Y\right)-{\log }_{10}{\left(N_X/N_Y\right)}_{\odot },$$

$$f(\mathbf{v})=\frac{1}{{\left(2\pi \right)}^{3/2}{\sigma }_r{\sigma }_{\theta }{\sigma }_{\varphi }}\exp \left[-\frac{v_r^2}{2{\sigma }_r^2}-\frac{v_{\theta }^2}{2{\sigma }_{\theta }^2}-\frac{v_{\varphi }^2}{2{\sigma }_{\varphi }^2}\right],$$

$$\frac{dR}{d{E}_{\text{nr}}}=\frac{{\rho }_{\chi }}{2m_{\chi }\mu }\sigma (q^2)F(q){\int }_{v_{\text{min}}}^{\infty }\frac{f\left(\text{v}+{\text{v}}_{\text{obs}}(t)\right)}{v}{d}^{3}v,$$

$${\psi }_A{\overline{\psi }}_A\to {\psi }_B{\overline{\psi }}_B.$$

$${\psi }_{B}X\to {\psi }_B{X}^{(\prime )},$$

$$-\frac{1}{4\Lambda^{2}}\Bigl{(}\overline{\psi}_{A}\psi_{A}\overline{\psi}_{B}\psi_{B}+\overline{\psi}_{A}\gamma^{\mu}\psi_{A}\overline{\psi}_{B}\gamma_{\mu}\psi_{B}+\frac{1}{2}\overline{\psi}_{A}\Sigma^{\mu\nu}\psi_{A}\overline{\psi}_{B}\Sigma_{\mu\nu}\psi_{B}+\overline{\psi}_{A}\gamma^{5}\psi_{A}\overline{\psi}_{B}\gamma^{5}\psi_{B}-\overline{\psi}_{A}\gamma^{\mu}\gamma^{5}\psi_{A}\overline{\psi}_{B}\gamma_{\mu}\gamma^{5}\psi_{B}\Bigr{)},$$

$${\mathcal{L}}_{\mathrm{int}}=\frac{1}{{\Lambda }^2}{\overline{\psi }}_A{\psi }_B{\overline{\psi }}_B{\psi }_A.$$

$$m_A>m_B>m_{{\gamma }^{\prime }}.$$

$$\mathcal{L}\supset -\frac{ϵ}{2}{F}_{\mu \nu }^{\prime }{F}^{\mu \nu },$$

$$\{m_A,m_B,m_{{\gamma }^{\prime }},\Lambda ,g^{\prime },ϵ\}.$$

$$m_A\simeq \mathcal{O}(10\text{GeV}),m_B\simeq \mathcal{O}(100\text{MeV}),m_{{\gamma }^{\prime }}\simeq \mathcal{O}(10\text{MeV}).$$

$$⟨{\sigma }_{A\overline{A}\to B\overline{B}}v{⟩}_{v\to 0}=\frac{1}{8\pi {\Lambda }^4}{\left(m_A+m_B\right)}^2\sqrt{1-\frac{m_B^2}{m_A^2}}.$$

$${\Omega }_A\simeq 0.2\left(\frac{5\times 10^{-26}{\text{cm}}^3/\text{s}}{⟨{\sigma }_{A\bar{A}\to B\bar{B}}v⟩}\right).$$

$$⟨{\sigma }_{A\overline{A}\to B\overline{B}}v⟩\approx 5\times 10^{-26}{\text{cm}}^3/\text{s}{\left(\frac{m_A}{20\text{GeV}}\right)}^2{\left(\frac{250\text{GeV}}{\Lambda }\right)}^4.$$

$$⟨{\sigma }_{B\overline{B}\to {\gamma }^{\prime }{\gamma }^{\prime }}v{⟩}_{v\to 0}=\frac{g^{\prime 4}}{2\pi }\frac{\left(m_B^2-m_{{\gamma }^{\prime }}^2\right)}{(m_{{\gamma }^{\prime }}^2-2m_B^2{)}^2}\sqrt{1-\frac{m_{{\gamma }^{\prime }}^2}{m_B^2}}.$$

$$\frac{{\Omega }_B}{{\Omega }_A}\simeq \frac{m_B}{m_A}\sqrt{\frac{⟨{\sigma }_{A\overline{A}\to B\overline{B}}v⟩}{⟨{\sigma }_{B\overline{B}\to {\gamma }^{\prime }{\gamma }^{\prime }}v⟩}}.$$

$$m_A=20\text{GeV},m_B=200\text{MeV},m_{{\gamma }^{\prime }}=20\text{MeV},g^{\prime }=0.5,ϵ=10^{-3},$$

Full text

\department

Department of Physics

\degree

Doctor of Philosophy

\degreemonth

June \degreeyear2017 \thesisdateMay 19, 2017

\supervisor

Jesse ThalerAssociate Professor

\chairman

Nergis MavalvalaAssociate Department Head of Physics

Boosting (In)direct Detection of Dark Matter

Lina Necib

{abstractpage}

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 $\lesssim 10$ 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.

Acknowledgments

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.

Contents

1. 1 Introduction

2. 1.1 Dark Matter

3. 1.1.1 Introduction

4. 1.1.2 Simulations

5. 1.1.3 WIMP Freeze-Out

6. 1.2 Spherical Cows of Dark Matter Indirect Detection

7. 1.3 Empirical Determination of Dark Matter Velocity Using Metal Poor Stars

8. 1.4 Boosted Dark Matter

9. 1.5 Additional Topics

10. 2 Spherical Cows of Dark Matter

11. 2.1 Introduction

12. 2.2 Methodology

13. 2.2.1 Illustris Simulation

14. 2.2.2 Computing $J$-factors

15. 2.2.3 Asymmetry Parameterization

16. 2.3 Galactic Analysis

17. 2.3.1 Observed Axis Ratio

18. 2.3.2 Quadrant Analysis

19. 2.3.3 Correlation with Baryon Disk

20. 2.3.4 Milky-Way-like Halos

21. 2.4 Extragalactic Analysis

22. 2.4.1 Axis Ratio

23. 2.4.2 Quadrant Analysis

24. 2.4.3 Correlation with Baryon Disk

25. 2.4.4 Halo Mergers / Subhalos

26. 2.5 Comparison to Photon Data

27. 2.5.1 Fermi data

28. 2.5.2 Cluster Data

29. 2.6 Conclusions

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

31. 3.1 Introduction

32. 3.2 The Eris Simulation

33. 3.3 Stellar Tracers for Dark Matter

34. 3.4 Empirical Velocity Distribution.

35. 3.5 Conclusions

36. 4 (In)Direct Detection of Boosted Dark Matter

37. 4.1 Introduction

38. 4.2 Two Component Dark Matter

39. 4.3 Thermal Relic Abundances and Present-Day Annihilation

40. 4.4 Detecting Boosted Dark Matter

41. 4.4.1 Flux of Boosted Dark Matter

42. 4.4.2 Detection of Boosted Dark Matter

43. 4.4.3 Backgrounds to Boosted Dark Matter

44. 4.4.4 Impact of Earth Attenuation

45. 4.5 Detection Prospects for Present and Future Experiments

46. 4.5.1 Event Selection

47. 4.5.2 Signal Rates

48. 4.5.3 Background Rates

49. 4.5.4 Estimated Experiment Reach

50. 4.6 Summary of Existing Constraints

51. 4.7 Conclusions and Other Possibilities

52. 5 Extending Boosted Dark Matter to DUNE

53. 5.1 Introduction

54. 5.2 Boosted Dark Matter

55. 5.2.1 Features of Boosted Dark Matter

56. 5.2.2 Flux of Boosted Dark Matter from Annihilation

57. 5.2.3 Implications of Forward Scattering

58. 5.2.4 Constant Amplitude Limit

59. 5.3 Event Selection

60. 5.3.1 Cherenkov Detectors: Super-K

61. 5.3.2 Argon-Based Detectors: DUNE

62. 5.3.3 Detector summary

63. 5.4 Background Modeling

64. 5.4.1 Cherenkov Detectors

65. 5.4.2 LArTPC Detectors

66. 5.5 Reach at neutrino experiments

67. 5.6 Dwarf Spheroidal Analysis

68. 5.6.1 J-factor of Dwarf Galaxies

69. 5.6.2 Event Reach

70. 5.7 Conclusions

71. 6 Conclusions

72. A Spherical Cows

73. A.1 Analysis of Decaying Dark Matter

74. A.1.1 Galactic Analysis

75. A.1.2 Extragalactic Analysis

76. A.2 Mass Correlation

77. A.3 Comparison of Asymmetry Parameterization to Axis Ratio

78. B Boosted Dark Matter

79. B.1 Analytic Approximations to Relic Abundances

80. B.2 Direct Detection of Non-Boosted DM

81. B.3 Boosted DM Scattering Off Hadrons

82. B.3.1 Elastic Scattering Off Hadrons

83. B.3.2 Deep Inelastic Scattering Off Hadrons

84. B.4 Understanding Forward Scattering

85. B.5 Comparing the Full Analysis with a Concrete Model

List of Figures

1. 1.1 Evolution of the abundance of a WIMP DM $\chi$ of mass $m_{\chi}=20$ GeV, and an annihilation cross section to SM of $\sigma=5\times 10^{-26}$ cm3/s.
2. 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 $=16$ (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 $3.2\times 10^{14}M_{\astrosun}$.
3. 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).
4. 2.3 The distribution of the asymmetry parameters $R_{\mathrm{opp}}$ (left) and $R_{\mathrm{adj}}$ (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. $J$-factors are computed over all halos (blue) as well as when omitting the inner $10^{\circ}$ disk (red), and through the inner $30^{\circ}$ only (green). All regions are centered on the halo center. The inset shows a zoom-in of the region of small $R_{\text{opp}}$/$R_{\text{adj}}$.
5. 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. $\overline{N}$ is the mean of the distribution (a flat distribution in $\cos\theta$).
6. 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.
7. 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).
8. 2.7 The distribution of the asymmetry parameters $R_{\mathrm{opp}}$ (left) and $R_{\mathrm{adj}}$ (right) for DM annihilation as observed from a point outside the halo, taken from DM-only and DM+baryons simulations. $J$-factors are computed over all halos (blue) and through the inner $5^{\circ}$ only (green).
9. 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. $\overline{N}$ is the mean of the distribution (a flat distribution in $\cos\theta$).
10. 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 $=512$ (see Eq. 2.4). The observer is located at a distance $2R_{200}=478$ kpc of the center. The halo mass is $2.72\times 10^{12}M_{\astrosun}$.
11. 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).
12. 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.
13. 2.12 Histogram of the observed axis ratio for annihilation and decay of DM in cluster-like halos with masses larger than $2\times 10^{12}M_{\odot}$. 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 $M>5\times 10^{13}M_{\astrosun}$. We find 22 our of 352 halos that satisfy that criteria.
14. 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 \text{[\alpha/Fe]}>0.2 and $\text{[Fe/H]}<-1,-2,-3$ 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, $r_{\odot}\sim 8$ kpc.
15. 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 $z$-axis is oriented along the stellar angular momentum vector. The stellar distributions are shown separately for different metallicities, with \text{[\alpha/Fe]}>0.2 and iron abundance varying from $\text{[Fe/H]}<-1$ (dotted brown) to $\text{[Fe/H]}<-3$ (solid orange). The distribution for all stars—dominated primarily by the disk—is also shown (solid blue). All distributions are shown for $|r-r_{\odot}|\leq 2$ 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 $\text{[Fe/H]}<-3$ distributions. The discrepancy in the $v_{\phi}$ 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 [1, 2].
16. 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 [3] for the SDSS and SHM velocity distributions.
17. 4.1 (Left) Production of boosted $\psi_{B}$ particles through $\psi_{A}$ annihilation in the galactic center: $\psi_{A}\overline{\psi}_{A}\rightarrow\psi_{B}\overline{\psi}_{B}$. This process would be considered “indirect detection” of $\psi_{A}$. (Right) Scattering of $\psi_{B}$ off terrestrial electron targets: $\psi_{B}e^{-}\rightarrow\psi_{B}e^{-}$. This process would be considered “direct detection” of $\psi_{B}$.
18. 4.2 (a) Ratio of the abundances $\Omega_{B}/\Omega_{A}$ as a function of $g^{\prime}$, fixing $m_{A}=20~{}\text{GeV}$, $m_{B}=0.2~{}\text{GeV}$, and $m_{\gamma^{\prime}}=20~{}\text{MeV}$. The solid line is the numerical solution of the Boltzmann equation in Eq. (B.1), the dotted line is the analytic estimate from assuming independent thermal freeze-out of $\psi_{A}$ and $\psi_{B}$ (naive freeze-out), and the dashed line is the analytic estimate from Eq. (4.12) (balanced freeze-out). (b) Evolution of the co-moving abundances $Y_{A}$ and $Y_{B}$ as a function of $x=m_{B}/T$ for the benchmark in Eq. (4.13). The solid lines show the actual densities per unit entropy, while the dashed lines are the equilibrium curves.
19. 4.3 Detection channels for boosted $\psi_{B}$ in neutrino experiments. (a) Elastic scattering on electrons. (b) Elastic scattering on protons (or nuclei). (c) Deep inelastic scattering on protons (or nuclei). For Cherenkov experiments, we find that the most promising channel is electron scattering.
20. 4.4 (a) Normalized recoil electron spectrum for different benchmark scenarios. Also indicated is the maximum scattered electron energy, given by Eq. (4.25) as well as the experimental threshold of Super-K in the Sub-GeV category (See Eq. (4.42)). (b) Recoil electron angular distribution for the same signal benchmarks, assuming a $\psi_{B}$ particle coming directly from the GC. The cutoff angle $\theta^{\prime\text{thresh}}_{e}$ is obtained by substituting the 100 MeV energy threshold into Eq. (4.35).
21. 4.5 Energy spectrum of signal and background events, normalized to the expected event yield over one year. The blue dashed line corresponds to the naive formula in Eq. (4.31) for the number of signal events in a $10^{\circ}$ search cone. The solid blue line is spectrum obtained from the convolution in Eq. (). The background spectrum of CC $\nu_{e}$ and $\overline{\nu}_{e}$ events comes from Super-K [4], scaled by a factor $\pi(10^{\circ})^{2}/(4\pi)$ to account for the nominal $10^{\circ}$ search cone. Note that data is available only for $E_{e}>100~{}\text{MeV}$, which is the same experimental threshold given in Eq. (4.42). Also indicated is the maximum scattered electron energy, given by Eq. (4.25).
22. 4.6 Angles involved in boosted DM detection. When a $\psi_{B}$ particle arrives at an angle $\theta_{B}$ from the GC, it scatters to produce an electron at angle $\theta_{e}$ with respect to $z$ ($\theta_{e}^{\prime}$ and $\phi^{\prime}_{e}$ with respect to $z^{\prime}$). To better isolate the signal from the uniform atmospheric neutrino background, we impose a search cone of half-angle $\theta_{C}$.
23. 4.7 Number of $\psi_{B}e^{-}\rightarrow\psi_{B}e^{-}$ signal events in Super-K, Hyper-K, PINGU, and MICA in the $m_{A}$/$m_{B}$ plane, for $m_{\gamma^{\prime}}=20~{}\text{MeV}$ (top) and $m_{\gamma^{\prime}}=50~{}\text{MeV}$ (bottom). The indicated regions are for 1 (left), 10 (center), 100 (right) detected events in a one year period, normalized to the couplings $\epsilon=10^{-3}$ and $g^{\prime}=0.5$. We have imposed the angular criteria of $\theta_{C}=10^{\circ}$ and the electron energy range of $\{100~{}\text{MeV},100~{}\text{GeV}\}$ ($\{1.33~{}\text{GeV},100~{}\text{GeV}\}$ for PINGU). Also shown are model-dependent constraints on the relic $\psi_{B}$ population from Sec. 4.6: the solid gray lines are from CMB heating (shown only for $g^{\prime}=0.5$), and the dashed gray lines are from DAMIC direct detection (which are independent of $g^{\prime}$, but can be eliminated by adding an inelastic splitting). The red star indicates the benchmark in Eq. (4.13).
24. 4.8 Yearly signal significance in the Sub-GeV category for our benchmark in Eq. (4.13) as a function of the search cone angle $\theta_{C}$. The peak around $10^{\circ}$ is seen for other parameter choices as well.
25. 4.9 Signal significance at Super-K, Hyper K, PINGU and MICA on the $m_{A}$/$m_{B}$ plane, for $m_{\gamma^{\prime}}=20~{}\text{MeV}$ (left) and $m_{\gamma^{\prime}}=50~{}\text{MeV}$ (right), fixing $\epsilon=10^{-3}$ and $g^{\prime}=0.5$. Shown are the $2\sigma$ reaches with 10 years of data, taking $\theta_{C}=10^{\circ}$ and adding the significances of the $E_{e}\in\{100~{}\text{MeV},1.33~{}\text{GeV}\}$ and $E_{e}\in\{1.33~{}\text{GeV},100~{}\text{GeV}\}$ categories in quadrature (only the latter for PINGU). Also shown is the current $2\sigma$ exclusion using all-sky data from Super-K, where we assume a 10% uncertainty on the background. The grey model-dependent limits are the same as in Fig. 4.7: the solid gray lines are constraints on $\psi_{B}$ from CMB heating and the dashed gray lines are from DAMIC. The red star is the benchmark from Eq. (4.13).
26. 4.10 Same as Fig. 4.9, but on the $g^{\prime}$/$m_{B}$ plane, for $m_{A}=20~{}\text{GeV}$ (left) and $m_{A}=50~{}\text{GeV}$ (right), fixing $\epsilon=10^{-3}$ and $m_{\gamma^{\prime}}=20~{}\text{MeV}$. The spikes in the CMB heating bounds (solid gray lines) are from Sommerfeld resonances.
27. 5.1 Scattering process of BDM $B$ off of electrons.
28. 5.2 Annihilation process that produces $B$ with a Lorentz boost.
29. 5.3 Geometry of a search cone for incoming $B$ particles originating at the GC and scattering off electrons at a neutrino experiment [5].
30. 5.4 $95\%$ limits on parameter space for BDM annihilation for Super-K, Hyper-K and DUNE for 10 and 40 kton in volume. The gray region is excluded by the fact that no excess has been detected in Super-K in the past 11 year data set.
31. 5.5 Map of the dSphs’ locations in Galactic coordinates used in this analysis. The center of the figure is the GC.
32. 5.6 $95\%$ limits on parameter space for BDM annihilation in a stacked analysis of dSphs. The gray region is excluded by the fact that no excess has been detected in Super-K in the past 11 year data set.
33. A.1 As Fig. 2.3, except for decay rather than annihilation.
34. A.2 As Fig. 2.7, except for decay rather than annihilation.
35. A.3 Histograms of the observed axis ratio for annihilation (left) and decay (right) for different mass bins: $M_{200}>2\times 10^{12}M_{\odot}$, $10^{10}M_{\odot}<M_{200}<2\times 10^{12}M_{\odot}$ and $M_{200}<10^{10}M_{\odot}$.
36. A.4 Histograms of the angle between the halo’s minor axis and the angular momentum vector of the baryons in annihilation (left) and decay (right) for three different mass bins: $M_{200}>2\times 10^{12}M_{\odot}$, $10^{10}M_{\odot}<M_{200}<2\times 10^{12}M_{\odot}$ and $M_{200}<10^{10}M_{\odot}$.
37. B.1 Left: Direct detection mechanism for $\psi_{A}$ via a $\psi_{B}$-$\phi$ loop. Right: Scattering cross section of $\psi_{A}$ on nucleons, sweeping $m_{B}=0.1~{}\text{GeV}$–$3~{}\text{GeV}$ and fixing $g^{\prime}=0.5$ and $\epsilon=10^{-3}$. Also shown are the current LUX limit (gray hashes).
38. B.2 Maximum observed angle of the electron $\theta_{e}$ as a function of the initial angle at which the Boosted particle $B$ was emitted for different values of the boost factor $\gamma_{B}$.
39. B.3 Feynman diagrams for the production and detection of DM particles. (Left) Diagram that controls the abundance of $A$ in the early universe as well as today’s production of $B$ with a Lorentz boost through $A$ annihilation. (Middle) Annihilation of $B$ to $\gamma^{\prime}$, diagram that contributes to CMB limits. (Right) Signal diagram of $B$ scattering off electrons.
40. B.4 Super-K, Hyper-K and DUNE limits for the model from Ref. [5] with an exposure of 13.6 years.

List of Tables

1. 2.1 The particle masses and softening lengths for the Illustris-1 and Illustris-1-Dark simulations [6]; “DM” subscripts label DM, while “b” subscripts label baryons.
2. 4.1 List of experiments studied in this chapter, their angular resolutions $\theta_{e}^{\rm res}$ on the Cherenkov-emitted electron direction, and the typical minimum energy threshold $E_{e}^{\rm thresh}$ of the detected electron. We note here that IceCube has too high of an energy threshold for our analysis, but we are interested in its future low-energy extensions such as PINGU and MICA. For PINGU, we have scaled the nominal volume (1 MTon) down by a factor of 2 to estimate particle identification efficiency. The MICA values are speculative at present, since there is not yet a technical design report.
3. 5.1 Detectors included in this analysis. We use the exposure time of Super-K as a reference for comparison with the rest of the experiments.
4. 5.2 A summary of the production frequency of free hadrons in collisions between atmospheric electron (anti)neutrinos and argon-40.
5. 5.3 Kinetic energy thresholds for DUNE to be able to detect various hadrons [7].
6. 5.4 Table of dSphs’s locations, distances and J-factors, compiled in Refs. [8, 9]

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 [10]. 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 [11, 12, 13]. The evolution of cosmology only helped solidify the existence of DM, especially with the extremely accurate measurements of the cosmic microwave background [14, 15, 16].

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 [17]. 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 $10\%$ 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

$N-$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 $\sim 10^{3}-10^{5}$ 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) [5]; 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 $N-$body simulations that also include star and gas particles [18, 19, 20, 21, 22]. Including baryonic physics furthers our understanding of the physics of DM, as the latter might be influenced by the kinematics of the baryons.

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 $\chi\chi\rightarrow\text{SM}~{}\text{SM}$, where $\chi$ denotes the DM particle. The evolution of the DM abundance as a function of inverse temperature $T$ of the universe is shown in Fig. 1.1. Initially, this process occurs along with SM SM$\rightarrow\chi\chi$, and therefore equilibrium is maintained between the DM and SM sectors. As the universe expands and cools down, the process (SM SM$\rightarrow\chi\chi$) is less and less efficient, and therefore the number density of DM is slowly depleted, as shown in Fig. 1.1 for $x=m_{\chi}/T\lesssim 20$, where $m_{\chi}$ 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 $x\sim 20$ 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 [23]. 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 [24]. I use the cosmological simulation Illustris [25, 26, 27] 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 [28, 29, 30, 31, 32, 33, 34, 35]. 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 [36]. 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 [37, 38, 39, 40, 41], and many have theorized about the possibility of a DM origin for this excess [42, 43, 44, 45, 46, 47, 48]. I therefore compared the full gamma ray Fermi sky map in energies $\sim 1-10$ 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 [49, 50, 51]. A few studies attributed these X-rays to DM being a sterile neutrino [52, 53]. 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 [23, 54], 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) [55, 56, 2], 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 [57]. 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 [58, 59] to study the local DM velocity distribution, as well as DM substructure [60].

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 [5]. It consists of two DM components $A$ and $B$, where $A$ is heavier and dominates the DM abundance. $A$ annihilates to $B$ 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 [61, 62, 63]. We are therefore directly detecting the $B$ component, while indirectly detecting the $A$ DM particle.

The boosted $B$ 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 [64]. 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) [65]. 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 [64]. 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 $W$ and $Z$ bosons as well as top quarks from backgrounds.

In Ref. [47], 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 [38].

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) [66]. This is of particular interest since an excess of high energy neutrinos has recently been detected, the origin of which is still unknown [67]. 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 [68, 69] 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 [70], 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 $\phi\nu\nu$, inducing temporal variations on neutrino parameters [71]. The scalar $\phi$ has a local field value at the spacetime coordinate ($t$, $\vec{x}$), and that can be written as

[TABLE]

where $\eta_{\phi}$ is a quantity set by the coupling of $\phi$ to neutrinos, the mass of $\phi$, and the local DM density. $\vec{v}\sim 10^{-3}$ denotes the virial velocity of DM. The oscillations of the scalar field, although faster than $\sim\mathcal{O}(\text{s})$, can introduce distorted neutrino oscillations by shifting the square mass differences as well as the neutrino angles by a quantity proportional to $\eta_{\phi}$. From precise measurements of neutrino parameters, we can set bounds on the mass and coupling of the scalar $\phi$, over decades in masses, from $\sim 10^{-23}$ eV to $\sim 10^{-10}$ eV, expecting better improvement by the future experiments DUNE [65] and JUNO [72].

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 [73, 74, 75, 76, 77]. Some of these are based on energy correlation functions, which are functions of energy fractions and angles of the different components of the jet [78]. 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 [79]. We used power counting techniques [80, 81] to develop a new set of observables, $N_{2}$, $M_{2}$, and $N_{3}$ that show excellent discrimination power, but also stability under changes in mass and $p_{T}$ cuts. Such a property is desirable in highly desirable in experimental searches [82]. These observables are currently being investigated by the CMS collaboration of the Large Hadron Collider, with a recent analysis published just days ago [83].

Chapter 2 Spherical Cows of Dark Matter

2.1 Introduction

Gravitational evidence for dark matter (DM) is well established [10, 11, 84], yet DM still evades all other means of detection [85, 86, 8, 48, 87]. 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 [88, 89, 18, 90, 35], and can be used to predict the properties of indirect signals from DM [91, 92, 93]. Hydrodynamic simulations include baryonic matter as well as DM, and thus can probe the impact of baryonic feedback on the DM distribution [35]. 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 [94, 95, 29, 96, 97, 98, 99, 30, 100, 101]. 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 [25, 102]. 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 $\sim$3.5 keV has been found in a stacked analysis of 73 galaxy clusters [49] and in other regions [52, 53, 50, 51]. Analysis of data from the Fermi Gamma-Ray Telescope (hereafter Fermi) [36] has shown an unexplained spherically symmetric excess of $\mathcal{O}$(GeV) gamma rays at the center of the Galaxy [37, 103, 38, 40, 39, 48, 41]. 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 $\sim 3.5$ 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 $J$-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 available111http://www.illustris-project.org/ hydrodynamic simulation that traces the evolution of DM particles, as well as gas, stars and black holes across redshifts from $z=127$ to today $z=0$ [104, 105, 25, 102, 26]. 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 [25]. In this work, we focus on the last snapshot at $z=0$, which reflects the simulated state of today’s Universe [27]. 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 $(106.5{\rm~{}Mpc})^{3}$. 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 $m_{\rm{b}}$. Gravity is included with softening of the potential at small scales to avoid numerical two-body particle scattering [106]. The softening lengths for both DM and baryons, $\epsilon_{\rm{DM}}$ and $\epsilon_{\rm{b}}$, are also shown in Table 2.1.

The Illustris simulation used the friends-of-friends (FOF) algorithm to identify DM halos [107]. 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 $1.6\times 10^{5}$ and $1.5\times 10^{5}$ halos for Illustris-1 and Illustris-1-Dark respectively. The halo mass function of the Illustris simulation described in Ref. [25] 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 $\sim 5\times 10^{9}M_{\astrosun}$ to $3\times 10^{14}M_{\astrosun}$.

2.2.2 Computing $J$-factors

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

[TABLE]

where $\rho$ 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 [108, 109, 20]. More explicitly, for a field $F(\vec{r})$, one can define a smoothed interpolated version $F_{s}(\vec{r})$, related to $F$ through a kernel function $W$

[TABLE]

where $d$ defines the length of the smoothing. The kernel function $W$ approaches a delta function as $d\rightarrow 0$. We use a cubic spline to compute the local density $\rho$:

[TABLE]

The smoothing kernel in this case is $W(r,d)=w(r/2d)$, where $d$ is chosen to be the distance to the 33rd nearest neighbor of the point considered. In order to compute the density along a particular line of sight defined by the galactic coordinates $(l,b)$, 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 $d$. 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 $1/r^{2}$.

We proceed to construct sky maps of annihilation/decay $J$-factors for each halo, by placing an observer at $R_{\odot}=8.5$ kpc from the center of the halo along the $x$-axis of the simulation, and compute the $J$-factor for different values of galactic coordinates $(l,b)$. The center of the halo is defined as the location of the gravitational potential minimum. We use the package HEALPix222http://healpix.sourceforge.net to divide the sky into equal area pixels [110]. The total number of pixels in a map is defined by

[TABLE]

where nside is an input parameter that defines the pixelation.

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, $R_{200}$ 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/$h$, where $h=0.704$ is the reduced Hubble constant so that $H=h\times 100$ km/sec/Mpc. The cosmological parameters used are $\Omega_{m}=0.2726$, and $\Omega_{\Lambda}=0.7274$ [25]. is 1659 kpc$/h=2356$ kpc, much larger than the observer distance from the center of the halo $R_{\odot}=8.5$ 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 [111, 112, 94, 97, 99, 100], 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 [113]

[TABLE]

where the sum is taken over the DM particles $n$ of the halo, and $i,j$ correspond to the coordinates of the particle $n$ projected on the plane perpendicular to the observer.444See Ref. [114] 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 $x$ 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), $x_{i},x_{j}$ run over the four combinations of the $y$ and $z$ 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 $J$-factor at that location. The new inertia tensor that we call the $\mathcal{J}$-tensor is therefore

[TABLE]

where the coordinates $z_{n,i}$ 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 $J$ is given by Eq. 2.1 at a point in the sky given by the coordinates $z_{n,i}$. 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 $\mathcal{J}$-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 $J$-factor associated with each quadrant as discussed in Sec. 2.2.2

[TABLE]

where $k\in\{1,2,3,4\}$ labels the different quadrants, and $R_{k}$ is the list of pixels in quadrant $k$. $J_{i}$ is the value of the $J$-factor found at pixel $i$. 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

[TABLE]

For signals that appear spherical from the point of view of the observer, the ratios defined in Eqs. 2.8 and 2.9 $R_{\text{opp}}=R_{\text{adj}}=0$. For signals that appear strongly elongated or asymmetric, $R_{\mathrm{opp}}\rightarrow 1$ or $R_{\mathrm{adj}}\rightarrow 1$ depending on which quadrants dominate the DM signal.

2.3 Galactic Analysis

In this section, we perform a statistical analysis of the sphericity of annihilation/decay signals as observed from a location $R_{\odot}=8.5$ kpc from the halo center, similar to the Earth’s separation from the center of the Milky Way [115, 116]. 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 $\approx 0.85$, 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 [117]. 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.

Moving to the $\mathcal{J}$-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 $\approx 0.9$. 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 [28, 29, 30, 31, 32, 34, 33, 35, 118].

2.3.2 Quadrant Analysis

In this section, we compute the distributions of $R_{\text{opp}}$ and $R_{\text{adj}}$ 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. [119, 38]). For every halo, we build a map of the $J$-factors at different points in the sky in annihilation and decay following the procedure in Sec. 2.2.2. We picked nside $=16$ (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 $10^{\circ}$ and finally we look at the inner cone of half angle $30^{\circ}$. 10 degrees from the perspective of an observer at a distance of 8.5 kpc covers one softening length $\epsilon_{\mathrm{DM}}=1.4$ kpc . We show the results for DM annihilation in Fig. 2.3 as a probability distribution of $R_{\text{opp}}$ and $R_{\text{adj}}$. 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 ($\Delta J/J_{\text{total}}\gtrsim 0.3$)666By $\Delta J/J_{\text{total}}$ we mean $R_{\text{opp}}$ and $R_{\text{adj}}$. 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 $R_{\text{opp}}$ and $R_{\text{adj}}$ 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 $10^{\circ}$ with the distribution across the whole halo. The distributions of $\Delta J/J_{\text{total}}$ 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 $\theta$ 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 $\vec{L}_{3}$ for a particular halo as

[TABLE]

where $\vec{r}_{i}=\vec{x}_{i}-\vec{x_{0}}$, with $\vec{x}_{i}$ ($\vec{x_{0}}$) the coordinates of the particle $i$ (the center of the halo), and $m_{\rm{gas}}$ ($\vec{v}$) is the mass (3D velocity) of the gas cell. We then project the angular momentum vector $\vec{L}_{3}$ on the plane perpendicular to the line between the observer and the halo, and label the new 2D angular momentum vector $\vec{L}$.

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 $\mathcal{J}$-tensor in annihilation and (3) computing the $\mathcal{J}$-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 $\vec{M}$.

In order to compute the angle between the minor halo axis and the angular momentum vector, we consider the normalized inner product $\vec{M}\cdot\vec{L}/(|\vec{M}||\vec{L}|)=\cos\theta$. In Fig. 2.4, we compare the newly found distribution of the angle $\theta$ to that of a flat distribution in $\cos\theta$ 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. $|\cos\theta|\sim 1$. This implies a slight preference for the halo’s major axis to be aligned with the baryonic disk.

Our results found in projection are consistent with previous 3D analyses [120, 121, 122, 123, 124, 114, 125, 21, 126, 127, 128, 129, 130] (Ref. [113] shows a 2D projected misalignment angle of order $10^{\circ}$).

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. [131])

[TABLE]

where $M_{200}$ is the mass of the halo enclosed in a sphere with a mean density 200 times the critical density of the Universe today. $M_{\odot}$ 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 [132, 93]

[TABLE]

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 $\approx 0.8-0.9$, although with a more peaked distribution.

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 $R_{\odot}$ 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 [36], 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

[TABLE]

where $R_{200}$ 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 $r>R_{200}$. 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 $r$, the halos cover $\sim 30^{\circ}$ of the map, which is higher than most extragalactic signals, but we do so in order to resolve the inner structure.

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 $0.1-0.4$. 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 $5^{\circ}$. 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.

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

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 $10\%$ ($1\%$) 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, $\sim 0.1-0.5$, 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 $1\%$ of the total mass of the halo, this tail is removed almost completely.

Second, we perform the quadrant analysis on the inner $5^{\circ}$ 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 $R_{\text{opp}}$ and $R_{\text{adj}}$ are peaked closer to zero (sphericity) when considering only pixels within the inner $5^{\circ}$ 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. [133]). Alternatively, one can also try to understand the virialization of the halos through a virialization parameter such as the one given in Ref. [134], though we do not do so in this work as it is computationally intensive.

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 [36, 135]. We employ the recommended data quality cuts: zenith angle $<90^{\circ}$, instrumental rocking angle $<52^{\circ}$, DATA_QUAL $>$ 0, LAT_CONFIG=1. We use the Ultraclean event class and select the top quartile of events by point spread function [136]. 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 $\sim 2-12$ 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 [37, 103, 38, 39, 48].

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 $=128$, 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 $\mathcal{J}$-tensor defined in Eq. 2.6, we find an average axis ratio of $0.54$, 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 $650\times 12=7800$ tests.

Second, we isolate the residual Fermi signal in the energy bin that dominates the signal $1.89-2.38$ 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. [41]. The region of interest in this analysis is $1^{\circ}<|b|<15^{\circ}$ and $|l|<15^{\circ}$, as the diffuse background templates are optimized to this region. We utilize standard template fitting methods (as in [41] 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 $0.99$, confirming previous results [38, 39] 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 [102]). 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.

2.5.2 Cluster Data

As an example of a potential extragalactic DM signal, we use X-ray images of 78 clusters taken by the telescope XMM-Newton [137, 138]. In a recent analysis [49, 50], the stacked spectrum of 73 of these 78 galaxy clusters has shown a line at $E=3.55$ 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 $10^{5}$ counts if the redshift $z<0.1$ and $10^{4}$ counts if the redshift is $0.1<z<0.4$. This sample finally yielded clusters with low redshift (less than 0.35) and masses larger than $5\times 10^{13}M_{\odot}$ and therefore we compare them to Illustris maps computed at $z=0$.

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 $\vec{x}_{0}$ of the cluster to be the center of mass, where pixel brightness is the mass equivalent. We then compute the $\mathcal{J}$-tensor given by Eq. 2.6 across a rectangle of pixels centered around $\vec{x}_{0}$.

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 $2\times 10^{12}M_{\odot}$ 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 $98\%$ ($60\%$) 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 [26]. 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 [139].

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 $\mathcal{J}$-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 $\sim 0.9$. 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 $0.1-0.5$. 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 $4-10\%$ more halos with $\theta<0.1\pi$ than expected from the uncorrelated case, where $\theta$ 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 $\sim 62\%$ over the flat distribution for an angle $\theta<0.1\pi$ 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 $\sim 0.5$, 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 $99\%$ 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 $\Lambda$CDM paradigm, the Milky Way’s DM halo forms from the hierarchical merger of smaller subhalos [140]. 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 [141, 142, 143, 144, 145, 146, 147, 148, 149, 150].

Knowledge of the DM velocity distribution is required to interpret results from direct detection experiments [151, 152], 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 [153, 154]. 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 [152]. Deviations from these assumptions can be important for certain classes of DM models (see [154] for a review).

$N$-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 [155, 156, 157]. 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 [158, 159, 160, 118, 161].

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 $\Lambda$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 [162, 163, 164, 165, 166, 167]. For example, the chemical abundance patterns of the stellar halo can be explained by the accretion—nearly 10 Gyr ago—of a few \sim$$5\times 10^{10} M*⊙* DM halos hosting dwarf-irregular galaxies [168, 169, 170]. The stars from these accreted galaxies would have characteristic chemical abundances.

A star’s abundance of iron, Fe, and $\alpha$-elements (O, Ca, Mg, Si, Ti) depends on its host galaxy’s evolution. Core-collapse supernova (SN), like Type II, result in greater $\alpha$-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 $\alpha$ 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 $\sim$20 kpc, is metal-poor, with an iron abundance of $\text{[Fe/H]}\sim-1.5$ and $\alpha$-enhancement of $\text{[}\alpha\text{/Fe]}\sim 0.3$ [171, 172, 173, 174, 175, 56].111The stellar abundance of element $X$ relative to $Y$ is defined as:

$[X/Y]=\log_{10}\left(N_{X}/N_{Y}\right)-\log_{10}\left(N_{X}/N_{Y}\right)_{\odot}\,,$

where $N_{i}$ is the number density of the $i^{\text{th}}$ 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 [23]. 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 $\sim$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 $z=90$ to today [23, 176]. It employs the TreeSPH code Gasoline [177] to simulate the evolution of the galaxy in a WMAP cosmology [178]. The mass resolution is $9.8\times 10^{4}$ and $2\times 10^{4}$ M⊙ for each DM and gas ‘particle,’ respectively. An overview of the simulation is provided in Refs. [23, 176, 54, 179, 180], and we summarize the relevant aspects for our study here.

The Eris DM halo has a virial mass of $M_{\text{vir}}=7.9\times 10^{11}$ M⊙ and radius $R_{\text{vir}}=239$ kpc, and experienced no major mergers after $z=3$. Within $R_{\text{vir}}$, there are $7\times 10^{6}$, $3\times 10^{6}$, and $8.6\times 10^{6}$ DM, gas, and star particles, respectively. At $z=0$, 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 [23, 54]

A star ‘particle’ of mass $6\times 10^{3}$ M⊙ is produced if the local gas density exceeds $5$ atoms/cm3. The star formation rate depends on the gas density, $\rho_{\text{gas}}$, as $d\rho_{*}/dt=0.1\,\rho_{\text{gas}}/t_{\text{dyn}}\propto\rho_{\text{gas}}^{1.5}$, where $\rho_{*}$ is the stellar density and $t_{\text{dyn}}$ is the dynamical time. Metals are redistributed by stellar winds and Type Ia and Type II SNe [179, 180]. 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 [181].

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 [54]. 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.

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 $\alpha$-element abundances. Figure 3.1 illustrates what happens when progressively stronger cuts are placed on [Fe/H], while keeping \text{[\alpha/Fe]}>0.2. As the cut on iron abundance varies from $\text{[Fe/H]}<-1$ to $\text{[Fe/H]}<-3$, the density fall-off becomes noticeably more shallow.

Because the focus of this work is the DM distribution in the Solar neighborhood, we consider galactocentric radii in the range $|r-r_{\odot}|\leq 2$ kpc, where $r_{\odot}=8$ kpc is the Sun’s position. In this range, the DM distribution falls off as $\rho(r)\propto r^{-2.07\pm 0.01}$, which is essentially consistent with the best-fit power-law for the most metal-poor subsample, which falls off as $\rho(r)\propto r^{-2.24\pm 0.12}$. This illustrates that the stars with lower iron abundance are adequate tracers for the underlying DM density distribution (see also Ref. [182]). The correspondence between the density distributions breaks down above $r\gtrsim 20$ kpc, indicating a transition from the inner to the outer halo that is consistent with observations [1].

Figure 3.2 compares the velocity distribution of candidate halo stars in Eris with that of the DM.222Throughout, we define the $z$-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 $v_{\phi}\simeq 220$ km/s and narrow dispersions in the radial and vertical directions. All distributions are shown for $|r-r_{\odot}|\leq 2$ kpc. Because direct detection experiments are only sensitive to DM within the Solar neighborhood, we restrict its vertical displacement from the disk to be $|z_{\text{DM}}|\leq 2$ 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 $v_{\rho}$ and $v_{z}$ 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 $v_{\rho}$ and $v_{z}$ probability distributions. The null hypothesis that the DM and stars share the same parent distribution is rejected at 95% confidence if the $p$-value is less than $0.05$. The $p$-values for the $(v_{\rho},v_{z})$ distributions are $(0.9,0.1)$ for $\text{[Fe/H]}<-3$, 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 $\sim$9% of all the DM in the Solar neighborhood in Eris [54]. 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 [183]. The effect on the stars is similar, and—indeed—more pronounced due to dissipative interactions between halo stars and the disk [54]. 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 [1, 2], and constrain the possible contributions from a dark disk [184]. 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 $v_{\phi}$ distributions that match just as well as those in the $v_{\rho}$ and $v_{z}$ cases.

We have verified that the results presented in Fig. 3.2 are robust even as the spatial and [$\alpha$/Fe] cuts are varied. We consider \text{[\alpha/Fe]}\in[0.2,0.4], remove the [$\alpha$/Fe] cut altogether, and study the region where $|r-r_{\odot}|\leq 1$ 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 [55, 56, 2]. The sample includes stars with $r$-band magnitude $r<20$ and heliocentric distances of $\sim$100 pc to 10 kpc that cover 6500 deg2 of sky at latitudes $|b|>20^{\circ}$ [2]. The metallicity of the halo stars is well-modeled by a Gaussian with mean $\text{[Fe/H]}=-1.46$ and standard deviation 0.30 dex [55]. The Galactic velocity distribution is provided for candidate halo stars with $\text{[Fe/H]}<-1.1$:

[TABLE]

where $\{\sigma_{r},\sigma_{\phi},\sigma_{\theta}\}=\{141,85,75\}\pm 5$ 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 $1\sigma$ errors. The peak velocity is located at $\sim$130 km/s. For comparison, the Eris DM speed distribution is shown for the region $|r-r_{\odot}|\leq 2$ kpc and $|z_{\text{DM}}|\leq 2$ kpc. ErisDark is a DM-only simulation generated with the same initial conditions as Eris and described in Ref. [159]; its DM speed distribution, plotted for $|r-r_{\odot}|\leq 2$ 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 ($\sigma=\sigma_{r,\phi,\theta}$), Eq. 3.1 simplifies to the Maxwell-Boltzmann distribution $f(v)\propto e^{-v^{2}/2\sigma^{2}}$. This corresponds to a collisionless isothermal distribution with density $\rho\left(r\right)\propto r^{-2}$, and yields a flat rotation curve with circular velocity $v_{c}^{2}=2\sigma^{2}$, where $v_{c}\sim 220$ km/s.

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 $m_{\chi}$ scatters off a nucleus with momentum transfer $q$ and effective cross section $\sigma(q^{2})$, the scattering rate is

[TABLE]

where $E_{\text{nr}}$ is the recoil energy of the nucleus, $\rho_{\chi}$ is the local DM density, $\mu$ is the DM-nucleus reduced mass, $F(q)$ is the exponential nuclear form factor [153], $v_{\text{min}}$ is the minimum velocity needed to scatter, and $\textbf{v}_{\text{obs}}(t)$ is the velocity of the lab frame relative to the Galactic frame. Assuming the exposure of the LUX experiment, with $3.35\times 10^{4}$ kg days and a minimum energy threshold of $1.1$ keV [3], 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 $\text{[Fe/H]}<-1.1$, 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 $\Lambda$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 $\sim$10 GeV. The wealth of data from Gaia [185] 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 [60].

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) [10, 186, 187]. 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 $v_{\rm DM,0}\simeq\mathcal{O}(10^{-3})$. 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 $\frac{\mu^{2}}{m_{N}}v_{\rm DM,0}^{2}$ ($\mu$ is the reduced mass of the DM-nucleus system, $m_{N}$ 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, $\psi_{A}$ and $\psi_{B}$ (which need not be fermions), with masses $m_{A}>m_{B}$. Species $\psi_{A}$ 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 $\psi_{A}\overline{\psi}_{A}\to\psi_{B}\overline{\psi}_{B}$ to set the relic abundance of $\psi_{A}$ appears in the assisted freeze-out scenario [188]. As an interesting side note, we will find that assisted freeze-out of $\psi_{A}$ can lead to a novel “balanced freeze-out” behavior for $\psi_{B}$. In App. B.1, we show that the relic abundance can scale like $\Omega_{B}\propto 1/\sqrt{\sigma_{B}}$ (unlike $\Omega_{B}\propto 1/\sigma_{B}$ for standard freeze-out). In this chapter, of course, we are more interested in the boosted $\psi_{B}$ population, not the thermal relic $\psi_{B}$ population.

[TABLE]

At the present day, non-relativistic $\psi_{A}$ particles undergo the same annihilation process in the galactic halo today, producing relativistic final state $\psi_{B}$ particles, with Lorentz factor $\gamma=m_{A}/m_{B}$. 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 [61, 189], IceCube/PINGU/MICA [190, 191, 192], KM3NeT [193], and ANTARES [194], as well as recent proposals based on liquid Argon such as LAr TPC and GLACIER [195, 196], and liquid scintillator experiments like JUNO [197, 72]. In such experiments, boosted DM can scatter via the neutral-current-like process

[TABLE]

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 $\psi_{i}\psi_{j}\to\psi_{k}\psi_{\ell}$ [198, 199, 188], semi-annihilation $\psi_{i}\psi_{j}\to\psi_{k}\phi$ (where $\phi$ is a non-DM state) [198, 200, 201, 202, 203], $3\rightarrow 2$ self-annihilation [204, 205, 206], or decay transition $\psi_{i}\to\psi_{j}+\phi$.

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 $\psi_{A}$ is isolated from the SM, one might think that $\psi_{B}$ 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 [207, 208, 209, 210]. So if we want to preserve the most attractive feature of the WIMP paradigm—namely, that the thermal relic abundance of $\psi_{A}$ is determined by its annihilation cross section, insensitive to other details—then $\psi_{B}$ must have efficient enough interactions with the SM to keep $\psi_{A}$ in thermal equilibrium at least until $\psi_{A}\overline{\psi}_{A}\to\psi_{B}\overline{\psi}_{B}$ freezes out. Such $\psi_{B}$-SM couplings then offer a hope for detecting the dark sector even if the major DM component $\psi_{A}$ has no direct SM couplings.

As a simple proof of concept, we present a two-component DM model of the above type, with $\psi_{A}$/$\psi_{B}$ now being specified as fermions. The dominant DM component $\psi_{A}$ has no (tree-level) interactions with the SM, such that traditional DM searches are largely insensitive to it. In contrast, the subdominant DM component $\psi_{B}$ has significant interactions with the SM via a dark photon $\gamma^{\prime}$ that is kinetically-mixed with the SM photon. The two processes related to the (in)direct detection of the $\psi_{A}$/$\psi_{B}$ dark sector are illustrated in Fig. 4.1. In the early universe, the process on the left, due to a contact interaction between $\psi_{A}$ and $\psi_{B}$, sets both the thermal relic abundance of $\psi_{A}$ as well as the production rate of boosted $\psi_{B}$ in the galactic halo today. The resulting boosted $\psi_{B}$ population has large scattering cross sections off nuclei and electrons via dark photon exchange, shown on the right of Fig. 4.1. Assuming that $\psi_{B}$ 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 $\psi_{B}$ via (indirect) detection of $\psi_{A}$ annihilation would offer the best non-gravitational probe of the dark sector.222Because $\psi_{A}$ has no direct coupling to the SM, the $\psi_{A}$ solar capture rate is suppressed. By including a finite $\psi_{A}$-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 [211], though not with the large boost factors we envision here which enable detection via Cherenkov radiation. Note, however, that $\psi_{B}$ 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.

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 $\psi_{A}$ annihilate into light stable $\psi_{B}$ particles (i.e. assisted freeze-out [188]) 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. [212, 213, 214, 215, 216, 217]. Such a feature enables this model to satisfy the increasingly severe constraints from DM detection experiments. A key lesson from secluded DM scenarios [212] is that it is often easier to detect the “friends” of DM (in this case $\psi_{B}$) rather than the dominant DM component itself [218]. 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 $B$- and $L$-number symmetries). Earlier work along these lines includes, for instance, the possibility of a mirror DM sector [219, 220, 208, 221]. Recently, multi-component DM scenarios have drawn rising interest motivated by anomalies in DM detection experiments [222, 223, 224] and possible new astrophysical phenomena such as a “dark disk” [225]. 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 $\psi_{A}$ and $\psi_{B}$, 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 $\psi_{A}$ and $\psi_{B}$ with Dirac masses $m_{A}>m_{B}$, which interact via a contact operator444Via a Fierz rearrangement, we can rewrite this operator as

-\frac{1}{4\Lambda^{2}}\Bigl{(}\overline{\psi}_{A}\psi_{A}\overline{\psi}_{B}\psi_{B}+\overline{\psi}_{A}\gamma^{\mu}\psi_{A}\overline{\psi}_{B}\gamma_{\mu}\psi_{B}+\frac{1}{2}\overline{\psi}_{A}\Sigma^{\mu\nu}\psi_{A}\overline{\psi}_{B}\Sigma_{\mu\nu}\psi_{B}+\overline{\psi}_{A}\gamma^{5}\psi_{A}\overline{\psi}_{B}\gamma^{5}\psi_{B}-\overline{\psi}_{A}\gamma^{\mu}\gamma^{5}\psi_{A}\overline{\psi}_{B}\gamma_{\mu}\gamma^{5}\psi_{B}\Bigr{)},

where $\Sigma^{\mu\nu}=\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$.

[TABLE]

This operator choice ensures an $s$-wave annihilation channel [226], $\psi_{A}\overline{\psi}_{A}\to\psi_{B}\overline{\psi}_{B}$ as in Fig. 4.1, which is important for having a sizable production rate of boosted $\psi_{B}$ today. A UV completion for such operator is shown in Fig. B.1a in App. B.2. Other Lorentz structures are equally plausible (as long as they lead to $s$-wave annihilation).

As an extreme limit, we assume that Eq. (4.3) is the sole (tree-level) interaction for $\psi_{A}$ at low energies and that $\psi_{A}$ is the dominant DM component in the universe today. We assume that both $\psi_{A}$ and $\psi_{B}$ are exactly stable because of separate stabilizing symmetries (e.g. a $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$).

The subdominant species $\psi_{B}$ is charged under a dark $U(1)^{\prime}$ gauge group, with charge $+1$ for definiteness. This group is spontaneously broken, giving rise to a massive dark photon $\gamma^{\prime}$ with the assumed mass hierarchy

[TABLE]

We will take the gauge coupling $g^{\prime}$ of the dark $U(1)^{\prime}$ to be sufficiently large (yet perturbative) such that the process $\psi_{B}\overline{\psi}_{B}\to\gamma^{\prime}\gamma^{\prime}$ efficiently depletes $\psi_{B}$ and gives rise to a small thermal relic abundance (see Eq. (4.12) below).

Via kinetic mixing with the SM photon [227, 228, 229] (strictly speaking, the hypercharge gauge boson),

[TABLE]

$\gamma^{\prime}$ acquires $\epsilon$-suppressed couplings to SM fields. In this way, we can get a potentially large cross section for $\psi_{B}$ to scatter off terrestrial SM targets, in particular $\psi_{B}e^{-}\to\psi_{B}e^{-}$ from $\gamma^{\prime}$ exchange (with large $g^{\prime}$ and suitable $\epsilon$) as in Fig. 4.1. In principle, we would need to account for the possibility of a dark Higgs boson $H^{\prime}$ 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 $U(1)^{\prime}$ [230, 231] or negligible couplings of $H^{\prime}$ to matter fields.

The parameter space of this model is defined by six parameters

[TABLE]

Throughout this chapter, we will adjust $\Lambda$ to yield the desired DM relic abundance of $\psi_{A}$, assuming that any DM asymmetry is negligible. Because the process $\psi_{B}e^{-}\to\psi_{B}e^{-}$ has homogeneous scaling with $g^{\prime}$ and $\epsilon$, the dominant phenomenology depends on just the three mass parameters: $m_{A}$, $m_{B}$, and $m_{\gamma^{\prime}}$. To achieve a sufficiently large flux of boosted $\psi_{B}$ particles, we need a large number density of $\psi_{A}$ particles in the galactic halo. For this reason, we will focus on somewhat low mass thermal DM, with typical scales:

[TABLE]

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 $\psi_{B}$ population. In addition, $\psi_{A}$ can acquire couplings to $\gamma^{\prime}$ through a $\psi_{B}$-loop, thus yielding constraints from direct detection of $\psi_{A}$, 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 $\psi_{A}$ and/or $\psi_{B}$ could have small Majorana masses which lead to mass splittings within each multiplet (for $\psi_{B}$ this would appear after $U(1)^{\prime}$ breaking) [232, 233]. As discussed in Refs. [234, 235, 236], 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 $\psi_{A}/\psi_{B}$ 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 $\psi_{A}$/$\psi_{B}$, we need to write down their coupled Boltzmann equations. In App. B.1, we provide details about this Boltzmann system (see also Refs. [188, 237, 42]), 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 $\psi_{A}\overline{\psi}_{A}\rightarrow\psi_{B}\overline{\psi}_{B}$ not only determines the thermal freeze-out of the dominant DM component $\psi_{A}$ but also sets the present-day production rate for boosted $\psi_{B}$ particles in Milky Way. Considering just the operator from Eq. (4.3), the thermally-averaged cross section in the $s$-wave limit is:

[TABLE]

As discussed in App. B.1, the Boltzmann equation for $\psi_{A}$ effectively decouples from $\psi_{B}$ when $\langle\sigma_{B\bar{B}\rightarrow\gamma^{\prime}\gamma^{\prime}}v\rangle\gg\langle\sigma_{A\bar{A}\rightarrow B\bar{B}}v\rangle$. In this limit, the relic density $\Omega_{A}$ takes the standard form expected of WIMP DM (assuming $s$-wave annihilation):

[TABLE]

Notice that in order to get the observed DM relic abundance $\Omega_{A}\approx 0.2$, the thermal annihilation cross section is around twice the “standard” thermal cross section $3\times 10^{-26}~{}\text{cm}^{3}/\text{s}$ where a Majorana fermion DM with $\simeq 100$ GeV mass is assumed. The slight discrepancy is because our $\Omega_{A}$ is the sum of the abundances of both Dirac particles $\psi_{A}$ and $\overline{\psi}_{A}$, and the $\psi_{A}$ we are interested in has lower mass $\lesssim 20$ GeV (see, e.g., Ref. [238]).

In the limit that $m_{B}\ll m_{A}$, we have

[TABLE]

Note that $m_{A}\ll\Lambda$ for our benchmark mass $m_{A}=20~{}\text{GeV}$, so it is consistent to treat the annihilation of $\psi_{A}$ as coming just from the effective operator in Eq. (4.3).

The thermal relic abundance of $\psi_{B}$ is more subtle. In the absence of $\psi_{A}$, the relic abundance of $\psi_{B}$ would be determined just by the annihilation process $\psi_{B}\overline{\psi}_{B}\rightarrow\gamma^{\prime}\gamma^{\prime}$, whose thermally-averaged cross section in the $s$-wave limit is

[TABLE]

However, the process $\psi_{A}\overline{\psi}_{A}\rightarrow\psi_{B}\overline{\psi}_{B}$ is still active even after $\psi_{A}$ freezes out with a nearly constant $\psi_{A}$ abundance well above its equilibrium value, which can have impact on the relic abundance of $\psi_{B}$. Let $x_{f,B}=m_{B}/T_{f,B}$, $T_{f,B}$ being the temperature at $\psi_{B}$ freeze-out. As explained in App. B.1, when $\frac{\sigma_{B}}{\sigma_{A}}(\frac{m_{B}}{m_{A}})^{2}\gg(x_{f})^{2}$ (i.e. large $g^{\prime}$), a good approximation to the relic abundance $\Omega_{B}$ is

[TABLE]

This $\Omega\propto 1/\sqrt{\sigma}$ behavior is very different from the usual DM abundance relation $\Omega\propto 1/\sigma$. It arises because in this limit, there is a balance between depletion from $\psi_{B}$ annihilation and replenishment from $\psi_{A}\overline{\psi}_{A}\rightarrow\psi_{B}\overline{\psi}_{B}$ conversion. To our knowledge, this “balanced freeze-out” behavior has not been discussed before in the DM literature.

In Fig. 4.2a, we show numerical results for $\Omega_{B}$ as a function of $g^{\prime}$: for small $g^{\prime}$, $\psi_{B}$ freezes out in the standard way with $\Omega_{B}\propto 1/\sigma_{B}$, while for large $g^{\prime}$, $\Omega_{B}$ exhibits the $1/\sqrt{\sigma_{B}}$ scaling from balanced freeze-out. Thus, as long as $g^{\prime}$ is sufficiently large, then $\psi_{B}$ will be a subdominant DM component as desired. In Fig. 4.2b, we show the full solution to the coupled Boltzmann equations for $\psi_{A}$ and $\psi_{B}$ (see Eq. (B.1)) for the following benchmark scenario:

[TABLE]

where we have adjusted $\Lambda=250~{}\text{GeV}$ to yield the cross section $\langle\sigma_{A\overline{A}\rightarrow B\overline{B}}v\rangle=5\times 10^{-26}\text{cm}^{3}/\text{s}$ needed to achieve $\Omega_{A}\simeq\Omega_{\rm DM}\approx 0.2$. For this benchmark, $\psi_{B}$ has a much smaller abundance $\Omega_{B}\simeq 2.6\times 10^{-6}\,\Omega_{\rm DM}$. We have chosen the reference masses to be safe from existing constraints but visible with a reanalysis of existing Super-K data, and we have chosen the reference value of $g^{\prime}$ to be comparable to hypercharge in the SM. The values of $m_{\gamma^{\prime}}$ and $\epsilon$ are also interesting for explaining the muon $g-2$ anomaly [239, 240].

This model, though simple, exhibits a novel $\psi_{B}$ freeze-out behavior, and the “balancing condition” behind Eq. (4.12) may be interesting to study in other contexts. For much of parameter space of our interest in this chapter, the $\Omega_{B}\propto 1/\sqrt{\sigma_{B}}$ scaling affects the CMB and direct detection constraints on $\psi_{B}$. As discussed in Sec. 4.6, this scaling implies that the constraints from CMB heating on $\psi_{B}$ annihilation are largely independent of $g^{\prime}$. Similarly, unless there is some kind of inelastic splitting within the $\psi_{B}$ multiplet, there is a firm direct detection bound on $m_{B}$ that is also largely independent of $g^{\prime}$. Note that the benchmark scenario in Eq. (4.13) indeed satisfies these bounds (see the star in Fig. 4.9).

4.4 Detecting Boosted Dark Matter

With $\psi_{A}$ being the dominant DM species, the annihilation process $\psi_{A}\overline{\psi}_{A}\to\psi_{B}\overline{\psi}_{B}$ is active in the galactic halo today, producing boosted $\psi_{B}$ particles. To compute the flux of $\psi_{B}$ incident on the earth, we can recycle the standard formulas from indirect detection of WIMP DM. Roughly speaking, the (in)direct detection of boosted $\psi_{B}$ particles from $\psi_{A}$ annihilation is analogous to the familiar process of indirect detection of neutrinos from WIMP annihilation. For this reason, the natural experiments to detect boosted DM are those designed to detect astrophysical neutrinos. As we will see, $\psi_{B}$ typically needs to have stronger interactions with the SM than real neutrinos in order to give detectable signals in current/upcoming experiments.

We also want to comment that, due to the small mass and suppressed thermal abundance, the non-relativistic relic $\psi_{B}$ particles can be difficult to detect through conventional direct and indirect DM searches, even with efficient interaction between $\psi_{B}$ and SM states. (See Sec. 4.6 for existing bounds on $\psi_{B}$.) Therefore, detecting boosted $\psi_{B}$ particles may be the only smoking gun from this two-component $\psi_{A}$/$\psi_{B}$ system.

4.4.1 Flux of Boosted Dark Matter

The flux of $\psi_{B}$ from the GC is

[TABLE]

where $r_{\text{Sun}}=8.33~{}\text{kpc}$ is the distance from the sun to the GC and $\rho_{\text{local}}=0.3~{}\text{GeV}/\text{cm}^{3}$ is the local DM density. Since the $\psi_{A}\overline{\psi}_{A}\to\psi_{B}\overline{\psi}_{B}$ annihilation process yields two mono-energetic boosted $\psi_{B}$ particles with energy $m_{A}$, the differential energy spectrum is simply

[TABLE]

The quantity $J$ is a halo-shape-dependent dimensionless integral over the line of sight,

[TABLE]

where $s$ is the line-of-sight distance to the earth, the coordinate $r(s,\theta)=(r_{\text{Sun}}^{2}+s^{2}-2r_{\text{Sun}}s\cos{\theta})^{1/2}$ is centered on the GC, and $\theta$ is the angle between the line-of-sight direction and the earth/GC axis. Assuming the NFW halo profile [241], we use the interpolation functions $J(\theta)$ provided in Ref. [242] and integrate them over angular range of interest. In particular, when trying to mitigate neutrino backgrounds in Sec. 4.4.3, we will require the $\psi_{B}e^{-}\to\psi_{B}e^{-}$ process to give final state electrons within a cone of angle $\theta_{C}$ from the GC.

To illustrate the scaling of the flux, we integrate over a $10^{\circ}$ cone around the GC and obtain

[TABLE]

For completeness, the flux over the whole sky is:

[TABLE]

These estimates are subject to uncertainties on the DM profile; for example, an Einasto profile would increase the flux by an $\mathcal{O}(1)$ factor [242].

Note that this GC flux estimate is the same as for any mono-energetic DM annihilation products.555Up to factors of 2 if the particles considered are Majorana or Dirac, and the number of particles created in the final state. Therefore we can estimate the expected bound on the boosted DM-SM cross section by reinterpreting neutrino bounds on DM annihilation. Looking at Ref. [243], the anticipated Super-K limit on 1-100 GeV DM annihilating in the Milky Way exclusively to monochromatic neutrinos is $10^{-21}-10^{-22}\simeq\text{cm}^{3}/\text{sec}$. This is four to five orders of magnitude weaker than a typical thermal annihilation cross section ($\simeq 10^{-26}~{}\text{cm}^{3}/\text{sec}$). Assuming thermal relic $\psi_{A}$ DM exclusively annihilates to boosted $\psi_{B}$ particles, we can estimate the bound on the $\psi_{B}$-SM cross section by scaling down the charged current neutrino scattering cross section ($10^{-38}~{}\text{cm}^{2}$, see Eq. (4.32)) by the corresponding factor. This gives an estimated bound of

[TABLE]

which is consistent with the cross section derived later in Eq. (4.29) for a benchmark model that is on the edge of detectability.666Our numbers are less consistent with Super-K bounds shown in conference proceedings in Ref. [244], which are two orders of magnitude more constraining than expected from Ref. [245]. However, the details of the Super-K analysis are not available for direct comparison.

4.4.2 Detection of Boosted Dark Matter

The flux of boosted $\psi_{B}$ particles estimated from Eq. (4.17) is rather small.777For comparison, the flux of non-relativistic relic $\psi_{B}$ particles incident on earth is approximately

$\Phi_{\rm local}\simeq\frac{\rho_{\rm local}v_{\rm 0}}{m_{B}}\frac{\Omega_{B}}{\Omega_{\rm DM}}=2.25\times 10^{3}~{}\text{cm}^{-2}\text{s}^{-1}\left(\frac{200~{}\text{MeV}}{m_{B}}\right)\left(\frac{\Omega_{B}}{10^{-5}}\right).$

where $v_{0}\simeq$ 220 km/sec. Therefore, in order to detect boosted $\psi_{B}$, one needs a large volume, small background detector sensitive to the (quasi-)elastic scattering process

[TABLE]

where $X$ and $X^{\prime}$ are SM states (possibly the same). Because the $\gamma^{\prime}$ is kinetically-mixed with the photon, $\psi_{B}$ can scatter off any SM state $X$ with electromagnetic couplings via $t$-channel exchange of $\gamma^{\prime}$.888There are also subdominant scatterings from weak charges as well. A large scattering cross section favors light $m_{\gamma^{\prime}}$, large $\epsilon$, and large $g^{\prime}$; the values of $m_{\gamma}^{\prime}\gtrsim 10~{}\text{MeV}$ and $\epsilon\sim 10^{-3}$ in the benchmark in Eq. (4.13) are (marginally) consistent with current limits on dark photons [246].

Existing neutrino detectors such as Super-K, IceCube, and their upgrades can be employed to detect boosted DM via Eq. (4.20). The strategy is to detect Cherenkov light from the final state charged particles, so the energy of outgoing $X^{\prime}$ must be above the Cherenkov threshold. In terms of a Lorentz factor, the threshold is

[TABLE]

where there is typically a stricter analysis threshold $E^{\rm thresh}$ on $X^{\prime}$ as well, depending on experimental specifics. Furthermore, one needs to distinguish $\psi_{B}$ scattering from the large background of neutrino scattering events, which we discuss more in Sec. 4.4.3.

As shown in Fig. 4.3, there are three detection channels for boosted $\psi_{B}$ at a neutrino detector: elastic scattering off electrons, elastic scattering off protons (or nuclei), and deep inelastic scattering (DIS) off protons (or nuclei). As discussed in more detail in App. B.3, although the total $\psi_{B}$ scattering cross section off protons and nuclei can be sizable, the detectable signal strengths in these channels are suppressed relative to scattering off electrons.999The reason is that $\psi_{B}$ scattering proceeds via $t$-channel exchange of the light mediator $\gamma^{\prime}$, so the differential cross section peaks at small momentum transfers, while achieving Cherenkov radiation (or DIS scattering) requires large momentum transfers. For elastic scattering, this logic favors electrons over protons in two different ways: an $\mathcal{O}(1~{}\text{GeV})$ $\psi_{B}$ can more effectively transfer momentum to electrons compared to protons because of the heavier proton mass, and protons require a larger absolute momentum transfer to get above the Cherenkov threshold. Compounding these issues, protons have an additional form-factor suppression, identifying proton tracks is more challenging than identifying electron tracks [247, 248, 249], and the angular resolution for protons is worse than for electrons at these low energies [249]. We note that liquid Argon detectors are able to reconstruct hadronic final states using ionization instead of Cherenkov light, so they may be able to explore the (quasi-)elastic proton channels down to lower energies, even with smaller detector volumes [195, 196]. Thus, we focus on the elastic scattering off electrons

[TABLE]

as the most promising detection channel, though we present signal studies for the other channels in App. B.3. At detectors like Super-K, the signal would appear as single-ring electron events coming from the direction of the GC.

We start by discussing the kinematics of scattering off electrons (the same logic would hold for protons). In the rest frame of an electron target with mass $m_{e}$, the momenta of incoming and outgoing particles are:

[TABLE]

For $\psi_{B}$ coming from nearly-at-rest $\psi_{A}$ annihilation,

[TABLE]

The maximum scattered electron energy occurs when $\vec{p}$ and $\vec{p}^{\,\prime}$ are parallel:

[TABLE]

The minimum detectable energy is set by the analysis threshold (assumed to be above the Cherenkov threshold),

[TABLE]

Of course, to have any viable phase space, $E_{e}^{\rm max}\geq E_{e}^{\rm min}$. From Eqs. (4.25) and (4.26), we can also express the viable kinematic region in terms of boost factors $\gamma_{e}$ and $\gamma_{B}$ (taking $m_{A}\gg m_{B}\gg m_{e}$):

[TABLE]

The differential cross section for $\psi_{B}$ elastic scattering off electrons is:

[TABLE]

where $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2xz-2yz$, $s=m_{B}^{2}+m_{e}^{2}+2E_{B}m_{e}$, $t=q^{2}=2m_{e}(m_{e}-E_{e})$, and one should make the replacement $E_{B}=m_{A}$ for our scenario. To give a numerical sense of the Cherenkov electron signal cross section, integrating Eq. (4.28) over the allowed kinematic region for the benchmark in Eq. (4.13) yields

[TABLE]

for an experimental threshold of $E_{e}^{\rm thresh}=100~{}\text{MeV}$. The approximate scaling is derived in the limit $m_{e}E_{e}^{\rm thresh}\ll m_{\gamma^{\prime}}^{2}\ll m_{e}E_{e}^{\rm max}$, where the dependance on $E_{B}$, $m_{B}$, and $E_{e}^{\rm thresh}$ is weaker than polynomial, which holds in the vicinity of the benchmark point but not in general. For completeness, the full cross section for $\psi_{B}$-electron scattering without an energy threshold cut is

[TABLE]

Since this cross section is rather high, we have to account for the possibility that $\psi_{B}$ particles might be stopped as they pass through the earth. In Sec. 4.4.4, we find that the attenuation of $\psi_{B}$ particles is mild, so we will treat the earth as transparent to $\psi_{B}$ particles in our analysis.

In Fig. 4.4a, we show the normalized, logarithmic electron spectrum for different benchmarks, including the one from Eq. (4.13). The electron energy $E_{e}$ peaks at relatively low values due to the $t$-channel $\gamma^{\prime}$, as discussed further in footnote 9. We note that the position of the peak depends both on $m_{B}$ and $m_{\gamma^{\prime}}$, though the dominant effect of $m_{\gamma^{\prime}}$ is to change the overall signal cross section (not visible in this normalized plot). The angular distribution of the recoil electron is shown in Fig. 4.4b. The signal is very forward peaked, as expected from $m_{B}\gg m_{e}$. This is advantageous when looking for boosted DM from the GC, since the recoil electrons’ direction is tightly correlated to that of the $\psi_{B}$’s.

In Fig. 4.5, we compare the energy profile of the signal to the observed background electron events at SK-I [4]. Using the benchmark model in Eq. (4.13), we plot the (logarithmic) energy spectrum of the yearly signal event yield within a cone of $10^{\circ}$ around the GC:

[TABLE]

where $\Phi_{\text{GC}}^{10^{\circ}}$ is defined in Eq. (4.17), $\Delta T$ is a year, and $N_{\text{target}}$ is the number of targets (electrons) at Super-K. Anticipating the analysis of Sec. 4.5.2, we also plot a more realistic spectrum obtained by convolving the signal scattering cross section $\psi_{B}e^{-}\rightarrow\psi_{B}e^{-}$ with the shape of the DM halo. This convolved spectrum matches nicely to the naive spectrum from Eq. (4.31), as expected given the peaked nature of the angular spectrum in Fig. 4.4b, with signal losses at low energies arising because less energetic electrons can be more easily deflected outside the search cone. Once the background from Ref. [4] is scaled by the appropriate factor of $(\pi(10^{\circ})^{2})/(4\pi)\approx 8\times 10^{-3}$, the signal for this benchmark is visible above the background, though the peak location is (accidentally) at a similar location.

4.4.3 Backgrounds to Boosted Dark Matter

The major background to the boosted DM signal comes from atmospheric neutrinos, which are produced through interactions of cosmic rays with protons and nuclei in the earth’s atmosphere. Atmospheric neutrino energy spectrum peaks around 1 GeV and follows a power law $E^{-2.7}$ at higher energies [250]. The scattering process $\psi_{B}e^{-}\to\psi_{B}e^{-}$ with an energetic outgoing electron faces a large background from charged-current (CC) electron-neutrino scattering $\nu_{e}n\rightarrow e^{-}p$ when the outgoing proton is not detected, as well as $\overline{\nu}_{e}p\rightarrow e^{+}n$ since Cherenkov-based experiments cannot easily distinguish electrons from positrons. For $\mathcal{O}(1~{}\text{GeV})$ neutrinos, the CC cross sections are [251]

[TABLE]

While smaller than the expected signal cross section in Eq. (4.29), the atmospheric neutrino flux is much higher than the boosted $\psi_{B}$ flux. The neutral current process $\nu_{e}e^{-}\rightarrow\nu_{e}e^{-}$ can also mimic the signal but it is subdominant to the CC interaction due to $m_{e}/m_{p}$ suppression [251].

There are a number of discriminants one could use to (statistically) separate our signal from the neutrino background.

• •

Angular restriction: Boosted $\psi_{B}$ particles have a definite direction because they come from the GC. In galactic coordinates, the atmospheric neutrino background has no preferred direction. Therefore, one can impose that the detected electron falls within a cone of half-opening angle $\theta_{C}$ with respect to the GC. As shown in Fig. 4.6, there are two relevant axes to consider: the $z$-axis connecting the earth to the GC and the $z^{\prime}$-axis in the direction that the $\psi_{B}$ travels along. Through $\psi_{B}e^{-}\to\psi_{B}e^{-}$ scattering, a $\psi_{B}$ particle coming from an angle $\theta_{B}$ ($\theta^{\prime}_{B}=0$) will yield a final state electron with scattering angle $(\theta^{\prime}_{e},\phi^{\prime}_{e})$, with

[TABLE]

and $\phi^{\prime}_{e}$ uniformly distributes between [math] and $2\pi$. To the extent that the electron energy is large and $m_{A}\gg m_{B}\gg m_{e}$, we have $\cos\theta_{e}\approx\cos\theta_{B}$. As we will see in Sec. 4.5.4, the optimum angle $\theta_{C}$ to maximize the signal acceptance while minimizing the neutrino background is around $10^{\circ}$, assuming perfect angular resolution.

• •

Energy restriction: Boosted $\psi_{B}$ particles have a mono-energetic spectrum ($E_{B}=m_{A}$), compared to the continuous atmospheric neutrino energy spectrum. This implies a correlation between the measured $E_{e}$ and $\cos\theta_{e}$. That said, we suspect that the typical angular resolution of neutrino experiments is not fine enough to make use of this feature. In fact, more important than energy resolution is to have a low energy threshold, since as shown in Fig. 4.4, the signal cross section peaks at small $E_{e}$.

• •

Absence of muon excess: The process $\psi_{B}e^{-}\to\psi_{B}e^{-}$ does not have a corresponding muon signature, whereas the neutrino CC process $\nu_{e}n\rightarrow e^{-}p$ is always accompanied by $\nu_{\mu}n\rightarrow\mu^{-}p$. So an electron excess from boosted DM should not have a correlated excess in muon events. One can also require fully-contained events to reduce the cosmic ray muon background.

• •

Anti-neutrino discrimination: The anti-neutrino background $\overline{\nu}_{e}p\rightarrow e^{+}n$ is in principle reducible since it involves a final state positron instead of an electron. Super-K cannot perfectly distinguish $\nu_{e}$ from $\overline{\nu}_{e}$ events, but as of the SK-IV analyses, they have used likelihood methods to separate these two categories by studying the number of decay electrons in each process. The purity of the $\nu_{e}$ sample is $62.8\%$ and that of $\overline{\nu}_{e}$ is $36.7\%$ [252, 253]. We have not used this feature in our current analysis. Adding gadolinium to Super-K would help tagging neutrons from the $\overline{\nu}_{e}$ CC process, and thus might improve the purity of these samples [254].

• •

Multi-ring veto: The process $\psi_{B}e^{-}\to\psi_{B}e^{-}$ leads to electron-like single-ring events only, without correlated multi-ring events. In contrast, neutrino CC process $\nu_{e}n\rightarrow e^{-}p$ can lead to multi-ring events when the outgoing proton energy is above Cherenkov threshold [249], or when the scattering is inelastic so that other charged hadronic states such as $\pi^{\pm}$ are produced. Argon-based detectors could improve the background discrimination since they can detect the hadronic final states from neutrino scattering better than water-based experiments. We note that for some extreme parameters (increasing $g^{\prime}$ or $\epsilon$), it is possible for $\psi_{B}$ to interact twice (or more) within the detector, also creating a potential multi-ring signal (or a lightly-ionizing track in a scintillator detector). That said, for such high cross sections, the signal would be heavily attenuated while traversing the earth (see Sec. 4.4.4). Another potential disadvantage of a multi-ring veto is that we might miss out on interesting signals such as $\psi_{B}$ scattering accompanied by $\gamma^{\prime}$ bremsstrahlung ($\psi_{B}e^{-}\to\psi_{B}e^{-}\gamma^{\prime}\to\psi_{B}e^{-}e^{+}e^{-}$).

• •

Solar neutrino/muon veto: Solar neutrinos dominate the background under around 20 MeV [255], though one can of course preform an analysis in solar coordinates and exclude events from the sun. In addition, there is a background from muons that do not Cherenkov radiate but decay to neutrinos in the detector volume; these are relevant in the range of 30–50 MeV and can be mitigated through fiducial volume cuts [256]. To avoid both of these complications, we will use a cut of $E_{e}>100~{}\text{MeV}$ in our analysis below. Of course the threshold of 100 MeV in Super-K could be brought down as low as 50 MeV (where solar neutrino backgrounds start to dominate). The potential advantage of looking in the 50-100 MeV range is that the backgrounds from atmospheric neutrinos are lower. The main disadvantage is the degradation of the angular resolution of the detector [257].

The first two points favor detectors with excellent angular resolution and low energy thresholds on the outgoing electron. The next three points mean that one could distinguish the boosted DM signal from neutrinos coming from WIMP DM annihilation in the GC; boosted DM only gives a single-ring electron signal whereas neutrinos from WIMPs would give equal contributions to an electron and muon signal, both single- and multi-ring events, and equal contributions to a neutrino and anti-neutrino signal. The last point suggests the interesting possibility of looking for boosted DM from the sun due to DM solar capture, though in the particular model we study in this chapter, the solar capture rate is too small to be visible, and any boosted DM particles from the sun would face considerable solar attenuation (see Sec. 4.4.4 below). The above criteria can be thought of as a general algorithm for background rejection, while specifics can be tailored to a particular experiment. For instance, “multi-ring veto” does not apply to PINGU where Cherenkov rings cannot be reconstructed and all non-$\mu$-like events are classified as “cascade events”.

4.4.4 Impact of Earth Attenuation

As seen in Eq. (4.30), the signal cross section $\sigma_{Be^{-}\rightarrow Be^{-}}$ is relatively high, so as they cross the earth, the $\psi_{B}$ particles might get deflected and lose energy. This is a particularly important effect for Northern Hemisphere experiments like Super-K, where a typical $\psi_{B}$ would have to traverse through $\sim 10^{5}$ km (75% of the earth’s diameter). The dominant cause for energy loss is (minimum) ionization of atoms. While not relevant for detection, the main source of angular deflection is scattering off nuclei. In the following, we base our discussion on the standard analysis of particle propagation through matter as developed in the PDG [258].

First we estimate the $\psi_{B}$’s energy loss. Just as for a heavy charged particle traversing the earth (see, e.g., Ref. [259]), the main energy loss mechanism is through ionization. For $\beta\gamma$ factors of 10–100, a muon loses $\approx$ 1 GeV of energy per meter of rock [258]. A muon scatters off nuclei via a $t$-channel $\gamma$ exchange, while a $\psi_{B}$ scatters off nuclei via the exchange of a $\gamma^{\prime}$. We can approximate the length required for a $\psi_{B}$ to lose 1 GeV by scaling the couplings and the propagator of the $\psi_{B}$-$e^{-}$ scattering process to those of the $\mu$-$e^{-}$ scatterings:

[TABLE]

where $t=2m_{e}(m_{e}-E_{e})\approx-10^{-4}~{}\text{GeV}^{2}$ for our key benchmark in Eq. (4.13). In this case, $\psi_{B}$ loses $\approx 1~{}\text{GeV}$ per $9\times 10^{8}~{}\text{cm}$, giving rise to a total expected loss of 1 GeV per trip through the earth ($R_{\oplus}=6.4\times 10^{8}~{}\text{cm}$). Since 1 GeV of energy loss is never more than $\simeq 10\%$ of the $\psi_{B}$’s initial energy in the parameter space of interest, we will assume the earth is transparent to $\psi_{B}$’s for the rest of the analysis. Accounting for the energy loss is approximately equivalent to shifting the plots in Figs. 4.7 and 4.9 by the energy loss on the $m_{A}$ axis. The parameter space of small $m_{A}$ is the most affected, but that region is already constrained by CMB bounds as shown in Sec. 4.6.

Turning to the angular distribution, the dominant source of deflection is from elastic scattering off of nuclei. Note that $\psi_{B}$-$e^{-}$ scattering processes lead to very small angles of deflection because of the mass hierarchy $m_{B}\gg m_{e}$; indeed for the key benchmark in Eq. (4.13), the maximum possible deflection per scatter is 0.14∘. In contrast, Coulomb-like scattering of $\psi_{B}$ particles of nuclei can give rise to a more substantial deflection (including full reversal). The mean-square change in angle per collision process is

[TABLE]

where $\sigma_{BN\rightarrow BN}$ is the scattering cross section of $\psi_{B}$’s off a nucleus $N$ (see Eq. (B.28)), and $E_{N}^{\text{max}}$ is defined similarly to Eq. (4.25). In the last step of Eq. (4.37), we have inserted the benchmark value from Eq. (4.13). Treating the deflection of $\psi_{B}$ particles as a random walk through the earth, the total deflection is

[TABLE]

where the quantity under the square root is the number of steps (interactions). The mean free path to interact with a nucleus of charge number $Z$ and atomic number $A$ is

[TABLE]

where $n$ is the number density. Under the conservative assumption that the earth is entirely made of iron (the benchmark $A$ and $Z$ values above), and taking the mass density of earth to be $\rho=5.5\text{ g/cm}^{3}$, the number of scatters is $\approx 64$ for the benchmark in Eq. (4.13), giving a total deflection of:

[TABLE]

We checked that for different values of $m_{A}$, $m_{B}$, and $m_{\gamma^{\prime}}$, the total deflection does not vary much compared to Eq. (4.40). Since this deflection is small compared with the search cone of 10∘ that is used in Sec. 4.5.2, we neglect the angular deflection of $\psi_{B}$’s in our analysis.

Interestingly, if a signal of boosted DM is found, we could potentially use the earth attenuation to our advantage by correlating candidate signal events with the position of the GC with respect to the experiment. Indeed, with high enough statistics, the effect of earth shadowing would give rise to time-dependent rates, energies, and angles for $\psi_{B}$ scattering. As mentioned in footnote 2, solar attenuation would have an adverse effect on possible boosted DM signals from the sun. Since the radius of the sun is 100 times larger than that of the earth, the $\psi_{B}$ particles would lose a factor of 100 more energy, so we would need $m_{A}\gtrsim\mathcal{O}\text{(100 GeV-1 TeV)}$ for $\psi_{B}$ particle to escape the sun. Alternatively, for a smaller scattering cross section of $\psi_{B}$ particles with the SM, the sun might then be a viable source of signal [260].

4.5 Detection Prospects for Present and Future Experiments

We now assess the detection prospects for boosted DM at present and future detectors for neutrinos and/or proton decay. In Table 4.1, we summarize the (approximate) capacities/sensitivities of some of the representative relevant experiments, given in terms of the detector volume $V_{\rm exp}$, electron energy threshold $E_{e}^{\rm thresh}$, and angular resolution $\theta_{e}^{\rm res}$. From this table, we can already anticipate which experiments are going to be best suited for boosted DM detection.

Due to the relatively small flux of boosted DM, a larger volume detector, such as IceCube, KM3NeT, or ANTARES would be favored in order to catch more signal events. However, the energy threshold for the original IceCube are much too high for our purposes (and similarly for KM3NeT and ANTARES), since the energy transferred to the outgoing electron is suppressed due to the $t$-channel $\gamma^{\prime}$ (see Fig. 4.4). Even the $\simeq 1$ GeV threshold of PINGU is not ideal, though it will have some sensitivity.

So although Super-K/Hyper-K have smaller detector volumes, their low energy threshold is better matched to the boosted DM signal. In addition, Super-K/Hyper-K have excellent angular resolution,101010More accurately, the angular resolution of fully contained 1 ring Multi-GeV electrons is $1.5^{\circ}$ while that of fully contained 1 ring Sub-GeV electrons is less than $3.3^{\circ}$ as shown in Ref. [264]. which makes it possible to optimize the $\theta_{C}$ search cone criteria. Ultimately, MICA would offer better coverage in the energy range of interest. It is also worth mentioning that the proposed experiments for proton decay based on large scale liquid Argon detectors [195, 196] can also be sensitive to boosted DM due to their low thresholds and large volume. As mentioned in footnote 9, liquid Argon detectors may also have sensitivity to the proton scattering channel as well.

In the following subsections, we discuss event selection, signal/background rates, and expected signal significance in the above experiments. For signal-only studies of the subdominant channels involving $\psi_{B}$ scattering off protons/nuclei, see App. B.3.

4.5.1 Event Selection

As discussed in Secs. 4.4.2 and 4.4.3, the leading boosted DM signal comes from elastic scattering off electrons ($\psi_{B}e^{-}\to\psi_{B}e^{-}$) and the leading background is from atmospheric neutrinos (mostly $\nu_{e}n\rightarrow e^{-}p$). In principle, one could use the full multivariate information about the kinematics of the outgoing electron to separate signal and background. In order to keep the analysis simple, we will do a cut-and-count study to estimate the sensitivity.

To isolate events coming from the GC, we will use the search cone $\theta_{C}$ described in Fig. 4.6. The dominant background from CC $\nu_{e}$ scattering of atmospheric neutrinos is assumed to be uniform across the sky, so the background in a search cone of half-angle $\theta_{C}$ scales proportional to $\theta_{C}^{2}$. Of course, one cannot take $\theta_{C}$ to be too small, otherwise the signal acceptance degrades. To optimize for the signal significance in Sec. 4.5.4, we will convolve the angular dependence of halo $J$-factor and the angular dependence of the $\psi_{B}e^{-}\to\psi_{B}e^{-}$ cross section to figure out the optimum $\theta_{C}$. Anticipating that result, we will find

[TABLE]

where $10^{\circ}$ applies to the high resolution experiments (Super-K/Hyper-K), and the other experiments are limited by their angular resolutions.

From Eqs. (4.25) and (4.26), we have minimum and maximum electron energies $E_{e}^{\rm min}$ and $E_{e}^{\rm max}$ for the signal. Ideally, one would adjust the energy selection for a given value of $m_{A}$ and $m_{B}$, and try to push the analysis threshold $E_{e}^{\rm min}$ to be as low as possible. To be conservative, we will take the standard Super-K event categories for fully-contained single-ring electron events (see e.g. Ref. [253])

[TABLE]

without attempting to do finer energy binning. For Super-K, Hyper-K, and MICA, we will use both categories as separate event selections; for the Sub-GeV category we will choose only zero-decay events. PINGU has a higher energy threshold and cannot reconstruct Cherenkov rings nor efficiently separate $\mu$-like and $e$-like events near threshold, so we will only use the Multi-GeV category, while also adding in backgrounds from multi-ring events and $\mu$-like events; we will also scale the PINGU effective volume down by a factor of 2 to account for an estimated reconstruction efficiency of $\sim 50\%$ [191]. Note that the 100 MeV lower bound of the Sub-GeV category is above the nominal 10 MeV threshold of Super-K, so there is room for improved signal acceptance. Similarly, when the 1.33 GeV upper bound of the Sub-GeV category is above $E_{e}^{\rm max}$, then we are overestimating the background.

4.5.2 Signal Rates

Imposing the $\theta_{C}$ and energy range requirements, the number of signal electron events is:

[TABLE]

where $\Delta T$ is the time duration of the observation, $N_{\text{target}}$ is the number of target electrons, $\Phi_{\text{GC}}$ is the DM flux from the GC, and $\sigma_{Be^{-}\rightarrow Be^{-}}$ is the $\psi_{B}$-electron scattering cross section (which depends on the energy integration range in Eq. (4.28)). The factor of 10 in the second line is the number of electrons per molecule of water. The DM flux and scattering cross section have to be convolved in order to isolate events that pass the $\theta_{C}$ requirement, and the angles in the last line are the same as in Fig. 4.6 with $\theta_{e}$ given in Eq. (4.34). The integration limits $\theta^{\prime}_{e}\in\{\theta^{\prime}_{\rm min},\theta^{\prime}_{\rm max}\}$ are given by Eq. (4.35) by requiring $E_{e}\in\{E_{e}^{\rm max},E_{e}^{\rm min}\}$ (note the reversal of the limits, and that $\theta^{\prime}_{\rm min}=0$ if Eq. (4.25) is more restrictive than the energy categories above).

To get a sense of the expected signal rate, we consider the number of signal events for $\theta_{C}=10^{\circ}$ in the combined categories:

[TABLE]

broken down by $21.1/\text{year}$ for Sub-GeV and $4.0/\text{year}$ for Multi-GeV, and the reference cross sections are based on the benchmark in Eq. (4.13). In our analysis below, we will always assume that $\langle\sigma_{A\overline{A}\rightarrow B\overline{B}}v\rangle$ takes on the thermal relic reference value.

Because $\sigma_{Be^{-}\rightarrow Be^{-}}$ scales homogeneously with $g^{\prime}$ and $\epsilon$, the number of signal events does as well, so the only non-trivial dependence is on the mass parameters $m_{A}$, $m_{B}$, and $m_{\gamma^{\prime}}$. In Fig. 4.7, we set two benchmark values $m_{\gamma^{\prime}}=20~{}\text{MeV}$ and $m_{\gamma^{\prime}}=50~{}\text{MeV}$, and show what part of the $m_{A}-m_{B}$ parameter space yields

[TABLE]

for $x=1,10,100$. These reference values for $m_{\gamma^{\prime}}$ have been chosen such that the $t$-channel scattering processes are not overly suppressed by the dark photon mass, and the reference $\epsilon$ is close to the maximum allowed by dark photon constraints. In the triangular regions in Fig. 4.7, the top edge is set by $m_{A}$ which controls the DM number density (and therefore the annihilation rate), the left edge is set by the requirement that $m_{B}>m_{\gamma^{\prime}}$, and the diagonal edge is set by the electron energy threshold.

In these figures, we have included model-dependent constraints from CMB heating and direct detection, discussed in the later Sec. 4.6.111111The bump around $m_{B}=10~{}\text{GeV}$ in the CMB heating bound is due to a Sommerfeld resonance. It is worth emphasizing that both of these constraints are due to the thermal relic $\psi_{B}$ population, and are independent of the boosted DM phenomenon. Indeed, as discussed at the end of Sec. 4.2, we could give $\psi_{B}$ a small Majorana mass splitting, which would eliminate the bound from (elastic) direct direction experiments while not affecting very much the kinematics of boosted $\psi_{B}e^{-}\rightarrow\psi_{B}e^{-}$ detection. The CMB constraints are more robust since they mainly depend on $\psi_{B}$ being in thermal contact with the SM via $\psi_{B}\overline{\psi}_{B}\rightarrow\gamma^{\prime}\gamma^{\prime}$, though the CMB constraints could potentially softened if $\gamma^{\prime}$ somehow decays to neutrinos (or to non-SM states).

4.5.3 Background Rates

The atmospheric neutrino backgrounds have been measured by Super-K over a 10.7 year period, during runs SK-I (1489 days), SK-II (798 days), SK-III (518 days) and SK-IV (1096 days), and the final results are summarized in Ref. [253]. In the Sub-GeV category, a total of 7755 fully-contained single-ring zero-decay electron events were seen the 100 MeV to 1.33 GeV energy range, giving a yearly background rate of

[TABLE]

In the Multi-GeV category, 2105 fully-contained single-ring electron events were seen in the 1.33 GeV to 100 GeV energy range [253, 265], yielding

[TABLE]

To estimate the background for PINGU (which lacks the ability to reconstruct Cherenkov rings), we add in multi-ring and $\mu$-like events in the Multi-GeV category, changing $197~{}\text{year}^{-1}$ to $634~{}\text{year}^{-1}$, which then has to be scaled by the effective PINGU detector volume.

For the boosted DM search, the background is reduced by considering only events where the electron lies in the search cone $\theta_{C}$. We assume a uniform background distribution from the entire sky, so the background within a patch in the sky of angle $\theta_{C}$ is:

[TABLE]

For $\theta_{C}=10^{\circ}$ relevant for Super-K, we have

[TABLE]

Ideally, we would use the full energy dependence of the background in order to optimize the signal/background separation, but given the rather low background rate, we will make the conservative choice to consider the whole Sub-GeV energy range.

Since one can estimate the background by looking at a side-band away from the GC, the background uncertainties in a $\theta_{C}$ cone should be dominated by Poisson fluctuations. For the all sky background, we note that Super-K saw a $\simeq 10\%$ mismatch between the measured atmospheric background and the Monte Carlo estimate [265, 253], so there is in fact a bit of room beyond Poisson fluctuations to accommodate a boosted DM signal in the current Super-K data.121212Associated with the published search for DM from the GC via upward going muons [266], there is also unpublished electron data from SK-I, -II, and -III [267, 268]. For $\cos\theta>0.8$ ($\theta_{C}\simeq 37^{\circ}$), around 600 Sub-GeV fully-contained single-ring zero-decay electron events were observed in a 7.7 year period. This number has subsequently been updated to around 850 events in the full 10.7 year data set [269]. In principle, these could be used to set a stronger bound than we show in this chapter, since no statistically significant excess is seen.

In order to have a fair comparison of the sensitivities at different experiments, for Hyper-K and MICA, we use the same event selection assuming the same exposure time, based on the available Super-K $\sim 10$ year data set, and simply scale up the background rate proportional to the detector volume $V_{\text{exp}}$ (and adjust $\theta_{C}$ for MICA). As already mentioned, since PINGU has a higher energy threshold and an inability to reconstruct Cherenkov rings, we rescale the full Multi-GeV category (single-ring supplemented by the multi-ring and $\mu$-like events).

4.5.4 Estimated Experiment Reach

Given the signal and background rates above, we can find the optimal search cone $\theta_{C}$ to maximize the significance

[TABLE]

In Fig. 4.8, we plot the significance as a function of search angle for our benchmark model in Eq. (4.13); we checked that other parameter choices show similar behavior. We see that the significance peaks at around $10^{\circ}$, and falls off somewhat slowly after that. For Super-K/Hyper-K with $3^{\circ}$ resolution, we can effectively ignore experimental resolution effects and take $\theta_{C}$ at the optimal value. For PINGU and MICA, we approximate the effect of the experimental resolution by taking $\theta_{C}=\theta_{e}^{\rm res}$; a more sophisticated treatment would be to apply Gaussian smearing to the electrons. This is the logic behind Eq. (4.41) above.

In Figs. 4.9 and 4.10, we show the $2\sigma$ sensitively possible with the 10.7 years of Super-K data, using the optimal $\theta_{C}=10^{\circ}$ selection criteria, as well as the estimated reach for Hyper K, PINGU, and MICA for the same period of time. We treat the Sub-GeV and Multi-GeV categories separately and report the overall significance as the quadrature sum of the significances obtained from the two categories. We also show the current bounds from Super-K that one can place without the $\theta_{C}$ selection (i.e. using the all-sky background), taking $\delta N_{\text{bkgd}}/N_{\text{bkgd}}=10\%$ to account for systematic uncertainties in the all-sky background. Here, we are only allowing for the two energy categories, and further improvements are possible if one adjusts the energy range as a function of $m_{A}$ and $m_{B}$.

Taken together, these experiments have substantial reach for boosted DM. The prospects for Super-K to find single-ring electron events from the GC are particularly promising, given that the data (with angular information) is already available [270] and one simply needs to change from lab-coordinates to galactic coordinates (as in Refs. [267, 268]).

4.6 Summary of Existing Constraints

Apart from the measured neutrino fluxes discussed in Sec. 4.5.3, we know of no model-independent constraints on the boosted DM phenomenon. There are, however, constraints on the particular model described here, and we summarize those constraints in this section. The most relevant bounds are due mainly to the relic $\psi_{B}$ population left over from thermal freeze-out, which leads to bounds from “Direct detection of non-relativistic $\psi_{B}$” and “CMB constraints on $\psi_{B}$ annihilation” described below and seen in Figs. 4.9 and 4.10.

• •

Limits on the dark photon $\gamma^{\prime}$. As discussed earlier, dark photon searches have set limits of $m_{\gamma^{\prime}}\gtrsim\mathcal{O}(10~{}\text{MeV})$ and $\epsilon\lesssim 10^{-3}$, assuming the dominant decay mode is $\gamma^{\prime}\to e^{+}e^{-}$ [246]. For $m_{\gamma^{\prime}}<\mathcal{O}(100~{}\text{MeV})$, beam dump experiments place a bound of roughly $\epsilon\gtrsim 10^{-5}$ [271]. We have used $m_{\gamma^{\prime}}=20~{}\text{MeV}$ and $\epsilon=10^{-3}$ as a benchmark in this chapter, which yields a detectable boosted DM signal while satisfying the current dark photon bounds. Our benchmark is also within the region of interest for explaining the muon $g-2$ anomaly [240, 239].

• •

Direct detection of non-relativistic $\psi_{A}$. Thermal relic $\psi_{A}$ particles are subject to constraints from conventional DM direct detection experiments (e.g. XENON, LUX, and CDMS) via their scattering off nuclei. As discussed in more detail in App. B.2, the constraints on $\psi_{A}$ are rather weak since $\psi_{A}$ has no tree-level interactions with the SM. That said, $\psi_{A}$ can scatter off nuclear targets via a $\psi_{B}$-loop. Since we have approximated the $\psi_{A}\overline{\psi}_{A}\psi_{B}\overline{\psi}_{B}$ interaction as a contact operator, this loop process is model-dependent. In Fig. B.1a and Eq. (B.15), we give an example UV completion involving an extra scalar $\phi$ that allows us to estimate the $\psi_{A}$-nucleon scattering cross section. Due to the loop factor and the mass suppression from $m_{\phi}\gg m_{A}$, the limits on $\psi_{A}$ are safe for most values of the parameter space, as shown in Fig. B.1b. As already mentioned, one could introduce inelastic splitting within the $\psi_{A}$ multiplet to further soften direct detection constraints [234, 235, 236].

• •

Direct detection of non-relativistic $\psi_{B}$. Despite the small relic abundance of $\psi_{B}$, it has a large $\psi_{B}$-nucleon scattering cross section, as calculated in App. B.3.1.

[TABLE]

where the scaling assumes $m_{B}\ll m_{p}$. Thus, direct detection experiments essentially rule out any elastic $\psi_{B}$-nucleon scattering above the detector threshold. Of course, in the parameter space of our interest, the $\psi_{B}$ mass is $\leq\mathcal{O}(1~{}\text{GeV})$, which is close to or below the threshold of LUX [272] and the low CDMS threshold analysis [273], and the most constraining limits come from CDMSLite [274] and DAMIC [275]. Because of this, light $\psi_{B}$ particles can evade existing direct detection bounds.

In Fig. 4.10, we demonstrate the constraints on the $(g^{\prime},m_{B})$ plane from the DAMIC experiment (which has a lower threshold than CDMSLite), using the effective nuclear cross section

[TABLE]

Essentially, the allowed parameter space is independent of $g^{\prime}$ and $(m_{\gamma^{\prime}})^{-4}$, since the expected $\psi_{B}p\rightarrow\psi_{B}p$ cross section is so large that any events above the energy threshold of the experiment would be seen. There is also the fact that when $g^{\prime}$ is $\mathcal{O}(10^{-2})$ and higher, the abundance scales as $g^{\prime-2}$ (see Eq. (4.12)), which cancels with the $g^{\prime 2}$ scaling of $\sigma_{Bp\rightarrow Bp}$, yielding a $g^{\prime}$-independent bound.131313The $\psi_{B}p\rightarrow\psi_{B}p$ cross section scales as $(m_{\gamma^{\prime}})^{-4}$, but we have checked that these bounds do not soften until $m_{\gamma^{\prime}}$ is higher than $\mathcal{O}(1~{}\text{GeV})$, which is not the regime we are studying in this chapter. Of course, as with $\psi_{A}$, the direct detection bound on $\psi_{B}$ could be alleviated by introducing inelastic splittings.

It has been recently pointed out that sub-GeV DM might be better constrained by scattering off electrons rather than off nuclei [276], as in recent XENON10 bounds [277]. In our case, these bounds are subsumed by CMB heating bounds discussed below. Note that for $\psi_{B}e^{-}\rightarrow\psi_{B}e^{-}$, the conventional direct detection process and the boosted DM detection process have very different kinematics, so one should not be surprised that the XENON10 bounds do not influence the boosted DM signal regions.

• •

Indirect detection of non-relativistic $\psi_{B}$. The annihilation process $\psi_{B}\overline{\psi}_{B}\to\gamma^{\prime}\gamma^{\prime}$ and the subsequent $\gamma^{\prime}$ decay to two $e^{+}/e^{-}$ pairs gives rise to a potential indirect detection signal in the positron and diffuse $\gamma$-ray channels. The recent constraint on DM annihilation in positron channel from AMS-02 is demonstrated in Refs. [278, 279], where the bound is strongest for 2-body final state and weaker when there are more particles in the final state like in our case. The suppressed relic abundance of $\Omega_{B}$ relative to $\Omega_{\rm DM}$ helps relieve the constraints on our model. In addition, at the sub-GeV mass which we are interested in, the background uncertainty of the above indirect detection limit is large due to solar modulation. The CMB considerations below give stronger constraints for the parameter range of our interest. The diffuse $\gamma$-ray signal from e.g. inverse Compton scattering of $e^{\pm}$ produced from $\psi_{B}$ annihilation has a smaller cross section and also faces large background uncertainty in the sub-GeV region. In fact, the $\gamma$-ray search for DM at Fermi, for instance, has a lower energy cutoff at $\sim 4$ GeV [280]. Indirect detection signals from $\psi_{A}$ annihilation have to go through higher order or loop processes, and are much suppressed.

• •

CMB constraints on $\psi_{B}$ annihilation. With a light mass of $m_{B}\lesssim\mathcal{O}(1~{}\text{GeV})$, thermal $\psi_{B}$ annihilation in the early universe may be subject to bounds from CMB heating [281]. The CMB constrains the total power injected by DM into ionization, heating, and excitations. For the dominant DM component with relic density $\Omega_{\rm DM}\approx 0.2$, the bound is directly imposed on the quantity:

[TABLE]

where $f_{\text{eff}}$ is the fraction of the annihilation power that goes into ionization, which depends on the annihilation channel and its energy scale. Though $\psi_{B}$ is a small fraction of total DM, it does annihilate into $\gamma^{\prime}$ which subsequently decays via $\gamma^{\prime}\to e^{+}e^{-}$. Therefore, the CMB spectrum constrains

[TABLE]

where the last relation is obtained using Eq. (4.12) for $\Omega_{B}/\Omega_{A}$, which is valid for large values of $g^{\prime}$ (typically for $g^{\prime}\gtrsim 0.1$) as explained in the App. B.1. These limits are illustrated in Fig. 4.10 for $f_{\text{eff}}=1$, which is a conservative assumption. Due to the presence of a light $\gamma^{\prime}$, there can be an extra Sommerfeld enhancement factor to the $\langle\sigma_{B\overline{B}\rightarrow\gamma^{\prime}\gamma^{\prime}}v\rangle$ in Eq. (4.56). For the parameter space we consider, we expect that this enhancement saturates at CMB time, which leads to an extra factor of [282]

[TABLE]

This enhancement contributes at low velocities, so we do not expect it to change the picture at freeze out, but it would be relevant in the CMB era where $v\approx 10^{-3}$. For our current parameter space, $S\approx 1$ until high values of $g^{\prime}=1$ where it becomes $\mathcal{O}(10)$. We incorporate the enhancement in the calculation of our CMB limits, as can be seen from the resonance peaks in Fig. 4.10.

• •

BBN constraints on $\psi_{B}$ annihilation. The energy injection from $\psi_{B}$ annihilation in the early universe can also alter standard BBN predictions [283, 284]. The constraints from hadronic final states are the most stringent, comparable to or even somewhat stronger at $\mathcal{O}(1~{}\text{GeV})$ than those from the CMB heating as discussed above [283]. However, as we focus on $m_{\gamma^{\prime}}$ of $\mathcal{O}(10~{}\text{MeV})$, the production of hadronic final states ($n,p,\pi$) from the leading annihilation channel $\psi_{B}\bar{\psi}_{B}\rightarrow\gamma^{\prime}\gamma^{\prime}$ followed by $\gamma^{\prime}$ decay are not kinematically possible. The subleading channel $\psi_{B}\overline{\psi}_{B}\rightarrow q\bar{q}$ is $\epsilon^{2}$ suppressed. Thus, the major energy injection to BBN is mostly electromagnetic from $\gamma^{\prime}\rightarrow e^{+}e^{-}$, and the related constraint in this case are much weaker than the CMB bound we have considered above [283].

• •

Dark matter searches at colliders. By crossing the Feynman diagrams in Fig. 4.3, we see that $\psi_{B}$ can be produced at colliders such as LEP, Tevatron, and the LHC. If $\psi_{B}$ were to interact with SM electrons or quarks via a heavy mediator, then collider searches would provide a stronger bound than direct detection at these low DM masses. However, this complementarity is lost when the interaction is due to a light mediator [285, 286, 287], which applies to our case where $\psi_{B}$ interacts with SM states via an $\mathcal{O}(10~{}\text{MeV})$ dark photon. In addition, compared to the irreducible main background from electroweak processes, e.g. $e^{+}e^{-}\rightarrow Z^{(*)}\rightarrow\nu\bar{\nu}$, the production cross section of $\psi_{B}$ is suppressed by $\epsilon^{2}\lesssim 10^{-6}$, so the collider constraints on our model are rather weak.

4.7 Conclusions and Other Possibilities

In this chapter, we presented a novel DM scenario which incorporates the successful paradigm of WIMP thermal freeze-out, yet evades stringent constraints from direct and indirect detection experiments, and predicts a novel signal involving boosted DM. The example model features two DM components, $\psi_{A}$ and $\psi_{B}$. The heavier particle $\psi_{A}$ (which is the dominant DM component) experiences assisted freeze-out [188] by annihilating into the lighter particle $\psi_{B}$ (which is the subdominant DM component). The whole dark sector is kept in thermal contact with the SM in the early universe via kinetic-mixing of a dark photon with the SM photon. Only $\psi_{B}$ couples directly to the dark photon (and hence to the SM), so the dominant DM component $\psi_{A}$ can largely evade current DM detection bounds. If such a scenario were realized in nature, then the leading non-gravitational signal of DM would come from annihilating $\psi_{A}$ particles in the galactic halo producing boosted $\psi_{B}$ particles that could be detected on earth via neutral-current-like scattering via the dark photon. In large volume neutrino or proton-decay detectors, the smoking gun for this scenario would be an electron signal pointing toward the GC, with no corresponding excess in the muon channel. Liquid argon detectors could potentially detect boosted DM through (quasi-)elastic proton scattering, as well as improve the rejection of the dominant neutrino CC backgrounds by vetoing on hadronic activity. Future experiments that use LArTPC technology for tracing the particle paths [288, 289] will provide both directionality and better background discrimination.

This phenomenon of boosted DM is generic in scenarios with multiple DM components. In fact, models with a single component DM could also potentially give rise to the same signature. If the stabilization symmetry is $\mathbb{Z}_{3}$, then the semi-annihilation process $\psi\psi\to\overline{\psi}\phi$ (where $\phi$ is a non-DM state) is allowed [290, 291, 292]. For $m_{\phi}=0$, the outgoing $\overline{\psi}$ would have energy $E_{\psi}=(5/4)m_{\psi}$. In the limit $m_{\psi}\gg m_{e}$, $\gamma_{\psi}=1.25$ implies a maximum $\gamma_{e}^{\rm max}=2\gamma_{\psi}^{2}-1=2.125$, which is above the Cherenkov threshold in water (and ice). Of course, the $\mathbb{Z}_{3}$ symmetry is not consistent with $\psi$ being charged under a $U(1)^{\prime}$, so additional model building would be necessary to get a sufficiently large scattering with the SM. But this example shows why non-minimal dark sectors tend to have some production cross section for boosted DM.

It is intriguing to consider other scenarios where DM mostly annihilates to other stable states in the dark sector. For example, if both $\psi_{A}$ and $\psi_{B}$ are charged under the $U(1)^{\prime}$ and the mass hierarchy is

[TABLE]

then the annihilation $\psi_{A}\overline{\psi}_{A}\to\gamma^{\prime}\gamma^{\prime}$ would be followed by the decay $\gamma^{\prime}\to\psi_{B}\overline{\psi}_{B}$, and the boosted $\psi_{B}$ particles could again be detected via $t$-channel $\gamma^{\prime}$ exchange with the SM. Of course, now $\psi_{A}$ itself has tree-level $\gamma^{\prime}$ exchange diagrams with the SM, but if $\psi_{A}$ has a Majorana mass splitting (allowing it to evade direct detection bounds), boosted DM would again be the dominant mode for DM discovery.141414There would also be interesting signals for $\psi_{B}$ in DM production/detection experiments [293].

The above scenario is particularly interesting in light of the gamma ray excess recently seen in the GC [38]. In the context of DM, this signal could be explained via cascade decays $\psi_{A}\overline{\psi}_{A}\to\gamma^{\prime}\gamma^{\prime}$ followed by $\gamma^{\prime}\to\text{SM}\,\text{SM}$ [45, 43, 46, 44]. Boosted DM could be produced in the same cascade process, since the dark photon could easily have comparable branching ratios for $\gamma^{\prime}\to\text{SM}\,\text{SM}$ and $\gamma^{\prime}\to\psi_{B}\overline{\psi}_{B}$ when $m_{\gamma^{\prime}}>m_{B}$. More generally, it is interesting to contemplate scenarios where $\psi_{A}$ partially annihilates to boosted $\psi_{B}$ and partially to SM states. For example, the bremsstrahlung process of $\psi_{A}\psi_{A}\rightarrow\psi_{B}\psi_{B}\gamma^{\prime}$, where the $\gamma^{\prime}$ decays to an electron-positron pair, can be a source of positrons that can be detected in experiments like AMS-02 [294] or indirectly in Gamma ray telescopes [295]. Of course, if the $\psi_{B}$ states are not too depleted, then they could give indirect detection signals of their own.

Finally, it is worth considering the broader experimental signatures possible in the paradigm of DM annihilating to stable dark sector states [214, 296, 225, 297, 217], with simple extensions/variations based on our current model. If $m_{B}\ll m_{e}$, then $\psi_{B}$ acts effectively like dark radiation, which may leave signatures in CMB observables such as $N_{\rm eff}$ [217]. If $\psi_{A}$ has a non-negligible solar capture cross section, then boosted DM could emerge from the sun. If $\psi_{B}$ takes up sizable fraction of the total DM abundance (perhaps via a leading asymmetric component), then the fact that $\psi_{B}$ has strong self-interactions may have implications for small scale structure of DM halos including the known anomalies such as cusp-core and too-big-to-fail problems [298, 299]. The potentially rich structure of the dark sector motivates a comprehensive approach to DM searches.

Chapter 5 Extending Boosted Dark Matter to DUNE

5.1 Introduction

Gravitational evidence for dark matter (DM) is overwhelming [10, 12, 17], but all nongravitational means of DM detection have not yet resulted in a definitive discovery. It is therefore essential to expand DM searches to encompass as many possible DM signals. Previous work [5] has proposed a new class of DM models called boosted dark matter (BDM) with novel experimental signatures at neutrino experiments. BDM search strategies are complementary to existing indirect detection searches for DM at neutrino detectors.

BDM expands the weakly interacting massive particle (WIMP) paradigm to a multicomponent dark sector that includes a component with a large Lorentz boost obtained today due to decay or annihilation of another dark particle at a location dense with DM. In this class of models, the boosted component can scatter off standard model (SM) particles similarly to neutrinos, and can thus be detected at neutrino experiments. Various extensions built on the BDM model [260, 300, 301, 302] have studied the potential reach at large volume neutrino detectors and even direct detection experiments.

In this paper, we present BDM searches assuming a constant scattering amplitude, which highlight the reach of different neutrino technologies with different experimental features, and in particular electron energy thresholds. Focusing on scenarios in which BDM scatters off electrons (and leaving scattering off protons to future work [303]) the scattering process of interest, shown in Fig. 5.1, is

[TABLE]

where $B$ is a subdominant DM component with a Lorentz boost due to the annihilation of another heavier dominant state $A$

[TABLE]

as shown in Fig. 5.2.

We present the potential reach for two searches for BDM, one where the boosted particle $B$ originates at the galactic Center (GC) and one where $B$ originates at dwarf galaxies (dSphs). Although dSphs are a great source for DM since they are low in astrophysical backgrounds, their DM density is lower than that of the GC, so we perform a stacked analysis to increase statistics and improve sensitivity.

We take advantage of $B$’s large Lorentz boost in reducing background as the emitted electrons scatter in the forward direction and therefore point to the origin of the BDM particle. This is different from the omnidirectional atmospheric neutrino background, dominated by the charged current processes

[TABLE]

Experiments of particular interest are Cherenkov detectors like Super-Kamiokande (Super-K) [61] and Hyper-Kamiokande (Hyper-K) [304], and liquid argon time projection chambers (LArTPCs) like the upcoming Deep Underground Neutrino Experiment (DUNE) [65]. Argon-based detectors utilize a new technology that has not previously been thoroughly investigated within the context of DM searches. We explore LArTPCs’ excellent angular resolution and particle identification in this paper, and emphasize the discrimination power of LArTPC experiments even with smaller volumes than their Cherenkov counterparts. We show the overall sensitivity of Super-K, Hyper-K and DUNE in setting limits on the DM-SM scattering cross section for the case of annihilation from another heavier component $A$. The decay case can be worked out in a similar fashion.

The rest of this paper is organized as follows: In Sec. 5.2, we introduce a simplified parametrization that captures BDM’s main features, and set up the framework to relate the expected number of detected events to the general properties of BDM. We then study event selection in Sec. 5.3 and background rejection in Sec. 5.4 for the Cherenkov and argon-based technologies. We finally show the experimental reach at current and future neutrino experiments to BDM originating in the GC in Sec. 5.5 and in dSphs in Sec. 5.6, and conclude in Sec. 5.7.

5.2 Boosted Dark Matter

5.2.1 Features of Boosted Dark Matter

One of the most studied paradigms of DM is that of WIMPs in which DM is a single cold thermal particle that froze out early in the Universe’s history. Various detection methods have been used to search for WIMP DM: direct detection in which nonrelativistic DM particles scatter off heavy nuclei [305, 272, 274, 85], and indirect detection in which SM particles resulting from DM annihilation/decay are detected (see for example, [278, 295, 38, 8]). Indirect detection signals originate in DM-dense regions, two of which are the GC and dSphs.

BDM is a class of multicomponent models in which a component of the dark sector has acquired a Lorentz boost today. Let the DM sector be composed of a dominant component $A$ and a subdominant component $B$.111$A$ and $B$ can be the same particle as in the case of a $Z_{3}$ symmetry for example [306, 307, 308, 309, 310, 198], and $A$ can correspond to more than one particle in the case of a more complex dark sector.

• •

The particle $B$ is boosted due to either annihilation or decay of a second state $A$, as shown in Fig. 5.2. Other processes that would boost the $B$ particle can be easily derived from the subsequent formalism, such as semiannihilation $AA\rightarrow B\phi$ [198] for example, with the energy of $B$ satisfying $E_{B}\gg m_{B}$.

• •

The boosted particle $B$ interacts with the SM through a scattering process. In this work, we focus on the case of $B$ scattering off electrons $Be^{-}\rightarrow Be^{-}$, as in Ref. [5]. We leave the case of $B$ scattering off protons [260, 300] to future work [303].222Proton scattering is more important for scenarios where DM, and in this case $A$, is captured in the Sun. This case depends on the capture scenario rather than the initial DM density, and therefore it is not incorporated in this work.

Searching for BDM therefore involves a hybrid approach, as one would directly detect the $B$ particle scattering off SM particles, and at the same time indirectly detect the $A$ component. In the following , we present a simplified parametrization of BDM in order to compare the reach of different neutrino detector technologies.

5.2.2 Flux of Boosted Dark Matter from Annihilation

The flux of $B$ produced in $A$ annihilation (see Fig. 5.2) within a region of interest (ROI) of a particular source is

[TABLE]

The annihilation J-factor $j_{\text{ann}}$ is obtained by integrating over the DM density squared along the line of sight at a particular position in the sky,

[TABLE]

The thermally averaged cross section $\langle\sigma_{\overline{A}A\rightarrow\overline{B}B}v\rangle$ is the annihilation cross section of the process that produces the $B$ particles, taken as a reference to be equal to the thermal cross section $\langle\sigma_{\overline{A}A\rightarrow\overline{B}B}v\rangle=3\times 10^{-26}\text{cm}^{3}/\text{sec}$. Any deviation is an overall rescaling of the flux.

As was previously argued in Ref. [5], the optimal choice of ROI for the GC analysis is $\approx 10^{\circ}$ around the GC for the case of annihilation.333The value of the optimal opening angle for decay ($A\rightarrow B\overline{B}$) cannot be taken as $10^{\circ}$ without a proper analysis. The initial value of the opening angle depends largely on the DM distribution. The fact that annihilation signals scale as the DM density squared while decay signals scale linearly with DM density means that DM will be less localized near the center, and that leads to a larger optimal choice of ROI. We therefore adopt the same ROI in this analysis. We define $J_{\text{ann}}$ as the integrated J-factor $j_{\text{ann}}(\Omega)$ over a patch of the sky, assuming an NFW profile [241], as

[TABLE]

where the numerical value corresponds to a cone of half angle $10^{\circ}$ around the GC [242].

The spectrum of $B$ is $dN_{B}/dE_{B}$, which in the case of the $A\overline{A}\rightarrow B\overline{B}$ process is

[TABLE]

Therefore, the integrated flux over a patch of the sky is

[TABLE]

The numerical values of the flux of DM integrated over the whole sky and over a cone of half angle $10^{\circ}$ for $\overline{A}A\rightarrow\overline{B}B$ are

[TABLE]

5.2.3 Implications of Forward Scattering

In the energy range of $\mathcal{O}(10\text{ MeV})-\mathcal{O}(100~{}\text{GeV})$, the dominant background for any neutrinolike signal is atmospheric neutrinos [255, 311, 253].444Solar neutrinos dominate below energies of $30$ MeV. Although we know the location of the Sun and can thereby veto solar neutrinos, we avoid this parameter space in order to be conservative as it is hard to estimate the ability of photomultipliers to trigger on events with such low energies. The key aspect in discriminating the background, which is omnidirectional, from the signal, which originates at a location dense in DM, is adopting a search cone strategy. As shown in Fig. 5.3, we veto all electrons that are emitted at an angle larger than $\theta_{C}$ around a particular source. This strategy takes advantage of forward scattering of the electron, emitted in the same direction as the incoming $B$.

As was computed in Ref. [5], the expected number of electron events $N_{\text{signal}}^{\theta_{C}}$ is obtained by convolving the initial DM distribution over the electron scattering angle of the $Be^{-}\rightarrow Be^{-}$ process, such that the emitted electron is scattered at angles smaller than $\theta_{C}$ around a particular source.

[TABLE]

where $\Delta T$ is the exposure time, and $N_{\text{target}}$ is the number of target electrons in the experiment considered. The angle $\theta_{B}$ is the polar angle of $B$ with respect to the source (GC or dSphs). The angle $\theta^{\prime}_{e}$ is the polar angle of $e^{-}$ with respect to the incoming direction of the $B$ (see Fig. 5.3). $f_{B}(\theta_{B})$ is the flux of the incoming $B$ particles as a function of the polar angle, integrated over the azimuthal angle. For a particular source, the total flux is related to $f_{B}$ by

[TABLE]

This is equal to Eq. (5.11) when $\alpha=10^{\circ}$.

As we show in App. B.4, in the limit where the energy of the BDM particle is much higher than the electron mass ($E_{B}\gg m_{e}$), highly boosted DM (with a Lorentz boost factor $\gamma_{B}\gg 1$) scatters off electrons which are then emitted in the forward direction ($\theta_{e}^{\prime}=0$). We can therefore use the electron scattering angle to infer the BDM’s origin. In this limit, the convolution of Eq. (5.12) can be simplified as

[TABLE]

It is important to note that the cross section $\sigma_{Be^{-}\rightarrow Be^{-}}^{\text{measured}}$, hereafter labeled $\mathcal{I}$, is not the total cross section, but rather the measured one, as the energy threshold of the experiment introduces an energy cutoff.

We write the measured cross section $\mathcal{I}$ as a function of the energy threshold $E_{\text{thresh}}$ in order to facilitate the comparison among experiments with different characteristics. Assuming that the limiting experimental factor is the energy threshold rather than the angular resolution, and this is a good approximation that follows from scatterings being in the forward direction and the excellent angular resolution of neutrino experiments, we write the measured cross section as a function of the measured energy of the emitted electron $E_{e}$

[TABLE]

The upper limit of integration is

[TABLE]

which is the maximum allowed by the kinematics of the scattering process.

5.2.4 Constant Amplitude Limit

In order to compare the reach of different experiments, we extract the dependence on the energy threshold while assuming a constant scattering amplitude. This simplifies the parameter space in order to better illustrate the reach of different experiments.

Let $\sigma_{0}$ be the total cross section for the process $Be^{-}\rightarrow Be^{-}$,

[TABLE]

If we assume a flat amplitude $|\mathcal{M}|^{2}=$ constant, we can then relate $\mathcal{I}$ defined in Eq. (5.15) with $\sigma_{0}$ defined in Eq. (5.17) by

[TABLE]

Below, we estimate limits on the quantity $\sigma_{0}$. The expected number of events given by Eq. (5.14) is

[TABLE]

5.3 Event Selection

The backgrounds to the signal process $Be^{-}\rightarrow Be^{-}$ are all processes in which an electron in the appropriate energy range is emitted from neutrino-induced scatterings. The processes with the highest cross sections are charged current neutrino scatterings $\nu_{e}+n\rightarrow e^{-}+p$ and $\overline{\nu_{e}}+p\rightarrow e^{+}+n$. For the energies of interest in $\mathcal{O}(10~{}\text{MeV})-\mathcal{O}(100~{}\text{GeV})$, the dominant background is atmospheric neutrinos. Neutrinos scattering in detectors produce both electrons and muons while the signal is present only in electron events. Therefore, an important feature of this BDM model is an excess in the electron channel over the muon channel. We now study the features of the signal that are used to discriminate against the background in Cherenkov and LArTPCs detectors separately.

5.3.1 Cherenkov Detectors: Super-K

We study Super-K as an example of Cherenkov detectors in this analysis. Super-Kamiokande is a large underground water Cherenkov detector, with a fiducial volume of 22.5 kton of ultrapure water. It has collected over 10 years of atmospheric data, which would be the target data set for this analysis [312, 4, 265].

The atmospheric neutrino backgrounds, as well as signal events in Super-K, are single-ring electrons, detected with the following properties.

• •

Energy range: for the electron to be detected in a Cherenkov experiment, the electron energy $E_{e}$ has to be above the Cherenkov limit $\gamma_{\text{water}}m_{e}$, with $\gamma_{\text{water}}=1.51$. The experimental threshold for the atmospheric neutrino analysis is, however, $E_{\text{thresh}}=100$ MeV, which is higher than $\gamma_{\text{water}}m_{e}$ and it is what sets the threshold on the electron detectability. This energy threshold is set such as to avoid Michel electrons which are the electrons produced in muon decay [313].

• •

Directionality: As we have previously argued, signal electrons are emitted in the forward direction, and therefore are a good tool to point at the origin of BDM. The angular resolution of Super-K improves as a function of the electron energy up to a point where all photomultipliers saturate, in which case it degrades and it gets harder to infer the direction of the electron. We therefore take a conservative value of the angular resolution as $5^{\circ}$ across all energies studied ($E_{e}\in[100\text{ MeV}-100\text{ GeV}]$). A more detailed study is required by the Super-K collaboration to find the appropriate resolution for this analysis. This conservative resolution is smaller than the full extent of the GC in the sky, so it will not impact the results. For the dSphs searches, we only trigger on electrons within $5^{\circ}$ from a particular source location.

• •

Gadolinium: Gadolinium has one of the highest neutron capture rates. Tests have been conducted for its use in Super-K. When added to Super-K, gadolinium captures emitted neutrons in the $\overline{\nu_{e}}+p\rightarrow e^{+}+n$ process and emits a distinctive 8 MeV photon, and therefore triggers on the $\overline{\nu_{e}}$ background [254, 314, 315, 316, 317, 318]. A full Super-K study will be able to estimate the reduction in background events when gadolinium is used, but it will not be included in this analysis.

Hyper-K is the future Super-K upgrade but with 25 times the fiducial mass,555The Hyper-K detector design might be modified for greater photomultiplier coverage and smaller mass [324], but we assume the volume used in the initial letter of intent for this study [189]. and thus will improve the sensitivity of Cherenkov detectors to BDM. In the following we assume it has the same properties as Super-K, from angular resolution to energy threshold [189, 63, 319, 320, 304, 321].

5.3.2 Argon-Based Detectors: DUNE

We now turn to the event selection at DUNE. DUNE is a planned LArTPC experiment which will be located at the Sanford Underground Research Lab. It will serve as the far detector for the long baseline neutrino facility and will be performing off-beam physics. It will include four 10-kton detectors. In the following, we study the sensitivity of 10 and 40 kton volume experiment to BDM [65]. The BDM features that we use to select potential signal events are the following:

• •

Energy range: To avoid being overwhelmed by the solar neutrino background, and to be conservative with the capability of the photodetector system to trigger on these events, we focus on the emitted electrons of energies $E_{e}>30~{}\text{MeV}$. This is a factor of 3 lower than a similar analysis at Super-K. Unlike Cherenkov detectors, Michel electrons are clearly associated with the parent muon track in LArTPCs. It is therefore easy to distinguish Michel electrons from electrons produced in charged current scatterings, and thus, the energy threshold can be lowered from 100 to 30 MeV.

• •

Absence of hadronic processes: The signal does not include any hadrons in the final state, and therefore, we can veto events with extra hadrons. The advantage of argon-based detectors over water/ice Cherenkov detectors is their ability to identify hadronic activity to low energies. We explore the details of the DUNE experiment in background discrimination in Sec. 5.4.

• •

Directionality: A feature of the LArTPC technology is its good angular resolution. With an estimated $1^{\circ}$ resolution of low energy electrons, the DUNE experiment will be able to reduce the background for the dSphs searches as the search cone can be as small as the resolution. This resolution has been studied for energies $\mathcal{O}(1~{}\text{GeV})$, and further study from liquid argon experiments should be carried out for a more accurate value for sub-GeV electron energies.

5.3.3 Detector summary

In Table 5.1, we summarize the experiments studied: Super-K and its upgrade Hyper-K for Cherenkov detectors, and two proposed volumes for DUNE as a LArTPC detector. Another detector with a potential of setting some limits on BDM is ICARUS [322, 323] as it ran 5 years deep underground with no cosmic contamination, but we expect Super-K with its present data set to set stronger limits on BDM. As a point of reference, we use the current Super-K exposure of 13.6 years for all experiments in order to estimate limits on BDM.

5.4 Background Modeling

We estimate the number of atmospheric neutrino background events in each experiment in turn.

5.4.1 Cherenkov Detectors

For Super-K and by extension Hyper-K, atmospheric neutrino data are already available, and help estimate the number of neutrino background events expected per year. Since we are not provided the electron spectrum, we use the full data set of events shown in Ref. [253] as the background. We use the fully contained single-ring electron events over the four periods of Super-K, SK-I (1489 days), SK-II (798 days), SK-III (518 days) and SK-IV (1096 days), or for a total of 10.7 years. We estimate the number of background events per year over all energies (provided in two categories sub-GeV and multi-GeV events) to be

[TABLE]

where $V_{\text{exp}}$ is the experimental volume. The number of background events can be scaled up for estimates of Hyper-K. For the BDM search within a cone of angle $\theta_{C}$ around a source, the number of expected background events is then

[TABLE]

which in the case of the GC analysis666Although the optimal value for the opening angle of the search cone depends largely on the DM distribution (J-factor), it also depends on the angular distribution of the scattering process, and has to be optimized separately given a particular scattering. and $\theta_{C}=10^{\circ}$ is

[TABLE]

A proper Super-K analysis can lower these estimates for the background by the use of the full background energy spectrum, and can thus improve the limits on BDM.

5.4.2 LArTPC Detectors

Previous studies have estimated the expected number of fully contained electron events to be 14053 per 350 kton year [65]. Therefore, we take the total number of electron events at DUNE to be 400 events per 10-kton-year. In order to optimize the analysis cuts, we generate a sample of 40,000 simulated atmospheric electron (anti)neutrino scattering events. The reactions inside the pure 40Ar target volume are simulated using the GENIE neutrino Monte Carlo software (v2.10.6) [325]. We model the atmospheric neutrinos with the Bartol atmospheric flux [255]. Since the flux varies slightly with geographic location and with altitude, we use the atmospheric flux available for the nearby MINOS far detector located in the Soudan Mine [326]. We use the neutrino flux that occurs at solar maximum777One expects the most conservative limit to occur at solar minimum. Indeed the flux is higher at solar minimum, but it is dominated by lower energy neutrinos which produce pions. The detection threshold for pions is low enough to improve background rejection at this limit. to provide the most conservative limit. Although charged current processes dominate the background in this energy range, neutral current processes are also simulated.

The dominant primary scattering processes are $\nu_{e}+n\rightarrow p+e^{-}$ and $\overline{\nu}_{e}+p\rightarrow n+e^{+}$. However, due to secondary intranuclear processes the final observable state can, and generally will, include additional hadrons. These are comprised almost entirely of protons, neutrons, pions, and kaons. Table 5.2 summarizes the frequency of different hadrons to be produced in the final state.

Approximately 99.72$\%$ of the simulated interactions contain a free hadron in the final state. This is a useful discriminant as a DM event would not produce a hadron in the final state. So, contingent on detectability, we are able to use these hadrons as a veto on charged current events.

To detect the emitted hadrons, DUNE is able to resolve hadronic activity down to low energy thresholds, provided in Table 5.3 [7]. Neutrons are harder to detect, and to be conservative, we assume that all neutrons escape detection, although future simulations of argon detectors might prove otherwise. Implementing the hadronic veto to the simulated dataset, we find that less than $32\%$ of simulated background processes pass the cut based on hadron tagging alone. We therefore estimate the number of background events over the whole sky to be

[TABLE]

For the searches within $10^{\circ}$ around the GC, the number of background events is

[TABLE]

Using the angular information of the events found by looking up to 10 degrees around the DM sources such as the GC, the background is about 1 event per year.

5.5 Reach at neutrino experiments

We now estimate the experimental sensitivity for BDM searches in the GC, leaving the analysis of dSphs to Sec. 5.6. We compare Cherenkov detectors’ large volume with the LArTPC’s ability to reduce background events through particle identification and explore key experimental features such as low energy thresholds and excellent angular resolution for both technologies.

To measure the sensitivity of an experiment, we define the signal significance as

[TABLE]

where $S$ is the number of signal events, and $B$ the number of background events. In the following, we estimate limits on the region of parameter space defined in Sec. 5.3.3 for a $2\sigma$ significance, using the exposure time shown in Table 5.1.

In Fig. 5.4, we show the $95\%$ limits of Super-K, Hyper-K, and DUNE to the effective cross section $\sigma_{0}$, defined in Eq. (5.17), as a function of $E_{\text{max}}$, defined in Eq. (5.16), in the constant amplitude limit. In the case of light BDM ($2m_{e}E_{B}\gg m_{B}^{2}$), $E_{\text{max}}\approx E_{B}$, while in the case of a heavy BDM ($2m_{e}E_{B}\ll m_{B}^{2}$), $E_{\text{max}}\approx 2m_{e}\gamma_{B}^{2}$. We plot the combination $\sigma_{0}/m_{A}^{2}$ since the number of signal events scales with the number density squared of DM in the case of annihilation.

We also show in Fig. 5.4 as the gray region, the bounds set currently by Super-K without any angular information, having assumed a systematic deviation in the number of events $\delta N_{\text{bkgd}}/N_{\text{bkgd}}=10\%$. This excludes cross sections per mass squared above $\sim 10^{-34}\text{cm}^{2}/\text{GeV}^{2}$. We find that DUNE with 10 kton is almost equally sensitive to BDM signals as Super-K is, for the same exposure, even though DUNE is three times smaller. This is due to its improved background rejection. DUNE can also explore lower electron energies at a comparable angular resolution and therefore lighter BDM.

Although different detector technologies can probe different features, DUNE can test for lighter BDM, while Super-K/Hyper-K can explore lower cross sections due to their large volumes. It is crucial that there is an overlapping region between both experiments; it allows the two experiments to cross-check possible signals and limits, which is especially interesting when comparing different technologies. Detecting a signal in both experiments would be one step towards confirming a DM detection.

5.6 Dwarf Spheroidal Analysis

Dwarf spheroidals are Milky Way satellite galaxies which are dense in DM and low in baryons; they are therefore good candidates for indirect detection searches, with low backgrounds [327, 328]. Although dSphs are less dense in DM than the GC, we can increase the sensitivity to BDM by stacking dSphs. In order to do so, we plot the direction of detected electron events in galactic coordinates, and correlate them with known sources within the experimental angular resolution, such as the dSphs as shown in Fig. 5.5.

5.6.1 J-factor of Dwarf Galaxies

Over the past few years, many dSphs have been found in large surveys [329, 330]. We list in Table 5.4 the locations of the brightest dSphs (in J-factors), the separating distance from the Earth, as well as their found J-factors in decay and annihilation, assuming a NFW profile.

The J-factors listed are integrated over a cone of half angle $0.5^{\circ}$ due to their small extent in the sky. Therefore, in detecting these sources, the search cone (see Fig. 5.3) has to be as small as possible and we therefore choose it to be the experimental angular resolution.

Although individually the J-factors of dSphs are 2 orders of magnitude lower than that of the GC, one can perform a stacked analysis of the dSphs which would effectively sum over the J-factors of all the dSphs considered to set more constraining limits. Such analysis is interesting as it can be a confirmation that a signal is potentially that of DM if it is detected in both the GC and dSphs.

5.6.2 Event Reach

We compute the number of background events as in Sec. 5.4, but here we limit the search angle to the experimental resolution. We find

[TABLE]

where $N_{\text{dSphs}}$ is the number of dSphs considered in the analysis.

Similarly to the GC analysis, we show in Fig. 5.6 the different experimental sensitivities. Although the reach is not as deep as that of the GC analysis, the dSphs analysis would be an excellent confirmation that any potential signal found in the GC is indeed consistent with a DM interpretation. Also, with future surveys, one might be able to push further the dSphs analysis sensitivity by finding more dSphs.

We also point out in this analysis that DUNE with only 10-kton will be able to outperform Super-K due to its excellent background rejection enabled by $1^{\circ}$ angular resolution. One caveat of this analysis is that when reducing the search cone to only 1 degree and 5 degrees for DUNE and Super-K respectively, we are only able to set limits reliably on BDM with a high boost factor $\gamma_{B}$ as the events have to be extremely forward (see App. B.4).

5.7 Conclusions

In this work, we have studied the experimental signatures of a class of DM models called boosted dark matter, in which one component has acquired a large Lorentz boost today and can scatter off electrons in neutrino experiments. Our analysis compared two neutrino technologies: Liquid argon detectors like DUNE and Cherenkov detectors like Super-K and Hyper-K.

We compared the excellent particle identification of LArTPC detectors by simulating neutrino events in argon, with the large volume of Cherenkov experiments to help further reduce the atmospheric neutrino background. Building a search strategy tuned for each experiment extends the physics reach of neutrino detectors from classic DM indirect detection to BDM direct detection, enabled by the ability to tag BDM particle on almost an event-by-event basis, especially in liquid argon experiments.

If the BDM component has a much higher energy than the electron mass, the electron is emitted in the forward direction, and can thus be used to trace back the origin of DM. Such a feature, coupled with a good angular resolution in neutrino experiments can help establish limits on BDM. The angular resolution can also help point back to the origin of DM; constructing a map of the origin of these sources can help correlate signals from neutrino detectors with other experiments, for example gamma rays at Fermi [36].

If a signal is detected, some BDM properties can be extracted. For example, the maximum Lorentz boost for an electron is related to that of the $B$ particle by

[TABLE]

We can therefore extract $E_{e}^{\text{max}}$ from the electron spectrum and obtain the boost factor of $B$. In the case of a monoenergetic signal, where all particles $B$ have energy $E_{B}$, we obtain a single value of $\gamma_{B}$. As we expect low statistics, we can only bound the Lorentz factor from below.

We performed two analyses, one for BDM originating from the GC and one in which we stacked signals from dSphs. We found that DUNE with 10 kton can perform as well as Super-K in the case of the GC analysis, and can outperform it in the dSphs analysis due to its superior angular resolution. In both analyses, we adopted a conservative strategy, in particular by using all atmospheric data across a wide range of energies as background. A dedicated experimental search from the Super-K and DUNE collaborations is able to properly estimate the background and improve the limits on BDM.

The largest constraints affecting the parameter space studied in this work are from our analysis of published Super-K data, where Super-K has not detected any excess of electron events over muon events above statistical fluctuations. Such limits are set without any angular information, and thus can be extended by the Super-K, Hyper-K and DUNE collaborations through a similar analysis to the one described in this work. Other limits, although not discussed above, are model specific and need to be taken into consideration when building a BDM model. These limits include direct detection bounds on any thermal component of a particle interacting with electrons and/or quarks: Direct detection limits on electron scattering are set by the same process that enables the $B$ particle detection at neutrino experiments, but affect the thermal $B$ component instead of the relativistic one [277]. Direct detection limits from proton scattering would affect $B$ particles with masses larger than $\mathcal{O}(1)$ GeV, making the ability of the DUNE experiment to lower the energy detection threshold of utmost importance [275, 339, 85]. Other possible limits include cosmic microwave background (CMB) constraints on the power injected by the thermal $B$ component into SM particles at early redshifts [281]. All these limits need to be studied properly when discussing a particular model of BDM. An example of such study has been implemented in Ref. [5].

DUNE is an excellent detector to cross-check present Cherenkov detectors and extend the reach of neutrino detectors in DM searches. Having multiple technologies for the hunt of DM is key in its eventual detection.

Chapter 6 Conclusions

To summarize, in Chapter 2, I studied the morphology of DM in indirect detection signals using the hydrodynamic simulation Illustris, and found that Galactic signals tend to be symmetric while extragalactic signals are not, due to mergers and DM substructure. In Chapter 3, I showed, using the zoom-in simulation Eris, that metal poor stars from the stellar halo and DM have similar kinematics, and used the velocity dispersions of metal poor halo stars found by SDSS to extract the DM velocity distribution. When this newly found velocity distribution is used, direct detection bounds are found to be almost an order of magnitude weaker at lower DM masses. In Chapter 4, I introduced a new class of DM models called BDM with hybrid direct and indirect detection signals that can be detected in neutrino experiments. Finally, in Chapter 5, I constructed new search strategies for BDM in neutrino experiments adapted to two separate detector technologies, Cherenkov detectors and liquid argon-based detectors.

Throughout this thesis I explored the difference of scales spanned by the problem that is DM: the largest scales (Hubble distance, clusters, and galaxies), and the smallest scales (particle and short-range interaction scales). I used hydrodynamic cosmological simulations to predict properties of DM discernable at direct and indirect detections experiments, and constructed models of DM with interesting experimental signatures. It is crucial to keep bridging the gap between particle physics and astrophysics to solve the problem of DM, and as the resolution as well as our understanding of the baryonic physics of hydrodynamic simulations improves, we will be able to test the cosmological implications of more DM models, and ultimately match observations with simulations. Only then could we discover and understand DM.

Appendix A Spherical Cows

A.1 Analysis of Decaying Dark Matter

In the text, we have performed the quadrant analysis for the case of annihilating DM for an observer located at $R_{\odot}=8.5~{}$ kpc as well as an observer situated well outside the halo. Here we perform a similar analysis for the case of decaying DM.

A.1.1 Galactic Analysis

Similarly to the analysis discussed in Sec. 2.3, we show in Fig. A.1 the distribution of the variables $R_{\text{opp}}$ and $R_{\text{adj}}$ in the case of decaying DM. The differences between the three studied regions (the entire halo, the omitted inner cone of half angle 10∘, and the region up to $30^{\circ}$ from the center of the halo) are less pronounced as the $J$-factor in the case of decay compared to annihilation signals, but the results are consistent with the previous analysis of Sec. 2.3.2.

A.1.2 Extragalactic Analysis

As in Sec. 2.4.2, We show in Fig. A.2 the probability distributions of $R_{\text{opp}}$ and $R_{\text{adj}}$ for the decay signals. The behavior is as expected from previous analyses; asphericity is less pronounced in decay signals compared to annihilation, and the data within a cone of half angle $5^{\circ}$ appears more symmetric. This is due to the off centered subhalos, as discussed in Sec. 2.4.4.

A.2 Mass Correlation

The morphology of halos is highly mass-dependent [97, 100, 130, 340]. We therefore categorize the masses of the halos of the Illustris simulation as follows:

• •

$M_{200}>2\times 10^{12}~{}M_{\odot}$: This subset corresponds to the cluster-sized halos of the simulation.111We have increased this range of masses for larger statistics. This subset is used to compare to the cluster X-ray data.

• •

$10^{10}~{}M_{\odot}<M_{200}<2\times 10^{12}~{}M_{\odot}$: This subset encompasses MW-way like halos as well as slightly less massive halos.

• •

$M_{200}<10^{10}~{}M_{\odot}$: This is the subset for the least massive halos.

In Fig. A.3, we plot the axis ratio for the different mass categories. We find consistent results that the more massive halos are the least expected to be spherical. In Fig. A.4, we plot the angular correlation between the angular momentum vector and the halo’s minor axis, as discussed in Sec. 2.4.3, but now broken by mass category. We find that the most massive halos, which are the least spherical show indeed the most correlation with the baryonic axis.

A.3 Comparison of Asymmetry Parameterization to Axis Ratio

We define the following ratios as they relate to Eqs. 2.8 and 2.9.

[TABLE]

and the related quantities

[TABLE]

In the case that the annihilation/decay signal profiles are modeled as perfect ellipses (as done in e.g. [38, 39]), the axis ratio is simply related to the parameter $r_{\mathrm{opp}}$, in the case where the quadrant boundaries lie at a $45^{\circ}$ degree angle to the major axis.

Consider an arbitrary intensity function of the form $I(\sqrt{(x/a)^{2}+(y/b)^{2}})$ (Note that real annihilation/decay profiles will not in general have this precise form). Let us assume the signal is sufficiently localized that neglecting the curvature of celestial sphere is a reasonable approximation (which will be true if the halo is sufficiently distant, and for a peaked intensity profile is true even for our own Galaxy), so we can define the quadrant boundaries as simply $|x|=|y|$. Then:

[TABLE]

(There are two possible definitions of $r_{\mathrm{opp}}$ in this case, one of which is the reciprocal of the other. This is an arbitrary choice, so we expect the distributions of $r_{\mathrm{opp}}$ and $1/r_{\mathrm{opp}}$ to be identical. In this case, we will arbitrarily choose the quadrants in the numerator to be those lying along the $x$-axis.)

We can restrict ourselves to the region with $x>0,~{}y>0$, and perform the integrals by the substitutions $X=x/a$, $Y=y/b$, followed by $X=R\cos\theta$, $Y=R\sin\theta$. Within this region, this procedure yields:

[TABLE]

Similarly,

[TABLE]

Independent of the boundaries on the integral over $R$ or the details of the function $I(R)$, we thus obtain:

[TABLE]

This result allows us to estimate limits on $r_{\mathrm{opp}}$ or $R_{\mathrm{opp}}$, when provided with limits on the axis ratio for a potential signal modeled as an ellipse, or vice versa.

In the limit where the ratio of major to minor axes is large (either $a/b\ll 1$ or $a/b\gg 1$), $r_{\mathrm{opp}}$ approaches $2(a/b)/\pi$ if $a/b\ll 1$, and $\pi(a/b)/2-1$ if $a/b\gg 1$. In the limit where $a/b\approx 1$, $r_{\mathrm{opp}}\approx 1+\frac{4}{\pi}(a/b-1)$, and $R_{\mathrm{opp}}\approx(2/\pi)|a/b-1|$. Thus we may consider $r_{\mathrm{opp}}$ with a specific choice of quadrants as a rough proxy for axis ratio (while being more general, and well-defined for cases where the signal is not actually elliptical), with the approximation being most accurate for near-spherical intensity profiles. (Note that there will be corrections to Eq. A.6 associated with the spherical coordinate system of the sky, and with any boundaries on the region of interest that are not only functions of $R$; if greater accuracy is desired in this conversion, $r_{\mathrm{opp}}$ should be computed numerically for the intensity profile and region of interest under study.)

Note that the parameter $r_{\mathrm{adj}}$ is identically 1 in the context of perfectly elliptical signal models, as the signal profile will be evenly bisected by any axis passing through its center.

Appendix B Boosted Dark Matter

B.1 Analytic Approximations to Relic Abundances

The coupled Boltzmann equations for the evolution of the $\psi_{A}$/$\psi_{B}$ abundances are

[TABLE]

where the factor of $\frac{1}{2}$ arises because $\psi_{A}$ and $\psi_{B}$ are Dirac fermions, and $n_{A}$ refers to the sum of the abundances for $\psi_{A}$ and $\overline{\psi}_{A}$ (and similarly for $n_{B}$). In terms of the comoving abundance $Y_{i}=n_{i}/s$, where $s$ is the entropy of the universe, and $x\equiv m_{B}/T$, we can rewrite the Boltzmann equations as

[TABLE]

where we have introduced the shorthand notations:

[TABLE]

and used the fact the total DM number is not changed by the $\psi_{A}\overline{\psi}_{A}\rightarrow\psi_{B}\overline{\psi}_{B}$ reaction, i.e.

[TABLE]

Obtaining accurate solutions requires solving the above coupled equations numerically. In much of the parameter space of interest, however, it is possible to obtain good analytic approximations based on two effectively decoupled equations. When $m_{B}<m_{A}$ and $\lambda_{B}\gg\lambda_{A}$, $\psi_{B}$ typically freezes out of equilibrium well after $\psi_{A}$ does. Therefore, the evolution of $Y_{A}$ in Eq. (B.2) becomes the conventional Boltzmann equation for one species of DM by taking $Y_{B}\approx Y^{\rm eq}_{B}$ at least up until the $\psi_{A}$ freeze-out time.111After $\psi_{B}$ freezes out, $Y_{B}\approx Y^{\rm eq}_{B}$ is invalid, so the two equations formally “re-couple”. Since $Y_{A}$ has approached its asymptotic value by then, though, it is insensitive to late-time details. In the case of $s$-wave annihilation of our interest, the relic abundance of $\psi_{A}$ can be well approximated by the familiar result [341] (with an extra factor of 2 to account for both $\psi_{A}$ and $\overline{\psi}_{A}$)

[TABLE]

where $T_{f,A}=m_{B}/x_{f,A}$ is the freeze-out temperature for $\psi_{A}$, and in the last step we used $sx^{3}/2H(m_{B})=0.132(g_{*s}/g_{*}^{1/2})M_{pl}m_{B}$.

The solution for $Y_{B}$ is more subtle, but can also be greatly simplified when the freeze-out times of $\psi_{A}$ and $\psi_{B}$ are well separated. If $x_{f,B}\gg x_{f,A}$, then we can drop terms suppressed by $(Y_{A}^{\rm eq}/Y_{B}^{\rm eq})^{2}$ in Eq. (B.3), and we can treat the effect of $\psi_{A}$ on $\psi_{B}$ freeze-out by taking $Y_{A}(x_{f,B})\simeq Y_{A}(x_{f,A})\simeq Y_{A}(\infty)$. Defining $\Delta\equiv Y_{B}-Y_{B}^{\rm eq}$, we rewrite Eq. (B.3) as:

[TABLE]

Focussing on the epoch when $\psi_{B}$ starts to deviate from equilibrium, we can apply the ansatz $\Delta=c\,Y_{B}^{\rm eq}$, where $c$ is $\mathcal{O}(1)$. The equilibrium distribution for $x\gg 1$ is

[TABLE]

where we only keep the leading power term in $x$ in the second line. Combining all these, we can rewrite Eq. (B.7) as a quadratic equation for $Y_{B}^{\rm eq}$,

[TABLE]

whose real positive solution is

[TABLE]

We can then equate this equation with $Y_{B}^{\rm eq}(x)\simeq x^{3/2}e^{-x}$ to solve numerically for $x_{f,B}$.

We can see that by removing the contribution from $\psi_{A}$ (i.e. the term $\propto\lambda_{A}Y_{A}^{2}(\infty)$) in Eq. (B.11), $\psi_{B}$ freezes out in the standard way. In particular, we have the approximate relation $x_{f,B}\simeq\log\lambda_{B}-\frac{1}{2}\log x_{f,B}$ which yields

[TABLE]

in analogy with Eq. (B.6). We also see that Eq. (B.11) approaches the standard freeze-out solution when $\lambda_{B}$ decreases and approaches $\lambda_{A}$, such that $\psi_{B}\overline{\psi}_{B}\rightarrow\gamma^{\prime}\gamma^{\prime}$ freezes out at temperatures comparable to $\psi_{A}\overline{\psi}_{A}\rightarrow\psi_{B}\overline{\psi}_{B}$; in that regime, the effect of $\psi_{A}$ on the $\psi_{B}$ evolution is subdominant since $Y_{A}^{\rm eq}<Y_{B}^{\rm eq}$ for $m_{A}>m_{B}$. Standard freeze-out of $\psi_{B}$ continues to hold when $\lambda_{B}\ll\lambda_{A}$, though the approximate solution Eq. (B.11) would not be valid in that regime, since $\Omega_{B}>\Omega_{A}$, in contradiction to our ansatz that $\psi_{A}$ constitutes the major DM component.

More surprising is the case of large $\lambda_{B}$. The $Y_{A}^{2}(\infty)$ term in Eq. (B.11) dominates when

[TABLE]

where we have estimated $x_{f,A}/x_{f,B}\simeq m_{A}/m_{B}$. Taking $Y_{B}^{\rm eq}(x_{f,B})\simeq Y_{B}(\infty)$, Eq. (B.11) reduces to

[TABLE]

This behavior is very strange from the point of view of standard freeze-out, since the abundance of $\psi_{B}$ scales like $1/\sqrt{\sigma_{B}}$ (instead of like the expected $1/\sigma_{B}$). A naive quick way of understanding this behavior is by setting $dY_{B}/dx\approx 0$ in Eq. (B.3) and dropping all $Y_{i}^{\rm eq}$ terms at late times, which immediately leads to Eq. (B.14). We call this “balanced freeze-out”, since the abundance of $\psi_{B}$ is set by the balance between a depleting term ($\propto\lambda_{B}Y_{B}^{2}$) and a replenishing term ($\propto\lambda_{A}Y_{A}^{2}$). Unlike in ordinary freeze-out where the expansion of the universe plays a key role in setting the abundance, in balanced freeze-out the main effect of the Hubble expansion is simply to drive $Y_{i}^{\rm eq}$ to zero at late times.

B.2 Direct Detection of Non-Boosted DM

In this paper, we have largely assumed that $\psi_{A}$ has no couplings to the SM. Given the contact interaction in Eq. (4.3), though, $\psi_{A}$ can interact with the dark photon via $\psi_{B}$ loops. In this appendix, we consider the direct detection bounds on $\psi_{A}$ from these loop processes. Of course, as with $\psi_{B}$, one can relax direct detection limits by giving $\psi_{A}$ an inelastic mass splitting.

The loop-induced couplings of $\psi_{A}$ to the dark photon depend on the UV completion of Eq. (4.3), and we consider exchange of a complex scalar $\phi$ with $U(1)^{\prime}$ charge as a concrete example. The Lagrangian for this system is

[TABLE]

where $D_{\mu}=\partial_{\mu}-ig^{\prime}A^{\prime}_{\mu}$. Integrating out $\phi$ yields the contact interaction in Eq. (4.3) with

[TABLE]

Through $\psi_{B}$-$\phi$ loops, $\psi_{A}$ acquires a coupling to the dark photon. In the limit $m_{\phi}\gg m_{A}\gg m_{B}$, the $\psi_{B}$-$\phi$ loop generates the effective dimension six operator

[TABLE]

which can lead to $\psi_{A}$-nucleon scattering as in Fig. B.1a. As discussed in the appendix of Ref. [342], the standard dimension five dipole operator $\overline{\psi}_{A}\Sigma^{\mu\nu}\psi_{A}F^{\prime}_{\mu\nu}$ does not appear after integrating out $\psi_{B}$ and $\phi$, because the interactions in Eq. (B.15) respect a chiral symmetry acting on $\psi_{A}$.

Similar to Ref. [342] (but replacing the photon with a dark photon), the dominant effect of Eq. (B.17) is to give rise to a charge-charge interaction between DM and a nucleus $N$. The spin-independent $\psi_{A}N\rightarrow\psi_{A}N$ cross section is

[TABLE]

where $m_{N}$ is the nucleus mass, $E_{R}$ is the nucleus recoil energy, $t=-2m_{N}E_{R}$ is the momentum-transfer-squared, $v$ is the DM velocity, $Z$ is the nucleus charge number, and $F^{2}(E_{R})$ is the nucleus charge form factor. The numerator in the expression above corresponds just to the lowest term in an expansion in $t$ (i.e. small momentum transfer). Spin-independent bounds on DM typically assume equal couplings to neutrons and protons, and can be expressed in terms of an effective nucleon cross section $\sigma_{n}$, with

[TABLE]

where $\mu$ is the DM-nucleon reduced mass, and $A$ is the nucleus mass number. Thus, we have

[TABLE]

Note that this cross section is momentum dependent, but for simplicity, we will take $E_{R}\simeq 10~{}\text{keV}$ to determine the typical value of $t$.

Near the benchmark in Eq. (4.13), $m_{A}$ is heavier than the proton so $\mu\simeq m_{p}$. The DM-nucleon cross section scales roughly as

[TABLE]

where we have set $\lambda=1$, ignored the logarithmic dependence on $m_{\phi}$ and $m_{B}$, and ignored the $m_{\gamma^{\prime}}$ dependence since $m_{\gamma^{\prime}}$ is comparable to the typical values of $t$. Since we adjust $\Lambda$ (equivalently $m_{\phi}/\lambda$) to get the right abundance of DM, and since $\Lambda^{4}\approx m_{A}^{2}/\langle\sigma_{A\bar{A}\rightarrow B\bar{B}}v\rangle$ from Eq. (4.8), we can rewrite this dependence as:

[TABLE]

In Fig. B.1b, we show the limits of the LUX experiment [272] on the direct detection of $\psi_{A}$ for different values of $(m_{A},m_{B},m_{\gamma^{\prime}})$, and see that these constraints are easily satisfied, though future direct detection experiments would have sensitivity.

B.3 Boosted DM Scattering Off Hadrons

In Sec. 4.4.2, we focused on the $\psi_{B}e^{-}\rightarrow\psi_{B}e^{-}$ mode for boosted DM detection. Here, we summarize the signal event rate for boosted DM scattering off protons or nuclei. Since the number of signal events is rather small, we have not pursued a background study, though we remark that the angular pointing for hadronic scattering is rather poor at the low energies we are considering. As discussed in the main text, both event rate and angular resolution for scattering off proton can be improved at liquid Argon detectors.

B.3.1 Elastic Scattering Off Hadrons

The elastic scattering $\psi_{B}N\rightarrow\psi_{B}N$ has similar kinematics to electron scattering (with the obvious replacement $m_{e}\to m_{N}$), except we have to include the electromagnetic form factor. We will express the cross section as a function of the electric and magnetic Sachs form factors $G_{E}$ and $G_{M}$. For protons, we can use the dipole approximation

[TABLE]

To compute the cross section, we use the Rosenbluth formula in the lab frame as cited in [343]

[TABLE]

The energies and momenta are defined the same way as Eq. (4.23), $M$ is the mass of the proton, and $\theta$ is the scattering angle of $\psi_{B}$.

The lowest momentum for a proton to Cherenkov radiate is $1.2~{}\text{GeV}$, and for our benchmark in Eq. (4.13), the proton cross section above this threshold is

[TABLE]

yielding an all-sky event rate of

[TABLE]

Due to the presence of the Cherenkov cutoff and the proton form factor, the elastic scattering rate in Eq. (B.25) varies little within the mass range of interest, as given in Eq. (4.7). When the transferred energy is above $2.5~{}\text{GeV}$, the elastic scattering cross section is rather small, and protons instead typically produce secondary hadronic showers [248]. In that case, one should transition to the DIS calculation below.

As mentioned in footnote 9, large volume liquid Argon detectors are able to detect scattered protons with energies much below the Cherenkov threshold using ionization signals, where the total elastic scattering cross section off protons would be relevant. We find the total elastic $\psi_{B}\,p\to\psi_{B}\,p$ scattering cross section for boosted $\psi_{B}$ with $m_{A}\gtrsim 1$ GeV to be

[TABLE]

which is insensitive to $m_{A}$, $m_{B}$ in the boosted $\psi_{B}$ regime due to the proton form factor. We see that the total elastic scattering rate off proton is much larger than the one with a Cherenkov cutoff, so $\psi_{B}\,p\to\psi_{B}\,p$ could potentially be the leading signal detectable at a liquid Argon detector.

Generalizing the previous calculation to a coherent nucleus of charge number Z:

[TABLE]

This same $\psi_{B}\,p\to\psi_{B}\,p$ calculation is relevant for direct detection of non-relativistic relic $\psi_{B}$. Taking the $q^{2}\to 0$ limit and integrating over all angles, we have the cross section

[TABLE]

where $\mu_{p}=m_{p}m_{B}/(m_{p}+m_{B})$ is the reduced mass of the dark matter and the proton. This is the basis for Eq. (4.53) shown earlier.

B.3.2 Deep Inelastic Scattering Off Hadrons

At sufficiently high energies, $\psi_{B}$ scattering off hadrons will behave more like deep inelastic scattering (DIS), where the final state is a hadronic shower. The DIS cross section is a convolution of the parton-level cross section with parton distribution functions (PDFs). The parton-level cross section $\hat{\sigma}$ is given by

[TABLE]

For the $\psi_{B}$-parton system: $\hat{s}+\hat{u}+\hat{t}=2m_{B}^{2}$, $\hat{t}=-Q^{2}$, and $\hat{s}=(1-x)m_{B}^{2}+xs$. We define $x$ by $p\equiv xP$ where $P$ is the 4-momentum of the initial proton at rest. We define $y\equiv\frac{2P\cdot q}{2P\cdot k}=\frac{-\hat{t}}{\hat{s}-m_{B}^{2}}$, which characterizes the fraction of the energy transferred from $\psi_{B}$ to the parton, since $y=\frac{q^{0}}{k^{0}}=1-\frac{E^{\prime}}{E}$ in the rest frame of the initial proton.

From these relations, we get the transferred momentum $Q^{2}=xy(s-m_{B}^{2})$, and $dx\,dQ^{2}=\frac{dQ^{2}}{dy}dx\,dy=x(s-m_{B}^{2})\,dx\,dy$. Including parton distribution functions, and using $x$/$y$ as variables, we obtain the resulting DIS cross section:

[TABLE]

where $f_{f}(x,Q)$ are PDFs with f indicating different flavor of fermion. For numerical evaluation, we use the MSTW2008 NNLO PDFs from Ref. [344]. The integration limits of Eq. (B.31) are $0\leq x\leq 1$ and $0\leq y\leq y_{\rm max}$, where applying the condition $\cos\theta\leq 1$ we obtain

[TABLE]

with $s=m_{B}^{2}+M_{p}^{2}+2M_{p}E_{B}$ and $E_{B}=m_{A}$. Unlike the familiar case of DIS initiated by nearly massless incoming particles, for the massive $\psi_{B}$ we consider here, $y_{\rm max}$ is not trivially $1$.

Since the PDFs are only reliable for transferred energies over $\sim 1$ GeV, we impose $Q^{2}\geq(1~{}\rm GeV)^{2}$ as a default cut for numerical integration. Analogous to the discussion for elastic scattering signals, for a particular experiment, there may be harder cut on phase space due to detector threshold energy. For our benchmark in Eq. (4.13), the DIS cross section above the $1~{}\text{GeV}$ threshold is

[TABLE]

yielding

[TABLE]

for the all-sky event rate at Super-K.

B.4 Understanding Forward Scattering

In Sec. 5.3, we assumed that when the energy $E_{B}$ of the boosted particle is greater than the electron mass

[TABLE]

the final state electron of the elastic scattering $Be^{-}\rightarrow Be^{-}$ is emitted in the forward direction. This is crucial as the observed electron can then point back to the origin of the $B$ particle. From kinematics, the scattering angle of the emitted electron relative to the incoming $B$, labeled $\theta_{e}^{\prime}$ as shown in Fig. 5.3, is

[TABLE]

where the energy of the emitted electron is $E_{e}$. Applying the assumption of Eq. (B.35), Eq. (B.36) becomes

[TABLE]

where

[TABLE]

with $i\in\{B,e\}$ being the $B$ and electron boost factors. We have expanded in large $\gamma_{e}$ and $\gamma_{B}$ in Eq. (B.4).

In the cases where $\gamma_{B},\gamma_{e}\gg 1$, we find to a good approximation that $\cos\theta^{\prime}_{e}\approx 1$ and $\sin\theta^{\prime}_{e}\approx 0$. The angle of the recoiled electron relative to the DM source $\theta_{e}$ is related to $\theta^{\prime}_{e}$ by

[TABLE]

where $\phi^{\prime}_{e}$ is the azimuthal angle of the recoiled electron with respect to the incoming $B$ as shown in Fig. 5.3 and is uniformly distributed between [math] and $2\pi$.

In order to estimate the error on the measured angle $\theta_{e}$ compared to the incoming $B$ angle $\theta_{B}$, we study the deviations in Eq. (B.39) from $\cos\theta_{e}=\cos\theta_{B}$. Taylor expanding around $\theta^{\prime}_{e}=0$, we find

[TABLE]

From Eq. (B.4), and in terms of the boost factors $\gamma_{e}$ and $\gamma_{B}$,

[TABLE]

The last approximation is found from the kinematics relation $\gamma_{e}^{\text{max}}=2\gamma_{B}^{2}-1$ and therefore $\gamma_{e}<2\gamma_{B}^{2}$. Taking $\sin\phi^{\prime}_{e}=1$ as its maximum value, we find that the deviation from the forward approximation is

[TABLE]

We show the results of the ratio of the observed electron angle by the incoming $B$ angle $\theta_{e}/\theta_{B}$, as a function of the $B$ angle $\theta_{B}$ in Fig. B.2. For every value of $\gamma_{e}$ found, there exists a minimal gamma factor of the original particle $B$ such that $\gamma_{e}=2(\gamma_{B}^{\text{min}})^{2}-1$. The solid curves in Fig. B.2 correspond to the ratio $\theta_{e}/\theta_{B}$ with

[TABLE]

for different values of $\gamma_{B}^{\text{min}}$. We find that values of $\gamma_{B}>20$ are suitable within the forward scattering approximation, with errors less than $20\%$. We also study the largest value of $1/\gamma_{e}$, which occurs at the experimental threshold $E_{\text{thresh}}$

[TABLE]

As discussed in Sec. 5.3, the experiment thresholds considered are $E_{\text{thresh}}=30$ MeV and $E_{\text{thresh}}=100$ MeV, which lead to a gamma factor of $\gamma_{e}^{\text{min}}=60-200$. We show the measured angle of the electron off the source as a function of the initial BDM angle $\theta_{B}$ for the events right at the energy threshold in dashed lines in Fig. B.2. This study can be properly incorporated within the experimental framework to estimate the systematics as a function of the emitted electron’s energy.

B.5 Comparing the Full Analysis with a Concrete Model

In this section, we summarize the model explored in Ref. [5], based on Ref. [188], and show the reach of the DUNE experiments in the appropriate parameter space. We start with a multicomponent DM model with two particle species $A$ and $B$, such that $A$ is the dominant DM component that interacts solely with $B$, and $B$ is the subdominant component that couples to the standard model. If a mass hierarchy exists such that $m_{A}\gg m_{B}$, the annihilation process $A\overline{A}\rightarrow B\overline{B}$ leads to particles $B$’s with energies $E_{B}=m_{A}$ and thus a high boost factor $\gamma_{B}=m_{A}/m_{B}$.

We further take the $B$-SM couplings to be through the kinetic mixing of a dark photon $\gamma^{\prime}$ with the photon. The mixing term is

[TABLE]

where $F^{\prime\mu\nu}$ is the dark photon field, $F^{\mu\nu}$ is the photon field, and $\epsilon$ is the coupling of the interaction. We take the coupling of $B$ to the dark photon to be $g^{\prime}$, which is large but perturbative. The model parameters are therefore:

[TABLE]

The cross section of the $A-\overline{A}$ annihilation (see the left diagram of Fig. B.3) is set such that we obtain the right abundance of $A$’s today, which brings the value of the cross section close to the thermal cross section. The abundance of $B$ particles is controlled by both the annihilation of the $A$ diagram as well as the annihilation of the $B$ diagram (middle diagram of Fig. B.3).

Finally, the scattering of $B$ particles off electrons is set by the right diagram of Fig. B.3. The same diagram with a nucleon instead of an electron is the one that sets direct detection bounds on the thermal component of $B$. This study focuses however on $B$ particles with masses below the ones studied so far in direct detection experiments. Of course higher $B$ masses can be evaded by the introduction of inelastic scattering [232, 233].

For Fig. B.4, we use the following benchmark (while varying $m_{A}$ and $m_{B}$), where the limits on the dark photon are consistent with those in Ref. [345].

[TABLE]

In Fig. B.4, we show the estimated limits of DUNE as well as Super-K and Hyper-K in the $m_{A}-m_{B}$ space, first presented in Ref. [5]. We find consistent results with Fig. 5.4, as Super-K and DUNE with 10 kton have similar sensitivity, with DUNE able to probe lower electron recoils. This is shown by the diagonal line in the triangular range of Fig. B.4 which can be thought of as the difference between $m_{A}$ and $m_{B}$, a quantity that is related to the energy of the emitted electron.