We employ the HOD approach suitably calibrated to the [OII] emission-line galaxies (ELGs) observed by DESI (Rocher et al., 2024). We assume a Gaussian HOD model for the central galaxies with the expected occupation 〈 N cent M 〉 = A c 2 π σ M ⋅ exp − log 10 M − log 10 M c 2 2 σ M 2 ≡ 〈 N cent GHOD M 〉 and a power law model for the satellite galaxies with ⟨ N sat M ⟩ = A s M − M 0 M 1 α . Here, A_c and A_s set the number density to match those of the observed galaxies; M_c is the typical halo mass to host the [OII] emitter as a central galaxy; σ M is the width of the Gaussian distribution. The low-mass cutoff M₀ determines whether a halo possesses a satellite galaxy, and α controls the halo mass dependence of the number of satellite galaxies. In this model, M₁ was introduced for normalization and is fixed at 1 0 13 h − 1 M ⊙ . The central galaxies are populated such that their positions and velocities match those of the halo centers. The satellite galaxies are distributed within a halo with a Navarro–Frenk and White profile (Navarro et al., 1996; Navarro et al., 1997) ρ r = ρ s r r s 1 + r r s 2 , where r is the distance from the halo center and r _s is the scale radius. The HOD parameters used in the present study are listed in Table 2.

We first generated 3D LIM mock catalogs using the COLA halo catalogs at redshift 1.0. It is expected that the star-forming galaxies at this redshift have strong CO(3-2) emissions as excellent targets for CONCERTO. We also included the interloper of CO(4-3) at redshift 1.67, which is redshifted to the same wavelength. Here, we briefly describe the method for assigning the CO(3-2) and CO(4-3) luminosities to each halo. First, we calculated the stellar mass of a halo using a double power-law stellar-to-halo mass ratio (SHMR) model proposed by Moster et al. (2010) as follows: M ∗ M h z = 2 A z M h M A z − β z + M h M A z γ z − 1 .

We then use the best-fit parameters derived in Girelli et al. (2020) (A = 0.0353, log(M _A) = 12.05, β = 0.88, and γ = 0.599) with a 0.2 dex scatter. Next, we use an empirical relationship between the stellar mass and star formation rate (SFR) to calculate the infrared luminosity (Béthermin et al., 2022); in this procedure, we split the galaxies into star-forming and passive galaxies using the star-forming fraction given by Equation 1: f SF M ∗ , z = 1 − f Q , 0 z 1 − erf log 10 M ∗ − log 10 M t z σ SF z 2 (1) where f Q , 0 ( z ) , M t ( z ) , and σ SF ( z ) denote the passive fractions of low-mass galaxies ( M ∗ ≪ M t ( z ) ) , stellar mass threshold separating the star-forming and passive galaxies, and width of the transition as a function of the redshift.

The SFR of a galaxy is calculated using Equation 2: log 10 SFR M ⊙ / yr = log 10 M ∗ 1 0 9 M ⊙ − m 0 + a 0 ⁡ log 10 z − a 1 max 0 , log 10 M ∗ 1 0 9 M ⊙ − m 1 − a 2 ⁡ log 10 1 + z 2 . (2)

The infrared luminosity L _IR is calculated from the SFR using the Kennicutt (1998) conversion factor of 1.0 × 1 0 − 10 M ⊙ / yr / L ⊙ . The CO(1-0) luminosity is then calculated from the L _IR using the empirical relationship of Equation 3 derived in Sargent et al. (2014): log 10 L CO 1 − 0 ′ K   k m   s − 1 p c 2 = 0.81 ⁡ log 10 L IR L ⊙ + 0.54 (3)

The other transitions are calculated using the spectral line energy distribution template suggested by Bournaud et al. (2015). For the CO(4-3) interlopers, we generated a halo catalog at redshift 1.67 and calculated the CO(4-3) line luminosity using the steps described above. Since we expect that the CO(4-3) contaminants at z = 1.67 are not spatially correlated with the distribution of the CO(3-2) emitters at z = 1.0, we used a different random seed to generate the initial condition for the z = 1.67 halo catalog. We then downselect 20% of the halos by assuming that 80% of the interlopers can be removed using the methods discussed in Section 5. Although the contamination fraction is uncertain and arbitrary, we choose this value to quantify the impacts of the interlopers on void statistics.

We list the basic observational parameters used to generate the LIM catalog in Table 3; here, the parameters were set to closely follow the specifics of CONCERTO. We conducted the following analysis on a single slice of the LIM mock corresponding to an observer frame frequency of 173 GHz. Since we set the spectral resolution to 1 GHz and angular resolution to 32 arcsec, a single pixel in the 2D image included emissions from cubes of approximately 0.36 h^–1 Mpc width and 19.5 h^–1 Mpc depth for CO(3-2) emitters as well as 0.50 h^–1 Mpc width and 18.1 h^–1 Mpc depth for CO(4-3) emitters in the comoving space. We applied a minimum intensity threshold of 1 kJy/sr (corresponding to 0.52 K in brightness temperature) to remove low-intensity pixels with low signal-to-noise ratios in the actual observations and treated the pixels above this threshold as “signals” in the 2D maps.

We located voids in the galaxy catalogs using a popular void finder called VIDE, whose full technical description is provided in Sutter et al. (2014); its applications in cosmology include cosmological parameter estimations using SDSS galaxies (Contarini et al., 2023) and constraints on neutrino masses (Bayer et al., 2024). Briefly, VIDE generates the Voronoi tessellation for a given particle (galaxy) distribution, where the density of each Voronoi cell is calculated as the inverse of its volume (or area in the 2D case). Cells that flow toward the same local density minimum are grouped together to form a zone; for each zone, the zone center is calculated as the volume-weighted center of all the particles in the zone, as shown in Equation 4: r center = ∑ i r i V i ∑ i V i . (4)

The initial radius R _ini of each zone is calculated as the radius of a sphere (or circle in the 2D case) using the corresponding volume (area) of the detected zone, as shown in Equation 5: 3 4 π R ini 3 = ∑ i V i , π R ini 2 = ∑ i V i . (5)

The central density of each zone is then calculated as the mean density within a spherical (circular) region around its center with a radius of 0.25R _ini. For direct comparisons with theoretical or analytic models, the zones are cleaned and rescaled as described in Jennings et al. (2013). The procedures are as follows:

• Zones with central densities greater than ( 1 + δ ν NL ) ρ ̄ will be removed from the void candidates. Here, ρ ̄ is the mean particle density and δ ν NL is the non-linear density threshold for void detection. In our study, we used δ ν NL = − 0.7 but note that the threshold can be set arbitrarily; however, the same threshold must be set for the theoretical prediction.

The void catalogs thus generated were used in the study subsequently.

One of the fundamental properties of a void is its size. Analogous to the halo mass function that characterizes the abundance of halos with different masses, the void size function describes the abundance of voids with different sizes and serves as a crucial statistical measure of the voids present in a large volume. A theoretical prediction of the void size function was provided by Sheth and van de Weygaert (2004) using the excursion-set approach (Press and Schechter, 1974), according to which the void size function in the linear regime is expressed by Equation 6: d n d ⁡ ln ⁡ r lin = f ln ⁡ σ σ V r d ⁡ ln σ − 1 d ⁡ ln ⁡ r (6) where f ln ⁡ σ corresponds to the total area occupied by the voids, as defined by Equation 7: f ln ⁡ σ σ = 2 ∑ j = 1 ∞ j π x 2 ⁡ sin j π D exp − j π x 2 2 (7) Here, σ is the square root of the mass variance, and x and D are defined by Equations 8, 9, respectively: x = D | δ ν L | (8) D = | δ ν L | δ c L + | δ ν L | , (9) Here, δ ν L is the under-density threshold to form voids, and δ c L is the over-density threshold to form halos, according to the excursion-set theory. The density threshold evolves with redshift according to Equation 10 : δ L z = δ L 0 D z D 0 , (10) where D ( z ) is the linear growth factor.

To predict the size function in a non-linear regime, we assume a volume-conserving model in which the total volume occupied by the voids is conserved within a single void-size bin. Under this assumption, the void size function can be rewritten as in Equation 11 : d n d ⁡ ln ⁡ r = d n d ⁡ ln ⁡ r lin V r L V r d ⁡ ln r L d ⁡ ln ⁡ r , (11) where r L is the void radius defined in linear theory. The non-linear density threshold can be converted to its linear equivalent (Equation 12 ) using the fitting function derived in Bernardeau (1994). Thus, δ ν L = C 1 − 1 + δ ν NL − 1 / C (12) with C = 1.594 .

For biased tracers, the density threshold will also be biased, as shown in Equation 13 : δ ν , tr NL = b ν δ ν NL , (13) where b ν corresponds to the tracer bias in the under-dense region. Its value can be directly calculated from the linear bias b eff using the empirical relation (Equation 14 ) derived in Contarini et al. (2019): b ν = 0.854 ± 0.007 b eff + 0.420 ± 0.010 . (14)

Since we set our threshold as δ ν , tr NL = − 0.7 , the linear density threshold used for the linear prediction is given by Equation 15 : δ ν L = 1.594 1 − 0 . 3 − 1 / 1.594 / b ν (15)

Since the formation and evolution of voids are governed by the nature of dark energy and dark matter, the abundance of voids with different sizes is sensitive to the cosmological parameters (Li et al., 2012; Verza et al., 2019). Herein, we explore the size function of a void using different tracers and their dependence on the cosmological differences.

The radial profile of the voids is another important characteristic of the cosmic structure, which may be sensitive to clustering scales smaller than the free-streaming lengths of dark matter and relic neutrinos. In our study, we defined the radial profile of the voids in terms of particle density as a function of distance from the void center; here, the void center was defined as the centroid of the void region.