#  Amplified J-factors in the Galactic Center for velocity-dependent dark
matter annihilation in FIRE simulations

Daniel McKeown1, James S. Bullock1, Francisco J. Mercado 1, Zachary Hafen1,
Michael Boylan-Kolchin2, Andrew Wetzel3, Lina Necib4, Philip F. Hopkins4,
Sijie Yu1
1Center for Cosmology, Department of Physics and Astronomy,4129 Reines Hall,
University of California Irvine, CA 92697, USA
2Department of Astronomy, The University of Texas at Austin, 2515 Speedway
Stop C1400, Austin, TX 78712, USA
3Department of Physics, University of California, Davis, CA 95616, USA
4TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA
91125, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ; Accepted XXX. Received YYY;
in original form ZZZ)

###### Abstract

We use FIRE-2 zoom cosmological simulations of Milky Way size galaxy halos to
calculate astrophysical J-factors for dark matter annihilation and indirect
detection studies. In addition to velocity-independent (s-wave) annihilation
cross sections $\langle\sigma v\rangle$, we also calculate effective J-factors
for velocity-dependent models, where the annihilation cross section is either
either p-wave ($\propto v^{2}/c^{2}$) or d-wave ($\propto v^{4}/c^{4}$). We
use 12 pairs of simulations, each run with dark-matter-only (DMO) physics and
FIRE-2 physics. We observe FIRE runs produce central dark matter velocity
dispersions that are systematically larger than in DMO runs by factors of
$\sim 2.5-4$. They also have a larger range of central ($\sim 400$ pc) dark
matter densities than the DMO runs ($\rho_{\rm FIRE}/\rho_{\rm DMO}\simeq
0.5-3$) owing to the competing effects of baryonic contraction and feedback.
At 3 degrees from the Galactic Center, FIRE J-factors are $5-50$ (p-wave) and
$15-500$ (d-wave) times higher than in the DMO runs. The change in s-wave
signal at 3 degrees is more modest and can be higher or lower ($\sim 0.3-6$),
though the shape of the emission profile is flatter (less peaked towards the
Galactic Center) and more circular on the sky in FIRE runs. Our results for
s-wave are broadly consistent with the range of assumptions in most indirect
detection studies. We observe p-wave J-factors that are significantly enhanced
compared to most past estimates. We find that thermal models with p-wave
annihilation may be within range of detection in the near future.

###### keywords:

galaxies: disc – glaxies: formation – cosmology: dark natter – cosmology:
theory

††pubyear: 2020††pubyear: 2021††pagerange:  Amplified J-factors in the
Galactic Center for velocity-dependent dark matter annihilation in FIRE
simulations –B

##  1 Introduction

There is significant astrophysical and cosmological evidence showing that non-
baryonic dark matter dominates the mass in the Universe (e.g. Einasto et al.,
1974; Rubin et al., 1978; Trimble, 1987; Wittman et al., 2000; Mandelbaum et
al., 2013; Kwan et al., 2017). Cosmological constraints have pinpointed that
the dark matter mass density relative to critical is $\Omega_{\rm dm}\simeq
0.27$ today (see Planck Collaboration et al., 2020, and references therein).
For a comprehensive historical perspective on the observational and
theoretical motivations for dark matter see Bertone & Hooper (2018).

One popular theory for dark matter suggests that it is made up of Weakly
Interacting Massive Particles (WIMPs) (WIMPs, see Jungman et al., 1996;
Dodelson, 2003; Feng, 2004; Bertone et al., 2005; Mo et al., 2010, for
discussions). In the standard WIMP scenario, where dark matter particles are
their own antiparticles, WIMPs self-annihilate and recombine in equilibrium
when the Universe is young, hot, and dense. As the Universe cools and expands,
annihilation rates become too low to maintain equilibrium, and the co-moving
particle abundance “freezes out.” The resultant abundance is set directly by
the interaction cross section during freeze-out, and this gives us a way to
relate a macroscopic observable ($\Omega_{\rm dm}\simeq 0.27$) to microscopic
properties of the particle. Specifically, if the thermally-averaged cross-
section during freeze out is
$\langle\sigma_{A}v\rangle\equiv\langle\sigma_{A}v\rangle_{T}\simeq 2.3\times
10^{-26}$ cm3 s-1, then the thermal relic density is naturally of the right
order of magnitude to match the observed abundance (see Steigman et al., 2012,
for a precise treatment).

The same annihilations that set the thermal abundance of WIMPs in the early
Universe should be occurring again today in regions where the dark matter has
become dense in dark matter halos. If those annihilations produce Standard
Model particles, this provides a means for indirect dark matter detection.
Specifically, an observed flux of Standard Model particles from such a
location could provide evidence for dark matter. One region of particular
interest is the Galactic Center (Bergström et al., 1998). Not only is the
Galactic Center expected to be dense in dark matter but its relative proximity
to Earth has made it a subject of significant study for indirect messengers of
annihilation, including cosmic rays and neutrinos (see Gaskins, 2016, for a
review). If annihilation to quarks and charged lepton states happens for dark
matter particles, this will ultimately produce photons with energies of order
$\sim 10\%$ of the dark matter particle mass, making gamma-ray observations of
particular interest for indirect searches for WIMP dark matter with
$m_{\chi}\sim 100$ GeV (Daylan et al., 2016).

An observed excess in gamma-ray emission from the Galactic Center based on
Fermi Large Area Telescope observations has sparked considerable interest as a
potential indirect detection signal (Hooper & Goodenough, 2011; Abazajian &
Kaplinghat, 2012; Martin et al., 2014; Ajello et al., 2016; Ackermann et al.,
2017). The basic excess has been confirmed by multiple groups (see Murgia,
2020, for a review) and is consistent with expectations for a dark matter
particle with mass $m_{\chi}\sim 10-100$ GeV annihilating with a velocity-
averaged cross section that matches the thermal WIMP expectation. A different
signal from the Andromeda galaxy halo is potentially consistent with this
interpretation (Karwin et al., 2020). Although a dark matter origin of the
Galactic Center Excess (GCE) is the most intriguing possibility, astrophysical
sources, including gamma-ray emitting pulsars (e.g. Abazajian, 2011; Bartels
et al., 2016) and supernovae remnants (e.g. Carlson & Profumo, 2014) are
plausible alternatives.

The case for an astrophysical interpretation has been strengthened over the
last several years, with analyses showing that the morphology of the excess
traces the flattened "boxy" stellar over-density of the Galactic bulge, rather
than the more spherical distribution one would expect for a dark matter signal
(Macias et al., 2018; Bartels et al., 2018). Based on this realization,
Abazajian et al. (2020) used templates for the galactic and nuclear stellar
bulges to show that the GC shows no significant evidence for DM annihilation
and used this to place strong constraints on the s-wave cross section. In
particular, in the case of a pure b-quark annihilation channel, assuming a
range of DM profiles consistent with numerical simulations, Abazajian et al.
(2020) ruled out s-wave cross sections $\gtrsim
0.015\langle\sigma_{A}v\rangle_{T}$ for dark matter masses $m_{\chi}\sim 10$
GeV.

One way that this thermal limit could be avoided is if dark matter
annihilation is velocity-dependent (Robertson & Zentner, 2009; Giacchino et
al., 2013; Choquette et al., 2016; Petač et al., 2018; Boddy et al., 2018;
Johnson et al., 2019; Argüelles et al., 2019; Board et al., 2021).
Specifically, in some models, symmetries forbid the s-wave contribution to the
annihilation cross-section, and the leading contribution to DM annihilations
could be p-wave $\sigma v\propto(v/c)^{2}$ or even d-wave $\sigma
v\propto(v/c)^{4}$. In the Milky Way, typical DM velocities are usually
thought to be $\sim 3\times 10^{-4}c$ near the Galactic Center, while at
thermal freeze out $v\sim 0.1c$. This means that for p-wave, the cross section
is expected to be suppressed by a factor of $\sim 10^{-5}$ compared to the
value during freeze-out.

For a fixed particle physics model with s-wave annihilation, the expected
annihilation signal depends on the square of the dark matter density along the
line-of-sight from the observer. This "astrophysical J-factor" is therefore
critical to the interpretation of any indirect dark matter search (Stoehr et
al., 2003; Diemand et al., 2007; Springel et al., 2008; Kamionkowski et al.,
2010; Grand & White, 2021). It is common for interpretive analyses to adopt
analytic priors for the dark matter profile shape inferred from cosmological
simulations and to normalize the profiles so that the local dark matter
density near the Sun matches observationally-inferred values (Necib et al.,
2019). For velocity-dependent models, the J-factor is generalized to include
the local velocity distribution (Boddy et al., 2018). Recently, Board et al.
(2021) used several cosmological zoom hydrodynamical simulations to
investigate the generalized J-factor for velocity dependent models for Milky-
Way size galaxies. They found that J-factors were enhanced for hydrodynamic
runs in p and d-wave cases and that the relative velocity distribution of dark
matter particles was well-described by a Maxwell-Boltzmann form. They also
concluded that the J-factor in all models was strongly correlated with the
local Solar dark matter density.

In this paper, we perform a similar analysis to that in Board et al. (2021)
utilizing 12 Milky Way mass zoom simulations with $\sim 10$ times better mass
resolution done as part of the FIRE-2 collaboration (Wetzel et al., 2016;
Garrison-Kimmel et al., 2017a; Hopkins et al., 2018a; Garrison-Kimmel et al.,
2019a; Lazar et al., 2020). These simulations represent some of the highest
resolution simulations of Milky-Way mass halos with full galaxy formation
physics yet run, and allow us to resolve the J-factor within $\sim 3$ degrees
of the Galactic Center. For each halo, we have a dark-matter-only version, and
this allows us to explore the differential effect of galaxy formation physics
on J-factor predictions. One major difference from Board et al. (2021) is that
we re-normalize every halo so that the local dark matter density at mock solar
locations are identical. This allows us to mimic what is done in indirect
detection analyses and to explore how differences in the shape of the dark
matter density and velocity profile will affect J-factor predictions from
simulation to simulation. The outline of this paper is as follows. In section
2 we describe our simulations. In section 3 we provide our nomenclature for
the effective J-factor and describe our analysis. In section 4 present our
results and conclude in section 5.

Table 1: (1) Simulation name. The suffix “DMO" stands for “Dark Matter Only" and refers to the same simulation run with no hydrodynamics or galaxy formation physics. (2) Factor $f$ by which dark matter particle masses have been multiplied ($m_{\rm dm}\rightarrow fm_{\rm dm}$) in order to normalize the dark matter density at $\rho(r=R_{\odot})=10^{7}$ $M_{\odot}$${\rm kpc}^{-3}=0.38$GeV cm-3 for a mock solar location $R_{\odot}=8.3$ kpc. (3) Stellar mass $M_{\star}$ of the central galaxy. (4) Virial mass (of raw simulation, not including the $f$ factor) defined by Bryan & Norman (1998). The following quantities are derived after normalizing (by $f$) to the local dark matter density: (5) Cumulative s-wave J-factor within 3 degrees of the Galactic Center. (6) Cumulative s-wave J-factor within 10 degrees of the Galactic Center. (7) Cumulative s-wave J-factor integrated over the sky. (8) Cumulative p-wave J-factor within 3 degrees of the Galactic Center. (9) Cumulative p-wave J-factor within 10 degrees of the Galactic Center. (10) Cumulative p-wave J-factor integrated over the sky. (11) Cumulative d-wave J-factor within 3 degrees of the Galactic Center. (12) Cumulative d-wave J-factor within 10 degrees of the Galactic Center. (13) Cumulative d-wave J-factor integrated over the sky.  Simulation | f | M⋆ | Mvir | Js(< $3^{\circ}$) | Js(< $10^{\circ}$) | J${}_{s}^{\rm tot}$ | Jp(< $3^{\circ}$) | Jp(< $10^{\circ}$) | J${}_{p}^{\rm tot}$ | Jd(< $3^{\circ}$) | Jd(< $10^{\circ}$) | J${}_{d}^{\rm tot}$
---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  $10^{10}$ M⊙ |  $10^{12}$ M⊙ |  | ($10^{22}$ | GeV2 cm-3) |  | ($10^{16}$ | GeV2 cm-3) |  | ($10^{10}$ | GeV2 cm-3)
M12i | 1.28 | 6.4 | 0.90 | 0.964 | 5.08 | 12.5 | 4.62 | 21.1 | 42.2 | 121 | 477 | 806
M12iDMO | 1.59 | - | 1.3 | 0.295 | 1.32 | 5.21 | 0.133 | 0.881 | 5.34 | 0.390 | 3.52 | 31.2
M12c | 1.26 | 6.0 | 1.1 | 1.03 | 5.63 | 15.5 | 4.32 | 20.9 | 46.6 | 98.5 | 423 | 790
M12cDMO | 1.83 | - | 1.3 | 1.26 | 4.62 | 12.3 | 0.545 | 2.93 | 11.3 | 1.42 | 10.9 | 59.6
M12m | 0.885 | 11 | 1.2 | 0.603 | 3.94 | 13.6 | 2.25 | 13.3 | 38.2 | 45.5 | 243.0 | 594
M12mDMO | 1.42 | - | 1.4 | 1.05 | 3.90 | 9.41 | 0.490 | 2.42 | 8.10 | 1.31 | 8.60 | 40.3
M12f | 1.01 | 8.6 | 1.3 | 0.844 | 5.31 | 14.5 | 3.65 | 20.0 | 44.4 | 86.5 | 411 | 765
M12fDMO | 1.82 | - | 1.6 | 0.620 | 2.38 | 6.67 | 0.303 | 1.61 | 7.02 | 0.890 | 6.46 | 44.2
M12w | 1.28 | 5.8 | 0.83 | 0.921 | 3.99 | 11.2 | 3.58 | 14.5 | 33.1 | 74.8 | 285 | 550
M12wDMO | 1.68 | - | 1.1 | 0.436 | 1.72 | 5.66 | 0.207 | 1.18 | 5.68 | 0.610 | 4.78 | 32.5
M12b | 0.990 | 8.1 | 1.1 | 1.04 | 6.49 | 16.6 | 5.82 | 29.7 | 58.0 | 182 | 759 | 1201
M12bDMO | 1.25 | - | 1.4 | 1.25 | 4.37 | 11.2 | 0.655 | 3.10 | 11.0 | 2.08 | 12.9 | 63.1
Juliet | 1.31 | 4.2 | 0.85 | 3.08 | 10.9 | 18.6 | 10.7 | 33.4 | 50.2 | 200 | 555 | 755
JulietDMO | 1.59 | - | 1.0 | 2.10 | 6.12 | 12.6 | 1.07 | 4.43 | 11.5 | 3.26 | 18.6 | 60.3
Romeo | 0.99 | 7.4 | 1.0 | 4.08 | 14.4 | 25.2 | 10.2 | 35.7 | 60.5 | 136 | 472 | 780
RomeoDMO | 1.26 | - | 1.2 | 1.74 | 5.73 | 12.2 | 0.881 | 3.83 | 11.3 | 2.56 | 14.7 | 61.2
Thelma | 1.17 | 7.9 | 1.1 | 0.304 | 2.18 | 8.86 | 1.22 | 8.12 | 26.7 | 26.9 | 188 | 449
ThelmaDMO | 1.70 | - | 1.3 | 0.721 | 2.52 | 6.69 | 0.310 | 1.62 | 6.87 | 0.825 | 6.31 | 42.3
Louise | 1.42 | 2.9 | 0.85 | 0.811 | 4.39 | 11.3 | 2.00 | 10.0 | 23.6 | 26.6 | 124 | 267
LouiseDMO | 1.41 | - | 1.0 | 1.01 | 3.77 | 10.4 | 0.557 | 2.85 | 10.3 | 0.825 | 6.30 | 42.3
Romulus | 1.00 | 10 | 1.53 | 7.30 | 18.4 | 28.6 | 26.8 | 63.7 | 92.2 | 529 | 1182 | 1610
RomulusDMO | 1.01 | - | 1.9 | 1.15 | 4.40 | 11.9 | 0.540 | 2.94 | 11.5 | 1.52 | 11.5 | 64.2
Remus | 1.10 | 5.1 | 0.97 | 1.67 | 7.48 | 16.6 | 4.65 | 19.5 | 39.0 | 69.4 | 272 | 500
RemusDMO | 1.18 | - | 1.3 | 1.57 | 5.45 | 13.3 | 0.814 | 3.87 | 12.5 | 2.53 | 16.0 | 68.9

##  2 Overview of Simulations

Our analysis relies on cosmological zoom-in simulations performed as part of
the Feedback In Realistic Environments (FIRE)
project111https://fire.northwestern.edu/ with FIRE-2 feedback implementation
(Hopkins et al., 2018b) with the gravity plus hydrodynamics code GIZMO
(Hopkins, 2015). FIRE-2 includes radiative heating and cooling for gas with
temperatures ranging from 10 K to $10^{10}\,{\rm K}$, an ionising background
(Faucher-Giguère et al., 2009), stellar feedback from OB stars, AGB mass-loss,
type Ia and type II supernovae, photoelectric heating, and radiation pressure.
Star formation occurs in gas that is locally self-gravitating, sufficiently
dense ($>1000$ cm-3), Jeans unstable, and molecular (following Krumholz &
Gnedin 2011). Locally, the star formation efficiency is set to $100\%$ per
free-fall time, though the global efficiency of star formation within a giant-
molecular cloud (or across larger scales) is self-regulated by feedback to
$\sim$1-10% per free-fall time (Orr et al., 2018).

In this work, we analyse 12 Milky-Way-mass galaxies (Table 1). These zoom
simulations are initialised following the approach outlined in Oñorbe et al.
(2014) using the MUSIC code (Hahn & Abel, 2011). Six of these galaxies were
run as part of the Latte suite (Wetzel et al., 2016; Garrison-Kimmel et al.,
2017b, 2019a; Samuel et al., 2020; Hopkins et al., 2018a) and have names
following the convention m12*. The other six, with names associated with
famous duos, are set in paired configurations to mimic the Milky Way and M31
(Garrison-Kimmel et al., 2019a, b). Analysis has shown these are good
candidates for comparison with the Milky Way (Sanderson et al., 2020). Gas
particles for the M12* runs have initial masses of $m_{\rm g,i}=7070\,$
$M_{\odot}$. The ELVIS on FIRE simulations have roughly two times better mass
resolution ($m_{\rm g,i}\simeq 3500-4000$ $M_{\odot}$). Gas softening lengths
are fully adaptive down to $\simeq$0.5$-$1 pc. The dark matter particle masses
are $m_{\rm dm}=3.5\times 10^{4}$ $M_{\odot}$for the Latte simulations and
$m_{\rm dm}\simeq 2\times 10^{4}$ $M_{\odot}$for the ELVIS runs. Star particle
softening lengths are $\simeq$4 pc physical and a dark matter force softening
is $\simeq$40 pc physical.

Lastly, each FIRE simulation has an analogous dark matter only (DMO) version.
The individual dark matter particle masses in the DMO simulations are larger
by a factor of $(1-f_{\rm b})^{-1}$ in order to keep the total gravitating
mass of the Universe the same, where $f_{\rm b}=\Omega_{\rm b}/\Omega_{\rm m}$
is the cosmic baryon fraction. The initial conditions are otherwise identical.
DMO versions of each halo are referred to with the same name as the FIRE
version with the added suffix “DMO."

As can be seen in Table 1, the stellar masses of the main galaxy in each FIRE
run (second column) are broadly in line with the Milky Way:
$M_{\star}\approx(3-11)\times 10^{10}$ $M_{\odot}$. The virial masses (Bryan &
Norman, 1998) of all the halos in these simulation span a range generally in
line with expectations for the Milky Way: $M_{\rm vir}\approx(0.9-1.8)\times
10^{12}$ $M_{\odot}$. In every case, the DMO version of each pair ends up with
a higher virial mass. This is consistent with the expectation that halos will
have lost their share of cosmic mass by not retaining all baryons in
association with feedback. As we discuss in the next section, in our primary
analysis we re-normalize all halos (both FIRE and DMO runs) so that they have
the same “local" dark matter density at the Solar location (by the factor $f$
listed in the table).

Figure 1:  Upper Left: Simulated dark matter halo density profiles for all DMO
(dashed) and FIRE (solid) runs (solid). Upper Right: Three-dimensional dark
matter velocity dispersion profiles for all DMO (dashed) and FIRE (solid)
simulations. Lower Left: Dark matter density profiles after re-normalizing.
The mass per particle in each simulation has been multiplied by a factor $f$
(ranging from $f=0.89-1.8$, see Table 1) to give them the same density at a
mock solar location: $\rho(r=R_{\odot})=10^{7}$ $M_{\odot}$${\rm kpc}^{-3}$
with $R_{\odot}=8.3$ kpc (vertical dotted line). Bottom Right: Velocity
dispersion profiles for each halo after re-scaling the particle velocities by
a factor $\sqrt{f}$. This roughly accounts for the re-scaling of the
mass/density profile. Note that after re-scaling, the DMO velocity dispersion
profiles become similar for $r\lesssim R_{\odot}$. Even when normalized at the
solar radius, there is almost an order of magnitude scatter in the inner
($\sim 400$pc) density for the FIRE simulations and they all have higher
central dark matter velocities than would have been expected from DMO
simulations.

Figure 2: Left: The ratio of the FIRE dark matter density profiles shown in
Figure 1 to its DMO counterpart. We see that in some cases the FIRE runs are
less dense than the DMO runs though most halos get denser at small radius in
response to galaxy formation. Right: The ratio of the FIRE dark matter halo
velocity dispersion profiles shown in Figure 1 divided by the DMO version for
each halo. We see that in every case, the process of galaxy formation heats
the dark matter at small radius.

###  2.1 Dark Matter Density and Velocity Dispersion Profiles

Lazar et al. (2020) provide an extensive discussion of the dark matter halo
density profiles for the simulations we analyze here. Every system has more
than $\sim 1000$ dark matter particles within the inner 400pc and is converged
outside of this radius according to the criteria discussed Hopkins et al.
(2018b). Though some systems are even better converged, for simplicity we
adopt the same convergence radius, $r_{\rm conv}=400$pc, for each halo and
only present values that depend on quantities outside of this radius. For an
adopted Solar location at radius $R_{\odot}=8.3$ kpc, our $400$ pc convergence
radius corresponds to an angle $\psi=2.75^{\circ}$ in projection at the
Galactic Center.

Figure 1 shows the spherically-averaged density and velocity dispersion
profiles of the simulations in our sample. The upper two panels show raw
simulation results, with differential density profiles on the left and
velocity-dispersion profiles on the right. The DMO simulations are in black
and FIRE runs are in blue. Note that the density profiles of the FIRE runs are
systematically steeper for $r\gtrsim 3$ kpc than the DMO runs. This is
consistent with the expectations that baryonic contraction makes halos more
concentrated at this stellar-mass scale (see Lazar et al., 2020). At smaller
radii ($r\lesssim 3$ kpc), however, the FIRE halos have a larger range of
central densities; sometimes feedback produces a core-like profile and
sometimes the halo remains fairly steep (see Mercado et al., 2021, for an
investigation into the origin of this variation).

The FIRE halos also have systematically higher central dark matter velocity
dispersion than the DMO runs. A similar result was reported by Robles et al.
(2019), which studied the velocity dispersion profiles of cold dark matter
halos using zoom in dark matter simulations that included a slowly-grown
Milky-Way disk potential. In these simulations, even without feedback, they
found that the central velocity dispersion of the dark matter was much higher
in runs with disk potentials compared to those without. Board et al. (2021)
also found that the central dark matter velocity dispersion was higher in
simulations that included full galaxy formation physics (though with a
different implementation than our own). Taken together, these results suggest
that the dark matter velocity dispersion at the center of Milky-Way mass halos
should be significantly higher than would be expected from DMO simulations,
irrespective of galaxy formation model.

The pair of panels on the bottom of Figure 1 show the profiles after we have
re-scaled them to the defaults we will use in the rest of the analysis. Our
aim here is to normalize each run to have the same local dark matter density
at the solar location $r=R_{\odot}$. We are motivated to do this because it is
customary in indirect-detection analyses to normalize the assumed profile at
$R_{\odot}$ and to marginalize about the local density range inferred by
observations. While our halos are Milky-Way like in virial mass and stellar
mass, they are not precise replicas of the Milky Way. By re-normalizing at the
solar location, our results become primarily about profile shape rather
normalization, and can be scaled appropriately as observational estimates of
the local density become more precise. We assume $R_{\odot}=8.3$ kpc (vertical
dotted lines in the left panels) and set the density there to be
$\rho(r=R_{\odot})=10^{7}$ $M_{\odot}$${\rm kpc}^{-3}=0.38$ GeV cm-3 (Guo et
al., 2020). We do this by scaling the particle masses in each simulation (post
process) by a factor $f$: m${}_{\rm dm}\rightarrow f\,$mdm. The values of $f$
for each simulation are given in Table 1 and range from $f=1.8$ to $f=0.89$.
We also re-scale the particle velocities in each simulation by a factor222This
assumes $v\propto\sqrt{fM/r}$. We have checked that the dark matter velocity
dispersion in our simulations does roughly scale with the local dark matter
density as $\sqrt{f}$. $\sqrt{f}$ in order to roughly account for the re-
scaling of the total mass: $v\rightarrow\sqrt{f}v$. Note that after re-
scaling, the DMO velocity dispersion profiles become similar for $r\lesssim
R_{\odot}$, as expected.

Even when normalized at the Solar radius, there is considerable scatter in the
inner density and the FIRE simulations display more variance than the DMO
simulations. The difference between DMO and FIRE is most systematic in the
inner velocity dispersion (bottom right of Figure 1). While the normalized DMO
simulations all have $\sigma\simeq 140\hskip 2.84544pt{\rm km\,s}^{-1}$ at
$r=400$ pc, the FIRE runs have $\sigma\simeq 350-550\,{\rm km\,s}^{-1}$ at the
same radius. While this scatter is interesting to note, giving a precise and
detailed answer as to why it occurs while require more analysis and is the
topic for future analysis. For this present paper, we note that it exists, and
that is has a significant impact on the magnitude of the J-factor signal, so
that the total magnitude of the J-factor signal varies quite signficantly from
halo to halo.

Figure 2 shows FIRE to DMO ratios for the normalized density profile (left)
and velocity dispersion profile (right) of each halo pair. On the left we see
that galaxy formation has generally made the halos less dense at large radii,
corresponding to steeper (contracted) density profile. The effect of galaxy
formation on the inner density is quite varied, with some systems (e.g.
Romulus and Romeo) displaying higher central densities in the FIRE runs, while
others (e.g. Thelma and m12m) have lower densities. There is no clear trend
with stellar mass or virial mass associated with these differences. Thelma and
Romeo, for example, have very similar galaxy masses and virial masses but
galaxy formation seems to have had an opposite effect on their relative
density profiles.This is likely an artifact of some of the important baryonic
differences between halos at late times, which are studied in other papers:
some have late-occurring minor mergers or strong stellar bars (which tend to
push DM outwards and lower central densities; e.g. Sanderson et al., 2017;
Debattista et al., 2019), others have strong torques or early multiple-mergers
which produce inflows and dense bulges and more compact disks (e.g. Garrison-
Kimmel et al., 2018; Ma et al., 2017). For the purposes of this paper, it is
noted that it impacts the calculated J-factor and affects variance from halo
to halo.

The right panel of Figure 2 shows again that the effect of galaxy formation on
the dark matter velocity dispersion is systematic. In every case the FIRE runs
are hotter, with $\sim 3-4$ times higher velocity dispersion than their DMO
counterparts at $r=400$ pc. In Appendix B we show that these halos become
baryon dominated within 3-8 kpc from the center (Figure 2. As discussed next,
this enhancement in central velocity dispersion has a systematic effect on the
dark-matter annihilation J-factors for velocity-dependent cross sections.

##  3 Astrophysical J-Factors

###  3.1 Definitions

If dark matter particle of mass $m_{\chi}$ is its own antiparticle with an
annihilation cross section $\sigma_{A}$, the resulting differential particle
flux produced by annihilation in a dark matter halo can be written as the
integral along a line of sight $\ell$ from the observer (located at the solar
location in our case) in a direction $\vec{\theta}$ in the plane of the sky
over pairs of dark matter particles with velocities $\vec{v}_{1}$ and
$\vec{v}_{2}$:

$\frac{d^{2}\Phi}{dEd\Omega}=\frac{1}{4\pi}\frac{dN}{dE}\int{\rm{d}}\ell\leavevmode\nobreak\
{\rm{d}}^{3}v_{1}\,{\rm{d}}^{3}v_{2}\frac{f(\vec{r},\vec{v}_{1})}{m_{\chi}}\frac{f(\vec{r},\vec{v}_{2})}{m_{\chi}}\,\frac{(\sigma_{A}v_{\rm
rel})}{2}\ .$ (1)

Here, $\vec{r}=\vec{r}(\ell,\vec{\theta})$ is the 3D position, which depends
on the distance along the line of sight $\ell$ and sky location
$\vec{\theta}$. The dark matter velocity distribution $f(\vec{r},\vec{v})$ is
normalized such that $\int d^{3}vf(\vec{r},\vec{v})=\rho(\vec{r})$, where
$\rho$ is the dark matter density at that location. The symbol $v_{\rm
rel}=|\vec{v}_{1}-\vec{v}_{2}|$ represents the relative velocity between pairs
of dark matter particles. The quantity $m_{\chi}$ is the dark matter particle
mass and $dN/dE$ is the particle energy spectrum ultimately produced by a
single annihilation.

Following Boddy et al. (2018), we parameterize the velocity-dependence of the
dark matter annihilation cross section as

$\sigma_{A}v_{\rm rel}=[\sigma v]_{0}\,Q(v_{\rm rel}),$ (2)

where $[\sigma v]_{0}$ is the overall amplitude and the function $Q(v)$
parameterizes the velocity dependence. For $s$-wave, $p$-wave, and $d$-wave
annihilation, $Q(v)=1$, $(v/c)^{2}$, and $(v/c)^{4}$, respectively. We can
then rewrite the differential particle flux as

$\frac{d^{2}\Phi}{dE{\rm{d}}\Omega}=\frac{(\sigma_{A}v)_{0}}{8\pi
m_{X}^{2}}\frac{dN}{dE_{\gamma}}\left[\frac{dJ_{Q}}{d\Omega}\right]\,.$ (3)

Here, the term in brackets absorbs all of the astrophysics inputs and defines
the astrophysical "J-factor"

$\frac{dJ_{Q}}{d\Omega}(\vec{\theta})=\int{\rm{d}}\ell\int{\rm{d}}^{3}v_{1}{f(\vec{r},\vec{v}_{1})}\int{\rm{d}}^{3}v_{2}{f(\vec{r},\vec{v}_{2})}\,Q(v_{\rm
rel})\,.$ (4)

In principle, the $\ell$ integral above sums pairs along the line-of-sight
from the observer ($\ell=0$) to infinity. In practice, we are focusing on
J-factors arising from an individual “Milky Way" halo, and truncate our
integrals at the halo’s edge (see below).

It is often useful to quote the cumulative J-factor within a circular patch of
sky of angular radius $\psi$ centered on the Galactic Center. In this case,
the patch defined by $\psi$ subtends a solid angle
$\Omega_{\psi}=4\pi\sin^{2}({\psi}/2)$ and we have:

$J_{Q}(<\psi)=\int_{0}^{\Omega_{\psi}}\frac{dJ_{Q}}{d\Omega}(\vec{\theta})\,{\rm{d}}\Omega\,.$
(5)

Figure 3: All-sky Hammer projections of J-factors ($dJ/d\Omega$) for JulietDMO
(left) and Juliet (right) in Galactic coordinates as viewed from mock solar
locations 8.3 kpc from the halo centers. Maps utilize bins of roughly 1.3
square degrees on the sky. From top to bottom, the rows assume s-wave, p-wave,
and d-wave annihilation. The color map in each pair of panels is fixed for
each type of assumed annihilation and is logarithmic in $dJ/d\Omega$, as
indicated by the bar along the top of each image. Note that FIRE runs (right)
produce systematically rounder J-factor maps on the sky. The p-wave and d-wave
maps are significantly brighter and more extended for the FIRE runs as well,
owing to the effects of galaxy formation in enhancing dark matter velocity
dispersion in the center of each halo.

Figure 4: Same as figure 3, but for M12cDMO (left) and M12c (right).

###  3.2 Approach

In what follows we aim to determine the astrophysical J-factors for each of
our simulated halos for s-wave, p-wave, and d-wave annihilation. In doing so
we approximate the dark matter distribution $f(\vec{r},\vec{v})$ as a
separable function:

$f(\vec{r},\vec{v}_{1})=\rho(\vec{r})\,g(\vec{v}(\vec{r})),$ (6)

with the dark matter density $\rho$ estimated using direct particle counts in
the simulation. In this estimate, we everywhere rely on the nearest 32 local
particles to determine the density. 333 As discussed in Board et al. (2021),
the inclusion of subhalos that can be resolved at this resolution is not
expected to increase global J-factors by more than $\sim 20\%$. We assume that
$g(\vec{v}(\vec{r}))$ at location $\vec{r}$ can be approximated as a local
Maxwellian, $g_{r}$. Board et al. (2021) have shown that the distribution of
velocities in their simulations are well approximated by this. Assuming a
local 3D velocity dispersion $\sigma=\sigma(r)$ we can write

$g_{r}(v){\rm{d}}^{3}v=\left(\frac{3}{2\pi\sigma^{2}(r)}\right)^{3/2}\,\exp{\left[-\frac{3v^{2}}{2\sigma^{2}(r)}\right]}\,{\rm{d}}^{3}v.$
(7)

We measure $\sigma$ using direct particle counts at a given radius in the
simulation. The Maxwellian assumption gives us
$\int{\rm{d}}^{3}v\,g_{r}(v)=1.0$, $\langle v^{2}\rangle=\sigma^{2}$, and
$\langle v^{4}\rangle=15/9\,\sigma^{4}$.

For standard $s$-wave annihilation, $Q(v)=1$, and the effective J-factor (Eq.
4) reduces to a simple integral over the density squared:

$\displaystyle\frac{dJ_{s}}{d\Omega}(\vec{\theta})$ $\displaystyle=$
$\displaystyle\int{\rm{d}}\ell\,\rho^{2}(\vec{r})\int{\rm{d}}^{3}v_{1}\,g_{r}(v_{1})\int{\rm{d}}^{3}v_{2}\,g_{r}(v_{2})$
(8) $\displaystyle=$ $\displaystyle\int{\rm{d}}\ell\,\rho^{2}[\ell(r)].$

For $p$-wave annihilation, $Q(v)=(v/c)^{2}$, and Eq. 4 becomes

$\displaystyle\frac{dJ_{p}}{d\Omega}(\vec{\theta})$ $\displaystyle=$
$\displaystyle\int{\rm{d}}\ell\,\rho^{2}(\vec{r})\int{\rm{d}}^{3}v_{1}\,g_{r}(v_{1})\int{\rm{d}}^{3}v_{2}\,g_{r}(v_{2})\,\frac{|\vec{v}_{1}-\vec{v_{2}}|^{2}}{c^{2}}$
(9) $\displaystyle=$
$\displaystyle\frac{2}{c^{2}}\int{\rm{d}}\ell\,\rho^{2}[\ell(r)]\,\sigma^{2}[\ell(r)].$

In the second line we have used
$\langle(\vec{v_{1}}-\vec{v_{2}})^{2}\rangle=2\langle v^{2}\rangle$, under the
assumption that the cross-terms vanish via local spherical symmetry. For
$d$-wave annihilation, $Q(v)=(v/c)^{4}$, which implies

$\displaystyle\frac{dJ_{d}}{d\Omega}(\vec{\theta})$ $\displaystyle=$
$\displaystyle\int{\rm{d}}\ell\,\rho^{2}(\vec{r})\int{\rm{d}}^{3}v_{1}\,g_{r}(v_{1})\int{\rm{d}}^{3}v_{2}\,g_{r}(v_{2})\,\frac{|\vec{v}_{1}-\vec{v_{2}}|^{4}}{c^{4}}$
(10) $\displaystyle=$
$\displaystyle\frac{48}{9c^{4}}\int{\rm{d}}\ell\,\rho^{2}[\ell(r)]\,\sigma^{4}[\ell(r)].$

The co-factor $48/9$ comes from
$\langle(\vec{v_{1}}-\vec{v_{2}})^{4}\rangle=2\langle v^{4}\rangle+2\langle
v^{2}\rangle^{2}$, again assuming that the cross terms vanish and that we have
a Maxwellian distribution ($\langle v^{4}\rangle=15/9\,\sigma^{4}$).

###  3.3 Geometric Setup

For each halo in our sample, we calculate J-factors as defined in Equation 4,
integrating from a mock Solar location (setting $\ell=0$) to the edge of the
halo, which we define as a sphere of radius $r=300$ kpc from the center of
each halo in every case. While the virial radii (Bryan & Norman, 1998) of our
halos range from $300-335$ kpc, we fix $300$ kpc as the halo boundary for
consistency. Since most of the J-factor signal comes from the inner halo,
changing the outer radius by $10\%$ has no noticeable affect on our results.

For the DMO runs, we assume that the Galactic Center corresponds to the halo
center and fix the observer location to be at a distance 8.3 kpc from the halo
center along the x-axis of the simulation. For FIRE runs, we position the
observer in the galaxy disk plane at a radius of 8.3 kpc from the halo center.
We define the disk plane to be perpendicular to the angular momentum vector of
all the stars within 20 kpc of the central galaxy. We note that while we have
placed the observer within the disk for physical consistency, in practice
galaxy formation renders the dark matter halos quite spherical within $\sim
10$ kpc (Zhu et al., 2016), so the location of the observer in the disk plane
is not a major factor on the results.

Figure 5: Differential J-factor profiles (Equation 4) as a function of angle
$\psi$ from the Galactic center for s-wave, p-wave, and d-wave annihilation
(from left to right). The FIRE simulation halos are shown as solid blue lines
and their DMO counterparts are shown as dashed black lines. For the s-wave
case (left) we see that the central J-factor values are similar for FIRE and
DMO cases, though the FIRE profiles are flatter at small angle, and more
extended on the sky. The FIRE p-wave profiles (middle) are a noticeably
amplified compared to the DMO cases. Their shapes are also significantly
different – with a flatter inner profile and sharper fall-off at angles beyond
10 degrees. The d-wave case (right panel) demonstrates the starkest difference
between FIRE and DMO runs.

Figure 6: Ratio of FIRE to DMO J-factors for each halo as a function of angle
$\psi$ from the Galactic center for s-wave, p-wave, and d-wave annihilation
(from left to right). Each halo pair has a unique color, as indicated. For the
case of s-wave (left), the ratio is of order unity, ranging from a factor of
$\sim 3$ higher to a factor of $\sim 0.3$ lower at small angles. For p-wave,
the FIRE runs have J-factors as much as $\sim 30$ times higher, though the
amplification can be as small as a factor of $\sim 3$. In the case of d-wave,
amplification factors as large as $\sim 100-400$ at small angle are seen.

Figure 7:  Cumulative J-factor within annular angle $\psi$ from the Galactic
Center (Eq. 5) for DMO (left) and FIRE (right) for s-wave (top), p-wave
(middle), and d-wave (bottom) annihilation. For the s-wave case, FIRE runs
have J-factors of the same order of magnitude as the DMO runs, while for
p-wave and d-wave the J-factors are considerably higher.

Figure 8: Schematic illustration of how the cross section versus particle mass
constraints from Abazajian et al. (2020) (black solid lines, upper panels)
would shift for s-wave (top), p-wave (middle), and d-wave (bottom)
annihilation based on the relative J-factors (blue) from each of our DMO
(left) and FIRE runs (right). The region above the lines is ruled out. The
dotted lines show the $[\sigma v]_{0}$ required for thermal dark matter to
match the observed abundance. FIRE results suggest that current constraints
for p-wave are much closer to the thermal cross section than would have been
expected from DMO halos, potentially within a factor of $\sim 10$ for $10$ GeV
WIMPS.

##  4 Results

Figures 3 and 4 illustrate graphically our results for two example halo pairs,
Juliet and m12c, respectively. We show all-sky Hammer projection maps of
$dJ/d\Omega$ for s-wave (top), p-wave (middle), and d-wave (bottom) for the
DMO (left) and FIRE (right) run of each halo. We are viewing the Galactic
Center (middle of each image) from mock solar locations as defined in the
previous section. The color bars are are mapped to J-factor amplitude as
indicated at the top of each image. Note that every row has the same color
mapping, so that the relative difference between DMO and FIRE runs can be seen
clearly for each assumed velocity dependence. The binning in these maps is
$1.3$ square degrees.

The first takeaway from these images is that the FIRE runs are significantly
brighter (with amplified J-factors) than the DMO runs for the p-wave and
d-wave cases. This is a direct result of the FIRE halos having enhanced dark
matter velocity dispersion compared to the DMO halos (e.g., the right panel of
Figure 2). The maps are also more extended from the Galactic Center. The
s-factor maps are not as different, given the modest differences in central
densities for these particular halos (see the left panel of Figure 2), though
substructure is significantly reduced in the FIRE runs, as expected from the
destructive effects of the central galaxy (Garrison-Kimmel et al., 2017a;
Kelley et al., 2019). Note that in all cases, including s-wave, the FIRE maps
are rounder on the sky – this is a result of galaxy formation tending to
sphericalize the dark matter distributions compared to DMO runs in halo
centers (Debattista et al., 2008; Bernal et al., 2016; Chua et al., 2019;
Kelley et al., 2019; Shen et al., 2021; Sameie et al., 2021). The fact that we
expect annihilation signals to be even rounder than in the DMO case should in
principle make it easier to detect or exclude dark matter annihilation in the
face of astrophysical backgrounds, which are expected to track more closely
the shape of the Galaxy (e.g. Abazajian et al., 2020). We will not focus on
quantifying this difference in shape on the sky here because it will be the
subject of future work. We will instead focus on azimuthally-averaged results
in what follows.

Figure 5 provides a summary of J-factor results for all of our FIRE halos
(solid blue) and DMO halos (dashed black). Plotted are $dJ/d\Omega$ profiles
(Equation 4) as a function of angle $\psi$ with respect to the Galactic
Center. Results for s-wave, p-wave, and d-wave are shown in separate panels,
from left to right. As expected, the FIRE runs are generally amplified
compared to DMO runs, especially for the p-wave and s-wave cases. The shapes
of the profiles are also significantly different in character. While the DMO
runs show a trend for J-factor profiles to be more peaked at small angle for
s-wave, and to become flatter and more extended on the sky as we progress to
p-wave and d-wave, the FIRE profiles are more similar in shape. In all cases
(s-wave, p-wave, and d-wave) the "emission" profile is fairly constant out to
$\sim 5-10$ degrees in the FIRE runs, with a steep fall-off towards larger
angles beyond that point.

Figure 6 shows, for each halo pair, the ratio of $dJ/d\Omega$ in FIRE to the
DMO case as a function of angle from the Galactic Center, $\psi$. Each halo
pair has a unique color, as indicated. For s-wave annihilation (left panel) we
see that the FIRE runs sometime produce higher J-factors (up to a factor of
$\sim 6$) at small angle and sometimes give decreased J-factors (as small as
$\sim 0.3$ of the DMO value). Because for s-wave annihilation the J-factor
depends only on the density, this behavior tracks what seen for the density
profiles (Figure 2). Sometimes feedback has produced a cored-out central
density profile, leading to a lower central J-factor; sometimes baryonic
contraction is more important, and this creates higher central densities and
higher J-factors at small angle.

The middle and right panels of Figure 6 show that the J-factor ratios are
higher, in all cases, for the FIRE runs for p-wave and d-wave annihilation.
This is because the dark matter velocities are always enhanced (see Figure 1)
by factors of $2.5-4$, which is enough to boost the p-wave and d-wave
annihilation with respect to DMO runs, even in the cases where the central
density is slightly smaller. Typical amplification factors are $\sim 10$ for
p-wave and $\sim 100$ d-wave. In cases like M12i and Romulus, where the dark
matter density is also higher in the FIRE runs, the p-wave and d-wave ratios
can very large (factors of $\sim 40-60$ and $\sim 400-500$, respectively).

Another way to see the difference between DMO and FIRE models is to compare
the cumulative J-factors within annular angle $\psi$ from the Galactic Center.
Figure 7 shows this quantity for each DMO (left) and FIRE (right) run. Each
halo pair has a unique color, as indicated. For s-wave annihilation (top left
and right panels) we see that the FIRE runs generally do have larger
integrated J-factors within 20 degrees of the Galactic Center, even in cases
(like Thelma) that have somewhat smaller signals within 10 degrees. The middle
and bottom panels of Figure 7 show that the J-factor cumulative totals are
higher, in all cases, for the FIRE runs for p-wave and d-wave annihilation.
Though the cumulative totals are more enhanced within $\sim 10$ degrees owing
the the same of the profile.

In the next section we briefly explore implications of our results for dark
matter indirect detection.

##  5 Implications

One of our primary results is that the DM velocities in our full galaxy
formation runs are significantly higher than would be expected from DMO runs;
this elevates the expected signal for fixed cross section in p-wave and d-wave
models (see, e.g. Figure 7). In what follows we aim to provide a schematic
illustration of how our results may impact attempts to constrain dark matter
models with thermal abundance cross sections, especially those with velocity
dependence. We use the published results of Abazajian et al. (2020) in this
illustrative example, and adopt a simple scaling of their published limits to
provide a first-order sketch of how our results may impact future attempts to
constrain s, p, and d-wave annihilation.

As discussed in the introduction, the realization that the observed Galactic
Center gamma-ray excess has a non-circular, boxy shape that traces the
Galactic Bulge (Macias et al., 2018; Bartels et al., 2018) has allowed
Abazajian et al. (2020) to rule out a number of thermal-abundance WIMP models
with s-wave annihilation channels for $m_{\chi}\lesssim 500$ GeV. Such a
result motivates the exploration of p-wave (and d-wave models). There is a
general expectation that velocity-suppressed p-wave and d-wave annihilation
will be far from detectable in Milky Way (though see Johnson et al., 2019).
This is because typical DM velocities in the Galactic Center are usually
thought to be $v\sim 100$ km s-1 (based on DMO simulations), compared values
$\sim 10^{3}$ times higher during thermal freeze out. While, from the point of
view of a model builder, such a suppression is “good" because it evades
direct-detection bounds, from the point of view of an observer or
experimentalist, this level of suppression is a potential nightmare: how can
we detect such a signal?

The thick black lines in the upper set of panels Figure 8 reproduces the
s-wave constraints published by Abazajian et al. (2020). The horizontal axis
shows the WIMP mass and vertical axis is the velocity-averaged cross section.
In our generalized language, the vertical axis specifically corresponds to the
normalization $[\sigma v]_{0}$, defined by $\langle\sigma v\rangle=[\sigma
v]_{0}Q(v_{\rm rel})$ in Equation 11, where $Q(v)=1$ is the s-wave case. Cross
sections above the black line are excluded. In deriving this constraint, they
assumed a $b\bar{b}$ annihilation channel and a plausible range of Milky Way
dark matter profiles (their "NFW" case) as expected from DMO simulations. The
dashed line shows the required cross section to produce the correct thermal
abundance of dark matter observed (Steigman et al., 2012).

The blue lines in Figure 8 provide schematic estimates for how the Abazajian
et al. (2020) limit would shift for s- (top), p- (middle), and d-wave (bottom)
annihilation for halos that match our simulation results. Here we have made
the simplistic assumption that limit will scale in direct proportion to the
integrated J-factor within $10^{\circ}$ of the Galactic Center. The range of
NFW profiles considered in Abazajian et al. (2020) have central densities
quite similar to our own M12wDMO case, and we use this to set the reference
J-factor for the constraint: $J_{s}(<10^{\circ})\equiv J_{\rm ref}=1.7\times
10^{22}$ GeV2 cm-3. Note that this reference Js-factor is at the lower range
of those from our DMO halos (See Table 1). This is mainly due to the fact that
while we normalized each halo to a local, solar dark matter density of $0.38$
GeV cm-3, they assumed a median normalization of $0.28$ GeV cm-3 at the solar
radius. Some of our DMO halos also deviate somewhat from a strict NFW shape,
effectively making them slightly denser at $\sim 1$ kpc than that shape would
predict for a fixed solar normalization (e.g. Lazar et al., 2020).

Each blue line is a scaled version of the black line. Specifically, for each
halo in our suite, we determine the ratio of the reference J-factor from
Abazajian et al. to the measured J-factor for s, p, and d-wave cases: $J_{\rm
ref}/J_{q}(<10^{\circ})$ (for $q=$ s, p, and d), and multiply the Abazajian
limit by that ratio to estimate the implied limit.

The top pair of panels shows how the implied limit scales for each of our DMO
runs (left) and FIRE runs (right) in the s-wave case. As mentioned above, even
in the DMO case, our halos tend to have larger J-factors than the halos used
in Abazajian et al. (2020) because of our chosen local-density normalization.
The spread in lines comes about because of the halo-to-halo scatter.
Interestingly, for the FIRE cases, which we regard as more realistic, the
limit lies above all lines, though is within $\sim 25\%$ of the upper
envelope. One way to interpret this is that the Abazajian limit is
conservative in comparison to our expectations, but not exceedingly so,
especially considering the sensitivity to local dark matter density
normalization. In this sense, our results are unlikely to affect current
constraints for the s-wave cross section significantly.

The middle panels show the implied constraints for the case of p-wave
annihilation, with DMO halos on the left and FIRE halos on the right. The
dotted line shows the required cross section normalization for the thermal
abundance in the p-wave annihilation case. The value is about $50$ times
higher than the s-wave thermal abundance normalization to make up for the fact
that the total cross section is suppressed by a factor $Q=(v/c)^{2}$ with
$v/c\sim 0.15$ during freeze-out. We work out the thermal abundance
normalization explicitly for p-wave and d-wave dark matter in Appendix A.

Note that the scaled limits in the DMO p-wave panel are more than two orders
of magnitude above the thermal cross section, which would suggest indirect
detection for such a model is unlikely. The right middle panel of Figure 8
unitizes what we believe to be more realistic p-wave J-factors from our FIRE
runs. The blue lines in this case suggest that a much more powerful constraint
is possible for p-wave annihilation than would have been expected from DMO
halos alone. In particular, we see that for low-mass WIMPS, we may be within a
factor of $\sim 10-20$ of detecting a p-wave annihilation signal if the Milky
Way resembles halos like Romulus, Juliet, and Romeo, which have among our
highest p-wave J-factors.

Finally, the bottom panels show the implied constraints for the case of d-wave
annihilation. The required cross section normalization for the thermal
abundance (dotted line) in the d-wave annihilation case is $\sim 2000$ times
higher than the s-wave thermal abundance normalization (Appendix A), though it
is still orders of magnitude below any of the scaled constraints. Inferred
constraints from the DMO halos (left) are five to six orders of magnitude
above the thermal abundance normalization. For the FIRE cases, the situation
is slightly better (roughly four orders of magnitude) though still far out of
reach.

Of course, realistic constraints will require a careful analysis of Fermi-LAT
Galactic Center observations, including templates for the stellar galactic and
nuclear bulges, variations in the Galactic diffuse emission models, and a
careful consideration of the shape of p-wave J-factor models of the kind show
in Figures 3 and 4, which tend to be even less boxy in the hydrodynamic runs
than would be expected in DMO. Based on the rough estimates presented here,
such an analysis is certainly warranted.

##  6 Conclusions

We have explored how galaxy formation affects predictions for the
astrophysical J-factors of Milky Way size dark matter halos. For a fixed
particle physics model, astrophysical J-factors are directly proportional to
the expected flux of Standard Model particles sourced by dark matter
annihilation, and therefore provide a crucial input for dark-matter indirect
detection searches in the Milky Way (see Eq. 4). In particular, we have used
twelve FIRE zoom simulations of Milky Way-type galaxies along with dark-
matter-only (DMO) versions of the same halos and worked out implications for
both velocity-independent (s-wave) and velocity-dependent (p-wave and d-wave)
annihilation cross sections.

One significant result is that the central dark matter velocity dispersion in
FIRE halos is systematically amplified by factors of $\sim 2.5-4$ compared to
their DMO counterparts (right panel of Figure 2). The effect of galaxy-
formation on the central dark matter density in the same halos is less
systematic, sometimes increasing and sometimes decreasing the central density,
with ratios ranging from $\sim 0.3-2.5$ (left panel of Figure 2). For p-wave
($\propto(v/c)^{2})$ and d-wave ($\propto(v/c)^{4}$) models, our FIRE-derived
J-factors are amplified by factors of $\sim 3-60$ and $\sim 10-500$ compared
to DMO runs (see Figure 6). FIRE halos generally produce J-factor profiles
that are flatter (less peaked) towards the Galactic Center (see Figure 5) and
rounder on the sky (see images 3 and 4). Note that these differences occur
despite the fact that we have normalized all of our halos to have the same
local (solar location) dark matter density. That is, these results are driven
by differences in the shape of the underlying dark matter density and velocity
dispersion profiles brought about by galaxy formation processes.

One basic implication of our results is that we expect p-wave and d-wave dark
matter annihilation to produce more easily detectable signals than would have
been expected from DMO halos. For example, while it is typical to suspect that
p-wave annihilation ($\propto(v/c)^{2})$) is suppressed to undetectable levels
444Though see Johnson et al. (2019), who have investigated the role of the
central black hole in altering the dark matter velocity dispersion. in the
Milky Way today (where $v\ll c$), we showed in section 5 that this may not be
the case. With the amplified velocities we see in our FIRE runs, the detection
of (or interesting constraints on) thermal-relic p-wave dark matter may not be
too far out of reach. In particular, by scaling the s-wave constraints from
Fermi-LAT derived by Abazajian et al. (2020), we showed that a similar
analysis could bring p-wave constraints to within a factor of $\sim 10$ of the
naive thermal cross section (see the right middle panel of Figure 8). Future
analyses that include detailed simulation-inspired priors on the shape of the
annihilation signal in these models, could potentially approach the thermal
value.

Another result worth highlighting is that we see significant scatter in the
density profiles (and associated s-wave J-factor profiles) in our FIRE runs.
As we begin to explore joint-constraints from multiple galaxies (e.g. M31 and
the Milky Way) it will be important to allow for the possibility that even
halos with similar halo masses are expected to have a large scatter in
J-factor normalization. Figure 5, for example, shows that the variance in our
FIRE s-wave J-factors is larger than one order of magnitude for our sample at
1 degree from the Galactic Center, despite the fact that our halos have been
fixed to have the same local density of dark matter at the solar location.

## Acknowledgements

We dedicate this paper to the memory of our dear friend and colleague José
Antonio Florez Velázquez. We would like to thank Tyler Kelly for helpful
suggestions and assistance in the analysis. We thank Kev Abazajian and Louis
Strigari for useful discussions.

DM, JSB, and FJM were supported by NSF grant AST-1910346. ZH is supported by a
Gary A. McCue postdoctoral fellowship at UC Irvine. PH is supported by NSF
Research Grants 1911233 & 20009234, NSF CAREER grant 1455342, NASA grants
80NSSC18K0562, HST-AR-15800.001-A. Numerical calculations were run on the
Caltech compute cluster “Wheeler,” allocations FTA-Hopkins/AST20016 supported
by the NSF and TACC, and NASA HEC SMD-16-7592.

AW received support from: NSF grants CAREER 2045928 and 2107772; NASA ATP
grant 80NSSC20K0513; HST grants AR-15809 and GO-15902 from STScI; a Scialog
Award from the Heising-Simons Foundation; and a Hellman Fellowship.

MBK acknowledges support from NSF CAREER award AST-1752913, NSF grant
AST-1910346, NASA grant NNX17AG29G, and HST-AR-15006, HST-AR-15809, HST-
GO-15658, HST-GO-15901, HST-GO-15902, HST-AR-16159, and HST-GO-16226 from the
Space Telescope Science Institute, which is operated by AURA, Inc., under NASA
contract NAS5-26555.

## Data Availability

The data supporting the plots within this article are available on reasonable
request to the corresponding author. A public version of the GIZMO code is
available at http://www.tapir.caltech.edu/ phopkins/Site/GIZMO.html.

## References

  * Abazajian (2011) Abazajian K. N., 2011, J. Cosmology Astropart. Phys., 2011, 010
  * Abazajian & Kaplinghat (2012) Abazajian K. N., Kaplinghat M., 2012, Phys. Rev. D, 86, 083511
  * Abazajian et al. (2020) Abazajian K. N., Horiuchi S., Kaplinghat M., Keeley R. E., Macias O., 2020, Phys. Rev. D, 102, 043012
  * Ackermann et al. (2017) Ackermann M., et al., 2017, ApJ, 840, 43
  * Ajello et al. (2016) Ajello M., et al., 2016, ApJ, 819, 44
  * Argüelles et al. (2019) Argüelles C. A., Diaz A., Kheirandish A., Olivares-Del-Campo A., Safa I., Vincent A. C., 2019, arXiv e-prints, p. arXiv:1912.09486
  * Bartels et al. (2016) Bartels R., Krishnamurthy S., Weniger C., 2016, Phys. Rev. Lett., 116, 051102
  * Bartels et al. (2018) Bartels R., Storm E., Weniger C., Calore F., 2018, Nature Astronomy, 2, 819
  * Bergström et al. (1998) Bergström L., Ullio P., Buckley J. H., 1998, Astroparticle Physics, 9, 137
  * Bernal et al. (2016) Bernal N., Necib L., Slatyer T. R., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 030
  * Bertone & Hooper (2018) Bertone G., Hooper D., 2018, Reviews of Modern Physics, 90, 045002
  * Bertone et al. (2005) Bertone G., Hooper D., Silk J., 2005, Phys. Rep., 405, 279
  * Board et al. (2021) Board E., Bozorgnia N., Strigari L. E., Grand R. J. J., Fattahi A., Frenk C. S., Marinacci F., Navarro J. F., 2021, arXiv e-prints, p. arXiv:2101.06284
  * Boddy et al. (2018) Boddy K. K., Kumar J., Strigari L. E., 2018, Phys. Rev. D, 98, 063012
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Carlson & Profumo (2014) Carlson E., Profumo S., 2014, Physical Review D, 90, 023015
  * Choquette et al. (2016) Choquette J., Cline J. M., Cornell J. M., 2016, Phys. Rev. D, 94, 015018
  * Chua et al. (2019) Chua K. T. E., Pillepich A., Vogelsberger M., Hernquist L., 2019, MNRAS, 484, 476
  * Daylan et al. (2016) Daylan T., Finkbeiner D. P., Hooper D., Linden T., Portillo S. K., Rodd N. L., Slatyer T. R., 2016, Physics of the Dark Universe, 12, 1
  * Debattista et al. (2008) Debattista V. P., Moore B., Quinn T., Kazantzidis S., Maas R., Mayer L., Read J., Stadel J., 2008, ApJ, 681, 1076
  * Debattista et al. (2019) Debattista V. P., Gonzalez O. A., Sanderson R. E., El-Badry K., Garrison-Kimmel S., Wetzel A., Faucher-Giguère C.-A., Hopkins P. F., 2019, MNRAS, 485, 5073
  * Diemand et al. (2007) Diemand J., Kuhlen M., Madau P., 2007, ApJ, 657, 262
  * Dodelson (2003) Dodelson S., 2003, Modern cosmology. Elsevier
  * Einasto et al. (1974) Einasto J., Kaasik A., Saar E., 1974, Nature, 250, 309
  * Faucher-Giguère et al. (2009) Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, 1416
  * Feng (2004) Feng J. L., 2004, arXiv e-prints, pp hep–ph/0405215
  * Garrison-Kimmel et al. (2017a) Garrison-Kimmel S., et al., 2017a, Monthly Notices of the Royal Astronomical Society, 471, 1709
  * Garrison-Kimmel et al. (2017b) Garrison-Kimmel S., et al., 2017b, MNRAS, 471, 1709
  * Garrison-Kimmel et al. (2018) Garrison-Kimmel S., et al., 2018, MNRAS, 481, 4133
  * Garrison-Kimmel et al. (2019a) Garrison-Kimmel S., et al., 2019a, MNRAS, 487, 1380
  * Garrison-Kimmel et al. (2019b) Garrison-Kimmel S., et al., 2019b, MNRAS, 489, 4574
  * Gaskins (2016) Gaskins J. M., 2016, Contemporary Physics, 57, 496
  * Giacchino et al. (2013) Giacchino F., Lopez-Honorez L., Tytgat M. H. G., 2013, J. Cosmology Astropart. Phys., 2013, 025
  * Grand & White (2021) Grand R. J. J., White S. D. M., 2021, MNRAS, 501, 3558
  * Guo et al. (2020) Guo R., Liu C., Mao S., Xue X.-X., Long R. J., Zhang L., 2020, Monthly Notices of the Royal Astronomical Society, 495, 4828
  * Hahn & Abel (2011) Hahn O., Abel T., 2011, MNRAS, 415, 2101
  * Hooper & Goodenough (2011) Hooper D., Goodenough L., 2011, Physics Letters B, 697, 412
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hopkins et al. (2018a) Hopkins P. F., et al., 2018a, Monthly Notices of the Royal Astronomical Society, 480, 800
  * Hopkins et al. (2018b) Hopkins P. F., et al., 2018b, MNRAS, 480, 800
  * Johnson et al. (2019) Johnson C., Caputo R., Karwin C., Murgia S., Ritz S., Shelton J., Collaboration F.-L., et al., 2019, Physical Review D, 99, 103007
  * Jungman et al. (1996) Jungman G., Kamionkowski M., Griest K., 1996, Phys. Rep., 267, 195
  * Kamionkowski et al. (2010) Kamionkowski M., Koushiappas S. M., Kuhlen M., 2010, Phys. Rev. D, 81, 043532
  * Karwin et al. (2020) Karwin C., Murgia S., Moskalenko I., Fillingham S., Burns A.-K., Fieg M., 2020, arXiv e-prints, p. arXiv:2010.08563
  * Kelley et al. (2019) Kelley T., Bullock J. S., Garrison-Kimmel S., Boylan-Kolchin M., Pawlowski M. S., Graus A. S., 2019, MNRAS, 487, 4409
  * Kolb (1990) Kolb T., 1990, The early universe. Addison-Wesley
  * Krumholz & Gnedin (2011) Krumholz M. R., Gnedin N. Y., 2011, ApJ, 729, 36
  * Kwan et al. (2017) Kwan J., et al., 2017, MNRAS, 464, 4045
  * Lazar et al. (2020) Lazar A., et al., 2020, Monthly Notices of the Royal Astronomical Society, 497, 2393
  * Lisanti (2017) Lisanti M., 2017, in Polchinski J., et al. eds, New Frontiers in Fields and Strings (TASI 2015. pp 399–446 (arXiv:1603.03797), doi:10.1142/9789813149441_0007
  * Ma et al. (2017) Ma X., Hopkins P. F., Wetzel A. R., Kirby E. N., Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Quataert E., 2017, MNRAS, 467, 2430
  * Macias et al. (2018) Macias O., Gordon C., Crocker R. M., Coleman B., Paterson D., Horiuchi S., Pohl M., 2018, Nature Astronomy, 2, 387
  * Mandelbaum et al. (2013) Mandelbaum R., Slosar A., Baldauf T., Seljak U., Hirata C. M., Nakajima R., Reyes R., Smith R. E., 2013, MNRAS, 432, 1544
  * Martin et al. (2014) Martin A., Shelton J., Unwin J., 2014, Phys. Rev. D, 90, 103513
  * Mercado et al. (2021) Mercado F. J., et al., 2021, Monthly Notices of the Royal Astronomical Society, 501, 5121
  * Mo et al. (2010) Mo H., Van den Bosch F., White S., 2010, Galaxy formation and evolution. Cambridge University Press
  * Murgia (2020) Murgia S., 2020, Annual Review of Nuclear and Particle Science, 70, 455
  * Necib et al. (2019) Necib L., Lisanti M., Garrison-Kimmel S., Wetzel A., Sanderson R., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2019, The Astrophysical Journal, 883, 27
  * Oñorbe et al. (2014) Oñorbe J., Garrison-Kimmel S., Maller A. H., Bullock J. S., Rocha M., Hahn O., 2014, MNRAS, 437, 1894
  * Orr et al. (2018) Orr M. E., et al., 2018, MNRAS, 478, 3653
  * Petač et al. (2018) Petač M., Ullio P., Valli M., 2018, J. Cosmology Astropart. Phys., 2018, 039
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Robertson & Zentner (2009) Robertson B. E., Zentner A. R., 2009, Phys. Rev. D, 79, 083525
  * Robles et al. (2019) Robles V. H., Kelley T., Bullock J. S., Kaplinghat M., 2019, MNRAS, 490, 2117
  * Rubin et al. (1978) Rubin V. C., Ford W. K. J., Thonnard N., 1978, ApJ, 225, L107
  * Sameie et al. (2021) Sameie O., et al., 2021, MNRAS, 507, 720
  * Samuel et al. (2020) Samuel J., et al., 2020, Monthly Notices of the Royal Astronomical Society, 491, 1471
  * Sanderson et al. (2017) Sanderson R., Wetzel A., Sharma S., Hopkins P., 2017, Galaxies, 5, 43
  * Sanderson et al. (2020) Sanderson R. E., et al., 2020, The Astrophysical Journal Supplement Series, 246, 6
  * Shen et al. (2021) Shen X., Hopkins P. F., Necib L., Jiang F., Boylan-Kolchin M., Wetzel A., 2021, MNRAS, 506, 4421
  * Springel et al. (2008) Springel V., et al., 2008, Nature, 456, 73
  * Steigman et al. (2012) Steigman G., Dasgupta B., Beacom J. F., 2012, Phys. Rev. D, 86, 023506
  * Stoehr et al. (2003) Stoehr F., White S. D. M., Springel V., Tormen G., Yoshida N., 2003, MNRAS, 345, 1313
  * Trimble (1987) Trimble V., 1987, ARA&A, 25, 425
  * Wetzel et al. (2016) Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, ApJ, 827, L23
  * Wittman et al. (2000) Wittman D. M., Tyson J. A., Kirkman D., Dell’Antonio I., Bernstein G., 2000, Nature, 405, 143
  * Zhu et al. (2016) Zhu Q., Marinacci F., Maji M., Li Y., Springel V., Hernquist L., 2016, Monthly Notices of the Royal Astronomical Society, 458, 1559

##  Appendix A Thermal cross section for velocity-dependent annihilation

In this subsection we provide an estimate for the required thermal cross
section to match the observed dark matter abundance for models with a
velocity-dependent annihilation cross section. We are assuming that

$\langle\sigma v\rangle=[\sigma v]_{0}Q(v_{\rm rel})$ (11)

where $[\sigma v]_{0}$ is a normalization constant and $Q(v)=1$,
$v^{2}/c^{2}$, and $v^{4}/c^{4}$, for s, p, and d-wave annihilation,
respectively. For the standard s-wave case ($Q=1$), the normalization required
for thermal abundance was worked out carefully by Steigman et al. (2012), who
find

$[\sigma v]_{0}^{T\,s}\simeq 2.3\times 10^{-26}\rm{cm}^{3}\,{\rm s}^{-1}$ (12)

for WIMP masses above $\sim 10$ GeV (Steigman et al., 2012). Our goal here is
to determine how the required normalization changes for p-wave and d-wave.

We will follow the textbook treatment of WIMP freeze-out (see, e.g. Kolb,
1990; Dodelson, 2003; Mo et al., 2010; Lisanti, 2017), with the standard
assumption that WIMPs are their own anti-particle and initially in thermal
equilibrium with Standard Model particles in the early universe with photon
temperature $T$. In this case the dark matter number density $n$ can be
tracked using a simplified version of the Boltzmann equation

$\frac{dn}{dt}+3H(t)n=\langle\sigma v\rangle(n_{eq}^{2}-n^{2}),$ (13)

where $H(t)$ is the Hubble parameter. Here, $n_{eq}$ is the equilibrium number
density in physical coordinates (Kolb, 1990), which scales with the expansion
factor $a$ as $n_{eq}\propto T^{3}\propto a^{-3}$ for relativistic species and
becomes thermally suppressed for non-relativistic species $\propto e^{-m/T}$,
where $m$ is the dark matter particle mass. It is useful to rewrite this
equation in terms of a new variable $Y\equiv n/T^{3}$, which scales out the
expansion of the universe and eliminates the $3H(t)$ term on the left-hand
side of Eq. 13:

$\frac{dY}{dt}=T^{3}\langle\sigma
v\rangle(Y_{eq}^{2}-Y^{2})\>\longrightarrow\>\frac{dY}{dx}=\frac{\lambda(x)}{x^{2}}(Y_{eq}^{2}-Y^{2}),$
(14)

where $Y_{eq}\equiv n_{eq}/T^{3}$ and the arrow points to an equivalent
equation that uses a new time variable $x=m/T$. The variable $\lambda(x)\equiv
m^{3}\langle\sigma v\rangle/H(m)$, where $H(m)$ in the denominator is the
Hubble parameter evaluated at the time when the temperature $T=m$. In the
standard (s-wave) treatment, $\lambda$ is a constant. More generally
$\lambda(x)\propto[\sigma v]_{0}\,x^{-n}$ with $n=0$, $1$, and $2$ for s-wave,
p-wave, and d-wave, respectively.

While there is no analytic solution to Equation 14, we can estimate how the
ultimate abundance will scale with input parameters. At early times ($x\ll
1$), the coefficient $\lambda/x^{2}$ will be very large and $Y$ will track the
equilibrium abundance. As time progresses, the $\lambda/x^{2}$ term will
become smaller than unity as the annihilation rate drops below the expansion
rate. After this point, the dark matter particles will no longer track the
equilibrium abundance but instead “freeze out". For typical models, the
freeze-out time is well into the non-relativistic regime $x_{f}\simeq 25$.
After this time ($x\gtrsim x_{f}$), we expect $Y\gg Y_{eq}$ because
$Y_{eq}\propto\exp(-x)$. In this limit, Equation 14 simplifies to

$\frac{dY}{dx}=\frac{-\lambda(x)}{x^{2}}Y^{2}.$ (15)

We can estimate $Y$ today by solving the above equation as an integral over
$x$ from the freeze-out time $x_{f}$ to today $x_{0}=m/T_{0}\simeq\infty$. The
value of $Y$ today, $Y_{\infty}$, maps directly to the dark matter density
today: $\Omega_{dm}\propto mY_{\infty}/\rho_{\rm crit}$. This allows us to
write

$\Omega_{dm}\propto\left[\int_{x_{f}}^{\infty}\frac{[\sigma
v]_{0}}{x^{n+2}}dx\right]^{-1}=\frac{(n+1)x_{f}^{n+1}}{[\sigma v]_{0}}.$ (16)

Equation 16 allows us to estimate how $[\sigma v]_{0}$ must change as a
function of velocity-dependence $n$ in order to keep $\Omega_{\rm dm}$ fixed:

$[\sigma v]_{0}\propto(n+1)x_{f}^{n+1}.$ (17)

Since we expect $x_{f}$ to vary only logarithmically with the cross section,
let us assume $x_{f}\simeq 25$ in all cases. According to Equation 17, scaling
the s-wave cross section normalization (12) by a factor $2x_{f}^{2}$ will give
us a thermal cross section estimate for p-wave ($n=1$):

$[\sigma v]_{0}^{T\,p}\approx 50[\sigma v]_{0}^{T\,s}=1.2\times
10^{-24}\rm{cm}^{3}\,{\rm s}^{-1}.$ (18)

For d-wave ($n=2$), the scaling is $3x_{f}^{3}$ and we have

$[\sigma v]_{0}^{T\,d}\approx 1875[\sigma v]_{0}^{T\,s}=4.3\times
10^{-23}\rm{cm}^{3}\,{\rm s}^{-1}.$ (19)

We have checked the above two numbers by solving Equation 14 numerically and
find them to be good approximations.

##  Appendix B Comparison of baryonic matter densities to dark matter in the
Galactic Centers of FIRE -2 halos

We examined the ratios of the densities of ordinary matter (gas and stars) to
that of dark matter in regions near the galactic center of FIRE -2 halos. Our
results are shown in Figure 9. As anticipated, the halos are dominated by
baryonic matter near their most central regions, and then fall below an even
1:1 at 3-8 kpc from the center. Interestingly, the ratios vary substatially
from halo to halo. For example, Louise and a few others have a baryonic
content roughly 10 times larger than their dark matter content within the
innermost resolved region (400 pc), but other halos, such as Thelma have
ratios of up to 20 to 30 times greater at the same radius. Understanding the
run-to-run variance in baryonic to dark matter density ratios will be a goal
of future work.

Figure 9: Ratio of baryonic density to dark matter density at as a function of
radius within FIRE halos. Note all halos are baryon dominated at inner radii
until falling to roughly 1:1 near a mock solar location $\sim 8$ kpc.