Why Cosmic Voids Matter: Nonlinear Structure&Linear Dynamics

We use the Magneticum suite of state-of-the-art hydrodynamical simulations to identify cosmic voids based on the watershed technique and investigate their most fundamental properties across different resolutions in mass and scale. This encompasses the distributions of void sizes, shapes, and content, as well as their radial density and velocity profiles traced by the distribution of cold dark matter particles and halos. We also study the impact of various tracer properties, such as their sparsity and mass, and the influence of void merging on these summary statistics. Our results reveal that all of the analyzed void properties are physically related to each other and describe universal characteristics that are largely independent of tracer type and resolution. Most notably, we find that the motion of tracers around void centers is perfectly consistent with linear dynamics, both for individual, as well as stacked voids. Despite the large range of scales accessible in our simulations, we are unable to identify the occurrence of nonlinear dynamics even inside voids of only a few Mpc in size. This suggests voids to be among the most pristine probes of cosmology down to scales that are commonly referred to as highly nonlinear in the field of large-scale structure.


Introduction
Observations of the cosmic microwave background (CMB) radiation have revealed the presence of tiny fluctuations in temperature and density when the Universe was only 370 000 years old [1,2]. Since then, gravity amplified these fluctuations by many orders of magnitude, building up what is known as the cosmic web in the process. What starts off as a Gaussian random field that is fully specified by its two-point function, evolves into an intricate system of gravitationally bound structures. These are composed of a multi-scale network of sheets, connected by filaments, which in turn are connected by dense nodes [3]. At the same time the remaining space is evacuated and expands at an increasing pace, creating enormous voids [4][5][6][7][8][9].
Today we can only infer this web of structures indirectly via large surveys that map out its luminous tracers, such as galaxies. The dominant fraction of the cosmic web's matter content is assumed to be composed of some form of cold dark matter (CDM), whose nature is yet to be determined. On top of that, the observed late-time acceleration in the expansion of space has revealed the presence of dark energy, so far consistent with a cosmological constant Λ [10,11]. Despite these big unknowns, simulations have allowed us to investigate structure formation in a ΛCDM universe in great detail [e.g., 12,13].
In addition, the study of radial velocity profiles has become important in the context of modeling redshift-space distortions (RSD) for the observable shapes of voids in redshift space [e.g., [58][59][60][61][62][63][64][65][66][67]. However, tracer velocities can hardly be measured directly, so these models make use of local mass conservation to relate velocities to densities. Although voids represent nonlinear objects with significant density contrasts with respect to the mean background value, the relationship between density and velocity revealed itself to be remarkably well described by the continuity equation at merely linear order in perturbation theory [22]. This insight has been key for the success of RSD models for voids. In contrast, such models are challenged by complex nonlinearities and non-Gaussian statistics appearing in the context of galaxy clustering, which considerably limit the extraction of cosmological information to relatively large scales [68].
The aim of this paper is to study the fundamental properties of voids in more detail with the help of hydrodynamical simulations. We compare different estimators of the density and velocity profiles of stacked voids and explore the use of weighting schemes. Furthermore, we examine how mass cuts and subsamplings of the underlying tracer distributions affect these profiles and test how well linear mass conservation holds up around voids. The outline of the paper is as follows: section 2 provides details on the simulation suite and section 3 introduces the void finder applied to our tracer catalogs, as well as estimators for void profiles. 3 × 10 9 M ⊙ /h in ultra-hr to select halos for void finding. These mass cuts are above the resolution limit of each corresponding box at redshift z = 0.29, as can be seen in figure 1. Halos with masses below these mass cuts will not be considered in any further analysis in this work. In all the aforementioned boxes we analyze voids found in halos from the hydrodynamical simulations and voids found in subsamplings of the underlying CDM that formed these halos.

VIDE void finding
To identify voids in both the halo, as well as the underlying CDM distribution, we make use of the Void IDentification and Examination toolkit vide 2 [97]. vide implements an enhanced version of the ZOnes Bordering On Voidness algorithm zobov [98], which is a watershed algorithm [99] that identifies local basins in the three-dimensional density field estimated from the positions of tracer particles. This density field is constructed via Voronoi tessellation, where each tracer particle j is assigned a unique Voronoi cell of volume V j .
Voronoi cells are defined as the volume of space that is closer to its associated tracer particle than to any other particle in the entire simulation box. Hence, the volumes of all Voronoi cells combined make up the whole simulation box. The density anywhere inside the cell of particle j is then simply given by n j = 1/V j . Starting from the local minima in the density field, the watershed algorithm searches for neighbouring cells with monotonically increasing density to find extended basins of density depressions, our cosmic voids. Additionally, vide allows for a density-based merging threshold within zobov, which is a free parameter. Adjacent basins will be added to a void only if the lowest density along their common ridge line is below this threshold in units of the mean tracer densityn. A low threshold prevents voids from extending deeply into overdense structures, thereby limiting the depth of the void hierarchy [97,98]. In case two adjacent basins are merged, there will be two voids with one encompassing the other, creating a 'parent' and a 'child' void, hence the original number of voids is maintained even after merging. The difference is that the child-level void is now considered as a sub-void of the parent-level void, which encompasses the volume of both voids. In this way, every void can only have one parent, but multiple children (sub-voids). The default value is set to a very low number to prevent any merging of adjacent basins, which results in a sample of parent voids without any sub-voids. A higher value approaching values of order one and above creates a void hierarchy with potentially many levels of sub-voids, but the total number of voids remains independent of the merging threshold. Unless stated otherwise, in this paper we use a default merging threshold of 10 −9n (no merging).
In the literature it is also common to use a threshold value of 0.2n for merging voids, which has a special physical significance. It derives from the spherical expansion model for an inverted top-hat perturbation in an Einstein-de Sitter universe, where the value of 0.2n marks the characteristic density inside the top hat when shell crossing occurs in its boundaries [98,100,101]. However, this only applies when the threshold is defined in the full matter density field and when spherical symmetry is assumed. For voids defined in the number density field of tracers one additionally has to account for the tracer bias [28,[102][103][104]. Moreover, the sparsity of tracer particles effectively smooths the density field on scales below their typical separation, which can affect density ridges between voids [105].
When investigating the effects of merging, we therefore restrict ourselves to the two extreme thresholds for merging, namely a value close to zero for no merging (our default) and a value of infinity for a fully merged void hierarchy. The former case will be referred to as an isolated, the latter as a merged void catalog, respectively. The resulting catalogs encompass all cases with a finite merging threshold in between. The merged catalog consists of one parent void with a deep hierarchy of children, whereas the isolated catalog only contains separate, non-overlapping voids.
Performing the void finding with vide results in catalogs of non-spherical voids with various properties. The center of each void is defined as the volume-weighted barycenter of all of its constituent Voronoi cells at the comoving tracer locations x j : We can think of this barycenter as the geometric center of a void, since its position is constrained by its boundary, which contains most of the void's tracer particles. This definition makes the position of the center very robust against Poisson fluctuations and preserves information about the void topology. Note that the barycenter does not necessarily coincide with the position of the lowest density Voronoi cell due to the lack of spherical symmetry in voids. For the same reason we can only define an effective radius r v , which corresponds to the radius of a sphere of identical volume as the void. It is calculated as the sum over all its associated Voronoi cell volumes: For quantifying the shapes of voids, vide calculates their inertia tensor, defined via: where we sum over all member particles with comoving coordinates x j , y j , and z j relative to the void's center. The other components of the inertia tensor are defined accordingly. The void ellipticity is then given in terms of the smallest (J 1 ) and largest (J 3 ) eigenvalues of the inertia tensor [97]: Further void properties of interest include the core densityn min and the compensation ∆ t [48]. The index 't' for tracers is either 'h' in case of finding voids in the halo distribution or 'm' for CDM, respectively. The core density is the density inside the largest Voronoi cell of a void, thus the cell of minimal density, expressed in units of the mean tracer density: The compensation is a measure of whether a void contains more or less tracer particles N t than an average patch of the Universe of same volume V . It conveys information about the environment the void is located in, whether it be one of higher or lower local average densitŷ n avg in units of the meann. Voids with ∆ t > 0 are referred to as being overcompensated and voids with ∆ t < 0 as being undercompensated [17,22]. The compensation is defined as:

Void profiles
The density profile n (i) v (r)/n − 1 of an individual void i is defined as the spherically averaged density contrast around the center of a void from its mean valuen in the Universe. When using tracer particles, the density in radial shells of thickness 2δr at a comoving distance r from the void center at the origin can be written down in its general form as: where Θ(r j ) ≡ ϑ [r j − (r − δr)] ϑ [−r j + (r + δr)] uses two Heaviside step functions ϑ to define the radial bins of the profile, with r j being the distance of tracer j from the void center. Here we sum over all tracers j in the vicinity of the void, up to our desired maximal distance. The w j in equation (3.7) are a placeholder for optional weights andw = 1 Nt k w k is the average of weights over all N t tracers. For the usual (unweighted) density profile we simply set w j = 1 for every tracer. Another option is to include mass weighting for halos as tracers, which we will explore in section 5.2.
Besides the individual density profiles of voids, we are also interested in stacked profiles. These stacks are simply an average over the individual profiles in ranges of different void properties, typically their radius r v . In order to maintain characteristic void features, such as their compensation wall, at a constant location from the void center in a given stack, and in order to better compare different stacks with each other, we calculate the profile of every individual void using constant radial bin sizes when expressed in units of their void radius, i.e. δr/r v = const. In this way every void's compensation wall is centered around r = r v , and not scattered among different distances when using constant units in physical scale. In this manner the stacked density profile is given by an average over equation (3.7): When calculating density profiles of voids using CDM tracers, we will denote these profiles by ρ v (r)/ρ − 1 instead of n v (r)/n − 1 for halo tracers, and n (w) v (r)/n (w) − 1 for mass-weighted density profiles based on halo tracers.
To investigate the radial movement of tracer particles around a void, we calculate its velocity profile around the barycenter. In our definition, positive velocities correspond to an outflow of tracers from the void, whereas velocities are negative when tracers move towards the void center. The velocity profile of each individual void i can be estimated by calculating: .

(3.9)
Here u j is the peculiar velocity vector of a given tracer j,r j = r j /r j the unit vector connecting the void center with the tracer particle j, and V j its Voronoi cell volume. By weighting the individual particle velocities with their Voronoi volumes V j , we ensure that the velocity profile is a volumetric representation of the true underlying velocity field [22]. This takes into account that a uniform sampling of velocity fields from an uneven tracer distribution is biased high in denser and biased low in emptier regions [67]. Moving on to stacked velocity profiles, there are two different ways of calculating them. The first option is equivalent to the way of stacking density profiles by simply averaging over the individual velocity profiles of each void in the sample: We will refer to stacks calculated with this method as individual stacks. Alternatively, we can average the denominator and numerator of equation (3.9) separately for all voids of the stack before taking their ratio: This method will be referred to as global stacks. There are typically many voids, in particular small ones, which contain no or very few tracer particles near their centers. Hence, as no velocity can be estimated there, it is set to zero by default. Even stacking many such individual profiles via equation (3.10) will result in no better velocity estimate. However, using global stacking one first gathers the tracers within the entire sample of voids and then divides by the normalization in equation (3.11). This guarantees better statistics for tracer counts in shells, but statistically favours larger voids, whose shell volumes are bigger in physical scale. Both methods come with advantages and disadvantages. Depending on the applications they are used for, one may be preferred over the other. We will discuss this issue in more detail in the following sections. Finally, local mass conservation allows us to relate the density profile with the velocity profile via the linearized continuity equation [106,107]: where f (z) = Ω γ m (z) is known as the linear growth rate of density perturbations. Here γ ≃ 0.55 is the growth index of matter perturbations in GR, Ω m (z) the matter density parameter, H(z) the Hubble parameter and b t the bias of tracer particles, with b t = 1 for CDM. Lastly, ∆(r) is the integrated density contrast, defined as: Combining equations (3.12) and (3.13) then yields: from which we can see the direct relation between density and velocity profiles. In section 6 we will investigate how well this equation holds up in the environments around cosmic voids.

Tracers
As void finding with vide only requires the positions of tracer particles of any kind, we will use both the distributions of halos and CDM in our simulations to identify voids in. We do not identify voids in the total distribution of matter (including baryons) in this work. From here on, we will refer to the voids identified in the CDM distribution simply as CDM voids and correspondingly voids identified in the halo distribution will be referred to as halo voids.
The corresponding numbers and number densities of halos with masses above the mass cuts of 1.0×10 12 M ⊙ /h in midres, 1.0×10 11 M ⊙ /h in highres and 1.3×10 9 M ⊙ /h in ultra-hr are shown in table 2. The midres simulations most closely match the expected tracer densities attainable with state-of-the-art galaxy surveys, such as Euclid [108].
In both midres and highres boxes we additionally perform subsamplings of the CDM particles which closely match the total number of halos used for the void finding. This is done in order to have identical mean tracer separationsr t and more specifically to find haloand CDM-defined voids of roughly similar ranges in void radii. Moreover, these matched subsamplings eliminate differences in the statistics derived from CDM and halo voids that are merely caused by different tracer number densities. In midres this equates to subsampling to around 0.066% of the total CDM particles, whereas in highres we subsample to around 0.034% of all CDM particles. The CDM voids of ultra-hr will be of no further relevance. While voids found in the three-dimensional distribution of matter are not directly accessible to observations, we nevertheless study these CDM voids to compare them to the ones identified using halos.

Voids
Here we present the isolated and merged void catalogs extracted from the distribution of halos and CDM in the midres and highres boxes at redshift z = 0.29. For details on the number of tracers used for the void finding, we refer to table 2 and section 4.1. Table 2 additionally contains the numbers of voids that were found using mass cuts and subsamplings, with no difference in void numbers between isolated and merged catalogs. We will not investigate the void properties and profiles from the ultra-hr simulation further due to the extremely low number of voids, instead we only use ultra-hr for a resolution study of linear mass conservation in section 6.
In figure 2 we present the void abundance of isolated and merged voids as a function of their radius (up to around 80 Mpc/h), ellipticity, and core density (up to a value of 1.5). The former is also known as the void size function [e.g., 18,103,[109][110][111][112]. We first note that whether or not voids are merged or isolated only significantly affects their distribution in radius, but their ellipticities and core densities are almost identical. Voids from the highres simulation reach smaller radii than in midres, because the better resolution provides more CDM particles and low-mass halos, and thus a higher overall tracer density.
Note that although the total number of halos and CDM particles are matched for each box, we observe almost twice as many CDM voids (see table 2), which implies smaller voids in the CDM for a fixed simulation volume. This difference is due to the bias of halos, which are preferentially located in regions of high matter density. CDM particles sample the density field without this bias, which enables the detection of voids in less dense regions as well. Therefore, establishing a one-to-one correspondence between halo voids and CDM voids from the same simulation is not feasible [105].
However, the void size functions of merged halo and CDM voids in both midres and highres all agree within their error bars for radii larger than about 45 Mpc/h. Figure 2 suggests that the abundance of merged voids above a certain size is much less dependent on tracer bias, subsampling fraction, mass cuts, and resolution. We additionally ran vide on a subsample of 200 million CDM tracers in midres, which confirmed this result for merged voids. No such convergence can be found in the void size function of isolated voids, which  continuously fragment into smaller ones as the density of tracers increases. However, at the smallest radii the void size function is not affected by merging. These are voids in the isolated catalog, which get relabelled to sub-voids in the merged catalog. As noted before, the total number of voids in merged and isolated catalogs is identical, since no new voids are 'created', only some voids are merged, resulting in the hierarchical structure of parent and child voids. From the lower left plot in figure 2 it is evident that the ellipticity distribution is more or less identical for voids in all cases, the curves in highres are simply shifted vertically towards higher values due to the larger number density of voids that can be found when resolving CDM particles at higher density, or halos of lower mass. All distributions consistently peak around a value of ε ≃ 0.15. There is a mild indication for merged voids to be slightly more elliptical, likely due to their increased amounts of substructure, which increases the possible complexity of void shapes. In contrast, the distributions in core density on the lower right panel exhibit clear differences across tracer type and resolution. Halo voids typically exhibit deeper core densities than their CDM counterparts due to halo bias, which amplifies fluctuations in the density field [28,105]. Moreover, the core densities show clear differences between midres and highres beyond a vertical shift of the distributions. Their maxima move towards lower density values as the resolution increases. This is because more nonlinear density fluctuations can be resolved on smaller scales. On the other hand, void merging does not affect the core densities, because the local density minima remain the same by construction.
The joint two-dimensional distributions of void properties are shown in figure 3 for the case of isolated halo voids of the midres simulation. Voids from other tracer types, merging thresholds, and resolutions have qualitatively similar distributions. Their ellipticity only slightly depends on void radius, but the distribution becomes narrower towards larger voids. Among the smallest voids we find a few highly elliptical ones with ε ≈ 1, which are most likely spurious due to the effects of tracer sparsity (see section 6.1). The core density exhibits a stronger correlation with void radius, as evident from the upper right plot of figure 3. Larger voids tend to have lower core densities with a more narrow distribution [113]. Towards small voids the distribution widens considerably, featuring some voids with minimum densities even above the mean background value. This is an expected consequence of a parameter-free void definition based on the watershed technique, which does not impose any density threshold on the interior of voids. Because it is purely topological, this method is able to identify local basins above the mean density. Finally, the distribution of void compensation with radius is depicted in the bottom plot. Its shape falls in between the ones for ellipticity and core density, exhibiting an anti-correlation between compensation and void radius. Small voids can be either overcompensated when found in very dense environments, or undercompensated when found inside larger voids. Because their distribution is very skewed towards positive compensation, however, on average small voids are overcompensated. On the contrary, large voids are preferentially undercompensated [17].

Profiles
In this section we focus on the stacked density and velocity profiles of voids. The profiles are calculated individually for each void out to five times its effective radius in bins of width 0.1× r v and then stacked (averaged) in bins of different void properties. We use isolated halo voids as our default, but investigate the impact of merging and CDM voids as well. Additionally, we calculate density profiles of halo voids using the underlying distribution of subsampled CDM. We then distinguish between their number density profiles and matter density profiles, depending on whether halos or CDM particles are used for the profile calculation, respectively. The latter case is relevant in weak lensing studies, where voids are identified via luminous tracers (e.g., galaxies) and the matter density field around them is probed by the gravitational shear of background objects [27,30]. The density around isolated voids gradually increases with smaller void size, reaching values above even the mean background density in the center of the smallest voids. Since their size is close to the mean tracer separation, one might classify these small voids as spurious. However, they are mostly embedded in environments of relatively high density [17], so the local mean tracer separation is smaller as well. They also exhibit clearly defined compensation walls around r = r v , which topologically defines them as voids. The slight increase in density     towards the center of the smallest voids is caused by the sparse sampling statistics on scales below the mean tracer separation. As they are only defined via a few particles, the density estimate from equation (3.7) returns biased results. In particular, the density is biased high when the shell volume in the denominator is very small near the void center, but happens to contain a tracer particle. When comparing CDM voids with halo voids, we observe that the latter have higher compensation walls and are slightly deeper near their center. This is consistent with the voids' core density distributions in figure 2, and the impact of halo bias amplifying the CDM density fluctuations [28,105]. However, we emphasize that a comparison of CDM voids with halo voids of the same size is not necessarily meaningful, due to their significantly different void size functions. For merged voids we find a somewhat different behavior, as shown in the lower panels of figure 4. While the density profiles of the smallest voids remain virtually   unchanged, merging gives rise to larger voids with slightly higher compensation walls and shallower cores than isolated voids, an effect already seen in reference [22]. This happens because merged voids contain sub-voids and hence sub-structures defining those. When the density profiles of merged voids are stacked, these sub-structures effectively get smoothed out, which leads to shallower cores and slightly higher compensation walls in the profiles. The density profiles of voids from the highres simulation are depicted in figure 5. In order to obtain a sufficient number of voids in each stack, we now adapt the bin widths in radius, covering a range from 4 Mpc/h to 30 Mpc/h for isolated voids on the top, and from 4 Mpc/h to 60 Mpc/h for merged voids on the bottom. We find consistent results in comparison to the midres simulation shown in figure 4, but resolve even smaller voids (as evident from figure 2). Again, we caution against directly comparing voids of a given size from different resolution simulations, similarly as for voids defined by different tracer types. Only in the regime where their void size functions converge, typically for larger merged voids, this may be meaningful.     Figure 6 presents the matter density profiles, ρ v (r)/ρ − 1, of isolated halo voids in both midres and highres, i.e. instead of using halos as tracers to calculate their profiles (as in figures 4 and 5), we now use subsampled CDM particles for this. The differences to the number density profiles of halo voids calculated using halos, n v (r)/n − 1, are depicted on the bottom panel of each plot. We notice similar trends as for the CDM void profiles, although less pronounced. The profiles are generally less dense in their centers, while at the same time their compensation walls are attenuated [28]. This is due to the CDM particles being more evenly distributed in space than the halos (i.e., halo bias), making their matter density profiles more similar to the ones of CDM voids.

Density profiles
So far we have exclusively focused on density profiles in bins of void radius. Figure 7 additionally presents the number density profiles of isolated halo voids from the midres simulation in bins of ellipticities ε, core densitiesn min , and compensations ∆ h . As before, we limit ourselves to voids with radii between 5 Mpc/h and 50 Mpc/h for the sake of comparability. Density profiles in bins of void ellipticity are shown in the top right panel of figure 7. As already seen in figure 3, the ellipticity of voids is only weakly correlated with their radius, especially in the presented range of ε between 0.0 and 0.26. Because the void size function peaks at small radii, small voids dominate in every bin of ellipticity, so their stacked profiles exhibit relatively high compensation walls. The least elliptical voids feature the sharpest walls and the deepest cores in their stacked density profiles, and vice versa. This is a simple consequence of spherical averaging in shells around the void center: if regions of similar density in its vicinity are non-spherical, they overlap with shells at different distances from it, which effectively smooths out the density profile. To account for this effect, an alternative stacking method has been proposed by the authors of reference [26] to follow the geometry of void boundaries when constructing void density profiles. This also leads to sharper compensation walls and deeper cores, in agreement with our findings.
The lower left panel of figure 7 depicts density profile stacks in bins of core density    figure 7, covering values from −1 (undercompensated) to 2 (overcompensated). Undercompensated voids tend to have the deepest cores and no clear compensation wall. In fact, the stacked density profile of the most undercompensated voids gradually increases outwards and only approaches the mean background density at a large distance from the void center. Yet, there is a hint for a void boundary at r = r v , where the slope in density becomes flatter. The higher the compensation of voids, the more pronounced their compensation wall becomes. Here the most overcom-pensated voids barely reach below the mean density in their cores, resembling the behavior of small voids. These on average exhibit higher compensation values, while large voids have a tendency to be undercompensated, as already evident from figure 3. The compensation is a measure of the environment a void is located in, since it depends on its size and total number of member particles. Undercompensated voids happen to be within underdense environments, whereas overcompensated voids in regions of higher local density. These two phenomena have been coined "void-in-void" and "void-in-cloud" scenario, respectively [101]. This defines an interesting transition point for exactly compensated voids with ∆ h = 0 [17]. Their density profile converges to the background density at the closest distance from the void center among all stacks in bins of compensation, around r ≃ 2r v . Apart from the imposed mass cuts on our halo catalogs, we have so far neglected the masses of individual halos for the calculation of void density profiles. When used as weights w j = M j for every halo j in equation (3.7), one can estimate the mass-weighted density profile around voids from all the matter contained in halos. This is different from using CDM particles for the profile estimation, because not all of the CDM is confined inside halos. Yet, it yields a special type of matter density profile around voids, which additionally probes the spatial distribution of halos depending on their mass. While individual halo masses are difficult to obtain in large-scale structure surveys, the magnitudes of their hosted galaxies can provide an observational proxy for them [114,115].

Mass weighting
In figure 8 we present the mass-weighted density profiles of halo voids from midres and highres boxes, with differences to their unweighted number density profiles on the bottom. The profiles are stacked in void radius bins identical to the ones used in figures 4 and 5.
Mass-weighted profiles are amplified with respect to number density profiles, featuring higher compensation walls and deeper cores, most notably for small voids. This implies that the  least massive halos tend to reside closer to the void centers, whereas more massive halos are rather located in the more clustered regions at the void boundary, as expected from previous studies [e.g., 49,104]. For larger voids, however, this effect diminishes. Alternatively, one can interpret the effect of mass weighting as a boost of the halo bias [115].

Velocity profiles
Having investigated the spatial distribution of CDM and halos inside voids, we now turn to their movements. Void dynamics are characterized by the coherent, radially directed flow of matter around their centers, which can be quantified via stacked velocity profiles. As described in section 3.2, we distinguish two different ways of calculating these, using individual stacks via equation (3.10) and using global stacks via equation (3.11). Figures 9 and 10 present the two methods for isolated halo voids from the midres and highres simulations with identical void radius bins as those used in figures 4 and 5 for number density profiles. The velocities of both halos, as well as the CDM are shown, with differences highlighted in the bottom panels. First of all, there is an exquisite agreement between the velocities of halos and CDM particles around halo voids for both stacking methods. This is as expected from the equivalence principle, stating that test particles in a common gravitational potential fall with the same speeds, irrespective of their mass and composition. Differences are becoming more visible when approaching the void centers, where the tracer statistics are sparser. Large voids are characterized by outflows, which steadily increase from close to zero at their center until a maximum near the compensation wall, then drop again to approach zero in the large distance limit. This is consistent with the deep and extended underdensity of matter in the vicinity of those voids, as seen in figure 6. As the compensation wall becomes more pronounced for voids of smaller size, the direction of motion features a turning point and matter flows radially inward towards the compensation wall [22]. This influx is most extreme around the smallest voids, which may eventually overcome their interior expansion and eliminate them [101]. One may argue this to be the case for the smallest voids in figure 9, but our results from the highres simulation in figure 10 show that voids of even smaller size experience internal outflows. While the most salient features in the velocity profiles are similar in both stacking methods, there are some important differences to notice in proximity of the void center. Individual stacks closely converge towards zero velocity towards the void centers, while global stacks reach finite values there. For global stacks the differences between halo and CDM velocities decrease with increasing void size, whereas the opposite trend is manifest in individual stacks, which yield the best agreement for the smallest voids. It is mainly the halo velocities that are affected by the choice of stacking method, CDM velocities are more consistent with each other. Since the same halos are used for the identification of voids, this suggests that the residuals between halo and CDM velocities are caused by the same sparse sampling effects already found in the density profiles of small voids in figures 4 and 5, as discussed in section 5.1.
In fact, the velocity profiles of CDM voids, which can be found in figures 17 and 18 of section 6, experience similar sampling artifacts, because the same tracers are used for the identification of voids and the calculation of their profiles. This becomes most severe for the voids defined by the fewest tracer particles. We should of course bear in mind that a direct comparison between CDM voids and halo voids has a limited scope, as their size functions are very different. One could also argue for other potential causes for a discrepancy between halo and CDM velocities, such as a velocity bias from halos, or some kind of nonlinear dynamics [116,117]. However, when increasing the simulation resolution in figure 10 the discrepancy between the velocities of halos and CDM diminishes, so we cannot find any evidence in favor of such effects. Moreover, we have checked that void merging does not impact any of these conclusions either.
Analogously to the density profiles in figure 7, we can choose void properties other than their radius to stack velocity profiles. Figure 11 presents this for global stacks of isolated  halo voids from the midres simulation in bins of void ellipticity ε, core densityn min , and compensation ∆ h , in addition to the bins in void radius r v . As before, we only select voids with radii between 5 Mpc/h and 50 Mpc/h for the sake of comparability, but now also compare both methods for stacking velocity profiles in the lower panels. When considering ellipticity bins in the upper right plot of figure 11, we notice that all bins experience the characteristic outflow of halos from the void center, with the maximum velocity steadily decreasing with increasing ellipticity, in correspondence with the shape of density profiles from figure 7. This is true for both stacking methods, although individual stacks result in generally lower velocities than global stacks, even at large distances from the void center. In global stacks the most elliptical voids are dominated by outflows of halos, only the more spherical voids exhibit inflows towards their compensation wall. In contrast, individual stacks feature inflows in all bins of ellipticity. This difference can be understood considering the definition of the two stacking methods in equations (3.10) and (3.11) together with figure 3, which shows that void ellipticities and radii are largely independent from each other, such that each ellipticity bin covers a wide range of void sizes. Therefore, an average over individual velocity profiles from equation (3.9) is biased towards small voids, which are more numerous than large ones. On the other hand, equation (3.11) sums up the volumeweighted tracer velocities around the full void sample before normalizing by the total Voronoi volume of all tracers per shell, which is biased towards shells containing more tracers and hence large voids. The dependence of velocity profiles on void radius, as shown in figure 9, then translates into the differences of the two stacking methods appearing in the upper right panel of figure 11.
The velocity stacks for bins in core density are depicted in the lower left panel of figure 11. The correspondence with the associated number density profiles from figure 7 is evident: the most underdense voids experience the strongest outflows and vice versa. Differences between the two stacking methods are present, but do not affect the general trends already manifest in figures 7 and 9. This is because the core densities of voids are more strongly correlated with their radii than ellipticities are, so the averaging effects discussed above are of lower importance here.
For the velocity profiles in compensation bins, as shown in the lower right panel of figure 11, the general trends are similar to the previous case, except that compensation has a stronger impact on the void environment than its core density, as expected from figure 7. Overcompensated voids are dominated by influx of matter, while undercompensated voids expand out to large distances from their center. At the boundary between these two regimes, the peculiar motions around compensated voids with ∆ h = 0 vanish at the smallest distance from the void center. Test particles near a region of average background density experience no net force and simply move with the Hubble flow.
In summary, an important conclusion from figure 11 is the fact that different estimators for stacked velocity profiles around voids can be biased in different ways, depending on the diversity of the considered void sample and the selection property used to bin the stacks. Individual stacks are biased towards the more numerous small voids, whereas global stacks are biased towards larger voids that are sampled with more tracer particles. The issue can be partially mitigated by limiting the range of void sizes per stack, but a comparison of both stacking methods is helpful in revealing residual biases the estimators may encounter.

Sampling effects
A comparison of the different profile estimators in the previous two subsections revealed some of their limitations in certain regimes. We expect these limitations to arise from the sparse statistics of tracers or voids, which are unavoidable when approaching scales near the average tracer separation at any resolution. However, we can artificially amplify the impact of sparse statistics by removing the lowest-mass halos, or by randomly subsampling the tracers used in the profile calculations. Nevertheless, we do not repeat the void identification step on modified tracer samples, because this would render a direct comparison between different void samples impossible. For example, this can be seen when comparing the number density profiles of isolated halo voids in the midres and highres simulations from figures 4 and 5. The distribution of void sizes is different in these two cases, and voids of the same size do not necessarily share the same properties. Without loss of generality we restrict ourselves to isolated halo voids from the highres simulation with radii ranging from 16 Mpc/h to 20 Mpc/h, and subsequently apply mass cuts and subsamplings to the halo sample for estimating the void profiles. Figure 12 presents the stacked number density profiles after selecting halos above different mass thresholds M min and after random subsamplings that match the previously selected number of halos. Mass cuts affect the density profiles in a similar way as mass weights (c.f. figure 8), revealing the tendency of more massive halos to be distributed within the voids' compensation wall, but being scarcer in the voids' core. On the contrary, random subsamplings have no significant effect on the number density profiles of voids, except for an increase of the error bars, as expected. This certifies that our density estimator from equation (3.7) is not biased for sparse tracer samples.
The corresponding velocity profiles are presented in figure 13. Evidently, individual stacks and global stacks are affected very differently by mass cuts and subsampling. While the velocity profiles from individual stacks continuously diminish when reducing the number of tracers, the ones from global stacks remain stable within their error bars. In this regard the global stacking method to estimate velocity profiles is preferred, because it does neither generate a bias from sparser tracer samples, nor depend on the tracer mass. A dependence on tracer mass would violate the equivalence principle adopted in the gravity solver of the simulation, so it must be spurious. The issue arises whenever tracers become too scarce to faithfully sample the velocity field, yielding a too low velocity estimate. It is particularly severe near the void center in individual stacks, where massive halos and other tracer particles are scarcest.

Linear mass conservation
After examining the stacked density and velocity profiles of voids separately in the previous section, we now want to investigate their interrelation via the linear continuity equation (3.14). For a given density profile, we refer to the velocity profiles calculated via that equation as linear theory profiles. We have to determine one free parameter in equation (3.14) when using halos as tracers: their bias b t . We achieve this by fitting the linear theory profiles to the velocity profiles estimated via equation (3.9) for individual voids, and via equations (3.10) and (3.11) for stacked voids, treating the bias as a free (inverse) amplitude.

Individual voids
We begin with testing the validity of equation ( midres simulation, while their velocity and linear theory profiles (with b t = 1) are shown on the right. These five voids are a random draw from our catalog, with the only condition to sample a range of different void radii. Although the individual density profiles are subject to high sample variance, one can perceive the characteristic underdensity near the voids' center, a compensation wall around r v , and the trend towards the mean background density at large scales.
Their velocity profiles mirror this fluctuating structure, with some voids dominated by outward motion on all scales and others experiencing infall towards their compensation wall from large distances. One particular void with r v = 32.06 Mpc/h (in yellow) exhibits a density peak in its very center, which serves as a good example for the sparsity effects that can occur in the estimator from equation (3.7) at the smallest inner shells. Nevertheless, the linear theory profiles match the measured velocity profiles with a remarkable accuracy in every single case. Due to the integral over the density profiles appearing in equation (3.14), the linear theory profiles become smoother than the measured velocity profiles. Near the void centers the differences are more significant due to the tracer sparsity inside voids, but still consistent with the scatter of the measurement. We emphasize that no free parameters have been adjusted in this procedure, the tracer bias of CDM is simply fixed to b t = 1 here.
When using halos as tracers we can no longer assume b t = 1, but have to determine their bias separately. The individual number density profiles of five randomly selected isolated halo voids of the midres simulation are presented in figure 15. We additionally revisit mass-weighted density profiles on the bottom panel for the same five voids, as introduced in section 5.2. The two types of number density profiles follow similar shapes, but the mass-weighted ones have a tendency for higher peaks and deeper troughs, because more massive halos tend to reside in regions of higher density, as discussed before. One void with r v = 23.7 Mpc/h (in red) is embedded within a larger underdensity located between r = 3r v and r = 4r v , a nice example for the void-in-void scenario. The individual velocity profiles reflect the density structure of all five voids. We now leave the bias b t in equation (3.14) as a free parameter to fit the linear theory profiles to the measured velocity profiles. Again, the match between the two is striking, even for the mass-weighted case. The only difference from the latter is an enhanced value of the best-fit bias parameter, which corroborates our conclusion from section 5.2.
For completeness we present additional profiles from five isolated halo voids of the highres simulation in figure 16. This time we intentionally selected four voids of similar size to point out that their profiles can still be very diverse. The void with r v = 19.71 Mpc/h (in red) could be ascribed to a "void-in-void-in-void" scenario, exemplifying how deeper hierarchies with multiple levels of sub-voids can occur. This is substantiated by its velocity profile, which exhibits three distinct peaks of enhanced outflow near the boundaries of each sub-void in the density profile. No matter how complex the void structure is, linear theory successfully reproduces the dynamics within each individual void. In the highres simulation the only difference to midres is that lower bias values are obtained, due to the lower mass cut of 10 11 M ⊙ /h. To summarize, mass conservation based on the linear continuity equation (3.14) provides an extremely accurate description for the dynamics operating around individual voids that are well resolved, independently of whether CDM or halos are used as tracer particles for void finding. To the best of our knowledge, this unique behavior has not been found in any environment of large-scale structure before.

Stacked voids
In the same fashion as performed for individual voids, we may apply the linear continuity equation (3.14) for stacked voids. To this end we can make use of the previously examined matter density profiles of CDM voids, number density profiles of halo voids (with or without With all these variants there are various options to test linear dynamics with stacked void profiles. Reference [22] already presented such a test based on individually stacked profiles of merged CDM voids, which we repeat here with isolated CDM voids for both stacking methods in figure 17. We use identical void radius bins as used in the upper left panel of figure 4. Furthermore, we treat b t as a single free parameter even for the case of CDM. In that case, deviations from b t = 1 can indicate additional biases of estimators. We confirm a very good agreement between the linear prediction and the velocity profiles for both individual and global stacks, especially for larger voids. The profiles from global stacks exhibit a more linear rise near their center, which is slightly better in line with linear theory than in the case of individual stacks. Deviations become stronger for small voids, especially close to their center due to tracer sparsity, as discussed in section 5.3. However, also at larger distances from the void center a constant offset remains for small voids. The best-fit values of the tracer bias parameter converge towards 1 for large voids, as expected for CDM tracers. Smaller voids exhibit values exceedingly higher than unity, in accordance with the results of reference [28]. These conclusions are not affected by the choice of merging threshold and remain valid for halo voids and their number density and velocity profiles, which yield higher values for the tracer bias, as expected from halos. With the methodology developed for individual void profiles in section 6.1, we can now introduce two additional tests for linear dynamics around stacked voids. The first option is to apply equation (3.14) to each individual density profile of a given void sample, fit b t for every void profile separately, and average the linear theory profiles and best-fit bias parameters in the end. We will refer to this method as individual fits. In the second option we also apply equation (3.14) to the individual void density profiles, but fix b t to 1. Then we average the resulting linear theory profiles and fit for a single 'global' tracer bias parameter in the end.  This method will be referred to as global fits. Both methods have the advantage that they provide an entire posterior distribution of predictions for the linear theory profile, which can be used to quantify uncertainties (error bars). Figure 18 presents the results from these two methods (individual and global fits) for isolated CDM voids. In all cases the measured velocity profiles and linear theory predictions match closely, with stronger differences again occurring around small voids. Moreover, differences between data and model gradually diminish towards larger voids in global stacks, whereas in individual stacks those are of the same magnitude for all void radius bins. In general, the results from global fits are very similar to the previous case from figure 17. In contrast, individual fits feature a slightly better agreement between data and model, particularly at larger distance from the void center. On the other hand, the best-fit values of the  Figure 19. Same as figure 18, but for isolated halo voids in the midres simulation.
tracer bias are somewhat closer to unity in the case of global fits. Note that for individual fits the obtained bias values are identical in both stacking methods, since the fitting is performed before stacking. For global fits the bias values slightly vary between the two stacking methods. Subsequent figures 19 and 20 present individual and global fits to the profiles of isolated halo voids in both the midres and highres simulation. Essentially all our conclusions from CDM voids remain valid for halo voids. Only the tracer bias parameters are higher, as expected, but also decrease towards larger void sizes. A comparison across different resolutions allows us to explore a wider range in void size and potentially identify a limit for the validity of linear theory for voids below some characteristic size. For the highres simulation, we use the following void radius bin edges:  Figure 20. Same as figure 19, but for the highres simulation.
12 Mpc/h exhibit smaller residuals between data and model than in midres. In order to reach the same level of agreement in midres, one has to look at much larger voids of around 22 Mpc/h in radius. Voids of that size are among the largest in highres and exclusively feature outflows, similar to the largest voids above 40 Mpc/h in midres. This is a strong indication for the fact that residual mismatches are due the sparsity of tracers in the midres simulation, rather than a limit of linear theory below some fundamental scale. As a final test we exchange the number density profiles of halo voids with their massweighted versions. This is presented in figure 21 for the midres simulation. The theoretical model still produces the correct shape of stacked velocity profiles, but residuals increase somewhat in comparison to using number density profiles. In individual stacks, the differences now increase with void size, although the mass-weighted density profiles of large voids are the most similar compared to their number density profiles (see figure 8). In either stacking method we observe that individual fits experience the largest discrepancies inside voids, where the maximal velocity of linear theory profiles considerably falls short of the measured  Figure 21. Same as figure 19, but using mass-weighted density profiles in equation (3.14).
velocities. On larger scales the profiles align again, except for the smallest voids. Global fits instead reproduce the maxima of the velocity profiles more accurately, but disagree more near the compensation walls. In all cases, the retrieved tracer bias values increase, as expected from mass weighting, but the agreement with linear theory is generally worse than for the unweighted case. This can be explained by the fact that this weighting scheme amplifies the impact of high-mass halos, which are much scarcer than the bulk of the halo population. In turn, this leads to a stronger sensitivity to the impact of tracer sparsity in the profile estimators, most notably close to the void centers. We have repeated the same analysis on the matter density profiles of halo voids and merged voids. Furthermore, we tested higher mass cuts when selecting halos for our void finding and using those halos for calculating the profiles. We also performed tests at different redshifts (z = 0 to z ≃ 4) and for a variety of cosmological parameter values in midres simulations of smaller box size, available within the Magneticum suite. We find identical conclusions, so we refrain from presenting this in additional figures here. One also has the  option to stack voids based on properties other than their radius, like in figures 7 and 11. In that case the linear theory profiles often yield a less accurate prediction due to the large ranges in void size per stack. However, this can be remedied using the methods of individual fits and individual stacks, because then voids are modeled independently from each other.

Resolution study
Our tests on the linear continuity equation (3.14) in the previous subsections show that velocity profiles calculated with this model agree remarkably well with measured profiles, both for individual, as well as for stacked profiles of both CDM and halo voids. A comparison between midres and highres boxes has not revealed any limitations for the validity of the model, apart from the impact of tracer sparsity. Here we extend this resolution study by one more layer, making use of the ultra-hr simulation. It covers a box of merely 48 Mpc/h a side, but features a much higher resolution than highres. We find a total of 346 isolated halo voids inside this simulation. This is too small of a sample to consider for realistic cosmological applications, but still very suited for testing linear dynamics at very small scales. Figure 22 shows both individual and stacked density profiles of voids on the left. The latter clearly exhibit larger error bars and therefore more noise than the previously analysed profiles, owed to the much smaller sample size. Nevertheless, the characteristic features are consistent with our previous findings, except that here we are examining voids of only a few Mpc in size, with a maximum radius of 4.5 Mpc/h. The right side of this figure presents the corresponding velocity profiles of the same voids, with individual stacks on the bottom. Most of these small voids experience infall from their environment, but there is one example of a clearly undercompensated void with r v = 4.15 Mpc/h (in green).
Based on their density profiles on the left side of figure 22, we repeat the application of equation (3.14) to predict the velocity profiles on the right. From the individual profiles on the top we find that even the smallest depicted void with radius 1.38 Mpc/h (in red) is consistent with linear dynamics within the scatter. This is still in line with the individual void profiles in both midres and highres from section 6.1. Despite their huge variety in shape, these individual profiles are all described remarkably well by the linear continuity equation. The stacked profiles on the bottom fully confirm this picture: even for the smallest bin in void radius between 1.0 Mpc/h and 2.5 Mpc/h, linear theory profiles (based on individual fits) align with stacked velocity profiles within their error margins on all scales around the void center.
Such scales, reaching 1 Mpc/h and below, are typically considered as highly nonlinear in the field of large-scale structure. Nevertheless, we find no evidence for the onset of nonlinear dynamics in void environments of that size. It is possible that this is a general characteristic of voids and holds irrespective of scale. At the same time, voids dominate the volume fraction of the Universe, which implies that linear dynamics should prevail within the largescale structure. It is conceivable that for practical reasons the tradition of studying the brightest galaxies (such as luminous red galaxies) and the most massive objects (such as galaxy clusters), which are located in the densest environments of the cosmos, may have concealed this insight from cosmologists so far.

Conclusion
In this paper we have explored the interrelation between a variety of general void properties across a substantial dynamic range in mass and scale by analyzing the hydrodynamical simulation suite Magneticum. This enabled us to reveal and inspect a number of universal characteristics of voids that persisted in all considered settings. Our main results can be summarized as follows: • Merging creates larger voids with a hierarchy of sub-voids. Stacking merged voids leads to shallower central density profiles, due to sub-structure inside these voids, and slightly higher compensation walls (figures 4 and 5). The smallest voids are not affected by this, because they host no sub-voids at a given tracer density. They are identical to isolated voids. In contrast to isolated voids, the size function of merged voids converges on large scales, irrespective of resolution and tracer type (figure 2).
• The stacked number density profiles of halo voids exhibit the same characteristics as the matter density profiles of CDM voids, albeit with an enhanced amplitude caused by the bias of halos (figures 4 and 5). Using their masses as weights in the density profile calculation even enhances this effect, as halo masses correlate with the density of their environment ( figure 8). However, the CDM density profiles of halo voids are very similar to the ones of CDM voids ( figure 6).
• The ellipticity of voids affects the shape of stacked density and velocity profiles via spherical shell averaging. More elliptical voids exhibit smoother density profiles with wider compensation walls and shallower cores (figure 7).
• The minimum and average density inside voids is anti-correlated with the height of their compensation wall: the shallowest voids feature the biggest walls and vice versa. Overcompensated voids, which live in overdense environments, are typically much smaller than undercompensated ones, which can be interpreted as sub-voids inside larger parent voids. These void properties are reflected in the general shape of their density and velocity profiles (figures 7 and 11).
• The stacked velocity profiles of voids mirror their density structure: large and undercompensated voids are dominated by coherent outflow, while small and overcompensated ones are dominated by infall towards their compensation wall. Halos and CDM move at the same speed around voids, as expected from the equivalence principle in general relativity (figures 9 and 10). Their radial velocity profiles accurately obey linear mass conservation. Residual deviations are caused by sampling artifacts arising from sparse tracer statistics. These become increasingly significant inside voids with extensions close to the resolution limit / mean tracer separation and depend on the type of estimator used for the velocity profile calculation (figures 17-21).
• Mass cuts and subsampling affect different estimators for void profiles in different ways. Density profiles are stable against subsampling, while mass cuts have a similar effect as mass weighting. Global stacks for velocity profiles are more independent of tracer mass and subsampling than individual stacks (figures 12 and 13).
• All previous conclusions on stacked profiles remain valid even for individual voids, although their profiles are more affected by sampling variance and noise. Nevertheless, the interrelations between density and velocity profiles of individual voids obey linear mass conservation with an exquisite accuracy (figures [14][15][16]. This fact can be exploited to explain and improve the linear relationship between the corresponding average profiles of stacked voids, which plays an important role in cosmological applications. • The apparent breakdown of linear dynamics inside the voids we analyze at moderate resolution is caused by sparse sampling of tracers close to their mean separation scale. We have confirmed this by comparing voids of a given size in simulations of increasing resolution and find no sign for the onset of nonlinear dynamics down to scales of order 1 Mpc/h ( figure 22). At a given resolution, velocity profiles are more severely affected by sparse sampling than density profiles, so the linearized continuity equation applied to densities is better suited for the modeling of velocities and observable RSD around voids than using more biased velocity profile templates measured in simulations.
These findings have a number of important implications for observational studies on cosmic voids and their use as cosmological probes. While it is not yet feasible to identify voids in the full three-dimensional matter distribution of the Universe, various tracers of the latter have already been considered for void finding, such as galaxies [e.g., [118][119][120][121][122], galaxy clusters [123], the Ly-α forest [124][125][126][127], or the 21cm emission from neutral Hydrogen [128,129]. We cannot foresee a compelling reason as to why our conclusions above would not apply for the voids identified within any of these tracers. This opens up a vast observational window to conduct cosmological experiments, allowing us to make use of voids ranging from a few (or possibly less) to hundreds of Mpc in diameter. As linear dynamics holds up in all voids irrespective of their size, void catalogs can not only be extended via enlarging survey volumes towards higher redshifts (as planned for a number of next-generation surveys, such as Euclid [108,130,131]), but also by conducting deeper observations that provide denser tracer samples and hence smaller sub-voids (e.g., as planned for Roman [132] and 4MOST [133]), even at low redshift. This would open up the possibility to maximize the number of observable linear modes of the density field of large-scale structure that can be exploited for the purpose of cosmological inference, for example via the Alcock-Paczynski effect and RSD [60,66,108], far beyond the previously imposed limits. A similar conclusion has been reached in reference [31], where voids with r v > 5Mpc/h have been found to behave very linearly in the sense that they are well-described by the Zel'dovich approximation. Moreover, the upcoming data sets will contain hundreds of thousands of voids with numerous individual properties, providing a rich playground for the latest machine learning applications [52,134,135].