# Got plenty of nothing: cosmic voids as a probe of particle dark matter

S. Arcari  E. Pinetti  N. Fornengo

###### Abstract

The search for a particle dark matter signal in terms of radiation produced by
dark matter annihilation or decay has to cope with the extreme faintness of
the predicted signal and the presence of masking astrophysical backgrounds. It
has been shown that using the correlated information between the dark matter
distribution in the Universe with the fluctuations of the cosmic radiation
fields has the potential to allow setting apart a pure dark matter signal from
astrophysical emissions, since spatial fluctuations in the radiation field due
to astrophysical sources and dark matter emission have different features. The
cross-correlation technique has been proposed and adopted for dark matter
studies by looking at dark matter halos (overdensities). In this paper we
extend the technique by focusing on the information on dark matter
distribution offered by cosmic voids, and by looking specifically at the
gamma-ray dark matter emission: we show that, while being underdense and
therefore producing a reduced emission as compared to halos, nevertheless in
voids the relative size of the cross-correlation signal due to decaying dark
matter vs. astrophysical sources is significantly more favourable, producing
signal-to-background ratios $S/B$ (even significantly) larger than 1 for decay
lifetimes up to $2\times 10^{30}$ s. This is at variance with the case of
halos, where $S/B$ is typically (even much) smaller than 1. We show that
forthcoming galaxy surveys such as Euclid combined with future generation
gamma-ray detectors with improved specifications have the ability to provide a
hint of such a signal with a predicted significance up to $5.7\sigma$ for
galaxies and $1.6\sigma$ for the cosmic shear. The bound on the dark matter
lifetime attainable exploiting voids is predicted to improve on current bounds
in a mass range for the WIMP of $25\div 900$ GeV.

##  1 Introduction

The evidence for the presence of dark matter in the Universe, while being
overwhelming, still lacks of an understanding of the nature of what we call
dark matter. A natural explanation relies on the existence of one (or more)
new type(s) of elementary particle(s), which would form the dark matter. A
test of the particle physics interpretation of dark matter in terms of a new
elementary particle is expected to produce a variety of signals which are
possible due to its particle physics nature. Among them, the production of
cosmic radiation in terms of radiation (from radio to gamma rays, depending on
the mass of the particle, which sets the maximal energy of the produced
signal) or neutrinos is one of the most investigated channels.

The search for a particle dark matter signal in terms of radiation produced by
dark matter annihilation or decay, while offering a wide set of opportunities
since the signal can be produced in every structure where dark matter is
present (galaxies, clusters of galaxies, filamentary structures connecting
them), nevertheless it is typically faint for most of the particle physics
models, which makes it often dominated by masking astrophysical backgrounds.
In order to attempt to extract the dark matter signal from the backgrounds, in
[1, 2] it was proposed to look for the cross-correlation between a
gravitational tracer of dark matter (like the cosmic shear or the galaxy
distribution) and cosmic radiation fields (like gamma-rays, for heavy dark
matter like WIMPs - see also [3]). The correlated information between where
dark matter is and the ensuing fluctuations induced on the radiation fields
could give a handle to separate the two signals, since spatial fluctuations in
the radiation field due to astrophysical sources (which are essentially
pointlike) and dark matter emission (more diffuse) have different features.

This has been further elaborated and explored on data for the cross-
correlation between dark matter gamma-ray emission with galaxies [4, 5, 6, 7,
8, 9, 10, 11, 12], clusters [13, 14, 15, 16], CMB lensing [17, 18] and cosmic
shear [19, 20, 21, 22, 23], discussed in different energy bands like X-rays
[24, 25] and the NIRB [26, 27, 28] or extended to different gravitational
tracers like the HI intensity mapping [29].

In the previous literature, the cross-correlation technique has been proposed
and adopted for dark matter studies by looking at dark matter halos
(overdensities). In this paper we investigate the possibility to extend the
cross-correlation technique by using the information on dark matter
distribution offered by cosmic voids, for which catalogs start to become
available [30, 31, 32]. While the signal that originates in cosmic voids is
expected to be weaker than the one produced in dark matter halos, due to the
fact that voids are underdense and therefore the electromagnetic emission
intensity produced by dark matter annihilation or decay is smaller in size
than the one produced in dark matter halos, nevertheless we will show in the
following that the relative size of the dark matter signal vs. the cross-
correlation signal due to astrophysical sources like active galactic nuclei
(AGN) or star-forming galaxies (SFG), which represent the relevant background
for our observable of interest, is significantly more favourable in voids than
in halos. Therefore, by selecting cosmic voids as large-scale structure (LSS)
tracers can potentially offer a cleaner signal as compared to the signal due
to overdensities: this makes the cross-correlation signal in voids an
interesting counter-part to the signal in halos, with a trade off between
weaker but cleaner vs. stronger but with higher background signal.

The paper is organized as follows: in Sec. 2 we summarize the formalism to
derive the cross-correlation angular power spectrum of any given pair of
intensity fields, and the associated variance, extended to a halo and void
model of LSS [33], which we briefly review in Sec. 3, along with the
prescriptions for the halo and void statistics and density distribution. Sec.
4 serves as a repertory of the cross-correlation 3D power spectra needed for
each case study, whilst Sec. 5 reports the corresponding window functions. In
Sec. 6 we show our results and finally conclusions are drawn in Sec. 7.
Appendix A and B discuss the gamma-ray luminosity function (GLF) of unresolved
astrophysical sources and their mass to luminosity function, respectively.
Appendix C shows all contributions to the total 3D power spectrum of each
source field considered in this work and its variation upon the choice of the
void profile. Appendix D outlines the correlation between the free parameters
of our analysis.

Throughout the paper we assume a flat $\Lambda$CDM cosmology with cosmological
parameters as derived by the _Planck_ satellite in 2018 [34]:
$H_{0}=67.66\;\text{km s}^{-1}\text{ Mpc}^{-1}$, $\Omega_{m}=0.3111$,
$\Omega_{\Lambda}=0.6889$.

##  2 The cross-correlation signal for particle dark matter

The dark matter signal we consider is the statistical cross-correlation
between the unresolved gamma-ray emission from particle dark matter
annihilation or decay and the distribution of mass in the Universe traced
through cosmic shear or through the distribution of galaxies in the Universe.
In this Section we briefly review the basic elements of the formalism to
cross-correlate two density fields by means of the angular power spectrum
(APS), while in the next Sections we will specify in details the relevant
ingredients that determine the APS for the density fields of interest. We
closely follow the formalism introduced in Refs. [1, 3], although we extend it
to include cosmic voids.

The source intensity of an observable $i$, taken along a given direction
$\vec{n}$, can be written as:

$I_{i}(\vec{n})=\int{\textrm{d}}\chi\;g_{i}(\chi,\;\vec{n})\tilde{W}(\chi)\;,$
(2.1)

where $\chi$ is the comoving radial distance, $g_{i}(\chi,\;\vec{n})$
represents the density field of the source $i$ and $\tilde{W}(\chi)$ is a
window function characterising the average intensity field as a function of
distance (or equivalently of redshift). It is convenient to define a
normalized window function $W(\chi)=\langle
g_{i}(\chi,\;\vec{n})\rangle\tilde{W}(\chi)$ such that $\langle
I_{i}\rangle=\int{\textrm{d}}\chi\;W(\chi)$. By expanding the intensity
fluctuations of two given fields $i$ and $j$ in spherical harmonics, we can
compute the cross-correlation angular power spectrum (CAPS):

$C_{\ell}^{ij}=\int\frac{{\textrm{d}}\chi}{\chi^{2}}\;W_{i}(\chi)W_{j}(\chi)P_{ij}\left(k=\frac{\ell}{\chi},\;\chi\right)\;,$
(2.2)

where the window function $W_{i}(\chi)$ describes how the observable $i$ is
distributed in redshift and its shape strongly depends on the physics behind
the chosen signal. The 3D power spectrum is defined through $\langle
f_{i}(\chi,\;\bm{k})f_{j}(\chi,\;\bm{k})\rangle=(2\pi)^{3}\delta_{D}(\bm{k}-\bm{k^{\prime}})P_{ij}(k,\;\chi)$
where $f_{i}=g_{i}-\langle g_{i}\rangle$ is the fluctuation of the density
field. In Eq. (2.2) we adopt the Limber approximation [35, 36, 37], which is
typically valid for the relevant scales explored in cross-correlation studies
involving gamma rays [3, 2, 29].

In the halo model [38], the 3D power spectrum can be split into two terms,
taking into account the correlations between source fields either from the
same halo (1Halo) or from two different halos (2Halo). In this work we extend
the CAPS formalism in a Halo-Void Model of LSS [33], allowing matter to lie
within halos and voids. The model can be further extended to include a dust
component [33], but we leave this further level of complexity to future works.
Introducing more structures in the model leads to new correlation terms,
therefore the power spectrum will be decomposed not only into 1Halo (1H) and
2Halo (2H) terms, but it will also include 1Void (1V) and 2Void (2V)
components as well as the Halo-Void (HV) mixed term. A detailed derivation of
the relevant terms will be given in the next sections.

Under the hypothesis of gaussianity, the variance of the predicted CAPS is
(see e.g. [29]):

$(\Delta
C^{ij}_{\ell})^{2}=\frac{1}{(2\ell+1)f_{\text{sky}}}\left[(C^{ij}_{\ell})^{2}+\left(C^{ii}_{\ell}+\frac{N^{i}}{(B_{\ell}^{i})^{2}}\right)\bigg{(}C^{jj}_{\ell}+\frac{N^{j}}{(B_{\ell}^{j})^{2}}\bigg{)}\right]\;,$
(2.3)

where $f_{\text{sky}}$ is the observed fraction of the sky, $C^{ii}_{\ell}$
and $C^{jj}_{\ell}$ represent the auto-correlation angular power spectra
associated to the observable $i$ and $j$, respectively, whereas $N^{i}$ and
$N^{j}$ are their corresponding noises. Their beam functions $B_{\ell}^{i}$
and $B_{\ell}^{j}$ in harmonic space refer to the angular resolution of the
chosen detector. Eq. (2.3) represents the uncertainty on the predicted CAPS,
and will be used to determine whether the signal is detectable.

All quantities in Eqs. (2.2) and (2.3) will be explicitly specified below. In
the following, $i$ labels the gamma-ray intensity, either from astrophysical
sources or annihilating/decaying particle dark matter, while $j$ represents
the gravitational tracer under consideration in our analysis, that is cosmic
shear and galaxy distribution.

##  3 The model for halos and voids

An overwhelming amount of theoretical and observational evidences favors the
idea that structure in the Universe has risen out of a nearly homogeneous
primordial cosmos through gravitational instability. The process of structure
formation is hierarchical (larger clustering are formed through the continuous
merging of smaller structures), as backed by most existing theories, and
therefore strongly depends on initial conditions, whose knowledge is one of
the primary questions in cosmology.

Tracing back the evolution of structures retains the potential to broaden our
understanding on the primordial Universe and much effort has been put into
developing both analytical and numerical schemes to explain such phenomena. On
one hand, linear and higher order perturbation theory descriptions of
gravitational clustering [39, 40, 41, 42] from Gaussian initial conditions
explain the evolution and mildly non-linear clustering of dark matter, but
break in highly non-linear regimes [43] (i.e. at scales smaller than few
megaparsecs) and do not provide a rigorous framework to describe the
clustering of galaxies. Also, very large scales are troublesome to work with,
from the observational side, because of the small amount of data available. On
the other hand, smaller non-linear scales can only be described by numerical
dark matter simulations of the LSS clustering [44, 45]. The latter show that
an initially smooth matter distribution evolves into a complex web of knots,
sheets and filaments. These numerical simulations provide detailed information
on the distribution of mass within these structures [46, 47] when performed at
high resolution but relatively small volume, but they are also useful to
constrain the abundance and spatial distribution of structures in the Universe
when performed at lower resolution with large volume [48, 49]. The drawback is
that simulations are usually computationally expensive and unable to provide
an analytical description of the initial conditions, on which they are highly
dependent.

Data from forthcoming large-area imaging and redshift surveys of galaxies [50]
and weak lensing [51, 52] will provide constraints on the dark matter
distribution on large scales as well as on the galaxy formation history. At
the same time, the Sunyaev–Zel’dovich effect [53] probes the distribution of
the pressure on large scales, and can be observed through wide-field surveys
[54]. Moreover, dark and/or baryonic matter leave their imprints on CMB in the
form of secondary temperature fluctuations on small scales [53, 55, 56].

The halo model (HM) [38] has been, up to today, the uttermost successful
analytical description of non-linear scales and provides a self-consistent
explanation of the observations discussed above. This model relies on the
central assumption that all matter in the Universe lies within dark matter
halos and offers a simple framework to explain the transition between non-
linear and linear scales, which are dominated by the 1Halo and 2Halo terms,
respectively. However, the predicted power spectrum is in accordance with
N-body simulations only within 20% around the transition and, at the same
time, it either requires a normalisation of the 2Halo term or for the halo
abundance to be integrated down to very low and untested masses (M $\ll
10^{4}\;h^{-1}\;\text{M}_{\odot}$), in order to account for all matter in the
Universe. As discussed in [33], there have been several attempts to modify the
HM [57, 58, 59, 60]. These alternative models typically either introduce new
free parameters that cannot be fitted using only halo properties, or do not
significantly improve the HM predictions.

Recently, a successful self-consistent modification of the HM has been
proposed by Voivodic et al. [33], namely the halo-void model (HVM). This model
is based on relaxing the central assumption of the HM and allowing matter to
lie not only within halos but also within cosmic voids.

Including voids as building blocks of the Universe leads to new terms in the
prediction of the cross-correlation signal: the 1Void, 2Void, and Halo-Void
terms, in addition to the 1Halo and 2Halo contributions (already present in
the HM). In the past decades, the scientific community has mostly focused on
modeling the dark matter halos, while cosmic voids have been largely
unappreciated. Nevertheless, voids constitute the dominant volume fraction of
the Universe and can be used as powerful, independent probes for our theories
of structure formation. Underdense regions are indeed ideal environments to
constraint dark energy [61, 62, 63, 64] and modified gravity [65, 66, 67, 68].
In addition, thanks to their extreme sensitivity to background cosmological
changes, voids are more closely related to initial conditions [69, 70, 71] and
can be valuable in order to constrain cosmology [72, 73, 74, 75, 76, 77, 78].
Also, they are the perfect case-study for the excursion set theory111The
excursion set approach provides a useful framework to describe the formation
histories of gravitationally bound structures such as virialized halos or
cosmic voids. [79] given their rather spherical symmetry222Unlike the
evolution of density peaks, primordial asphericity of negative density
perturbations is quickly lost as they expand [80]. and the existence of a
refined repertory of fitting prescriptions for their density profile, mass
function and linear bias.

The HVM requires the following ingredients: the halo and void density profile,
the halo and void mass function and the halo and void linear bias. The model
computes the halo and void statistics through the excursion set formalism [79]
with two barriers [68, 81, 80, 82]. This approach guarantees a fully self-
consistent model that takes into account the void-in-cloud, void-in-void and
cloud-in-void effects, in addition to the cloud-in-cloud effect that appears
in the HM. Moreover, there is no need for a normalization on large scales or
for the abundances to be integrated down to very low halo masses, as the
matter within smaller halos is taken into account in larger voids. Voivodic et
al. [33] also show how considering both halos and voids improves the
transition between the 2Halo and the 1Halo term.

Following the prescriptions of the HVM, the total matter density field is
given by the sum of the halo and void contributions:

$\rho(\textbf{x})=\sum_{i}^{\text{halos}}\rho_{h}(\textbf{x}-\textbf{x}_{i}\,|\,M_{i})+\sum_{j}^{\text{voids}}\rho_{v}(\textbf{x}-\textbf{x}_{j}\,|\,M_{j})\;,$
(3.1)

where $\rho_{h}(\textbf{x}-\textbf{x}_{i}\,|\,M_{i})$ is the density profile
of a halo with mass $M_{i}$ centered at $x_{i}$ and
$\rho_{v}(\textbf{x}-\textbf{x}_{j}\,|\,M_{j})$ is the density profile of a
void with mass $M_{j}$ centered at $x_{j}$. Note how the HM is recovered when
the last term in Eq. (3.1) is neglected, i.e. if we set the matter density in
voids to zero across the whole Universe. We can rewrite Eq. (3.1) as:

$\begin{split}\rho(\textbf{x})=\int{\textrm{d}}M\int{\textrm{d}}^{3}x^{\prime}\Bigg{[}&\sum_{i}^{\text{halos}}\delta_{D}(M-M_{i})\delta_{D}(\textbf{x}^{\prime}-\textbf{x}_{i})\rho_{h}(\textbf{x}-\textbf{x}^{\prime}\,|\,M)\\\
+&\sum_{j}^{\text{voids}}\delta_{D}(M-M_{j})\delta_{D}(\textbf{x}^{\prime}-\textbf{x}_{j})\rho_{v}(\textbf{x}-\textbf{x}^{\prime}\,|\,M)\Bigg{]}\;.\end{split}$
(3.2)

The 2-point correlation function reads

$\xi(\textbf{r})=\frac{1}{\overline{\rho}_{m}^{2}}\langle\rho(\textbf{x})\rho(\textbf{x}+\textbf{r})\rangle-1\;,$
(3.3)

where $\overline{\rho}_{m}$ is the average matter density in the Universe.
Plugging Eq. (3.2) into Eq. (3.3), we get that Eq. (3.3) can be written as the
sum of three terms: a pure halo term (containing the correlation between
particles within one or two halos), a pure void term (containing the
correlation between particles within one or two voids) and a mixed term
(containing the correlation between two particles within a halo and a void,
respectively). In particular, the pure halo correlation reads:

$\begin{split}\xi_{h}(\textbf{r})&=\dfrac{1}{\overline{\rho}_{m}^{2}}\int{\textrm{d}}M_{1}{\textrm{d}}M_{2}{\textrm{d}}^{3}x_{1}{\textrm{d}}^{3}x_{2}\;\rho_{h}(\textbf{x}-\textbf{x}_{1}\,|\,M_{1})\rho_{h}(\textbf{x}-\textbf{x}_{2}+\textbf{r}\,|\,M_{2})\\\
&\times\left\langle\sum_{i,\,j}\delta_{D}(M_{1}-M_{i})\delta_{D}(M_{2}-M_{j})\delta_{D}(\textbf{x}_{1}-\textbf{x}_{i})\delta_{D}(\textbf{x}_{2}-\textbf{x}_{j})\right\rangle\;,\end{split}$
(3.4)

where we can distinguish two terms: the 1Halo term for particles residing
within the same halo and the 2Halo term for particles residing in different
halos:

$\displaystyle
I^{1H}=\delta_{D}(M_{1}-M_{2})\delta_{D}(\textbf{x}_{1}-\textbf{x}_{2})\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{1}}$
(3.5) $\displaystyle
I^{2H}=\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{1}}\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{2}}[1+\xi_{hh}(\textbf{x}_{1}-\textbf{x}_{2}\,|\,M_{1},\,M_{2})]\;,$
(3.6)

where $\textrm{d}n_{h}/\textrm{d}M$ denotes the halo mass function, which
measures the differential number density of halos in the mass range
$[M,\;M+{\textrm{d}}M]$:

$\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{1}}=\left\langle\sum_{i}\delta_{D}(M_{1}-M_{i})\delta_{D}(\textbf{x}_{1}-\textbf{x}_{i})\right\rangle$
(3.7)

and $\xi_{hh}(\textbf{x}_{1}-\textbf{x}_{2}\,|\,M_{1},\,M_{2})$ is the halo-
halo 2-point correlation function of halos with mass $M_{1}$ and $M_{2}$. A
similar discussion applies to the other terms of the 2-point correlation
function. However, for the cross-correlation the only non vanishing term is
naturally the one with particles residing in two distinct structures. The
terms appearing in Eq. (3.3) can be summarised as:

$\displaystyle\xi^{1H}(\textbf{r})$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\int{\textrm{d}}M\,\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M}\int{\textrm{d}}^{3}y\,\rho_{h}(\textbf{y}|M)\rho_{h}(\textbf{y}+\textbf{r}|M)$
(3.8) $\displaystyle\xi^{2H}(\textbf{r})$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\int{\textrm{d}}M_{1}\,\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{1}}b_{h}(M_{1})\int{\textrm{d}}M_{2}\,\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{2}}b_{h}(M_{2})$
(3.9)
$\displaystyle\times\int{\textrm{d}}^{3}y_{1}\,\rho_{h}(\textbf{y}_{1}|M_{1})\int{\textrm{d}}^{3}y_{2}\,\rho_{h}(\textbf{y}_{2}|M_{2})\xi^{L}(\textbf{y}_{1}-\textbf{y}_{2})$
$\displaystyle\xi^{1V}(\textbf{r})$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\int{\textrm{d}}M\,\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M}\int{\textrm{d}}^{3}y\,\rho_{v}(\textbf{y}|M)\rho_{v}(\textbf{y}+\textbf{r}|M)$
(3.10) $\displaystyle\xi^{2V}(\textbf{r})$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\int{\textrm{d}}M_{1}\,\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M_{1}}b_{v}(M_{1})\int{\textrm{d}}M_{2}\,\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M_{2}}b_{v}(M_{2})$
(3.11)
$\displaystyle\times\int{\textrm{d}}^{3}y_{1}\,\rho_{v}(\textbf{y}_{1}|M_{1})\int{\textrm{d}}^{3}y_{2}\,\rho_{v}(\textbf{y}_{2}|M_{2})\xi^{L}(\textbf{y}_{1}-\textbf{y}_{2})$
$\displaystyle\xi^{HV}(\textbf{r})$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\int{\textrm{d}}M_{1}\,\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{1}}b_{h}(M_{1})\int{\textrm{d}}M_{2}\,\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M_{2}}b_{v}(M_{2})$
(3.12)
$\displaystyle\times\int{\textrm{d}}^{3}y_{1}\,\rho_{h}(\textbf{y}_{1}|M_{1})\int{\textrm{d}}^{3}y_{2}\,\rho_{v}(\textbf{y}_{2}|M_{2})\xi^{L}(\textbf{y}_{1}-\textbf{y}_{2})\;,$

where we have used the linear approximation (tree-level) for the structure-
structure 2-point correlation ($\xi_{xy}$ for $x,y=h,\,v$). The density
contrast $\delta_{x}(\textbf{x}\,|\,M)$ is then obtained from the linear
matter density contrast $\delta_{m}^{L}(\textbf{x})$ as (see e.g. [38, 33]):

$\delta_{x}(\textbf{x}\,|\,M)=b^{L}_{x}(M)\delta_{m}^{L}(\textbf{x})\;,$
(3.13)

where $b^{L}_{x}(M)$ denotes the linear bias of $x$. For any combination of
structures $xy$ we have:

$\xi_{xy}(\textbf{r})=b_{x}(M_{x})\,b_{y}(M_{y})\,\xi^{L}(\textbf{r})\;,$
(3.14)

where we have dropped the $L$ apex on the linear bias to ease the notation.
Two conditions have to hold: the total matter density of the Universe has to
be equal to the sum of the total matter density in halos and voids, and matter
does not have to be biased with respect to itself. These two conditions lead
to the two following constraints:

$\displaystyle\overline{\rho}_{h}+\overline{\rho}_{v}=\overline{\rho}_{m}$
(3.15) $\displaystyle 1-\overline{b}_{h}-\overline{b}_{v}=0\;,$ (3.16)

where

$\displaystyle\overline{\rho}_{h}=\int{\textrm{d}}M\,M\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M}$
(3.17)
$\displaystyle\overline{\rho}_{v}=\int{\textrm{d}}M\,M\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M}\;,$
(3.18)

and

$\displaystyle\overline{b}_{h}=\frac{1}{\overline{\rho}_{m}}\int{\textrm{d}}M\,M\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M}b_{h}(M)$
(3.19)
$\displaystyle\overline{b}_{v}=\frac{1}{\overline{\rho}_{m}}\int{\textrm{d}}M\,M\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M}b_{v}(M)\;.$
(3.20)

Assuming spherical symmetry ($\rho_{x}(\textbf{r}\,|\,M)=\rho_{x}(r\,|\,M)$
for $x=h,\,v$) and Fourier transforming Eqs. (3.8)-(3.12) we have:

$\displaystyle P^{1H}(k)$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\int{\textrm{d}}M\;\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M}|\rho_{h}(k\,|\,M)|^{2}$
(3.21) $\displaystyle P^{2H}(k)$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\left[\int{\textrm{d}}M\;\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M}\rho_{h}(k\,|\,M)b_{h}(M)\right]^{2}P^{L}(k)$
(3.22) $\displaystyle P^{1V}(k)$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\int{\textrm{d}}M\;\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M}|\rho_{v}(k\,|\,M)|^{2}$
(3.23) $\displaystyle P^{2V}(k)$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\left[\int{\textrm{d}}M\;\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M}\rho_{v}(k\,|\,M)b_{v}(M)\right]^{2}P^{L}(k)$
(3.24) $\displaystyle P^{HV}(k)$
$\displaystyle=\frac{1}{\overline{\rho}^{2}_{m}}\int{\textrm{d}}M_{1}\;\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{1}}\rho_{h}(k\,|\,M_{1})b_{h}(M_{1})\int{\textrm{d}}M_{2}\;\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M_{2}}\,\rho_{v}(k\,|\,M_{2})b_{v}(M_{2})\,P^{L}(k)\;,$
(3.25)

where $P^{L}(k)$ is the linear matter power spectrum, given by the Fourier
transform of the linear 2-point correlation function $\xi^{L}(r)$. Note how
the HM is easily recovered when setting $\rho_{v}(k\,|\,M)=0$. The total
matter power spectrum is now given by:

$P(k)=P^{1H}(k)+P^{2H}(k)+P^{1V}(k)+P^{2V}(k)+2\,P^{HV}(k)$ (3.26)

and, using Eqs. (3.16) and (3.15), it can be shown that the total matter power
spectrum reduces to the linear matter power spectrum on very large scales
($k\ll 1$ Mpc-1). However, for the HM, these constraints return
$\overline{b}_{h}=1$ and $\overline{\rho}_{h}=1$, which, for a standard halo
mass function and bias, have a very slow convergence, requiring integration
down to tiny and untested masses. The HVM solves this shortcoming by taking
into account the matter within smaller halos in larger voids (see Sec. 3 and
[33] for a brief discussion on this point).

In the following sections we discuss the ingredients required to compute the
power spectrum: halo/void density profiles, mass functions and linear biases.
For the last two, we adopt the recipes provided by Voivodic et al. [33],
derived from the excursion set theory with two static barriers [81, 68, 33].
Hence, the obtained mass functions and linear biases naturally incorporate the
exclusion of voids inside halos and vice versa.

###  3.1 Density profiles

In this work, we adopt for halos the standard Navarro-Frenk-White (NFW)
profile [47], and for voids the Hamaus-Sutter-Wandelt (HSW) profile [83]. Both
assume a spherical density distribution of matter in the structure.

### Halos

The NFW density profile for a halo of virial mass $M$ can be parameterized as
follows [47]:

$\rho_{\text{NFW}}(r\,|\,M)=\frac{\rho_{s}}{c(M)\dfrac{r}{r_{\rm
vir}}\left(1+c(M)\dfrac{r}{r_{\rm vir}}\right)^{2}}\;,$ (3.27)

where $\rho_{s}$ is the characteristic density for which the volume-integrated
profile returns the virial mass, $r_{\rm vir}$ is the virial radius333The
virial radius is defined through
$M=\frac{4\pi}{3}\overline{\rho}_{m}\Delta_{\rm vir}r_{\rm vir}^{3}$, where
$\overline{\rho}_{m}$ is the average matter density of the Universe and
$\Delta_{\rm vir}\sim 330$ is the halo density contrast at the time of
virialization for the fiducial cosmology. and $c(M)$ is the concentration
parameter [84]. Its Fourier transform, truncated at $r_{\rm vir}$ is given by:

$\rho_{h}(k\,|\,M)=4\pi\int^{r_{\rm
vir}}_{0}{\textrm{d}}r\;r^{2}\,\frac{\sin{kr}}{kr}\,\rho_{\text{NFW}}(r\,|\,M)\;.$
(3.28)

### Voids

Figure 1:  (_Left_): Hamaus-Sutter-Wandelt void density contrast [74], at
$z=0$, for different void radii. Smaller voids are more underdense in their
central region but show higher compensation walls in their outskirt, i.e. they
tend to be overcompensated. Larger voids are less underdense in their central
region but show lower compensation walls in their outskirt, i.e. they tend to
be undercompensated. In the small box we show the integrated overdensity, as a
function of $r/r_{v}$, for a normalized HSW profile, so that the void
overdensity $\Delta_{v}\simeq 0.2047$ is reached, for any void, at its radius.
This allows us to consider spherical voids with fixed void overdensity.
(_Right_): Void density contrast for the profile proposed in Ref. [33], based
on empty centers and no compensation walls.

Given their low density, vastness and unsphericity at some scales, voids are
troublesome to study both with observations and simulations. Most void finders
show them as deeply underdense in their interiors444Voids of smaller size show
emptier central regions, while larger voids tend to be slightly denser in
their centers [74]., and the profiles exhibit overdense compensation walls
with a maximum located slightly outside their effective radii, shifting
outwards for larger voids. The height of the compensation wall decreases with
void size, causing the inner profile slope to become shallower and the wall to
widen. This trend divides all voids into being either overcompensated or
undercompensated, depending on whether the total mass within their
compensation wall exceeds or falls behind their missing mass in the center,
respectively. Ultimately, at sufficiently large distances to the void center,
all profiles approach the mean background density. The HSW profile [74] is a
simple empirical formula that can accurately capture the properties described
above:

$\frac{\rho_{v}(r\,|\,r_{v})}{\overline{\rho}_{m}}-1=\delta_{c}\frac{1-\left(\dfrac{r}{r_{s}}\right)^{\alpha}}{1+\left(\dfrac{r}{r_{v}}\right)^{\beta}}\;,$
(3.29)

where $r_{s}$ is the characteristic radius for which
$\rho_{v}=\overline{\rho}_{m}$, $r_{v}$ is the effective radius of the
considered void and $(\alpha,\;\beta,\;\delta_{c})$ are free parameters. In
particular, $\delta_{c}$ represents the central density contrast and $\alpha$
and $\beta$ are the inner and outer slopes of the compensation wall. These
parameters are usually determined through best-fits of N-body dark matter
simulations but Hamaus et al. [74] also provide empirical functions, which we
use in our analysis, relating them to $r_{s}$ and $r_{v}$, therefore
decreasing the degrees of freedom.

Nevertheless, the choice of the profile strongly depends on the void finder,
which N-body simulations should be calibrated on. Other works have used void
finders that neither rely on any central particle, but on Voroni vertices, nor
show any compensation wall (see e.g. [33]). These properties can be summarized
in the one-parameter empirical formula proposed in [33]:

$\frac{\rho_{v}(r\,|\,r_{v})}{\overline{\rho}_{m}}-1=\frac{1}{2}\left[1+\tanh{\left(\frac{y-y_{0}}{s(r_{v})}\right)}\right]-1\;,$
(3.30)

where $y=\ln{(r/r_{v})}$ and $y_{0}=\ln{(r_{0}/r_{v})}$. The radius $r_{0}$
can be parameterized in terms of $s$, which remains the only free quantity.
Ref. [33] shows, through N-body simulations, that its dependence on $r_{v}$ is
very weak and can be safely fixed to $s=0.75$, for all void radii.

The left and right panel of Fig. 1 show the void density contrast described by
Eq. (3.29) and Eq. (3.30), respectively, for multiple bins in void radii, and
follow efficiently all of the related properties outlined above.

The Fourier transform of the void profile is given by:

$\rho_{v}(k\,|\,M(r_{v}))=4\pi\int^{r_{v}}_{0}{\textrm{d}}r\;r^{2}\,\frac{\sin{kr}}{kr}\,\rho_{v}(r\,|\,M(r_{v}))\;,$
(3.31)

where, in principle,
$M(r_{v})=4\pi\int_{0}^{r_{v}}{\textrm{d}}r\;r^{2}\,\rho_{v}(r\,|\,r_{v})$.
However, in the following discussion we will consider spherical voids with an
overdensity $\Delta_{v}\simeq 0.2047$ (see e.g. [80]), whose mass can be
computed through
$M(r_{v})=\frac{4\pi}{3}\,\overline{\rho}_{m}\Delta_{v}r^{3}_{v}$. Hence, the
profile needs to be properly normalized in order for voids of any radius to
reach the chosen overdensity $\Delta_{v}$:

$\rho_{v}(r\,|\,r_{v})\longrightarrow\frac{\Delta_{v}}{\Delta(r_{v})}\rho_{v}(r\,|\,r_{v})\;,$
(3.32)

where

$\Delta(r_{v})=\frac{3}{r_{v}^{3}}\int_{0}^{r_{v}}{\textrm{d}}r\;r^{2}\,\frac{\rho_{v}(r\,|\,r_{v})}{\overline{\rho}_{m}}\;.$
(3.33)

Note how, by doing so, Eq. (3.31) can be rewritten in terms of the normalized
profile, but the integration should be performed up to a void radius extracted
from a different $M(r_{v})$ relation, in order to obtain the same value for
the Fourier transform of the void profile. Calculating the latter through the
normalized profile is only an expedient to ease the computation, nonetheless,
the physical profile remains that shown in either panels of Fig. 1, which
follow the void properties previously discussed. The small box in the left
panel of Fig. 1 shows the evolution of Eq. (3.33) as a function of the
integration upper limit, for the HSW normalized profile. Independently of the
void size, the required void overdensity $\Delta_{v}$ is reached when
integrating up to the void radius.

Studies of voids’ properties are yet at their early stages and the choice of
the void profile usually depends on the used void finder. However, while for
definiteness in the analyses shown below we adopt the HSW profile,
nevertheless we have verified that the actual choice of the profile has only a
sub-percent effect on the power spectra (see Sec. 4), and we leave further
discussion on the use of different profiles, in the context of our work, to
future studies.

###  3.2 Mass functions

The halo/void mass functions used in the HVM are computed through the
excursion set theory with two static linear barriers (2SB) [68, 33] to avoid
double counting of matter. We show here a comparison of these mass functions
with those computed following Sheth-Tormen [85], for halos, and standard
Press-Schechter theory [86], for voids.

A central ingredient is the variance of the linear density field smoothed at
some scale $R$ (with a top-hat window function $\widehat{W}$):

$\sigma^{2}(R)=\int\frac{{\textrm{d}}k}{2\pi^{2}}\;k^{2}P^{L}(k)\widehat{W}^{2}(kR)\;,$
(3.34)

where for any structure we have:

$M=\frac{4\pi}{3}\,\overline{\rho}_{m}R^{3}\;.$ (3.35)

### Halos

Figure 2:  (_Left_): Halo mass function as a function of halo mass. (_Right_):
Void mass function as a function of void radius/mass. Both results are
computed at $z=0$. We compare two different models: Sheth-Tormen for halos or
Press-Schecther for voids (purple) and 2SB (green). The insets show how, at
very small masses, the excursion set prediction decreases rapidly, differently
from the standard Sheth-Tormen or Press-Schecther abundances.

The halo mass function can be written as:

$\frac{{\textrm{d}}n_{h}}{{\textrm{d}}\ln{M}}=f_{h}(\sigma)\frac{\overline{\rho}_{m}}{M}\frac{{\textrm{d}}\ln{\sigma^{-1}}}{{\textrm{d}}\ln{M}}\;,$
(3.36)

where $f_{h}(\sigma)$, the multiplicity function, determines the model. In
this work, we consider the multiplicity functions:

$\displaystyle f_{h}^{\text{sheth99}}(\sigma)$
$\displaystyle=2A\nu\left(1+\frac{1}{(a\nu)^{p}}\right)\sqrt{\frac{a}{2\pi\nu}}\;\text{exp}\left(-\frac{a\nu}{2}\right)$
(3.37) $\displaystyle f_{h}^{\text{2SB}}(\sigma)$
$\displaystyle=2\,\sum_{n}\frac{n\pi}{\delta_{T}^{2}}\,\sigma^{2}\,\sin{\left(\frac{n\pi\delta_{c}}{\delta_{T}}\right)}\,\text{exp}\left[-\frac{n^{2}\pi^{2}}{2\delta_{T}^{2}}\,\sigma^{2}\right]\;,$
(3.38)

where $\delta_{T}=\delta_{c}+|\delta_{v}|$ with $\delta_{c}\simeq 1.686$ and
$\delta_{v}\simeq-2.717$ the linear critical (over and under) densities for
halo and void formation, respectively. Here $f_{h}^{\text{sheth99}}$ is the
Sheth-Tormen multiplicity function [85], with
$\nu=(\delta_{c}^{2}/\sigma^{2})$ and $(A,a,p)$ free parameters.
$f_{h}^{2\text{SB}}$ is the excursion set prediction for a model with two
static barriers.

Note how $f_{h}^{\text{sheth99}}$ is normalized to unity, whereas
$f_{h}^{2\text{SB}}$ is not, since it already accounts for the existence of
voids. As a direct consequence, the excursion set predictions with two
barriers tend to decrease rapidly at small masses, as shown by the left panel
of Fig. 2. This directly reflects on no need to integrate Eq. (3.17) down to
very low and untested masses.

### Voids

The void mass function is usually expressed in terms of its radius, rather
than its mass:

$\frac{{\textrm{d}}n_{v}}{{\textrm{d}}\ln{r}}=\Delta_{v}\left[\frac{f_{v}(\sigma)}{V(r_{L})}\frac{{\textrm{d}}\ln{\sigma^{-1}}}{{\textrm{d}}\ln{r_{L}}}\right]_{r_{L}=\,r/1.7}\;,$
(3.39)

where $\Delta_{v}=({r_{L}}/{r_{v}})^{1/3}\simeq 0.2$ is the void overdensity
and $f_{v}(\sigma)$ determines the model. The last equation follows from
Eq.(3.36) and derives from the conservation requirement of the volume density,
as discussed by Jennings et al. [82].

In this work, we consider the following multiplicity functions:

$\displaystyle f_{v}^{\text{PS}}(\sigma)$
$\displaystyle=\sqrt{\frac{2}{\pi}}\frac{|\delta_{v}|}{\sigma}\;\text{exp}^{-\frac{\delta^{2}_{v}}{2\sigma^{2}}}$
(3.40) $\displaystyle f_{v}^{\text{2SB}}(\sigma)$
$\displaystyle=2\,\sum_{n}\frac{n\pi}{\delta_{T}^{2}}\,\sigma^{2}\,\sin{\left(\frac{n\pi|\delta_{v}|}{\delta_{T}}\right)}\,\text{exp}\left[-\frac{n^{2}\pi^{2}}{2\delta_{T}^{2}}\,\sigma^{2}\right]\;,$
(3.41)

where $\delta_{c}$, $\delta_{v}$ and $\delta_{T}$ are the same parameters
discussed for the halo mass function. Here $f_{v}^{\text{PS}}$ is the Press-
Schechter multiplicity function [86] and $f_{v}^{2\text{SB}}$ is the excursion
set prediction for a model with two static barriers.

In the right panel of Fig. 2 we compare the Press-Schechter and 2SB void mass
functions. Once again, we note how the excursion set prediction decreases
rapidly at small masses (differently from Press-Schechter theory), avoiding
the need to integrate Eq. (3.18) down to tiny untested masses.

###  3.3 Linear bias

We present here the linear biases used in this work. In particular, we show
the excursion set predictions [81, 80, 82, 68] proposed by Voivodic et al.
[33]. Using these biases makes the model fully self-consistent with the 2SB
mass functions. Moreover, the combination of 2SB bias and mass function avoids
double counting of matter, excluding the overlap of structures.

### Halos

The relevant functions for halos that we use are:

$\displaystyle b_{h}^{\text{sheth99}}(\sigma)$
$\displaystyle=1+\frac{a\nu-1}{\delta_{c}}+\frac{2p}{\delta_{c}\left(1+(a\nu)^{p}\right)}$
(3.42) $\displaystyle b^{2\text{SB}}_{h}(\sigma)$
$\displaystyle=1-\frac{\sum_{n}\frac{n\pi}{\delta_{T}^{2}}\,\sin{\left(\frac{n\pi\delta_{c}}{\delta_{T}}\right)}\,\text{exp}\left[-\frac{n^{2}\pi^{2}}{2\delta_{T}^{2}}\,\sigma^{2}\right]\left[\textrm{cotan}{\left(\frac{n\pi\delta_{c}}{\delta_{T}}\right)\frac{n\pi}{\delta_{T}}}\right]}{\sum_{n}\frac{n\pi}{\delta_{T}^{2}}\,\sin{\left(\frac{n\pi\delta_{c}}{\delta_{T}}\right)}\,\text{exp}\left[-\frac{n^{2}\pi^{2}}{2\delta_{T}^{2}}\,\sigma^{2}\right]}\;.$
(3.43)

Here $b_{h}^{\text{sheth99}}$ is the Sheth-Tormen linear bias [87], with
$\nu=\delta_{c}^{2}/\sigma^{2}$ and $(a,\,p)$ are the same free parameters of
the Sheth-Tormen mass function (Eq. (3.37)). The quantity $b_{h}^{2\text{SB}}$
is the prediction for a model with two static barriers.

In the left panel of Fig. 3, we compare the linear bias predictions for the
functions above. Voivodic et al. [33] show how the 2SB linear bias is in very
good agreement with simulations.

Figure 3:  (_Left_): Halo linear bias as a function of halo mass. (_Right_):
Void linear bias as a function of void radius/mass. Both results are computed
at $z=0$. We compare two different models: Sheth-Tormen for halos or Press-
Schechter for voids (purple) and 2SB (green). For voids, we show the absolute
value of the bias, as the latter is usually negative for larger voids. This is
strictly connected to the absence of compensation walls for very large voids.

### Voids

In this work we consider the following functions:

$\displaystyle b_{v}^{\text{PS}}(\sigma)$
$\displaystyle=1+\frac{\nu-1}{\delta_{v}}$ (3.44) $\displaystyle
b^{2\text{SB}}_{v}(\sigma)$
$\displaystyle=1+\frac{\sum_{n}\frac{n\pi}{\delta_{T}^{2}}\,\sin{\left(\frac{n\pi|\delta_{v}|}{\delta_{T}}\right)}\,\text{exp}\left[-\frac{n^{2}\pi^{2}}{2\delta_{T}^{2}}\,\sigma^{2}\right]\left[\frac{n\pi}{\delta_{T}}\textrm{cotan}{\left(\frac{n\pi|\delta_{v}|}{\delta_{T}}\right)}\right]}{\sum_{n}\frac{n\pi}{\delta_{T}^{2}}\,\sin{\left(\frac{n\pi|\delta_{v}|}{\delta_{T}}\right)}\,\text{exp}\left[-\frac{n^{2}\pi^{2}}{2\delta_{T}^{2}}\,\sigma^{2}\right]}\;,$
(3.45)

where $\delta_{c}$, $\delta_{v}$ and $\delta_{T}$ are the same discussed for
the halo mass function. Here, $b_{v}^{\text{PS}}$ is the Press-Schechter
linear bias [86], with $\nu={\delta_{v}^{2}}/{\sigma^{2}}$. The quantity
$b_{v}^{2\text{SB}}$ is the excursion set prediction for a model with two
static barriers.

In the right panel of Fig. 3, we compare the linear bias predictions for the
functions above. Note that we show the absolute value of the bias since the
latter is negative for larger voids in the 2SB model. This anti-correlation
with the matter field is tightly connected with the compensation walls, as
larger voids tend to have smaller compensation walls (see Sec. 3.1) and,
therefore, to be under-compensated, resulting in a negative sign for the
linear bias.

### Convergence at small masses

As we know, the total power spectrum should follow linear theory on large
scales. This is naturally achieved when the constraints of Eqs. (3.17) and
(3.16) are satisfied. The former ensures that the total matter contained in
all structures matches the total matter in the Universe, while the latter
ensures that matter is not biased with respect to itself.

These constraints can be rewritten as:

$\displaystyle I^{\rho}$
$\displaystyle=\frac{\overline{\rho}}{\overline{\rho}_{m}}=\frac{1}{\overline{\rho}_{m}}\int_{0}^{\infty}{\textrm{d}}M\;M\,\frac{{\textrm{d}}n}{{\textrm{d}}M}=1$
(3.46) $\displaystyle I^{b}$
$\displaystyle=\overline{b}=\frac{1}{\overline{\rho}_{m}}\int_{0}^{\infty}{\textrm{d}}M\;M\,\frac{{\textrm{d}}n}{{\textrm{d}}M}\,b(M)=1\;,$
(3.47)

where both integrals are given by the sum of the contributions from all of the
considered structures.

Ref. [33] studies the problem by varying the lower bound of the integrals in
Eqs. (3.46) and (3.47) for mass values $M_{\rm min}>0$. They show that the
standard HM (with Sheth-Tormen-like halo mass function and linear bias), does
not converge to unity even down to very low masses ($M_{\rm
min}<10^{4}\;\text{M}_{\odot}$), pointing out the difficulty of the HM to
recover efficiently the linear matter power on large scales and implying the
need to either integrate Eqs. (3.46) and (3.47) down to very low masses or
suitably normalize the 2Halo term. On the contrary, the sum of contributions
from halos and voids, within the HVM, to Eqs. (3.46) and (3.47) converges to
unity already for $M_{\rm min}\simeq 10^{9}\;\text{M}_{\odot}$. Hence,
incorporating voids in a model of LSS eliminates the need for an exotic re-
normalization and saves up computational-time.

In conclusion, allowing matter to lie within multiple structures, through the
excursion set theory with two barriers [81, 68, 33], offers as a fully self-
consistent model, that correctly recovers the matter power on large scales,
and as an efficient tool to a variety of cosmological and astrophysical
studies. The foundation of such a success lies in considering the mass of
smaller halos within voids of larger size.

##  4 Power spectra

Early decoupling of dark matter allowed for its perturbations to grow
undisturbed and form deep potential wells that accreted baryonic matter after
recombination. As a consequence, galaxies are tightly coupled with dark matter
and offer as a powerful tracer of its distribution. At the same time, galaxies
host astrophysical sources that emit in several energy bands and dark matter
structures bend the light, from distant sources, in the weak lensing regime.
We then expect a correlation between gravitational tracers, such as i) the
galaxy distribution and ii) the cosmic shear, and gamma-ray emitters, such as
a) annihilating/decaying dark matter and b) unresolved astrophysical sources.
The latter, while being an interesting signal by itself in the study of the
unresolved component of the gamma-ray sky, nevertheless represent an
irreducible background for pure dark matter studies.

In this section we generalize the formalism introduced in the previous
sections to adapt the 3D power spectrum to the cross-correlation signal
between two different source fields. In particular, we are interested in
correlating i-ii) with a-b) within the HVM. For any couple of observables $i$
and $j$, the power spectra are:

$\displaystyle P^{1H}_{ij}(k)$
$\displaystyle=\int{\textrm{d}}M\;\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M}f_{i}^{h*}(k\,|\,M)f_{j}^{h}(k\,|\,M)$
(4.1) $\displaystyle P^{2H}_{ij}(k)$
$\displaystyle=\int{\textrm{d}}M_{1}\;\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{1}}f_{i}^{h*}(k\,|\,M_{1})b_{h}(M_{1})\int{\textrm{d}}M_{2}\;\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{2}}f_{j}^{h}(k\,|\,M_{2})b_{h}(M_{2})P^{L}(k)$
(4.2) $\displaystyle P^{1V}_{ij}(k)$
$\displaystyle=\int{\textrm{d}}M\;\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M}f_{i}^{v*}(k\,|\,M)f_{j}^{v}(k\,|\,M)$
(4.3) $\displaystyle P^{2V}_{ij}(k)$
$\displaystyle=\int{\textrm{d}}M_{1}\;\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M_{1}}f_{i}^{v*}(k\,|\,M_{1})b_{v}(M_{1})\int{\textrm{d}}M_{2}\;\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M_{2}}f_{j}^{v}(k\,|\,M_{2})b_{v}(M_{2})P^{L}(k)$
(4.4) $\displaystyle P^{HV}_{ij}(k)$
$\displaystyle=\int{\textrm{d}}M_{1}\;\frac{{\textrm{d}}n_{h}}{{\textrm{d}}M_{1}}f_{(i}^{h*}(k\,|\,M_{1})b_{h}(M_{1})\int{\textrm{d}}M_{2}\;\frac{{\textrm{d}}n_{v}}{{\textrm{d}}M_{2}}\,f_{j)}^{v}(k\,|\,M_{2})b_{v}(M_{2})\,P^{L}(k)\;,$
(4.5)

where we have generalized Eqs. (3.21)-(3.25) for any pair of source fields
$f_{i},f_{j}$. The Halo-Void term in Eq. (4.5) is given by the sum of the
symmetric permutations of the considered fields. In the notation, we have
omitted the redshift dependence of the 3D power spectrum for simplicity, but
it has been considered throughout the computation.

The main idea of this work relies on the observational ability to identify
cosmic voids: this would then allow to perform the cross-correlations between
the relevant gravitational tracers (related to halos or voids) and the
corresponding gamma-ray emission (due to dark matter annihilation/decay or the
astrophysical sources) from those structures. In this case, we can split the
CAPS (2.2) into two terms, each one depending on the 3D power spectrum
generated only from halos or voids:

$\displaystyle P_{ij}^{h}(k,\;z)$
$\displaystyle=P_{ij}^{1H}(k,\;z)+P_{ij}^{2H}(k,\;z)$ (4.6) $\displaystyle
P_{ij}^{v}(k,\;z)$ $\displaystyle=P_{ij}^{1V}(k,\;z)+P_{ij}^{2V}(k,\;z)\;.$
(4.7)

In order to calculate the power spectra of Eqs. (4.6), (4.7) and the window
functions of Eq. (2.2), we need to specify the relation between gravitational
tracers and the underlying large-scale structure as well as how the DM
particles and astrophysical sources hosted in halos and voids produce the
gamma-ray emission observed by our telescopes. In our work we consider cosmic
shear and galaxy catalogs as our reference gravitational tracers. Clearly, in
order to calculate the total signal, we need to add also the HV term.

Concerning the relation between the fields that define our signals with the
underlying halo or void mass distribution, let us start with weak lensing,
which directly depends on the distribution of dark matter across the line of
sight; therefore the cosmic shear source field is intuitively proportional to
the density profile of dark matter structures. Similarly, decaying dark matter
is sourced by the distribution of dark matter itself, as it requires a single
particle to occur in the decaying process, and therefore it directly traces
the mass density. In both cases we have:

$f^{x}_{l}(k,\;z\,|\,M)=f^{x}_{d}(k,\;z\,|\,M)=\frac{\mathcal{F}[\rho_{x}](k,\;z\,|\,M)}{\overline{\rho}_{m}(z)}\;,$
(4.8)

where $l$ stands for "lensing", $d$ for "decaying dark matter", $x=(h,\;v)$
denotes the considered structure and $\mathcal{F}[\varphi]$ denotes the
Fourier transform of the field $\varphi$, which in our case can be either the
density $\rho_{x}$ for decaying dark matter or the density squared
$\rho^{2}_{x}$ for annihilating dark matter.

Also the galaxy distribution is sourced by the density of dark matter
structures, however the latter must be weighted through the halo occupation
distribution (HOD) [38, 88, 89, 90, 91, 5], namely
$\langle\tilde{N}_{g}\rangle$, which is an indicator of the number of galaxies
in each region of a dark matter halo:

$f^{x}_{g}(k,\;z\,|\,M)=\frac{\langle\tilde{N}_{g}^{x}\rangle}{\overline{n}_{g}^{x}(z)}=\frac{\langle
N_{\text{cen}}(M)\rangle+\langle
N_{\text{sat}}(M)\rangle\mathcal{F}[\rho_{x}](k,\;z\,|\,M)/M}{\overline{n}_{g}^{x}(z)}\;,$
(4.9)

where
$\overline{n}_{g}^{x}=\int{\textrm{d}}M\;\dfrac{{\textrm{d}}n_{x}}{{\textrm{d}}M}\left(\langle
N_{\text{cen}}\rangle+\langle N_{\text{sat}}\rangle\right)$ the average number
of galaxies. Following [5], and references therein, the average number of
central galaxies $\langle N_{\text{cen}}\rangle$ and satellite galaxies
$\langle N_{\text{sat}}\rangle$ in halos can be modeled as:

$\displaystyle\langle N_{\text{cen}}\rangle$
$\displaystyle=\frac{1}{2}\left[1+\textrm{erf}\left(\frac{\log M-\log M_{\rm
th}}{\sigma_{\log M}}\right)\right]$ (4.10) $\displaystyle\langle
N_{\text{sat}}\rangle$
$\displaystyle=\left(\frac{M}{M_{*}}\right)^{\alpha}\textrm{exp}\left[-\frac{M_{\rm
cut}}{M}\right]\;,$ (4.11)

where $M_{\rm th}$ denotes the approximate halo mass required to populate the
halo with the considered type of galaxies and $\sigma_{\log M}$ governs the
width of the transition between $0$ and $1$ for the central galaxy. The
satellite distribution is described by a power law (of order $\alpha$) with an
exponential cutoff $M_{\rm cut}$ at low masses. In this paper we consider the
galaxy distribution from the 2MASS catalog [92], as it is one of the most
extended almost-all-sky catalogs. Exploiting the HOD results of Ref. [93], we
adopt a step function for central galaxies ($\langle N_{\text{cen}}\rangle=0$
for $M<10^{12.1}\;\text{M}_{\odot}$ and $\langle N_{\text{cen}}\rangle=1$ for
$M\geq 10^{12.1}\;\text{M}_{\odot}$) and take the following parameters for
satellite galaxies: $\alpha=1.2$, $M_{*}=10^{13.5}\;\text{M}_{\odot}$ and
$M_{\rm cut}=0$.

Figure 4:  Comparison of the total cross-correlation 3D power spectrum between
gamma-ray sources and gravitational tracers (at $z=0.5$): cosmic shear
(_Left_) and the 2MASS galaxy catalog (_Right_). The shown gamma-ray sources
are DM decay (solid black), DM annihilation (dashed black), BL Lac (magenta),
FSRQ (yellow), mAGN (cyan) and SFG (crimson red).

In the case of voids, we need to adapt the HOD to a void occupation
distribution (VOD). In this regard, we assume that the total number of
galaxies is the same within halos and voids of the same mass. Intuitively, if
the local density fractions of the Universe components remain unchanged in
voids, we expect the fraction of matter contained in galaxies to be the same
as in halos. Ref. [94] shows that voids usually contain bluer galaxies,
however, the general properties (such as color distribution, bulge to total
ratios, and concentrations) are remarkably similar to those in halos. The main
difference resides in their distribution within the structure; in particular,
the number of galaxies tends to increase with the radius of the void (see e.g.
[74]), while void centers are expected to be extremely underdense, leaving
arguably no space for central galaxies555Indeed, void finders are usually
based on tessellation methods around void centers, found as the minima in the
density distribution, thus representing points sensibly far from the closest
galaxies.. This can be accounted for by re-scaling the void density profile
with a parameterization of the galaxy number density contrast (in voids) shown
in [94]. We define the number density contrast of galaxies as
$1+\delta_{g}(r)=n_{g}(r)/\overline{n}_{v}$ (where $n_{g}(r)$ is the radial
number density of galaxies and $\overline{n}_{v}$ is the average density of
the considered void). The density of galaxies within voids can be rewritten as

$\rho_{v}^{\rm
gal}(r)=\rho_{v}(r)\frac{n_{g}(r)}{\overline{n}_{v}}\equiv\rho_{v}(r)(1+\delta_{g}(r))\;,$
(4.12)

where the values for $\delta_{g}(r)$ have been extrapolated by the data of
Ref. [94]. By plugging this expression in the HOD formula shown above and
setting $\langle N_{\text{cen}}(M)\rangle=0$ for any void we obtain the needed
VOD.

In the case of annihilating dark matter, the process requires both a particle
and an antiparticle to occur. The source field is therefore proportional to
the square of the density profile of dark matter structures:

$f^{x}_{a}(k,\;z\,|\,M)=\frac{1}{\Delta_{x}^{2}(z)}\frac{\mathcal{F}[\rho_{x}^{2}](k,\;z\,|\,M)}{\overline{\rho}^{2}_{m}(z)}\;,$
(4.13)

where $a$ stands for "annihilating dark matter" and $\Delta_{x}^{2}(z)$ is the
clumping factor, which acts as a transfer function between
$\langle\rho^{2}\rangle$ and $\langle\rho\rangle^{2}$:

$\Delta^{2}_{x}(z)=\frac{\langle\rho_{x}^{2}\rangle}{\overline{\rho}^{2}_{m}}=\int{\textrm{d}}M\;\frac{{\textrm{d}}n_{x}}{{\textrm{d}}M}\int{\textrm{d}}^{3}x\;\frac{\rho_{x}^{2}(\textbf{x}\,|\,M)}{\overline{\rho}^{2}}\;.$
(4.14)

When considering astrophysical sources we have to replace the structure mass
$M$ with the luminosity $L$ of the source, the mass function
${\textrm{d}}n_{x}/{\textrm{d}}M$ with the GLF
$\phi_{x}(L,\;z)={\textrm{d}}n_{x}/{\textrm{d}}L$ and the linear bias with
$b_{s}^{x}(M(L),\;z)$, following [3]. In this work we consider four different
classes of unresolved gamma-ray astrophysical sources: BL Lacertae objects (BL
Lacs), flat-spectrum radio quasars (FSRQs), misaligned AGN (mAGN) and star-
forming galaxies (SFGs). The GLFs and mass-to-luminosity functions $M(L)$ of
each astrophysical class are shown in Appendix A and B, respectively. The
source field can be expressed as:

$f_{s}^{x}(z\,|\,L)=\frac{L}{\langle g^{x}_{s}\rangle}\;,$ (4.15)

where $\langle g^{x}_{s}(z)\rangle=\int{\textrm{d}}L\;L\,\phi_{x}(L,\;z)$ is
the average luminosity density of sources.

Gamma-ray astrophysics has deeply improved in the last decade, with thousands
of sources detected, allowing for a better understanding of the most violent
phenomena of the Universe. However, the unresolved gamma-ray background (UGRB)
still remains to be fully understood. Moreover, previous works mostly focus on
unresolved sources within dark matter halos, while their presence and
properties within the most underdense regions of the Universe remain somewhat
foggy and require further studies. Works on the properties of galaxies [94]
and AGN [95, 96] in voids show that the environment makes very little impact
on general properties. AGN are slightly more common in underdense regions, but
only for the most luminous galaxies. At the same time, voids are generally
younger than halos, due to the hierarchy in structure formation, showing bluer
galaxies and increased rates of structure formation (the effect being evident
only in the rare very massive galaxies). Nevertheless, the accretion rates do
not show particular changes and we expect mostly negligible variations between
the luminosity properties of sources hosted by galaxies in halos and those in
voids of the same mass. Given these considerations, we assume the number of
sources within voids and halos of the same mass to be constant, which is
intuitively similar to what we have done for the void occupation distribution.
Under this assumption, we can write the void GLF as the one of halos re-scaled
by the volume fraction at fixed mass $M_{*}$666Dimensionally speaking, the GLF
is a number density per unit of luminosity; assuming no dependence of the
luminosity properties on the environment and the number of sources to be
constant among different structures of the same mass, the GLF goes as the
inverse of a volume.:

$\phi_{v}(L,\;z)=\frac{V_{h}(M_{*})}{V_{v}(M_{*})}\,\phi_{h}(L,\;z)\sim
6.2\times 10^{-4}\,\phi_{h}(L,\;z)\;,$ (4.16)

where $V_{h}(M_{*})/V_{v}(M_{*})=\Delta_{v}/\Delta_{h,\,\text{vir}}$.

In the left panel of Fig. 4, we show the total 3D power spectrum of the cross-
correlation between cosmic shear and gamma-ray emitters (decaying/annihilating
DM and the four classes of unresolved astrophysical sources discussed above)
at $z=0.5$, computed within the HVM through the sum of contributions from Eqs.
(4.1)-(4.5). We discuss, with more detail, the importance of each term in the
Appendix. The same result is shown in the right panel of Fig. 4 for the cross-
correlation of the 2MASS galaxy catalog and the unresolved gamma-ray sky.

##  5 Window functions of gamma-ray emission, galaxies and cosmic shear

| $\Gamma$ |  $L_{\rm min}$ $[$erg s${}^{-1}]$ |  $L_{\rm max}$ $[$erg s${}^{-1}]$
---|---|---|---
BL Lacs [97] | $2.11$ | $7\times 10^{43}$ | $10^{52}$
FSRQ [98] | $2.44$ | $10^{44}$ | $10^{52}$
mAGN [99, 100] | $2.37$ | $10^{40}$ | $10^{50}$
SFG [101] | $2.7$ | $10^{37}$ | $10^{42}$
Table 1: Spectral index $\Gamma$, minimum and maximum luminosities and related
references for the classes of unresolved gamma-ray astrophysical sources
considered in this paper.

The last needed ingredient to compute the CAPS in Eq. (2.2) is the window
function for the different observables considered here.

The window function of weak lensing takes the form (see e.g. [102]):

$W_{l}(\chi)=\frac{3}{2}\frac{H_{0}^{2}}{c^{2}}\Omega_{m}(1+z)\chi\int_{\chi}^{\infty}{\textrm{d}}\chi^{\prime}\;\frac{\chi^{\prime}-\chi}{\chi^{\prime}}\frac{{\textrm{d}}N}{{\textrm{d}}\chi^{\prime}}(\chi^{\prime})\;,$
(5.1)

where ${\textrm{d}}N/{\textrm{d}}\chi$ denotes the redshift distribution of
the background sources (see [103]), normalized to unit area.

For galaxies, the window function simply reduces to their redshift
distribution (e.g. see [4]):

$W_{g}(\chi)=\frac{{\textrm{d}}N_{g}}{{\textrm{d}}z}\frac{{\textrm{d}}z}{{\textrm{d}}\chi}=\frac{{\textrm{d}}N_{g}}{{\textrm{d}}z}\frac{H(z)}{c}\;.$
(5.2)

As mentioned above, we consider the galaxy distribution from the 2MASS catalog
[92].

In the case of decaying dark matter we have:

$W_{d}(E,\;z)=\frac{1}{4\pi}\frac{\Omega_{\text{DM}}\rho_{c}}{m_{\text{DM}}\tau_{d}}\frac{{\textrm{d}}N_{d}}{{\textrm{d}}E}\left[E(1+z)\right]\text{e}^{-\tau\left[E(1+z),\;z\right]}\;,$
(5.3)

where $\Omega_{\text{DM}}$ is the cosmological abundance of DM and $\rho_{c}$
the critical density of the Universe today. $m_{\text{DM}}$ and $\tau_{d}$
denote the mass and decay lifetime of the DM particle (here we consider
$m_{\text{DM}}=100\;\text{GeV}$ and $\tau_{d}=3\times 10^{27}\;\text{s}$, from
[104, 105], for definiteness, about at its conservative lower bound) and
${\textrm{d}}N_{d}/{\textrm{d}}E$ is the number of photons emitted by the
decay in the energy band $[E,\;E+{\textrm{d}}E]$, here for the decaying
channel $b\overline{b}$ (see e.g. [106]). Finally, $\tau$ is the optical depth
for absorption on the line of sight (taken here from [107]). The latter is
relevant for gamma-rays mainly due to pair production on the extra-galactic
background light emitted by galaxies in the ultraviolet, optical, and infrared
bands.

Figure 5:  (_Left_): Comparison between the window functions of
decaying/annihilating dark matter with those of the 2MASS galaxy distribution
and weak lensing. (_Right_): Comparison between the window functions of each
class of gamma-ray astrophysical sources and those of gravitational tracers.
In both cases the functions are computed for gamma-rays with $E_{\gamma}=5$
GeV, and normalized to the redshift integrated intensity.

For the case of annihilating dark matter, since the source field depends on
$\rho^{2}$, the clumping factor of Eq. (4.14) enters the window function,
leading to a dependence on the structure $x$, differently from the decaying
case:

$W_{a}^{x}(E,\;z)=\frac{\left(\Omega_{\text{DM}}\rho_{c}\right)^{2}}{4\pi}\frac{\langle\sigma_{a}v\rangle}{2m_{\text{DM}}^{2}}(1+z^{3})\Delta_{x}^{2}(z)\frac{{\textrm{d}}N_{a}}{{\textrm{d}}E}\left[E(1+z)\right]\text{e}^{-\tau\left[E(1+z),\;z\right]}\;,$
(5.4)

where $\langle\sigma_{a}v\rangle$ is the velocity averaged annihilating cross
section (we consider the thermal value of $3\times
10^{-26}\;\text{cm}^{3}\text{s}^{-1}$ for a WIMP dark matter candidate with a
mass around $100\;\text{GeV}$, see [104]) and
${\textrm{d}}N_{a}/{\textrm{d}}E$ is the number of photons emitted by
annihilations in the energy band $[E,\;E+dE]$, we consider the channel
$b\overline{b}$ [106].

The window function of astrophysical sources, within a structure $x$, yields
(see e.g. [29]):

$W_{s}^{x}(E,\;z)=\left(\frac{d_{L}(z)}{1+z}\right)^{2}\int_{L_{\rm
min}}^{L_{\rm
max}(z)}{\textrm{d}}L\;\frac{{\textrm{d}}F}{{\textrm{d}}E}(E,\;L,\;z)\phi_{x}(L,\;z)\;,$
(5.5)

where $\phi_{x}$ is the GLF, $d_{L}(z)=(1+z)\,\chi(z)$ the luminosity distance
and ${\textrm{d}}F/{\textrm{d}}E$ the spectral energy distribution (SED)
which, assuming a power law for the number of gamma-rays in the energy
interval $(E,\;E+{\textrm{d}}E)$, reads777The normalization factor in front of
the energy power law comes from defining the luminosity in the energy range
$(0.1,\;100)$ GeV.:

$\frac{{\textrm{d}}F}{{\textrm{d}}E}(E,\;L,\;z)=\frac{L}{4\pi
d_{L}^{2}(z)}(1+z)(2-\Gamma)\left[100^{2-\Gamma}-0.1^{2-\Gamma}\right]^{-1}\left(\frac{E}{\text{GeV}}\right)^{-\Gamma}\;\text{GeV}^{-2}\;,$
(5.6)

where $\Gamma$ is the spectral index of the considered source class (see Table
1). The minimum and maximum luminosities in Eq. (5.5) depend on the intrinsic
properties of the source class and are shown in Table 1. However, since we are
dealing with the unresolved sky, $L_{\rm max}$ shall never be larger than the
luminosity corresponding to the sensitivity of the detector. In this work we
consider the _Fermi_ -LAT telescope for gamma-ray surveys and assume a
detector flux sensitivity $F_{\rm
sens}=10^{-10}\;\text{cm}^{-2}\text{s}^{-1}$, for photons in the energy band
$1-100$ GeV, which is well compatible with 8 years of data taking888The
constraint on $L_{\rm max}$ is more relevant at low redshift, where $L_{\rm
sens}<L_{\rm max}$ and the window function tends to fall rapidly. and a
slightly better sensitivity for a future improved detector, discussed below.
Note that Eqs. (5.3), (5.4) and (5.5) must be integrated in the considered
energy band, before being plugged into the CAPS of Eq. (2.2).

In Fig. 5 we show the normalized (to the average intensity $\langle
I\rangle=\int{\textrm{d}}\chi\;W(\chi)$) window functions for gamma-rays with
$E_{\gamma}=5$ GeV. We compare the gravitational tracers (galaxy distribution
from the 2MASS catalog and the cosmic shear) with the gamma-ray emitters:
annihilating/decaying dark matter in the left panel and astrophysical sources
in the right panel. We see how the latter peaks around redshifts $0.5-1$,
while the former, being completely unresolved, peaks at very low redshifts and
then rapidly decays.

##  6 Results

Bin | $E_{\text{min}}\;[\text{GeV}]$ | $E_{\text{max}}\;[\text{GeV}]$ | $N^{\gamma}\;[\text{cm}^{-4}\text{s}^{-2}\text{sr}^{-1}]$ | $f_{\text{sky}}$ | $\sigma_{0}^{\text{Fermi}}\;[\text{deg}]$ | $E_{b}\;[\text{GeV}]$
---|---|---|---|---|---|---
1 | $0.5$ | $1.0$ | $1.056\times 10^{-17}$ | $0.134$ | $0.87$ | $0.71$
2 | $1.0$ | $1.7$ | $3.548\times 10^{-18}$ | $0.184$ | $0.50$ | $1.30$
3 | $1.7$ | $2.8$ | $1.375\times 10^{-18}$ | $0.398$ | $0.33$ | $2.18$
4 | $2.8$ | $4.8$ | $8.324\times 10^{-19}$ | $0.482$ | $0.22$ | $3.67$
5 | $4.8$ | $8.3$ | $3.904\times 10^{-19}$ | $0.549$ | $0.15$ | $6.31$
6 | $8.3$ | $14.5$ | $1.768\times 10^{-19}$ | $0.574$ | $0.11$ | $11.0$
7 | $14.5$ | $22.9$ | $6.899\times 10^{-20}$ | $0.574$ | $0.09$ | $18.2$
8 | $22.9$ | $39.8$ | $3.895\times 10^{-20}$ | $0.574$ | $0.07$ | $30.2$
9 | $39.8$ | $69.2$ | $1.576\times 10^{-20}$ | $0.574$ | $0.07$ | $52.5$
10 | $69.2$ | $120.2$ | $6.205\times 10^{-21}$ | $0.574$ | $0.06$ | $91.2$
11 | $120.2$ | $331.1$ | $3.287\times 10^{-21}$ | $0.597$ | $0.06$ | $199.5$
12 | $331.1$ | $1000$ | $5.094\times 10^{-22}$ | $0.597$ | $0.06$ | $575.4$
Table 2: Gamma-ray energy bins used in this analysis, adherent with 8 years of
data taking from _Fermi_ -LAT Pass 8 (see [108]). $N^{\gamma}$ is the auto-
correlation noise, $f_{\text{sky}}$ the observed fraction of the sky outside
the combined Galactic and point-source masks and $\sigma_{0}^{\text{Fermi}}$
the $68\%$ containment angle of the PSF, referred to the geometric center of
each energy bin $E_{b}=\sqrt{E_{\text{min}}\,E_{\text{max}}}$.

As previously discussed, we underline the potential of cross-correlating
gamma-ray emission from particle dark matter (either through annihilations or
decays) with the most underdense regions of the Universe. Previous works have
proposed [3, 1, 29] and used [4, 13, 20] the cross-correlation formalism
between dark matter emission and gravitational tracers within dark matter
halos, and have shown how the gamma-ray signal generated by unresolved
astrophysical sources typically surpasses that of dark matter, making the
signal-to-noise ratio somehow unfavourable. As discussed in Sec. 4, the
astrophysical signal depends on the GLF which can be more than three orders of
magnitude lower in voids than in halos, assuming the number of sources to be
constant for halos and voids with the same mass and no relevant changes in
their general properties. For this reason, we show here that, while the dark
matter signal from voids is reduced in size as compared to halos, nevertheless
in comparison to its astrophysical background the cross-correlation signal can
be favored with respect to astrophysical sources, ticking the potential
usefulness of exploiting voids to research on the dark sector.

On the observational side, the most efficient gamma-ray detector today is the
_Fermi_ -LAT telescope, which has contributed strongly in the advancement of
our knowledge on extremely powerful events in the Universe. With its excellent
angular and energy resolutions and more than 10 years of service, it has been
used to determine the composition of the UGRB [108, 109, 110, 6]. In this work
we use, as a point of reference, the specifications adopted for the analysis
of the UGRB performed in [108], based on 8 years of data taking and a
selection of events with optimal angular resolution and background rejection.
Ref. [108] is currently the most up-to-date analysis of the statistical
fluctuations of the UGRB, and we therefore adopt the specifications
(sensitivity, energy binning, angular resolution, noise) as a reference for
this analysis. We will then forecast the reach of a future gamma-ray detector
(called Fermissimo for definiteness and in continuity with previous analyses
[2, 29]), as specified in [29] and which is modeled on improved specifications
(like e.g. [111]).

We compute the cross-correlations on the 12 energy bins of [108] and reported
in Table 2. For each energy bin we display the measured photon noise
$N^{\gamma}$, the observed portion of the sky $f_{\text{sky}}$ and the $68\%$
containment angle of the _Fermi_ -LAT point spread function (PSF) around the
geometric center of the bin. These parameters are necessary to calculate the
variance on the CAPS (see Eq. (2.3)). The beam function of the telescope
depends on the photon event class and the energy spectrum and is available
through the _Fermi_ tools. An overall good analytical approximation is given
by (see [29]):

$B_{\ell}^{\gamma}=\text{exp}\left[-\frac{\sigma_{b}^{2}(\ell,\;E)\ell^{2}}{2}\right]\;,$
(6.1)

where the dispersion angle of each energy bin evolves as (see [29]):

$\sigma_{b}(\ell,\;E)=\sigma_{0}^{\text{Fermi}}(E)\left[1+0.25\,\sigma_{0}^{\text{Fermi}}(E)\ell\right]^{-1}\;,$
(6.2)

for which $\sigma_{0}^{\text{Fermi}}(E)$ is the $68\%$ containment angle of
the _Fermi_ -LAT PSF and can be defined in terms of a referenced value
$\sigma_{0}^{\text{Fermi}}(E_{\rm ref}=0.5\;\text{GeV})=1.2$ deg (see [29]):

$\sigma_{0}^{\text{Fermi}}(E)=\sigma_{0}^{\text{Fermi}}(E_{\rm
ref})\left[\frac{E}{E_{\rm ref}}\right]^{-0.95}+0.05\;\text{deg}\;.$ (6.3)

Figure 6:  (_Top_): CAPS of decaying (solid)/annihilating (dashed) dark matter
and cosmic shear both in halos (green) and voids (orange). (_Bottom_):
Relative difference between the cross-correlation signal in halos and voids,
for the two production cases. Note that in the case of voids, the annihilating
signal is largely disfavored, being its relative difference much higher than
for decaying dark matter. The dark matter particle is $m_{\text{DM}}=100$ GeV,
the production channel is $b\bar{b}$, the annihilation rate is $\langle\sigma
v\rangle=3\times 10^{-26}\;\text{cm}^{3}\text{s}^{-1}$, the decay lifetime is
$\tau_{d}=3\times 10^{27}\;\text{s}$.

On the cosmic shear side the currently operating survey is DES [51], which
provides weak lensing shape catalogs with more than $10^{8}$ galaxies [112,
113]. However, Euclid [52] will be the next generation galaxy survey and its
forecasts [103] are extremely promising, especially for cosmic shear as the
telescope is expected to be very sensitive to this effect. Galaxy surveys have
achieved an incredibly high angular resolution and their beam function in
harmonic space $B_{\ell}$ can be considered equal to unity in the multipole
range of interest here, while the noise associated to the auto-correlation is,
for Euclid: $N^{\text{Euclid}}=\sigma^{2}_{\epsilon}/\overline{\rm N}_{g}$,
where $\sigma_{\epsilon}=0.3$ is the intrinsic ellipticity and $\overline{\rm
N}_{g}=30\;\text{arcmin}^{-2}$ is the average number of galaxies per steradian
(see e.g [2]). Similarly, for galaxies the beam function can be set to unity
and the noise simply reduces to the average number of galaxies.

In the gaussian error estimates on the signals of Eq. (2.3), the cross-
correlation term is largely sub-dominant, while the gamma-ray auto-correlation
term is dominated by the noise999This is only true when considering the
specifications of _Fermi_ -LAT, while when considering an improved future
detector, namely Fermissimo, the noise and beam function are largely reduced
and their contribution to the error is comparable with that of the gamma-ray
auto-correlation, and therefore we include them in the error estimate. (see
e.g. [29]). The variance of the cross-correlation between gamma-rays and
gravitational tracers (namely "$t$") can therefore be safely approximated by:

$(\Delta
C^{t\gamma}_{\ell})^{2}\simeq\frac{1}{(2\ell+1)f_{\text{sky}}}\left[\frac{N^{\gamma}}{(B_{\ell}^{\gamma})^{2}}\left(C^{tt}_{\ell}+N^{t}\right)\right]\;.$
(6.4)

For definiteness, as a representative case of WIMP, in the following we will
consider a dark matter particle of 100 GeV mass and annihilating or decaying
into the $b\bar{b}$ channel with a canonical-thermal-relic annihilation rate
of $\langle\sigma_{a}v\rangle=3\times 10^{-26}\;\text{cm}^{3}\,\text{s}^{-1}$
or a decay lifetime of $\tau_{d}=3\times 10^{27}\;\text{s}$ (close to its
conservative lower bound [104, 105]), respectively.

Let us turn first to the discussion of the cross-correlation signal between
dark matter gamma-ray emission and gravitational tracers. In particular, we
first compare the decaying signal with the annihilating one (we chose to use
only the cosmic shear as a reference for gravitational tracers, for
simplicity). The CAPS of Eq. (2.2) depends on the spectra and window functions
for dark matter and cosmic shear, whose prescriptions are shown in Secs. 4 and
5, respectively. In the top panel of Fig. 6 we show the CAPS of dark matter
decay (solid lines) and annihilation (dashed lines), separately for halos
(green) and voids (orange). The signal produced by decaying dark matter is
higher than that of annihilating dark matter across the whole range of
multipoles and for both types of structures. While both annihilating and
decaying signals can be modified by changing $\langle\sigma v\rangle$ or
$\tau$ (which we have set here close to their bounds), nonetheless, as a
general property, the relative difference between the signal in voids and
halos is highly disfavored for annihilations. In fact, as shown by the bottom
panel of Fig. 6, in the case of decay the signal in voids is about a factor of
50 smaller than the corresponding signal in halos, while for annihilations the
ratio is significantly depressed. This behavior is intuitively connected to
the dependence of the source function on the square of the density, enhancing
the already conspicuous ratio between the density of halos and voids. We
conclude that underdense regions are not an efficient probe of annihilating
dark matter and we will focus only on decaying dark matter for the following
discussion.

Remaining in the realm of decaying dark matter, the relevant angular power
spectra are shown in Fig. 7 for the cross-correlation between cosmic shear and
gamma-ray emitters, and in Fig. 8 for that between the galaxy distribution and
gamma-ray emitters. In each figure, we display the contribution from the four
classes of unresolved astrophysical sources (BL Lac, FSRQ, mAGN and SFG) and
from decaying dark matter, both for halos (in the left panel) and voids (in
the right panel). The results refer to the sum of contributions from all the
energy bins of Table 2. Among astrophysical sources, for the cosmic shear
case, BL Lacs and SFGs provide for the dominant terms, while FSRQs are
strongly sub-leading. In the case of galaxies, however, BL Lacs give way to
mAGN, which climb up to the dominant contribution, together with SFGs, for
high multipoles. This derives from a complex interplay between the redshift
and energy dependence of the CAPS. Intuitively, mAGN have a greater impact
when cross-correlating with galaxies as the 2MASS window function is peaked at
lower redshift (where blazars tend to have extremely low power, due to the
detector sensitivity cutoff) than the cosmic shear one, which, on the
contrary, peaks right around the maximum of the window function for BL Lacs,
that is strongly suppressed at low redshift (these behaviors are trivial by
looking at the right panel of Fig. 5).

Figure 7:  Cross-correlation angular power spectrum between cosmic shear and
the unresolved gamma-ray sky. The total signal (solid black) is given by the
sum of contributions from all gamma-ray sources: astrophysical sources and
dark matter decays (we consider a DM mass $m_{\text{DM}}=100\;\text{GeV}$ and
decay lifetime $\tau_{d}=3\times 10^{27}\;\text{s}$). In the left panel we
show the signals computed within halos and in the right panel those within
voids. The result refers to the sum of contributions from all energy bins of
Table 2 and the error bars are obtained from the Gaussian estimate of the
variance of the signals. Error bars are calculated for the Fermi-LAT
configuration discussed in the text. Figure 8:  The same as Fig. 7 for
galaxies $\times$ gamma-rays.

The main observation relative to Fig. 7 and 8 is that in the case of halos,
the dark matter signal is either completely dominated by astrophysical sources
or too close to their signal to be effectively distinguished by the present-
day detectors. However, as shown by the right panel of Figs. 7 and 8, in voids
decaying DM provides the dominant contribution, even almost completely setting
the total signal, dominating astrophysical sources by more than three orders
of magnitude for the cosmic shear and from three to one order of magnitude
(for increasing multipoles) for galaxies. This behavior is due to the expected
much lower density of astrophysical sources in underdense environments, as
discussed in Sec. 4 when examining the difference between the GLF in halos and
voids.

In summary, though the UGRB signal is expected to be much lower in voids than
in halos, the signal from decaying dark matter is expected to dominate over
the astrophysical background by a sizeable amount, contrary to the case of
halos, making this option (if accessible, given the sensitivity of gamma-ray
detectors) a background-free signal for dark matter lifetimes 3-4 orders of
magnitude smaller than current bounds.

Figure 9:  Signal-to-background ratio (at $\ell=100$) as a function of the
dark matter decay lifetime. In the left panel we show the result for the
cross-correlation with cosmic shear, in the right panel that with galaxies. In
both cases we report values for four different values of the WIMP mass. The
signal-to-background ratio is shown in black for the all-sky cross-correlation
(given by the sum of contributions from halos and voids) and in oranges for
the correlation signal computed within cosmic voids only.

This can be seen by looking at Fig. 9, where we report the ratio between the
dark matter signal and the astrophysical background as a function of the decay
lifetime and for different WIMP masses. In the left panel we show the result
for cosmic shear, in the right panel that of galaxies. The grey lines refer to
the total contribution from both halos and voids and lead to a signal-to-
background ratio (even significantly) lower than unity, if not for small
values of the decay lifetime (basically already outside current bounds). On
the contrary, the signal-to-background ratio within voids (orange lines) is
way larger than unity also for very large values of the decay lifetime. This
"background-free" situation for a relatively large fraction of the dark matter
particle parameter space is analogous to the case of the galactic antideuteron
signal, where a large signal-to-background ratio is present at low $\bar{D}$
energies [114]. These two channel ($\bar{D}$ and voids cross-correlation)
share the feature of being potentially offering promising opportunities,
though requiring high sensitivities [115].

The two gravitational tracers used here provide similar results. However,
Figs. 7 and 8 are helpful to determine which of them might be a better probe
for decaying dark matter. In particular, the relative difference between the
dark matter signal and that of astrophysical sources in voids, as discussed
above, is on average higher for the cross-correlation with the cosmic shear,
while being greater for galaxies for central multipoles as shown in Fig. 9. At
the same time, the difference between the total signal in halos and in voids
is lower for galaxies, especially at low multipoles where is set around one
order of magnitude, compared to the two orders of magnitude for the cosmic
shear.

Figure 10:  The same as Fig. 7 using the Fermissimo configuration to compute
error bars. Figure 11:  The same as Fig. 8 using the Fermissimo configuration
to compute error bars.

Concerning detectability of a signal, we need to consider that the errors on
the signals computed through Eq. (2.3) with the Fermi-LAT configuration are
large at single multipoles, as well for the individual energy bins of Figs. 7
and 8. However, we can exploit the full range of multipoles accessible with
the Fermi-LAT PSF, i.e. ranges between about $\ell\simeq 50$ (below which the
large-scale galactic foreground is a limiting factor) and $\ell\simeq$
(200$\div 1000)$ depending on the energy bin [108]. This can be traced to the
strong domination of the gamma-ray noise/beam term from _Fermi_ -LAT, as
mentioned above (see Eq. (6.4)).

In order to investigate the potentiality of the cross-correlation signals to
probe dark matter, we consider a future gamma-ray detector with improved
specifications, following [29]. First, we assume the exposure of the detector
to be larger by a factor of 2 as compared with the current _Fermi_ -LAT
specification here considered. Second, we assume that the detector PSF can be
improved, adopting the same behaviour of the beam function expressed in Eq.
(6.1) but with a better angular resolution (see [29]):

$\sigma_{0}(E)=\alpha_{\sigma}\times\sigma_{0}^{\text{Fermi}}(E_{\rm
ref})\left[\frac{E}{E_{\rm ref}}\right]^{-0.95}+0.001\;\text{deg}\;,$ (6.5)

for which we assume $\alpha_{\sigma}=0.2$, for definiteness. Finally, thanks
to the better angular resolution, which means smaller mask, we adopt a
(somewhat optimistic) larger sky-fraction coverage of $f_{\text{sky}}=0.8$,
allowing to slightly reduce the impact of noise. Using the same energy bins of
Table 2, we scale the noise as the inverse of the exposure, as in Ref. [29].
We refer to this configuration as _Fermissimo_ and show its effect of reducing
the variance on the cross-correlation signal in Figs. 10 and 11 (for the
cosmic shear and the 2MASS galaxy distribution, respectively), when compared
to the _Fermi_ results of Figs. 7 and 8. With the statistical technique
discussed below, we find that the fiducial model we are adopting in the
analysis ($m_{\textrm{DM}}=100$ GeV and $\tau_{d}=3\times 10^{27}$ s) would
lead to a statistical significance of $5.7\sigma$ in the case of galaxies
(being $1.6\sigma$ in the case of cosmic shear). The corresponding numbers for
the Fermi-LAT configuration are $1.3\sigma$ and $0.6\sigma$, respectively.

Let us in fact now quantify the reach in terms of detectability of a signal or
setting bounds on the DM properties that the cross-correlation technique with
voids can lead to. To this aim we adopt the Fisher matrix formalism [116, 117,
118], which is a likelihood-based statistics suitable to forecast the
capabilities of future experiments. Let us therefore consider our observable,
namely the $C_{\ell}$’s, modeled through a defined set of free parameters
$\theta_{a\,=\,1,\,\dots\,,\,n}$ and with an associated covariance
$\Gamma_{\ell\ell^{\prime}}$; we can write the Fisher matrix (a representation
of the covariance matrix for the model parameters, associated to a maximum-
likelihood estimate) as:

$F_{ab}=\sum_{\ell\ell^{\prime}}\frac{\partial
C_{\ell}}{\partial\theta_{a}}\,\Gamma_{\ell\ell^{\prime}}^{-1}\,\frac{\partial
C_{\ell^{\prime}}}{\partial\theta_{b}}\;.$ (6.6)

For any free parameter $\theta_{a}$, its error can be estimated as
$\sigma(\theta_{a})=\sqrt{(F^{-1})_{aa}}$ and therefore we can set a bound at
$n$-$\sigma$ significance level through:

$\theta_{a}^{\text{bound}}=n\times\sigma(\theta_{a})\;.$ (6.7)

In our case, we model the total CAPS at fixed WIMP mass with seven free
parameters:

$\begin{split}C_{\ell}&=p\,\tilde{C}_{\ell}^{\text{DM}}+A_{\text{1V}}\,\tilde{C}_{\ell,\,\text{1V}}^{\text{astro}}(\Gamma_{i})+A_{\text{2V}}\,\tilde{C}_{\ell,\,\text{2V}}^{\text{astro}}(\Gamma_{i})\;,\\\
\tilde{C}_{\ell,\,\text{1V/2V}}^{\text{astro}}(\Gamma_{i})&=\sum_{\text{i=1}}^{4}\tilde{C}_{\ell,\,\text{1V/2V}}^{i}(\Gamma_{i})\;,\end{split}$
(6.8)

where $p=(3\times 10^{27}\;\text{s}/\tau_{d})$ is the decay lifetime-
normalization with respect to the fiducial cross-correlation for DM computed
with $\tau_{d}=3\times 10^{27}$ s and $A_{\text{1V}}$ and $A_{\text{2V}}$ are
free normalizations for the 1V and 2V contributions to the total astrophysical
signal, computed with our fiducial model (see Appendix A). The latter are
given by the sum of the contributions from the four unresolved astrophysical
sources (BL Lac, FSRQ, mAGN and SFG), each computed with the fiducial value
for the related spectral index $\Gamma_{i}$ from Table 1. The two
normalizations account for variations in shape related to the multipoles,
while the spectral indexes carry all of the information about the photon’s
energy dependence of the spectra (see Sec. 5). With our gaussian theoretical
assumptions, the covariance matrix can be assumed to be diagonal:
$\Gamma_{\ell\ell^{\prime}}(E)=\delta_{\ell\ell^{\prime}}(\Delta
C_{\ell}(E))^{2}=\delta_{\ell\ell^{\prime}}\sigma^{2}$. Therefore the Fisher
matrix takes the form:

$F=\begin{pmatrix}\sum_{\ell,E}\frac{(\tilde{C}_{\ell}^{\text{DM}})^{2}}{\sigma^{2}}&\sum_{\ell,E}\frac{\tilde{C}_{\ell}^{\text{DM}}\tilde{C}_{\ell,\,\text{1V}}^{\text{astro}}}{\sigma^{2}}&\cdots\\\\[8.53581pt]
\sum_{\ell,E}\frac{\tilde{C}_{\ell}^{\text{DM}}\tilde{C}_{\ell,\,\text{1V}}^{\text{astro}}}{\sigma^{2}}&\sum_{\ell,E}\frac{(\tilde{C}_{\ell,\,\text{1V}}^{\text{astro}})^{2}}{\sigma^{2}}&\cdots\\\\[8.53581pt]
\vdots&\vdots&\ddots\\\\[8.53581pt] \end{pmatrix}\;,$ (6.9)

where the sums are over all considered multipoles and the energy bins of Table
2. It is a 7x7 symmetric square matrix since we are considering seven
parameters: the DM normalization $p$, the two astrophysical normalizations
$A_{\text{1V}}$ and $A_{\text{2V}}$ and the four spectral indexes
$\Gamma_{\text{BL}}$, $\Gamma_{\text{FSRQ}}$, $\Gamma_{\text{mAGN}}$ and
$\Gamma_{\text{SFG}}$. From Eq. (6.7) we can then compute the bound on $p$ and
consequently on the decay lifetime $\tau_{d}$, with respect to varying values
of the WIMP mass. The Fisher formalism allows to efficiently consider all
uncertainties on the chosen models and, in particular, to consider the impact
of the astrophysical background noise on the detection of the dark matter
signal.

Figure 12:  Forecast of the bounds on the DM properties (decay lifetime
$\tau_{d}$ vs. WIMP mass $m_{DM}$) attainable through the study of the cross-
correlation between the unresolved gamma-ray sky and either cosmic shear
(green) or the galaxy distribution (red) within cosmic voids. We show both the
$2\sigma$ bound and the $5\sigma$ detection reach, considering a highly
efficient gamma-ray detector, namely Fermissimo. We also report previously
derived bounds from cosmic-ray antiprotons [104] and from the isotropic gamma-
ray background [105]. (_Left_): conservative bound. (_Right_): the bound
accounts for priors on the non-DM parameters; $\pm\;60\%$ for astrophysical
normalizations and $\pm\;0.4$ for spectral indexes.

In Fig. 12 we show the $2\sigma$ (95% C.L. bound reach) and $5\sigma$
(detection limit) reach for the Fermissimo configuration, attainable with
voids (cosmic shear in green and galaxies in red). We compare these forecasts
with bounds obtained with cosmic antiprotons [104] and the isotropic gamma-ray
background (IGRB) [105]. In the left panel we report the bound for the
marginalized case (conservative). In the right panel we show the bounds
obtained by accounting for priors on the non-DM parameters (determined from
the measured uncertainties on the fiducial model parameters): for the
astrophysical normalizations $A_{\text{1V}}$ and $A_{\text{2V}}$ we consider
variations of $\pm\;60\%$, while for spectral indexes we take a prior of
$\pm\,0.4$ on the fiducial values of $\Gamma_{i}$. In general, the galaxies
would allow to reach tighter bounds than cosmic shear. The forecasted bounds
are competitive with previous works only for galaxies in the conservative
scenario and for both galaxies and shear in the non-conservative scenario,
when accounting for priors on the parameters. In particular, the current
bounds could be improved (in the case of galaxies) for $m_{\textrm{DM}}\sim
25\div 900$ GeV in the conservative case (and in the whole range of considered
masses for the non-conservative case), while being generally worse
(conservative) or slightly worse (non-conservative) for cosmic shear. Let us
notice that the drop in sensitivity toward lower masses is due to the fact
that we are considering photon energies starting from 500 MeV (to conform to
the specifications of Ref. [108], from which we adopt relevant information,
like the photon noise, see Table 2): by decreasing the lower energy around 100
MeV would allow us to improve the forecasted bounds for lighter dark matter,
making the lines flatter, since an energy threshold at 500 MeV cuts out a
relevant fraction of the signal for dark matter masses close to 10 GeV,
especially for the $b\bar{b}$ decay channel we are considering here. The
detector PSF at those energies would be worse, though.

In conclusion, although voids lead to a much weaker cross-correlation signal
than halos, their signal-to-background ratio is much larger than 1, for decay
lifetimes up to $2\times 10^{30}$ s. This represents a situation in which a
detection would point toward a dark matter interpretation, contrary to the
case of a detection in halos, for which the signal-to-background ratio is
unfavourable. The signal in voids is nevertheless difficult to achieve and
requires improved sensitivities of gamma-ray detectors: with improved
specifications like those discussed in the text, the reach on bounds on the
dark matter lifetime can exceed current bounds obtained through different
techniques in a mass range between 25 GeV and 900 GeV.

##  7 Conclusion

In this paper we have discussed the idea of using cosmic voids as probes for
particle dark matter. We considered the cosmic shear and the galaxy
distribution as gravitational tracers of the matter distribution in the
Universe to investigate the nature of the unresolved gamma-ray background,
which can hide a signal due to dark matter annihilation or decay. Large-scale
structures are responsible for the bending of light in the weak lensing regime
as well as for gamma-ray emission either from unresolved astrophysical sources
or particle dark matter annihilations and decays. At the same time, the
distribution of galaxies is tightly connected to the distribution of dark
matter as well as to that of astrophysical sources hosted by them. We thus
expect a positive correlation between these gravitational tracers and the
unresolved gamma-ray emission.

We have therefore computed the cross-correlation angular power spectrum
between the fluctuations due to the inhomogeneous distribution of matter in
the Universe, traced by either the cosmic shear or the galaxy distribution,
and those induced in the unresolved gamma-ray emission, by both astrophysical
sources and dark matter annihilations/decays. We separated the signals into
the contributions coming from voids and halos, respectively, to investigate
the main differences among their properties within the two types of
structures. We show that in voids, the cross-correlation signal generated by
dark matter annihilation is strongly disfavored with respect to dark matter
decay, as it depends on the square of the structure’s density, leading to a
much greater difference between the signal in halos and the signal in voids
than the case of decaying dark matter.

On the other hand, in the case of voids and for decaying dark matter, we find
that the dark matter cross-correlation signal can largely exceed the
astrophysical counterpart, for an interestingly large section of the dark
matter parameter space. For decay lifetimes up to $2\times 10^{30}$ s the
signal-to-background ratio is (even significantly) larger than 1. This is at
variance with the cross-correlation from halos, where the astrophysical signal
dominates over the dark matter signal. This makes the signal from voids a
potentially "background-free" option for dark matter searches for decay
lifetimes about 3 orders of magnitude larger than current bounds.

The size of signal is nevertheless small, which makes observational
opportunities to require next generation detectors. We found in fact that the
combination of forthcoming galaxy surveys, such as Euclid, and the gamma-ray
_Fermi_ -LAT telescope is not able yet to detect the signal. However,
potential detectability can be achieved by considering an improved gamma-ray
detector, with better angular resolution and slightly larger exposure: we
assumed a factor 2 larger exposure as compared to the Fermi-LAT, and an
angular resolution a factor of 5 better than Fermi-LAT. In this case, the
cross-correlation which uses galaxies as a dark matter tracer, would reach a
5.7$\sigma$ stastistical significance for the fiducial dark matter models we
adopted in the analysis ($m_{\text{DM}}\sim 100\;\text{GeV}$ and
$\tau_{d}=3\times 10^{27}\;\text{s}$) and will be able to improve on current
bounds on the dark matter decay lifetime for masses in the range
$m_{\textrm{DM}}\sim 25\div 900$ GeV.

## Acknowledgments

We acknowledge support from: Departments of Excellence grant awarded by the
Italian Ministry of Education, University and Research (MIUR); Research grant
The Dark Universe: A Synergic Multimessenger Approach, Grant No. 2017X7X85K
funded by the Italian Ministry of Education, University and Research (MIUR);
Research grant The Anisotropic Dark Universe, Grant No. CSTO161409, funded by
Compagnia di Sanpaolo and University of Torino; Research grant TAsP
(Theoretical Astroparticle Physics) funded by Istituto Nazionale di Fisica
Nucleare (INFN); Fermi Research Alliance, LLC under Contract No. DE-
AC02-07CH11359 with the U.S. Department of Energy, Office of High Energy
Physics.

##  Appendix A Gamma-ray luminosity functions

The astrophysical sources are characterised by their GLF
$\phi_{\gamma}={\textrm{d}}N_{\gamma}/({\textrm{d}}L\,{\textrm{d}}V)$, which
specifies the number of gamma-ray sources per unit of luminosity $L$ and
comoving volume $V$. In this appendix we provide the GLFs for each
astrophysical class considered in our analysis.

###  A.1 Blazars

Following the LDDE model of Ref. [97], the GLF of BL Lacs and FSRQs can be
parameterized as a broken power-law in luminosity and redshift:

$\displaystyle\phi_{\gamma}\left(L,\,z\right)=\dfrac{A}{\ln(10)\,L}\left[\left(\dfrac{L}{L_{\star}}\right)^{\gamma_{1}}+\left(\dfrac{L}{L_{\star}}\right)^{\gamma_{2}}\right]^{-1}\times\left[\left(\dfrac{1+z}{1+z_{c}(L)}\right)^{-p_{1}}+\left(\dfrac{1+z}{1+z_{c}(L)}\right)^{-p_{2}}\right]^{-1}\,,$
(A.1)

where $z_{c}=z_{\star}\,(L/10^{48}\textrm{erg s}^{-1})^{\beta}$ and the other
parameters are specified in Table 3.

|  $A$ [Mpc-3] |  $L_{\star}$ [erg s-1] | $\gamma_{1}$ | $\gamma_{2}$ | $p_{1}$ | $p_{2}$ | $z_{\star}$ | $\beta$
---|---|---|---|---|---|---|---|---
BL Lacs | $9.20\times 10^{-11}$ | $2.43\times 10^{48}$ | $1.12$ | $3.71$ | $4.50$ | $-12.88$ | $1.67$ | $4.46\times 10^{-2}$
FSRQs | $3.06\times 10^{-9}$ | $0.84\times 10^{48}$ | $0.21$ | $1.58$ | $7.35$ | $-6.51$ | $1.47$ | $0.21$
Table 3: Parameters of the gamma-ray luminosity function for BL Lacs and
FSRQs.

###  A.2 Misaligned active galactic nuclei

$\beta_{l}$ |  $L_{l\star}$ [W/Hz] |  $\rho_{l\star}$ [Mpc-3] | $z_{l\star}$ | $k_{l}$ | $\beta_{h}$ |  $L_{h\star}$ [W/Hz] |  $\rho_{h}$ [Mpc-3] | $z_{h\star}$
---|---|---|---|---|---|---|---|---
$0.586$ | $10^{26.48}$ | $10^{-7.523}$ | $0.71$ | $3.48$ | $2.42$ | $10^{27.39}$ | $10^{-6.757}$ | $2.03$
Table 4: Parameters of the radio luminosity function for mAGN.

The GLF for mAGN can be derived from the radio luminosity function (RLF)
through

$\phi_{\gamma}(L,\,z)=\dfrac{k\,\eta}{(1+z)^{2-\Gamma}}\,\dfrac{1}{\ln(10)\;L^{\textrm{151
MHz}}_{\textrm{tot}}}\dfrac{{\textrm{d}}L^{\textrm{151
MHz}}_{\textrm{tot}}}{{\textrm{d}}L}\;\rho_{r}(L^{\textrm{151
MHz}}_{\textrm{tot}}(L),\,z)\;,$ (A.2)

where $k=3.05$, $\Gamma=2.37$ and

$\eta=\dfrac{{\textrm{d}}^{2}V_{W}/{\textrm{d}}z{\textrm{d}}\Omega}{{\textrm{d}}^{2}V/{\textrm{d}}z{\textrm{d}}\Omega}\;.$
(A.3)

The comoving volume
$\dfrac{{\textrm{d}}^{2}V_{W}}{{\textrm{d}}z\,{\textrm{d}}\Omega}$ used by
[99] and the one in the standard $\Lambda$CDM cosmology
$\dfrac{{\textrm{d}}^{2}V}{{\textrm{d}}z\,{\textrm{d}}\Omega}$ are

$\displaystyle\dfrac{{\textrm{d}}^{2}V_{W}}{{\textrm{d}}z\,{\textrm{d}}\Omega}$
$\displaystyle=$
$\displaystyle\dfrac{c^{3}\,z^{2}\,(2+z)^{2}}{4\,H^{3}_{0,W}\,(1+z)^{3}}\;,$
(A.4)
$\displaystyle\dfrac{{\textrm{d}}^{2}V}{{\textrm{d}}z\,{\textrm{d}}\Omega}$
$\displaystyle=$
$\displaystyle\dfrac{c\,d^{2}_{L}(z)}{H_{0}\,(1+z)^{2}\sqrt{(1-\Omega_{\Lambda}-\Omega_{m})(1+z)^{2}+(1+z)^{3}\Omega_{m}+\Omega_{\Lambda}}}\;.$
(A.5)

with $H_{0,W}=50$ $\textrm{km}\,\textrm{s}^{-1}\,\textrm{Mpc}^{-1}$. The
relation between core radio luminosity and gamma-ray luminosity is provided in
[100], while [119] derived the correlation between core and total
luminosities:

$\displaystyle\log L$ $\displaystyle=$ $\displaystyle 2+1.008\log
L_{r,\textrm{core}}$ (A.6) $\displaystyle\log L^{\textrm{5
GHz}}_{r,\textrm{core}}$ $\displaystyle=$ $\displaystyle 4.2+0.77\log
L^{\textrm{1.4 GHz}}_{r,\textrm{tot}}\;.$ (A.7)

The reference radio frequency in Eq. (A.2) is 151 MHz thus, following [120],
we consider the power-law scaling

$\dfrac{L_{r}}{\nu}\propto v^{-\alpha_{r}}$ (A.8)

with $\alpha_{r}=0.80$. The radio luminosity function [99] can be expressed as
the sum of two components

$\rho_{r}(L_{r},z)=\rho_{l}(L_{r},z)+\rho_{h}(L_{r},z)\;,$ (A.9)

where

$\displaystyle\begin{cases}\rho_{l}=\rho_{l\star}\left(\dfrac{L_{r}}{L_{l\star}}\right)^{-\beta_{l}}\,\exp\left(-\dfrac{L_{r}}{L_{l\star}}\right)\,\left(1+z\right)^{k_{l}}\quad\quad\;\;&\textrm{for}\quad
z<z_{l\star}\\\\[14.22636pt]
\rho_{l}=\rho_{l\star}\left(\dfrac{L_{r}}{L_{l\star}}\right)^{-\beta_{l}}\,\exp\left(-\dfrac{L_{r}}{L_{l\star}}\right)\,\left(1+z_{l\star}\right)^{k_{l}}\quad&\textrm{for}\quad
z\geq z_{l\star}\end{cases}$ (A.10)

and

$\rho_{h}=\rho_{h\star}\left(\dfrac{L_{r}}{L_{h\star}}\right)^{-\beta_{h}}\exp\left(-\dfrac{L_{h\star}}{L}\right)\,\exp\left\\{-\dfrac{1}{2}\left(\dfrac{z-z_{h\star}}{z_{h0}}\right)^{2}\right\\}\,.$
(A.11)

For $z<z_{h\star}$ we adopted $z_{h0}=0.568$, while for $z\geq z_{h\star}$ we
used $z_{h0}=0.956$. All the parameters included in Eqs.(A.10) and (A.11) are
specified in Table 4.

###  A.3 Star-forming galaxies

From the infrared luminosity function (ILF)
$\phi_{\textrm{IR}}=\dfrac{{\textrm{d}}N_{\textrm{IR}}}{{\textrm{d}}\log_{10}(L_{\textrm{IR}})\,{\textrm{d}}V}$
we can obtain the GLF of SFGs:

$\phi_{\gamma}(L,\,z)=\phi_{\textrm{IR}}\,\dfrac{{\textrm{d}}\log_{10}(L_{8-1000\,\mu\textrm{m}})}{{\textrm{d}}L_{0.1-100\,\textrm{GeV}}}\,.$
(A.12)

The luminosity $L_{0.1-100\,\textrm{GeV}}$ between 0.1 GeV and 100 GeV and the
luminosity $L_{8-1000\,\mu\textrm{m}}$ between 8 $\mu$m and 1000 $\mu$m are
related via [121]

$\log_{10}\,\left(\dfrac{L_{0.1-100\textrm{GeV}}}{\textrm{erg
s}^{-1}}\right)=\alpha_{\textrm{IR}}\;\log_{10}\,\left(\dfrac{L_{8-1000\,\mu\textrm{m}}}{10^{10}\,L_{\odot}}\right)+\beta_{\textrm{IR}}$
(A.13)

with coefficients $\alpha_{\textrm{IR}}=1.09$ and $\beta_{\textrm{IR}}=39.19$.
Following [101], the ILF can be written as the sum of quiescent spiral
galaxies, starburst galaxies and SFG hosting a concealed or low-luminosity
AGN:

$\phi_{\textrm{IR}}=\phi_{\textrm{spiral}}+\phi_{\textrm{starburst}}+\phi_{\textrm{SF-
AGN}}\;.$ (A.14)

The ILF of each sub-class can be modelled as

$\phi_{i}=\phi_{0,\,i}(z)\left(\dfrac{L_{8-1000\,\mu\textrm{m}}}{L_{0,\,i}}\right)^{1-\gamma_{i}}\,\exp\left[-\dfrac{1}{2\sigma^{2}_{i}}\log^{2}_{10}\left(1+\dfrac{L_{8-1000\,\mu\textrm{m}}}{L_{0,\,i}}\right)\right]\;,$
(A.15)

where $i$ = {spiral, starburst, SF-AGN}. The radio luminosity can be written
as a function of the gamma-ray luminosity using Eq. (A.13) and the
normalization $\phi_{0,i}$ reads

$\phi_{0,i}=\begin{cases}\phi_{\star,\,{i}}\left(\dfrac{1+z}{1.15}\right)^{k_{R1,\,i}}&\textrm{for}\;z\leq
z_{\star,i}\\\\[8.53581pt]
\phi_{\star,\,{j}}\left(\dfrac{1+z_{\star,i}}{1.15}\right)^{k_{R1,\,i}}\,\left(\dfrac{1+z}{1+z_{\star,i}}\right)^{k_{R2},\,\textrm{i}}&\textrm{for}\;z>z_{\star,i}\,.\end{cases}$
(A.16)

All the parameters in Eqs. (A.15) and (A.16) are specified in Table 5.

| $\gamma$ | $\sigma$ | $\log_{10}(L_{\star}/L_{\odot})$ | $\log_{10}(\phi_{\star}/\textrm{Mpc}^{-3})$ | $k_{L}$ | $k_{R1}$ | $k_{R2}$ | $z_{\star}$
---|---|---|---|---|---|---|---|---
spiral | $1.0$ | $0.50$ | $9.78$ | $-2.12$ | $4.49$ | $-0.54$ | $-7.13$ | $0.53$
starburst | $1.0$ | $0.35$ | $11.17$ | $-4.46$ | $1.96$ | $\,3.79$ | $-1.06$ | 1.1
SF-AGN | $1.2$ | $0.40$ | $10.80$ | $-3.20$ | $3.17$ | $\,0.67$ | $3.17$ | $1.1$
Table 5: Parameters of the infrared luminosity function for the three sub-
classes of SFGs. Figure 13:  3D power spectrum of the cross-correlation
between cosmic shear and gamma-ray emitters (decaying/annihilating DM and
unresolved astrophysical sources), at $z=0.5$, divided into the contributions
from a single halo (1H), two halos (2H), a single void (1V), two voids (2V)
and an halo and a void (HV). Figure 14:  The same as Fig. 13 for the cross-
correlation of gamma-ray emitters with the 2MASS galaxy catalog.

##  Appendix B Mass-to-luminosity relations

The astrophysical sources are better characterised by their luminosity, as
opposed to the DM halos which are better characterised by the halo mass.
Therefore, the power spectra involving the astrophysical components include an
integral over the luminosity and the relation $M(L)$ between the luminosity of
the astrophysical class and the host DM halo is required. Following [2], the
mass-to-luminosity relations for blazars (BL Lacs and FSRQs), mAGN and SFGs
read

$\displaystyle M_{i}(L)$ $\displaystyle=$ $\displaystyle
10^{13}\,M_{\odot}\;\left(\dfrac{M_{\star,i}}{10^{8.8}\,(1+z)^{1.4}}\right)^{0.645}\quad
i=\textrm{blazar, mAGN}$ (B.1) $\displaystyle M_{\textrm{SFG}}(L)$
$\displaystyle=$
$\displaystyle\dfrac{10^{12}\,M_{\odot}}{(1+z)^{1.61}}\,\left(\dfrac{L}{6.8\cdot
10^{39}\textrm{ erg/s}}\right)^{0.92}\;,$ (B.2)

with

$\displaystyle M_{\star,\textrm{blazar}}$ $\displaystyle=$ $\displaystyle
10^{9}\,\left(\dfrac{L}{10^{48}\textrm{ erg/s}}\right)^{0.36}\;$ (B.3)
$\displaystyle M_{\star,\textrm{mAGN}}$ $\displaystyle=$ $\displaystyle
4.6\cdot 10^{9}\,\left(\dfrac{L}{10^{48}\textrm{ erg/s}}\right)^{0.16}\;.$
(B.4)

##  Appendix C Terms of the 3D power spectrum and variations on the void
profile

When computing the 3D power spectrum within the HVM, three additional terms,
with respect to the HM, arise: the correlation between points within the same
void (1Void term), within two different voids (2Void term) and between points
belonging to a halo and a void, respectively (Halo-Void term). In Figs. 13 and
14 we show the 3D power spectrum of the cross-correlation between gamma-ray
emitters (decaying/annihilating DM, BL Lacs, FSRQs, mAGN and SFGs) and either
cosmic shear or the 2MASS galaxy catalog, respectively (calculated at
$z=0.5$). For each panel of the figures, we report the contributions from the
terms of Eqs. (4.1)-(4.5). In general, the linear power on large scales is
well recovered and the transition between linear and non-linear scales (around
$k\sim 1\;h\,\text{Mpc}^{-1}$) is clearly visible in all cases.

Note how the 1Halo and 1Void terms in the cross-correlation with galaxies
(Fig. 14) act as a shot noise, since the integrand of Eqs. (4.1) and (4.3) are
fairly constant with the scale $k$, for the relevant source functions.

The 3D power spectrum depends on the considered structure’s profile for all
cross-correlations used in this work (see Sec. 4). However, while the halo
density profile has been studied for a long time and the NFW profile is
typically considered the standard choice, the void density profile is debated
in the literature and strongly depends on the used void finder (as discussed
in Sec. 3.1). Here, we considered both voids with a radius-dependent central
density, showing compensation walls, and voids with empty centers and no
compensation walls. The former are well reproduced by the Hamaus-Sutter-
Wandelt profile [74] of Eq. (3.29) and the latter by the profile proposed in
Voivodic et al. [33], shown in Eq. (3.30). Despite having different
properties, Fig. 15 shows how the choice of the void profile leads to sub-
percent variations of the 3D power spectrum, for all cross-correlations.

Figure 15:  Relative variations on the 3D power spectrum upon using the HSW
profile [74] or the tanh-profile [33]. We show them for the cross-correlation
of the 2MASS galaxy catalog (green) or the cosmic shear (blue) with the
decaying dark matter signal (solid) and the total unresolved astrophysical
signal (dashed). Differences are confined well below the few per mille level.

##  Appendix D DM and non-DM parameters’ correlation in the Fisher analysis

Figure 16:  Bivariate correlations between the DM free parameter of the Fisher
matrix analysis $p$ and the six non-DM parameters. All are referred to a WIMP
with $m_{\text{DM}}=100$ GeV.

In Sec. 6 we showed the forecasted bound on the DM decay lifetime, computed
through a Fisher matrix analysis with seven parameters: one normalization for
the DM signal, two normalizations for the astrophysical signal and the four
spectral indexes of unresolved astrophysical sources, accounting for energy-
dependent variations of the $C_{\ell}$’s. The Fisher matrix in Eq. (6.9) can
be seen as the inverse covariance of the model’s free parameters. Thus, it can
be used to perform a multivariate study on the degeneration of the parameter.
In particular, we can consider the six submatrixes related to the correlations
between the DM parameter $p$ and each of the non-DM parameters
$A_{\text{1V}}$, $A_{\text{2V}}$, $\Gamma_{\text{BL}}$,
$\Gamma_{\text{FSRQ}}$, $\Gamma_{\text{mAGN}}$ and $\Gamma_{\text{SFG}}$. In
Fig. 16 we show the bivariate correlations for a WIMP with a mass of 100 GeV.
In general, no strong correlation is present and the non-DM parameters have
little impact on the constraining power of $p$ (and consequently on that of
the DM decay lifetime).

## References

  * [1] S. Camera, M. Fornasa, N. Fornengo and M. Regis, _A Novel Approach in the Weakly Interacting Massive Particle Quest: Cross-correlation of Gamma-Ray Anisotropies and Cosmic Shear_ , _Astrophys. J. Lett._ 771 (2013) L5 [1212.5018].
  * [2] S. Camera, M. Fornasa, N. Fornengo and M. Regis, _Tomographic-spectral approach for dark matter detection in the cross-correlation between cosmic shear and diffuse $\gamma$-ray emission_, _JCAP_ 06 (2015) 029 [1411.4651].
  * [3] N. Fornengo and M. Regis, _Particle dark matter searches in the anisotropic sky_ , _Front. Physics_ 2 (2014) 6 [1312.4835].
  * [4] J.-Q. Xia, A. Cuoco, E. Branchini and M. Viel, _Tomography of the Fermi-lat $\gamma$-ray Diffuse Extragalactic Signal via Cross Correlations With Galaxy Catalogs_, _Astrophys. J. Suppl._ 217 (2015) 15 [1503.05918].
  * [5] A. Cuoco, J.-Q. Xia, M. Regis, E. Branchini, N. Fornengo and M. Viel, _Dark Matter Searches in the Gamma-ray Extragalactic Background via Cross-correlations With Galaxy Catalogs_ , _Astrophys. J. Suppl._ 221 (2015) 29 [1506.01030].
  * [6] M. Fornasa et al., _Angular power spectrum of the diffuse gamma-ray emission as measured by the Fermi Large Area Telescope and constraints on its dark matter interpretation_ , _Phys. Rev. D_ 94 (2016) 123005 [1608.07289].
  * [7] M. Regis, J.-Q. Xia, A. Cuoco, E. Branchini, N. Fornengo and M. Viel, _Particle dark matter searches outside the Local Group_ , _Phys. Rev. Lett._ 114 (2015) 241301 [1503.05922].
  * [8] J.-Q. Xia, A. Cuoco, E. Branchini, M. Fornasa and M. Viel, _A cross-correlation study of the Fermi-LAT $\gamma$-ray diffuse extragalactic signal_, _Mon. Not. Roy. Astron. Soc._ 416 (2011) 2247 [1103.4861].
  * [9] M. Shirasaki, S. Horiuchi and N. Yoshida, _Cross-Correlation of the Extragalactic Gamma-ray Background with Luminous Red Galaxies_ , _Phys. Rev. D_ 92 (2015) 123540 [1511.07092].
  * [10] A. Cuoco, M. Bilicki, J.-Q. Xia and E. Branchini, _Tomographic imaging of the Fermi-LAT gamma-ray sky through cross-correlations: A wider and deeper look_ , _Astrophys. J. Suppl._ 232 (2017) 10 [1709.01940].
  * [11] S. Ammazzalorso, N. Fornengo, S. Horiuchi and M. Regis, _Characterizing the local gamma-ray Universe via angular cross-correlations_ , _Phys. Rev. D_ 98 (2018) 103007 [1808.09225].
  * [12] D. Hashimoto, O. Macias, A.J. Nishizawa, K. Hayashi, M. Takada, M. Shirasaki et al., _Constraining dark matter annihilation with HSC Low Surface Brightness Galaxies_ , _JCAP_ 01 (2020) 059 [1906.06701].
  * [13] E. Branchini, S. Camera, A. Cuoco, N. Fornengo, M. Regis, M. Viel et al., _Cross-correlating the $\gamma$-ray sky with Catalogs of Galaxy Clusters_, _Astrophys. J. Suppl._ 228 (2017) 8 [1612.05788].
  * [14] D. Hashimoto, A.J. Nishizawa, M. Shirasaki, O. Macias, S. Horiuchi, H. Tashiro et al., _Measurement of redshift dependent cross correlation of HSC clusters and Fermi $\gamma$ rays_, 1805.08139.
  * [15] M. Colavincenzo, X. Tan, S. Ammazzalorso, S. Camera, M. Regis, J.-Q. Xia et al., _Searching for gamma-ray emission from galaxy clusters at low redshift_ , _Mon. Not. Roy. Astron. Soc._ 491 (2020) 3225 [1907.05264].
  * [16] X. Tan, M. Colavincenzo and S. Ammazzalorso, _Bounds on WIMP dark matter from galaxy clusters at low redshift_ , _Mon. Not. Roy. Astron. Soc._ 495 (2020) 114 [1907.06905].
  * [17] N. Fornengo, L. Perotto, M. Regis and S. Camera, _Evidence of Cross-correlation between the CMB Lensing and the $\Gamma$-ray sky_, _Astrophys. J. Lett._ 802 (2015) L1 [1410.4997].
  * [18] C. Feng, A. Cooray and B. Keating, _Planck Lensing and Cosmic Infrared Background Cross-Correlation with Fermi-LAT: Tracing Dark Matter Signals in the Gamma-Ray Background_ , _Astrophys. J._ 836 (2017) 127 [1608.04351].
  * [19] M. Shirasaki, S. Horiuchi and N. Yoshida, _Cross-Correlation of Cosmic Shear and Extragalactic Gamma-ray Background: Constraints on the Dark Matter Annihilation Cross-Section_ , _Phys. Rev. D_ 90 (2014) 063502 [1404.5503].
  * [20] T. Tröster et al., _Cross-correlation of weak lensing and gamma rays: implications for the nature of dark matter_ , _Mon. Not. Roy. Astron. Soc._ 467 (2017) 2706 [1611.03554].
  * [21] M. Shirasaki, O. Macias, S. Horiuchi, S. Shirai and N. Yoshida, _Cosmological constraints on dark matter annihilation and decay: Cross-correlation analysis of the extragalactic $\gamma$-ray background and cosmic shear_, _Phys. Rev. D_ 94 (2016) 063522 [1607.02187].
  * [22] M. Shirasaki, O. Macias, S. Horiuchi, N. Yoshida, C.-H. Lee and A.J. Nishizawa, _Correlation of extragalactic $\gamma$ rays with cosmic matter density distributions from weak gravitational lensing_, _Phys. Rev. D_ 97 (2018) 123015 [1802.10257].
  * [23] DES collaboration, _Detection of Cross-Correlation between Gravitational Lensing and $\gamma$ Rays_, _Phys. Rev. Lett._ 124 (2020) 101102 [1907.13484].
  * [24] F. Zandanel, C. Weniger and S. Ando, _The role of the eROSITA all-sky survey in searches for sterile neutrino dark matter_ , _JCAP_ 09 (2015) 060 [1505.07829].
  * [25] A. Caputo, M. Regis and M. Taoso, _Searching for Sterile Neutrino with X-ray Intensity Mapping_ , _JCAP_ 03 (2020) 001 [1911.09120].
  * [26] P. Serra, G. Lagache, O. Doré, A. Pullen and M. White, _Cross-correlation of cosmic far-infrared background anisotropies with large scale structures_ , _Astron. Astrophys._ 570 (2014) A98 [1404.1933].
  * [27] Y. Gong, A. Cooray, K. Mitchell-Wynne, X. Chen, M. Zemcov and J. Smidt, _Axion decay and anisotropy of near-IR extragalactic background light_ , _Astrophys. J._ 825 (2016) 104 [1511.01577].
  * [28] A. Caputo, A. Vittino, N. Fornengo, M. Regis and M. Taoso, _Searching for axion-like particle decay in the near-infrared background: an updated analysis_ , _JCAP_ 05 (2021) 046 [2012.09179].
  * [29] E. Pinetti, S. Camera, N. Fornengo and M. Regis, _Synergies across the spectrum for particle dark matter indirect detection: how HI intensity mapping meets gamma rays_ , _JCAP_ 07 (2020) 044 [1911.04989].
  * [30] D.C. Pan, M.S. Vogeley, F. Hoyle, Y.-Y. Choi and C. Park, _Cosmic Voids in Sloan Digital Sky Survey Data Release 7_ , _Mon. Not. Roy. Astron. Soc._ 421 (2012) 926 [1103.4156].
  * [31] K.A. Douglass, D. Veyrat and S. BenZvi, _Updated void catalogs of the SDSS DR7 main sample_ , 2202.01226.
  * [32] Q. Mao et al., _A Cosmic Void Catalog of SDSS DR12 BOSS Galaxies_ , _Astrophys. J._ 835 (2017) 161 [1602.02771].
  * [33] R. Voivodic, H. Rubira and M. Lima, _The Halo Void (Dust) Model of Large Scale Structure_ , _JCAP_ 10 (2020) 033 [2003.06411].
  * [34] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [35] D. Limber, _The analysis of counts of the extragalactic nebulae in terms of a fluctuating density field_ , _Astrophys. J._ 134 (1953) A6.
  * [36] N. Kaiser, _Weak gravitational lensing of distant galaxies_ , _Astrophys. J._ 388 (1992) 272.
  * [37] N. Kaiser, _Weak lensing and cosmology_ , _Astrophys. J._ 498 (1998) 26 [astro-ph/9610120].
  * [38] A. Cooray and R.K. Sheth, _Halo Models of Large Scale Structure_ , _Phys. Rept._ 372 (2002) 1 [astro-ph/0206508].
  * [39] F. Bernardeau, S. Colombi, E. Gaztanaga and R. Scoccimarro, _Large scale structure of the universe and cosmological perturbation theory_ , _Phys. Rept._ 367 (2002) 1 [astro-ph/0112551].
  * [40] D. Baumann, A. Nicolis, L. Senatore and M. Zaldarriaga, _Cosmological Non-Linearities as an Effective Fluid_ , _JCAP_ 07 (2012) 051 [1004.2488].
  * [41] J.J.M. Carrasco, M.P. Hertzberg and L. Senatore, _The Effective Field Theory of Cosmological Large Scale Structures_ , _JHEP_ 09 (2012) 082 [1206.2926].
  * [42] J.J.M. Carrasco, S. Foreman, D. Green and L. Senatore, _The Effective Field Theory of Large Scale Structures at Two Loops_ , _JCAP_ 07 (2014) 057 [1310.0464].
  * [43] T. Konstandin, R.A. Porto and H. Rubira, _The effective field theory of large scale structure at three loops_ , _JCAP_ 11 (2019) 027 [1906.00997].
  * [44] R. Teyssier, _Cosmological hydrodynamics with adaptive mesh refinement: a new high resolution code called ramses_ , _Astron. Astrophys._ 385 (2002) 337 [astro-ph/0111367].
  * [45] V. Springel, _The Cosmological simulation code GADGET-2_ , _Mon. Not. Roy. Astron. Soc._ 364 (2005) 1105 [astro-ph/0505010].
  * [46] B. Moore, T.R. Quinn, F. Governato, J. Stadel and G. Lake, _Cold collapse and the core catastrophe_ , _Mon. Not. Roy. Astron. Soc._ 310 (1999) 1147 [astro-ph/9903164].
  * [47] J.F. Navarro, C.S. Frenk and S.D.M. White, _The Structure of cold dark matter halos_ , _Astrophys. J._ 462 (1996) 563 [astro-ph/9508025].
  * [48] VIRGO collaboration, _Galaxy clusters in Hubble volume simulations: Cosmological constraints from sky survey populations_ , _Astrophys. J._ 573 (2002) 7 [astro-ph/0110246].
  * [49] R.E. Angulo, C.M. Baugh, C.S. Frenk and C.G. Lacey, _Extending the halo mass resolution of $N$-body simulations_, _Mon. Not. Roy. Astron. Soc._ 442 (2014) 3256 [1310.3880].
  * [50] LSST collaboration, _LSST: from Science Drivers to Reference Design and Anticipated Data Products_ , _Astrophys. J._ 873 (2019) 111 [0805.2366].
  * [51] DES collaboration, _The Dark Energy Survey_ , astro-ph/0510346.
  * [52] Euclid Theory Working Group collaboration, _Cosmology and fundamental physics with the Euclid satellite_ , _Living Rev. Rel._ 16 (2013) 6 [1206.1225].
  * [53] R.A. Sunyaev and Y.B. Zeldovich, _The Velocity of clusters of galaxies relative to the microwave background. The Possibility of its measurement_ , _Mon. Not. Roy. Astron. Soc._ 190 (1980) 413.
  * [54] ACT collaboration, _Atacama Cosmology Telescope measurements of a large sample of candidates from the Massive and Distant Clusters of WISE Survey - Sunyaev-Zeldovich effect confirmation of MaDCoWS candidates using ACT_ , _Astron. Astrophys._ 653 (2021) A135 [2105.00068].
  * [55] R.K. Sachs and A.M. Wolfe, _Perturbations of a cosmological model and angular variations of the microwave background_ , _Astrophys. J._ 147 (1967) 73.
  * [56] N. Lemarchand, J. Grain, G. Hurier, F. Lacasa and A. Ferté, _Secondary CMB anisotropies from magnetized haloes - I. Power spectra of the Faraday rotation angle and conversion rate_ , _Astron. Astrophys._ 630 (2019) A149 [1810.09221].
  * [57] A. Mead, J. Peacock, C. Heymans, S. Joudaki and A. Heavens, _An accurate halo model for fitting non-linear cosmological power spectra and baryonic feedback models_ , _Mon. Not. Roy. Astron. Soc._ 454 (2015) 1958 [1505.07833].
  * [58] F. Schmidt, _Towards a self-consistent halo model for the nonlinear large-scale structure_ , _Phys. Rev. D_ 93 (2016) 063512 [1511.02231].
  * [59] A.Y. Chen and N. Afshordi, _Amending the halo model to satisfy cosmological conservation laws_ , _Phys. Rev. D_ 101 (2020) 103522 [1912.04872].
  * [60] P. Valageas and T. Nishimichi, _Combining perturbation theories with halo models for the matter bispectrum_ , _Astron. Astrophys._ 532 (2011) A4 [1102.0641].
  * [61] R. Biswas, E. Alizadeh and B.D. Wandelt, _Voids as a Precision Probe of Dark Energy_ , _Phys. Rev. D_ 82 (2010) 023002 [1002.0014].
  * [62] G. Pollina, M. Baldi, F. Marulli and L. Moscardini, _Cosmic voids in coupled dark energy cosmologies: the impact of halo bias_ , _Mon. Not. Roy. Astron. Soc._ 455 (2016) 3075 [1506.08831].
  * [63] M. Sahlén, I.n. Zubeldía and J. Silk, _Cluster–void Degeneracy Breaking: Dark Energy, Planck, and the Largest Cluster and Void_ , _Astrophys. J. Lett._ 820 (2016) L7 [1511.04075].
  * [64] A. Pisani, P.M. Sutter, N. Hamaus, E. Alizadeh, R. Biswas, B.D. Wandelt et al., _Counting voids to probe dark energy_ , _Phys. Rev. D_ 92 (2015) 083531 [1503.07690].
  * [65] A. Barreira, M. Cautun, B. Li, C. Baugh and S. Pascoli, _Weak lensing by voids in modified lensing potentials_ , _JCAP_ 08 (2015) 028 [1505.05809].
  * [66] Y.-C. Cai, N. Padilla and B. Li, _Testing Gravity using Cosmic Voids_ , _Mon. Not. Roy. Astron. Soc._ 451 (2015) 1036 [1410.1510].
  * [67] E.L.D. Perico, R. Voivodic, M. Lima and D.F. Mota, _Cosmic voids in modified gravity scenarios_ , _Astron. Astrophys._ 632 (2019) A52 [1905.12450].
  * [68] R. Voivodic, M. Lima, C. Llinares and D.F. Mota, _Modelling Void Abundance in Modified Gravity_ , _Phys. Rev. D_ 95 (2017) 024018 [1609.02544].
  * [69] G. D’Amico, M. Musso, J. Norena and A. Paranjape, _Excursion Sets and Non-Gaussian Void Statistics_ , _Phys. Rev. D_ 83 (2011) 023521 [1011.1229].
  * [70] K.C. Chan, N. Hamaus and M. Biagetti, _Constraint of Void Bias on Primordial non-Gaussianity_ , _Phys. Rev. D_ 99 (2019) 121304 [1812.04024].
  * [71] H. Song and J. Lee, _The Mass Distribution of SDSS Galaxy Groups in Void Regions and Its Implication on the Primordial non-Gaussianity_ , _Astrophys. J. Lett._ 701 (2009) L25 [0811.1339].
  * [72] Y.-C. Cai, A. Taylor, J.A. Peacock and N. Padilla, _Redshift-space distortions around voids_ , _Mon. Not. Roy. Astron. Soc._ 462 (2016) 2465 [1603.05184].
  * [73] T. Chantavat, U. Sawangwit, P.M. Sutter and B.D. Wandelt, _Cosmological parameter constraints from CMB lensing with cosmic voids_ , _Phys. Rev. D_ 93 (2016) 043523 [1409.3364].
  * [74] N. Hamaus, B.D. Wandelt, P.M. Sutter, G. Lavaux and M.S. Warren, _Cosmology with Void-Galaxy Correlations_ , _Phys. Rev. Lett._ 112 (2014) 041304 [1307.2571].
  * [75] N. Hamaus, P.M. Sutter, G. Lavaux and B.D. Wandelt, _Probing cosmology and gravity with redshift-space distortions around voids_ , _JCAP_ 11 (2015) 036 [1507.04363].
  * [76] N. Hamaus, A. Pisani, P.M. Sutter, G. Lavaux, S. Escoffier, B.D. Wandelt et al., _Constraints on Cosmology and Gravity from the Dynamics of Voids_ , _Phys. Rev. Lett._ 117 (2016) 091302 [1602.01784].
  * [77] C.-H. Chuang, F.-S. Kitaura, Y. Liang, A. Font-Ribera, C. Zhao, P. McDonald et al., _Linear redshift space distortions for cosmic voids based on galaxies in redshift space_ , _Phys. Rev. D_ 95 (2017) 063528 [1605.05352].
  * [78] G. Lavaux and B.D. Wandelt, _Precision cosmography with stacked voids_ , _Astrophys. J._ 754 (2012) 109 [1110.0345].
  * [79] J.R. Bond, S. Cole, G. Efstathiou and N. Kaiser, _Excursion set mass functions for hierarchical Gaussian fluctuations_ , _Astrophys. J._ 379 (1991) 440.
  * [80] R.K. Sheth and R. van de Weygaert, _A Hierarchy of voids: Much ado about nothing_ , _Mon. Not. Roy. Astron. Soc._ 350 (2004) 517 [astro-ph/0311260].
  * [81] A. De Simone, M. Maggiore and A. Riotto, _Excursion Set Theory for generic moving barriers and non-Gaussian initial conditions_ , _Mon. Not. Roy. Astron. Soc._ 412 (2011) 2587 [1007.1903].
  * [82] E. Jennings, Y. Li and W. Hu, _The abundance of voids and the excursion set formalism_ , _Mon. Not. Roy. Astron. Soc._ 434 (2013) 2167 [1304.6087].
  * [83] N. Hamaus, P.M. Sutter and B.D. Wandelt, _Universal Density Profile for Cosmic Voids_ , _Phys. Rev. Lett._ 112 (2014) 251302 [1403.5499].
  * [84] A.R. Duffy, J. Schaye, S.T. Kay and C. Dalla Vecchia, _Dark matter halo concentrations in the Wilkinson Microwave Anisotropy Probe year 5 cosmology_ , _Mon. Not. Roy. Astron. Soc._ 390 (2008) L64 [0804.2486].
  * [85] R.K. Sheth and G. Tormen, _Large scale bias and the peak background split_ , _Mon. Not. Roy. Astron. Soc._ 308 (1999) 119 [astro-ph/9901122].
  * [86] W.H. Press and P. Schechter, _Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensation_ , _Astrophys. J._ 187 (1974) 425.
  * [87] R.K. Sheth, H.J. Mo and G. Tormen, _Ellipsoidal collapse and an improved model for the number and spatial distribution of dark matter haloes_ , _Mon. Not. Roy. Astron. Soc._ 323 (2001) 1 [astro-ph/9907024].
  * [88] A.A. Berlind and D.H. Weinberg, _The Halo occupation distribution: Towards an empirical determination of the relation between galaxies and mass_ , _Astrophys. J._ 575 (2002) 587 [astro-ph/0109001].
  * [89] A.A. Berlind, D.H. Weinberg, A.J. Benson, C.M. Baugh, S. Cole, R. Dave et al., _The Halo occupation distribution and the physics of galaxy formation_ , _Astrophys. J._ 593 (2003) 1 [astro-ph/0212357].
  * [90] Z. Zheng, A.L. Coil and I. Zehavi, _Galaxy Evolution from Halo Occupation Distribution Modeling of DEEP2 and SDSS Galaxy Clustering_ , _Astrophys. J._ 667 (2007) 760 [astro-ph/0703457].
  * [91] Z. Zheng, A.A. Berlind, D.H. Weinberg, A.J. Benson, C.M. Baugh, S. Cole et al., _Theoretical models of the halo occupation distribution: Separating central and satellite galaxies_ , _Astrophys. J._ 633 (2005) 791 [astro-ph/0408564].
  * [92] 2MASS collaboration, _The Two Micron All Sky Survey (2MASS)_ , _Astron. J._ 131 (2006) 1163.
  * [93] SDSS collaboration, _The Luminosity and color dependence of the galaxy correlation function_ , _Astrophys. J._ 630 (2005) 1 [astro-ph/0408569].
  * [94] S.G. Patiri, F. Prada, J. Holtzman, A. Klypin and J. Betancort-Rijo, _The Properties of Galaxies in Voids_ , _Mon. Not. Roy. Astron. Soc._ 372 (2006) 1710 [astro-ph/0605703].
  * [95] A. Constantin, F. Hoyle and M.S. Vogeley, _Active Galactic Nuclei in Void Regions_ , _Astrophys. J._ 673 (2008) 715 [0710.1631].
  * [96] A. Furniss, P.M. Sutter, J.R. Primack and A. Domínguez, _A Correlation Between Hard Gamma-ray Sources and Cosmic Voids Along the Line of Sight_ , _Mon. Not. Roy. Astron. Soc._ 446 (2015) 2267 [1407.6370].
  * [97] M. Ajello et al., _The Cosmic Evolution of Fermi BL Lacertae Objects_ , _Astrophys. J._ 780 (2014) 73 [1310.0006].
  * [98] M. Ajello et al., _The Luminosity Function of Fermi-detected Flat-Spectrum Radio Quasars_ , _Astrophys. J._ 751 (2012) 108 [1110.3787].
  * [99] C.J. Willott, S. Rawlings, K.M. Blundell, M. Lacy and S.A. Eales, _The radio luminosity function from the low-frequency 3crr, 6ce & 7crs complete samples_, _Mon. Not. Roy. Astron. Soc._ 322 (2001) 536 [astro-ph/0010419].
  * [100] M. Di Mauro, F. Calore, F. Donato, M. Ajello and L. Latronico, _Diffuse $\gamma$-ray emission from misaligned active galactic nuclei_, _Astrophys. J._ 780 (2014) 161 [1304.0908].
  * [101] C. Gruppioni et al., _The Herschel PEP/HerMES Luminosity Function. I: Probing the Evolution of PACS selected Galaxies to z~4_ , _Mon. Not. Roy. Astron. Soc._ 432 (2013) 23 [1302.5209].
  * [102] M. Bartelmann, _Gravitational Lensing_ , _Class. Quant. Grav._ 27 (2010) 233001 [1010.3829].
  * [103] Euclid collaboration, _Euclid preparation: VII. Forecast validation for Euclid cosmological probes_ , _Astron. Astrophys._ 642 (2020) A191 [1910.09273].
  * [104] N. Fornengo, L. Maccione and A. Vittino, _Constraints on particle dark matter from cosmic-ray antiprotons_ , _JCAP_ 04 (2014) 003 [1312.3579].
  * [105] C. Blanco and D. Hooper, _Constraints on Decaying Dark Matter from the Isotropic Gamma-Ray Background_ , _JCAP_ 03 (2019) 019 [1811.05988].
  * [106] J.A.R. Cembranos, A. de la Cruz-Dombriz, A. Dobado, R.A. Lineros and A.L. Maroto, _Photon spectra from WIMP annihilation_ , _Phys. Rev. D_ 83 (2011) 083507 [1009.4936].
  * [107] S. Razzaque, C.D. Dermer and J.D. Finke, _The stellar contribution to the extra-galactic background light and absorption of TeV gamma-rays_ , _Astrophys. J._ 697 (2009) 483 [0807.4294].
  * [108] Fermi-LAT collaboration, _Unresolved Gamma-Ray Sky through its Angular Power Spectrum_ , _Phys. Rev. Lett._ 121 (2018) 241101 [1812.02079].
  * [109] Fermi-LAT collaboration, _The spectrum of isotropic diffuse gamma-ray emission between 100 MeV and 820 GeV_ , _Astrophys. J._ 799 (2015) 86 [1410.3696].
  * [110] Fermi-LAT collaboration, _Anisotropies in the diffuse gamma-ray background measured by the Fermi LAT_ , _Phys. Rev. D_ 85 (2012) 083007 [1202.2856].
  * [111] N. Topchiev et al., _High-energy gamma-ray studying with GAMMA-400_ , _PoS_ ICRC2017 (2018) 802 [1707.04882].
  * [112] DES collaboration, _Dark Energy Survey Year 3 results: Cosmological constraints from galaxy clustering and weak lensing_ , _Phys. Rev. D_ 105 (2022) 023520 [2105.13549].
  * [113] DES collaboration, _Dark Energy Survey Year 3 results: Cosmological constraints from galaxy clustering and weak lensing_ , _Phys. Rev. D_ 105 (2022) 023520 [2105.13549].
  * [114] F. Donato, N. Fornengo and D. Maurin, _Antideuteron fluxes from dark matter annihilation in diffusion models_ , _Phys. Rev. D_ 78 (2008) 043506 [0803.2640].
  * [115] T. Aramaki et al., _Review of the theoretical and experimental status of dark matter identification with cosmic-ray antideuterons_ , _Phys. Rept._ 618 (2016) 1 [1505.07785].
  * [116] H. Jeffreys, _Theory of Probability_ , Oxford Classic Texts in the Physical Sciences, Oxford University Press,, Oxford, U.K. (1939).
  * [117] M.S. Vogeley and A.S. Szalay, _Eigenmode analysis of galaxy redshift surveys I. theory and methods_ , _Astrophys. J._ 465 (1996) 34 [astro-ph/9601185].
  * [118] M. Tegmark, A. Taylor and A. Heavens, _Karhunen-Loeve eigenvalue problems in cosmology: How should we tackle large data sets?_ , _Astrophys. J._ 480 (1997) 22 [astro-ph/9603021].
  * [119] L. Lara, G. Giovannini, W.D. Cotton, L. Feretti, J.M. Marcaide, I. Marquez et al., _A New sample of large angular size radio galaxies. 3. Statistics and evolution of the grown population_ , _Astron. Astrophys._ 421 (2004) 899 [astro-ph/0404373].
  * [120] Y. Inoue, _Contribution of the Gamma-ray Loud Radio Galaxies Core Emissions to the Cosmic MeV and GeV Gamma-Ray Background Radiation_ , _Astrophys. J._ 733 (2011) 66 [1103.3946].
  * [121] Fermi-LAT collaboration, _GeV Observations of Star-forming Galaxies with \textitFermi LAT_ , _Astrophys. J._ 755 (2012) 164 [1206.1346].