© The Authors 2023

1. Introduction

Cosmic voids are underdense regions in the large-scale structure (LSS) which together with halos, walls, and filaments make up the cosmic web. They have recently become a popular independent probe for cosmological parameters as they provide information about the parameter combination fσ₈(z) through the study of redshift space distortions (RSDs) from the imprint left on the quadrupole of the void-galaxy cross-correlation function (CCF; Cai et al. 2016; Nadathur & Percival 2019; Nadathur et al. 2019b). Further on, careful modelling of the RSDs around voids originating from the peculiar velocities of galaxies also opens up for the applied fiducial cosmology to be tested through the Alcock-Paczyǹski effect (Alcock & Paczyński 1979; Sutter et al. 2012; Hamaus et al. 2022). In addition to this, the empty nature of cosmic voids makes them sensitive to diffuse components of our Universe, such as massive neutrinos, and effects of modified gravities which scale inversely with density (Zivick et al. 2015; Cai et al. 2015; Massara et al. 2015; Falck et al. 2018; Kreisch et al. 2019; Fiorini et al. 2022).

The presence of massive neutrinos suppresses structure formation on scales smaller than the neutrino free-streaming length (Lesgourgues & Pastor 2006), while f(R)-modified gravity enhances it on scales smaller than the Compton wavelength of the scalaron (e.g., Cataneo et al. 2015). This results in observational degeneracies which have previously been studied in simulations (Baldi et al. 2014; Giocoli et al. 2018; Contarini et al. 2021). Both massive neutrinos and modified gravity are independently large research fields. Determining the absolute mass scale of neutrinos is not only an important quest within particle physics but also essential in order to understand LSS formation within cosmology. The massive neutrinos make up a fraction of the matter content of the Universe ranging from 0.5 − 2% based on the lower and upper limits of the sum of the neutrino masses, 0.06 eV ≲ ∑m_ν ≲ 2.4 eV, provided by particle physics experiments (Particle Data Group 2022; Aker et al. 2022). This range can be further constrained through cosmological observations. Effects of massive neutrinos on LSS and signatures left in the cosmic microwave background (CMB; Archidiacono et al. 2017) can be used to estimate the sum of the neutrino masses. New space missions, such as Euclid (Laureijs et al. 2011), are now sensitive enough to accurately pick up on the effects of massive neutrinos on the matter power spectrum. Euclid aims to measure the sum of the neutrino masses to a precision better than 0.03 eV. This will be achieved through a combined analysis of galaxy clustering and weak gravitational lensing. For small neutrino masses, ∑m_ν < 0.1 eV, this is precise enough to determine the neutrino mass hierarchy.

Although constraints from cosmology are typically tighter than the ones obtained from particle physics experiments, they do depend on an assumed cosmological model. Because of this, disentangling degeneracies between massive neutrinos and various cosmologies is important in order to fully determine the neutrino mass scale. Accounting for massive neutrinos when studying alternatives to ΛCDM (Lambda cold dark matter) is therefore also necessary to make sure that constraints on the model in question are properly estimated. For f(R)-modified gravity, the constraints on the model parameters might be slightly alleviated when massive neutrinos are included in the estimates, due to their degeneracy.

The theory of general relativity (GR) has been thoroughly investigated in high density regions, for instance through local solar system tests (Will 2014). To adhere to the resulting measurements, modified gravity theories typically need to be screened as a function of density. This contributes to the constraints on the parameters of these models and also makes deviations from ΛCDM hard to observe in highly screened environments. Cosmic voids are underdense regions sensitive to both massive neutrinos and f(R)-modified gravity. Because of this, voids have been proposed as a means to separate their known degenerate effects. The void size function at high redshifts for large voids is found by Contarini et al. (2021) to be a promising candidate for the task.

Before further exploring voids in the context of massive neutrinos and f(R)-modified gravity, we must first take a look at the theory behind the models developed to extract cosmological information from them. In our case, we investigate the void-galaxy CCF in redshift space, the void velocity profile, and a reconstruction method for putting galaxy positions observed from redshift space back into real space. These are all based on linear theory and, to a degree, a standard ΛCDM universe. In this paper, we explore the need for change in these models in order to use them to describe voids in a universe with massive neutrinos and f(R)-modified gravity. To do so, we compared them to N-body simulations containing massive neutrinos and f(R)-modified gravity, both independently and combined. We studied how well the model for the void velocity profile and the void-halo CCF in redshift space fit the simulation data in each case. The latter was done both by using all the real space information from the simulations directly as model input and also by treating the simulation data as observational data, using redshift space halo positions and a reconstruction method to gain an estimate of the corresponding real space positions of the halos. In addition, we performed Markov chain Monte Carlo (MCMC) fits for fσ₈(z) and the Alcock-Paczyǹski parameter, ϵ, in order to test how well we are able to recover the fiducial cosmology.

This paper is structured as follows: We start by presenting background theory for cosmic voids, massive neutrinos, and f(R)-modified gravity in Sect. 2 and then explain our simulation set-up in Sect. 3. In Sect. 4 we detail the methodical approach of the paper before reporting our results in Sect. 5. Finally, we conclude in Sect. 6.

2. Theory

2.1. Cosmic voids

In the following subsections, we recap some models developed to extract information from cosmic voids. This includes the modelling of the void-galaxy CCF in redshift space, the velocity profile, and a reconstruction process. In addition, we briefly explain how voids can be used to test cosmology through the Alcock-Paczyński effect. The numerical definition of voids in this work is addressed in Sect. 4.2.

2.1.1. Cross-correlation function

The void-galaxy CCF gives us information about how galaxies are distributed around voids as a function of the void-centre galaxy separation (Cai et al. 2016; Nadathur & Percival 2019; Woodfinden et al. 2022). In real space, we denote this separation vector by r, and in redshift space by s. As we observe in redshift space, deconstructing the separation vectors into components perpendicular and parallel to the line-of-sight (LOS) is advantageous, and gives the following relations

$$\begin{array}{cc}& s_{\perp }=r_{\perp },\hfill \end{array}$$(1)

$$\begin{array}{cc}& s_{\Vert }=r_{\Vert }+\frac{v_{\Vert }}{\mathit{aH}}.\hfill \end{array}$$(2)

Here, H is the Hubble rate, a the scale factor, and v_∥ the component of the galaxy peculiar velocity parallel to the LOS. The LOS can be defined as from the observer to the void centre, the galaxy or the midpoint between them. We use the observer to void-centre definition in this paper.

The void-galaxy CCF in redshift space, ξ^s(s), can be related to the void-galaxy CCF in real space, ξ^r(r), by the streaming model (Peebles 1980; Fisher 1995; Paillas et al. 2021):

$$\begin{array}{c}\hfill 1+{\xi }^s(\mathbf{s})={\displaystyle \int }(1+{\xi }^r(\mathbf{r}))P(v_{\Vert },\mathbf{r})\mathrm{d}{v}_{\Vert },\end{array}$$(3)

where P(v_∥, r) is the probability distribution function (PDF) of the galaxies’ peculiar velocities parallel to the LOS. This is a mapping from real to redshift space, based on the fact that void-galaxy pairs are conserved when moving between the two¹. The galaxy velocities can further be separated into two components, one describing the coherent outflow velocity away from the void centre, and the other an additional stochastic motion of the galaxies. The first term can be assumed spherically symmetric around and radially directed out from the void centre when averaging over a large number of voids, due to statistical isotropy in real space. Following this argument, we can write the LOS galaxy velocity component as

$$\begin{array}{c}\hfill v_{\Vert }(r,{\mu }_r)=v_r(r){\mu }_r+{\stackrel{\sim }{v}}_{\Vert },\end{array}$$(4)

where ${\stackrel{\sim }{v}}_{\Vert }$ is the random velocity component parallel to the LOS, μ_r = cos θ, and θ is the angle between the LOS and r so that v_r(r)μ_r is the part of the coherent outflow velocity profile projected along the LOS.

By performing a coordinate change from $v_{\Vert }\to {\stackrel{\sim }{v}}_{\Vert }$, as shown in Woodfinden et al. (2022), Eq. (3) can be written as

$$\begin{array}{c}\hfill 1+{\xi }^s(\mathbf{s})={\displaystyle \int }(1+{\xi }^r(\mathbf{r}))P({\stackrel{\sim }{v}}_{\Vert },\mathbf{r})J_{\mathbf{rs}}\mathrm{d}{\stackrel{\sim }{v}}_{\Vert },\end{array}$$(5)

where J_rs is the Jacobian resulting from the shift from real to redshift space and $P({\stackrel{\sim }{v}}_{\Vert },\mathbf{r})$ is the PDF of the random velocity component centred around zero. Writing out the Jacobian we get the expression

$$\begin{array}{c}\hfill J_{\mathbf{rs}}=[1+\frac{v_r}{\mathit{raH}}+\frac{(v_r^{\prime }-v_r/r)}{\mathit{aH}}{\mu }_r^2{]}^{-1},\end{array}$$(6)

with derivatives with respect to r denoted by a prime. Simulations show that $P({\stackrel{\sim }{v}}_{\Vert },\mathbf{r})$ is approximated well by a Gaussian (Nadathur & Percival 2019; Paillas et al. 2021) and we therefore consider

$$\begin{array}{c}\hfill P({\stackrel{\sim }{v}}_{\Vert },\mathbf{r})=\frac{1}{\sqrt{2\pi {\sigma }_{v_{\Vert }}^2(r)}}exp(-\frac{{\stackrel{\sim }{v}}_{\Vert }^2}{2{\sigma }_{v_{\Vert }}^2(r)}),\end{array}$$(7)

with the additional assumption of a spherically symmetric velocity dispersion, ${\sigma }_{v_{\parallel }}^2(r)$.

2.1.2. Velocity modelling

Equation (5) depends on the coherent mean outflow velocity profile around the void, $\mathbf{v}(r)=v_r(r)\widehat{\mathbf{r}}$. This quantity can be provided to the model through a template, for example, obtained from a simulation, or it can be modelled separately. Typically, using linear perturbation theory and the continuity equation (Peebles 1980, 1993), the expression

$$\begin{array}{c}\hfill v_r^{\mathrm{Lin}}(r)=-\frac{1}{3}faH\mathrm{\Delta }(r)r\end{array}$$(8)

is used for the radial outflow velocity profile. Here, f = dlnD/dlna is the linear growth rate, with D as the growth factor, and Δ(r) is the integrated density,

$$\begin{array}{c}\hfill \mathrm{\Delta }(r)=\frac{3}{r^3}{\displaystyle {\int }_0^r}\delta (x)x^2\mathrm{d}x,\end{array}$$(9)

with δ(r) the dark matter overdensity profile around the void. Nadathur et al. (2019b) show that Δ(r)∝σ₈, which leads to v_r(r) depending on the combination of the parameters fσ₈, where σ₈ is the amplitude of the linear matter power spectrum at a scale of 8 h⁻¹Mpc. It is also worth noting that in the process of deriving Eq. (8), the growth rate is assumed constant, which does not hold for cosmologies with massive neutrinos or modified gravity, where the growth rate is scale-dependent (e.g., Hernández 2017; Mirzatuny & Pierpaoli 2019).

Although Eq. (8) is commonly used to model the coherent outflow velocity, the expression in Eq. (5) should hold for any spherically symmetric velocity profile. One other example of this is the more general velocity profile from Paillas et al. (2021), which introduces an extra degree of freedom,

$$\begin{array}{c}\hfill v_r^{\mathrm{Gen}}(r)=-\frac{1}{3}\frac{faH\mathrm{\Delta }(r)r}{1+A_v{\xi }_0(r)}.\end{array}$$(10)

Here, ξ₀(r) is the void-galaxy CCF monopole and A_v is a free parameter. This expression is able to better fit the void velocity profile by adjusting A_v. However, the value of A_v must be estimated from, for example, simulations if we want to use Eq. (10) for void-galaxy CCF modelling.

2.1.3. Reconstruction

When working with observational data, we only have access to galaxy positions in redshift space. However, the model of the void-galaxy CCF in redshift space depends on the real space void-galaxy CCF. One way to obtain this quantity is by using a reconstruction method to estimate the positions of galaxies in real space.

In the context of voids, Nadathur et al. (2019a) propose one such method, based on Nusser & Davis (1994). It involves solving the Zel’dovich equation (Zel’dovich 1970),

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathrm{\Psi }+\frac{f}{b}\mathrm{\nabla }·(\mathrm{\Psi }·\widehat{\mathbf{r}})\widehat{\mathbf{r}}=-\frac{{\delta }_g}{b},\end{array}$$(11)

for the displacement field, Ψ. Here, f is the growth rate and b is the linear bias relating the galaxy overdensity field, δ_g, to the matter overdensity field, δ, in redshift space. The displacement field, Ψ, relates the Eulerian and Lagrangian positions of the galaxies. A part of this displacement, given by

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathrm{RSD}}=-f(\mathrm{\Psi }·\widehat{\mathbf{r}})\widehat{\mathbf{r}},\end{array}$$(12)

is due to the linear RSDs resulting from the galaxies’ peculiar velocities. By shifting the galaxy positions in redshift space by −Ψ_RSD, we can approximately obtain their real space positions. As seen above, this procedure depends on an input value of the growth rate, f, and the linear galaxy bias, b, which combined make the reconstruction parameter β = f/b. Consequentially, the reconstruction method depends on a fiducial cosmological model through the growth rate.

2.1.4. The Alcock-Paczyński effect

One way to inspect a chosen cosmology is through the Alcock-Paczyński (AP) effect (Alcock & Paczyński 1979). This test only involves the geometry of the Universe and is performed by studying the ratio of the observed angular and redshift size of an object with a known shape. An example is cosmic voids, which on average should have a spherical configuration. A correctly chosen fiducial cosmology should, when converting from observed redshifts to physical distances, reproduce the spherical shape of the average void. Deviations from this show up as anisotropies in the void-galaxy CCF in addition to the contribution from peculiar velocities. If we can model the RSDs from velocities accurately enough, the AP effect can be used on voids as a test of the fiducial cosmology.

We parameterise the AP effect through the parameters

$$\begin{array}{c}\hfill {\alpha }_{\Vert }=\frac{D_H(z)}{D_H^{\mathrm{fid}}(z)},{\alpha }_{\perp }=\frac{D_{\mathrm{A}}(z)}{D_{\mathrm{A}}^{\mathrm{fid}}(z)},\end{array}$$(13)

which are distance ratios parallel and perpendicular to the LOS between the true cosmology and the fiducial cosmological model². Here, D_H(z) = c/H(z), which is the Hubble distance at redshift z, and D_A(z) is the comoving angular diameter distance.

Taking into account the difference between the true cosmology and the fiducial model, the void-galaxy CCF in redshift space scales like

$$\begin{array}{c}\hfill {\xi }^s(s_{\Vert },s_{\perp })={\xi }^{s,\mathrm{fid}}({\alpha }_{\Vert }{s}_{\Vert }^{\mathrm{fid}},{\alpha }_{\perp }{s}_{\perp }^{\mathrm{fid}})\end{array}$$(14)

due to the void-galaxy pair separation, s^fid, dependency of the calculation. The two AP parameters can be rewritten into

$$\begin{array}{cc}\hfill \alpha & ={\alpha }_{\perp }^{2/3}{\alpha }_{\Vert }^{1/3},\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill ϵ& =\frac{{\alpha }_{\perp }}{{\alpha }_{\Vert }},\hfill \end{array}$$(16)

where α characterises the volume dilation and ϵ quantifies anisotropic distortions, which largely affect the quadrupole of the CCF, as Nadathur et al. (2019b) demonstrate. When performing the calculation of the CCF model, the input parameters ξ^r(r), Δ(r), and σ_v∥ are scaled by α in order to not use the absolute void size as a standard ruler (Nadathur et al. 2019a). This leaves only a dependence on ϵ. The underlying assumption behind the AP modelling in Eq. (14) is that the positions of voids in a fiducial cosmology where the AP parameters are not unity, is (in a statistical sense) just an AP stretching of the void positions in the true cosmology. This is not always the case, and the implications of this assumption for extracting constraints from the void-galaxy CCF in future surveys is the topic of an upcoming paper.

2.2. Massive neutrinos

From neutrino flavour oscillation experiments, we know that at least two of the three neutrino mass states have a non-zero mass (Fukuda et al. 1998; Ahmad et al. 2002; Araki et al. 2005). Currently, the flavour oscillation experiments of solar and atmospheric neutrinos put constraints on the differences between the neutrino mass states given by (Particle Data Group 2022)

$$\begin{array}{cc}& \mathrm{\Delta }{m}_{21}^2=(7.53\pm 0.18)\times {10}^{-5}{\mathrm{eV}}^2,\hfill \end{array}$$(17)

$$\begin{array}{cc}& \mathrm{\Delta }{m}_{32}^2=(-2.536\pm 0.034)\times {10}^{-3}{\mathrm{eV}}^2(\mathrm{IH}),\hfill \end{array}$$(18)

$$\begin{array}{cc}& \mathrm{\Delta }{m}_{32}^2=(2.453\pm 0.033)\times {10}^{-3}{\mathrm{eV}}^2(\mathrm{NH}).\hfill \end{array}$$(19)

These experiments are not yet sensitive enough to distinguish between the ordering of the mass states, and we, therefore, have two possibilities: the normal hierarchy (NH), where m₁ < m₂ ≪ m₃, and the inverted hierarchy (IH), where m₃ ≪ m₁ < m₂. Put together, these measurements provide a lower bound on the sum of the neutrino masses at ∑m_ν ≳ 0.06 eV and ∑m_ν ≳ 0.1 eV for the NH and IH respectively. An upper bound on the sum of the neutrino masses is also provided by particle physics and the most recent constraint, at ∑m_ν ≲ 2.4 eV, comes from the KATRIN experiment, which investigates the single β-decay of molecular tritium (Aker et al. 2022).

Massive neutrinos make up a fraction of the matter content of the Universe, Ω_ν, described through the relation (Lesgourgues & Pastor 2006)

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\nu }{h}^2\approx \frac{\sum m_{\nu }}{93.14\mathrm{eV}}.\end{array}$$(20)

Typically, this is taken out of the dark matter budget, Ω_c, so that the total matter content, Ω_m = Ω_c + Ω_ν + Ω_b, stays constant. Here, Ω_b and Ω_m denote the energy density of baryons and the total matter respectively. As the massive neutrinos are relativistic at early times, they free-stream out of local peaks in the density throughout the first stages of structure formation (Lesgourgues & Pastor 2006). This, in addition to massive neutrinos altering the background evolution, leads to a suppression of matter fluctuations on scales smaller than the neutrino free-streaming length. The strength of this scale-dependent suppression depends on the neutrino mass and shows up in the linear matter power spectrum at small scales as

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{P}_m}{P_m}\approx -8f_{\nu }.\end{array}$$(21)

This is a comparison between a cosmology with massless neutrinos and massive neutrinos, where P_m is the linear matter power spectrum considering the total matter and f_ν = Ω_ν/Ω_m. At non-linear scales, this suppression is known to be even stronger, with a maximum around k ∼ 1 h Mpc⁻¹, followed by a turn-around which results in a spoon-like shape of the matter power spectrum ratio (Hannestad et al. 2020).

The presence of massive neutrinos results in observable features in the LSS and CMB (Will 2014), which makes it possible to obtain an upper bound on the sum of the neutrino masses through cosmological observations. Recently, Di Valentino et al. (2021) find ∑m_ν ≲ 0.09 eV at a 95% confidence level by analysing a combination of data sets. This puts pressure on the IH, which might further be confirmed by the Euclid space mission (Laureijs et al. 2011). It should, however, be noted that although cosmology can provide a tighter upper bound on the sum of the neutrinos masses than current particle physics experiments, the analysis depends on the choice of a cosmological model.

2.3. f(R)-modified gravity

In the concordance ΛCDM model, we have the general theory of relativity (GR) with the Einstein-Hilbert action

$$\begin{array}{c}\hfill S={\displaystyle \int }\frac{R}{16\pi G}\sqrt{-g}{\mathrm{d}}^{4}x,\end{array}$$(22)

where R is the Ricci scalar and g is the determinant of the metric tensor, g_μν. To derive the Einstein field equations, we also need to insert the Lagrangian density of matter, ℒ_m, which describes the matter in the theory,

$$\begin{array}{c}\hfill S={\displaystyle \int }(\frac{R}{16\pi G}+{\mathcal{L}}_m)\sqrt{-g}{\mathrm{d}}^{4}x.\end{array}$$(23)

In f(R)-modified gravity theories, the Einstein-Hilbert action is modified by adding a function depending on the Ricci scalar,

$$\begin{array}{c}\hfill S={\displaystyle \int }(\frac{R+f(R)}{16\pi G}+{\mathcal{L}}_m)\sqrt{-g}{\mathrm{d}}^{4}x.\end{array}$$(24)

For the specific case of Hu-Sawicki f(R)-modified gravity, we have

$$\begin{array}{c}\hfill f(R)=-m^2\frac{c_1{(R/m^2)}^n}{c_2{(R/m^2)}^n+1},\end{array}$$(25)

where $m^2=H_0^2{\mathrm{\Omega }}_m$ and c₁, c₂, and n are constant, dimensionless, non-negative parameters that describe the model (Hu & Sawicki 2007). Requiring that the model gives dark energy in the form of a cosmological constant, these three parameters can be reduced to two, n and f_R0, where

$$\begin{array}{c}\hfill f_{R0}=-n\frac{c_1}{c_2^2}{\left(\frac{{\mathrm{\Omega }}_m}{3({\mathrm{\Omega }}_m+4{\mathrm{\Omega }}_{\mathrm{\Lambda }})}\right)}^{n+1}\end{array}$$(26)

and c₁/c₂ = 6Ω_Λ/Ω_m. Increasing the value of |f_R0| gives greater deviations from GR, which has an effect on structure formation. This was demonstrated in Hu & Sawicki (2007), where raising the value of |f_R0| showed an enhancement in the matter power spectrum on small scales.

The f(R) model has become the ‘fiducial’ model used to investigate cosmological signatures of modified gravity. It is well studied in the literature and constraints have already been produced based on various cosmological probes. For example, an analysis by Cataneo et al. (2015) including cluster, CMB, supernova, and baryon acoustic oscillation (BAO) data obtains an upper bound on the Hu-Sawicki f(R) theory given by log₁₀|f_R0|<  − 4.79 for n = 1 at 95.4% confidence level. This familiarity is the reason why we chose Hu-Sawicki f(R) as our modified gravity model for this work. We set the fiducial value to |f_R0| = 10⁻⁵, as this is the value that the LSS is currently able to probe (Koyama 2016). In addition, we occasionally include some analysis with |f_R0| = 10⁻⁴ to enhance the effects of the modified theory, even though this value has been ruled out by some probes (e.g., Cataneo et al. 2015).

2.4. Massive neutrinos and f(R) gravity degeneracy

Modifications of gravity are often described as adding a fifth force which, in addition to gravity, works attractively. This fifth force is in f(R)-modified gravity carried by the scalaron, the scalar degree of freedom, quantified by df/dR. The Compton wavelength of the scalaron, λ_C, determines the range of the fifth force, which also establishes the scale on which Hu-Sawicki f(R)-modified gravity enhances structure formation (e.g., Llinares et al. 2014; Cataneo et al. 2015),

$$\begin{array}{c}\hfill {\lambda }_C^0=3\sqrt{\frac{(n+1)}{{\mathrm{\Omega }}_m+4{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}\sqrt{\frac{|f_{R0}|}{{10}^{-6}}}{h}^{-1}\mathrm{Mpc}.\end{array}$$(27)

For |f_R0| = 10⁻⁵ and n = 1, this corresponds to scales k ≳ 0.1 h Mpc⁻¹ at redshift zero.

Massive neutrinos change the time of matter-radiation equality, leading to changes in the LSS. In addition, they suppress structure formation on scales smaller than the neutrino free-streaming length (Lesgourgues & Pastor 2006),

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{FS}}=7.7\frac{1+z}{\sqrt{{\mathrm{\Omega }}_{\mathrm{\Lambda }}+{\mathrm{\Omega }}_m{(1+z)}^3}}(\frac{1\mathrm{eV}}{\sum m_{\nu }})h^{-1}\mathrm{Mpc},\end{array}$$(28)

due to their inability to cluster as relativistic particles. This corresponds to scales of k ≳ 0.02 h Mpc⁻¹ for ∑m_ν = 0.15 eV at redshift zero.

The scales where structure formation is enhanced by f(R) gravity coincide with the neutrino free-streaming length. The suppression effect of the massive neutrinos can therefore, depending on neutrino mass and the value of f_R0, counteract the structure enhancement of f(R)-modified gravity. This degeneracy has been shown to affect observables such as the matter power spectrum, the halo mass function (HMF) and the halo bias (Baldi et al. 2014).

3. Simulations

A variety of simulations incorporating massive neutrinos and modified gravity in different ways already exist (e.g., Adamek et al. 2022; Winther et al. 2015), although fewer take both into account simultaneously (e.g., Baldi et al. 2014; Giocoli et al. 2018). We approach this in our own way, by developing a RAMSES-based code, ANUBISIS, which offers the option of both massive neutrinos and several different modified gravity theories. This code is a merger between ANUBIS³ and ISIS, as explained in further detail below.

3.1. ANUBIS

For massive neutrino simulations, we have implemented neutrino particles in the N-body and hydrodynamical code RAMSES (Teyssier 2002). RAMSES implements the particle mesh (PM) method with adaptive mesh refinement (AMR), which enables higher resolution in denser regions of the simulation box. As RAMSES is a Newtonian N-body code, the most straightforward way to add neutrinos is to implement them as their own particle family and alter the equations of motion (EOM) to handle arbitrarily high momenta. Originally, the gravitational potential is calculated by solving the Poisson equation, but we instead solved the geodesic equation written in terms of canonical momentum, as commonly done for N-body codes including relativistic particles (Ma & Bertschinger 1994; Adamek et al. 2017, 2022). As a first approach, we ignored frame-dragging, scattering of particles on gravitational waves, and anisotropic stress, leaving us with the EOMs

$$\begin{array}{cc}& \frac{\partial q_i}{\partial \tau }=-\frac{2q^2+m^2{a}^2}{\sqrt{q^2+m^2{a}^2}}\frac{\partial {\varphi }_i}{\partial x^i},\hfill \end{array}$$(29)

$$\begin{array}{cc}& \frac{\partial x_i}{\partial \tau }=\frac{q_i}{\sqrt{q^2+m^2{a}^2}},\hfill \end{array}$$(30)

where q is the canonical momentum, a is the scale factor, m is the cold dark matter (CDM) particle mass, ∂ϕ/∂xⁱ is the gravitational force (which for modified gravity models includes additional fifth-forces), τ is conformal time, and x_i is the coordinate three-position vector.

In addition to including the massive neutrinos as a new particle family, we also incorporated radiation into the background evolution through the Hubble function, with an additional option to read the Hubble function, H(z), from file. As of now, we have made no changes to how the time-stepping scheme already incorporated in RAMSES operates, which makes for quite slow calculations for the lighter neutrino masses. This could be implemented at a later stage by, for example, detaching the neutrino particles from the grid and allowing them to travel further, letting the CDM particles solely determine the size of the time steps, or by requiring that the neutrino runs use the same refinements and level structure as an equivalent massless neutrino simulation. The latter alternative would allow a more direct comparison between ΛCDM and massive neutrino simulations, especially when studying various ratio properties.

The performance of the ANUBIS code is tested against other codes incorporating massive neutrinos in the massive neutrino code comparison project (MNCCP) within Euclid (Adamek et al. 2022), although with a somewhat lower resolution than what is expected for convergence within 1% between the RAMSES and GADGET-3 codes (Schneider et al. 2016). The upper panel of Fig. 1 displays the neutrino suppression on the CDM + baryon (cb) matter power spectrum for different values of the sum of the neutrino masses at z = 0 for both ANUBIS and GADGET-3. The plot shows good agreement between the two codes, especially for lower neutrino masses. At higher masses, the codes start to deviate at smaller scales. This is most likely due to internal differences in resolution for RAMSES, between the ΛCDM and massive neutrino simulations. As the neutrino mass increases, structure formation is more suppressed, leading to fewer refinements created by the RAMSES AMR-scheme. This results in a slightly lower resolution for the neutrino simulations, compared to the ΛCDM one, which shows up in the ratio. This can be resolved by a higher particle density. The lower panel of Fig. 1 displays the neutrino suppression on the HMF for the same neutrino masses as for the matter power spectrum ratios in the upper panel. The magnitude of the suppression is in agreement for both GADGET-3 and ANUBIS, and increases with the halo mass, corresponding to the findings of Brandbyge et al. (2010).

Fig. 1. Comparison of the neutrino suppression on the cb power spectrum and HMF for ANUBIS and GADGET-3 at z = 0 for various values of the sum of the neutrino masses. Top: the neutrino suppression on the cb power spectrum. The greyed-out area indicates the parts of the ratios lying outside of the Nyquist frequency of the simulations. Bottom: the neutrino suppression on the HMF. Only halos with Npart ≥ 11 are included.

3.2. ISIS

ISIS is a cosmological N-body code incorporating scalar-tensor gravitational theories including screening mechanisms into RAMSES (Llinares et al. 2014). This is done by implementing a non-linear implicit solver for a generic scalar field, which can treat various scalar-tensor-modified gravity theories. This also includes f(R) theories, which can be rewritten into the scalar-tensor format. In particular, ISIS contains Hu-Sawicki f(R)-modified gravity, as described in Sect. 2.3.

The upper panel of Fig. 2 shows the ratio of the cb matter power spectrum for a ΛCDM (lcdm_small) and f(R) (fofr_small) simulation run with ISIS, with the parameter |f_R0| = 10⁻⁴. Here, we can clearly see the enhancement of structure at small scales, with a peak at similar scales to the trough in the neutrino suppression in Fig. 1. For the HMF in the bottom panel of Fig. 2, the amount of large mass halos is increased, compared to ΛCDM, as opposed to the ANUBIS simulations with massive neutrinos that show a decrease. This enhancement is due to a boost of gravity resulting from the fifth force present in f(R)-modified gravity. Summed up, these figures show the opposite behaviour of what is observed in Sect. 3.1, as expected.

Fig. 2. CDM + baryon (cb) matter power spectrum and HMF of the fofr_small simulation, compared to lcdm_small. The power spectrum and HMF ratios are displayed in the top and bottom panels respectively. For the power spectrum ratio, the different σ8-values of the simulations are taken into account through scaling, and the greyed-out area indicates the parts of the ratios lying outside of the Nyquist frequency of the simulations.

The two simulations introduced here, dubbed lcdm_small and fofr_small, are occasionally used as additions to the simulations performed for this work. This is further explained in the following section, along with the simulation details.

3.3. ANUBIS + ISIS = ANUBISIS

To obtain a RAMSES-based code which includes both the effects of massive neutrinos and f(R) gravity, we have merged the ANUBIS and ISIS codes into one code, ANUBISIS. This provides us with the opportunity to run simulations with massive neutrinos and modified gravity both independently and simultaneously. For this paper, we have performed a suite of such simulations with properties as presented in Table 1. Summarised, we have six simulations, one with ΛCDM cosmology, one with Hu-Sawicki f(R)-modified gravity, two with ΛCDM and massive neutrinos, and two with massive neutrinos and Hu-Sawicki f(R) gravity combined. In addition to these six simulations run with the ANUBISIS code, we also have two simulations, one ΛCDM and one f(R), which was previously run with the ISIS code, as presented in Sect. 3.2. These are included as the f(R) simulation was run with |f_R0| = 10⁻⁴, and therefore better demonstrates the effects of modified gravity in specific cases where we are interested in studying this further.

In this section, we present some general properties of the ANUBISIS simulations. To better help distinguish the results, we introduce a general linestyle guide where pure ΛCDM is always shown as a full line, f(R) as a dashed line, massive neutrino runs as a dotted line, and the combination of f(R)-modified gravity and massive neutrinos as a dash-dotted line. The upper panel of Fig. 3 shows the ratio of the CDM+baryon (cb) matter power spectrum between the various ANUBISIS simulations and the ΛCDM case. At large scales, there is originally a slight excess in the power spectrum for the massive neutrino and f(R) + massive neutrino runs. This appears as a result of the grid settings used in the simulations, which is further explained in Appendix A. Essentially, it amounts to a σ₈-scaling (2 − 3% difference), which is accounted for in the figure. On smaller scales, we observe the expected suppression of structure for the massive neutrino simulations and an enhancement of structure for the f(R) simulation. The f(R) + 0.15 eV simulation is slightly suppressed, compared to the pure f(R) run, but it is still dominated by the modified gravity effects due to the low neutrino mass. The f(R) + 0.6 eV simulation, on the other hand, is mostly dominated by the massive neutrinos due to the high neutrino mass, but it is still slightly enhanced, compared to the pure 0.6 eV run.

Fig. 3. CDM + baryon (cb) matter power spectrum and HMF of the different ANUBISIS simulations, compared to ΛCDM. The power spectrum and HMF ratios are displayed in the top and bottom panels respectively. For the power spectrum, the greyed-out area indicates the parts of the ratios lying outside of the Nyquist frequency of the simulations.

In the bottom panel of Fig. 3, the HMF of the various simulations, compared to ΛCDM, is shown. In general, we again see that the massive neutrinos suppress the formation of the most massive halos. This is mostly due to the massive neutrinos reducing the maximum cluster mass on a linear level (Brandbyge et al. 2010). Because of the fifth force contribution to gravity, we again expect the opposite effect for the f(R) simulation. This is, however, not very clear for the |f_R0| = 10⁻⁵ simulation. If we instead look back at the smaller box simulations from ISIS in Fig. 2 with |f_R0| = 10⁻⁴, this effect is much more prominent.

Another interesting property is the scale-dependent halo bias. Here, we define it as the ratio between the cross-power spectrum of halos and CDM+baryons and the auto matter power spectrum of CDM+baryons,

$$\begin{array}{c}\hfill b(k)=\frac{P{(k)}_{\mathrm{halo},\mathrm{cb}}}{P{(k)}_{\mathrm{cb}}}.\end{array}$$(31)

This is defined by the cold species instead of the total matter, as it yields a more universal and scale-independent result in the presence of massive neutrinos (Castorina et al. 2014). In Fig. 4, we see the halo bias for the ANUBISIS simulations at the top, and the ratio of the biases, compared to ΛCDM, at the bottom. The ratios show a bump at the same scales where we see a trough for the massive neutrino power spectrum ratios. In our case, this shows that our halo power spectrum is less sensitive to the neutrino mass than the cb power spectrum, as discussed by Hassani et al. (2022). In Fig. 5 we also show the scale-dependent bias for the pure ISIS simulations. In general, we see that the halo bias increases with neutrino mass and decreases with the f_R0 parameter, compared to ΛCDM, in line with previous findings in the literature (e.g., Arnold et al. 2019; Chiang et al. 2019).

Fig. 4. Scale-dependent halo bias for the various ANUBISIS simulations. The upper panel shows the value of the bias for each simulation, while the lower panel shows the ratio with ΛCDM for each model. The greyed-out area indicates the parts lying outside of the Nyquist frequency of the simulations.

Fig. 5. Scale-dependent halo bias for the fofr_small and lcdm_small simulations. The upper panel shows the value of the bias for each simulation, while the lower panel shows the ratio with ΛCDM. The greyed-out area indicates the parts lying outside of the Nyquist frequency of the simulations.

Looking at the bias values on linear scales, we can calculate an estimate of the linear halo bias. This is presented in Table 2, together with β and σ₈. The two latter are simulation parameters, β being the reconstruction parameter of Sect. 2.1.3 and σ₈ the amplitude of the linear matter power spectrum at a scale of 8 h⁻¹Mpc. For the ΛCDM and massive neutrino simulations, σ₈ is a known parameter obtained from the linear power spectrum calculated by CLASS⁴ as a part of the initial condition set-up (see Sect. 3.4). For the simulations including modified gravity, an estimation of the linear σ₈-value was found by observing that ${\sigma }_{8,\mathtt{CLASS}}^{\nu }/{\sigma }_{8,\mathtt{CLASS}}^{\mathrm{\Lambda }\mathrm{CDM}}\approx {\sigma }_{8,\mathtt{ANUBISIS}}^{\nu }/{\sigma }_{8,\mathtt{ANUBISIS}}^{\mathrm{\Lambda }\mathrm{CDM}}$ for the linear values calculated from CLASS and the non-linear values calculated from ANUBISIS. We then assumed that this relation also holds in the modified gravity context.

3.4. Initial conditions

The initial conditions used for the ANUBISIS simulations were generated by following the same procedure as outlined in Adamek et al. (2022). In short, the linear matter power spectra and transfer functions were generated by CLASS (Lesgourgues 2011; Blas et al. 2011), and rescaled to z = 127 through the REPS⁵ code (Zennaro et al. 2017), which takes the effects of massive neutrinos into account through a two-fluid description and by including radiation in the background evolution and implementing a scale-dependent growth rate. Positions and velocities for the CDM and neutrino particles were then generated by a version of the N-GenIC code⁶ that also has been modified to include the scale-dependence of the growth rate and growth factor in the presence of massive neutrinos.

For the initial conditions, we used Ω_b = 0.049 and Ω_c ≃ 0.27. We kept the total matter density fixed at Ω_m ≃ 0.319 by adjusting the ratio of Ω_c and the neutrino density parameter, Ω_ν = ∑m_ν/(93.14h² eV). Besides this, we used the Hubble constant, h = 0.67, the scalar spectral index, n_s = 0.9619, the CMB temperature, T_CMB = 2.7255 K, and the amplitude of the primordial power spectrum, A_s = 2.215 × 10⁻⁹, at the pivot scale k_p = 0.05 Mpc⁻¹. This is the same cosmology as applied by Adamek et al. (2022).

We used the method outlined above to generate initial conditions for our ΛCDM and massive neutrino simulations. For the simulations including f(R)-modified gravity, we used the outputs from the other simulations at z = 4 as initial conditions. The deviations from GR are small at this stage, and only become more important at later times, as demonstrated for z = 3 and |f_R0| = 10⁻⁵ by Zhao et al. (2011). For the two independent ISIS runs, the initial conditions were generated by Grafic2 (Bertschinger 2001), with the parameters Ω_c = 0.267, Ω_b = 0.045, h = 0.719, and n_s = 1.0.

4. Method

In this section, we detail the various codes and packages used to study our simulation data. We also provide an overview of the full analysis process.

4.1. Halo finder

To identify halos in our simulations, we used the ROCKSTAR⁷ halo finder (Behroozi et al. 2012). ROCKSTAR locates halos in phase space by applying the 3D friends-of-friends (FOF) method to pinpoint overdense regions. For each of the groups created in the FOF procedure, the linking length (the characteristic length scale grouping the particles together) is reduced progressively, so that subgroups emerge in a hierarchical structure. Seeds are then placed in the lowest substructures and particles are assigned to the halo seed within the shortest phase-space distance. From this, the relationship between host and subhalos is computed and unbound particles are removed. Finally, halo properties are calculated.

When applying ROCKSTAR to our simulations, we read the gadget files directly to obtain the CDM particle properties and the cosmological parameters. The force resolution was set to be approximately the smallest distance resolved by the simulation, and the linking length was set to l = 0.28 h⁻¹Mpc, which is 0.2 of the mean inter-particle separation. In the end, we removed the subhalos and were left with a host halo catalogue. As shown in Adamek et al. (2022), ANUBIS underestimates the number of halos, especially at low mass. This can be improved upon by increasing the particle density but was not done for the simulations presented in this paper. Because of this, we made no further mass cuts to the halo catalogue.

4.2. Void finder

When identifying voids within the halo catalogue of a simulation, we must first define a void numerically. In this paper, we apply two different void definitions, spherical voids and voxel voids.

For the spherical void detection, we used the void finder module implemented in Pylians⁸, which is based on the method detailed in Banerjee & Dalal (2016). Here, the algorithm is provided with a range of radii corresponding to spherical regions of various sizes to be investigated, along with an overdensity threshold given by $\mathrm{\Delta }=\rho /\overline{\rho }-1$. Starting with the largest radius, the simulation box is divided into a grid and the density profile is calculated in each voxel and smoothed on a scale corresponding to the given radius. Voxels with densities below the threshold are identified and sorted. Starting with the lowest density voxel, all the voxels around it within the given radius are considered and added as a part of the void. This is then repeated for the next to lowest density voxel and so on. If any of the voxels within the corresponding radius are already assigned to another void, the new void is rejected. Once this is done for all the voxels below the density threshold, the algorithm moves on to the next to largest radius provided, smooths the density field at that scale, and once again proceeds as explained above. This is repeated for all the radii given initially, and the resulting voids identified are spherical regions of these specific sizes, with the centre at the lowest density voxel. In our case, we set Δ = −0.75 and provided the algorithm with 47 void radii between 16 − 82 h⁻¹Mpc.

For detecting voids by using the voxel void definition, we employed the Revolver⁹ void finder (Nadathur et al. 2019c). Here, the voxel voids are defined by using a watershed method. The simulation box is divided into a grid and local minima of the density field are located. Around each of these minima, surrounding voxels with increasing overdensity are added to the pool up until a voxel with a lower overdensity than the one previously added is discovered. Using this definition, the identified voids may have any shape, as opposed to the spherical void definition. The centre of the void again lies within the voxel with the lowest density. The Revolver void finder can also perform reconstruction when provided with a halo catalogue given in redshift space. This is based on linear theory, as described in Sect. 2.1.3, and requires the halo bias, b, and growth rate, f, as input parameters. For a more detailed description of the spherical and voxel void definitions, along with other void-identifying algorithms, see Massara et al. (2022).

4.3. Covariance matrix

Each of our simulation boxes only has one realisation. To obtain an estimate of the uncertainty in the redshift space void-halo CCF calculated from the simulation data, we attained the covariance matrix through a Jackknife method. We applied the Jackknife estimator as implemented in Pycorr¹⁰, following Mohammad & Percival (2022).

The ANUBISIS simulation boxes were equally divided into n = 512 sub-boxes¹¹. Individual Jackknife realisations were then made by calculating the correlation function for the full volume with one of the sub-boxes removed at a time. This lead to n Jackknife realisations, each with a volume fraction (n − 1)/n of the original volume. Based on this, each element in the covariance matrix, C_ij, is given by

$$\begin{array}{c}\hfill {\mathsf{C}}_{\mathit{ij}}=\frac{n-1}{n}{\displaystyle \sum _{k=1}^n}[{\xi }_i^k-{\overline{\xi }}_i][{\xi }_j^k-{\overline{\xi }}_j],\end{array}$$(32)

where

$$\begin{array}{c}\hfill {\overline{\xi }}_i=\frac{1}{n}{\displaystyle \sum _{k=1}^n}{\xi }_i^k\end{array}$$(33)

is the mean estimate from all the n Jackknife realisations.

It is important to note that since we only have one simulation of each individual case, and we used the data from the simulations to calculate ξ^r, v_r(r), and σ_v∥, the model and data ξ^s are correlated. Ideally, the input parameters should be calculated from the average of many mock simulations, as is done in for example Woodfinden et al. (2022) and Nadathur et al. (2019b), to avoid this problem, along with the issue of the uncertainty associated with only having one measurement of ξ^r. One way to deal with the correlated errors is to calculate the covariance matrix for the CCF of the difference between the model and the data, ${\xi }_{\text{diff}}^s={\xi }_{\text{model}}^s-{\xi }_{\text{data}}^s$. A more detailed explanation and demonstration of the effects of this approach can be found in Appendix A of Radinović et al. (2023). We did not include this step in our analysis as we focus more on the comparison between the simulations as opposed to reducing the statistical error for a single simulation.

4.4. Analysis pipeline

For the ANUBISIS simulation data, the CDM and neutrino particles were given separate identifiers. It is the CDM particle information that goes into the procedure detailed below unless otherwise stated.

First, the particle data were given to the ROCKSTAR halo finder, which identified halos in the simulation boxes as detailed in Sect. 4.1. This provided both the positions and velocities of the halos, enabling us to put the halos into redshift space if needed. The halo catalogues were then given to the Revolver void finder, which identified voxel voids in between the halos, as explained in Sect. 4.2. In order to test how well reconstruction (Sect. 2.1.3) performs in the case of massive neutrinos and modified gravity, this step was not only performed directly on the real space halo catalogue but also with the redshift space halo catalogue as input. We performed the reconstruction step by applying Revolver, which solves Eq. (11) through an iterative fast Fourier transform algorithm (Burden et al. 2014, 2015). This is the same algorithm used for BAO reconstruction (e.g., Gil-Maŕin et al. 2020; Bautista et al. 2021). We then had two void catalogues, one directly identified in real space, and one in reconstructed real space, both found by the voxel void definition. We then performed the same void finding with the spherical void definition, as detailed in Sect. 4.2, both with the real space and reconstructed real space halo positions. After obtaining void catalogues in real and reconstructed real space, for the voxel and spherical void definitions, we performed a radius cut at 21 h⁻¹Mpc for all catalogues, to eliminate void discoveries that coincide with the average spacing between halos. For each catalogue, the remaining voids with their coinciding halos were then stacked on top of each other to obtain an ‘average’ void with more statistically robust properties and spherical symmetry. From the halo distribution in and around these voids, we calculated the density profile, mean velocity profile, and LOS velocity dispersion. The latter was computed as detailed by Radinović et al. (2023). These are all ingredients needed to calculate the void-halo CCF in redshift space as shown in Eq. (5).

We want to compare the modelled void-halo CCF in redshift space with the void-halo CCF calculated directly from the simulation data. In order to obtain the latter we used the correlation function estimation wrapper, Pycorr, which currently utilises the Corrfunc¹² CCF engine (Sinha & Garrison 2019, 2020). Provided the redshift space halo catalogue, real space void catalogue, and randomly distributed reference catalogues, Pycorr calculates the void-halo CCF of the simulation data as we would observe it. This is done using the Landy-Szalay estimator (Landy & Szalay 1993),

$$\begin{array}{c}\hfill \xi (r,\mu )=\frac{D_1{D}_2(r,\mu )-D_1{R}_2(r,\mu )-D_2{R}_1(r,\mu )+R_1{R}_2(r,\mu )}{R_1{R}_2(r,\mu )},\end{array}$$(34)

where D is the simulated data catalogues, R is the random catalogues, and 1 and 2 are the halos and voids. We used 50 bins between 0 − 150 h⁻¹ Mpc for r, the spatial distance between the pairs, and 200 bins from −1 to 1 for μ, the angular separation. For the random catalogues, we used 50 times as many halos and voids as the simulations. Pycorr also has an inbuilt Jackknife estimator, which allowed us to calculate the covariance matrix of the CCF, as described in Sect. 4.3. We calculated the void-halo CCF using both the voids identified in real space and in reconstructed real space.

Having the simulated CCF, we used Victor¹³ to calculate the theoretical model and compare it to the data, both in the real and reconstructed case. To do so, Victor requires ξ^s from the simulations, a covariance matrix, and ξ^r, v_r(r), Δ_DM, and σ_v∥ as input to the model. Victor also assesses the goodness of fit between the data and the model, and provides a χ² value upon request.

Through an interface with Cobaya (code for bayesian analysis¹⁴, Torrado & Lewis 2019, 2021), we can also use Victor to perform MCMC fits of the parameters in the void-halo CCF model. We assumed a Gaussian form of the likelihood

$$\begin{array}{c}\hfill log\mathcal{L}=-\frac{1}{2}({\xi }_{\mathrm{model}}^s-{\xi }_{\mathrm{data}}^s){\mathsf{C}}^{-1}({\xi }_{\mathrm{model}}^s-{\xi }_{\mathrm{data}}^s{)}^{\mathrm{T}},\end{array}$$(35)

along with flat priors for the parameters fσ₈, σ_v∥, β, and ϵ, encompassing the fiducial values. In the case where we used the actual real space data from the simulations, β was not included. In the case of reconstruction, β was allowed to vary. When comparing CCF data to theory, reconstruction was only performed for the fiducial β-value of each simulation. However, to allow for β to vary for the MCMC fits, we performed reconstruction, void finding, and the CCF calculation repeatedly for 11 different β-values, which were provided to Victor to perform the fit. Ideally, covariance matrices should also be computed in each reconstruction case, as it depends on β. This is a time-consuming task, and as a first approach, we instead assumed a fixed β = β^fiducial for all the covariance matrices. We also kept the input ξ^r equal to the one calculated for β^fiducial in all cases, meaning that the β-dependence only shows up in ξ^s. This is the same approach taken by Radinović et al. (2023).

5. Results

In this section, we present the results of our analysis. In addition, we discuss the implications of our results.

5.1. Void abundance

First, we take a look at the general void population identified in all the ANUBISIS simulations. In Fig. 6 we show the abundance of voids depending on the effective void radius. At the top, the voxel void definition has been used to find the voids, and at the bottom, the spherical void definition. For the voxel voids, the difference in the number of voids is clearly apparent in the approximate range 10 − 40 h⁻¹ Mpc. For the spherical voids, however, this is not the case. If we zoom in, the same ordering is present, but this is not apparent at first glance due to the high number of voids identified in the given range. The spherical void finding algorithm, as detailed in Sect. 4.2, identifies voids by smoothing the field for a given top-hat radius and declaring spheres with density below a certain threshold as voids. For small radii, this results in a large amount of voids, some of which might be the result of shot noise, for all simulations.

Fig. 6. Histogram of effective void size for the various ANUBISIS simulations. The voxel voids display a notable difference in the number of voids in the ∼10 − 40 h−1 Mpc range. Both void definitions also show a slight opposite abundance difference amongst the largest voids in the simulation catalogues.

For the f(R) and f(R) dominated f(R)+0.15 eV simulations, we see an increase in the number of voids within the 10 − 40 h⁻¹ Mpc effective radius range, compared to ΛCDM. This is due to the fifth force enhancing gravity in these regions, resulting in a more effective ‘evacuation’ of the areas. Li et al. (2012) also report a higher number of the large voids in their f(R) simulations, although it should be noted that their maximum void size is around 15 h⁻¹ Mpc due to their smaller, but higher resolution, simulation boxes.

The most massive neutrino simulation, 0.60 eV, shows the opposite behaviour. Here, the amount of voids in the given size range is suppressed, compared to ΛCDM. This is a result of the neutrinos slowing down the clustering and thereby the evolution of the voids towards lower densities. Massara et al. (2015) also report fewer large voids in their massive neutrino simulations, compared to ΛCDM, where large in their case is in the range 20 − 40 h⁻¹ Mpc. For very large voids, R_eff ≳ 45 h⁻¹Mpc, we see a turn-around of the ordering for both void definitions. Cai et al. (2015) also found a higher amount of very large radii (≳25 h⁻¹ Mpc) voids in their regular GR simulations, compared to f(R). Based on investigation they suggest that this is a result of the largest voids being less empty in the f(R) case, as the enhanced gravity inside makes it easier for halos to form, compared to ΛCDM. When looking at a big region in ΛCDM versus a big region in the f(R) simulation, the number density of halos would be lower for ΛCDM, making it easier to pass the void identification criteria. The opposite could then be argued for the massive neutrino simulations. It should be noted that the differences observed for very large voids in Fig. 6 are enhanced due to the logarithmic scale, and could also be affected by the small sample size.

For the ΛCDM simulation, we find in total N_halo ≈ 1.3 × 10⁶. From this, we chose a radius cut for the voids as $R_{\mathrm{eff}}^{\mathrm{cut}}=1.5{(N_{\mathrm{halo}}/V_{\mathrm{box}})}^{-1/3}\approx 21h^{-1}\mathrm{Mpc}$. This was applied to the void catalogues to make sure that we exclude what are simply empty regions in between halos in the simulations, and not actual voids.

5.2. Velocity profile

When using galaxy surveys to identify voids, the velocity profile needed for the CCF model is typically modelled by the linear velocity profile or estimated from simulations. Because of this, we want to study how the simulated velocity profiles for various cosmologies compare to the linear model. In Fig. 7, the mean radial outflow void velocity profile for each ANUBISIS simulation is shown for both the voxel and spherical void definitions. The individual profiles are compared to the theoretical linear velocity profile, as shown in Eq. (8), and the more general velocity profile, as shown in Eq. (10). For the linear velocity profile, the growth rate, $f={\mathrm{\Omega }}_m^{0.55}$, corresponding to the expected ΛCDM value, was applied for all simulations. In doing so, we ignore the scale dependence of the growth rate in the f(R) and massive neutrino cosmologies and investigate whether or not this leads to biased results. For the general profile, the best-fit values of the A_v-parameter are given in each panel of Fig. 7. These were obtained through the least squares method.

Fig. 7. Radial void velocity profiles for all ANUBISIS simulations together with the linear velocity profile and a fit to the more general velocity profile, presented in Sect. 2.1.2. The voxel and spherical void definitions are shown in the left and right columns respectively. The Av-parameter in each panel shows the best-fit value of the general profile.

The general velocity profile is, with the best fit A_v-parameter values, by construction, a good match to the velocity profiles found in the simulation data. For the voxel voids, the best-fit parameter value is consistently higher for f(R) and decreases with increasing neutrino mass, compared to ΛCDM. This is expected as a result of higher velocities in the f(R) case, and lower in the massive neutrino case. The mixed simulations lie somewhere in between. For the spherical voids, we do not see the same pattern for the neutrino masses and the value of A_v. However, for the spherical voids, the general velocity profile does not fit as well as for the voxel voids at higher radii. When estimating A_v, this must be taken into account, in addition to the peak of the velocity profile. If we could observe the velocity profile directly over a large void sample, the best-fit A_v-value could be an interesting parameter to explore as a cosmological probe. Observing the velocity profile directly is not possible in traditional galaxy surveys, which is why it is currently modelled or estimated from simulations.

From visual inspection, it is clear that the theoretical linear velocity profile is not a good fit for the data close to the void centre. This is due to a sparse tracer sample in this region and also depends on the void centre definition (Massara et al. 2022). It affects the voxel void definition more than the spherical one. In the voxel void case, the linear model overestimates the clustering. We can see that this is somewhat compensated for in the f(R) simulations, due to higher velocities. This is even more visible for the pure ISIS simulations, as shown in the lower left panel of Fig. 8. For the spherical voids, this is less of an issue, although the linear model slightly underestimates the clustering. In the latter case, this could be compensated by adopting a different value of the growth rate, f, in the linear theory model, as demonstrated in Fig. 9. For the voxel voids, increasing the value of f in line with expectations (e.g., Mirzatuny & Pierpaoli 2019) only increases the discrepancies.

Fig. 8. Radial void velocity profiles for the fofr_small and lcdm_small simulations together with the linear velocity profile and a fit to the more general velocity profile. The voxel and spherical void definitions are shown in the left and right columns respectively. The Av-parameter shows the best-fit value for the general profile.

Fig. 9. Linear velocity profile, compared to the fofr_small simulation, for various values of the growth rate, f. The regular linestyle and colour show the velocity profile for the simulation. For ΛCDM we have $f={\mathrm{\Omega }}_m^{0.55}$ ≈ 0.53.

Altogether, it is evident that the linear velocity profile is not a drastically worse fit for the f(R) or massive neutrino simulations. Before we can address any differences due to cosmology, however, we must make sure that a sparse tracer sample does not affect our model-data comparison. As it stands now, applying the linear velocity profile to the CCF modelling results in modifications similar to changes in fσ₈, due to the offset induced by the number of tracers (Massara et al. 2022). In addition, we should point out that even if the velocity modelling was improved to forgo this issue, the modifications expected from the |f_R0| = 10⁻⁵f(R) simulation and the most massive neutrino simulation are mostly on scales smaller than the average void size.

5.3. Velocity dispersion

The CCF theory presented in Sect. 2.1.1 depends on the velocity dispersion of the void velocity profile. For all the ANUBISIS simulations, Fig. 10 shows this quantity for both the voxel and spherical void definitions, scaled by the average amplitude beyond 75 h⁻¹Mpc, σ_*. The velocity dispersion profile has the same shape for all the various cosmologies, but a different amplitude, σ_*, as presented in Table 3. The relative values of the amplitudes coincide with expectations. The halos in and around the void in the f(R) simulation feel a stronger gravitational pull towards the edge, compared to the ΛCDM case. In general, there is also more clustering which leads to higher velocities. The relative values of the amplitudes in the f(R) and ΛCDM cases coincide well with the findings of Fiorini et al. (2022). The opposite behaviour is seen for the massive neutrino simulations, where less clustering results in lower velocities. However, although these differences can be calculated from the simulated data, they can not readily be observed.

Fig. 10. LOS velocity dispersion scaled by σ* as given in Table 3. The voxel and spherical void definitions are shown in the upper and lower panels respectively. All the simulated cosmologies result in similar profile shapes. The largest differences are located towards the void centres, which are dominated by shot noise.

5.4. Void-halo CCF: data versus model

We are interested in the void-halo CCF of the various simulated cosmologies and how well the theory presented in Sect. 2.1.1 reproduces the simulation data. Figure 11 shows the monopole, quadrupole, and hexadecapole of the void-halo CCF for the ANUBISIS simulations. The left and right columns display the results for the voxel and spherical void definitions respectively. The voids were identified from the real space halo catalogue, and the redshift space halo positions were calculated directly from the real space positions and corresponding LOS velocities.

Fig. 11. Void-halo CCF for the ANUBISIS simulation data. The voxel and spherical void definitions are shown in the left and right columns respectively. The error bars are estimated using a Jackknife technique and are only shown for the ΛCDM case for a tidier visual representation of the data.

For the monopole, the differences between the simulations are clearer for the voxel void definition, as was also the case earlier when looking at the void abundance. The shape of the monopole is closely related to the shape of the density profile around the void, as both in essence describe how halos are distributed around the void centre as a function of radius. The f(R) and f(R)+0.15 eV simulations show the narrowest profiles (most evolved) and the 0.6 eV simulation the broadest (least evolved), while the others fall somewhere in between. From the void abundance histogram in Fig. 6, we saw that the f(R) case has more voids in the 10 − 40 h⁻¹ Mpc effective radius range and the 0.6 eV has less, compared to ΛCDM. This further explains what we see here in terms of the shapes of the profiles, as the monopole is based on the overall averaged void from each simulation box. For the voxel void definition, the central value goes towards −1.0, while for the spherical void definition, we see a very flat core stabilising at around −0.88. This is again a result of the void-finding algorithm, where a sphere of a given radius with a smoothed density below −0.75 is classified as a void. As mentioned before, this leads to the detection of a large amount of small voids, and in addition, most of the detected voids have central densities that are larger than −1.0.

The quadrupole is most easily studied for the voxel void definition. Here, we see, in all cases, the expected smooth shape resulting from RSDs as a consequence of halos inside the void moving towards the small overdensity at the edge, but also halos outside the void moving back towards the void edge. For the spherical voids, we instead see a sharp feature proceeding the edge of the void. The difference between the quadrupoles of the two void definitions is due to the shapes of the void density profiles and the presence of a velocity dispersion. This is further explained through the use of a toy model in Appendix B. Within the error bars, the quadrupole for each simulation is not easily distinguishable. The most massive neutrino simulation shows signs of lower peaks due to overall lower velocities for the voxel void definition, but higher resolution simulations should be performed in order to confirm this.

The hexadecapole is, from the theoretical model, expected to be small (Cai et al. 2016). This is the case for all simulations and both void definitions, except for a little dip in the voxel void ${\xi }_4^s$ around s ∼ 35 h⁻¹ Mpc. This small signal in the hexadecapole shows why it can still be important to include this multipole when fitting the model to the data to obtain cosmological parameters, instead of assuming it to be zero.

Figures 12–14 respectively display the monopole, quadrupole, and hexadecapole of the void-halo CCF in redshift space as calculated directly from the ANUBISIS simulation data, compared to its theoretical value calculated by Victor as outlined in Sect. 2.1.1. The simulation data follow the regular colour and linestyle pattern, while the theory is presented with full brown and blue lines. The brown line is the relevant theoretical multipole as calculated by Victor, with the velocity profile from the simulation data given to the model as a velocity template, and the blue line is the same, only with the linear velocity profile as displayed in Eq. (8) given as input. For each simulation, void definition, and multipole, only the ΛCDM simulation data line has error bars. This is to illustrate their magnitude but otherwise reduce cluttering.

Fig. 12. Difference between the redshift and real space monopole of the void-halo CCF from the ANUBISIS simulations vs. theory calculated by Victor. The results are obtained by using the voxel and spherical void definitions in the left and right columns respectively. Two theory results from Victor are included, one where the velocity profile is given by linear theory and one where it is given by a template. The template used in this case is the velocity profile calculated from the simulation data.

Fig. 13. Quadrupole of the void-halo CCF from the simulations vs. theory calculated by Victor. Layout and input as explained for Fig. 12.

Fig. 14. Hexadecapole of the halo-void CCF from the simulations vs. theory calculated by Victor. Layout and input as explained for Fig. 12.

The monopole data-theory comparison in Fig. 12 is presented through the difference between the redshift space and real space monopole. This is done to better display the discrepancies between the data and the two models. For all simulations and both void definitions, we see good agreement. There is a preference towards inputting the velocity profile from the simulations in the model, compared to using the linear velocity profile. This is expected simply because we are providing the model with the actual data that we want to reproduce. We already saw, in Sect. 5.2, that the linear velocity profile does not reproduce the velocity profile directly calculated from the halos in the simulations very well. However, if this study was performed on observational data, we would not have had access to the exact mean outflow velocity profile. We could then either use the linear approximation or make estimates through simulations.

The quadrupole data-theory comparison in Fig. 13 is the most interesting to study, as it is the RSDs observed here that, through modelling, can give us an estimate of the fσ₈-value. It is again evident from visual inspection that the model reproduces the data well for all simulations. For the voxel void definition, applying the linear velocity profile as an approximation gives results within the error bars for all simulations. For the spherical void definition, the consequence of applying the linear velocity profile is more significant, even though the difference between the data and linear theory velocity profiles in Fig. 7 is larger for the voxel void definition. This could be partly because even though the linear and simulation velocity profiles deviate more towards the void centre for the voxel voids, compared to the spherical voids, they actually coincide better for r ≳ 30 h⁻¹Mpc. This is also where we see some of the larger discrepancies between the spherical void simulation quadrupole and the modelled quadrupole with the linear velocity as input. Even though we see a difference in the model fit for the two void definitions, there is no significant difference in the performance of the model between the various simulated cosmologies. Ideally, we need higher resolution simulations along with more accurate velocity modelling in order to quantify this.

Figure 14 shows that also the hexadecapole is well described by the model for both void definitions. This further encourages including it in analysis when it is available.

5.5. Reconstructed void-halo CCF: data versus model

The reconstruction process, detailed in Sect. 2.1.3, assumes linearity and a constant growth rate which might affect its performance when modified gravity or massive neutrinos are introduced. Figure 15 shows how reconstruction performs for each individual ANUBISIS simulation. The lines, following the regular colour and style pattern, show a histogram of the difference between the LOS coordinate (defined as the z-direction in our simulation boxes) in real space and reconstructed real space. The fainter line, in the same colour, shows a histogram of the difference in the LOS coordinate in real space and redshift space, before any reconstruction is performed. The reconstruction process illustrated in this figure is executed with the fiducial β-values presented in Table 2.

Fig. 15. Difference in the LOS coordinate between real space and reconstructed real space for the ANUBISIS simulations shown in the regular linestyle and colour. The weaker full line of the same colour shows the difference between the LOS coordinate in real and redshift space.

For each simulation, the histogram of the difference between real and redshift space is always broader than that for the difference between real and reconstructed real space. This tells us that performing reconstruction has indeed led us closer to the actual real space values. The method is tested on ΛCDM simulations by, for example, Nadathur et al. (2019a), Woodfinden et al. (2022) and Radinović et al. (2023), and we achieve similar results in our case. The non-ΛCDM simulations again obtain comparable results to the ΛCDM reference. Although the f(R) dominated ones (f(R) and f(R)+ 0.15 eV) seem to perform slightly worse, and the heaviest neutrino mass, 0.6 eV, slightly better. The reconstruction method is based on linear theory, and one possible explanation for this could therefore be that the massive neutrino simulations, with their suppressed clustering and lower velocities, are closer to the linear approximation on the scales where the reconstruction is performed. For the f(R) dominated simulations, the opposite would then be the case. This is further highlighted for the pure ISIS simulations, as shown in Fig. 16, where the deviations from GR are larger. As a test, we also ran the reconstruction process for this simulation when varying the growth rate, f, by 10 − 20% (as for Fig. 9). This did not lead to a significant difference in the histogram, but it did show up in the quadrupole as expected, and explained, by Nadathur et al. (2019a).

Fig. 16. Difference in the LOS coordinate between real and reconstructed real space and real and redshift space for the fofr_small and lcdm_small simulations.

Figure 17 displays the quadrupole of the void-halo CCF for both the simulation data and the theory when the reconstruction step has been performed. In practice, this means that we performed the analysis as before, but instead of using the actual real space positions to calculate ξ^r, we instead used the reconstructed real space information. This mimics a possible procedure for performing the quadrupole analysis when using observational data. The model does, however, also require the density profile of the CDM, the mean outflow velocity profile, and the velocity dispersion. For this, we used the same as before, which was calculated from the actual real space data. This works as an approximation, as these quantities are only templates given to the model. For the velocity dispersion, it was also necessary to use the real space simulation data, as we did not have velocity information for the halos in reconstructed real space.

Fig. 17. Quadrupole of the reconstructed void-halo CCF from the simulations vs. theory calculated by Victor. The real space CCF used as theory input is calculated from the reconstructed real space halo catalogue and the voids identified in it. The layout follows Fig. 12.

Comparing Figs. 13 and 17, it is clear that the model better fits the data when the reconstruction step is not involved. This is expected, as reconstruction only approximates the real space halo positions, and can not fully recreate the real space CCF. The deviation between the model and data does, however, appear to be comparable for the different ANUBISIS simulations.

5.6. MCMC fits

A way to test the CCF theory for the different simulated cosmologies is to check if they all reproduce the original parameter values of the simulations in a similar manner. Figure 18 displays the posterior distributions of the fσ₈ and ϵ parameters based on the void-halo CCF simulation data and model calculations for all ANUBISIS simulations. Each triangle plot shows the result for one simulation and both void definitions and follows the regular colour and linestyle pattern. The contours display the 68% and 95% confidence intervals and the grey dashed lines show the fiducial values, although with an assumed value of $f\approx {\mathrm{\Omega }}_m^{0.55}$ for all simulations. We only display the fσ₈ and ϵ parameters, as the LOS velocity dispersion, σ_v∥, is not an observable. The corresponding parameter values are listed in Table 4.

Fig. 18. MCMC fits for the ANUBISIS simulations. Each triangle plot shows the fσ8 and ϵ estimates of one simulation for both the voxel and spherical void definitions. The grey dashed lines show the fiducial values and the contours display the 68% and 95% confidence intervals.

It is evident from Fig. 18 that the voxel voids more closely reproduce the fiducial parameter values overall. This could possibly be because the voxel void definition results in a smoother CCF that is easier to fit than the very sharply varying spherical void CCF. We do, however, see that the model prediction is able to follow the sharp features (e.g., Fig. 13). The spherical void catalogue also contains a larger amount of voids with effective radii at the lower end of the included radius range, compared to the voxel void catalogue (Fig. 6). It should be further explored whether or not performing the fits in different radius bins would alleviate the discrepancies between the spherical and voxel voids’ abilities to reproduce the fiducial parameter values. This investigation is saved for future work.

In general, we find that the spherical voids lead to higher fσ₈-values, compared to the voxel voids. This might be due to a large number of small voids with deep density profiles identified by the spherical void finder. This gives, on average, a higher velocity peak, as can be observed in Figs. 7 and 8, and could lead to larger inferred values of fσ₈. The ϵ-values, on the other hand, seem to be quite consistent between the two definitions for all simulations. The fact that the voxel and spherical voids coincide the best in their prediction for the f(R)  +  0.15 eV simulation is most likely a coincidence, which illustrates how the different model components and void definitions can vary. In general, even though there is a difference between the void definitions, there does not seem to be a clear distinction for how well the fiducial values are reproduced for the various simulated cosmologies.

In Fig. 19, we see the results of the MCMC fits when the reconstruction step is included and β is allowed to vary. The numerical values are given in Table 4. Again, the voxel void definition better reproduces the fiducial values of the ANUBISIS simulations. In fact, the differences between the voxel and spherical void definitions have increased. We still see quite consistent ϵ-values between both void definitions, with a slightly lower value for the spherical voids except for the 0.15 eV and 0.6 eV cases. The fσ₈-values, however, are now consistently lower for the spherical voids, compared to the voxel voids, opposite of before. This might lie in the ability of the reconstruction method to correctly predict the high velocities found in the smaller spherical voids. Yet again, there is no clear difference between the simulated cosmologies and how accurately the fiducial parameter values are recovered. This underlines that the reconstruction method should be further studied for various void definitions.

Fig. 19. MCMC fits for the ANUBISIS simulations with reconstruction. Each triangle plot shows the fσ8 and ϵ estimates of one simulation for both the voxel and spherical void definitions. The grey dashed lines show the fiducial values and the contours display the 68% and 95% confidence intervals.

It should be noted that for the fits that include reconstruction and a free β-value, we encounter the issue of a multimodal likelihood surface for β. This issue is more prominent for the voxel void definition than for the spherical one. It does, however, seem to be an artefact of the reconstruction model and, as a first approach, we dealt with it by constraining the range of the flat prior. The new range was chosen to include the most likely values resulting from the parameter fit performed with the original prior.

6. Conclusions

Cosmic voids are promising and independent probes of gravity and cosmology expected to provide stringent constraints on fσ₈ and the ratio between the Hubble distance and transverse comoving distance in upcoming space missions such as Euclid (Hamaus et al. 2022; Radinović et al. 2023). Voids are also proposed as grounds for studying the effects of massive neutrinos and modified gravity, both separately and simultaneously, and as a way to break known degeneracies (Li et al. 2012; Massara et al. 2015; Voivodic et al. 2017; Kreisch et al. 2019; Perico et al. 2019; Contarini et al. 2021). In this paper, we investigate the performance of models describing the void-halo CCF, the radial mean void velocity profile, and a reconstruction method using ANUBISIS simulations with both massive neutrinos and Hu-Sawicki f(R) gravity. Our simulation suite consists of six simulations, one reference ΛCDM cosmology, one with Hu-Sawicki f(R)-modified gravity (|f_R0| = 10⁻⁵), two with massive neutrinos (∑m_ν = 0.15 eV,  0.6 eV), and two with both massive neutrinos and modified gravity (|f_R0| = 10⁻⁵ + ∑m_ν = 0.15 eV,  0.6 eV). Occasionally, we supplement with a ΛCDM and f(R) simulation (|f_R0| = 10⁻⁴) run with ISIS, when we want to investigate the effects of modified gravity further.

We find that the linear void velocity model (Eq. (8)) fits the velocity profiles calculated directly from the data similarly for all ANUBISIS simulations. This is the case both when using the voxel void and spherical void definitions to locate voids in the halo catalogues. In fact, the increased velocities in the f(R) dominated simulations happen to partially compensate for the bias between the velocity profile calculated from the data and the one calculated from linear theory, arising from a sparse tracer sample (Massara et al. 2022). This shows that the linear velocity model needs to be improved further before differences due to massive neutrinos or modified gravity can be seen clearly. We also fit the velocity profile from the simulation data to a more general velocity model (Eq. (10)), with one free parameter. This model fits the data well and has a different best-fit parameter value for the various simulations. If the radial velocity profile could be observed directly with good precision from galaxy surveys, the value of this free parameter would be an interesting observable. The value is, however, degenerate for massive neutrinos and f(R)-modified gravity. Similarly, we find that the LOS velocity dispersion profiles have the same shape for all the simulated cosmologies, but different amplitudes. This amplitude is an interesting parameter, but it has the same issue of degeneracy, in addition to not being readily observable.

We compare the monopole, quadrupole, and hexadecapole in redshift space calculated from the simulation data with two theoretical model outcomes calculated by Victor. The only difference in the two models is the velocity profile input, which in one case is the radial linear velocity profile and in the other case, the velocity profile calculated from the simulation data. For all the multipoles, we see that the model fits the data with similar accuracy for all the different simulations. As expected, there is a preference towards the model that has the velocity profile from the data given as input, both for the voxel and spherical void definitions. As it stands, we need more accurate modelling and high-resolution surveys if we hope to use the void-galaxy CCF to put constraints on f(R)-modified gravity or massive neutrinos. However, our simulations show that modified gravity and massive neutrinos affect the monopole and quadrupole in opposite ways, making the void-galaxy CCF yet another degenerate observable.

When calculating the model prediction of the void-halo CCF in redshift space, we have in the above case used the real space CCF information from the simulation boxes as model input. For data obtained from a galaxy survey, this information is unknown and can be approximated through a reconstruction process where the redshift space positions are put back into real space by solving for the displacement field (Nadathur et al. 2019a). When performing reconstruction for all simulations, we find that the process performs similarly for all simulated cosmologies, although slightly better for the most massive neutrino case and slightly worse for the modified gravity case, compared to ΛCDM. This could be due to the reduced clustering and lower velocities in the massive neutrino simulations better adhering to the linear approximations of the reconstruction model, and the opposite for the modified gravity simulation. When fitting the CCF model to the data when the reconstruction step is included, the model provides a slightly worse fit to the data. This is expected as we are now providing the model with information that only approximates the actual real space CCF that is needed as model input. The model-data fit is, however, comparable for all the ANUBISIS simulations also in this scenario.

During our analysis, we kept the growth rate, f, constant in all cases, effectively ignoring the scale dependence for modified gravity and massive neutrino cosmologies. We used the expected ΛCDM value, $f={\mathrm{\Omega }}_m^{0.55}$, for all simulations to keep the study consistent, and instead investigated if this choice gives visible biases in our model fits. As recapped above, all model-data fits perform similarly for the different simulations. The linear velocity profile does, however, show a slight improvement for the f(R) simulations, especially visible for the pure ISIS run. We, therefore, tried varying the value of f in this model to see if this could further improve the fit. For the voxel void definition, which already overestimates the clustering, increasing f to values more in line with expectations for f(R)-modified gravity, only contributes to further mismatch between the model and the simulation data. For the spherical void definition, which suffers less from the effect of a sparse tracer sample, increasing f can give a better fit. This only shows that the choice of void finding algorithm and the current sparse tracer sample issue in the void velocity modelling is important and must be consistently handled before possible effects of a different growth rate in the various simulated cosmologies can be resolved.

For all the ANUBISIS simulations, we performed MCMC fits for fσ₈ and the Alcock-Paczyǹski parameter, ϵ. We did this both with and without reconstruction. In both cases, we find that the voxel void definition better recovers the fiducial values of the fσ₈-parameter for the various simulations. The ϵ-values are more consistently recovered between the two void definitions. The discrepancies are most prominent when reconstruction is included, where the fσ₈-values recovered for the spherical voids are consequently underestimated for all simulated cosmologies, as opposed to slightly overestimated without the reconstruction step. Still, for both void definitions, the ability to recover the fiducial parameter values is similar for the different simulations. There is no indication that the parameter estimations are more accurate for a ΛCDM cosmology, compared to f(R)-modified gravity or massive neutrinos.

Our investigations suggest that the current limitations of the void theories must be dealt with in order for the models to be accurate enough to clearly showcase differences between the various simulated cosmologies. This is particularly clear for the velocity profile modelling and the reconstruction method. In addition, higher resolution simulations are necessary to lower the uncertainties in the data. Once this is in place, cosmic voids are a promising ground for studying cosmologies with massive neutrinos and f(R)-modified gravity.

Acknowledgments

We thank the Research Council of Norway for their support. Parts of our computations were performed on resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway. We also wish to thank Slađana Radinović for her contributions in the form of valuable input, discussions and instructions regarding the void component of this paper.

References

Appendix A: ANUBISIS resolution

In the upper panel of Fig. 3, we present the matter power spectrum ratio for the non-ΛCDM simulations, compared to the reference ΛCDM run. We mention that at larger scales there is originally an excess in the power spectrum for the massive neutrino and f(R) + massive neutrino runs. This is attributed to the coarse base grid settings used in the simulations. With the resources available to us at the time of running ANUBISIS, we made appropriate adjustments in order to constrain our memory and time usage. To accomplish this, we chose a base grid of 512³ for all the simulations, which results in a lower resolution at early times and at large scales where there is less clustering. This is due to the nature of the RAMSES AMR-scheme, which creates more refinements as structure grows. In a ratio comparison, such an internal consequence of the algorithm should cancel out, but as the massive neutrinos contribute when calculating the time steps, the simulations containing these particles automatically operate with smaller time steps. This leads to a more detailed time development of the CDM particles in the massive neutrino and f(R) + massive neutrino simulations, compared to the ΛCDM and pure f(R) runs. Effectively, we resolve the simulations containing massive neutrinos better at large scales, resulting in the refinement of more clustering in these cases. This evens out as time passes and more refinements are made, but for large neutrino masses, this also happens for small scales, where the suppressed structure formation leads to fewer refinements within the code, compared to the ΛCDM case. To help even out the issues at large scales, we ran the ΛCDM and pure f(R) simulations with one-fourth of the time step originally calculated by the code. The effect of adjusting the time step is illustrated in Fig. A.1. Ideally, the simulations should be run with a coarser base grid and with a higher particle density, which would also contribute to better resolution overall. This is not done in our case, due to computational costs.

Fig. A.1. Matter power spectrum ratio between a 0.15 eV massive neutrinos simulation and a ΛCDM simulation with varying timestep, performed by ANUBIS. The simulation box has Lbox = 1024 h−1Mpc, Ncdm = Nν = 10243, and a fixed coarse base grid of 5123. The time step is calculated by the RAMSES algorithm for both runs, but for the ΛCDM case it is multiplied by a reduction factor, A. We see here that one way to compensate for the coarse base grid, especially at large scales, is to shorten the time step. This requires a longer run time for the ΛCDM simulation, but not more memory as would be the result of adopting a finer grid.

In the end, this issue is practically equivalent to having a slightly different σ₈-value (2 − 3% lower) at large scales. We have not accounted for this in our analysis, as the resulting changes are small. For our MCMC fits, the true σ₈-value is still encompassed in the flat prior interval, and the only difference in Figs. 18 and 19 would be a slight adjustment of the fiducial fσ₈-line towards lower values for ΛCDM and f(R). Simulations with the same issue still showed good agreement with other codes in the MNCCP (Adamek et al. 2022).

Appendix B: Quadrupole shape

If the redshift space void-halo CCF model (Eq. 3) did not include a velocity dispersion for the mean outflow velocity profile, the shape of the quadrupole in redshift space would be qualitatively similar for all void definitions. This would mean a trough inside the void due to outflow towards the overdensity at the edge, and a peak outside due to infall towards the same overdensity, akin to the shape we see for the voxel voids in Fig. 11. Adding a velocity dispersion to the model effectively corresponds to ‘smoothing’ the RSD quadrupole.

If we take the real space radial coordinate, r, to be equal to the redshift space coordinate, s, in the theory model, then from Eq. (5) it can be shown that for v_r/(raH)≪1 we roughly have Λ. This means that it is deviations from linearity of the void velocity profile (equivalent to deviations from a flat density profile) that lead to a quadrupole inside the void. To better understand the shape of the quadrupole for different void definitions, it is useful to consider a simple toy model. We define a density profile that is perfectly flat inside the void, and equal to zero outside: δ = δ₀ for r < R_void and δ = 0 for r > R_void. This then gives a velocity profile v_r ∝ r inside the void, and $v_r\propto R_{\text{void}}^3/r^2$ outside, if we apply Eq. (8). This again gives ${\xi }_2^s=0$ inside the void, and ${\xi }_2^s\propto R_{\text{void}}^3/r^4$ outside. At the void edge, we have a sharp feature. For real voids, the density profile is never perfectly flat on the inside, and the deviations break the linearity of the outflow velocity, resulting in a trough in the quadrupole inside the void. The flatter the density profile is inside the void, the sharper the peak of the velocity profile is, leading to sharper features in the quadrupole. The monopole, quadrupole, and hexadecapole, in the absence of a velocity dispersion, are illustrated by the green line in Fig. B.1. The sharp features at R_void are visible both for the quadrupole and hexadecapole.

Fig. B.1. Monopole, quadrupole, and hexadecapole for the toy model of the void density profile both with (blue dashed line) and without (green line) velocity dispersion. The flat density profile inside the void results in a sharp feature in the quadrupole at the void radius, which smoothes out and gives a peak right before the trough inside the void.

The above scenario is the case when there is no velocity dispersion. If we add this to the model, the imposed smoothing effect simply reduces the amplitude of the no-dispersion quadrupole for the voxel void definition. This means that the smoothing length given by the size of the velocity dispersion is smaller than the features seen in the no-dispersion quadrupole. On the other hand, for spherical voids, where the density profile is very flat in the void centre and the quadrupole has sharp features as a consequence, this smoothing results in a sharp peak followed by a sharp trough close to the void edge (as seen in Fig. 11). This shows up at the point where the density profile stops being flat, which is also the point where the velocity profile has a sharp peak. The blue dashed line in Fig. B.1 shows this effect clearly for the quadrupole. The smaller the average void is, the larger this first peak is, as the smoothing effect is more dramatic. This shape is also similar to what previous studies have found for Zobov voids (Nadathur et al. 2019a,b; Nadathur & Percival 2019; Massara et al. 2022; Woodfinden et al. 2022), which is not surprising, as these voids also have a fairly flat density profile, δ ∼ −1.0, at the void centre.

All Tables

All Figures

Fig. 1. Comparison of the neutrino suppression on the cb power spectrum and HMF for ANUBIS and GADGET-3 at z = 0 for various values of the sum of the neutrino masses. Top: the neutrino suppression on the cb power spectrum. The greyed-out area indicates the parts of the ratios lying outside of the Nyquist frequency of the simulations. Bottom: the neutrino suppression on the HMF. Only halos with Npart ≥ 11 are included.

In the text

Fig. 2. CDM + baryon (cb) matter power spectrum and HMF of the fofr_small simulation, compared to lcdm_small. The power spectrum and HMF ratios are displayed in the top and bottom panels respectively. For the power spectrum ratio, the different σ8-values of the simulations are taken into account through scaling, and the greyed-out area indicates the parts of the ratios lying outside of the Nyquist frequency of the simulations.

In the text

	Fig. 3. CDM + baryon (cb) matter power spectrum and HMF of the different ANUBISIS simulations, compared to ΛCDM. The power spectrum and HMF ratios are displayed in the top and bottom panels respectively. For the power spectrum, the greyed-out area indicates the parts of the ratios lying outside of the Nyquist frequency of the simulations.
In the text

	Fig. 4. Scale-dependent halo bias for the various ANUBISIS simulations. The upper panel shows the value of the bias for each simulation, while the lower panel shows the ratio with ΛCDM for each model. The greyed-out area indicates the parts lying outside of the Nyquist frequency of the simulations.
In the text

	Fig. 5. Scale-dependent halo bias for the fofr_small and lcdm_small simulations. The upper panel shows the value of the bias for each simulation, while the lower panel shows the ratio with ΛCDM. The greyed-out area indicates the parts lying outside of the Nyquist frequency of the simulations.
In the text

	Fig. 6. Histogram of effective void size for the various ANUBISIS simulations. The voxel voids display a notable difference in the number of voids in the ∼10 − 40 h−1 Mpc range. Both void definitions also show a slight opposite abundance difference amongst the largest voids in the simulation catalogues.
In the text

	Fig. 7. Radial void velocity profiles for all ANUBISIS simulations together with the linear velocity profile and a fit to the more general velocity profile, presented in Sect. 2.1.2. The voxel and spherical void definitions are shown in the left and right columns respectively. The Av-parameter in each panel shows the best-fit value of the general profile.
In the text

	Fig. 8. Radial void velocity profiles for the fofr_small and lcdm_small simulations together with the linear velocity profile and a fit to the more general velocity profile. The voxel and spherical void definitions are shown in the left and right columns respectively. The Av-parameter shows the best-fit value for the general profile.
In the text

	Fig. 9. Linear velocity profile, compared to the fofr_small simulation, for various values of the growth rate, f. The regular linestyle and colour show the velocity profile for the simulation. For ΛCDM we have $f={\mathrm{\Omega }}_m^{0.55}$ ≈ 0.53.
In the text

	Fig. 10. LOS velocity dispersion scaled by σ* as given in Table 3. The voxel and spherical void definitions are shown in the upper and lower panels respectively. All the simulated cosmologies result in similar profile shapes. The largest differences are located towards the void centres, which are dominated by shot noise.
In the text

	Fig. 11. Void-halo CCF for the ANUBISIS simulation data. The voxel and spherical void definitions are shown in the left and right columns respectively. The error bars are estimated using a Jackknife technique and are only shown for the ΛCDM case for a tidier visual representation of the data.
In the text

Fig. 12. Difference between the redshift and real space monopole of the void-halo CCF from the ANUBISIS simulations vs. theory calculated by Victor. The results are obtained by using the voxel and spherical void definitions in the left and right columns respectively. Two theory results from Victor are included, one where the velocity profile is given by linear theory and one where it is given by a template. The template used in this case is the velocity profile calculated from the simulation data.

In the text

	Fig. 13. Quadrupole of the void-halo CCF from the simulations vs. theory calculated by Victor. Layout and input as explained for Fig. 12.
In the text

	Fig. 14. Hexadecapole of the halo-void CCF from the simulations vs. theory calculated by Victor. Layout and input as explained for Fig. 12.
In the text

	Fig. 15. Difference in the LOS coordinate between real space and reconstructed real space for the ANUBISIS simulations shown in the regular linestyle and colour. The weaker full line of the same colour shows the difference between the LOS coordinate in real and redshift space.
In the text

	Fig. 16. Difference in the LOS coordinate between real and reconstructed real space and real and redshift space for the fofr_small and lcdm_small simulations.
In the text

	Fig. 17. Quadrupole of the reconstructed void-halo CCF from the simulations vs. theory calculated by Victor. The real space CCF used as theory input is calculated from the reconstructed real space halo catalogue and the voids identified in it. The layout follows Fig. 12.
In the text

	Fig. 18. MCMC fits for the ANUBISIS simulations. Each triangle plot shows the fσ8 and ϵ estimates of one simulation for both the voxel and spherical void definitions. The grey dashed lines show the fiducial values and the contours display the 68% and 95% confidence intervals.
In the text

	Fig. 19. MCMC fits for the ANUBISIS simulations with reconstruction. Each triangle plot shows the fσ8 and ϵ estimates of one simulation for both the voxel and spherical void definitions. The grey dashed lines show the fiducial values and the contours display the 68% and 95% confidence intervals.
In the text

Fig. A.1. Matter power spectrum ratio between a 0.15 eV massive neutrinos simulation and a ΛCDM simulation with varying timestep, performed by ANUBIS. The simulation box has Lbox = 1024 h−1Mpc, Ncdm = Nν = 10243, and a fixed coarse base grid of 5123. The time step is calculated by the RAMSES algorithm for both runs, but for the ΛCDM case it is multiplied by a reduction factor, A. We see here that one way to compensate for the coarse base grid, especially at large scales, is to shorten the time step. This requires a longer run time for the ΛCDM simulation, but not more memory as would be the result of adopting a finer grid.

In the text

	Fig. B.1. Monopole, quadrupole, and hexadecapole for the toy model of the void density profile both with (blue dashed line) and without (green line) velocity dispersion. The flat density profile inside the void results in a sharp feature in the quadrupole at the void radius, which smoothes out and gives a peak right before the trough inside the void.
In the text