© The Authors 2025

1. Introduction

The connection between galaxies and the dark matter halos they inhabit is a cornerstone of our understanding of galaxy formation and evolution within the large-scale structure of the Universe (White & Rees 1978). The halo occupation distribution (HOD), which statistically describes the number of galaxies residing in halos of a given mass, has become a powerful tool for exploring this connection (e.g., Jing et al. 1998; Ma & Fry 2000; Peacock & Smith 2000; Seljak 2000; Scoccimarro et al. 2001; Berlind & Weinberg 2002; Cooray & Sheth 2002; Berlind et al. 2003; Zheng et al. 2005; Yang et al. 2007; Rodriguez et al. 2015).

The halo mass is the primary factor in determining a galaxy’s characteristics and population, and initial HOD models thus often assumed a primary dependence on halo mass. That said, recent studies show that the properties of the large-scale environment in which these halos reside also exert an influence (e.g., Shi & Sheth 2018; Montero-Dorta & Rodriguez 2024), including on the way galaxies populate dark matter halos (Zehavi et al. 2018; Artale et al. 2018; Alfaro et al. 2020, 2021, 2022; Perez et al. 2024, among others). This environmental influence becomes more pronounced when focusing on regions with matter densities that significantly deviate from the mean. Specifically, marked deviations in the HOD are found in extreme environments, such as the densest regions and cosmic voids (Alfaro et al. 2020, 2021; Perez et al. 2024).

Cosmic voids are characterized by their extremely low densities and reduced gravitational interactions and provide unique environments for studying galaxy evolution (Rojas et al. 2004; Hoyle et al. 2005, 2012; Patiri et al. 2006; Ceccarelli et al. 2008). Previous studies have demonstrated that galaxies in voids are typically bluer, fainter, and have younger stellar populations with intense star formation activity (Rojas et al. 2004, 2005; Hoyle et al. 2005, 2012; Jian et al. 2022; Rodríguez-Medrano et al. 2023).

Focusing on HOD in voids, Alfaro et al. (2020) found that halos within these regions are systematically populated by fewer galaxies. This difference becomes more pronounced for higher-mass halos. These results were consistent across two different types of simulations: a semi-analytic model and a hydrodynamical simulation. Building on this work, Alfaro et al. (2022) conducted a similar analysis using observational data from the Sloan Digital Sky Server Data Release 12 (SDSS-DR12) spectroscopic catalog (Alam et al. 2015) and the galaxy group catalog of Rodriguez & Merchán (2020). Despite the inherent observational biases, the results supported the predictions from simulations, demonstrating the influence of these vast under-dense regions on how galaxies populate dark matter halos.

The aforementioned studies provide evidence of the impact of voids on the evolutionary history of halos, but these differences were analyzed at redshift z = 0. A natural question that arises from this context is how the HOD in these environments evolves over time. Using the semi-analytic MultiDark Planck 2 Semi-Analytic Galaxies (MDPL2-SAG) catalog (Knebe et al. 2018), we explored the environmental dependence of the HOD in greater depth by focusing on its temporal evolution within cosmic voids. By studying the HOD at different epochs, we aim to understand the underlying mechanisms driving this process. We constrained our analysis to between z = 0 and z = 0.5 for two reasons. First, we did not expect a remarkable evolution in the void size function of our different redshift samples, which would ensure consistent region identification. Second, observational limitations increasingly complicate void and galaxy group detection at higher redshifts because of declining tracer density and increasing selection effects. These constraints motivated our focus on the local Universe and provided a robust baseline for understanding void galaxy populations.

This paper is organized as follows: In Sect. 2 we describe the data used in this work, including the MDPL2-SAG catalog and our void identification methodology. Section 3 presents our analysis of how cosmic voids affect the HOD across different redshifts, examining both the global impact of voids and their specific effects on z = 0 void galaxies. In Sect. 4 we analyze the formation and environmental evolution of void galaxies at z = 0. We examine the formation times of their host halos, their spatial distribution within voids, and the evolution of their local density. These results provide insight into how the void environment influences galaxy formation and help explain the distinct HOD trends observed in these regions. Finally, in Sect. 5 we summarize our main findings and present our conclusions.

2. Data

In this section we describe the datasets used in this work to study the temporal evolution of the HOD in cosmic voids. For this, we employed a catalog of semi-analytic galaxies (SAGs) and their corresponding dark matter halos, derived from the MDPL2 simulation and the SAG model.

2.1. MDPL2-SAG catalog

The galaxy catalog used in this work is the MDPL2-SAG catalog (Knebe et al. 2018), publicly available in the database COSMOSIM¹ (Riebe et al. 2013). The catalog was built using the MDPL2 cosmological simulation (Prada et al. 2012; Klypin et al. 2016), which follows the evolution of 3840³ dark matter particles in a cubic box of 1000 h⁻¹ Mpc on side, with cosmological parameters consistent with Planck observations (Planck Collaboration XIII 2016). Dark matter halos and subhalos were identified using the ROCKSTAR algorithm (Behroozi et al. 2013a), and their merger trees were constructed with CONSISTENTREES (Behroozi et al. 2013b). Dark matter halos were populated with galaxies using the semi-analytic model of galaxy formation and evolution SAG (Cora et al. 2018). This model includes the main physical processes relevant to galaxy formation, such as radiative cooling of gas, star formation, supernova feedback, chemical enrichment, growth of supermassive black holes, active galactic nucleus feedback, and galaxy mergers. The resulting catalog, known as MDPL2-SAG, has also been used in our previous studies of galaxy formation and evolution (Alfaro et al. 2020, 2021). To study the evolution of galaxies and halos over time, we select 10 snapshots from the simulation, ranging from snapshot 125 to 107, which covers a redshift range from z = 0.0 to z = 0.49. The details of these snapshots, including the total number of galaxies and the identified voids in each, are presented in Table 1.

To ensure the resolution of the objects, we selected in each snapshot only those galaxies that reside in halos with M_200c greater than 5 × 10¹⁰ h⁻¹ M_⊙, where M_200c is the mass contained in a sphere enclosing 200 times the critical density of the Universe, ρ_c = 2.7755 × 10¹¹ h² M_⊙ Mpc⁻³. The MDPL2 simulation uses dark matter particles with a mass of 8.271 × 10⁸ h⁻¹ M_⊙, and the halos we analyzed contain a minimum of 50 particles, ensuring reliable mass estimates and structural properties. In addition, we only considered galaxies with stellar masses (M_⋆) greater than 5 × 10⁸ h⁻¹ M_⊙ and with absolute magnitudes in the r band brighter than M_r − 5log₁₀(h) = − 16. This selection guarantees a sample of well-resolved galaxies with reliable information along all the range of redshift selected. It is important to note that in our analysis we exclude type 2 galaxies, also known as orphan galaxies. These galaxies are those that have lost their subhalo and whose tracking in the simulation requires special treatment. In the MDPL2-SAG catalog, the galaxy type classification follows a standard hierarchical criterion: type 0 (central galaxies) are defined as the most massive galaxy within a given dark matter halo, residing at the halo’s center of mass. Type 1 (satellite galaxies) are all other galaxies within the same halo that are less massive than the central galaxy. This ensures that each halo contains exactly one central galaxy, with any additional galaxies classified as satellites. Since our main focus is the study of the HOD, we restricted ourselves to type 0 (central galaxies) and type 1 (satellite galaxies). Information on the type of galaxy was obtained directly from the MDPL2-SAG catalog. Our final samples are presented in detail in Table 1.

2.2. The cosmic voids catalogs

We identified spherical voids in each snapshot of the MDPL2-SAG catalog using the public algorithm SPARKLING² (Ruiz et al. 2015, 2019). This is the same algorithm used in Alfaro et al. (2020, 2022), which enables a direct comparison of the results. Briefly, the algorithm works as follows:

• First, the density field is estimated using a Voronoi tessellation of the galaxy distribution within each snapshot, using galaxies as tracers.

• For each Voronoi cell, a density is calculated as the inverse of the cell volume, and a density contrast (δ_cell) is defined relative to the mean density of the tracers ($\overline{\rho }$).

• All Voronoi cells that satisfy δ_cell < 0 are selected as under-dense regions and potential void centers.

• Around these centers, the spherical regions grow until the integrated density contrast Δ(r) = δ(< r) reaches a threshold value Δ_lim on a given scale r = R_void, defined as the void radius.

• The algorithm then performs a sampling procedure to refine the void center location, seeking to maximize the void radius. This involves making random displacements to the center, covering the volume of the corresponding sphere associated with the Voronoi cell, and accepting the new position if the new void radius is larger.

• Finally, overlapping void candidates are removed by keeping the largest voids that do not superpose with any other candidate.

For each snapshot, we used the brightest galaxies, specifically those with absolute magnitudes brighter than M_r − 5log₁₀(h) = − 20, as tracers for the density field. The voids were defined using a density contrast threshold of Δ_lim = −0.9, which means that the identified voids have an integrated density that is, at most, 10% of the mean density of the tracers. This choice is motivated by the fact that identified voids with this value of Δ_lim exhibit the most significant differences in their HODs compared to the general halo population (Alfaro et al. 2020).

By applying this algorithm to each snapshot, we obtained an initial sample of cosmic voids. To construct the final void sample, we examined the void radii distribution in each snapshot and identified the minimum radius at which the distribution of the void radii starts to decrease. This decrease indicates that the void sample is not complete for smaller radii due to resolution effects, something expected taking into account the shot-noise predominance in low density small scales. Therefore, for each snapshot, we only considered voids with radii larger than this identified minimum radius. As an example, Fig. 1 shows the void size distributions for snapshots 125 and 107 (z = 0 and z = 0.49, respectively) and the minimum radius in each case. Taking this radius cut in each void sample into account, we obtained our cosmic void catalogs to study the HOD redshift evolution within these low-density environments. The final number of voids identified in each snapshot and the total number of void galaxies are shown in Table 1.

Fig. 1. z = 0 (black lines) and z = 0.49 (red lines) void size distribution. The dashed vertical lines mark the radius up to which the sample can be considered complete.

3. Effects of cosmic voids in the HOD across redshifts

3.1. The global effect of the cosmic voids

To understand how the global effect of voids on the halo population changes over time, we analyzed void galaxies and halos at different redshifts. It is important to clarify that in this part of the work we were not tracking the same objects across snapshots using merger trees. Instead, we independently identified cosmic voids at each redshift and analyzed the properties of galaxies and halos that reside within these voids at each specific epoch.

We considered snapshots from z = 0 to z = 0.49 (snapshots 125 to 107; see Table 1) of the MDPL2-SAG simulation. For our samples in each snapshot, we selected all halos whose central galaxy is located within the boundaries of a cosmic void. We then included all the galaxies associated with these halos.

Following the same procedure used in Alfaro et al. (2020) to determine the HOD within the voids, we performed a stacking analysis of the selected halos in each snapshot. This allowed us to compute the average number of galaxies as a function of halo mass within the void environment, ⟨N|M_halo⟩. We then compared this measurement with the overall HOD measured in the entire simulation box. This provided us with a clear picture of how the presence of voids affects the halo occupation at different epochs.

Figure 2 shows the HOD measurements for galaxies with M_r − 5log₁₀(h) <  − 17 at redshifts z = 0, z = 0.304, and z = 0.49. We note that the calculations were performed for the complete sample of snapshots, but for simplicity, only these three redshifts are displayed. The same applies to the limiting values of M_r − 5log₁₀(h), where we explored all the results presented in this work for the magnitudes M_r, lim = −16, −18, −19, −20, and −21, consistently finding coherent values. We conclude that the detected effects impact the overall galaxy population regardless of their magnitude. Therefore, we chose to present the results only for M_r − 5log₁₀(h) <  − 17, as this sample contains a large number of well-defined objects that are not near the cutoff values used to define the catalogs. In this and all subsequent graphs that feature error bars, the bars represent the standard error estimated by the Jackknife resampling method. In each case, we used 50 subsamples, since the jackknife estimates stabilize once at least 30 subsamples are employed.

Fig. 2. Top panels: HOD measured for galaxies with Mr − 5log10(h) <  − 17 at redshifts z = 0, z = 0.3, and z = 0.5 for the MDPL2-SAG simulation. The blue lines represents the HOD of the complete catalogs, while the black and red lines show the HOD of halos inside the cosmic voids identified in each snapshot. Bottom panels: Ratio of the HOD inside voids to the HOD in the complete catalogs at each redshift. The error bars were estimated via the Jackknife resampling method and using 50 subsamples.

In each panel of Fig. 2, we compare the HOD measured for halos inside voids with the HOD measured for all halos in the corresponding snapshot. The bottom panels show the ratios of the HOD within voids to the overall HOD for the same three redshifts. The HOD within voids is systematically lower than the overall HOD at all redshifts, which agrees with previous results of our previous work showing that the number of galaxies inside voids is lower than in other environments. In all cases, up to ∼10¹² h⁻¹ M_⊙, the halos inside the voids do not show significant differences with the redshift compared to the global result. However, for halos more massive than ∼10¹² h⁻¹ M_⊙, the differences between the HOD in voids and overall increase with the halo mass for all redshifts, which are larger for larger masses. It is also worth noting that the number of halos comprising the two last mass bin samples decreases when we consider higher redshifts, to the point where the sample corresponding to z = 0.49 does not reach the last mass bin.

These results show that the effect of the void environment on the HOD is consistent across all analyzed redshifts, suggesting that the mechanisms that cause the HOD suppression within voids are similar through cosmic time.

3.2. The effect of the cosmic voids on z = 0 void galaxies

In the previous section, we analyzed the global effect of voids on the HOD, finding that this global effect persists with redshift. Now, we focus on how the voids influence the galaxy population residing within them at z = 0, and trace the evolution of this particular population over time.

To examine the impact of voids on the HOD of their member galaxies, we tracked the evolution of central galaxies identified within voids at z = 0. This was achieved using the GalaxyStaticID parameter from the MDPL2-SAG catalog, which remains static across all simulation snapshots and allows for consistent galaxy tracking throughout cosmic time. Unlike the previous section, where the HOD of halos was analyzed independently at each redshift, this approach focuses on a fixed population of galaxies traced across snapshots.

For each snapshot, we computed the HOD of halos hosting central galaxies that were identified within voids at z = 0 and remain central in subsequent snapshots. By fixing M_halo and M_r − 5log₁₀(h) to their z = 0 values, we isolated the HOD evolution from effects such as halo mass accretion and stellar evolution. This methodology provides a consistent view of how the HOD evolves over time for a single, defined galaxy population.

Figure 3 shows the HOD of central galaxies in voids at z = 0 across different redshifts for galaxies with M_r − 5log₁₀(h) <  − 17. The top panel presents the HOD curves in the ten selected snapshots between z = 0 and z = 0.49, using a black line for the sample at z = 0 and different shades of red that become lighter at higher redshifts. The bottom panel shows the relationship between the HOD at each redshift and the HOD at z = 0, highlighting the relative changes over time. The figure reveals distinct trends depending on the halo mass. For halos with M_halo ≲ 2 × 10¹¹ h⁻¹ M_⊙, the HOD remains nearly constant across all redshifts, indicating that these low-mass halos do not significantly accrete satellite galaxies over time. However, beyond this mass threshold, the HOD curves start to grow as the redshift decreases, reflecting the progressive accretion of satellite galaxies. This growth becomes more pronounced for halos with higher masses, highlighting the increasing efficiency of massive halos in the capture of satellites over cosmic time. This is a consequence of the hierarchical structure formation scenario, where massive halos grow predominantly through mergers and accretion at later times (Cooray & Sheth 2002).

Fig. 3. Evolution of the HOD of central galaxies in voids at z = 0. Top panel: HOD for central galaxies in voids at z = 0 across different redshifts. Bottom panel: Ratio between the HOD at each redshift and the HOD at z = 0. The error bars were estimated via the Jackknife resampling method and using 50 subsamples.

To compare the HOD of central galaxies in voids with the general behavior of halos, we computed the overall HOD by considering the complete box of each snapshot. We selected all central galaxies at z = 0 that remain as central in each snapshot, fixed their M_halo and M_r − 5log₁₀(h) values to their z = 0 properties, and calculated the HOD for galaxies with M_r − 5log₁₀(h) <  − 17. Figure 4 shows the relationship between the HODs in the previous figure and the general HOD at each redshift. To avoid line overlap and highlight the main trends, we present only the results for three snapshots: z = 0, z = 0.304, and z = 0.49.

Fig. 4. Ratio between the HOD of central galaxies located in voids at z = 0 and the general HOD of the Universe at three different redshifts. The solid black line with diamonds corresponds to z = 0, the dashed red line with squares to z = 0.304, and the dashed pink line with triangles to z = 0.49. The HOD for all cases was calculated for galaxies with Mr − 5log10(h) <  − 17 at z = 0 with halo masses fixed to their z = 0 values. The error bars were estimated via the Jackknife resampling method and using 50 subsamples.

Figure 4 shows that, in general, galaxies residing in voids at z = 0 maintain a lower HOD than the average across the entire redshift range studied. No systematic evolution is observed in the relationship between the HOD of this particular sample of central galaxies and the general trend.

Finally, Fig. 5 illustrates the relative evolution of the HOD in cosmic voids compared to the general Universe in different redshift and halo mass ranges. We quantified this evolution through a differential metric (Δ) that compares the normalized HOD in the voids against the cosmic mean, using z = 0 as the reference epoch. The differential metric was computed as

Fig. 5. Variation in the HOD within voids and in the complete snapshot between each redshift compared to the HOD at z = 0. The errors were estimated via the Jackknife resampling method and using 50 subsamples.

$$\begin{array}{cc}\hfill \mathrm{\Delta }& =[(1-\frac{{\mathrm{HOD}}_{\mathrm{all}}}{{\mathrm{HOD}}_{\mathrm{all},z0}})-(1-\frac{{\mathrm{HOD}}_{\mathrm{voids}}}{{\mathrm{HOD}}_{\mathrm{all},z0}})]\times 100\hfill \\ \hfill & =(\frac{{\mathrm{HOD}}_{\mathrm{voids}}}{{\mathrm{HOD}}_{\mathrm{all},z0}}-\frac{{\mathrm{HOD}}_{\mathrm{all}}}{{\mathrm{HOD}}_{\mathrm{all},z0}})\times 100,\hfill \end{array}$$(1)

where each term in parentheses on the first line of the equation represents the corresponding percentage with respect to the HOD value at z = 0 for each mass bin. This normalization allows us to track relative changes in the galaxy occupation of halos while accounting for the natural evolution of the overall HOD. The uncertainties were derived through jackknife resampling and proper error propagation.

The heatmap reveals several notable trends. First, a clear mass-dependent evolution emerges, with stronger differential effects at higher halo masses across all redshifts. This finding further supports the idea that the void environment exerts a more significant influence on the galaxy occupation of massive halos. Secondly, a redshift evolution is evident, where the magnitude of Δ generally increases with redshift for halos with M_halo ≳ 10¹² h⁻¹ M_⊙, but remains practically unchanged for less massive halos. This indicates a growing disparity in galaxy occupation patterns between low- and high-mass halos, driven by the void environment as we look back in time.

4. Void galaxies at z = 0: Formation and evolution of the local environment

Building on the previous section’s investigation of the impact of voids on the HOD and the properties of galaxy populations whose central galaxy resides within a void at z = 0, this section offers a more detailed examination of their characteristics. Our goal is to understand how these galaxies have been shaped by the void environment and to explore the underlying reasons for their significantly lower HODs compared to the broader galaxy population.

We began our analysis by investigating the formation times of halos hosting central galaxies within voids at z = 0. This step is essential to understanding the evolutionary history of these systems and how their development might be influenced by the void environment. For halos within voids, we computed two key formation metrics: the formation redshift (z_form), defined as the redshift when half of the maximum halo mass was accreted, and the redshift of the last major merger (z_lastMM), defined as the redshift at which the last major merger occurred (mass ratio > 0.3). Both parameters are readily inferred from the halfmass_scale and lastmm_scale outputs provided by the ROCKSTAR algorithm. Figure 6 presents these results, showing the formation times as a function of the halo mass. The top panel displays the mean values of z_form, while the bottom panel illustrates the distribution of z_lastmm. A clear trend emerges from these results: the less massive halos within voids are, on average, younger than their counterparts in the general population, whereas the most massive halos in voids are older than the average.

Fig. 6. Halo formation redshifts as a function of the halo mass. Top panel: Results for zform, defined as the redshift when half of the maximum halo mass was accreted. Bottom panel: Results for zlastmm, the redshift at which the last major merger occurred. The blue lines represent the complete z = 0 snapshot and the black line the z = 0 void halo samples. The error bars were estimated via the Jackknife resampling method and using 50 subsamples.

This finding, when correlated with the HOD results from the previous section, suggests that the void environment exerts distinct influences on halos of different masses. Additionally, the formation history of halos appears to play a crucial role in shaping the galaxy populations they host. These observed trends in halo formation times and merger histories align with recent studies of galaxies in under-dense environments. For example, Rosas-Guevara et al. (2022) identified analogous mass-dependent patterns in the timing of merger events across different cosmic regions using hydrodynamical simulations, noting that major mergers in massive-void galaxies tend to occur more recently toward the inner void than in its outskirts. Similarly, Rodríguez-Medrano et al. (2023) find that in the last 2 Gyr, the fraction of accreted mass from high-mass halos in voids has declined relative to low-mass halos, suggesting that more massive systems in these regions have reached a dynamically quieter state. Our results corroborate these findings and extend them by elucidating how such formation histories specifically shape the HOD within void regions. Building on this framework, we proceeded to investigate the spatial arrangement and evolutionary pathways of halo densities in void environments.

To quantify the local density, we measured the distance from each central galaxy to its fifth nearest neighboring galaxy brighter than M_r − 5log₁₀(h) = − 20, regardless of whether the neighbor is a central or satellite galaxy, and normalized this quantity by the average fifth nearest bright neighboring of the complete snapshot. This parameter, referred as D₅, provides a reliable measure of the local galaxy density surrounding each halo. A value of D₅ < 1 indicates a high local density, while the values > 1 suggest a more isolated environment. The choice of M_r − 5log₁₀(h) = − 20 as a brightness threshold in each snapshot is consistent with our void identification methodology, ensuring the use of the same tracers of the distribution of the underlying field matter for environmental characterization. Table 2 shows the physical scale corresponding to D₅ = 1 in each snapshot, which typically represents environments with densities comparable to the mean cosmic density. As expected, the scales of D₅ = 1 tend to increase as we approach z = 0. Figure 7 illustrates the spatial distribution of central galaxies at z = 0, showing a slice of the simulation box with a thickness of 50 h⁻¹ Mpc. Each central galaxy is represented by a point, colored according to its D₅, providing an intuitive visualization of density variations throughout the large-scale structure.

Fig. 7. Spatial distribution of central galaxies at z = 0 within a simulation box slice spanning 100 Mpc h−1 < z < 150 Mpc h−1. The color of each point represents D5, the distance to the fifth nearest bright neighbor, normalized by the mean value of the complete snapshot, which serves as a proxy for the local galaxy density. Left panel: Area of 1000 Mpc h−1 per side. Right panel: Zoomed-in view of the highlighted region.

In addition to these parameters, we analyzed the spatial distribution of these central galaxies at each redshift. For each galaxy, we computed the distance (r) from the center of the void it will inhabit at z = 0. In this way, r/R_void serves as an indicator of the galaxy’s location within the region that will eventually evolve into a void at z = 0. Thus, the combined analysis of r/R_void and D₅ allows us to determine whether these galaxies have always resided in the inner regions of the area that will evolve into a void at z = 0. In addition, it enables us to quantify how much the local density surrounding them has changed as this region evolved.

Figure 8 visualizes these results. The left panel shows the relationship between the halo mass of central galaxies in voids at z = 0 and their normalized distance from the center of the void for five of the redshifts studied, which decreases as the redshift increases. The right panel illustrates how D₅ varies with the mass of these same halos.

Fig. 8. Left panel: Normalized distance from the void center of central galaxies in voids at z = 0 as a function of their halo mass for five redshift samples. Circle sizes decrease with increasing redshift, and the color scale represents the mean D5 for each mass bin. Right panel: Variation in D5 with halo mass for the same samples. The color scale indicates the mean r/Rvoid for each bin. The error lines were estimated via the Jackknife resampling method and using 50 subsamples.

The combined analysis of the two panels reveals how the halos in each mass bin evolve within the region that will form the void at z = 0. We observe a clear trend indicating a correlation between halo mass and their location within the void. Lower-mass halos tend to populate the innermost regions, while more massive halos are found closer to the void walls. This mass segregation within voids is consistent with previous findings by Gottlober et al. (2003) and Ricciardelli et al. (2014, 2017) for galaxies. Additionally, as expected, the local density of halos increases as we move toward the outer regions of the void, leading to a direct correlation between the halo mass and its local density.

Examining the evolution of these results with redshift, we find that, in general, galaxies that populate a void at z = 0 have always resided within the volume that defines these regions. The overall trend shows that these galaxies migrate from more central areas toward the outer parts of the void. This behavior aligns with the well-established dynamics of cosmic voids. Previous studies have demonstrated that voids behave similarly to under-dense regions in an expanding Universe, where the inner shells expand faster than the outer ones, creating a distinct velocity field characterized by coherent outflows (Dubinski et al. 1993; Padilla et al. 2005; Ceccarelli et al. 2006). This differential expansion causes galaxies within voids to gradually migrate toward the boundaries over cosmic time, consistent with our findings. At the same time, we see that the evolution of local density differs depending on the halo mass range, following a pattern similar to that observed in the HOD. Halos with M_halo ≲ 10¹² h⁻¹ M_⊙ show no significant variations in D₅. However, for halos with M_halo ≳ 10¹² h⁻¹ M_⊙, it becomes evident that the local density decreases as the redshift approaches zero. This effect is increasingly pronounced for more massive halos: halos with M_halo ≳ 10¹³ h⁻¹ M_⊙ transition from residing in locally overdense regions (D₅ < 1) to being surrounded by locally under-dense environments (D₅ > 1). In particular, the local density of halos in the highest mass bin remains above 1 even at z = 0, reinforcing the idea that these halos lie near the void boundaries, at the interface with higher-density structures such as sheets or filaments.

These findings are consistent with a scenario in which the voids are populated by halos with distinct formation histories driven by their halo masses, as we describe in Fig. 6. Low-mass halos are younger than average and predominantly occupy the inner regions of voids. Most of these halos have inhabited void interiors since high redshifts, evolving in sub-dense regions. In contrast, more massive void halos tend to be located on void walls and are older than average. Most of these halos transition from regions of high local density to sub-dense and average environments, potentially impacting their development. These differences in void halo populations may explain the variations detected in the HOD within these regions.

5. Summary and conclusions

We have investigated the evolution of the HOD within cosmic voids from redshift z = 0.5 to the present day, utilizing the MDPL2-SAG catalog and the SPARKLING void finder algorithm. We aimed to understand how the under-dense void environment influences the HOD over cosmic time and explored the underlying mechanisms responsible for the distinctive galaxy populations observed in these regions.

Our analysis of the HOD across different redshifts revealed several key findings:

• Persistently low HODs in voids: We confirm that the HOD within cosmic voids is consistently lower than the global HOD across the redshift range studied (z = 0 to z = 0.49). This is evident for halos with masses greater than ∼10¹² h⁻¹ M_⊙. This indicates that the void environment consistently hinders galaxy occupation in halos, especially in more massive systems, throughout cosmic history.

• Stable relative HOD for z = 0 void galaxies: When tracing the HOD evolution of galaxies that reside in voids at z = 0 (and were central galaxies over time), we found that their HOD remains consistently lower than the global HOD at all redshifts. Importantly, the ratio of their HOD to the global HOD does not show significant redshift evolution. This implies that the factors that determined the lower HOD of these specific void galaxies were already in place at higher redshifts and maintained a consistent relative effect over time.

• Mass-dependent environment in voids: The spatial distribution of void halos at z = 0 exhibits a clear mass dependence. Lower-mass halos tend to reside in the inner regions of voids, whereas more massive halos are preferentially found near the void boundaries, in agreement with previous studies (Gottlober et al. 2003; Ricciardelli et al. 2014, 2017). Moreover, the local density around massive void halos decreases significantly from high redshifts to z = 0, indicating a transition from initially denser environments to under-dense void regions. In contrast, lower-mass halos have consistently resided in under-dense environments since early times.

• Mass-dependent halo formation: The particular environmental conditions within voids create a relatively quiescent setting for halos, where mergers and interactions with other structures are less frequent compared to higher-density regions. Low-mass halos, which predominantly inhabit locally under-dense regions, may have experienced inhibited matter accretion, delaying their assembly compared to halos of the same mass in an average-density environment. Conversely, massive halos have transitioned from initially high-density regions to medium- or low-density void environments, which may have facilitated their earlier assembly compared to similar halos in other cosmic environments.

These results suggest a scenario in which the void environment exerts a complex, mass-dependent influence on galaxy formation and halo occupation. The persistently lower HOD in voids across redshifts suggests that fundamental environmental mechanisms that suppress galaxy formation within these under-dense regions are at play. The distinct formation histories and environmental evolution of low- and high-mass halos in voids likely contribute to the observed mass-dependent HOD trends. The mass scale of ∼10¹² h⁻¹ M_⊙ emerges as a natural threshold for the global void effects on the HOD, the specific lowering of the HOD for z = 0 void galaxies, the formation time of void halos, and the local density evolution of z = 0 void galaxies. Below this mass, the void environment appears to have a minimal impact on all these effects, suggesting that galaxy formation in these lower-mass halos is governed by processes largely independent of the void environment itself. Above this mass, the void environment becomes a significant factor, suppressing galaxy occupation, particularly in more massive halos. In such haloes, their location on void walls and their transition from denser to under-dense environments might influence their accretion history and satellite galaxy capture, leading to observed HOD suppression, particularly at later times. The distinct evolutionary paths oflow- and high-mass halos within voids, in terms of formation time and environmental history, likely underpin these mass-dependent HOD variations.

Future work with hydrodynamical simulations is crucial to disentangling the specific physical mechanisms driving the mass-dependent impact of voids on the HOD, including gas accretion, feedback, and mergers. They should also explore broader environmental factors, such as tidal forces and merger rates, and incorporate analysis of galaxy properties and their evolution across cosmic time to fully understand the interplay of environment, halo formation history, and galaxy occupation in these under-dense regions.

The results in this work are in conflict with observations. Therefore, it would be interesting to develop methodologies that allow us to estimate the HOD precisely at different redshifts using spectroscopic surveys or with the photometric redshifts that will be available in the coming years.

Acknowledgments

The authors thank Dra. Yamila Yaryura for her help with the MDPL2-SAG data. The COSMOSIM database used in this paper is a service of the Leibniz Institute for Astrophysics Potsdam (AIP). The MULTIDARK database was developed in cooperation with the Spanish MultiDark Consolider Project CSD2009-00064. FR acknowledges support from the ICTP through the Junior Associates Programme 2023-2028. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at the Leibniz Supercomputing Centre (LRZ, www.lrz.de). This work has been partially supported with grants from Agencia Nacional de Promoción Científica y Tecnológica, the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina) and the Secretaría de Ciencia y Tecnología de la Universidad Nacional de Córdoba (SeCyT-UNC, Argentina). Figures were developed using MATPLOTLIB (Hunter 2007) and some of them were post-processed with INKSCAPE (https://inkscape.org).

References

All Tables

All Figures

	Fig. 1. z = 0 (black lines) and z = 0.49 (red lines) void size distribution. The dashed vertical lines mark the radius up to which the sample can be considered complete.
In the text

Fig. 2. Top panels: HOD measured for galaxies with Mr − 5log10(h) <  − 17 at redshifts z = 0, z = 0.3, and z = 0.5 for the MDPL2-SAG simulation. The blue lines represents the HOD of the complete catalogs, while the black and red lines show the HOD of halos inside the cosmic voids identified in each snapshot. Bottom panels: Ratio of the HOD inside voids to the HOD in the complete catalogs at each redshift. The error bars were estimated via the Jackknife resampling method and using 50 subsamples.

In the text

	Fig. 3. Evolution of the HOD of central galaxies in voids at z = 0. Top panel: HOD for central galaxies in voids at z = 0 across different redshifts. Bottom panel: Ratio between the HOD at each redshift and the HOD at z = 0. The error bars were estimated via the Jackknife resampling method and using 50 subsamples.
In the text

Fig. 4. Ratio between the HOD of central galaxies located in voids at z = 0 and the general HOD of the Universe at three different redshifts. The solid black line with diamonds corresponds to z = 0, the dashed red line with squares to z = 0.304, and the dashed pink line with triangles to z = 0.49. The HOD for all cases was calculated for galaxies with Mr − 5log10(h) <  − 17 at z = 0 with halo masses fixed to their z = 0 values. The error bars were estimated via the Jackknife resampling method and using 50 subsamples.

In the text

	Fig. 5. Variation in the HOD within voids and in the complete snapshot between each redshift compared to the HOD at z = 0. The errors were estimated via the Jackknife resampling method and using 50 subsamples.
In the text

Fig. 6. Halo formation redshifts as a function of the halo mass. Top panel: Results for zform, defined as the redshift when half of the maximum halo mass was accreted. Bottom panel: Results for zlastmm, the redshift at which the last major merger occurred. The blue lines represent the complete z = 0 snapshot and the black line the z = 0 void halo samples. The error bars were estimated via the Jackknife resampling method and using 50 subsamples.

In the text

Fig. 7. Spatial distribution of central galaxies at z = 0 within a simulation box slice spanning 100 Mpc h−1 < z < 150 Mpc h−1. The color of each point represents D5, the distance to the fifth nearest bright neighbor, normalized by the mean value of the complete snapshot, which serves as a proxy for the local galaxy density. Left panel: Area of 1000 Mpc h−1 per side. Right panel: Zoomed-in view of the highlighted region.

In the text

Fig. 8. Left panel: Normalized distance from the void center of central galaxies in voids at z = 0 as a function of their halo mass for five redshift samples. Circle sizes decrease with increasing redshift, and the color scale represents the mean D5 for each mass bin. Right panel: Variation in D5 with halo mass for the same samples. The color scale indicates the mean r/Rvoid for each bin. The error lines were estimated via the Jackknife resampling method and using 50 subsamples.

In the text