# Exploring the cosmological synergy between galaxy cluster and cosmic void
number counts

D. Pelliciari,1,2 S. Contarini,1,3,4 F. Marulli,1,3,4 L. Moscardini,1,3,4 C.
Giocoli,3,4 G. F. Lesci,1,3 and K. Dolag 6,7

1 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum
Università di Bologna, via Piero Gobetti 93/2, I-40129 Bologna, Italy
2 Istituto Nazionale di Astrofisica (INAF) – Istituto di Radioastronomia
(IRA), via Gobetti 101, 40129 Bologna, Italy
3 Istituto Nazionale di Astrofisica (INAF) - Osservatorio di Astrofisica e
Scienza dello Spazio (OAS), via Piero Gobetti 93/3, I-40129 Bologna, Italy
4 Istituto Nazionale di Fisica Nucleare (INFN) - Sezione di Bologna, viale
Berti Pichat 6/2, I-40127 Bologna, Italy
6 Universitats-Sternwarte, Fakultat fur Physik, Ludwig-Maximilians-Universitat
Munchen, Scheinerstr.1, 81679 Munchen, Germany
7 Max-Planck-Institut fur Astrophysik, Karl-Scharzschildstr. 1, D-85748,
Garching, Germany E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Galaxy clusters and cosmic voids, the most extreme objects of our Universe in
terms of mass and size, trace two opposite sides of the large-scale matter
density field. By studying their abundance as a function of their mass and
radius, respectively, i.e. the halo mass function (HMF) and void size function
(VSF), it is possible to achieve fundamental constraints on the cosmological
model. While the HMF has already been extensively exploited providing robust
constraints on the main cosmological model parameters (e.g. $\Omega_{\rm m}$,
$\sigma_{8}$ and $S_{8}$), the VSF is still emerging as a viable and effective
cosmological probe. Given the expected complementarity of these statistics, in
this work we aim at estimating the costraining power deriving from their
combination. To this end, we exploit realistic mock samples of galaxy clusters
and voids extracted from state-of-the-art large hydrodynamical simulations, in
the redshift range $0.2\leq z\leq 1$. We perform an accurate calibration of
the free parameters of the HMF and VSF models, needed to take into account the
differences between the types of mass tracers used in this work and those
considered in previous literature analyses. Then, we obtain constraints on
$\Omega_{\rm m}$ and $\sigma_{8}$ by performing a Bayesian analysis. We find
that cluster and void counts represent powerful independent and complementary
probes to test the cosmological framework. In particular, the constraining
power of the HMF on $\Omega_{\rm m}$ and $\sigma_{8}$ improves drastically
with the VSF contribution, increasing the $S_{8}$ constraint precision by a
factor of about $60\%$.

###### keywords:

cosmological parameters – large-scale structure of Universe – methods:
numerical

††pubyear: 2022††pagerange: Exploring the cosmological synergy between galaxy
cluster and cosmic void number counts–A

##  1 Introduction

Thanks to the remarkable achievements obtained by experiments on the
cosmological microwave background (CMB, see e.g. Komatsu et al., 2011; Planck
Collaboration et al., 2020) as well as by spectroscopic and photometric galaxy
surveys (e.g. Alam et al., 2017), cosmology has recently entered a new
precision era, in which cosmological parameters are constrained with sub-
percent accuracy. The $\Lambda$-cold dark matter ($\Lambda$CDM) model, being
described by only six parameters, got strong empirical support by the current
cosmological observations. In this framework, the matter-energy budget of our
Universe is dominated by two dark components, i.e. dark energy (DE),
responsible for its present accelerated expansion (Riess et al., 1998;
Perlmutter et al., 1999), and cold dark matter (CDM), which is the dominant
matter component in galaxies and galaxy clusters.

However, statistically significant tensions between the cosmological
constraints derived by early-and late-time measurements have also arisen in
the last years (see e.g. Di Valentino et al., 2020a, b, c; Perivolaropoulos &
Skara, 2021, for recent reviews). Among those, the most relevant ones are the
discrepancies between the derived values of the Hubble constant, $H_{0}$, and
the tension on the growth of cosmic structures, often quantified in terms of
the parameter $S_{8}\equiv\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$, where
$\Omega_{\rm m}$ is the matter density parameter and $\sigma_{8}$ is related
to the amplitude of the matter power spectrum. Many works have been focused on
these tensions (see e.g. Abdalla et al., 2022, for a recent review of
astrophysical and cosmological tensions). In recent years numerous theoretical
solutions have been proposed, going from modified gravity theories to
alternative DE theories (see e.g. Yoo & Watanabe, 2012; Amendola et al., 2013;
Joyce et al., 2016, for an extensive review).

To shed light on these tensions it is fundamental to consider different and
complementary cosmological probes, whose synergy can improve the precision on
the parameter measurements. In fact, the highest precision is achieved by
considering probe combination, i.e. computing the joint constraints from
different analyses, such as e.g. combining CMB anisotropies and large-scale
strucures observations (Webster et al., 1998; Gawiser & Silk, 1998; Bridle et
al., 1999), as well as supernovae Type-Ia (SNIa), baryonic acoustic
oscillations, weak gravitational lensing and galaxy clustering (see e.g. DES
Collaboration, 2019). Moreover, Wu et al. (2022) recently discussed the
constraining power of the combination between 21 cm intensity mapping, fast
radio bursts, gravitational wave standard sirens and strong gravitational
lensing, which are novel probes that will be greatly developed in the near
future (see Moresco et al., 2022, for a recent review about alternative
cosmological probes).

In terms of mass and size, galaxy clusters and cosmic voids are the most rare
and extreme objects in the Universe. Very massive clusters and large voids
follow the exponential tail of their mass/size distribution functions, which
makes their abundance a probe very sensitive to the change of cosmological
parameters. It is well known that galaxy clusters, associated with massive DM
haloes, provide a powerful cosmological probe when exploiting their abundances
(Vikhlinin et al., 2009; Cunha et al., 2009; Balaguera-Antolínez & Porciani,
2013; Hasselfield et al., 2013; Planck Collaboration et al., 2014, 2016;
Bocquet et al., 2019; Costanzi et al., 2019; Lesci et al., 2022a), as well as
their clustering properties (Hong et al., 2012; Veropalumbo et al., 2016;
Marulli et al., 2018, 2021; Lesci et al., 2022b). Also, cosmological
constraints become even more stringent when cluster counts are combined with
clustering measurements (see e.g. Mana et al., 2013; Salvati et al., 2018),

Concerning cosmic voids, their number counts (Platen et al., 2008; Pisani et
al., 2015b; Contarini et al., 2019, 2022a), density profiles (Hamaus et al.,
2014; Ricciardelli et al., 2014; Nadathur et al., 2015), and redshift-space
distortions of void-galaxy cross-correlation (Nadathur et al., 2020; Hamaus et
al., 2020; Hamaus et al., 2022), are rapidly becoming competitive and
alternative cosmological probes. Indeed, their unique low density interiors
and large sizes make these objects powerful tools to study the properties of
the DE (Bos et al., 2012; Pisani et al., 2015b; Verza et al., 2019) and the
effects of modified gravity theories (Spolyar et al., 2013; Barreira et al.,
2015; Contarini et al., 2021), to which they are highly sensitive. Cosmic
voids are also well suited to test primordial non-Gaussianity (Chan et al.,
2019), and other physical phenomena beyond the Standard Model of particle
physics (Reed et al., 2015; Yang et al., 2015; Baldi & Villaescusa-Navarro,
2016), and to investigate elusive components such as massive neutrinos
(Villaescusa-Navarro et al., 2013; Massara et al., 2015; Schuster et al.,
2019; Kreisch et al., 2019; Contarini et al., 2021), which are predominant in
these regions.

The number counts of galaxy clusters and cosmic voids can be considered as
independent and complementary cosmological probes given the opposite regimes
in sides of the matter density field they map. We then expect some degrees of
orthogonality between the corresponding cosmological constraints (Sahlén et
al., 2016; Bayer et al., 2021; Kreisch et al., 2022). This important feature
represents a desired property to extract as much information as possible from
the joint combination of these two probes, and it has not been exploited in
great details so far.

In this work we investigate the constraining power of a joint analysis of the
number counts of galaxy clusters and cosmic voids. We consider redshift-space
cluster and void catalogues extracted from a large hydrodynamic cosmological
simulation in four redshift bins, spanning in the range $0.2\leq z\leq 1$.
Recent works have already studied the constraining power derived from the
combination of halo and void statistics obtained from cosmological simulations
(Bayer et al., 2021; Kreisch et al., 2022). Nevertheless, the analysis we
present in this work is directly applicable to real data sets, thanks to the
bias-dependent void size function parametrisation (Section 2.2), the inclusion
of the most relevant observational effects (Section 4.1), the model
calibrations on hydrodynamical simulations (Section 5.1), and the both
realistic and conserving mass and radius selections considered for clusters
and voids.

The paper is structured as follows. In Section 2 we introduce the theoretical
context of the analysis we carry out in this work, focusing on the considered
halo mass function and void size function models. In Section 3 we describe the
simulations used to construct cluster and void mock catalogues, along with the
prescriptions for a reliable data set preparation. In Section 4 we introduce
the observational effects considered in the simulation analysis, and we
present the Bayesian framework used to assess the posterior distributions of
the model parameters. We also estimate the correlation between the two
cosmological probes considered, verifying their statistical independence.
Subsequently, we outline in Section 5 the adopted calibration method
considered for the analysis, required for an accurate comparison between the
data and the theoretical models, and we show the results of the combined
analysis of cluster and void number counts. Finally, we draw our conclusions
in Section 6.

##  2 Cosmological probes

The abundance of both haloes and cosmic voids can be modelled with the
excursion-set formalism (Bond et al., 1991; Zentner, 2007), also known as the
extended Press-Schechter formalism (Press & Schechter, 1974), which can be
used to describe the evolution of spherical overdensities and underdensities.
In particular, haloes form from the growth of positive fluctuations, i.e.
overdensities, while cosmic voids form from the evolution of negative ones,
i.e. underdensities. In what follows, we describe the two theoretical models
used in this work, i.e. the halo mass function (HMF) and the void size
function (VSF), which describe the comoving number density of haloes and
voids, as a function of their mass and radius, respectively. We then compute
number counts from number density by taking into account the most relevant
observational effects (see Section 4.1).

###  2.1 Cluster counts

Following the Press-Schechter formalism, the comoving number density,
$n(M,z)$, of DM haloes having masses between $M$ and $M+\mathrm{d}M$ at
redshift $z$ can be written as (Press & Schechter, 1974; Cole, 1991; Bond et
al., 1991; Mo & White, 1996; Sheth & Tormen, 1999; Sheth et al., 2001):

$\frac{\mathrm{d}n(M,z)}{\mathrm{d}\ln
M}=\frac{\bar{\rho}}{M}f_{\ln\sigma}(M,z)\frac{\mathrm{d}\ln\sigma^{-1}(M,z)}{\mathrm{d}\ln
M}\textrm{ ,}$ (1)

when expressed through its logarithmic differential. Here $\bar{\rho}$ is the
mean cosmic background density and $f_{\ln\sigma}(M,z)$ is the so-called
multiplicity function, which depends on the mass definition adopted to
identify simulated haloes (see e.g. Tinker et al., 2008). For the Press-
Schechter mass function, the latter is defined as:

$f_{\ln\sigma}(M,z)=\frac{\mathrm{d}\sigma(M,z)}{\mathrm{d}\ln
M}=\sqrt{\frac{2}{\pi}}\frac{\delta_{\rm
c}}{\sigma_{M}}\exp\Bigl{(}-\frac{\delta_{\rm
c}^{2}}{2\sigma^{2}(M,z)}\Bigr{)}\textrm{ ,}$ (2)

where $\delta_{\rm c}$ is the critical overdensity threshold computed in the
linear theory. Here $\sigma_{M}$ is the square root of the variance of the
linear density field extrapolated to the redshift at which haloes are
identified, after smoothing with a spherical top-hat filter with a radius $R$
enclosing a mass $M$. This can be expressed in terms of the present-day matter
power spectrum $P(k)$:

$\sigma^{2}(M,z)=\frac{D^{2}(z)}{2\pi^{3}}\int\mathrm{d}^{3}\vec{k}P(\vec{k},z=0)\hat{W}^{2}[\vec{k};R(M)]\textrm{
,}$ (3)

where $D(z)$ is the growth factor of linear perturbations, normalised such
that $D(z=0)=1$ (Carroll et al., 1992), $\hat{W}[\vec{k};R(M)]$ is the
Fourier-transform of the window function and $R(M)$ is obtained by requiring
that $M/(4\pi R^{3}/3)=\bar{\rho}$. Here $\bar{\rho}$ is expressed in comoving
coordinates, so its value is independent of the redshift. The dependence of
the HMF on the cosmological parameters $\Omega_{\rm m}$ and $\sigma_{8}$ lies
both in the mass variance, $\sigma^{2}(M,z)$, and in the mean density of the
Universe, $\bar{\rho}$.

It can be shown that for a $\Lambda$CDM cosmology the HMF assumes a universal
functional form when the mass variance is parametrised by means of the
following scaled variable (Sheth & Tormen, 1999; Sheth et al., 2001; Courtin
et al., 2011; Despali et al., 2016):

$\nu\equiv\frac{\delta^{2}_{c}}{\sigma_{M}^{2}}\textrm{ .}$ (4)

With this parametrisation the HMF can be written as:

$\nu
f(\nu)=\frac{M^{2}}{\bar{\rho}}\frac{\mathrm{d}n}{\mathrm{d}M}\frac{\mathrm{d}\ln
M}{\mathrm{d}\ln\nu}\textrm{ ,}$ (5)

where:

$\nu f(\nu)=A\biggl{(}1+\frac{1}{\nu^{\prime
p}}\biggr{)}\biggl{(}\frac{\nu^{\prime}}{2\pi}\biggr{)}^{1/2}e^{-\nu^{\prime}/2}\textrm{
,}$ (6)

with $\nu^{\prime}=a\nu$. The parameters of this model are $(a,p,A)$ and can
be calibrated using cosmological simulations (Manera et al., 2010; Courtin et
al., 2011). In particular, $a$, $p$, and $A$ define the mass function cut-off
at high masses, its shape in the low mass range and its normalisation,
respectively. Despali et al. (2016) found that the HMF parametrisation of Eq.
(6) is universal if clusters are defined by means of the virial overdensity
$\Delta_{\rm vir}$. However, since both observed and simulated clusters are
usually defined via other overdensity criteria, an extended universal fitting
formula was provided by Despali et al. (2016) as follows:

$\displaystyle a$ $\displaystyle=0.4332x^{2}+0.2263x+0.7665$ (7)
$\displaystyle p$ $\displaystyle=-0.1151x^{2}+0.2554x+0.2488$ $\displaystyle
A$ $\displaystyle=-0.1362x+0.3292\textrm{ ,}$

where the quantity $x\equiv\log[\Delta(z)/\Delta_{\text{vir}}(z)]$ properly
rescales the spherical overdensity used to define the clusters, $\Delta(z)$,
with the virial overdensity, $\Delta_{\text{vir}}(z)$.

The HMF universality highlighted by Despali et al. (2016) is a fundamental
feature which can be exploited in cosmological analyses. Indeed, its universal
shape can be used to derive unbiased cosmological constraints by modelling
data at any redshift and for standard cosmological scenarios. In particular,
the coefficients in Eq. (7) have been calibrated considering DM-only
cosmological simulations, with haloes identified with a spherical overdensity
halo finder (e.g. Tormen, 1998; Tormen et al., 2004; Giocoli et al., 2008). A
re-calibration of the fitting relation parameters is therefore necessary when
using different type of simulations (e.g. hydrodynamic) and/or other halo
finder methods (see Section 3).

###  2.2 Void counts

As already mentioned, the abundance of cosmic voids as a function of their
radius can be modelled as well with the excursion-set formalism (Sheth & van
de Weygaert, 2004; Jennings et al., 2013). Given the analogies among
overdensity and underdensity, the linear VSF can be written as:

$\frac{\mathrm{d}n}{\mathrm{d}\ln R}\biggr{|}_{\rm
L}=\frac{f_{\ln\sigma}}{V(R_{\rm
L})}\frac{\mathrm{d}\ln\sigma^{-1}(R,z)}{\mathrm{d}\ln R_{\rm L}}\mathrm{,}$
(8)

where $V(R_{\rm L})=\frac{4}{3}\pi R_{\rm L}^{3}$ is the volume of a spherical
region of radius $R_{\rm L}$ and the subscript "L" indicates that the
quantities are obtained in linear theory. $\sigma^{2}(R,z)$ is the variance of
the matter power spectrum on a scale $R$ and, analogously to $\sigma^{2}(M,z)$
entering in the HMF, embeds the cosmological dependence of the model. Here the
multiplicity function can be expressed as an infinite series (Sheth & van de
Weygaert, 2004):

$f_{\ln\sigma}=2\sum_{j=1}^{\infty}j\pi
x^{2}\sin(j\pi\mathcal{D})\exp\biggl{[}-\frac{(j\pi
x)^{2}}{2}\biggr{]}\textrm{ ,}$ (9)

where:

$\mathcal{D}\equiv\frac{|\delta_{\rm v,L}|}{\delta_{\rm c}+|\delta_{\rm
v,L}|}\quad\textrm{,}\quad x\equiv\frac{\mathcal{D}}{|\delta_{\rm
v,L}|}\,\sigma(R,z)\textrm{ ,}$ (10)

and $\delta_{\rm v,L}$ is the linear underdensity threshold. Following the
spherical evolution, nonlinear voids have expanded by a factor of $\simeq 1.7$
with respect to their linear stage. Consequently, the nonlinear (subscript
"NL") void abundance becomes:

$\frac{\mathrm{d}n}{\mathrm{d}\ln R}\biggr{|}_{\rm
SvdW}=\frac{\mathrm{d}n_{\rm L}}{\mathrm{d}\ln R_{\rm L}}\biggl{|}_{R_{\rm
L}=R/1.7}\textrm{ .}$ (11)

The latter has been proposed by Sheth & van de Weygaert (2004) (SvdW,
hereafter) and takes into account the suppression impact of the halo formation
on the evolving population of voids (the void-in-cloud phenomenon) expressing
$\mathcal{D}$ as a function of the two threshold parameters $\delta_{\rm v}$
and $\delta_{\rm c}$. This translates into a shift of the linear model of Eq.
(8) towards larger radii, without any change in amplitude, which means that in
the SvdW model the comoving number density of voids is conserved in the
nonlinear regime (Sheth & van de Weygaert, 2004).

Jennings et al. (2013) pointed out that the latter is an incorrect assumption
especially for large voids. Following this model, the cumulative fraction of
volume occupied by voids exceeds the total volume of the Universe. Indeed, the
SvdW model overpredicts the abundance of voids at any radius (Jennings et al.,
2013; Ronconi & Marulli, 2017), unless the threshold $\delta_{\rm v,L}$ is
considered as a free parameter, which can be fine-tuned for different redshift
and cosmic tracers by means of cosmological simulations. However, this
severally affects the possibility of using the VSF as a cosmological probe.

In order to account for this problem, Jennings et al. (2013) proposed the
volume-conserving size function model, or Vdn, in which the volume fraction
occupied by cosmic voids is conserved in the transition from the linear to the
nonlinear regime. This VSF model, which has been tested recently in many works
(Jennings et al., 2013; Ronconi & Marulli, 2017; Ronconi et al., 2019; Verza
et al., 2019; Contarini et al., 2019, 2021, 2022a), is defined as:

$\frac{\mathrm{d}n}{\mathrm{d}\ln R}\biggr{|}_{\rm
Vdn}=\frac{\mathrm{d}n}{\mathrm{d}\ln R}\Biggr{|}_{\rm L}\frac{V(r_{\rm
L})}{V(R)}\frac{\mathrm{d}\ln R_{\rm L}}{\mathrm{d}\ln R}\textrm{ .}$ (12)

In the VSF model, voids are defined as non-overlapping spheres embedding an
internal underdensity contrast $\delta_{\rm v,NL}^{\rm DM}$ relative to the
unbiased tracer distribution (DM particles).

Figure 1: The spatial distribution of simulated galaxies (blue dots), galaxy
clusters (yellow dots) and voids (white circles) in the snapshots considered
in this work, extracted from the Magneticum simulations. The points
representing galaxy clusters have dimensions proportional to their mass, while
more massive galaxies have lighter colours. Void circles have have a radius
corresponding to their size. The different panels refer to redshifts z=1,
0.72, 0.52, 0.2, from left to right. Here the cosmic objects considered follow
the selections adopted in this work, i.e. $M_{*}\geq 10^{10}\ h^{-1}\
\rm{M}_{\odot}$ for galaxies, $M_{500c}\geq 10^{14}\ h^{-1}\rm{M}_{\odot}$ for
clusters and $R_{\rm eff}\geq 3.6\ \lambda_{\rm mgs}(z)$ for voids. The axes
are expressed in comoving coordinates, considering a $60$ $h^{-1}\ \rm Mpc$
thick slice in the third coordinate.

In this work we consider voids identified in the distribution of simulated
galaxies, which are biased tracers of the matter density field. In this case,
the radii of the spherical voids predicted by the Vdn model have to be re-
scaled in order to match the radii of voids identified in the nonlinear biased
tracer distribution (Jennings et al., 2013; Ronconi et al., 2019; Contarini et
al., 2019). In our analysis, this is achieved by applying the cleaning
algorithm developed by Ronconi & Marulli (2017) (see Section 3.1, for further
details). The matching of the radii is obtained by re-scaling the model
underdensity threshold $\delta_{\rm v,NL}^{\rm DM}$ by means of a bias factor
$b$:

$\delta_{\rm v,NL}^{\rm tr}=b\delta_{\rm v,NL}^{\rm DM}\rm,$ (13)

where $\delta_{\rm v,NL}^{\rm tr}$ is the tracer density contrast imposed to
be embedded by voids in the biased field and selected during the void cleaning
procedure. This threshold can be chosen from a wide range of negative values
(see Contarini et al., 2019, for further details), as long as the
corresponding linear density contrast is inserted in the Vdn model, as we will
describe shortly. Following the studies of Pollina et al. (2017); Pollina et
al. (2019), Contarini et al. (2019, 2021) showed that the bias $b$ does not
coincide with the one computed on large scales, $b_{\text{eff}}$, measurable
e.g. from the two-point correlation function (2PCF) of the tracer density
field, demonstrating that the bias measured inside cosmic voids is linearly
related to $b_{\text{eff}}$:

$\mathcal{F}(b_{\text{eff}})=B_{\text{slope}}b_{\text{eff}}+B_{\text{off}}\textrm{
,}$ (14)

with $B_{\text{slope}}$ and $B_{\text{off}}$ being the slope and the offset of
this relation, respectively. Consequently, the model nonlinear threshold
$\delta_{\rm v,DM}^{\rm NL}$ can be computed as:

$\delta_{\rm v,NL}^{\rm DM}=\frac{\delta_{\rm v,NL}^{\rm
tr}}{\mathcal{F}(b_{\text{eff}})}\textrm{ .}$ (15)

Moreover, Contarini et al. (2021) found a negligible dependence of this
calibration relation on the considered cosmological model, that is a
fundamental prerequisite to exploit the VSF as a cosmological probe.

Finally, the theoretical VSF model requires the corresponding linear
underdensity threshold. Having $\delta_{\rm v,NL}^{\rm DM}$, one can convert
it into its linear value via the following relation (Bernardeau, 1994):

$\delta_{\rm v,L}^{\rm DM}=\mathcal{C}[1-(1+\delta_{\rm v,NL}^{\rm
DM})^{-1/\mathcal{C}}]\text{ ,}$ (16)

with $\mathcal{C}=1.594$, the latter being exact for cosmologies with
$\Lambda=0$, and very precise also for other values of $\Lambda$, especially
when applied to the negative density contrasts.

##  3 The Magneticum simulations

In this work we make use of simulated galaxy and galaxy cluster catalogues
extracted from the Magneticum Pathfinder111http://www.magneticum.org/ (Dolag
et al., in preparation), a large set of cosmological, hydrodynamical
simulations, having box volumes going from $(12\ h^{-1}\ \text{Mpc})^{3}$
(Box6) to $(2688\ h^{-1}\ \text{Mpc})^{3}$ (Box0). These boxes have different
resolutions, and thus allow us to study the formation and evolution of both
large-scale structures as well as phenomena that occur on smaller scales, from
the motion of galaxies to the physics of gas inside them. The Magneticum
simulations have been run with the parallel code P-GADGET3, which is an
updated version of the TreeSPH GADGET-2 code presented in Springel (2005). In
this code, the gravitational forces are computed through a TreePM algorithm,
while the hydrodynamics is modelled through an updated SPH algorithm (Springel
& Hernquist, 2002; Beck et al., 2016), from which it is possible to properly
track the gas turbulence (Dolag et al., 2005).

Furthermore the main baryonic physics phenomena are implemented in these
hydrodynamic simulations, following the methods presented in Springel &
Hernquist (2003), such as the cooling of gas, the star formation and
supernovae feedback. In addition, black hole and active galactic nuclei
feedbacks are included (Di Matteo et al., 2008; Fabjan et al., 2010;
Hirschmann et al., 2014), as well as thermal conduction (Dolag et al., 2004),
stellar population and chemical enrichment models (Tornatore et al., 2003,
2007). Following the latter, metals are produced by SNII, SNIa, and by low and
intermediate mass stars in the asymptotic giant branch. The simulation assumes
a Chabrier initial stellar mass function (Chabrier, 2003) and takes into
account also the different release amount of metals and energy depending on
the stellar mass (Padovani & Matteucci, 1993). A more detailed description can
be found in Hirschmann et al. (2014) and Teklu et al. (2015).

The Magneticum simulations have been carried out with the seven-year Wilkinson
Microwave Anistropy Probe (WMAP7) cosmology (Komatsu et al., 2011), which we
will consider as fiducial. In this scenario the main cosmological parameters
are set to $\Omega_{\mathrm{m}}=0.272$, $\Omega_{\Lambda}=0.728$,
$\sigma_{8}=0.809$, $H_{0}=70.4\ \text{km}\ \text{s}^{-1}\text{Mpc}^{-1}$, and
$n_{\mathrm{s}}=0.963$. For our analysis we exploit the Box1a (for details see
Dolag et al., 2016), which is a large simulation that follows the evolution of
$2\cdot 1526^{3}$ particles in a comoving volume of $(896\ h^{-1}\
\text{Mpc})^{3}$, with which it is possible to make a detailed statistical
analysis of galaxy clusters and cosmic voids. This snapshot is built
considering a mass resolution of DM and gas particles of $1.3\times 10^{10}\
h^{-1}M_{\odot}$ and $2.6\times 10^{9}\ h^{-1}M_{\odot}$, respectively. The
simulation outputs, from which redshift-space cluster and void catalogues have
been constructed, are selected at four redshifts, $z=0.2,0.52,0.72,1$.

###  3.1 Data preparation

We construct redshift-space mock catalogues of galaxies and galaxy clusters
starting from their positions in the real space by following the same
methodology used in Bianchi et al. (2012) and Marulli et al. (2011); Marulli
et al. (2012a); Marulli et al. (2012b, 2017), which we briefly outline in the
following. We set a virtual observer at redshift $z=0$, and place the centre
of the different simulation snapshots at a comoving distance $D_{\rm c}$,
corresponding to their redshift:

$D_{\rm c}(z)=c\int_{0}^{z}\frac{\mathrm{d}z^{\prime}_{\rm
c}}{H(z^{\prime}_{\rm c})}\text{ ,}$ (17)

where $c$ is the speed of light, $z_{\rm c}$ is the cosmological redshift due
to the Hubble recession velocity and $H(z)$ is the Hubble parameter. The
latter, for a $\Lambda$CDM cosmology and considering the contribution of
radiation as negligible, becomes:

$H(z)=H_{0}\bigl{[}\Omega_{\rm
m}(1+z)^{3}+\Omega_{\Lambda}\bigr{]}^{1/2}\text{ .}$ (18)

We then transform the comoving coordinates of real-space mock sources into
angular positions and observed redshifts. For each galaxy in a given
catalogue, the latter is computed from $z_{\rm c}$ as:

$z_{\rm obs}=z_{\rm c}+\frac{v_{\parallel}}{c}(1+z_{\rm c})\text{ ,}$ (19)

with $v_{\parallel}$ being the galaxy centre of mass velocity projected along
the line-of-sight. As in Marulli et al. (2012b, 2017), we do not take into
account any systematic error in the observed redshift, neither we consider
geometrical distortions during the redshift-space mock catalogues generation,
i.e. during this procedure we assume the correct, true underlying cosmology of
the simulations.

In the Magneticum simulations, haloes are identified with a Friend-of-Friend
(FoF) algorithm (Davis et al., 1985), with a linking length $b=0.16$ applied
to the distribution of DM particles (Dolag et al., 2009) while galaxy clusters
are detected, starting from halo centres, as spherical overdensities. In our
case, we consider spheres containing a linear overdensity
$\Delta=500\rho_{c}$, where $\rho_{\rm c}\equiv 3H^{2}/8\pi G$ is the critical
density of the Universe at the considered redshift. We then select only
clusters having mass $M_{\text{500c}}\geq 10^{14}\ h^{-1}\rm{M}_{\odot}$,
choosing this mass selection in order to mimic the range of masses for galaxy
clusters in real current and future surveys, like the extended Roentgen Survey
with an Imaging Telescope Array (eROSITA) (Merloni et al., 2012). Also, for
smaller masses we verified the presence of incompleteness in the cluster
counts. We do not focus on any observable-mass relation in this work, and we
consider cluster masses as a direct observable.

We underline that the adopted halo mass definition is not unique and that it
is possible to switch to a different one under the assumption of a mean halo
profile. For example, the mass cut off $M_{\text{500c}}\geq 10^{14}\ h^{-1}\
M_{\odot}$ would approximately translate into $M_{\text{200c}}\geq 1.4\times
10^{14}\ h^{-1}\ M_{\odot}$ by considering a Navarro-Frenk-White profile with
a concentration $\rm c=4$ (for further details see e.g. Despali et al., 2016).
We therefore do not expect the halo mass definition to have a significant
impact on the results presented in our analysis.

Voids are instead detected in the redshift-space distribution of galaxies. In
the Magneticum simulations galaxies are identified using a modified version of
subfind, an algorithm aimed at extracting bound structures from a large
$N$-body simulation characterised by the presence of gas and star particles
(Springel et al., 2001; Dolag et al., 2009). We apply then a galaxy mass
selection of $M_{*}\geq 10^{10}\ h^{-1}\ \rm{M}_{\odot}$ following the choice
of Marulli et al. (2017), who analysed the same cosmological simulations. We
underline that there is not a direct link or correspondence between the galaxy
cluster and galaxy mass cuts examined in this work.

In each simulated box, redshift-space voids are identified with the public
Void IDentification and Examination
toolkit222https://bitbucket.org/cosmicvoids/vide_public (VIDE, Sutter et al.,
2015), which is based on the parameter-free void finding algorithm ZOnes
Bordering On Voidness code (ZOBOV, Neyrinck, 2008). The latter finds voids in
a 3D distribution of tracer particles, without introducing any free parameter
or assumptions about the void shape.

The extracted void catalogues were then cleaned with the algorithm implemented
inside the
CosmoBolognaLib333https://gitlab.com/federicomarulli/CosmoBolognaLib
C++/Python libraries (Marulli et al., 2016) and exploited in recent works
(Ronconi & Marulli, 2017; Ronconi et al., 2019; Contarini et al., 2019, 2021,
2022a; Contarini et al., 2022b, c). As already mentioned in Section 2.2, the
cleaning procedure is necessary because the VSF models described above
consider voids as underdense, spherical, non-overlapping regions with a
specific internal density contrast. Hence the aim of adopting this procedure
is to prepare catalogues matching the theoretical definition used to develop
the VSF model. The cleaning algorithm can be briefly described in three main
steps: it removes non-relevant objects, i.e. voids with a radius outside of a
given range and/or with a central contrast that is too high; it rescales the
effective void radii, $R_{\rm eff}$, in order to make the voids embedding a
specific density contrast; finally it checks for overlapping voids, eventually
rejecting the ones with the highest central density. The impact of the
cleaning procedure on the void abundances has been extensively tested in
Ronconi et al. (2019), which applied it to a set of cosmological simulations
at different redshifts having different spatial resolutions and box side
lengths.

Once the cleaning algorithm is applied, we select only voids having effective
radii $R_{\rm eff}\geq 3.6\ \lambda_{\rm mgs}(z)$, with $\lambda_{\rm mgs}(z)$
being the mean galaxy separation of the catalogues from which voids are
detected. This quantity corresponds to $\lambda_{\rm mgs}\simeq 6.6\ h^{-1}\rm
Mpc$ at $z=0.2$ and increases up to $\lambda_{\rm mgs}\simeq 8.9\ h^{-1}\rm
Mpc$ at $z=1$ because of the lower number of galaxies at higher redshifts. We
consider this extremely conservative radius cut in order to ensure reliable
constraints on cosmological parameters (see Section 5.2). Indeed, with this
selection we avoid considering the spatial scales characterised by
incompleteness, i.e. those voids affected by the limited resolution of the
simulation. Indeed, we verified that a less stringent radius selection would
lead to cosmological test constraints slightly in disagreement with the true
values of the Magneticum simulation parameters. We refer the reader to
previous works (i.e. Jennings et al., 2013; Ronconi et al., 2019; Contarini et
al., 2019) for further details on the impact of the spatial resolution on void
number counts.

Figure 1 shows the evolution of the large-scale distribution of galaxies (blue
dots), galaxy clusters (yellow dots) and voids (white circles) inside slices
of thickness 60 $h^{-1}\rm Mpc$, extracted from the Magneticum snapshots. We
also report simulated voids following the selection procedure described above,
i.e. by applying cleaning and selection procedures. As one can note from this
figure, the large-scale structure clustering increases with cosmic time,
giving rise to more massive collapsed objects (i.e. galaxies and galaxy
clusters) and, at the same time, to larger cosmic voids. In this
representation however, some voids seem to be located in regions of galaxy
crowding. This is just a projection effect given by the finite thickness of
the shown slices. Moreover, some underdensity regions do not have an
associated void: this can be explained in terms of the applied void selection,
for which we considered only very large voids (see Section 3.1).

##  4 Methods

###  4.1 From real-space to redshift-space

The theoretical framework outlined in Section 2 is only suitable for a real-
space analysis. Considering redshift-space galaxy and galaxy cluster
catalogues, it is necessary to correct for geometric and dynamic distortions.
The former arise when the fiducial cosmology assumed in the analysis is
different from the true one, while the latter are due to the peculiar motion
of galaxies with respect to the Hubble flow. We will start our discussion by
providing a model correction for the geometric distortions, which must be
taken into account during the Bayesian Markov Chain Monte Carlo (MCMC)
analysis. In fact, as the cosmological parameters change at each step of the
MCMC, the mapping between redshift and comoving distance will be different, as
well as the simulated box volume. The former is known as the Alcock–Paczynski
effect (AP, Alcock & Paczynski, 1979) and influences the observed shape of
large-scale structures, making them appearing elongated along the line of
sight. Instead, the variation of the simulated volume affects the predicted
number density of both clusters and voids.

Regarding clusters, we will include the effect of the varying simulated volume
only, since the AP effect has a negligible impact on their masses. Viceversa,
voids must be considered as spatially extended objects, hence the geometrical
distortions must be taken into account when modelling their observed sizes.

Voids identified assuming the true cosmological model can be modelled as
spherical objects having radius a $R$. However, it can be shown that assuming
a fiducial cosmology different from the true one444Measurements computed in
the fiducial cosmology are indicated here with primed symbols., the void
radius $R^{\prime}$ can be written as (see Contarini et al., 2022a, and
references therein):

$R=q_{\parallel}^{1/3}q_{\bot}^{1/3}R^{\prime}\rm,$ (20)

where $q_{\parallel}$ and $q_{\bot}$ are defined via the following relations:

$\displaystyle r_{\parallel}^{\prime}$
$\displaystyle=\frac{H(z)}{H^{\prime}(z)}r_{\parallel}=q_{\parallel}^{-1}r_{\parallel}\rm,$
(21) $\displaystyle r_{\bot}^{\prime}$ $\displaystyle=\frac{D_{\rm
A}(z)}{D_{\rm A}^{\prime}(z)}r_{\bot}=q_{\bot}^{-1}r_{\bot}\rm,$

i.e. as the ratio between the Hubble parameter, $H(z)$, and the ratio between
the comoving angular-distances, $D_{\rm A}(z)$, computed assuming the two
cosmologies, respectively. In Eq. (21), $r_{\parallel}$ and $r_{\bot}$ are the
parallel and perpendicular projections, with respect to the line of sight, of
the distance vector connecting two points at redshift $z$. In applying this
correction, we consider the mean redshift of the centers of voids included in
each redshift interval. This approximation is accurate enough since the
selected samples of voids cover small redshift ranges.

As already mentioned, the other effect concerns the variation of the simulated
box volume. This can be modelled as:

$V_{\rm box}(z)=L_{\rm box}^{2}\cdot\Delta D_{\rm C}(z)\rm,$ (22)

where $\Delta D_{\rm C}(z)=D_{\rm C}(z_{\rm max})-D_{\rm C}(z_{\rm min})$ is
the difference between the comoving distance of the furthest object and the
nearest object in the simulated box. Here the redshift value depends on the
simulated snapshot being considered.

With these corrections, the number counts of clusters in a given mass bin
$\Delta M\equiv M_{2}-M_{1}$ can be obtained as:

$N_{\rm cl}(\Delta M;V_{\rm box})=V_{\rm
box}\int_{M_{1}}^{M_{2}}\text{d}M\frac{\mathrm{d}n}{\mathrm{d}M}\rm,$ (23)

where $\mathrm{d}n/\mathrm{d}M$ is the HMF given by Eq. (6), and $V_{\rm box}$
is computed as in Eq. (22). Accordingly, the number counts of voids are
computed from the VSF model of Eq. (12) as:

$N_{\rm v}(\Delta R;V_{\rm box})=V_{\rm
box}\int_{R_{1}}^{R_{2}}\frac{\text{d}R}{R}\frac{\mathrm{d}n}{\mathrm{d}\ln
R}\biggr{|}_{\rm Vdn}\rm,$ (24)

where $\Delta R\equiv R_{2}-R_{1}$ is a given void radius bin, with $R$
computed from Eq. (20), i.e. taking into account the AP effect.

The inclusion of the geometric distortions in the HMF and VSF models, via the
simulated volume change and the AP effect (for cosmic voids only), introduces
a further dependence on the cosmological model. In particular, the total
matter density parameter $\Omega_{\mathrm{m}}$ enters in the computation of
the Hubble parameter $H(z)$ and in the comoving angular-distance
$D_{\mathrm{A}}$ (see Eq. 21), as well as in the comoving distance
$D_{\mathrm{C}}$ (see Eq. 22).

Finally, regarding dynamic distortions, we do not consider their impact on the
HMF since the effect on the observed cluster masses can be considered
negligible. Instead, it is fundamental to include the redshift-space
distortions in the VSF model because the observed void radii undergoes a
significant enlargement due to the coherent outflow of tracers from the void
centres (Pisani et al., 2015a). This effect can be encapsulated in the
nuisance parameters, $B_{\mathrm{slope}}$ and $B_{\mathrm{off}}$, of the
extended Vdn theory (Contarini et al., 2022a). In our analysis, this will be
performed in Section 5.1, during the calibration of the VSF with the measured
redshift-space void number counts.

###  4.2 Bayesian analysis

Having extracted our data sets $\mathcal{D}$, consisting of cluster and void
counts measured from the simulations, we follow a Bayesian approach to sample
the posterior distribution of two fundamental cosmological parameters, i.e.
the total matter density parameter, $\Omega_{\rm m}$, and the $z=0$ amplitude
of the matter power spectrum, $\sigma_{8}$. Let $\Theta=(\Omega_{\rm
m},\sigma_{8},\textbf{q})$ denote the parameter vector, with q being the set
of the internal HMF and VSF model parameters. Then the posterior distribution
of $\Theta$, given a cosmological model, can be computed from the Bayes’s
theorem:

$\mathcal{P}(\Theta|\mathcal{D})\propto\mathcal{L}(\mathcal{D}|\Theta)p(\Theta)\rm,$
(25)

where $\mathcal{L}(\mathcal{D}|\Theta)$ is the likelihood function and
$p(\Theta)$ the priors considered for the parameters.

We then make use of MCMC methods for the posterior sampling, taking into
account the correlation between different data sets where needed (see Section
4.3). In the general case, i.e. where the data set covariance matrix
$\mathcal{C}$ is non-diagonal, the following Gaussian likelihood is
considered:

$\mathcal{L}(\mathcal{D}|\Theta)=\prod_{i,j}\frac{1}{(2\pi)^{n/2}\sqrt{|\mathcal{C}|}}\exp{\Bigl{(}-\frac{1}{2}\chi_{ij}^{2}\Bigr{)}}\textrm{
,}$ (26)

where

$\chi_{ij}^{2}=(N_{ij}-\mu_{ij})^{\text{T}}\mathcal{C}^{-1}(N_{ij}-\mu_{ij})\textrm{
.}$ (27)

Here $n$ is the number of items in the data set, $N_{ij}\equiv
N(\xi_{i},z_{j}|\mathcal{D})$ is the number counts for the given observable
$\xi$ in the $i$-th bin (mass or radius) at redshift $z_{j}$, and
$\mu_{ij}\equiv N(\xi_{i},z_{j}|\Theta)$ is the expectation value in the
cosmological model considered. Note that the covariance matrix of a given
probe (i.e. HMF or VSF, in our case),
$\mathcal{C}\equiv\mathcal{C}(\xi_{i},z_{j})$, depends both on the measured
statistic in the given bin, $\xi_{i}$, and on the mean redshift of the bin,
$z_{j}$. In this work we exploit the combination between different posterior
distributions, obtained separately from our cosmological probes. This can be
achieved via the product of their likelihood functions, or more generally
their posterior product.

###  4.3 Correlations between cluster and void counts

To check for possible correlations between cluster and void counts, we compute
the cross-covariance matrix of the two probes considered, employing the
jackknife sampling technique (Parr & Schucany, 1980; Efron, 1982; Escoffier et
al., 2016). The latter is used to generate new mock samples that are
representative of an underlying population, i.e. the analysed catalogues in
this work, whose statistics we want to estimate. This technique works by
sequentially removing one sub-sample from the entire data set, computing the
desired statistics on the remaining samples.

Figure 2: Jackknife cross-correlation matrix for galaxy cluster and void
counts, for the mass and radius selections considered in the paper. The
colorbar on the right represents the Pearson correlation coefficient
associated to cluster mass and void radius bins. Data sets at different
redshifts are separated by grey solid lines, while mass bins and radius bins
are divided by solid black lines.

For each redshift considered, we divide the cluster and void catalogues into
$N_{\rm sub}=125$ equally sized sub-regions. A value for the number counts of
galaxy clusters $N_{\alpha}$ and cosmic voids $N_{\beta}$ can be associated
with $\alpha\sim(i_{M},i_{z})$, a generic bin in mass and redshift, and
$\beta\sim(i_{R},i_{z})$, a generic bin in radius and redshift. The jackknife
covariance element between two generic bins $\alpha$ and $\beta$ can then be
estimated as (see e.g. Efron, 1982):

$\mathcal{C}^{\text{jk}}_{\alpha\beta}=\frac{N_{\text{sub}}}{N_{\text{sub}}-1}\sum_{i=1}^{N_{\text{sub}}}\bigl{(}N^{\text{jk}}_{\alpha}(i)-\bar{N}^{\text{jk}}_{\alpha}\bigr{)}\bigl{(}N^{\text{jk}}_{\beta}(i)-\bar{N}^{\text{jk}}_{\beta}\bigr{)}\textrm{,}$
(28)

where:

$\bar{N}^{\text{jk}}_{\alpha}=\frac{1}{N_{\text{sub}}}\sum_{i=1}^{N_{\text{sub}}}N^{\text{jk}}_{\alpha}(i)$
(29)

is the average of the jackknife counts in a specific bin and
$N^{\text{jk}}_{\alpha}(i)=N^{\text{tot}}_{\alpha}-N_{\alpha}(i)$ is a
jackknife sample, defined as the difference between the number counts in the
total volume and the ones in the volume of one sub-sample. In this way, the
cross-covariance matrix between the number counts of galaxy clusters and
cosmic voids is a block matrix, and can be written as:

$\mathcal{C}^{\text{jk}}=\left(\begin{array}[]{c|c}\mathcal{C}_{\rm
HH}&\mathcal{C}_{\rm VH}\\\ \hline\cr\mathcal{C}_{\rm HV}&\mathcal{C}_{\rm
VV}\end{array}\right)\textrm{ ,}$ (30)

where $\mathcal{C}_{\rm HH},\mathcal{C}_{\rm VH},\mathcal{C}_{\rm HV}$ and
$\mathcal{C}_{\rm VV}$ are the covariance matrices between cluster mass bins,
cluster mass and void radius bins (and vice-versa) and void radius bins at all
redshifts ($z=0.2,0.52,0.72,1$), respectively. Since the covariance matrix has
been computed using a limited set of mock catalogues, we include the
statistical corrections suggested by Percival et al. (2014) to improve the
accuracy in the estimation of the model parameters. Additionally, we apply the
prescriptions by Hartlap et al. (2007) to correct for the numerical
uncertainties in the inversion of the covariance matrix, needed in Eq. (27).
We note some level of noise in the computed covariance matrix. We discuss the
stability of our results against the variation of $N_{\rm sub}$ in Appendix A.

In Figure 2 we present the cross-correlation matrix between our two probes,
computed for each bin as the Pearson correlation coefficients:

$\rho_{ij}=\frac{\mathcal{C}_{ij}}{\sqrt{\mathcal{C}_{ii}\times\mathcal{C}_{jj}}}\textrm{
.}$ (31)

We clearly note that the internal correlation matrix of galaxy cluster counts,
i.e. the upper left block matrix in Figure 2 presents non-vanishing off-
diagonal elements. These features are due to the fact that we are analysing
the redshift evolution of the same large-scale structure, so we are
considering the same objects evolving with time. We checked that these
correlation features in fact disappear when estimating again the covariance
matrix for cluster number counts, but resampling from different cosmological
sub-boxes at each redshift. Despite this approach would be more accurate in
reproducing the behaviour of real cosmic objects, it severely reduces the
measured cluster number counts and so would negatively affects the statistical
relevance of the analysis. For this reason we applied this methodology as a
consistency test only. On the contrary, we do not get evident correlation
features in the void count correlation matrix. We attribute this fact to both
the noisier data and the effect of the cleaning procedure, which reduces the
possibility of following the evolution of the same void sample with cosmic
time. Notably, we can appreciate how the level of correlation of the off-
diagonal block matrices $\mathcal{C}_{\rm VH}$ and $\mathcal{C}_{\rm HV}$ is
consistent with zero, confirming the high statistical independence of cluster
and void number counts, as also found in previous works (Sahlén et al., 2016;
Bayer et al., 2021; Kreisch et al., 2022).

##  5 Results

Before comparing the simulated data sets extracted from the Magneticum
simulations with the theoretical models described in Section 2, a self-
calibration (i.e. on the same simulation outputs, considering the fiducial
cosmology) of the latter is necessary, because of the differences in the halo
definitions between the simulations employed here and those assumed by the
theoretical HMF model, as already explained in Section 2.1. In this case, the
calibration will allow us to determine the posterior distribution for the
coefficients that appears in Eq. (7). Consequently, we will marginalise over
the latter in order to find the posterior distributions for the cosmological
parameters $\Omega_{\rm m}$ and $\sigma_{8}$ (see Section 5.2).

Figure 3: Posterior probability distribution of the coefficients of the
fitting relation parameters ($a,A$) reported in Eq. (7), for a cluster mass
selection of $M_{500c}\geq 10^{14}\ h^{-1}\ \rm{M}_{\odot}$. Dark and light
areas show the $68\%$ and $95\%$ confidence regions, respectively. On the top
of each columns we report the projected 1D marginal posterior distributions,
with the associated maximum posterior values and relative $1\sigma$
uncertainties. Figure 4: Posterior probability distribution of
$B_{\text{slope}}$ and $B_{\text{off}}$, coefficients of the
$\mathcal{F}(b_{\text{eff}})$ relation, reported in Eq. (14). Dark and light
areas show the 68% and 95% confidence regions, while at the top of each column
we report the projected 1D distributions, with the values associated to the
maximum of the posterior and the relative $1\sigma$ uncertainties.

Moreover, also the VSF model has to be calibrated (Contarini et al., 2019,
2021, 2022a). We make use of the Vdn model described in Section 2.2, modifying
properly the underdensity threshold that enters in the model, as reported in
Eq. (14). Specifically, the calibration we perform on the VSF model is to
assess the posterior distribution of the free parameters of the function
$\mathcal{F}(b_{\mathrm{eff}})$, i.e. $B_{\mathrm{slope}}$ and
$B_{\mathrm{off}}$.

For both the calibrations and constraints on cosmological parameters, we fit
the joint data set at all the redshift considered with the model presented in
Section 2. In particular, we compute the constraints by concatenating the data
sets concerning the two probes analysed at different redshifts into a unique
vector, constructing the likelihood with the correct covariance matrix for the
given cosmological probe.

###  5.1 Model self-calibrations

The calibration performed in Despali et al. (2016) considered DM-only
cosmological simulations, while the Magneticum simulations also follow the
baryonic component. Moreover, Despali et al. (2016) defined the galaxy
clusters with a halo finder method different from the one used in this work.
For these reasons, we have to re-calibrate the parameters of the HMF model, by
assuming the functional relation given by Eq. (7). This translates into
finding other coefficients for the relations $a=a(x),p=p(x)$ and $A=A(x)$. We
chose to consider $p$ as fixed to the value provided by Despali et al. (2016)
in the most general fitting case, that is $p=0.2536$. We made this choice
because $p$ significantly affects only the low mass end of the HMF, and we
cannot precisely constrain it due to the mass selection applied and the
simulation mass resolution.

We expressed the general form of the fitting formula Eq. (7) as:

$\displaystyle a$ $\displaystyle=a_{1}x^{2}+a_{2}x+a_{3}$ (32) $\displaystyle
A$ $\displaystyle=-A_{1}x+A_{2}$ $\displaystyle p$
$\displaystyle=0.2536\textrm{ ,}$

and we calibrate the new relations for the model parameters by performing a
Bayesian MCMC analysis on the simulated HMF, considering flat priors for the
five free parameters ($a_{1},a_{2},a_{3},A_{1},A_{2}$) and sampling their
posterior distributions. Since the correlation matrix associated to the HMF
(upper left block matrix in Figure 2) is non-diagonal, during the calibration
procedure we consider a Gaussian likelihood as in Eq. (26).

The result of the calibration is presented in Figure 3, which shows the
posterior distribution of the coefficients of Eq. (32). We find strong
correlations between the different calibration parameters. Moreover the
marginalised 1D posterior distributions present an approximately Gaussian
shape, with a low degree of skewness, the latter being more prominent for the
$a_{i}$ parameters. While $A_{1}$ and $A_{2}$ result statistically in
agreement with the values obtained from the calibration performed in Despali
et al. (2016) (see Eq. 7), the values for the $a_{i}$ coefficients are in
$>2\sigma$ disagreement with respect to their recovered values. This is an
expected outcome due to the different adopted halo definitions (see Despali et
al., 2016, for further details), and also to the effect of the baryonic
physics (see e.g. Bocquet et al., 2016; Ragagnin et al., 2021), which is
included in the simulations analysed in this work.

Figure 5: Measured number counts of galaxy clusters (red dots), from the
Magneticum (Box1a) at redshifts $z=0.2,0.52,0.72,1$, for ${M_{500\rm c}\geq
10^{14}\ h^{-1}\rm{M}_{\odot}}$. Red solid lines represents the HMF model from
Despali et al. (2016) after the parameter re-calibration. Lower sub-panels
report the residuals of the cluster counts, computed as the difference of the
measured data and the HMF model, in units of the data errors. The light grey
bands represent the $1\sigma$ intervals for the residuals. Figure 6: Measured
number counts of voids (blue dots) identified in the distribution of galaxies,
from the Magneticum simulation (Box1a) at redshifts $z=0.2,0.52,0.72,1$. The
blue lines represent the corresponding predictions of the Vdn model, extended
and calibrated according to the prescriptions provided in Section 2.2. The
lower sub-panels report the residuals of the void counts, computed as the
difference of the measured data and the VSF model, in units of the data
errors. The light blue bands represent the $1\sigma$ intervals for the
residuals.

Similarly to what we have done for the HMF parameters, here the VSF model
calibration consists in finding new values for the parameters
$B_{\text{slope}}$ and $B_{\text{off}}$ of Eq. (14), by sampling their
posterior distributions via a MCMC analysis, with the cosmological parameters
fixed to the simulation true values.

In order to do this, we first compute the effective large-scale bias,
$b_{\text{eff}}$, following the method delined in Marulli et al. (2013) and
Marulli et al. (2018). In particular, we construct a 10-time denser random
catalogue with the same properties of the Magneticum simulation boxes and
measure the 2PCF of real-space galaxies making use of the Landy-Szalay
estimator (Landy & Szalay, 1993). We then sample the posterior distribution of
the free parameter of the model, $b_{\text{eff}}$, following the same method
described in Contarini et al. (2019, Appendix A). We obtain $b_{\rm
eff}=1.321\pm 0.007,1.549\pm 0.008,1.714\pm 0.009,1.98\pm 0.01$ for the four
redshifts considered in our analysis, modelling the galaxy 2PCF in the range
$[20-40]\ h^{-1}\rm Mpc$. Then we calibrate the VSF model by fitting the
measured void number counts. In particular we consider flat priors for $B_{\rm
slope}$ and $B_{\rm off}$, marginalising over the values of $b_{\rm eff}$ by
assuming 1D Gaussian priors centred on the computed mean values and with
standard deviation equal to their reported errors. As for cluster counts, we
considered a Gaussian likelihood (Eq. 26) with the covariance matrix computed
in Section 4.3.

We present the result of this calibration in Figure 4, where we show the
posterior distribution of $B_{\rm slope}$ and $B_{\rm off}$, which have both
approximately Gaussian shape. Also in this case we found a strong correlation
between the parameters, with degeneracy direction consistent with recent works
(Contarini et al., 2021, 2022a). We find a discrepancy with the results found
in Contarini et al. (2022a), i.e. $B_{\rm slope}=0.96\pm 0.04$ and $B_{\rm
off}=0.44\pm 0.07$. This is due to the different mass tracer type and
selection used by these authors: in Contarini et al. (2022a) voids are traced
by H$\alpha$ galaxies, simulated by means of a halo occupation distribution
algorithm and selected according to a magnitude cut (see Contarini et al.,
2022a, for further details). In this work we analyse instead voids identified
in hydrodynamical mock samples of galaxies, selected with a threshold in
stellar mass.

###  5.2 Cosmological test constraints

We can now present the measurements of the number counts of both galaxy
clusters and cosmic voids extracted from redshift-space catalogues, and
compare these with the theoretical predictions given by the calibrated models.
We rejected voids that are too close to the simulation boundaries, considering
an edge of $15\ h^{-1}$ Mpc, since the cleaning algorithm cannot rescale them
accurately. We present the measured number counts of galaxy clusters and
cosmic voids in Figures 5 and 6, respectively, for the four redshift bins
considered. The bars represent the Poissonian errors related to the data while
the $\pm 1\sigma$ uncertainty of theoretical models are shown as shaded
regions. From the residuals, computed as the difference between the measured
data points and the theoretical models in units of the data error, one can
appreciate how the measured number counts are accurately fitted by the
theoretical models. Indeed, residuals are in all cases within $1\sigma$.

Once measured the HMF and VSF models from the simulated data sets, we extract
test constraints on the cosmological parameters $\Omega_{\rm m}$ and
$\sigma_{8}$. As already explained, we consider cluster and void abundances as
a function of the mass and radius, respectively, as statistically independent,
meaning that we do not take into account the cross-covariance matrix between
the two probes. We perform a Bayesian MCMC analysis, considering flat priors
on $\Omega_{\rm m}$ and $\sigma_{8}$ with ranges $[0.1-0.5]$ and $[0.5-1.0]$,
respectively, and marginalising over the non-cosmological parameters, i.e.
considering them as nuisances. Moreover, during the variation of $\Omega_{\rm
m}$ in the MCMC, we imposed the flat spatial geometry of the Universe by
rescaling the value of the CDM component as ${\Omega_{\rm cdm}=\Omega_{\rm
m}-\Omega_{\rm b}}$ in order to hold true the condition $\Omega_{\rm
m}+\Omega_{\Lambda}=1$.

The correlation between nuisance parameters of the models, i.e.
$(a_{i},A_{i})$ for the HMF and $(B_{\rm slope},B_{\rm off})$ for the VSF, is
taken into account by considering a multivariate normal distribution as
corresponding prior, centred on the parameter average values obtained from the
calibrations (Sections 5.1). As already mentioned, we use the jackknife
covariance matrix for both cluster and void counts. Then, we progressively
stack the cluster and void abundance data sets relative to different
simulation snapshots, sampling the joint posterior distribution of
$\Omega_{\rm m}$ and $\sigma_{8}$ for increasing redshift ranges. This
procedure is aimed at assessing the constraining power on the growth of cosmic
structures deriving from the synergy of cluster and voids number counts at
different cosmic times.

Figure 7 shows the resulting contours for the HMF model considering redshift-
space clusters with masses $M_{500\rm c}\geq 10^{14}\ h^{-1}\rm{M}_{\odot}$,
for the VSF model with voids having radius $R_{\rm eff}\geq 3.6\ \lambda_{\rm
mgs}(z)$, and for their combination, considering a gradually increasing
maximum redshift up to $z=1$. In all cases, the HMF and VSF contours present
different orientations in the $\Omega_{\rm m}$-$\sigma_{8}$ parameter space.
In particular they are almost perpendicular, which is a powerful feature when
we consider their combination. As expected, the precision of the constraints
increases by adding the cosmological information of different Universe epochs.
We underline indeed how the combined confidence contour shrinks by going from
$z\leq 0.2$ to $z\leq 1$, remaining however perfectly in agreement with the
true cosmological parameters of the simulations. We also notice that the VSF
gains more constraining power with the extension of the redshift range with
respect to the HMF. Thanks to the negligible correlation between the void
samples at different redshifts, the combination of the corresponding VSF
constraints results indeed very effective. For the HMF, instead, adding
cosmological information from higher redshifts is less advantageous mainly
because of the correlation between the galaxy clusters belonging to the
different simulation snapshots: as already explained, this is taken into
account by including in the modelling of cluster counts the full-covariance
matrix (see Section 4.3) and results in a milder increase of the constraining
power since it derives from the joint analysis of dependent data sets.

The final combined constraints on $\Omega_{\rm m}$ and $\sigma_{8}$, obtained
by considering all the data sets up to $z=1$, are the following:

$\displaystyle\Omega_{\rm m}=0.270\pm 0.007$ (33)
$\displaystyle\sigma_{8}=0.809\pm 0.005\rm,$

which, as well as those coming from the single probes, are well centered on
the truth values of the simulation, $\Omega_{\rm m}=0.272$ and
$\sigma_{8}=0.809$. We notice that, despite the maximum of the posterior
distribution deriving from the VSF is slightly shifted towards lower values in
the $\Omega_{\rm m}$-$\sigma_{8}$ space, the true simulation cosmological
parameters are widely included in the the corresponding 68$\%$ confidence
region.

To have an easier comparison between the combination constraints and the ones
from HMF and VSF alone, we derive 1D constraints also for
$S_{8}\equiv\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$, which are reported in Table
1. First of all we note how the truth value $S_{8}=0.77$ is well recovered
within the 68$\%$ confidence region for all the redshift ranges considered in
the table. Then, in order to evaluate the contribution of the VSF in the
performed probe combination, we show how much the precision on $S_{8}$
increases with respect to that computed from the HMF alone. This percentage
improvement, denoted here as $\mathcal{I}_{\rm HMF}$, is computed as the ratio
between the relative errors on $S_{8}$ from HMF and from the combined
analysis. We find an improvement on the constraining power of about $60\%$
compared to that obtained with the HMF alone, almost independent of the
maximum redshift considered in the analysis.

This result highlights the impressive contribution of the VSF in enhancing the
constraining power derived with standard probes like the HMF, moreover it
underlines the fundamental importance of exploiting complementary probes in
future redshift surveys to achieve precise and independent cosmological
constraints.

Table 1: Test constraints up to a given redshift (first column) on the $S_{8}$ parameter. The reported errors represents the $68\%$ confidence regions on $S_{8}$. In the last column $\mathcal{I}_{\rm HMF}$ represents the improvement on the constraining power achieved by the combination of HMF and VSF compared to the HMF alone, computed as the ratio between the relative errors on $S_{8}$ from HMF and from HMF+VSF. Redshift | $S_{8}(\rm HMF)$ | $S_{8}(\rm VSF)$ | $S_{8}(\rm HMF+VSF)$ | $\mathcal{I}_{\rm HMF}$
---|---|---|---|---
$z\leq 0.2$ | $0.763^{+0.019}_{-0.017}$ | $0.75^{+0.10}_{-0.09}$ | $0.77\pm 0.01$ | $67\%$
$z\leq 0.52$ | $0.760^{+0.017}_{-0.014}$ | $0.75\pm 0.07$ | $0.77\pm 0.01$ | $65\%$
$z\leq 0.72$ | $0.760^{+0.016}_{-0.015}$ | $0.74\pm 0.06$ | $0.767^{+0.009}_{-0.01}$ | $63\%$
$z\leq 1$ | $0.765^{+0.017}_{-0.013}$ | $0.76\pm 0.06$ | $0.77\pm 0.01$ | $67\%$
Figure 7: Posterior probability distribution on $\Omega_{\rm m}$ and
$\sigma_{8}$ from VSF (blue), HMF (red) and their combination (gray), for
clusters and voids extracted from Magneticum simulations, considering their
joint analysis up to a given maximum redshift, i.e. $z=0.2,0.52,0.72,1$. The
dashed black lines show the truth values of the cosmological parameters of the
Magneticum simulations.

##  6 Conclusions

In this work we evaluated the constraining power of the combination between
two large-scale cosmological probes, namely the number counts of galaxy
clusters and cosmic voids. We studied the potential constraints achievable by
sampling the posterior distribution of the parameters $\Omega_{\rm m}$ and
$\sigma_{8}$ through the analysis of the redshift-space cluster and void
catalogues extracted from the Magneticum Box1a simulation, in four redshift
snapshots ($0.2\leq z\leq 1$). Since this simulation takes into account the
baryonic physics, we re-calibrated the HMF and VSF parameter models to
minimise systematic uncertainties. We considered only galaxy clusters having
masses $M_{500c}\geq 10^{14}\ h^{-1}\rm{M}_{\odot}$, and voids with an
effective radii $R_{\rm eff}\geq 3.6\ \lambda_{\rm mgs}(z)$, with
$\lambda_{\rm mgs}(z)$ being the mean galaxy separation from which voids are
detected. We chose these conservative selections to make our mock catalogues
similar to those expected from large-scale surveys and, at the same time, to
not be influenced by the resolution of the simulation. We computed the cross-
covariance matrix between the two probes, finding no statistically relevant
correlations between cluster mass and void radius bins at different redshifts.
By looking at the cross-covariance matrix we found, instead, correlation
features within cluster mass bins at different redshifts.

The cross-correlation terms are dominated by noise, so we have decided not to
include their contribution in the calculation of the final constraints. It
would be interesting to repeat our analysis with a larger set of mock
catalogues. This would allow us to better investigate the cross-correlation
terms, thus possibly improve the accuracy of the results. Moreover, it would
be interesting to investigate the impact of cross-correlations between
different redshift bins. This is beyond the scope of the current analysis, and
it will be investigated in a forthcoming work.

The most remarkable result of our analysis concerns the high constraining
power derived from the combination of the HMF and VSF on the $\Omega_{\rm m}$
\- $\sigma_{8}$ space: being the confidence contours computed with these
probes highly perpendicular and their cross-correlation negligible, the
cluster and void number count joint analysis is of great effectiveness. In
particular, we found that the constraining power given by the cluster counts
improves by about $60\%$ when combined with void counts. This highlights the
important contribution that void counts can provide to reduce the relative
errors on the cosmological parameters related to the growth of large-scale
structures, i.e. $\Omega_{\rm m}$, $\sigma_{8}$ and the derived quantity
$S_{8}$. Notably, we found that the improvement remains approximately constant
regardless of the maximum redshift considered.

This work is based on some modelling simplifications. Firstly, we did not take
into account any uncertainties on observable-mass scaling relation and derived
redshifts for the considered simulated tracers, as if galaxy cluster masses
were perfectly measured. Also, in the building of the redshift-space
catalogues, we mapped the tracer comoving coordinate into the observed ones by
assuming the true cosmological parameters of the simulation, taking into
account only the geometrical distortions arising from the cosmology variation
acting in the MCMC analysis. Moreover we modelled the geometric distortions
for voids only, since these effects are expected to be negligible in the
counts of galaxy clusters, the latter being less-spatially extended objects
than voids. Our study was also limited to the constraints on $\Omega_{\rm m}$
and $\sigma_{8}$ and it will be fundamental to consider also the remaining
cosmological parameters on which the analysed probes strongly depend. Finally,
we neglected the super-sample covariance (SSC) (Hu & Kravtsov, 2003; Takada &
Bridle, 2007; Takada & Spergel, 2014; Krause & Eifler, 2017; Lacasa et al.,
2018), which concerns the fact that simulations and observations map a limited
portion of the Universe. Indeed, due to the missing modes larger than the
simulation box, neglecting the SSC might introduce non-negligible effects.
Recently, Bayer et al. (2022) estimated the SSC specifically for cluster and
void statistics. Although the SSC effect can be considered negligible for
cluster and void counts alone, this is no longer true when considering their
cross-correlation, potentially causing a change in the relative improvements
of the combination analysis (Takada & Bridle, 2007; Takada & Spergel, 2014;
Lacasa & Rosenfeld, 2016). We will address all these issues in future works.

In order to obtain well grounded statistical results from the analysis of
large-scale structures, it would be interesting to search for galaxy clusters
and cosmic voids in much larger cosmological volumes. Upcoming galaxy surveys
like the Euclid mission (Laureijs et al., 2011; Scaramella et al., 2014;
Amendola et al., 2018), the Dark Energy Spectroscopic Instrument (DESI) (Hang
et al., 2021; Besuner et al., 2021) and Vera C. Rubin Observatory LSST (LSST
Dark Energy Science Collaboration, 2012) will map a very large fraction of our
Universe, with volumes of the order of $10^{2}$ ($h^{-1}$ Gpc)3, giving even
more tight constraints on cosmological parameters and, mostly, improving our
knowledge about the underlying cosmological model. The tighter constraints
achieved by combining cluster and void counts, compared to those obtained by
exploiting cluster counts alone, will significantly contribute to break cosmic
degeneracies.

In this analysis we proved how the high complementarity of these cosmological
probes, together with the great error reduction given by their combination,
will provide a fundamental contribution in shedding light on tensions possibly
affecting cosmological observations. We thus conclude that the methods
presented in this work, including model calibrations on suitable mock
catalogues, will motivate the exploitation of the combined analysis of galaxy
cluster and cosmic void counts from future wide field galaxy surveys.

## Acknowledgements

We thank the referee for the useful suggestions, which helped us to improve
the paper. We acknowledge the use of computational resources from the parallel
computing cluster of the Open Physics
Hub555https://site.unibo.it/openphysicshub/en at the Physics and Astronomy
Department in Bologna. CG and LM acknowledge the support from the grant PRIN-
MIUR 2017 WSCC32 ZOOMING, and the support from the grant ASI n.2018-23-HH.0.
CG acknowledges funding from the Italian National Institute of Astrophysics
under the grant "Bando PRIN 2019", PI: Viola Allevato and from the HPC-Europa3
Transnational Access programme HPC17VDILO. KD acknowledges support by the
COMPLEX project from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program grant agreement
ERC-2019-AdG 882679 as well as support through the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094 - 390783311. The calculations for the
hydrodynamical simulations were carried out at the Leibniz Supercomputer
Center (LRZ) under the project pr83li.

## Data Availability

The simulation data underlying this article will be shared on reasonable
request to the corresponding author.

## References

  * Abdalla et al. (2022) Abdalla E., et al., 2022, Journal of High Energy Astrophysics, 34, 49
  * Alam et al. (2017) Alam S., et al., 2017, MNRAS, 470, 2617
  * Alcock & Paczynski (1979) Alcock C., Paczynski B., 1979, Nature, 281, 358
  * Amendola et al. (2013) Amendola L., Kunz M., Motta M., Saltas I. D., Sawicki I., 2013, Phys. Rev. D, 87, 023501
  * Amendola et al. (2018) Amendola L., et al., 2018, Living Reviews in Relativity, 21, 2
  * Balaguera-Antolínez & Porciani (2013) Balaguera-Antolínez A., Porciani C., 2013, J. Cosmology Astropart. Phys., 2013, 022
  * Baldi & Villaescusa-Navarro (2016) Baldi M., Villaescusa-Navarro F., 2016, arXiv e-prints, p. arXiv:1608.08057
  * Barreira et al. (2015) Barreira A., Cautun M., Li B., Baugh C. M., Pascoli S., 2015, J. Cosmology Astropart. Phys., 2015, 028
  * Bayer et al. (2021) Bayer A. E., et al., 2021, ApJ, 919, 24
  * Bayer et al. (2022) Bayer A. E., Liu J., Terasawa R., Barreira A., Zhong Y., Feng Y., 2022, arXiv e-prints, p. arXiv:2210.15647
  * Beck et al. (2016) Beck A. M., et al., 2016, MNRAS, 455, 2110
  * Bernardeau (1994) Bernardeau F., 1994, ApJ, 427, 51
  * Besuner et al. (2021) Besuner R., et al., 2021, arXiv e-prints, p. arXiv:2101.11794
  * Bianchi et al. (2012) Bianchi D., Guzzo L., Branchini E., Majerotto E., de la Torre S., Marulli F., Moscardini L., Angulo R. E., 2012, MNRAS, 427, 2420
  * Bocquet et al. (2016) Bocquet S., Saro A., Dolag K., Mohr J. J., 2016, MNRAS, 456, 2361
  * Bocquet et al. (2019) Bocquet S., et al., 2019, ApJ, 878, 55
  * Bond et al. (1991) Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440
  * Bos et al. (2012) Bos E. G. P., van de Weygaert R., Dolag K., Pettorino V., 2012, MNRAS, 426, 440
  * Bridle et al. (1999) Bridle S. L., Eke V. R., Lahav O., Lasenby A. N., Hobson M. P., Cole S., Frenk C. S., Henry J. P., 1999, MNRAS, 310, 565
  * Carroll et al. (1992) Carroll S. M., Press W. H., Turner E. L., 1992, ARA&A, 30, 499
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chan et al. (2019) Chan K. C., Hamaus N., Biagetti M., 2019, Phys. Rev. D, 99, 121304
  * Cole (1991) Cole S., 1991, ApJ, 367, 45
  * Contarini et al. (2019) Contarini S., Ronconi T., Marulli F., Moscardini L., Veropalumbo A., Baldi M., 2019, MNRAS, 488, 3526
  * Contarini et al. (2021) Contarini S., Marulli F., Moscardini L., Veropalumbo A., Giocoli C., Baldi M., 2021, MNRAS, 504, 5021
  * Contarini et al. (2022a) Contarini S., et al., 2022a, arXiv e-prints, p. arXiv:2205.11525
  * Contarini et al. (2022b) Contarini S., Pisani A., Hamaus N., Marulli F., Moscardini L., Baldi M., 2022b, arXiv e-prints, p. arXiv:2212.03873
  * Contarini et al. (2022c) Contarini S., Pisani A., Hamaus N., Marulli F., Moscardini L., Baldi M., 2022c, arXiv e-prints, p. arXiv:2212.07438
  * Costanzi et al. (2019) Costanzi M., et al., 2019, MNRAS, 488, 4779
  * Courtin et al. (2011) Courtin J., Rasera Y., Alimi J. M., Corasaniti P. S., Boucher V., Füzfa A., 2011, MNRAS, 410, 1911
  * Cunha et al. (2009) Cunha C., Huterer D., Frieman J. A., 2009, Phys. Rev. D, 80, 063532
  * DES Collaboration (2019) DES Collaboration 2019, Physical Review Letters, 122, 171301
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * Despali et al. (2016) Despali G., Giocoli C., Angulo R. E., Tormen G., Sheth R. K., Baso G., Moscardini L., 2016, MNRAS, 456, 2486
  * Di Matteo et al. (2008) Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki D., 2008, ApJ, 676, 33
  * Di Valentino et al. (2020a) Di Valentino E., et al., 2020a, arXiv e-prints, p. arXiv:2008.11283
  * Di Valentino et al. (2020b) Di Valentino E., et al., 2020b, arXiv e-prints, p. arXiv:2008.11285
  * Di Valentino et al. (2020c) Di Valentino E., et al., 2020c, arXiv e-prints, p. arXiv:2008.11286
  * Dolag et al. (2004) Dolag K., Jubelgas M., Springel V., Borgani S., Rasia E., 2004, ApJ, 606, L97
  * Dolag et al. (2005) Dolag K., Vazza F., Brunetti G., Tormen G., 2005, MNRAS, 364, 753
  * Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
  * Dolag et al. (2016) Dolag K., Komatsu E., Sunyaev R., 2016, MNRAS, 463, 1797
  * Efron (1982) Efron B., 1982, The Jackknife, the Bootstrap and Other Resampling Plans. Society for Industrial and Applied Mathematics (https://epubs.siam.org/doi/pdf/10.1137/1.9781611970319), doi:10.1137/1.9781611970319, https://epubs.siam.org/doi/abs/10.1137/1.9781611970319
  * Escoffier et al. (2016) Escoffier S., et al., 2016, arXiv e-prints, p. arXiv:1606.00233
  * Fabjan et al. (2010) Fabjan D., Borgani S., Tornatore L., Saro A., Murante G., Dolag K., 2010, MNRAS, 401, 1670
  * Gawiser & Silk (1998) Gawiser E., Silk J., 1998, Science, 280, 1405
  * Giocoli et al. (2008) Giocoli C., Pieri L., Tormen G., 2008, MNRAS, 387, 689
  * Hamaus et al. (2014) Hamaus N., Sutter P. M., Wandelt B. D., 2014, Phys. Rev. Lett., 112, 251302
  * Hamaus et al. (2020) Hamaus N., Pisani A., Choi J.-A., Lavaux G., Wandelt B. D., Weller J., 2020, J. Cosmology Astropart. Phys., 2020, 023
  * Hamaus et al. (2022) Hamaus N., et al., 2022, A&A, 658, A20
  * Hang et al. (2021) Hang Q., Alam S., Peacock J. A., Cai Y.-C., 2021, MNRAS, 501, 1481
  * Hartlap et al. (2007) Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399
  * Hasselfield et al. (2013) Hasselfield M., et al., 2013, J. Cosmology Astropart. Phys., 2013, 008
  * Hirschmann et al. (2014) Hirschmann M., Dolag K., Saro A., Bachmann L., Borgani S., Burkert A., 2014, MNRAS, 442, 2304
  * Hong et al. (2012) Hong T., Han J. L., Wen Z. L., Sun L., Zhan H., 2012, ApJ, 749, 81
  * Hu & Kravtsov (2003) Hu W., Kravtsov A. V., 2003, ApJ, 584, 702
  * Jennings et al. (2013) Jennings E., Li Y., Hu W., 2013, MNRAS, 434, 2167
  * Joyce et al. (2016) Joyce A., Lombriser L., Schmidt F., 2016, Annual Review of Nuclear and Particle Science, 66, 95
  * Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
  * Krause & Eifler (2017) Krause E., Eifler T., 2017, MNRAS, 470, 2100
  * Kreisch et al. (2019) Kreisch C. D., Pisani A., Carbone C., Liu J., Hawken A. J., Massara E., Spergel D. N., Wandelt B. D., 2019, MNRAS, 488, 4413
  * Kreisch et al. (2022) Kreisch C. D., Pisani A., Villaescusa-Navarro F., Spergel D. N., Wandelt B. D., Hamaus N., Bayer A. E., 2022, ApJ, 935, 100
  * LSST Dark Energy Science Collaboration (2012) LSST Dark Energy Science Collaboration 2012, arXiv e-prints, p. arXiv:1211.0310
  * Lacasa & Rosenfeld (2016) Lacasa F., Rosenfeld R., 2016, J. Cosmology Astropart. Phys., 2016, 005
  * Lacasa et al. (2018) Lacasa F., Lima M., Aguena M., 2018, A&A, 611, A83
  * Landy & Szalay (1993) Landy S. D., Szalay A. S., 1993, ApJ, 412, 64
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Lesci et al. (2022a) Lesci G. F., et al., 2022a, A&A, 659, A88
  * Lesci et al. (2022b) Lesci G. F., et al., 2022b, A&A, 665, A100
  * Mana et al. (2013) Mana A., Giannantonio T., Weller J., Hoyle B., Hütsi G., Sartoris B., 2013, MNRAS, 434, 684
  * Manera et al. (2010) Manera M., Sheth R. K., Scoccimarro R., 2010, MNRAS, 402, 589
  * Marulli et al. (2011) Marulli F., Carbone C., Viel M., Moscardini L., Cimatti A., 2011, MNRAS, 418, 346
  * Marulli et al. (2012a) Marulli F., Baldi M., Moscardini L., 2012a, MNRAS, 420, 2377
  * Marulli et al. (2012b) Marulli F., Bianchi D., Branchini E., Guzzo L., Moscardini L., Angulo R. E., 2012b, MNRAS, 426, 2566
  * Marulli et al. (2013) Marulli F., et al., 2013, A&A, 557, A17
  * Marulli et al. (2016) Marulli F., Veropalumbo A., Moresco M., 2016, Astronomy and Computing, 14, 35
  * Marulli et al. (2017) Marulli F., Veropalumbo A., Moscardini L., Cimatti A., Dolag K., 2017, A&A, 599, A106
  * Marulli et al. (2018) Marulli F., et al., 2018, A&A, 620, A1
  * Marulli et al. (2021) Marulli F., Veropalumbo A., García-Farieta J. E., Moresco M., Moscardini L., Cimatti A., 2021, ApJ, 920, 13
  * Massara et al. (2015) Massara E., Villaescusa-Navarro F., Viel M., Sutter P. M., 2015, J. Cosmology Astropart. Phys., 2015, 018
  * Merloni et al. (2012) Merloni A., et al., 2012, arXiv e-prints, p. arXiv:1209.3114
  * Mo & White (1996) Mo H. J., White S. D. M., 1996, MNRAS, 282, 347
  * Moresco et al. (2022) Moresco M., et al., 2022, arXiv e-prints, p. arXiv:2201.07241
  * Nadathur et al. (2015) Nadathur S., Hotchkiss S., Diego J. M., Iliev I. T., Gottlöber S., Watson W. A., Yepes G., 2015, MNRAS, 449, 3997
  * Nadathur et al. (2020) Nadathur S., et al., 2020, MNRAS, 499, 4140
  * Neyrinck (2008) Neyrinck M. C., 2008, MNRAS, 386, 2101
  * Norberg et al. (2009) Norberg P., Baugh C. M., Gaztañaga E., Croton D. J., 2009, MNRAS, 396, 19
  * Padovani & Matteucci (1993) Padovani P., Matteucci F., 1993, ApJ, 416, 26
  * Parr & Schucany (1980) Parr W. C., Schucany W. R., 1980, International Statistical Review / Revue Internationale de Statistique, 48, 73
  * Percival et al. (2014) Percival W. J., et al., 2014, MNRAS, 439, 2531
  * Perivolaropoulos & Skara (2021) Perivolaropoulos L., Skara F., 2021, arXiv e-prints, p. arXiv:2105.05208
  * Perlmutter et al. (1999) Perlmutter S., et al., 1999, ApJ, 517, 565
  * Pisani et al. (2015a) Pisani A., Sutter P. M., Wandelt B. D., 2015a, arXiv e-prints, p. arXiv:1506.07982
  * Pisani et al. (2015b) Pisani A., Sutter P. M., Hamaus N., Alizadeh E., Biswas R., Wandelt B. D., Hirata C. M., 2015b, Phys. Rev. D, 92, 083531
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A20
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A24
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Platen et al. (2008) Platen E., van de Weygaert R., Jones B. J. T., 2008, MNRAS, 387, 128
  * Pollina et al. (2017) Pollina G., Hamaus N., Dolag K., Weller J., Baldi M., Moscardini L., 2017, MNRAS, 469, 787
  * Pollina et al. (2019) Pollina G., et al., 2019, MNRAS, 487, 2836
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
  * Ragagnin et al. (2021) Ragagnin A., Saro A., Singh P., Dolag K., 2021, MNRAS, 500, 5056
  * Reed et al. (2015) Reed D. S., Schneider A., Smith R. E., Potter D., Stadel J., Moore B., 2015, MNRAS, 451, 4413
  * Ricciardelli et al. (2014) Ricciardelli E., Quilis V., Varela J., 2014, MNRAS, 440, 601
  * Riess et al. (1998) Riess A. G., et al., 1998, AJ, 116, 1009
  * Ronconi & Marulli (2017) Ronconi T., Marulli F., 2017, A&A, 607, A24
  * Ronconi et al. (2019) Ronconi T., Contarini S., Marulli F., Baldi M., Moscardini L., 2019, MNRAS, 488, 5075
  * Sahlén et al. (2016) Sahlén M., Zubeldía Í., Silk J., 2016, ApJ, 820, L7
  * Salvati et al. (2018) Salvati L., Douspis M., Aghanim N., 2018, A&A, 614, A13
  * Scaramella et al. (2014) Scaramella R., et al., 2014, Proceedings of the International Astronomical Union, 10, 375
  * Schuster et al. (2019) Schuster N., Hamaus N., Pisani A., Carbone C., Kreisch C. D., Pollina G., Weller J., 2019, J. Cosmology Astropart. Phys., 2019, 055
  * Sheth & Tormen (1999) Sheth R. K., Tormen G., 1999, MNRAS, 308, 119
  * Sheth & van de Weygaert (2004) Sheth R. K., van de Weygaert R., 2004, MNRAS, 350, 517
  * Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1
  * Spolyar et al. (2013) Spolyar D., Sahlén M., Silk J., 2013, Phys. Rev. Lett., 111, 241103
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel & Hernquist (2002) Springel V., Hernquist L., 2002, MNRAS, 333, 649
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Sutter et al. (2015) Sutter P. M., et al., 2015, Astronomy and Computing, 9, 1
  * Takada & Bridle (2007) Takada M., Bridle S., 2007, New Journal of Physics, 9, 446
  * Takada & Spergel (2014) Takada M., Spergel D. N., 2014, MNRAS, 441, 2456
  * Teklu et al. (2015) Teklu A. F., Remus R.-S., Dolag K., Beck A. M., Burkert A., Schmidt A. S., Schulze F., Steinborn L. K., 2015, ApJ, 812, 29
  * Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
  * Tormen (1998) Tormen G., 1998, MNRAS, 297, 648
  * Tormen et al. (2004) Tormen G., Moscardini L., Yoshida N., 2004, MNRAS, 350, 1397
  * Tornatore et al. (2003) Tornatore L., Borgani S., Springel V., Matteucci F., Menci N., Murante G., 2003, MNRAS, 342, 1025
  * Tornatore et al. (2007) Tornatore L., Borgani S., Dolag K., Matteucci F., 2007, MNRAS, 382, 1050
  * Veropalumbo et al. (2016) Veropalumbo A., Marulli F., Moscardini L., Moresco M., Cimatti A., 2016, MNRAS, 458, 1909
  * Verza et al. (2019) Verza G., Pisani A., Carbone C., Hamaus N., Guzzo L., 2019, J. Cosmology Astropart. Phys., 2019, 040
  * Vikhlinin et al. (2009) Vikhlinin A., et al., 2009, ApJ, 692, 1060
  * Villaescusa-Navarro et al. (2013) Villaescusa-Navarro F., Vogelsberger M., Viel M., Loeb A., 2013, MNRAS, 431, 3670
  * Webster et al. (1998) Webster M., Bridle S. L., Hobson M. P., Lasenby A. N., Lahav O., Rocha G., 1998, arXiv e-prints, pp astro–ph/9802109
  * Wu et al. (2022) Wu P.-J., Shao Y., Jin S.-J., Zhang X., 2022, arXiv e-prints, p. arXiv:2202.09726
  * Yang et al. (2015) Yang L. F., Neyrinck M. C., Aragón-Calvo M. A., Falck B., Silk J., 2015, MNRAS, 451, 3606
  * Yoo & Watanabe (2012) Yoo J., Watanabe Y., 2012, International Journal of Modern Physics D, 21, 1230002
  * Zentner (2007) Zentner A. R., 2007, International Journal of Modern Physics D, 16, 763

##  Appendix A Covariance matrix convergence

In this work, we derived cosmological constraints from hydrodynamical
simulations by estimating the cross-covariance matrix, $\mathcal{C}$, between
galaxy cluster and void number counts. It is important to verify that the
noise in the covariance does not have a significant impact on the results, as
it could potentially affect the estimation of the cosmological parameters in
the Bayesian analysis. Therefore, in this Appendix we test the stability of
our results using different sampling techniques.

In Section 4.3 we computed $\mathcal{C}$ by dividing the tracer catalogues
(galaxy clusters and voids) into $N_{\rm sub}=125$ equally sized sub-regions.
To test the covariance convergence and stability we compare now both jackknife
and bootstrap resampling methods (Efron, 1982; Norberg et al., 2009). When
adopting the former, we re-compute $\mathcal{C}$ as in Section 4.3, but
varying this time $N_{\rm sub}\in[8,1000]$. Then, we extract again the
combined constraints on $\Omega_{\rm m}$ and $\sigma_{8}$ for each considered
case. When adopting the bootstrap technique, we derive instead the
cosmological constraints as a function of the number of resamplings $N_{\rm
mock}$, considering a range of $N_{\rm mock}\in[100,2400]$. For these tests we
fixed the nuisance model parameters of HMF and VSF (i.e.
$a_{1},a_{2},a_{3},A_{1},A_{2},B_{\mathrm{slope}},B_{\mathrm{off}}$) to those
calibrated in Section 5.1, considering all data sets up to $z=1$.

Figure 8: Posterior median values (black solid lines) for $\Omega_{\rm m}$ and
$\sigma_{8}$ (top and bottom panels, respectively) as a function of the number
of jackknife sub-boxes (left panels) or the number of bootstrap resamplings
(right panels) used for the cross-covariance matrix estimation. The presented
constraints are obtained considering all data sets up to $z=1$. The shaded
gray regions represent the $68\%$ confidence intervals around the posterior
median of the presented cases. We report with blue shaded bands the $68\%$
confidence intervals of the constraints obtained by applying the jackknife
sampling technique with $N_{\mathrm{sub}}=125$ (see Section 5.2). The dotted-
dashed black lines represent the true cosmological values of the Magneticum
simulations.

In Figure 8 we report the median values for $\Omega_{\rm m}$ and $\sigma_{8}$
and the corresponding $68\%$ uncertainties obtained for the different cases
analysed. We note that the cosmological constraints are stable and consistent
with those reported in Eq. (33), for all the methodologies employed. We
conclude that the final constraints presented in this work, obtained by
adopting the jackknife sampling technique with $N_{\mathrm{sub}}=125$, can be
considered unaffected by relevant numerical issues related to the estimation
of the covariance matrix.