Ringing the universe with cosmic emptiness: void properties through a combined analysis of stacked weak gravitational and Doppler lensing

An essential aspect of cosmic voids is that these underdense regions provide complementary information about the properties of our Universe. Unlike dense regions, voids are avoided by matter and are less contaminated by baryonic processes. The first step to understanding the properties of cosmic voids is to correctly infer their mass profiles. In the literature, various techniques have been implemented. In this paper, we review them and implement a new technique that is based on Doppler lensing. We use a relativistic $N$-body code, \textsc{Gevolution}, to generate cosmological mass perturbations and implement a three-dimensional ray-tracing technique, which follows the evolution of a ray-bundles. We focus on the various properties of cosmic voids (e.g. void size function, 2-point correlation function, and the density profile of voids), and compare the results with their universal trends. We show that when weak-lensing is combined with the Doppler lensing we obtain even tighter constraints than weak-lensing alone. We also obtain better agreement between density profiles within central parts of voids inferred from lensing and density profiles inferred from halo tracers. The implication of the result relevant to the ongoing and prospective low-redshift spectroscopic surveys is briefly discussed.


INTRODUCTION
Cosmic voids -the underdense regions -form a significant portion of the cosmic web filling more than half of the volume of the observable Universe (Vogeley et al. 1994). Because of their peculiar nature (Bertschinger 1985;White et al. 1987;Peebles 2001), significant attention has been paid to understanding the physics of cosmic voids since their first discovery in the late 1970s (Gregory & Thompson 1978;Jõeveer et al. 1978). The physics of cosmic voids introduces a new cosmic probe in understanding the evolution of our Universe Moresco et al. 2022). Recent advancements in both galaxy redshift surveys (Dawson et al. 2013;The Dark Energy Survey Collaboration 2005;Eisenstein et al. 2011;De Jong et al. 2013) and cosmological simulations (Springel 2005;Dolag et al. 2016;Potter et al. 2017;Villaescusa-Navarro et al. 2020) have attracted scientists to focus on voids and their relevance in the evolution of large-scale structures (LSS) to gain insights into our Universe (Patiri et al. 2006;Hahn et al. 2007;Chan et al. 2014;Habouzit et al. 2020;Bayer et al. 2021).
The void abundance, also known as the void size function, is ★ E-mail: mhos5289@uni.sydney.edu.au one of the features of cosmic voids, that is particularly sensitive to the dark energy equation of state (Lee & Park 2009;Bos et al. 2012;Pisani et al. 2015) and offers a new framework for exploring dark energy and modified gravity (Lee & Park 2009;Bos et al. 2012;Spolyar et al. 2013;Cai et al. 2015;Pisani et al. 2015;Zivick et al. 2015;Hamaus et al. 2016;Falck et al. 2018;Perico et al. 2019;Verza et al. 2019). Sahlén et al. (2016) demonstrated that the combined analysis of the abundance of voids and clusters might break the degeneracy among cosmological parameters and provide tighter constraints on modified gravity using extreme value statistics. In recent decades, researchers have also concentrated their efforts on extracting cosmological information by examining various aspects of cosmic voids, including the void-galaxy auto/cross-correlation function (Paz et al. 2013;Hamaus et al. 2014aHamaus et al. , 2015Liang et al. 2016;Cai et al. 2016;Pollina et al. 2017;Achitouv et al. 2017;Nan & Yamamoto 2018;Nadathur & Percival 2019;Nadathur et al. 2019a;Panchal et al. 2020;Khoraminezhad et al. 2022), and the universal density profile (Lavaux & Wandelt 2012;Sutter et al. 2012a,b;Krause et al. 2013;Pisani et al. 2014;Ricciardelli et al. 2014;Hamaus et al. 2014b;Leclercq et al. 2015;; Barreira et al. 2015;Pollina et al. 2017;Panchal et al. 2020;Rodríguez Medrano et al. 2021;Wilson & Bean 2021), which could serve as a potential probe for studying cosmology and galaxy formation with greater precision. Weak gravitational lensing is defined as the distortion of the size (in the form of convergence, ) and shape (in the form of shear, ) of background galaxies by the mass distribution along the line of sight. Thus, it can be directly used to map the distribution of dark matter (Kaiser 1992;Wittman et al. 2000; Bartelmann & Schneider 2001;Hoekstra & Jain 2008). Weak-lensing by cosmic voids, also known as void lensing, has received a lot of attention in recent decades (Amendola et al. 1999;Krause et al. 2013;Higuchi et al. 2013;Cai et al. 2015; Barreira et al. 2015;Baker et al. 2018;Davies et al. 2018Davies et al. , 2021. The weak-lensing tangential shear signal around cosmic voids has recently been investigated (Melchior et al. 2014;Clampitt & Jain 2015;Sánchez et al. 2017) and reported to be moderately significant (4.4 − 7 ) while Fang et al. (2019) found the most significant detection (14 ). Subsequently, Gruen et al. (2016) and Brouwer et al. (2018) utilised a new approach to quantify the tangential shear signal with higher significance (10 − 15 ). Besides the weak-lensing signal by cosmic voids, the lensing imprint of cosmic voids on cosmic microwave background (Chantavat et al. 2016;Raghunathan et al. 2020;Vielzeuf et al. 2021;Kovács et al. 2022) and intrinsic alignments of galaxies around cosmic voids (D'Assignies D. et al. 2022) have recently been studied.
Doppler lensing is an alternative approach to gravitational lensing that originates from redshift-space distortion (arises due to the peculiar velocities of galaxies; Kaiser 1987;Paz et al. 2013;Cai et al. 2016;Hamaus et al. 2017Hamaus et al. , 2020Hawken et al. 2017;Achitouv 2019;Correa et al. 2022) and is represented in terms of weak-lensing (Bolejko et al. 2013;Bacon et al. 2014). As redshift distortions are caused by the underlying matter distribution so, it could be beneficial to examine Doppler lensing in conjunction with weak-lensing to improve the precision of lensing measurements. Bacon et al. (2014) measured the Doppler lensing signal in combination with the weak-lensing signal by using a Newtonian simulation. They compute the 2-point auto/cross-correlation function and observed tighter constraints for cosmological parameters. They also showed that the Doppler lensing dominates over gravitational lensing at medium-to-low redshifts, leading to a new probe to consider in conjunction with weak-lensing for greater precision to study the underlying matter distribution in and around cosmic structures. Subsequently, Bonvin et al. (2017) investigated the effect of Doppler magnification on the size of galaxies and constructed an estimator to measure the velocities using this effect. In our recent paper , we studied both weak-lensing and Doppler lensing signals in and around the cosmic voids/haloes and demonstrated that the most optimal technique for mapping the mass distribution that combines both gravitational and Doppler lensing effects should target the redshift range ≈ 0.3 − 0.4.
In this paper, we investigate further benefits of combining the gravitational and Doppler lensing methods. This is done by the means of numerical simulations. Numerical simulations are essential for developing a deeper understanding of the LSS in the universe as it becomes complex on the nonlinear scale, and scientists have published a large number of works on detecting weak-lensing signals by the mass distribution of LSS using Newtonian simulations (Jain et al. 2000;Takahashi et al. 2011;Valageas et al. 2012). But in the presence of relativistic sources, Newtonian simulation has some limitations, and one may need to consider relativistic simulations to capture all the relativistic effects (Green & Wald 2012). In this work, we will consider a relativistic approach to cosmological simulations developed by Adamek et al. (2016a,b).
In this paper, we measure the lensing signals around cosmic voids and explore various properties of cosmic voids (e.g. void size function, 2-point auto/cross-correlation function, and density profile of voids). We constrain the lensing model parameters from the theoretical predictions and numerical data to see how much constraining power we can enhance by considering solely weaklensing shear or magnification signal and combined lensing (weaklensing + Doppler lensing) signals. In the following, we outline the theoretical background and numerical recipes in Section 2, describe the results of the various properties of cosmic voids in Section 3. The implications of our findings of the lensing measurements are then discussed in Section 4 and the drawn conclusions are summarised in Section 5.

Weak gravitational lensing
Gravitational lensing is a unique technique for directly mapping the mass distribution of dark matter throughout the cosmos. Weak gravitational lensing systematically distorts the images of source galaxies by mapping them to new locations on the sky. The lens equation for a gravitationally lensed image can be expressed as where is the deflection angle, is the angular position of the source on the sky, and is the observed position of the lensed image. Following Bartelmann & Schneider (2001), the deformation matrix A is the equivalent Jacobian matrix of the lens mapping and can be written as This deformation matrix can be parameterised using more physically instructive notions such as convergence ( ) and shear ( ≡ 1 + 2 ), and it can be represented as And, the total lensing magnification is expressed by the determinant of the inverse matrix The convergence describes the magnification and de-magnification of source images, whereas the shear describes how much the image is stretched along its axis. Depending on the mass distribution of the universe, the convergence along the line of sight can be positive, implying a magnified image of the source, or negative, implying a de-magnified image of the source. The positive or negative value of the convergence indicates that the photon travels through an overdense or underdense area of the universe. The weak-lensing statistical quantities i.e. convergence, shear, and magnification are also related to the projected surface mass density (Σ) and differential surface mass density of the lens (ΔΣ) and can be expressed as where is the physical transverse distance on the lens plane and Σ cr is written as where are the angular diameter distances between the observer (o), lens (l), and source (s).

Doppler lensing
Doppler lensing is the apparent change in object size and magnitude due to peculiar velocities. There are two fundamental differences between standard weak gravitational lensing and Doppler lensing. Firstly, weak-lensing is the phenomenon in which the mass distribution along the line of sight influences both the size and shape of the source galaxies, whereas Doppler lensing affects only the size of the galaxies while leaving their shape unaltered. Secondly, the largest contribution to the weak-lensing signal comes from the matter that is approximately located in mid-distance between the observer and the observed galaxies, and consequently, the precise location of the source galaxies is less significant. In contrast, Doppler lensing is the local consequence of the observed structure, thus the source is at the lens. The reason why for the Doppler lensing the precise location of the source is important is that the effect is caused by the peculiar motion of galaxies (Bonvin 2008;Bolejko et al. 2013;Bacon et al. 2014;Bonvin et al. 2017). The convergence sourced by the Doppler lensing is given by following normalised equation 1 where is the Doppler convergence, is the source redshift, is the co-moving distance of the sources, is the Hubble parameter, is the velocity between source galaxies and the observer, is the speed of light, and is the unit vector from the source to the observer. There are two scenarios can be observed because of the propagation direction of the photon • · > 0, implying that the objects are moving towards us, thus the Doppler convergence < 0. This condition suggests that the observed structures will be measured smaller and dimmer than the usual objects at their observed redshift.
• · < 0, implying that the objects are moving away from us, thus the Doppler convergence > 0. This condition suggests that the observed structures will be measured larger and brighter than the usual objects at their observed redshift.
Due to the factor ( ) −1 in Eqn. 8, the Doppler lensing convergence decreases with redshift, which limits its application to low-redshift only.

−body simulation
In this paper, we use relativistic -body code G 2 (Adamek et al. 2016a,b) to generate the weak perturbations, in which the perturbed Friedmann-Lemaître-Robertson-Walker line element 1 we use the Hubble parameter with a unit of ℎ −1 Mpc, so the product of and is a dimensionless quantity. 2 https://github.com/gevolution-code/gevolution-1.2 of the metric has the form where represents the scale factor of the background, are the comoving Cartesian coordinates, and is the conformal time. Φ and Ψ are the scalar perturbations, is the vector perturbation and ℎ is the tensor perturbation. To minimise computational time and memory, we consider only scalar (Φ and Ψ) perturbations here since vector and tensor perturbations have the negligible effect on light propagation. For more details we refer the reader to these papers (Lepori et al. 2020;Ema et al. 2022).
The G simulation consists in a baseline of standard Λ cold dark matter (ΛCDM) model and we run it by using the following cosmological parameters: ℎ = 0.67556, Ω c = 0.2638, Ω b = 0.048275, and a radiation density that includes massless neutrinos with eff = 3.046. Our considered CDM simulation follows the evolution of 256 3 CDM particles in a cubic comoving volume (320 ℎ −1 Mpc) 3 from redshift ini = 127 to the present epoch, and the linear initial conditions are computed with CLASS (Blas et al. 2011) by assuming a primordial power spectrum with amplitude = 2.215 × 10 −9 and spectral index = 0.9619. This provides a spatial resolution of 1.25 ℎ −1 Mpc, and corresponds to a minimum resolved halo mass (with 20 particles) of 5 × 10 11 ℎ −1 M . We run G simulation for 21 realisations (by changing the initial conditions but keeping the identical cosmological parameters), which takes around 250 CPU hours to provide the output in the forms of particle snapshots and weak potentials.

Void finder
Cosmic voids are the under-dense patches of the universe surrounded by filaments and sheets of galaxies that make up the cosmic web. There are a variety of void-finding techniques in the literature, including Voronoi tessellation and watershed methods (Platen et al. 2007;Neyrinck 2008;Sutter et al. 2012a;Nadathur et al. 2019a), Delaunay density estimation methods (Schaap & van de Weygaert 2000;Zhao et al. 2016), 2D projections (Clampitt & Jain 2015), Hessian-based methods (Adermann et al. 2017(Adermann et al. , 2018, tunnel methods Davies et al. 2021), spherical underdensity estimation methods (Colberg et al. 2005;Padilla et al. 2005;Li 2011;Pan et al. 2012;Hoyle et al. 2012;Villaescusa-Navarro et al. 2018). We recommend the readers to these papers (Colberg et al. 2008;Cautun et al. 2018) for further information on the comparison of void finding methodologies. Here we employ R (REalspace VOid Locations from surVEy Reconstruction) 3 , a publicly available code, to find the cosmic voids. It is a three-dimensional (3D) void finder that is been utilised extensively in simulations as well as observed catalogues (Nadathur et al. 2019a,b;Raghunathan et al. 2020;Li et al. 2020;Jeffrey et al. 2021;Khoraminezhad et al. 2022). The methodology of the R is based on the parameterfree void-finding algorithm Z (ZOnes Bordering On Voidness; Neyrinck 2008; Sutter et al. 2015), which uses Voronoi tessellation to estimate local density minima from the input tracers.
The code R assumes the centre of the void is the centre of the largest empty sphere that can be inscribed in the void, and consequently the place of the lowest density of the galaxy distribution . The barycentre from watershed voids has been extensively used in literature because it is more stable against shot noise and is dependent on the entire geometry of the cosmic void (Hamaus et al. , 2016Cousinou et al. 2019;Hamaus et al. 2020). However, the main results (lensing results) of this work, which will be discussed in Section 4, are not affected by the various assumptions about the centre of the cosmic voids. The ray-tracing algorithm is not sensitive to such a definition (it is affected by the matter distribution along the line of sight, not by any artificial definition of the centre of the void). Though realistic void shapes are far from spherical (for a nice visualisation we refer the reader to Fig. 3 of Nadathur 2016) but for simplicity the Z algorithm computes the effective void radius void from the total volume of each void (which is effectively the sum of the volumes of its constituent Voronoi cells, void = ) as void = The Z algorithm also computes the minimum galaxy density contrast near cosmic voids using the equation ,min = ,min /¯−1, where ,min is the minimum galaxy number density and¯is the average matter density. It is also called as core density contrast of each void from the estimation of Voronoi tessellation field estimator (VTFE) reconstruction.
The R algorithm takes a discrete tracer distribution as an input, in this instance DM haloes or particles. Haloes are identified by using the publicly available code R (Behroozi et al. 2013), a phase-space friends-of-friends (FOF) halo finder. We feed the DM particle snapshots from 21 distinct G simulations at = 0.22 into R 4 , which provides on average 3 × 10 4 DM haloes per realisation. The redshift of = 0.22 was chosen on the basis of the relevance of the Doppler lensing, which dominates for < 0.3 (see Bacon et al. (2014)). A detailed explanation also given in the conclusion part of our recent paper ). To find DM haloes, we set the FOF refinement fraction = 0.7 and 3D-FOF linking length = 0.28, for more details we refer the reader to the original paper (Behroozi et al. 2013). Then we utilise DM haloes as an input tracer into R to identify voids, which provide on average 210 cosmic voids per realisation with a radius spanning from 10 ℎ −1 Mpc − 65 ℎ −1 Mpc.
Running the VTFE technique on a large number of DM particles is computationally challenging, so we utilise the downsampling routine to randomly sub-sample the DM particles from the output of 21 distinct G simulations to minimise the computational time and memory requirements, which is common in the literature Khoraminezhad et al. 2022). In the down-sampling routine for each realisation of G simulation, we keep the DM particles down to a constant average density of 1.68 × 10 6 particles per cubic box size (320 ℎ −1 Mpc) 3 , which equates to 10% of the DM particles. Then we use the sub-sample of the DM particles as an input tracer into R to identify voids, which provide on average 9315 cosmic voids per realisation with a radius spanning from 2 ℎ −1 Mpc − 22 ℎ −1 Mpc. It is worth noting that the void sizes and counts vary depending on the average density of the various input tracers. 4 https://bitbucket.org/gfcstanford/rockstar/src/main/
In this paper, we use 3D Ray Bundle Tracers (3D-RBT; Ema et al. 2022;Hossen et al. 2022), the 3D ray-tracing algorithm that relies on the design of the RBM (Fluke et al. 1999(Fluke et al. , 2002. The steps of the ray-tracing algorithm are briefly described below, however for more details on 3D-RBT, we refer the reader to our recent works Hossen et al. 2022): • Generating weak perturbations: the first step is to generate relativistic weak potentials using -body simulation. We run 21 G simulations with the same cosmological parameters but different initial conditions and save the weak potentials and particle snapshots for different redshifts.
• Finding cosmic structures: from the particle snapshots of G simulations at redshift = 0.22, we identify DM haloes using the publicly available code R . Then we utilise DM haloes as an input tracer into R to identify cosmic voids.
• Projection directions: we fix the observer at = 0 and project bundles of photons (initially circular shape) in the directions of cosmic voids with three distinct groups of voids radii by slightly changing the projection angles.
• Solving null geodesics: at each step of the integration we obtain the metric (9) from G . This is then used to evaluate the Christoffel symbols and obtain the direction of the path of the photon, and consequently the full bundle. The RBM implemented in this paper uses eight null geodesics and a central null geodesic per bundle. In addition to this at each step of integration we check the null condition. This is an essential sanity check to ensure that we are obtaining the correct photon path profile.
• Ellipse fitting and lensing statistical quantities: the initial circular shape of the bundle deforms as photons traverse through the cosmic web. The shape of the bundle is calculated by fitting an ellipse into the bundle's shapes. This is then used to compute the weak-lensing magnification and shear as where represents the area of the image, is the area of the source, a & d are the semi-major and semi-minor axes of an ellipse, respectively. The weak-lensing convergence ( ) is then calculated from Eqn. 4, and is It is worth noting that ray-tracing algorithms are computationally challenging, requiring significant amounts of time and memory. For technical details and challenges of implementation of these algorithms we refer the reader to our earlier works Hossen et al. 2022).

Doppler lensing algorithm
Doppler lensing is a tool that can be used to measure the lensing signal (in conjunction with the weak-lensing signal) at low redshifts (Bolejko et al. 2013;Bacon et al. 2014;Bonvin et al. 2017;Hossen et al. 2022). Unlike weak gravitational lensing, which is caused by the mass distribution along the line of sight, Doppler lensing is purely related to the velocity of the source galaxies. Doppler lensing algorithm is based on solving Eqn. 8 and the steps of computing Doppler convergence are as follows: • Generating particle snapshots: to generate particle snapshots at redshift = 0.22, we conduct 21 G simulations with different initial conditions but keeping the identical cosmological parameters. These are the same simulations and snapshots used in the weak-lensing algorithm part.
• Finding cosmic structures: by using R , we detect DM haloes and extract their information, such as positions, masses, velocities, and so on. Then, employing R , we identify cosmic voids and extract their information, such as positions, radii, volume, and so on.
• Computing impact parameter: then we compute the value of the impact parameter, which is the distance from the centre of the void to the centre of the haloes, modified by the factor of cos where is the angle between the direction pointing from the observer towards the centre of the void, and the direction pointing from the centre of the voids towards the galaxy, cf. Hossen et al. (2022).
• Computing and stacking the Doppler convergence: then we calculate the dot product of the normalised velocities ( / ) of the haloes around cosmic voids and the direction vectors ( ) between the observer and the observed haloes. The Doppler convergence value is then obtained by multiplying the remaining parts of Eqn. 8 with the dot product result. Finally, we stack Doppler convergence values for voids with three distinct subgroups of void radii. Sheth & van de Weygaert (2004) have modelled for the first time the statistical analysis of the Void Size Function (VSF), also known as the void abundance, which is simply characterised as the number of voids in a given radius bin at a particular redshift. It is a potential approach for probing dynamical dark energy and modified gravity (Bos et al. 2012;Cai et al. 2015;Pisani et al. 2015;Verza et al. 2019), constraining neutrino masses (Massara et al. 2015;Kreisch et al. 2019;Schuster et al. 2019;Zhang et al. 2020;Contarini et al. 2021), standard and coupled dark energy cosmologies (Pollina et al. 2016), cubic Galileon and nonlocal gravity cosmologies (Barreira et al. 2015), cluster-void degeneracy breaking (Sahlén et al. 2016), etc.

Void size function
Here we aim to focus on VSF for two different input tracers that are used to identify cosmic voids. Figure 2 shows the cumulative VSF as a function of the effective void radius (number density of voids with radii above void ) for 21 realisations of G simulations at redshift = 0.22. As expected the VSF decreases with the increasing values of effective void radius in both halo field and DM particle field voids. It is observed from Fig. 2 that the VSF is clearly dependent on the input tracers that are utilised to identify cosmic voids. In comparison to the voids identified in both fields, the DM particle field voids, for example, are smaller and more numerous than the halo field voids. The fundamental reason for such variability in void sizes and counts is that the distribution of haloes is sparser than the distribution of DM particles ).

Void 2-point correlation function
Here we aim to measure the 2-point correlation function (2PCF) for void-void, halo-void, and DM-void configurations. The voidgalaxy cross-correlation function is a useful tool since it encodes the same information as the stacked radial density profile of voids (for more details see Section 3.3). In recent decades, much emphasis has been paid to mapping the stacked profile of cosmic voids using the auto or cross-correlation function in both real and redshift space (Paz et al. 2013;Hamaus et al. 2014aHamaus et al. , 2015Liang et al. 2016;Cai et al. 2016;Pollina et al. 2017;Achitouv et al. 2017  compute the auto/cross-correlation function as where represents the data catalogue, is the synthetic random catalogue, ( ) and ( ) represent the pair counts with separation in the data and random catalogues respectively. It is worth noting that, in the context of uniform periodic randoms such as for periodic simulation boxes, N analytically calculates the random pairings ( ) to minimise computing time.
The variation of the halo-void and void-void 2PCF as a function of separation ( ) for 21 realisations of G simulations at redshift = 0.22 is shown in Fig. 3. These correlation functions are calculated using data from all voids, haloes, and DM particles (down-sampled) for each realisation of the simulation box at redshift = 0.22. Then for certain separation bins, we compute the mean auto/cross-correlation function and demonstrate it in Fig. 3 as a function of the separation. As expected, the characteristics of these correlation functions are quite similar in both halo field and DM particle field voids. However, for halo field voids, higher separation bins are required because the sparser distribution of haloes results in larger void sizes than DM particle field voids. Again, the higher 1 variation of the void-void correlation function for halo field voids is due to the smaller number of voids identified in the halo distribution than DM particle distribution. The profiles for these correlation functions are consistent with the existing literature Cai et al. 2016;Pollina et al. 2017;Achitouv et al. 2017;Nan & Yamamoto 2018;Khoraminezhad et al. 2022).
For the void-void correlation function in both halo and DM particle fields, we observed that the peak of the void-void correlation function is approximately two times the mean effective void radius (i.e. ∼ 2¯v oid ) of the sample, which is consistent with the existing literature (Hamaus et al. 2014c;Kreisch et al. 2021). It is because void walls typically contact at twice the average effective void radius since they can't overlap, a consequence of void exclusion scale (for more details, see Hamaus et al. (2014a); Chan et al. (2014)). Due to this, the sample will likely have two void centres, each located at a distance equal to twice the sample's mean effective void radius (Kreisch et al. 2021). Generally speaking, voids can't be separated by smaller sizes since they will inevitably merge. Additionally, void-void correlation function appears to have negligible systematic errors and biases from redshift-space distortions, making this statistic the cleanest one for identifying geometric distortions (Hamaus et al. 2014c).
For the tracer-void (halo-void or DM-void) correlation function in both halo and DM particle fields, we observed that the peak of the tracer-void correlation function is approximately equal to the mean effective void radius (i.e. ∼¯v oid ) of the sample. The compensation walls of voids, which are composed of high-density structures including sheets, filaments, and clusters of galaxies, are the most dynamic areas. Galaxies far from the centre of the void have a net overdensity and so exhibit consistent infalls rather than outflows on large scales. So, random motions of these structures are responsible to observe the peak of the correlation function within the compensation wall of cosmic voids, a consequence of galaxyvoid exclusion (Padilla et al. 2005;Chan et al. 2014). It is also observed that the halo-void correlation function provides a slightly higher amplitude than the DM-void correlation function in the DM particle field voids (see right side of Fig. 3). The reason for such variation is due to the tracer bias ( t ) around cosmic voids. Since the distribution of haloes is more sparser than the distribution of DM particles (which indicates halo distributions exhibit stronger fluctuations than the DM), one may obtain a higher tracer bias factor ( t > 1) for computing the halo-void correlation function than the DM-void correlation function (for further information on how tracer bias impacts void properties, see these papers (Pollina et al. 2017).

Void density profile
Besides the void sizes, the distribution of the spherically averaged stacked matter density profile in and around cosmic voids is a significant quantity to shape our knowledge about voids (Lavaux & Wandelt 2012;Sutter et al. 2012a,b;Krause et al. 2013;Pisani et al. 2014;Ricciardelli et al. 2014;Hamaus et al. 2014b;Leclercq et al. 2015;Barreira et al. 2015;Pollina et al. 2017;Panchal et al. 2020;Rodríguez Medrano et al. 2021;Wilson & Bean 2021). Hamaus et al. (2014b) introduced the following simple empirical function for capturing the spherically averaged stacked matter density profile in and around cosmic voids where is the central density contrast at which the radial distance from the void centre = 0, t is the average density of tracer, is a scale radius at which the void profile density vt = t , and the inner and outer slopes of the void profile are represented by and .
Cross-correlation (void-tracers), by definition, counts tracers at a given distance from the void centre per unit volume, implying that the void-tracers correlation function carries the same information as the stacked radial density profile of voids Pollina et al. 2017Pollina et al. , 2019. Following Pollina et al. (2017) the average number density of the tracers from the centre of the cosmic voids, also known as the average stacked density profile of voids, can be mathematically expressed as where r is the radial distance from the centre of the void, v is the number of voids, t is the number of tracers, v is the average void density, t is the average tracer density, is the total observed volume, c is the centre of the th void, t is the centre of the th tracer, vt is the void-tracer number density, is the Dirac delta function, and vt represents the void-tracer correlation function.
The identification of cosmic voids in the halo fields offers closer to real observational applications, so we are interested in the halo field voids from now. Another reason is that we will look at lensing measurements (which is related to the integrated stacked density profile of cosmic voids) in the context of ray-tracing algorithms later, and we will need to account for the complete picture of the density distribution of the simulation volume to capture the lensing characteristics adequately.
The halo-void correlation function (equivalent to the average density profile of cosmic voids) is then computed using N for three distinct subgroups of void radii, and the theoretical model data is generated using Eqn. 14. We then employ Markov Chain Monte Carlo (MCMC) method to estimate the posteriors of the model parameters from log-likelihood analysis. To perform this, we utilise a publicly accessible code Emcee (Foreman-Mackey et al. 2013). The following parameter vector summarises the free parameters of our model: and we consider the following uniform prior ranges for the values of our model parameters  (17) The results of the MCMC analysis in terms of the constraints on the parameters of the density profile are presented in Fig. 4. The best fit models are presented in Fig. 5, which shows the variation of the average density profile of cosmic voids as a function of impact parameter ( /¯v oid ) for 21 realisations of G simulations at redshift = 0.22. We find good fits of our simulation data with model. The different colour markers show three distinct groups of void radii with best fit parameter values. It is clear from Fig. 5 in the regime 0 ≤ /¯v oid < 1, the profiles of average stacked density of cosmic voids for three distinct groups of void radii increase from minimum value (∼ −0.8) towards the background value i.e. ∼ 0. The density profiles attain their maximal value in the regime 1 ≤ /¯v oid < 2, the region associated with the filaments and walls around the voids, and then start to drop with increasing the impact parameter. And, in the final regime /¯v oid ≥ 2 the density distribution reaches its average value from the centre of the cosmic voids and the density profile of voids follows a constant background value, ∼ 0. The properties of the stacked density profiles of voids presented here are consistent with existing literature (Ricciardelli et al. 2014;Hamaus et al. 2014b;Leclercq et al. 2015;Barreira et al. 2015;Pollina et al. 2017;Panchal et al. 2020;Rodríguez Medrano et al. 2021;Wilson & Bean 2021).

Model fits
Integrating the stacked density profile as described in Eqn. 14 along the line of sight provides the theoretical prediction for the statistical lensing measurement. The projected surface mass density and differential surface mass density are connected to the lensing statistical quantities (i.e. convergence, shear, and magnification). Following Krause et al. (2013), the differential surface mass density can be expressed as where the projected surface mass density can be written as and los is the line-of-sight distance, void is the distance to the centre of the void, and is the angle between the direction to the centre of the voids and direction of the line of sight; finally¯is the cosmological mean mass density. Therefore, the theoretical predictions for weak-lensing convergence ( theory ) and magnification ( theory ) can be written as where los is the density along the line-of-sight and can be represented as Eqn. 14. And, the weak-lensing tangential shear profile can be expressed as

MCMC analysis and likelihoods
We use a MCMC approach to estimate parameter posteriors and generate likelihood contours for the entire parameter space. We employ consistent priors, as described in Section 3.3 but for lensing measurement we fix and parameters (by taking the best fit parameter values listed in Table 1) and estimate the model parameters using log-likelihood analysis. We compute the likelihood functions for measuring weak-lensing shear and magnification, as well as Doppler convergence, using the corresponding theoretical predictions (as described in Section 2.1.2 and Section 4.1) and numerical data, but for combined lensing analysis, we compute the total log-likelihood as the sum of all log-likelihood functions as where L WL− , L WL− , and L DL− are the weak-lensing magnification, weak-lensing shear, and Doppler lensing convergence likelihoods, respectively. The log-likelihood of the weak-lensing shear and magnification, and Doppler convergence can be expressed as where can be weak-lensing shear ( ), magnification ( ), and Doppler convergence ( ) and data is the error of the numerical data of each bin. We first use the MCMC approach to constrain the parameters of a void using the gravitational lensing shear. Figure 6 shows the posterior distribution of the model parameters for weak-lensing tangential shear measurement by stacking voids having three distinct groups of void radii, and table 2 lists the best fit parameter values for computing stacked weak-lensing tangential shear around cosmic voids. It is clear from Fig. 6 that the bigger voids provide the smallest contours for weak-lensing shear measurements, whereas the smaller voids give the largest contours, which is consistent with the density contrast contours of cosmic voids in Fig. 4. It is likely because the average number density in cosmic voids is higher for smaller voids, causing fluctuations in the number density of void galaxies inside cosmic voids, leading to higher statistical uncertainties and weakening the constraining power as compared to bigger voids. However, it is worth noting that due to the effects of tracer bias around cosmic voids (Pollina et al. 2017, the shear constraints underestimate density and overestimate the parameter as compared to halo tracers. We next investigate the constraints from weak-lensing magnification. The contours of the model parameters for weak-lensing magnification measurement for three distinct groups of void radii is depicted in Fig. 7. The contour plots generated from the weaklensing magnification measurement, like the weak-lensing tangential shear measurement, are the largest for smaller voids as compared to the bigger ones. But the constraints of the lensing model parameters are tighter for the weak-lensing tangential shear measurement than those for the weak-lensing magnification. As in the case of the weak-lensing shear, when compared against the profile obtained based on the halo tracers, the weak-lensing magnification underes- timates density and overestimates . Due to the fluctuation of the density distribution, voids identified in the halo field provide a higher tracer bias factor than the DM particle field (Pollina et al. 2017. This tracer bias around cosmic voids is responsible for the underestimation of density when a halo tracer is used rather than a DM tracer. Figure 8 shows the contours of the model parameters for Doppler convergence measurement for three distinct groups of void radii. We observed that looks good but is not constrained by Doppler lensing.  Finally, we want to combine all the lensing signals (weaklensing shear + weak-lensing magnification + Doppler lensing). Figure 9 shows the posterior distribution of the model parameters for combined lensing measurement by stacking voids having three distinct groups of void radii at lens redshift = 0.22. Table 2 lists the best fit parameter values for combined lensing analysis. For combined lensing analysis, we found a better improvement of the contour size for all ranges of void radius, and the size of the contours is larger for smaller voids as compared to bigger ones. It is apparent from Figs. 9 and 10 the constraints of the lensing model parameters are significantly tighter for the combined lensing case than weak-lensing alone.

CONCLUSION AND OUTLOOK
Weak gravitational lensing and Doppler lensing signals have been examined numerically in this work within the framework of relativistic -body simulations. We began by examining the different features of cosmic voids, such as the VSF, the void 2PCF, and the void density profile, for two distinct sets of input tracers. Numerous studies have identified these void characteristics extensively, but this is the first time we have compared the findings for halo-field and DM particle field voids in the context of relativistic -body simulations. Then, using our ray-tracing technique (3D-RBT), we projected bundles of photons in the direction of cosmic voids with three distinct groups of void radii and extracted the statistical quantities associated with weak-lensing, namely convergence, shear, and magnification. We also constructed a 3D Doppler lensing algorithm to quantify the Doppler convergence around the cosmic voids. We investigated the impact of mass distribution on lensing signals around cosmic voids and ran MCMC samplings from our mock weak-lensing and Doppler lensing data to forecast the accuracy for measuring lensing signals.
For the analysis of cosmic void properties, we considered two tracer fields (such as DM particles and haloes) and then identified cosmic voids from those tracer fields. Due to the sparser halo distributions, the number of cosmic voids in the DM particle field is smaller and more numerous than those in the halo field. To begin, we examined the VSF characteristics as a function of the radius of cosmic voids and observed that they are highly dependent on the tracer used to identify cosmic voids (cf. Fig. 2). The characteristics of halo-void and void-void 2PCF are quite similar in both halo field and DM particle field voids (cf. Fig. 3). The void-void correlation function for halo field voids has a higher statistical uncertainty due to the smaller number of voids identified in the halo distribution than DM particle distribution. Then we analysed the density contrast of cosmic voids as a function of impact parameter and showed that the characteristic profiles of density contrast for three distinct groups of void radii (cf. Fig. 5) follow the universal trends (Krause et al. 2013;Ricciardelli et al. 2014;Hamaus et al. 2014b;Barreira et al. 2015).
The most important result of this paper was the analysis of constraining lensing model parameters from the theoretical predictions and numerical data, followed by the combination of weaklensing and Doppler lensing signals to increase constraining power. The weak-lensing tangential shear profile provides comparatively tighter constraints for the density contrast as compared to the weaklensing magnification for all three distinct groups of void radii. The combined lensing (weak-lensing shear + weak-lensing magnification + Doppler lensing convergence) analysis generated contours are smaller than any of the individual contours (e.g. weak-lensing shear or magnification), for all combinations of parameters. It shows that the combined weak-lensing and Doppler lensing signals improved the constraints, and also brings it to a better agreement with the profile extracted from the halo-tracers. The results of this study should be considered as a proof-of-concept that studied and demonstrated the capability of mapping the matter distribution inside cosmic voids using Doppler lensing in combination with gravitational lensing. The findings provided here are applicable to current and future low-redshift spectroscopic surveys, such as the DESI Bright Galaxy survey (Ruiz-Macias et al. 2020;Zarrouk et al. 2022 Recently, Davies et al. (2021) have shown that combining void abundance with weak-lensing shear substantially enhanced the constraining power and resulted in tighter constraints on the measurement of cosmological parameters. It would be interesting to assess the constraining power by combining void abundance, weak-lensing statistics, and Doppler lensing signal to precisely measure the lensing model parameters and constrain the cosmological parameters. The covariance effects on the joint measurement of the weak-lensing and Doppler lensing signals and taking shape noise into account would be interesting, but we leave this for future studies.