GCEshort = GCE, long = Galactic Centre Excess MCMCshort = MCMC, long = Markov
Chain Monte Carlo

# No evidence for $p$\- or $d$-wave dark matter annihilation from local large-
scale structure

A. Kostić<EMAIL_ADDRESS>Max Planck Institute for Astrophysics,
Karl-Schwarzschild-Straße 1, 85748 Garching, Germany  D. J. Bartlett
<EMAIL_ADDRESS>CNRS & Sorbonne Université, Institut d’Astrophysique
de Paris (IAP), UMR 7095, 98 bis bd Arago, F-75014 Paris, France
Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road,
Oxford, OX1 3RH, UK  H. Desmond  Institute of Cosmology & Gravitation,
University of Portsmouth, Dennis Sciama Building, Portsmouth, PO1 3FX, UK

###### Abstract

If dark matter annihilates into standard model particles with a cross-section
which is velocity dependent, then Local Group dwarf galaxies will not be the
best place to search for the resulting gamma ray emission. A greater flux
would be produced by more distant and massive halos, with larger velocity
dispersions. We construct full-sky predictions for the gamma-ray emission from
galaxy- and cluster-mass halos within $\sim 200{\rm\,Mpc}$ using a suite of
constrained $N$-body simulations (CSiBORG) based on the Bayesian Origin
Reconstruction from Galaxies algorithm. Comparing to observations from the
Fermi Large Area Telescope and marginalising over reconstruction uncertainties
and other astrophysical contributions to the flux, we obtain constraints on
the cross-section which are two (seven) orders of magnitude tighter than those
obtained from dwarf spheroidals for $p$-wave ($d$-wave) annihilation. We find
no evidence for either type of annihilation from dark matter particles with
masses in the range $m_{\chi}=2-500{\rm\,GeV}/c^{2}$, for any channel. As an
example, for annihilations producing bottom quarks with
$m_{\chi}=10{\rm\,GeV}/c^{2}$, we find $a_{1}<2.4\times
10^{-21}{\rm\,cm^{3}s^{-1}}$ and $a_{2}<3.0\times 10^{-18}{\rm\,cm^{3}s^{-1}}$
at 95% confidence, where the product of the cross-section, $\sigma$, and
relative particle velocity, $v$, is given by $\sigma v=a_{\ell}(v/c)^{2\ell}$
and $\ell=1,2$ for $p$-, $d$-wave annihilation, respectively. Our bounds,
although failing to exclude the thermal relic cross-section for velocity-
dependent annihilation channels, are among the tightest to date.

##  I Introduction

An unambiguous detection of non-gravitational interactions of dark matter
continues to evade us. Astrophysical objects are ideal targets in the search
for novel dark matter physics, since low interaction rates can be compensated
by the large quantities of dark matter, leading to potentially detectable
signals. For over a decade, the origin of the observed excess of gamma rays at
the centre of our galaxy – known as the GCE [1, 2, 3] – has been debated, with
explanations such as unresolved point sources or dark matter annihilation
being proposed [e.g. 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. The majority of work on this
excess has focused on $s$-wave annihilation, where the product of the self-
annihilation cross-section, $\sigma$, and relative velocity, $v$, is
independent of $v$. Such models can be strongly constrained by considering the
gamma ray flux from nearby dwarf galaxies [31, 32, 33], which are able to rule
out the thermal relic cross-section for dark matter particle masses relevant
to explain the GCE.

Not all dark matter annihilation models require $\sigma v$ to be independent
of velocity. For example, for fermionic dark matter annihilation to spin-0
particles, even parity final states cannot have an $s$-wave contribution in
parity conserving theories [34, 35]. $p$-wave annihilation – where $\sigma v$
is proportional to the square of the relative velocity – dominates when
Standard Model fermion-antifermion pairs are produced from Majorana fermions
in models with minimal flavour violation due to chirality suppression of
$s$-wave annihilation [36]. If dark matter is instead a real scalar, then
$d$-wave annihilation – where the fourth power of the relative velocity is
relevant – is dominant [37, 38].

If the GCE were due to velocity-dependent dark matter annihilation, then,
despite their close proximity, dwarf spheroidals cannot provide tight
constraints on the annihilation cross-section, since the low velocity
dispersions in these objects result in heavily suppressed annihilation rates.
One should therefore test whether these velocity-dependent models are
compatible with other observations of gamma ray fluxes, since these can evade
bounds from these objects. In these scenarios, one would expect a large signal
from extra-galactic halos, where the velocity dispersion should be larger. It
has therefore been recently suggested to use extra-galactic halos to search
for velocity-dependent dark matter annihilation signals in gamma ray data
[39].

In this paper we perform this test for $p$-wave and $d$-wave annihilation,
where the cross-section is proportional to the square and fourth power of the
relative velocity between dark matter particles, respectively. As in our
previous work constraining $s$-wave annihilation [40], we utilise the CSiBORG
[41, 42, 43] suite of constrained $N$-body simulations to produce full-sky
templates for the dark matter annihilation flux, which are then compared to
the observations of the Fermi Large Area Telescope via a MCMC algorithm. The
initial conditions for these simulations are inferred using the BORG (Bayesian
Origin Reconstruction from Galaxies) algorithm and are designed to produce
three-dimensional dark matter density fields which match the observed
positions of galaxies in the 2M++ galaxy catalogue [44, 45]. We marginalise
over uncertainties in this reconstruction as well as over non-dark-matter
contributions to the gamma ray sky. Our constraints on the self-annihilation
cross-section are two (seven) orders of magnitude tighter than those obtained
using dwarf spheroidal galaxies for $p$-wave ($d$-wave) annihilation.

In Section II we describe how we forward model the dark matter annihilation
flux, discuss the models used for competing astrophysical contributions to the
signal, and outline our inference procedure. These predictions are compared to
observations, presented and discussed in Sections III and IV, and we conclude
in Section V.

##  II Methods

###  II.1 Bayesian large-scale structure inference

In this work we study the expected gamma ray flux from dark matter
annihilation in extragalactic halos. To determine the masses and locations of
these halos, we use the CSiBORG suite of 101 dark matter-only constrained
simulations [41, 42, 43]. The initial conditions of these simulations produce
present-day density fields which match the observed number densities of
galaxies in the 2M++ galaxy compilation [44, 45] and are inferred using the
BORG algorithm [see e.g. 46, 47, 48, 49, 44]. This algorithm produces a Markov
Chain of plausible initial conditions via the application of a Bayesian
forward model and marginalises over galaxy bias parameters. The initial
conditions are constrained in a box of length $677.7h^{-1}\;\text{Mpc}$ with
$256^{3}$ voxels. For each CSiBORG simulation these initial conditions are
augmented with white noise to a resolution of $2048^{3}$ within
$155h^{-1}\;\text{Mpc}$ of the Milky Way (giving a mass resolution of
$4.4\times 10^{9}\;\mathrm{M}_{\odot}$) and run to $z=0$ using RAMSES [50].
The watershed halofinder PHEW [51] is applied on-the-fly to the dark matter
particles with the standard threshold of $200\rho_{c}$. We do not consider
sub-halos and the minimum halo mass used ($M_{200\text{m}}$) is $4.4\times
10^{11}\mathrm{M}_{\odot}$. The resulting halo catalogues are publicly
available at [52].

###  II.2 Velocity-dependent $J$ factor

In this work we consider dark matter annihilation of particles of mass
$m_{\chi}$, with a cross-section $\sigma$ which depends on the relative
velocity, $v$, between annihilating dark matter particles as

$\sigma v=a_{\ell}S_{\ell}\left(\frac{v}{c}\right),$ (1)

where $a_{\ell}$ is a constant and $S_{\ell}$ can be an arbitrary function of
$v/c$. At observed energy $E_{\gamma}$, the differential photon flux,
$\mathrm{d}\Phi_{\gamma}/\mathrm{d}E_{\gamma}$, is

$\frac{\mathrm{d}\Phi_{\gamma}}{\mathrm{d}E_{\gamma}}=\frac{a_{\ell}}{8\pi
m_{\chi}^{2}}J^{\ell}\sum_{i}{\rm
Br}_{i}\left.\frac{\mathrm{d}N_{i}}{\mathrm{d}E_{\gamma}^{\prime}}\right|_{E_{\gamma}^{\prime}=E_{\gamma}\left(1+z\right)},$
(2)

where we sum over annihilation channels with branching ratios ${\rm Br}_{i}$.
This equation is valid if dark matter is entirely comprised of a single
particle whose antiparticle is itself. We introduced the $J$ factor

$\begin{split}J^{\ell}_{\mu,p}=\int_{\mu,p}\mathrm{d}\Omega\mathrm{d}\mu\,\mathrm{d}^{3}v_{1}\,\mathrm{d}^{3}v_{2}\,&S_{\ell}\left(\frac{\left|\bm{v}_{1}-\bm{v}_{2}\right|}{c}\right)\times\\\
&f\left(\bm{r}(\mu),\bm{v_{1}}\right)f\left(\bm{r}(\mu),\bm{v_{2}}\right),\end{split}$
(3)

where $f\left(\bm{r}(\mu),\bm{v}\right)$ is the distribution function of DM
particles, which we integrate over the velocities of the two particles,
$\bm{v}_{1}$ and $\bm{v}_{2}$, and along the line of sight coordinate, $\mu$.
Since we define these quantities on a discretised map, we also integrate
across the HEALPix 111http://healpix.sf.net pixel area $p$. For simplicity, we
assume that the decays occur via one channel, so the branching ratio, ${\rm
Br}_{i}$ is either 0 or 1. Note that the photon distribution function,
$\mathrm{d}N_{i}/\mathrm{d}E_{\gamma}^{\prime}$ should be evaluated at the
redshift of the halo. Since we will be using low-redshift halos in this work
($z\lesssim 0.05$) we will neglect this redshift effect and thus assume that
$E_{\gamma}^{\prime}\approx E_{\gamma}$. Introducing redshift dependence would
require inserting a redshift dependent optical depth into Eq. 2, which, in our
case, will only marginally affect the results (see for example the discussion
in Section IV.B.5 in [40]).

In this work we will consider $p$-wave and $d$-wave annihilation, i.e.
$S_{\ell}(x)=x^{2\ell}$, where $\ell=1$ ($p$-wave) or $\ell=2$ ($d$-wave). In
the following two sections we describe how one can obtain the $J$ factor for
these models from the density profile of a dark matter halo.

####  II.2.1 $J$ factor for $p$-wave annihilation

For $p$-wave annihilation, the $J$ factor can be re-written as

$\begin{split}J^{(\ell=1)}_{\mu,p}&=2\int_{\mu,p}\mathrm{d}\Omega\mathrm{d}\mu\rho^{2}\left(\bm{r}(\mu)\right)\left<\frac{\bm{v}^{2}}{c^{2}}\right>(\bm{r}(\mu))\\\
&=\frac{2}{c^{2}}\int_{\mu,p}\mathrm{d}\Omega\mathrm{d}\mu\rho^{2}\left(\bm{r}(\mu)\right)\left(\sigma_{r}^{2}(\bm{r}(\mu))+2\sigma_{\theta}^{2}(\bm{r}(\mu))\right),\end{split}$
(4)

where $\rho\left(\bm{r}(\mu)\right)$ is the dark matter density at the
position $\mu$ along the line of sight and $\sigma_{r}$ and $\sigma_{\theta}$
are the corresponding radial and tangential velocity dispersion, respectively.
We keep the explicit dependence on the line of sight position in order to
emphasize that the $J$ factor needs to be evaluated with respect to the
observer. However, since the total $J$ factor within a pixel is a scalar
quantity, one can also first evaluate it within the corresponding volume with
respect to the halo centre and then make the appropriate projection. This is
precisely how the clumpy package – which we utilise here – works (see [40, 53,
54, 55]). Now, to obtain the dispersion for a given halo, we assume spherical
symmetry, and therefore must solve the Jeans equation

$\frac{1}{\rho}\frac{\mathrm{d}}{\mathrm{d}r}\left(\rho\sigma_{r}^{2}\right)+2\beta\frac{\sigma_{r}^{2}}{r}=-\frac{\mathrm{d}\Phi}{\mathrm{d}r},$
(5)

where we introduce an anisotropy parameter

$\beta\left(r\right)=1-\frac{\sigma_{\theta}^{2}\left(r\right)}{\sigma_{r}^{2}\left(r\right)},$
(6)

and $\Phi$ is the gravitational potential of the halo.

We assume that all our halos have Navarro-Frenk-White (NFW) [56] density
profiles

$\rho\left(r\right)=\frac{\rho_{0}}{r/r_{\rm s}\left(1+r/r_{\rm
s}\right)^{2}},$ (7)

with scale length $r_{\rm s}$, defined such that the concentration of the halo
is $c=r_{\rm vir}/r_{\rm s}$ for virial radius $r_{\rm vir}$. In this work we
use the mean mass-concentration relation of Prada  _et al._ [57], Sánchez-
Conde and Prada [58] and assume zero scatter in this relation. Since the $J$
factor traces the square of the density, one may expect that subhalos or
clumps may provide a significant contribution to the $J$ factor and should
also be included. Although this is true for $s$-wave annihilation, the small
velocity dispersion of these clumps means that the overall signal is dominated
by the smooth component and we therefore ignore unresolved substructure. This
approximation has been validated by explicitly calculating the $J$ factor for
halos from hydrodynamical simulations [59], although we note that this
conclusion does not hold for the case of Sommerfeld enhanced annihilation (see
for example [60, 61]). We do not consider this case in this paper and leave it
for future work.

Defining the dimensionless radial coordinate $s\equiv r/r_{\rm vir}$, one
finds that the radial velocity dispersion for circular orbits is (see Eq. (14)
in [62])

$\begin{split}\frac{\sigma_{r}^{2}}{V_{\rm
vir}^{2}}(s)&\stackrel{{\scriptstyle\beta=0}}{{=}}\frac{c^{2}s\left(1+cs\right)^{2}}{2\left(\ln\left(1+c\right)-c/\left(1+c\right)\right)}\left[\pi^{2}+6\operatorname{Li}_{2}\left(-cs\right)\right.\\\
&\left.-\ln\left(cs\right)+\left(1+\frac{1}{c^{2}s^{2}}-\frac{4}{cs}-\frac{2}{1+cs}\right)\times\ln\left(1+cs\right)\right.\\\
&\left.+3\ln^{2}\left(1+cs\right)-\frac{1}{cs}-\frac{6}{1+cs}-\frac{1}{\left(1+cs\right)^{2}}\right],\end{split}$
(8)

where $V_{\rm vir}$ is the virial velocity, and $\operatorname{Li}_{2}$ is the
dilogarithm function. Similarly, for the extreme cases of a constant
anisotropy parameter $\beta=0.5$ and $\beta=1.0$ one finds (see Eq. (15) and
(16) from [62])

$\begin{split}\frac{\sigma_{r}^{2}}{V_{\rm
vir}^{2}}(s)\stackrel{{\scriptstyle\beta=0.5}}{{=}}&\frac{c\left(1+cs\right)^{2}}{\left(\ln\left(1+c\right)-c/\left(1+c\right)\right)}\left[-\frac{\pi^{2}}{3}-2\operatorname{Li}_{2}(-cs)\right.\\\
&+\frac{2}{1+cs}+\frac{\ln(1+cs)}{cs}+\frac{\ln(1+cs)}{1+cs}\\\
&\left.-\ln^{2}(1+cs)-\frac{1}{2(1+cs)^{2}}\right],\\\
\frac{\sigma_{r}^{2}}{V_{\rm
vir}^{2}}(s)\stackrel{{\scriptstyle\beta=1}}{{=}}&\frac{\left(1+cs\right)^{2}}{s\left(\ln\left(1+c\right)-c/\left(1+c\right)\right)}\left[\frac{\pi^{2}}{6}+\operatorname{Li}_{2}(-cs)\right.\\\
&\left.-\frac{1}{1+cs}-\frac{\ln(1+cs)}{1+cs}+\frac{\ln^{2}(1+cs)}{2}\right].\end{split}$
(9)

In Fig. 1 we show how these different velocity dispersion profiles affect the
angular dependence of the normalised $J$ factor (labelled as $\tilde{J}$) for
a halo of mass $M=1.2\times 10^{15}M_{\odot}$ (one of the most massive halos
in CSiBORG simulation 9844) at a distance of $\sim 100\,\mathrm{Mpc}$ from the
observer. The quantity $\tilde{J}$ we calculate as

$\tilde{J}^{\ell}(\theta)=\frac{J^{\ell}(\theta)}{\int\mathrm{d}\Omega
J^{\ell}(\theta)},$ (10)

where $\theta$ is the angular distance from the halo centre as seen on the
sky. We perform this calculation from close to the halo centre,
$\theta\sim{0.01}^{\circ}$, all the way up to the two times the virial radius
of the halo. Similar results for the case of the Milky Way halo (modelled as a
NFW profile) were obtained through the Eddington inversion method in [63]. We
show the results for $p$-wave with different values of the anisotropy
parameter $\beta$ (Eqs. 8 and 9) and $d$-wave calculations. The $\beta=0.5$
(dashed) and $\beta=1.0$ (dotted-dashed) cases are more strongly peaked at the
centre than the corresponding $\beta=0$ case, reflecting the behaviour of the
velocity dispersion profiles (see also Fig. 1 in [62]). Furthermore, for
$\beta=1.0$, $\tilde{J}$ diverges at the centre and hence needs to be
regularised (see Eq. 20). The $d$-wave $\tilde{J}$ shows qualitatively similar
behaviour as the corresponding $p$-wave case (thick line). For more details on
the $d$-wave calculation, see Section II.2.2 and Appendix A.

It should be noted that, for a more realistic scenario, the anisotropy
parameter should be taken as a radially varying function (see [64, 65, 66]).
We do not consider these models here, given that the three extreme scenarios
for the anisotropy parameter we pick should bracket reality. In our fiducial
analysis we choose $\beta=0$, however, and show in Section IV.2.1 that this is
a conservative choice. The corresponding $J$ factor calculation then simply
consists of substituting Eq. 8 into Eq. 4, noting that in this case
$\sigma_{\theta}=\sigma_{r}$, and integrating the final expression for all our
halos. We do this numerically using the clumpy package at a HEALPix resolution
of $\texttt{nside}=2048$, which is then subsequently degraded to the coarser
resolution ($\texttt{nside}=256$) at which we perform the inference. This
allows a more faithful calculation of the $J$ factor than performing the
calculation at $\texttt{nside}=256$, since the dependence of $J$ on $\rho$ is
non-linear.

Figure 1: The angular distribution of the normalised $J$ factor, $\tilde{J}$
(Eq. 10), as a function of angle on the sky from the halo centre, $\theta$,
for $p$-wave ($\ell=1$) and $d$-wave ($\ell=2$) annihilation channels. The
results are shown relative to the $s$-wave ($\ell=0$) result and for different
choices of anisotropy parameter, $\beta$. For this example, we used one of the
most massive halos from the CSiBORG 9844 realisation with $M=1.2\times
10^{15}M_{\odot}$ at a distance of $96\,\mathrm{Mpc}$ from the observer.

####  II.2.2 $J$ factor for $d$-wave annihilation

For $d$-wave annihilation the corresponding expression for the $J$ factor
takes the form (see also [67])

$J^{(\ell=2)}_{\mu,p}=\int_{\mu,p}\mathrm{d}\Omega\mathrm{d}\mu\rho^{2}(\bm{r})\left(2\left<\frac{\bm{v}^{4}}{c^{4}}\right>(\bm{r})+\frac{10}{3}\left<\frac{\bm{v}^{2}}{c^{2}}\right>^{2}(\bm{r})\right).$
(11)

We see that the second term appearing can be computed from a power of the
velocity dispersion (as in $p$-wave; see Eq. 4), but the first term requires
the evaluation of the fourth moment of the velocity distribution. In general,
this would require for solving the Boltzmann hierarchy up to the fourth
moment. It will not be possible to do so, however, since this system of
equations is not closed for a general self-gravitating stellar system. In
order to solve the Jeans system of differential equations, one must impose
additional constraints. This can be done by requiring the stellar system has
certain symmetries. One of the possibilities is to assume an ergodic
distribution function ($\beta=0$, see Section 4 - Box 4.3 of [68]). This then
allows to write down the following equation, relating $\left<v_{r}^{2}\right>$
and $\left<v_{r}^{4}\right>$

$\frac{\mathrm{d}}{\mathrm{d}r}\left(\rho\left<v_{r}^{4}\right>\right)\stackrel{{\scriptstyle\beta=0}}{{=}}-3\rho\left<v_{r}^{2}\right>\frac{\mathrm{d}\Phi}{\mathrm{d}r},$
(12)

where by $v_{r}$ we have denoted the projection of the velocity to the radial
direction (fixing the coordinate system to the halo centre of mass). In
addition to calculating $\left<v_{r}^{4}\right>$, we need to also calculate
the other projected moments as well as the cross-terms in order to obtain the
total fourth moment of the velocity, as is necessary for evaluating the Eq.
11. This is straightforward, due to the ergodicity of the distribution
function, as shown in Appendix A. The final expression for the fourth velocity
moment is then given by

$\left<\bm{v}^{4}\right>\stackrel{{\scriptstyle\beta=0}}{{=}}5\left<v_{r}^{4}\right>.$
(13)

This, and the solution for $\left<\bm{v}^{2}\right>$ from Eq. 8, are then
enough to evaluate the total $J$ factor for $d$-wave annihilation. See
Appendix A for more details.

###  II.3 Gamma ray data

We compare the dark matter annihilation templates to the gamma ray
observations from Fermi Large Area Telescope between mission weeks 9 and 634.
We select photons in the upper quartile of angular resolution (PS3) and event
class SOURCEVETO in the energy range $500{\rm\,MeV}-50{\rm\,GeV}$. The maximum
zenith angle is chosen to be $90^{\circ}$. These are binned spatially onto
HEALPix [69, 70] maps with resolution $\texttt{nside}=256$ and into 9
logarithmically spaced energy bins. The processing of the photon files is
performed using the Fermi
Tools222https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ and FermiPy
[71]. We also use this software to convolve our templates with the Fermi point
spread function and to obtain the exposure maps. To reduce the impact of
galactic emission and the GCE, we mask the region with galactic latitude
$|\lambda|<30^{\circ}$.

Although some of the photons detected by Fermi could have been produced by
dark matter annihilation, there are other sources which will also contribute,
and thus it is important to incorporate these in the analysis. As in our
previous work [40], we consider three contributions: point sources (psc), an
isotropic background (iso) and the Milky Way (gal). Each contribution is
modelled with a fixed spatial template for energy bin $i$,
$\\{T_{i}^{t}\left(\hat{r}\right):t\in\\{{\rm iso,\,gal,\,psc}\\}\\}$, where
we assign a different normalisation, $A_{i}^{t}$, for each energy bin and
template. We jointly infer these parameters with the amplitude of the $J$
factor template. The point source template is comprised of all extended and
point sources in the Large Area Telescope 12-year Source Catalog (4FGL-
DR3)333https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermilpsc.html and the
galactic component is modelled using the templates described in [72]. The
scalings $\\{A_{i}^{t}\\}$ are defined such that they all equal unity if the
Fermi data is perfectly described by these templates, however jointly
inferring these parameters will allow us to capture any imperfections in the
spectral modelling. We thus require an initial guess for the isotropic
component’s spectral shape, for which we use the Fermi Isotropic Spectral
Template444https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.
All templates are convolved with the Fermi point spread function using the
Fermi Tools.

###  II.4 Likelihood model

As in [40], we constrain the dark matter annihilation cross-section in two
steps. First, we do not assume a functional form for the spectrum, but infer
the normalisation of each template in each energy bin for each CSiBORG
simulation separately. We then fit this inferred spectrum to a functional form
in the second step, which allows us to constrain $a_{\ell}$.

For pixel $p$ and CSiBORG simulation $j$, the mean number of counts in energy
bin $i$ is predicted to be

$\lambda_{ipj}=\mathcal{F}_{ip}\times\left(A_{i}^{\rm
J}J_{ipj}+\sum_{t}A_{i}^{t}T_{ip}^{t}\right),$ (14)

where the $\mathcal{F}_{ip}$ is the Fermi exposure, and we have defined

$J_{ipj}\equiv\frac{J_{pj}}{J_{0}}\Delta E_{i},\quad
T_{ip}^{t}\equiv\int_{p}T_{i}^{t}\left(\hat{r}\right)\mathrm{d}\Omega,$ (15)

where $J_{0}=10^{13}{\rm\,GeV^{2}\,cm^{-5}}$ sets the units. We assume that
the photon counts are Poisson distributed, such that the likelihood for
observing $n_{ip}$ counts in pixel $p$ and energy bin $i$ given the mean
$\lambda_{ipj}$ is

$\mathcal{L}\left(n_{ip}|\lambda_{ipj}\right)=\frac{\lambda_{ipj}^{n_{ip}}\exp\left(-\lambda_{ipj}\right)}{n_{ip}!}.$
(16)

Given this likelihood, we then apply Bayes’ identity

$\mathcal{P}\left(A_{i}^{\rm J},\\{A_{i}^{t}\\},j|\mathcal{D}_{i}\right)\\\
=\frac{\mathcal{L}\left(\mathcal{D}_{i}|A_{i}^{\rm
J},\\{A_{i}^{t}\\},j\right)P\left(A_{i}^{\rm
J}\right)P\left(\\{A_{i}^{t}\\}\right)P\left(j\right)}{\mathcal{Z}\left(\mathcal{D}_{i}\right)},$
(17)

to obtain the posterior distribution for the model parameters, $\mathcal{P}$,
where $\mathcal{Z}$ is the evidence. For our priors, $P$, we enforce the
improper, uniform priors $a_{\ell},A_{i}^{\rm J}\geq 0$, and a uniform prior
on our non-dark matter templates, $A_{i}^{t}$, between 0.5 and 1.5, where
$A_{i}^{t}\in\\{A_{i}^{\text{gal}},A_{i}^{\text{iso}},A_{i}^{\text{psc}}\\}$.
We run the No U-Turn Sampler implemented in numpyro [73, 74] for 5000 steps
after an initial 1000 warmup steps. We find this typically leads to at least
$\mathcal{O}\left(2000\right)$ effective samples in all parameters, and
Gelman-Rubin statistics consistent with unity within $10^{-2}$.

We then marginalise over all other parameters and average over the
$N_{\text{sim}}=101$ CSiBORG realisations, assigning them equal weights, to
obtain the one-dimensional posterior for the annihilation flux in a given
energy bin

$\mathcal{P}\left(A_{i}^{\rm J}|\mathcal{D}_{i}\right)\\\ =\frac{1}{N_{\rm
sim}}\sum_{j}\int\mathrm{d}\\{A_{i}^{t}\\}\mathcal{P}\left(A_{i}^{\rm
J},\\{A_{i}^{t}\\},j|\mathcal{D}_{i}\right).$ (18)

Assuming a given annihilation channel and particle mass, and that each energy
bin is independent, the set of $A_{i}^{\rm J}$ can be trivially utilised to
constrain $a_{\ell}$ by evaluating the predicted $A_{i}^{\rm J}$ for the model
and multiplying the resultant probabilities. We use the pre-computed spectra
calculated by Jeltema and Profumo [75] and provided by the Fermi
collaboration555https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html
which capture the energy spectra of intermediately produced standard model
particles and their final products.

##  III Results

Figure 2: Posterior distributions of the amplitudes of the $J$ factor
templates for each energy bin, $i$, for $p$-wave (left) and $d$-wave (right)
annihilation. The black points are the maximum posterior points and the error
bars indicate the $1\sigma$ confidence intervals. For reference we also plot
the expected $A_{i}^{\rm J}$ for dark matter annihilation via the $b\bar{b}$
channel for dark matter particles of mass $m_{\chi}=100{\rm\,GeV}/c^{2}$, with
$a_{1}=10^{-19}{\rm\,cm^{3}s^{-1}}$ ($p$-wave) and
$a_{2}=10^{-16}{\rm\,cm^{3}s^{-1}}$ ($d$-wave). All fluxes are consistent with
zero, indicating no detection of dark matter $p$-wave or $d$-wave
annihilation. Figure 3: 95% upper limit on the dark matter $p$-wave (upper,
$\ell=1$) and $d$-wave (lower, $\ell=2$) annihilation cross-section parameter,
$a_{\ell}$, where the annihilation cross-section is $\sigma
v=a_{\ell}(v/c)^{2\ell}$. The solid lines and points are from this work. Each
constraint is derived at a fixed dark matter particle mass $m_{\chi}$ as
indicated on the $x$-axis. The red dotted line is a constraint obtained from
dwarf spheroidal (dSph) satellites of the Milky Way [76, 77] for the
$b\bar{b}$ channel. We find that our constraints are approximately two orders
of magnitude tighter than those from dSphs due to the larger velocity
dispersions in massive extragalactic halos. The dSph constraint for $d$-wave
is off the top of the plot. The brown dashed line is the constraint for the
$e^{+}e^{-}$ channel, directly detecting the resulting electrons and
positrons, assuming they are produced by the galactic halo [78]. This is
slightly tighter than our $e^{+}e^{-}$ constraint. Figure 4: Constraints on
$p$-wave (upper, $\ell=1$) and $d$-wave (lower, $\ell=2$) annihilation cross-
section, $\sigma v=a_{\ell}\left(v/c\right)^{2\ell}$, from halos in a given
mass range from a single CSiBORG simulation (9844). Each mass bin has a width
$\Delta\log_{10}\left(M_{\rm h}/\mathrm{M}_{\odot}\right)=1$, and we use a
sliding window. We plot the 95% upper limits on $a_{\ell}$ for the $b\bar{b}$
and $e^{+}e^{-}$ channels; the dashed horizontal lines are the constraints
obtained from all halos together. We see that our constraints are driven by
the most massive halos and that this is more important for $d$-wave
annihilation due to the stronger dependence of annihilation rate on the
relative velocity between dark matter particles.

The first step of our inference infers the amplitudes of the $J$ factor
templates, $A_{i}^{\rm J}$, for each energy bin, $i$. In Fig. 2 we plot the
one-dimensional posterior distributions of these parameters and see that these
coefficients are consistent with zero at $1\sigma$ confidence for all energy
bins for both $p$\- and $d$-wave annihilation. In the second step we fit this
inferred spectrum to a dark matter annihilation model and hence can constrain
$a_{\ell}$. For reference, in Fig. 2 we plot example spectra for annihilation
via the $b\bar{b}$ channel for $m_{\chi}=100{\rm\,GeV}/c^{2}$ at a fiducial
$a_{\ell}$.

In Fig. 3 we convert these posteriors to 95% upper limits on the $p$\- and
$d$-wave self-annihilation cross-section ($\sigma v=a_{\ell}(v/c)^{2\ell}$) as
a function of dark matter particle mass, $m_{\chi}$, for the production of any
standard-model quark, charged lepton or gauge boson (except photons). We do
not attempt to constrain $m_{\chi}$, and instead find constraints on
$a_{\ell}$ conditioned on this variable for the range
$m_{\chi}\in[2,500]{\rm\,GeV}/c^{2}$. In all cases we find no evidence for
$p$-wave or $d$-wave dark matter annihilation, as all the fluxes in Fig. 2 are
consistent with zero. For dark matter annihilating to $b\bar{b}$, if
$m_{\chi}=10{\rm\,GeV}/c^{2}$, we constrain $a_{1}<2.4\times
10^{-21}{\rm\,cm^{3}\,s^{-1}}$ and $a_{2}<3.0\times
10^{-18}{\rm\,cm^{3}\,s^{-1}}$ at 95% confidence. For an annihilation which
forms lighter quarks, we obtain practically identical constraints, whereas
lepton production yields constraints which are approximately an order of
magnitude less stringent. For large $m_{\chi}\gtrsim 100{\rm\,GeV}/c^{2}$, our
upper limits for the production of top quarks or $W$ or $Z$ bosons are similar
to those for lighter quarks, but these processes cannot occur at lower
$m_{\chi}$ so we do not obtain constraints for smaller dark matter particle
masses.

Note that, in Fig. 3, the constraints weaken as $m_{\chi}$ approaches $m_{A}$,
if the annihilation product is $A\bar{A}$. Near this value, the process is
strongly suppressed and hence the constraints no longer lie within the range
of the plot, so we do not plot these points.

For comparison, for $p$-wave, in Fig. 3 we also plot the constraints one
obtains for the $b\bar{b}$ channel from dwarf spheroidals. As anticipated from
the much smaller velocity dispersion in dwarf galaxies, our constraints are
significantly tighter, by approximately two orders of magnitude for $p$-wave
annihilation, and seven orders of magnitude for $d$-wave. This demonstrates
the increased importance of large halo masses for velocity-dependent
annihilation, whereas for $s$-wave annihilation the close proximity of these
objects more than compensates for their small masses, leading to tighter
constraints [see e.g. Figure 7 of 40]. We perform a more detailed comparison
to the literature in Section IV.4.

##  IV Discussion

###  IV.1 Understanding the results

To determine which halos are responsible for our constraints, we produce $J$
factor maps for halos only within a given mass range and run the inference
using only those halos. For computational convenience, we perform this on one
representative CSiBORG simulation only (9844). The results of this inference
are plotted in Fig. 4, where we use moving mass bins of width
$\Delta\log_{10}\left(M_{\rm h}/\mathrm{M}_{\odot}\right)=1$ and consider dark
matter particles of mass $m_{\chi}=8{\rm\,GeV}/c^{2}$ annihilating via the
$b\bar{b}$ or $e^{+}e^{-}$ channels. We see that our constraints are dominated
by the most massive objects in our simulations, i.e. halos of mass $M_{\rm
h}\sim 10^{14-16}\,\mathrm{M}_{\odot}$. This is due to both the larger
densities and velocity dispersions in massive objects, as was predicted by
Baxter  _et al._ [39] (see their Fig. 3). If we continued Fig. 4 to larger
$M_{\rm h}$ we would see that the constraints would begin to weaken, due to
the lack of very massive objects in our simulation. We observe that the
constraining power as a function of mass for $p$-wave annihilation plateaus at
an earlier halo mass than for $d$-wave annihilation due to the weaker
dependence on the relative velocity between dark matter particles, which
increases with mass. Qualitatively similar results are seen for other channels
and dark matter particle masses. Again, this is consistent with the
expectations of Baxter  _et al._ [39]. We find that the total $J$ factors we
obtained, for the 128 most massive halos of one of our CSiBORG realisation
(9844), agree to within factor of $\sim 5$ when compared to the expression in
Eq. (3.3) of [39]. This is a consequence of more detailed modelling of the
$J$-factor in this work; for example, if Baxter  _et al._ [39] chose a
different mass definition, given the strong dependence of the $p$\- and
$d$-wave annihilation channels on total mass of the object, this difference
could be alleviated.

To assess which energy bins drive our constraints, for each channel and
$m_{\chi}$ we calculate the log-likelihood at the 95% upper limit on
$a_{\ell}$ for each energy bin and compare to the values for $a_{\ell}=0$. For
both $p$\- and $d$-wave annihilation, we find that the first three energy bins
are the most constraining for $m_{\chi}\lesssim 30{\rm\,GeV}/c^{2}$, whereas
the fourth energy bin becomes particularly important for higher masses. This
is the same behaviour as for $s$-wave annihilation [40], which is attributed
to the lower expected photon counts in the higher energy bins, and therefore
weaker constraining power. As one increases $m_{\chi}$, one expects more
photons at higher energy, and thus the relative importance of these bins
increases.

###  IV.2 Systematic uncertainties

####  IV.2.1 $J$ factor calculation

In Equation 18 we marginalised over uncertainties in the large-scale density
modes which are constrained from the BORG algorithm, as well as smaller-scale
unconstrained modes which were added as white noise in the initial conditions.
In our fiducial analysis we use the full suite of 101 CSiBORG simulations,
however to verify this is a sufficiently large number, we rerun our analysis
one hundred times using fifty randomly selected CSiBORG simulations (with
repeats) each time. This yields an uncertainty on the 95% upper limit of
$a_{\ell}$ of no more than 5.9% and 6.5% for $p$-wave and $d$-wave
annihilation, respectively, when considering all channels and particle masses
separately. The median change in the constraint across these runs at a given
mass and channel changes by no more than 1.2% compared to the fiducial
analysis. Thus, we conclude that this is an appropriate number of simulations.

For our $p$-wave analysis, we assumed circular orbits for dark matter
particles ($\beta=0$), leading to the radial velocity dispersion given in Eq.
8. Similar expressions exist for alternative values of $\beta$, thus to
determine the impact of this assumption on our constraints, we compute $J$
factor maps assuming $\beta=0.5$ or $\beta=1$ (although the latter case is
unrealistic as this corresponds to purely radial orbits) for a representative
CSiBORG simulation (9844) and rerun our analysis to obtain alternative
constraints on $a_{1}$.

For $\beta=1$, one finds that the $J$ factor diverges logarithmically towards
the centre, and thus simply computing Eq. 4 would predict an infinite
annihilation flux from each halo. In practice, if the annihilation rate
becomes sufficiently high, then it will disrupt the density profile of the
halo as the dark matter particles will be depleted towards the centre,
presumably forming a cored profile. Thus, Eq. 4 is only valid up to some core
radius, $r_{\rm c}$. We estimate this to be the radius at which the dynamical
time scale of the halo ($1/\sqrt{G\rho}$) equals the annihilation time scale
($\rho\sigma v/m_{\chi}$). Since we expect this to occur close to the centre
of the halo, we use the leading order term for $\sigma_{r}^{2}$,

$\frac{\sigma_{r}^{2}}{V_{\rm
vir}^{2}}\stackrel{{\scriptstyle\beta=1}}{{=}}\frac{\pi^{2}-9}{6\left(\log\left(1+c\right)-\frac{c}{1+c}\right)}\frac{1}{s}+\mathcal{O}\left(s^{0}\right),$
(19)

and assume $\rho\sim\rho_{0}r_{\rm s}/r$ (see Eq. 7). This leads to a core
radius of

$r_{\rm c}\sim\left(\frac{\pi^{2}-9}{6}\frac{a_{1}}{m_{\chi}}\frac{V_{\rm
vir}^{2}}{\log\left(1+c\right)-\frac{c}{1+c}}\sqrt{\frac{\rho_{0}}{G}}\right)^{\frac{2}{3}}r_{\rm
vir}.$ (20)

We have verified that, for the core radii that this implies, the
$\mathcal{O}(s^{0})$ term is negligible compared to the leading order term.

For simplicity, we compute Eq. 4 using $r_{\rm c}$ as the lower limit of the
integral and add the contribution from the cored region separately.
Furthermore, we note that Eq. 20 depends on $a_{1}$ and $m_{\chi}$, and thus a
fully-consistent analysis should recompute the $J$ factor at each step in the
MCMC and for each channel. However, this is computationally infeasible.
Instead, we set $a_{1}=2.4\times 10^{-21}\,{\rm cm^{3}\,s^{-1}}$ and
$m_{\chi}=100{\rm\,GeV}/c^{2}$ in this formula; this is a characteristic WIMP
mass with $a_{1}$ at approximately the 20th percentile for the $b\bar{b}$
channel when $\beta=0.5$. This typically leads to a core radius of $r_{\rm
c}\sim 0.1-1{\rm\,pc}$. We find that the constraints on $a_{1}$ for $\beta=1$
are much tighter than this, and thus our approximation yields a conservative
constraint.

For both $\beta=0.5$ and $\beta=1$, we find that the anisotropic $J$ factors
lead to tighter constraints for all masses and channels, by up to a factor of
$2.3$ for $\beta=0.5$ and up to two orders of magnitude for $\beta=1$. This
simply reflects the general behaviour of the $J$ factor for these two
different anisotropic profiles. As can be seen from Fig. 1, the $\beta=0.5$
and $\beta=1.0$ profiles are more strongly peaked towards the centre than the
corresponding $\beta=0$ case, with $\beta=1.0$ diverging more quickly towards
the centre and giving total $J$ factors between $50-100$ times larger than the
corresponding $\beta=0$ case. From this and Eq. 2, it is then easy to see that
the derived constraints for these respective cases will show similar trends.
Choosing $\beta=0$ is therefore the most conservative of the cases considered,
justifying our use of it in our fiducial analysis.

For all halos considered in this work, we have not attempted to add
corrections to the halo density profiles due to baryonic effects, and have
instead used profiles calibrated against dark-matter-only simulations. In the
presence of baryons, adiabatic contraction during galaxy formation [79, 80]
can steepen the central density profile, or it can be made shallower due to
the subsequent stellar feedback [81, 82]. The different resulting density
profiles would change our $J$ factor calculations and hence our constraints on
$a_{\ell}$. Indeed, cosmological zoom hydrodynamical simulations have
demonstrated that baryons do affect the $J$ factors for velocity-dependent
annihilation, typically increasing $J$ from a factor of a few up to factors of
several hundred [83, 84]. Assuming this increase in $J$ can be applied to our
halos, then our constraints on $a_{\ell}$ are conservative, however we leave
it to future work to implement more precise corrections for baryonic physics.

####  IV.2.2 Constrained simulation volume

As already argued in Section VI.B.5 of [40], the choice of HEALPix nside=256
for our inference pipeline allows for circumventing the issue of the limited
volume we have for our CSiBORG simulation suite ($z\lesssim 0.05$). This is
because any sources beyond the CSiBORG volume would most likely be unresolved
at this resolution and hence captured by the isotropic template.

Our discovery that the most massive halos in the local universe dominate the
$p$\- and $d$-wave constraints (running the inference using the $\sim 20$
largest halos gives constraints which differ by a factor only $\mathcal{O}(1)$
compared to the full result) raises the question of whether stronger
constraints could be obtained using a dedicated cluster catalogue extending to
higher redshift than CSiBORG . To address this we calculate the $J$ factor
produced by the Master Catalogue of X-ray Clusters (MCXC; [85]), a compilation
of clusters detected through their X-ray-emitting gas.

We find that the total $J$ factor of the CSiBORG catalogue is larger by factor
of $\sim 10$ than that from MCXC. This is a result of the incompleteness of
the MCXC catalogue as well as the various selection effects involved, which
outweigh its larger redshift range and are absent from CSiBORG . The CSiBORG
suite is therefore optimal in terms of its all-sky constraining power.

Nevertheless, since the $p$\- and $d$-wave annihilation channels are strongly
dependent on relative velocity of dark matter particles, and thus the total
halo mass, it might still be worthwhile to explore constrained simulations
with larger simulation volumes than CSiBORG to include the contributions from
more distant massive halos. Further constraining power may be extracted from
the angular dependence of the $s$-, $p$\- and $d$-wave annihilation signals,
which we do not fully exploit due to the relatively coarse resolution of our
all-sky maps. Useful future work would therefore be to use higher-resolution
simulations or observations of clusters to compare their predicted and
observed flux profiles in greater detail. Note that the angular resolution of
the Fermi data is approximately $\texttt{nside}\approx 1024$, so significant
improvements in resolution are possible. The positions and masses of local
dark matter clusters may be found in the public CSiBORG halo catalogues [52].

####  IV.2.3 Cluster masses

As demonstrated in Fig. 4 and indicated by Baxter  _et al._ [39], for $p$-wave
and $d$-wave annihilation, the constraints on cross-section are dominated by
the largest mass objects. It is therefore more important to correctly obtain
cluster masses than it is when considering $s$-wave models [cf. Fig. 9 of 40].
In this work we have used the masses obtained from constrained $N$-body
simulations, with initial conditions inferred using the BORG algorithm. This
is an alternative to traditional estimates of cluster masses, which typically
use one of the virial theorem [86], X-ray measurements [87], the Sunyaev
Zel’dovich effect [88, 89] or weak lensing [90, 91]. These measurements tend
to lead to significantly different measurements of masses for the same
clusters [e.g. Fig. 3 of 92], so it is not clear whether these traditional
methods offer a more appropriate route to obtaining $J$ factor maps as each
method would produce very different constraints.

It has been found that the halo mass function in CSiBORG is higher than the
average $\Lambda$CDM expectation at the massive end [43]. Although cosmic
variance is high in this regime, this still represents a $\sim 2\sigma$
effect. To determine the effect of this potential mass discrepancy, we obtain
a single new posterior sample of the initial conditions with a different
gravity model (using the COmoving Lagrangian Acceleration (COLA) method [93]),
different time-stepping in the forward-model, a different method for
generating power spectra (CLASS [94] instead of Eisenstein & Hu [95]) and a
more robust likelihood [96], designed to remove the effects of unknown
systematic uncertainties on scales larger than $1.4^{\circ}$. We then run a
simulation with these initial conditions but otherwise identical to CSiBORG .
This reconstruction produces 2- and 3-point statistics of the density field
and a halo mass function which is closer to that expected from a random
$\Lambda$CDM realisation. However, the re-simulated halo masses are
systematically smaller than those observed and some objects (e.g. the Perseus
cluster) are less prominent than in our fiducial reconstruction. These
modifications to the older BORG algorithm are studied in more detail in [97].

Due to the variety of changes, this simulation is useful to test the
robustness of our results to the reconstruction analysis. We produce $J$
factor maps for this simulation and rerun the end-to-end inference. We find
this new simulation produces weaker constraints than our fiducial analysis,
with changes in the constraints of up to a factor of $\sim 5$ and $\sim 12$
for $p$\- and $d$-wave, respectively. This is to be expected given the
systematically smaller cluster masses obtained with the updated initial
conditions and is much greater than the variation in constraints between
CSiBORG realisations, which is typically $\sim 40-50\%$. We note that this
variation is comparable to what one would obtain using measured clustered mass
because different measurements give highly variable results [92]. It is the
most significant uncertainty in our analysis.

####  IV.2.4 Astrophysical templates

In our fiducial analysis, we use the latest galactic diffuse model provided by
the Fermi Collaboration (gll_iem_v07). Although we find little degeneracy
between the normalisation of the astrophysical templates and the amplitude of
the $J$-factor maps (and hence $a_{\ell}$ constraints), we nonetheless re-run
our inference using an older template (gll_iem_v02) to estimate the systematic
effect of the choice of galactic template. As in Bartlett  _et al._ [40], we
find that the constraints usually weaken slightly when using the older model,
by up to 24% for $p$-wave and up to 31% for $d$-wave, although at low
$m_{\chi}$ the $d$-wave constraints can be up to 7% tighter with the older
model. This variation is small compared to that caused by the cluster masses.

###  IV.3 Thermal relics

In this work we have focused on dark matter particles of mass
$2-500{\rm\,GeV}/c^{2}$ since, if the thermally-average annihilation cross-
section is $\left<\sigma v\right>_{\rm th}\approx 3\times
10^{-26}{\rm\,cm^{3}s^{-1}}$ [98] and these particles are thermal relics, then
the current abundance of dark matter is reproduced. This cross-section is
comparable to the electroweak coupling of the Standard Model, motivating the
study of these dark matter candidates. In fact, the previously quoted
$\left<\sigma v\right>_{\rm th}$ is appropriate for $s$-wave annihilation
only. Hence, in this section we extend the analysis of Steigman  _et al._ [98]
to obtain the equivalent values for $p$\- and $d$-wave annihilation, and
compare to the $a_{\ell}$ we obtain.

We assume that $\sigma v=a_{\ell}v^{2\ell}$ and that dark matter particles are
non-relativistic and follow a Boltzmann distribution. The thermally-averaged
cross-section at temperature $T$ is then

$\left<\sigma
v\right>=\frac{\int\mathrm{d}^{3}k_{1}\mathrm{d}^{3}k_{2}\exp\left(-\frac{E_{1}+E_{2}}{T}\right)\sigma
v}{\int\mathrm{d}^{3}k_{1}\mathrm{d}^{3}k_{2}\exp\left(-\frac{E_{1}+E_{2}}{T}\right)}=\frac{4^{\ell}\Gamma\left(\ell+\frac{3}{2}\right)}{\Gamma\left(\frac{3}{2}\right)}\frac{a_{\ell}}{x^{\ell}},$
(21)

where $k_{1}$ and $k_{2}$ are the phase-space momenta and $x=m_{\chi}/T$. One
must then solve the equation governing the number density, $n$, of dark matter

$\frac{\mathrm{d}Y}{\mathrm{d}x}=\frac{s\left<\sigma
v\right>}{Hx}\left[1+\frac{1}{3}\frac{\mathrm{d}\left(\log g_{\rm
s}\right)}{\mathrm{d}\left(\log T\right)}\right]\left(Y_{\rm
eq}^{2}-Y^{2}\right),$ (22)

where $Y=n/s$ for entropy density $s$, $g_{\rm s}$ is the number of
relativistic degrees of freedom contributing to the entropy density, $H$ is
the Hubble parameter, and the equilibrium value of $Y$ is

$Y_{\rm
eq}=\frac{45}{2\pi^{4}}\left(\frac{\pi}{8}\right)^{1/2}\frac{g_{\chi}}{g_{\rm
s}}x^{3/2}\exp\left(-x\right),$ (23)

where the dark matter particle has $g_{\chi}$ degrees of freedom. Since this
is a numerically stiff differential equation, we follow Steigman  _et al._
[98] and solve for $W=\log Y$ instead. We also take $g_{\rm s}(T)$ from the
calculations of Laine and Schröder [99]. Using an implicit Runge-Kutta method
of the Radau IIA family of order 5 [100], we solve from $x=1$ to $x=10^{3}$,
with an initial condition set using the approximate solution [101]

$Y\left(x\right)\approx\frac{x^{\ell+2}}{2\lambda}+Y_{\rm eq}\left(x\right),$
(24)

which is valid long before freeze-out, where $\lambda\equiv 2.76\times
10^{35}m_{\chi}\left<\sigma v\right>_{0}$ for $\left<\sigma
v\right>=\left<\sigma v\right>_{0}x^{-\ell}$ and where $m_{\chi}$ is in ${\rm
GeV}$ and $\left<\sigma v\right>_{0}$ is in ${\rm cm^{3}\,s^{-1}}$. We then
convert our solution to a density parameter

$\Omega_{\chi}=\frac{8\pi G}{3H_{0}^{2}}m_{\chi}s_{0}Y_{0},$ (25)

where $H_{0}$, $s_{0}$ and $Y_{0}$ are the values of $H$, $s$ and $Y$
evaluated today. By comparing this value to $\Omega_{\chi}h^{2}=0.11$ [102],
for $H_{0}=100h{\rm\,km\,s^{-1}\,Mpc^{-1}}$, we repeat the analysis for
varying $a_{\ell}$ at fixed $m_{\chi}$ to obtain the thermal relic cross-
section which reproduces the present-day dark matter abundance.

This is plotted in Fig. 5 as a function of $m_{\chi}$ for $\ell=0,1,2$. We
find that the required $p$-wave coefficient is $a_{1}\sim(1-3)\times
10^{-25}{\rm\,cm^{3}\,s^{-1}}$, whereas for $d$-wave this is
$a_{2}\sim(4-9)\times 10^{-25}{\rm\,cm^{3}\,s^{-1}}$. In both cases, this is
several orders of magnitude smaller than the constraints we find in Fig. 3,
and thus we cannot rule out dark matter being a thermal relic with a velocity-
dependent cross-section using large scale structure. These conclusions are
consistent with those discussed in the context of the extragalactic gamma-ray
background [103, 104]; a detection of a $p$-wave annihilation signal would
suggest a non-thermal origin of the present-day dark matter abundance.

We note that approximate analytic solutions of Eq. 25 exist for $p$\- and
$d$-wave annihilation [e.g. 105, 37]. In these cases, if freeze-out occurs
when $x=x_{\rm f}$ and $x_{\rm f}$ is the same for all $\ell$, then to obtain
the $p$\- and $d$-wave $\left<\sigma v\right>_{\rm th}$ from the $s$-wave
analysis, one would multiply $a_{0,{\rm th}}$ by $x_{\rm f}/3\approx 6.7$ and
$x_{\rm f}^{2}/20\approx 20$, respectively, where we have used $x_{\rm
f}\approx 20$. These ratios give the dashed horizontal lines in Fig. 5, when
multiplied by the fiducial value of $a_{0}=3\times
10^{-26}{\rm\,cm^{3}s^{-1}}$.

Figure 5: Annihilation cross-section coefficient ($\sigma
v=a_{\ell}\left(v/c\right)^{\ell}$) as a function of dark matter particle
mass, $m_{\chi}$, required to give a present-day dark matter abundance of
$\Omega_{\chi}h^{2}=0.11$ for a thermal relic. We plot the results for $s$-,
$p$\- and $d$-wave annihilation and see that the required $a_{\ell}$ increases
by less than an order of magnitude for a unit change of $\ell$. The dashed
horizontal lines give the standard value of $a_{0}=3\times
10^{-26}{\rm\,cm^{3}s^{-1}}$ for $s$-wave, and this value multiplied by
$x_{\rm f}/3$ and $x_{\rm f}^{2}/20$ for $p$\- and $d$-wave, assuming a
freeze-out value of $x$ of $x_{\rm f}=20$.

###  IV.4 Comparison to the literature

An attempt to combine the relic abundance requirements with the GCE and lack
of detection in dwarf galaxies was proposed by Choquette  _et al._ [106]. They
suggested that a metastable particle undergoing $s$-wave annihilation can
produce the thermal abundance of dark matter particles, and these daughter
particles can then undergo $p$-wave annihilation in the present epoch. They
find that a dark matter particle mass $m_{\chi}\sim 90{\rm\,GeV}/c^{2}$ with
an annihilation cross-section given by $a_{1}\sim 10^{-20}{\rm\,cm^{3}s^{-1}}$
for the $b\bar{b}$ channel can reproduce the GCE and survive all other tests.
For this mass and channel, the $a_{1}$ required to explain the GCE is ruled
out at over $1\sigma$ confidence in our analysis, suggesting that large scale
structure disfavours this interpretation. In particular, they convert
constraints on $a_{0}$ from clusters to constraints on $a_{1}$ by multiplying
by multiples of the square of the velocity dispersion, and they find the
Fornax cluster to be the most constraining object, with constraints which are
$\sim 2$ orders of magnitude weaker than those given in Fig. 3. Since we
perform a more careful analysis and obtain tight constraints on $a_{1}$,
future work should be dedicated to repeating the analyses cited by Choquette
_et al._ [106] but with the correct velocity-dependent $J$ factors.

If dark matter annihilated via the $e^{+}e^{-}$ channel, then the resulting
cosmic rays from the Milky Way could be detectable in Voyager 1 and the AMS-02
data. Boudaud  _et al._ [78] considered this effect, and found constraints on
$p$-wave annihilation for this channel which are tighter than our constraints,
by up to approximately one order of magnitude (see Fig. 3). Their results are
found to be robust to choice of profile used for the galactic halo, although
are only applicable to the $e^{+}e^{-}$ channel.

This is supported by Johnson  _et al._ [107], who obtain a null detection of
$p$-wave annihilation as the cause of the GCE, since they do not find a spike
in emission around Sgr A*. Moreover, since $s$\- and $p$-wave annihilation
produce different spatial distributions of emission [63] for the same density
profile, it is argued that, to explain the GCE, $p$-wave annihilation would
require a slope of the density profile of the Milky Way much steeper than
normally considered ($\gamma\sim 1.5-1.7$ compared to $\gamma\sim 1.2-1.4$ for
$s$-wave) [108], further questioning whether velocity-dependent annihilation
is the cause of this detection.

Alongside the GCE, one of the most common targets for indirect dark matter
annihilation searches are dwarf galaxies in the Local Group due to their
proximity and low baryon fraction. In Fig. 3 we compare our $p$-wave
constraints to those obtained by these objects [76, 109, 77] for the
$b\bar{b}$ channel, and find that our constraints are up to two orders of
magnitude tighter. Similarly, we find that our $d$-wave constraints are
approximately seven orders of magnitude tighter than those previously found
from studying 25 dwarf spheroidal galaxies of the Milky Way [77], and hence
lie outside the range of Fig. 3. This is in contrast to the results of
$s$-wave annihilation, where dwarf galaxies give a constraint which is an
order of magnitude tighter than that from large scale structure [40]. Assuming
that clusters and dwarf galaxies have velocity dispersions of
$\mathcal{O}(100){\rm\,km\,s^{-1}}$ and $\mathcal{O}(1-10){\rm\,km\,s^{-1}}$,
respectively, then one could estimate the $p$\- and $d$-wave constraints by
multiplying the $s$-wave ones by $\sigma^{2\ell}$. Given this ratio of
velocity dispersions, one would anticipate that clusters can produce $p$\- and
$d$-wave constraints which are factors of $10^{2-3}$ and $10^{5-6}$ tighter,
approximately as we find. Despite being less constraining for velocity-
dependent annihilation, if one did detect such a signal, then dwarf galaxies
could still have an important role, since each object would have the same
spectrum, but the different density and velocity profiles would result in
different spatial distributions of the signal, which one could use to
determine $\ell$ [110].

If dark matter annihilated, then the energy injected into the intergalactic
medium would modify the Universe’s post-recombination ionisation history, and
hence the observed cosmic microwave background (CMB). The large velocity
dispersion in halos means that this is a predominantly low-redshift effect and
could be detectable in the lowest multipoles of the CMB. Combining results
from HST with BAO measurements, a joint analysis of data from Planck, WMAP9,
ACT and SPT found no evidence of $p$-wave dark matter annihilation, ruling out
cross-sections corresponding to $a_{1}\sim
10^{-16}-10^{-12}{\rm\,cm^{3}\,s^{-1}}$ for the production of leptons at
masses $m_{\chi}\sim 1-10^{3}{\rm\,GeV}/c^{2}$ at 95% confidence [111]. Our
constraints are approximately four to five orders of magnitude more stringent
than these results. One can similarly constrain dark matter annihilation via
this process using the Lyman-$\alpha$ forest, but these constraints are also
$\sim 4$ orders of magnitude weaker than those found in this work [112].

##  V Conclusions

Annihilation of dark matter particles with velocity-dependent cross-sections
can result in fluxes of gamma rays which are higher from extragalactic halos
than from dwarf galaxies in the Local Group. In this work, we use the CSiBORG
suite of constrained $N$-body simulations of the nearby Universe (within
$155h^{-1}\;\text{Mpc}$) to forward-model the expected emission from all halos
down to $M_{\text{vir}}=4.4\times 10^{11}M_{\odot}$ and therefore produce
full-sky, field-level, energy-dependent intensity predictions for $p$-wave and
$d$-wave dark matter annihilation. We combine these with templates describing
emission by baryonic processes in the Milky Way, an isotropic background and
point sources, and compare to observations from the Fermi Large Area
Telescope. Marginalising over the amplitudes of these templates and
uncertainties in the reconstructed dark matter density field, we find no
evidence for velocity-dependent dark matter annihilation to any standard model
particle.

We place upper limits on the annihilation cross-section as a function of dark
matter particles mass, $m_{\chi}$, in the range
$m_{\chi}=2-500{\rm\,GeV}/c^{2}$ and for different annihilation channels. For
the $b\bar{b}$ channel, defining $\sigma v=a_{\ell}\left(v/c\right)^{2\ell}$
where $\ell=1,2$ for $p$ and $d$-wave, respectively, we constrain
$a_{1}<2.4\times 10^{-21}{\rm\,cm^{3}\,s^{-1}}$ and $a_{2}<3.0\times
10^{-18}{\rm\,cm^{3}\,s^{-1}}$ at 95% confidence for
$m_{\chi}=10{\rm\,GeV}/c^{2}$. Our constraints are weaker for the $e^{+}e^{-}$
channel, with $a_{1}<1.5\times 10^{-20}{\rm\,cm^{3}\,s^{-1}}$ and
$a_{2}<1.9\times 10^{-17}{\rm\,cm^{3}\,s^{-1}}$ at 95% confidence for this
particle mass. We find that these constraints are dominated by the most
massive halos ($M_{\rm h}\sim 10^{14-16}\,\mathrm{M}_{\odot}$). Our
constraints are four to five orders of magnitude more stringent than those on
$p$-wave annihilation from cosmic microwave background analyses, and two
(seven) orders of magnitude tighter than $p$-wave ($d$-wave) constraints from
dwarf galaxies, as one would expect from the larger velocity dispersion in
higher mass objects. Our bounds are however still several orders of magnitude
above the value required to rule out dark matter being a thermal relic with a
velocity-dependent annihilation cross-section.

In this work we have demonstrated the constraining power of the full halo
field on velocity-dependent dark matter annihilation through a full-sky field-
level inference, and obtained constraints which are several orders of
magnitude tighter than alternative probes. Future work should be dedicated to
higher-resolution inferences using the most constraining (highest mass)
objects, incorporating other observations to constrain and marginalise over
the density and velocity profiles, as well as investigating the case of
Sommerfeld-enhanced annihilation channels as proposed in [60, 61].

The dominant systematic uncertainty in our analysis is the density field
reconstruction (in particular its ability to reproduce high halo masses
accurately), and hence future work should focus on understanding and improving
this modelling (see [97]). Since anisotropic motions and baryonic processes
can increase the $J$ factor and hence the predicted signal, more precise
modelling of these effects may achieve tighter constraints even with the same
data.

## Data availability

This work required modifications to the clumpy package in order to model
velocity-dependent dark matter annihilation. These modifications will be
included in a future clumpy release. Our $J$ and $D$ factor maps, which may be
useful for other indirect detection analyses of current or future data, will
be made available upon publication. As supplementary material on arXiv we
provide the Mathematica notebooks used in deriving Appendix A.

###### Acknowledgements.

We thank Céline Combet, Sten Delos, Moritz Hütten, Jens Jasche, Eiichiro
Komatsu, Guilhem Lavaux, David Maurin, Fabian Schmidt and Ewoud Wempe for
useful discussions. We thank Jonathan Patterson for smoothly running the
Glamdring Cluster hosted by the University of Oxford, where most of the data
processing was performed. AK acknowledges support from the Starting Grant
(ERC-2015-STG 678652) “GrInflaGal” of the European Research Council at MPA.
DJB is supported by the Simons Collaboration on “Learning the Universe” and
was supported by STFC and Oriel College, Oxford. HD is supported by a Royal
Society University Research Fellowship (grant no. 211046). This project has
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (grant agreement No
693024). This work was done within the Aquila Consortium (https://www.aquila-
consortium.org/). Some of the results in this paper have been derived using
the healpy and HEALPix [69, 70], clumpy [53, 54, 55], numpyro [73, 74], numpy
[113], scipy [114] and FermiPy [71] packages. For the purpose of open access,
the authors have applied a Creative Commons Attribution (CC BY) licence to any
Author Accepted Manuscript version arising.

## References

  * Goodenough and Hooper [2009] L. Goodenough and D. Hooper, arXiv e-prints , arXiv:0910.2998 (2009).
  * Ajello  _et al._ [2016] M. Ajello _et al._ , ApJ 819, 44 (2016).
  * Linden  _et al._ [2016] T. Linden, N. L. Rodd, B. R. Safdi, and T. R. Slatyer, Phys. Rev. D 94, 103013 (2016).
  * Hooper and Goodenough [2011] D. Hooper and L. Goodenough, Physics Letters B 697, 412 (2011).
  * Hooper and Linden [2011] D. Hooper and T. Linden, Phys. Rev. D 84, 123005 (2011).
  * Hooper and Slatyer [2013] D. Hooper and T. R. Slatyer, Physics of the Dark Universe 2, 118 (2013).
  * Zhou  _et al._ [2015] B. Zhou _et al._ , Phys. Rev. D 91, 123010 (2015).
  * Daylan  _et al._ [2016] T. Daylan _et al._ , Physics of the Dark Universe 12, 1 (2016).
  * Cholis  _et al._ [2021] I. Cholis, Y.-M. Zhong, S. D. McDermott, and J. P. Surdutovich, arXiv e-prints , arXiv:2112.09706 (2021).
  * Grand and White [2022] R. J. J. Grand and S. D. M. White, MNRAS 511, L55 (2022).
  * Abazajian [2011] K. N. Abazajian, J. Cosmology Astropart. Phys 2011, 010 (2011).
  * O’Leary  _et al._ [2015] R. M. O’Leary, M. D. Kistler, M. Kerr, and J. Dexter, arXiv e-prints , arXiv:1504.02477 (2015).
  * Petrović  _et al._ [2015] J. Petrović, P. D. Serpico, and G. Zaharijas, J. Cosmology Astropart. Phys 2015, 023 (2015).
  * Lee  _et al._ [2016] S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer, and W. Xue, Phys. Rev. Lett. 116, 051103 (2016).
  * Buschmann  _et al._ [2020] M. Buschmann, N. L. Rodd, B. R. Safdi, L. J. Chang, S. Mishra-Sharma, M. Lisanti, and O. Macias, Phys. Rev. D 102, 023023 (2020).
  * Gautam  _et al._ [2021] A. Gautam _et al._ , arXiv e-prints , arXiv:2106.00222 (2021).
  * Hooper  _et al._ [2013] D. Hooper, I. Cholis, T. Linden, J. M. Siegal-Gaskins, and T. R. Slatyer, Phys. Rev. D 88, 083009 (2013).
  * Cholis  _et al._ [2015a] I. Cholis, D. Hooper, and T. Linden, J. Cosmology Astropart. Phys 2015, 043 (2015a).
  * Leane and Slatyer [2019] R. K. Leane and T. R. Slatyer, arXiv e-prints , arXiv:1904.08430 (2019).
  * Abazajian and Kaplinghat [2012] K. N. Abazajian and M. Kaplinghat, Phys. Rev. D 86, 083511 (2012).
  * Gordon and Macías [2013] C. Gordon and O. Macías, Phys. Rev. D 88, 083521 (2013).
  * Abazajian  _et al._ [2014] K. N. Abazajian, N. Canac, S. Horiuchi, and M. Kaplinghat, Phys. Rev. D 90, 023526 (2014).
  * Calore  _et al._ [2015] F. Calore, I. Cholis, and C. Weniger, J. Cosmology Astropart. Phys 2015, 038 (2015).
  * Horiuchi  _et al._ [2016] S. Horiuchi, M. Kaplinghat, and A. Kwa, J. Cosmology Astropart. Phys 2016, 053 (2016).
  * Petrović  _et al._ [2014] J. Petrović, P. D. Serpico, and G. Zaharijaš, J. Cosmology Astropart. Phys 2014, 052 (2014).
  * Carlson and Profumo [2014] E. Carlson and S. Profumo, Phys. Rev. D 90, 023015 (2014).
  * Cholis  _et al._ [2015b] I. Cholis, C. Evoli, F. Calore, T. Linden, C. Weniger, and D. Hooper, J. Cosmology Astropart. Phys 2015, 005 (2015b).
  * Gaggero  _et al._ [2015] D. Gaggero, M. Taoso, A. Urbano, M. Valli, and P. Ullio, J. Cosmology Astropart. Phys 2015, 056 (2015).
  * Macias  _et al._ [2018] O. Macias, C. Gordon, R. M. Crocker, B. Coleman, D. Paterson, S. Horiuchi, and M. Pohl, Nature Astronomy 2, 387 (2018).
  * Bartels  _et al._ [2018] R. Bartels, E. Storm, C. Weniger, and F. Calore, Nature Astronomy 2, 819 (2018).
  * Ackermann  _et al._ [2015] M. Ackermann _et al._ , Phys. Rev. Lett. 115, 231301 (2015).
  * Albert  _et al._ [2017] A. Albert _et al._ , ApJ 834, 110 (2017).
  * Gammaldi  _et al._ [2021] V. Gammaldi _et al._ , Phys. Rev. D 104, 083026 (2021).
  * Kim and Lee [2007] Y. G. Kim and K. Y. Lee, Phys. Rev. D 75, 115012 (2007).
  * Lee  _et al._ [2008] K. Y. Lee, Y. G. Kim, and S. Shin, Journal of High Energy Physics 2008, 100 (2008).
  * Kumar and Marfatia [2013] J. Kumar and D. Marfatia, Phys. Rev. D 88, 014035 (2013).
  * Giacchino  _et al._ [2013] F. Giacchino, L. Lopez-Honorez, and M. H. G. Tytgat, J. Cosmology Astropart. Phys 2013, 025 (2013).
  * Toma [2013] T. Toma, Phys. Rev. Lett. 111, 091301 (2013).
  * Baxter  _et al._ [2022] E. J. Baxter, J. Kumar, A. D. Paul, and J. Runburg, J. Cosmology Astropart. Phys 2022, 026 (2022).
  * Bartlett  _et al._ [2022] D. J. Bartlett, A. Kostić, H. Desmond, J. Jasche, and G. Lavaux, Phys. Rev. D 106, 103526 (2022).
  * Bartlett  _et al._ [2021] D. J. Bartlett, H. Desmond, and P. G. Ferreira, Phys. Rev. D 103, 023523 (2021).
  * Desmond  _et al._ [2022] H. Desmond, M. L. Hutt, J. Devriendt, and A. Slyz, MNRAS 511, L45 (2022).
  * Hutt  _et al._ [2022] M. L. Hutt, H. Desmond, J. Devriendt, and A. Slyz, MNRAS 516, 3592 (2022).
  * Lavaux and Jasche [2016] G. Lavaux and J. Jasche, MNRAS 455, 3169 (2016).
  * Jasche and Lavaux [2019] J. Jasche and G. Lavaux, Astron. Astrophys. 625, A64 (2019).
  * Jasche and Wandelt [2012] J. Jasche and B. D. Wandelt, MNRAS 425, 1042 (2012).
  * Jasche and Wandelt [2013] J. Jasche and B. D. Wandelt, MNRAS 432, 894 (2013).
  * Jasche  _et al._ [2010] J. Jasche, F. S. Kitaura, B. D. Wandelt, and T. A. Enßlin, MNRAS 406, 60 (2010).
  * Jasche  _et al._ [2015] J. Jasche, F. Leclercq, and B. D. Wandelt, J. Cosmology Astropart. Phys 2015, 036 (2015).
  * Teyssier [2002] R. Teyssier, A&A 385, 337 (2002).
  * Bleuler  _et al._ [2015] A. Bleuler, R. Teyssier, S. Carassou, and D. Martizzi, Computational Astrophysics and Cosmology 2, 5 (2015).
  * Hutt  _et al._ [2022] M. L. Hutt, H. Desmond, J. Devriendt, and A. Slyz, “The effect of local universe constraints on halo abundance and clustering,”  (2022).
  * Charbonnier  _et al._ [2012] A. Charbonnier, C. Combet, and D. Maurin, Computer Physics Communications 183, 656 (2012).
  * Bonnivard  _et al._ [2016] V. Bonnivard, M. Hütten, E. Nezri, A. Charbonnier, C. Combet, and D. Maurin, Computer Physics Communications 200, 336 (2016).
  * Hütten  _et al._ [2019] M. Hütten, C. Combet, and D. Maurin, Computer Physics Communications 235, 336 (2019).
  * Navarro  _et al._ [1997] J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 490, 493 (1997).
  * Prada  _et al._ [2012] F. Prada, A. A. Klypin, A. J. Cuesta, J. E. Betancort-Rijo, and J. Primack, MNRAS 423, 3018 (2012), arXiv:1104.5130 [astro-ph.CO] .
  * Sánchez-Conde and Prada [2014] M. A. Sánchez-Conde and F. Prada, MNRAS 442, 2271 (2014).
  * Piccirillo  _et al._ [2022] E. Piccirillo _et al._ , J. Cosmology Astropart. Phys 2022, 058 (2022).
  * Facchinetti  _et al._ [2022] G. Facchinetti, M. Stref, T. Lacroix, J. Lavalle, J. Pérez-Romero, D. Maurin, and M. A. Sánchez-Conde, arXiv e-prints , arXiv:2203.16491 (2022), arXiv:2203.16491 [astro-ph.CO] .
  * Lacroix  _et al._ [2022] T. Lacroix _et al._ , J. Cosmology Astropart. Phys 2022, 021 (2022).
  * Łokas and Mamon [2001] E. L. Łokas and G. A. Mamon, MNRAS 321, 155 (2001).
  * Boddy  _et al._ [2018] K. K. Boddy, J. Kumar, and L. E. Strigari, Phys. Rev. D 98, 063012 (2018).
  * Osipkov [1979] L. Osipkov, Pisma v Astronomicheskii Zhurnal 5, 77 (1979).
  * Merritt [1985] D. Merritt, The astronomical journal 90, 1027 (1985).
  * Gerhard [1991] O. E. Gerhard, Monthly Notices of the Royal Astronomical Society 250, 812 (1991).
  * Boucher  _et al._ [2022] B. Boucher, J. Kumar, V. B. Le, and J. Runburg, Physical Review D 106, 023025 (2022).
  * Binney and Tremaine [2011] J. Binney and S. Tremaine,  _Galactic dynamics_ , Vol. 13 (Princeton university press, 2011).
  * Zonca  _et al._ [2019] A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset, E. Hivon, and K. Gorski, Journal of Open Source Software 4, 1298 (2019).
  * Górski  _et al._ [2005] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, ApJ 622, 759 (2005).
  * Wood  _et al._ [2017] M. Wood, R. Caputo, E. Charles, M. Di Mauro, J. Magill, J. S. Perkins, and Fermi-LAT Collaboration, in  _35th International Cosmic Ray Conference (ICRC2017)_ , International Cosmic Ray Conference, Vol. 301 (2017) p. 824.
  * Acero  _et al._ [2016] F. Acero _et al._ , ApJS 223, 26 (2016).
  * Phan  _et al._ [2019] D. Phan, N. Pradhan, and M. Jankowiak, arXiv preprint arXiv:1912.11554 (2019).
  * Bingham  _et al._ [2019] E. Bingham _et al._ , J. Mach. Learn. Res. 20, 28:1 (2019).
  * Jeltema and Profumo [2008] T. E. Jeltema and S. Profumo, J. Cosmology Astropart. Phys 2008, 003 (2008).
  * Zhao  _et al._ [2016] Y. Zhao, X.-J. Bi, H.-Y. Jia, P.-F. Yin, and F.-R. Zhu, Phys. Rev. D 93, 083513 (2016).
  * Boddy  _et al._ [2020] K. K. Boddy, J. Kumar, A. B. Pace, J. Runburg, and L. E. Strigari, Phys. Rev. D 102, 023029 (2020).
  * Boudaud  _et al._ [2019] M. Boudaud, T. Lacroix, M. Stref, and J. Lavalle, Phys. Rev. D 99, 061302 (2019), arXiv:1810.01680 [astro-ph.HE] .
  * Blumenthal  _et al._ [1986] G. R. Blumenthal, S. M. Faber, R. Flores, and J. R. Primack, ApJ 301, 27 (1986).
  * Gnedin  _et al._ [2004] O. Y. Gnedin, A. V. Kravtsov, A. A. Klypin, and D. Nagai, ApJ 616, 16 (2004).
  * Pontzen and Governato [2012] A. Pontzen and F. Governato, Mon. Not. Roy. Astron. Soc. 421, 3464 (2012).
  * Del Popolo and Pace [2016] A. Del Popolo and F. Pace, Astrophys. Space Sci. 361, 162 (2016), [Erratum: Astrophys.Space Sci. 361, 225 (2016)], 1502.01947 .
  * Board  _et al._ [2021] E. Board _et al._ , J. Cosmology Astropart. Phys 2021, 070 (2021).
  * McKeown  _et al._ [2022] D. McKeown _et al._ , MNRAS 513, 55 (2022).
  * Piffaretti  _et al._ [2011] R. Piffaretti, M. Arnaud, G. Pratt, E. Pointecouteau, and J.-B. Melin, Astronomy & Astrophysics 534, A109 (2011).
  * Merritt [1987] D. Merritt, ApJ 313, 121 (1987).
  * Evrard  _et al._ [1996] A. E. Evrard, C. A. Metzler, and J. F. Navarro, ApJ 469, 494 (1996).
  * Sunyaev and Zeldovich [1970] R. A. Sunyaev and Y. B. Zeldovich, Comments on Astrophysics and Space Physics 2, 66 (1970).
  * Sunyaev and Zeldovich [1980] R. A. Sunyaev and I. B. Zeldovich, ARA&A 18, 537 (1980).
  * Bonnet  _et al._ [1994] H. Bonnet, Y. Mellier, and B. Fort, ApJ 427, L83 (1994).
  * Fahlman  _et al._ [1994] G. Fahlman, N. Kaiser, G. Squires, and D. Woods, ApJ 437, 56 (1994).
  * Stopyra  _et al._ [2021] S. Stopyra, H. V. Peiris, A. Pontzen, J. Jasche, and P. Natarajan, MNRAS 507, 5425 (2021).
  * Tassev  _et al._ [2013] S. Tassev, M. Zaldarriaga, and D. J. Eisenstein, J. Cosmology Astropart. Phys 2013, 036 (2013).
  * Blas  _et al._ [2011] D. Blas, J. Lesgourgues, and T. Tram, Journal of Cosmology and Astro-Particle Physics 2011, 034 (2011).
  * Eisenstein and Hu [1998] D. J. Eisenstein and W. Hu, ApJ 496, 605 (1998).
  * Porqueres  _et al._ [2019] N. Porqueres, D. Kodi Ramanah, J. Jasche, and G. Lavaux, A&A 624, A115 (2019), arXiv:1812.05113 [astro-ph.CO] .
  * Stopyra  _et al._ [2023] S. Stopyra, H. V. Peiris, A. Pontzen, J. Jasche, and G. Lavaux, “Towards accurate field-level inference of massive cosmic structures,” (2023), arXiv:2304.09193 [astro-ph.CO] .
  * Steigman  _et al._ [2012] G. Steigman, B. Dasgupta, and J. F. Beacom, Phys. Rev. D 86, 023506 (2012).
  * Laine and Schröder [2006] M. Laine and Y. Schröder, Phys. Rev. D 73, 085009 (2006).
  * Hairer and Wanner [1996] E. Hairer and G. Wanner,  _Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems_, Vol. 14 (1996).
  * Kolb and Turner [1990] E. W. Kolb and M. S. Turner,  _The Early Universe_, Vol. 69 (1990).
  * Komatsu  _et al._ [2011] E. Komatsu _et al._ , ApJS 192, 18 (2011).
  * Campbell  _et al._ [2010] S. Campbell, B. Dutta, and E. Komatsu, Phys. Rev. D 82, 095007 (2010).
  * Campbell and Dutta [2011] S. Campbell and B. Dutta, Phys. Rev. D 84, 075004 (2011).
  * Griest and Seckel [1991] K. Griest and D. Seckel, Phys. Rev. D 43, 3191 (1991).
  * Choquette  _et al._ [2016] J. Choquette, J. M. Cline, and J. M. Cornell, Phys. Rev. D 94, 015018 (2016).
  * Johnson  _et al._ [2019] C. Johnson, R. Caputo, C. Karwin, S. Murgia, S. Ritz, J. Shelton, and Fermi-LAT Collaboration, Phys. Rev. D 99, 103007 (2019).
  * Kiriu  _et al._ [2022] K. Kiriu, J. Kumar, and J. Runburg, arXiv e-prints , arXiv:2208.14002 (2022).
  * Zhao  _et al._ [2018] Y. Zhao, X.-J. Bi, P.-F. Yin, and X. Zhang, Phys. Rev. D 97, 063013 (2018).
  * Baxter  _et al._ [2021] E. J. Baxter, J. Kumar, A. B. Pace, and J. Runburg, J. Cosmology Astropart. Phys 2021, 030 (2021).
  * Diamanti  _et al._ [2014] R. Diamanti, L. Lopez-Honorez, O. Mena, S. Palomares-Ruiz, and A. C. Vincent, J. Cosmology Astropart. Phys 2014, 017 (2014).
  * Liu  _et al._ [2021] H. Liu, W. Qin, G. W. Ridgway, and T. R. Slatyer, Phys. Rev. D 104, 043514 (2021), arXiv:2008.01084 [astro-ph.CO] .
  * Harris  _et al._ [2020] C. R. Harris _et al._ , Nature 585, 357 (2020).
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , Nature Methods 17, 261 (2020).
  * Dejonghe [1986] H. Dejonghe, Phys. Rep. 133, 217 (1986).

##  Appendix A Calculation of $J$ factor for $d$-wave

To calculate the $d$-wave $J$ factor at a given distance from the halo centre,
one needs to evaluate the corresponding fourth velocity moment at that
distance, $\langle\bm{v}^{4}\rangle$. To do so, we assume ergodicity of the
halo dark matter distribution function, which directly leads to the relation
Eq. 12. However, this equation solves for the fourth moment of the radial
velocity component, while we need the total $\langle\bm{v}^{4}\rangle$. One
can show that in the case of an ergodic distribution function the following
holds (see Problem 4.29 in [68] and Section 1.3 from [115])

$\langle v_{\theta}^{n-j}v_{\phi}^{j-k}v_{r}^{k}\rangle=\langle
v_{r}^{n}\rangle\frac{\Gamma(\tfrac{k+1}{2})\Gamma(\tfrac{n-j+1}{2})\Gamma(\tfrac{j-k+1}{2})}{\pi\Gamma(\tfrac{n+1}{2})}.$
(26)

From this, the following relations directly follow

$\displaystyle\langle v_{\theta}^{4}\rangle=\langle v_{\phi}^{4}\rangle$
$\displaystyle=\langle v_{r}^{4}\rangle$ (27) $\displaystyle\langle
v_{\theta}^{2}v_{\phi}^{2}\rangle,=\langle v_{\theta}^{2}v_{r}^{2}\rangle$
$\displaystyle=\langle v_{\phi}^{2}v_{r}^{2}\rangle=\frac{2}{3}\langle
v_{r}^{4}\rangle,$ (28)

as does Eq. 13. Now, the only remaining issue is finding a solution to Eq. 12,
since the second part of Eq. 11 can be directly obtained from Eq. 8. Note,
that since we are assuming ergodicity, $\beta=0$.

We assume NFW profiles (Eq. 7) for our halos, which we rewrite as

$\rho(s)=\frac{\rho_{0}}{cs(1+cs)^{2}},$ (29)

where $c$ is the concentration, and $s=r/r_{\mathrm{vir}}$. We can then solve
Eq. 12 using Eq. 8. Noting that the gravitational potential in this notation
is

$\Phi\left(s\right)=-4\pi G\rho_{0}r_{\rm
s}^{2}\frac{\log\left(1+cs\right)}{cs},$ (30)

one obtains

$\displaystyle\rho\langle v_{r}^{4}\rangle=$
$\displaystyle\frac{4}{3}\pi^{2}G^{2}\rho_{0}^{3}r_{s}^{4}\begin{aligned}
&\bigg{[}-\frac{6(cs(cs(8cs+21)+6)-1)\log^{2}(cs+1)}{c^{3}s^{3}}+\frac{6\text{Li}_{2}\left(c^{2}s^{2}\right)(10cs+9\log(cs+1))}{cs}\\\
&-\frac{6\log(cs+1)}{c^{2}s^{2}(cs+1)}\bigg{(}cs\left(-3\pi^{2}(cs+1)+cs(3cs+35)+25\right)+3cs(cs+1)\log(cs)+2\bigg{)}\\\
&+288\text{Li}_{3}(cs+1)+\text{Li}_{2}(cs)\left(-\frac{108\log(cs+1)}{cs}-120\right)+72\text{Li}_{2}(cs+1)(3-4\log(cs+1))\\\
&-\frac{52(\pi cs+\pi)^{2}+174cs+177}{(cs+1)^{2}}+\frac{6}{cs(cs+1)^{2}}\\\
&+\frac{6(8cs+9)\log^{3}(cs+1)}{cs}+72\log(cs)(3-2\log(cs+1))\log(cs+1)+18\log(cs)\bigg{]}\end{aligned}$
$\displaystyle+\mathrm{const},$ (31)

where $\operatorname{Li}_{2}$ and $\operatorname{Li}_{3}$ represent the
dilogarithm and trilogarithm, respectively. We use the implementations
provided in the polylogarithm library
666https://github.com/Expander/polylogarithm to evaluate these. We impose
$\lim_{s\to\infty}(\rho\langle v_{r}^{4}\rangle)(s)=0$, which fixes the
integration constant to $0$. This is a consequence of the fact that the
velocity dispersion itself goes to zero far away from the halo, for the case
of ergodic phase-space distribution function which we assumed throughout for
our halos. This in turn implies that the velocity distribution has a vanishing
second moment and hence, by definition, the fourth moment is zero too in this
limit.

This analytical solution was implemented within the clumpy package and tested
against a numerical solution. For all our results regarding $d$-wave, we
always utilise the analytical solution from Appendix A and use the numerical
solution only for cross-checking the results. We demonstrate the obtained $J$
factor for $d$-wave in Fig. 1 obtained by substituting the expression from
Appendix A into Eq. 11 and normalising by the total integrated $J$ factor
value and dividing by the corresponding quantity for $s$-wave (labelled as
$\tilde{J}^{\ell=2}/\tilde{J}^{\ell=0}$). One sees that the shape of the
angular distribution of the normalised $J$ factor is qualitatively very
similar to the $\beta=0$, $p$-wave case and shows the same behaviour towards
the halo centre and the edge of the halo.

Note that similar results have been obtained in Boddy  _et al._ [63], who used
Eddington inversion formula for obtaining the velocity moments. We chose to
directly integrate the Jeans equations instead, since we were also interested
in studying the anisotropic cases for the $p$-wave annihilation channel, and
hence for consistency used the same approach for the $d$-wave $J$ factor
calculation.