Signature of Massive Neutrinos from the Clustering of Critical Points. I. Density-threshold-based Analysis in Configuration Space

For this first work, we adopt a technique centered on the extraction and classification of critical points from smoothed density fields. The construction of such fields is obtained via a grid interpolation strategy, directly from the DM particle distribution of N-body simulations. In essence, in the first step DM particles are turned into cells of given density, which are properly smoothed up to a desired scale. Critical points are then identified and classified according to their type. Subsequently, those points are sorted in terms of density threshold (or rarity), and for a fixed rarity their clustering statistics is characterized as a function of M_ν and redshift. The advantage of this three-step procedure resides in its simplicity, ease of implementation, good efficiency, and effective numerical performance. In what follows, we briefly describe the three main building blocks of our pipeline.

We begin with the production of density fields, which are constructed from the DM particle output of selected QUIJOTE N-body realizations using the Python libraries for the analysis of numerical simulations (Pylians; Villaescusa-Navarro 2018). ⁷ Pylians are a set of Python, Cython, and C libraries that are helpful in the analysis of large-volume runs, providing a number of statistical tools—as well as routines to create and manipulate density fields. Specifically, at a fixed redshift, we build 3D density fields ρ( x ) from the spatial positions ( x ) of N-body particles using a triangular-shape cloud (TSC) mass-assignment scheme to deposit particle masses into the grid; in this process, no weights are associated with each particle. Density contrast fields, defined as $\delta ({\boldsymbol{x}})=\rho ({\boldsymbol{x}})/\bar{\rho }-1$, with $\bar{\rho }$ the corresponding average background density, are readily inferred from the global knowledge of ρ( x ). We adopt N_grid = 1024 voxels per dimension on a regular grid spacing—corresponding to 2N _p , with N _p = 512 the number of particles per dimension in a given QUIJOTE simulation snapshot. We also tested different mass-assignment schemes, as well as a number of voxel choices for N_grid, and found negligible dependence on our results.

In the next step, fields are subsequently smoothed before the extraction and classification of critical points. More explicitly, DM density fields created from QUIJOTE snapshots are smoothed with a Gaussian kernel W_G over R_G = 3 h⁻¹ Mpc scales, using Py-Extrema. ⁸ Smoothing is relevant for this density-threshold-based methodology: we provide more extensive details related to our choice of the smoothing scale in Appendix A. We will also return to this important point in a companion publication, where we confront the present technique with a (fundamentally different) persistence-based approach—which instead does not rely directly on smoothing. In essence, the effect of smoothing is analogous to a resolution cutoff, erasing small-scale structures: at a given z-interval, structures of different sizes and masses are present, and therefore selecting an explicit R_G would correspond to averaging over specific density regions. A suitable smoothing choice would then be relevant if we were to identify critical points with actual physical structures, as the effective typical size of such objects depends on the adopted smoothing scale of the density field: to this end, see the considerations presented in Section 4.1.

Once density fields have been properly smoothed, critical points are identified with the same detection algorithm adopted in Cadiou et al. (2020) and in Shim et al. (2021), based on a local quadratic estimator centered on a second-order Taylor expansion of the density field about critical points (see also Gay et al. 2012 for additional details). Such expansion yields

$$\begin{eqnarray}&&{\boldsymbol{x}}-{{\boldsymbol{x}}}_{c}\sim {({\rm{\nabla }}{\rm{\nabla }}{\rho })}^{-1}{\rm{\nabla }}{\rho },\end{eqnarray} \tag{ 1 }$$

with x _c the spatial position of the critical point. In essence, for each grid cell the gradient and Hessian of the density field are computed, Equation (1) is solved, and solutions found at a distance greater than 1 pixel are discarded via the requirement ${\max }_{i}(| {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{{\rm{c}}}| )\lt 1\,\mathrm{pixel}$, with i = 1, 2, 3—only the critical point closest to the center of the cell is retained. Subsequently, each cell enclosing a critical point is flagged, and an efficient loop is performed over flagged cells containing multiple critical points of the same type. In this process, critical points are classified according to their corresponding rank order. Specifically, the classification is based on the number of negative eigenvalues of the density Hessian—i.e., 0 for voids (${ \mathcal V }$ or minima), 1 for walls (${ \mathcal W }$ or "saddle 1" or wall-type saddles), 2 for filaments (${ \mathcal F }$ or "saddle 2" or filament-type saddles), and 3 for peaks (${ \mathcal P }$ or maxima). ⁹

At this stage of the pipeline, critical points have been properly identified and classified. Next, all of the points are sorted according to the chosen density threshold ν = δ/σ, with δ the overdensity previously defined (but now referred to the smoothed density field) and σ the corresponding rms fluctuation of the field. Following Shim et al. (2021), we adopt an identical ad hoc convention in defining rarity levels, namely, we always sample the population that provides the same abundance for a given type of critical point. This choice represents also a relevant point that should be kept in mind later on, as it has some impact on the interpretation of the results presented in Section 4 (i.e., the adopted rarity definition allows one to avoid density overlapping among critical points, thus enhancing more remarkable features in correlations).

Our rarity level definition can be formalized as follows. Let us first introduce the notation $\alpha \in \{0,1,2,3\}\equiv \{{ \mathcal V },{ \mathcal W },{ \mathcal F },{ \mathcal P }\}$ to specify the critical point type. We will also use the indices ℓ and k in place of α, whenever necessary (i.e., α ≡ ℓ, k). In essence, for peak- and filament-type (respectively, void- and wall-type) critical points, we select ensembles with ν higher (respectively, lower) than a given threshold ${\nu }_{\alpha ,{ \mathcal R }}$, where ${ \mathcal R }$ indicates the relative abundance cut. ¹⁰ The threshold is then fixed as the rarity for each type of critical point α yielding the same relative abundances defined by the ratios

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{{\rm{C}},\alpha }(\nu \geqslant {\nu }_{\alpha ,{ \mathcal R }})}{{n}_{{\rm{C}},\alpha }},\end{eqnarray} \tag{ 2 }$$

with $\alpha \in \{{ \mathcal P },{ \mathcal F }\}\equiv \{3,2\}$ for peaks and filaments, respectively, and

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{{\rm{C}},\alpha }(\nu \leqslant {\nu }_{\alpha ,{ \mathcal R }})}{{n}_{{\rm{C}},\alpha }},\end{eqnarray} \tag{ 3 }$$

with $\alpha \in \{{ \mathcal W },{ \mathcal V }\}\equiv \{1,0\}$ for walls and voids, respectively. In the previous expressions, n_C,α indicates the entire ensemble of critical points of type α, while ${n}_{{\rm{C}},\alpha }(\nu \geqslant {\nu }_{\alpha ,{ \mathcal R }})$ and ${n}_{{\rm{C}},\alpha }(\nu \leqslant {\nu }_{\alpha ,{ \mathcal R }})$ represent the subset containing critical points of type α above and below the selected threshold ${\nu }_{\alpha ,{ \mathcal R }}$, respectively. For ease of notation, in what follows we define $\gamma \equiv {\nu }_{\alpha ,{ \mathcal R }}$ and also indicate n_C,α (ν ≥ γ) or n_C,α (ν ≤ γ) generically with ${n}_{{\rm{C}},\alpha }^{\gamma }$. In our study we consider three cuts in rarity, namely ${ \mathcal R }\in \{5,10,20\}$—where the abundances are expressed in percentages.

We then compute the spatial clustering statistics of critical points above (below) the threshold—i.e., real space two-point auto- and cross-correlations ξ(r) as a function of the separation r—using the Halotools ¹¹ platform (Hearin et al. 2017), equipped with efficient algorithms for calculating clustering statistics (including cross-correlations). For our two-point calculations in configuration space, we adopt the widely used minimal-variance Landy−Szalay (LS) estimator (Landy & Szalay 1993) assuming periodic boundary conditions and utilize random catalogs always at least 20 times bigger in size than the considered data. ¹² Specifically, for a given catalog ${C}_{\alpha }^{\gamma }$ of critical points above (below) a selected density threshold γ having total size ${n}_{{\rm{C}},\alpha }^{\gamma }$ and assuming a corresponding random catalog ${R}_{\alpha }^{\gamma }$ of total size ${n}_{{\rm{R}},\alpha }^{\gamma }$ characterized by a random uniform probability distribution of points within the same volume, the LS estimator for the two-point correlation function ${{\xi }}_{k{\ell }}^{{\gamma },{\rm{L}}{\rm{S}}}$ above (below) γ reads as follows: ¹³

$$\begin{eqnarray}&&{{\xi }}_{k{\ell }}^{{\gamma },{\rm{L}}{\rm{S}}}=\sqrt{{\eta }{\beta }}\displaystyle \frac{{\unicode{x02329}}{C}_{k}^{{\gamma }}{C}_{{\ell }}^{{\gamma }}{\unicode{x0232A}}}{{\unicode{x02329}}{R}_{k}^{{\gamma }}{R}_{{\ell }}^{{\gamma }}{\unicode{x0232A}}}-\sqrt{{\eta }}\displaystyle \frac{\sqrt{{\unicode{x02329}}{C}_{k}^{{\gamma }}{R}_{{\ell }}^{{\gamma }}{\unicode{x0232A}}{\unicode{x02329}}{C}_{{\ell }}^{{\gamma }}{R}_{k}^{{\gamma }}{\unicode{x0232A}}}}{{\unicode{x02329}}{R}_{k}^{{\gamma }}{R}_{{\ell }}^{{\gamma }}{\unicode{x0232A}}}+1,\end{eqnarray} \tag{ 4 }$$

where the normalization factors are specified by

$$\begin{eqnarray}&&{\eta }=\displaystyle \frac{({n}_{{\rm{R}},k}^{{\gamma }}-1)({n}_{{\rm{R}},{\ell }}^{{\gamma }}-1)}{{n}_{{\rm{C}},k}^{{\gamma }}{n}_{{\rm{C}},{\ell }}^{{\gamma }}}\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{\beta }=\displaystyle \frac{{n}_{{\rm{R}},k}^{{\gamma }}{n}_{{\rm{R}},{\ell }}^{{\gamma }}}{({n}_{{\rm{C}},k}^{{\gamma }}-1)({n}_{{\rm{C}},{\ell }}^{{\gamma }}-1)}.\end{eqnarray} \tag{ 6 }$$

Note that all of the correlation measurements obtained via Equation (4) are reported as a function of the spatial separation r, conventionally expressed in units of h⁻¹ Mpc. Unless otherwise specified, we always use a bin size of 5 h⁻¹ Mpc: this choice is motivated by the study presented in Appendix A, where we also discuss smoothing and bin size effects on clustering measurements. Moreover, in order to accurately determine the spatial locations of the inflection points of the correlation functions at large scales (see Section 4.2), we also employ a "refinement" technique where we adopt a smaller bin size of 1 h⁻¹ Mpc in selected regions near such points.

We begin with an instructive visualization of the geometrical distribution of critical points identified and classified via the "density-threshold-based" procedure described in Section 3. Figure 2 shows the spatial location of critical points in QUIJOTE-simulated patches at z = 0, color-coded by type, for three different rarity thresholds. The squared patches in the figure are characterized by a side of 100 h⁻¹ Mpc in the x–y plane and a 50 h⁻¹ Mpc depth along the z-coordinate. The underlying normalized density fields are the results of sampling the density fields with Pylians assuming N_grid = 1024. Top panels refer to the baseline cosmology (massless neutrino scenario), middle panels are for a massive neutrino cosmology with M_ν = 0.1 eV, and bottom panels display results for M_ν = 0.4 eV. From left to right, the rarity threshold is increased in terms of relative abundance cut, corresponding to 5%, 10%, and 20%, respectively, via the selection criteria detailed in Section 3.4.

Zoom In Zoom Out Reset image size

As anticipated before (i.e., Section 3.3), care must be taken in directly identifying critical points derived from a "density-threshold-based" approach with cosmological structures. This is primarily because such methodology relies on smoothing, and therefore the effective typical size of the LSS cosmic web constituents depends on the smoothing scale of the density field, which acts as a resolution cutoff. To this end, see also the discussion in Shim et al. (2021), who adopted a smoothing scale of R_G = 6 h⁻¹ Mpc corresponding to ∼10¹⁵ M_⊙ average density regions. As noted by the same authors, with such a choice of R_G only large voids can be resolved, and only 5% of all density peaks represent virialized galaxy clusters at z = 0, while the rest is still in a collapse process. From Figure 2, one can deduce that decreasing the rarity threshold (namely, moving from the left to the right panels in the figure) provides a more detailed mapping of the corresponding underlying cosmological structures. This fact can be inferred from correlation function measurements: larger values of ${ \mathcal R }$ manifest into smaller dispersions/scatter in ξ(r), as we show in Sections 4.2 and 4.3. Moreover, it is also interesting to compare (even just at the visual level) the spatial distribution of critical points for a fixed rarity value with increasing neutrino mass: differences in the abundance and spatial location of critical points when M_ν ≠ 0.0 eV translate into morphological and topological differences that carry rich cosmological information.

In this regard, Table 2 reports the overall abundance of the entire set of critical points at z = 0, classified by type, in the three cosmological frameworks considered in this work. These measurements are the results of running Py-Extrema on the smoothed density fields and represent averages over 100 independent QUIJOTE realizations at fixed cosmology. Note that for a Gaussian field it is expected that the numbers of peaks (${ \mathcal P }$) and voids (${ \mathcal V }$) are the same, as well as the numbers of filaments (${ \mathcal F }$) and walls (${ \mathcal W }$). Moreover, in Gaussian fields, the ratio between filaments and peaks (${ \mathcal F }/{ \mathcal P }$) or between walls and voids (${ \mathcal W }/{ \mathcal V }$) is estimated to be ∼3.05, and also the total number of extrema and saddle points is preserved at the first non-Gaussian perturbative order. In addition, for sufficiently large volumes, the ratio between the number of peaks and walls over voids and filaments ($[{ \mathcal P }/{ \mathcal W }]/[{ \mathcal V }/{ \mathcal F }]$) in a Gaussian field is very close to unity throughout the entire z-evolution. This is because the genus topology, equal to the alternate sum of critical points, should be preserved, namely,

$$\begin{eqnarray}&&{n}_{\mathrm{cp}}^{{ \mathcal P }}-{n}_{\mathrm{cp}}^{{ \mathcal F }}+{n}_{\mathrm{cp}}^{{ \mathcal W }}-{n}_{\mathrm{cp}}^{{ \mathcal V }}=0\,,\end{eqnarray} \tag{ 7 }$$

with ${n}_{\mathrm{cp}}^{(i)}$ the mean number density of peaks, filaments, walls, and voids—respectively.

Table 2. Total Number of Critical Points at z = 0, Classified by Type, for the Three Cosmologies Considered in This Work

	Mν (eV)	z = 0
	0.0	228,094 ± 382
Minima (${ \mathcal V }$)	0.1	227,770 ± 355
	0.4	226,527 ± 402
	0.0	680,088 ± 823
Saddle 1 (${ \mathcal W }$)	0.1	678,915 ± 775
	0.4	674,157 ± 797
	0.0	669,208 ± 858
Saddle 2 (${ \mathcal F }$)	0.1	667,952 ± 826
	0.4	662,960 ± 766
	0.0	223,951 ± 337
Maxima (${ \mathcal P }$)	0.1	223,571 ± 319
	0.4	222,031 ± 325

Download table as:  ASCIITypeset image

However, nonlinear evolution breaks the symmetry between underdense and overdense regions. Table 3 summarizes all of the previous considerations, showing clear departures from Gaussianity—expected because of nonlinear structure formation, in addition to the presence of massive neutrinos. Note that the effect of a nonzero neutrino mass is represented by an overall decrement in the abundance of critical points, much more pronounced for higher neutrino masses. For peaks, this fact has repercussions at the level of the halo mass function, and it can be interpreted in the framework of the halo model (see, e.g., Rossi 2017, their Section 4.3).

Table 3. Relevant Abundance Ratios of Critical Points at z = 0 for the Three Cosmologies Adopted in This Study

	Mν (eV)	z = 0
	0.0	0.9818 ± 0.0022
${ \mathcal P }/{ \mathcal V }$	0.1	0.9816 ± 0.0021
	0.4	0.9802 ± 0.0023
	0.0	0.9840 ± 0.0017
${ \mathcal F }/{ \mathcal W }$	0.1	0.9839 ± 0.0017
	0.4	0.9834 ± 0.0016
	0.0	2.9888 ± 0.0059
${ \mathcal F }/{ \mathcal P }$	0.1	2.9877 ± 0.0056
	0.4	2.9859 ± 0.0056
	0.0	2.9816 ± 0.0062
${ \mathcal W }/{ \mathcal V }$	0.1	2.9807 ± 0.0058
	0.4	2.9761 ± 0.0063
	0.0	0.9661 ± 0.0028
(${ \mathcal P }/{ \mathcal W }$)/(${ \mathcal V }/{ \mathcal F }$)	0.1	0.9657 ± 0.0026
	0.4	0.9639 ± 0.0027

Note. Here we highlight the fact that nonlinear evolution and the effects of massive neutrinos (if present) break the symmetry between overdense and underdense regions—expected instead in a purely Gaussian field.

Download table as:  ASCIITypeset image

Next, we move to autocorrelation measurements. Figure 3 shows the clustering statistics of critical points in configuration space at z = 0 for three different rarity thresholds (${ \mathcal R }=5,10,20$), computed with the LS estimator via Equation (4). From left to right, minima, wall-type saddles, filament-type saddles, and maxima are displayed for the three cosmologies considered in our analysis. BAO peaks are clearly visible in all of the panels, with the vertical dashed gray lines marking their exact scale (expressed with r in h⁻¹ Mpc). Note that the positions of the BAO peaks are all the same, independently of critical point type, neutrino mass, and rarity—namely, 102.5 h⁻¹ Mpc, which precisely coincides with the BAO expected location in the reference cosmology. This remarkable aspect is a clear indication that critical points trace the BAO peak similarly to DM, halos, and galaxies and are faithful representations of their complementary structures.

Zoom In Zoom Out Reset image size

Several interesting features can be inferred from Figure 3. First, it is evident that a nonzero neutrino mass generally corresponds to higher BAO amplitudes, regardless of the specific critical point type, and those amplitudes increase as M_ν is augmented—namely, the larger the neutrino mass, the bigger the BAO amplitudes; this implies that it may be possible to infer neutrino mass signatures from such differential measurements, as we argue in this section. Moreover, the corresponding BAO amplitudes of minima and wall-type saddles are smaller than those of maxima and filament-type saddles. The similarity in shape between minima and saddles 1 (between maxima and saddles 2) can be explained by considering how critical points below (above) the threshold are selected—as detailed in Section 3.4. In fact, for minima and wall-type saddles the specific critical point abundance is determined from the lower limit of the density field, while the opposite is done for maxima and filament-type saddles. Generally, extrema (i.e., minima and maxima) show noisier and less smooth correlation function shapes when compared to saddle-point clustering: notice, in fact, that their corresponding error bars in Figure 3 are much larger. This is readily explained by the fact that the time evolution of extrema is more nonlinear than that of saddle points, as also reported by Cadiou et al. (2020) and Kraljic et al. (2022). In addition, more marked features in the BAO peaks and noisier correlation function shapes are seen with smaller rarity thresholds (see, e.g., the top panels of the figure, where ${ \mathcal R }=5$). This finding can be simply interpreted via linear bias, meaning that the more the tracer gets biased, the stronger the clustering (see, e.g., Kaiser 1984; Desjacques et al. 2018; Shim et al. 2021; Kraljic et al. 2022). Note that the effect of a rarity cut is similar to that of smoothing (see Appendix A, as well as the previous references): an increase in smoothing implies a decrease in the number of volume elements along with an increase in bias, and consequently the spatial correlation function gets noisier because of enlarged statistical uncertainties, and its features appear more enhanced.

We then focus on two key features that can be inferred from our two-point correlation measurements in configuration space (i.e., Figure 3): (1) the amplitudes of the various BAO peaks as a function of ${ \mathcal R }$, and (2) the spatial locations of their corresponding inflection points at large spatial separations. The first aspect is quantified via Figures 4 and 5, while the second one is characterized by Figures 6 and 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Specifically, Figure 4 shows the measured BAO amplitudes (r² ξ) at z = 0 along the x-axis, versus their corresponding spatial position marked by horizontal gray dashed lines, for the three cosmological scenarios and the three rarity thresholds examined. From top to bottom, the four panels refer to minima, saddles 1, saddles 2, and maxima, respectively. Note in particular that the differences in terms of BAO amplitudes between a massless neutrino framework and a cosmology with M_ν = 0.1 eV are rather small, especially for the case of minima—for which their values almost coincide when ${ \mathcal R }=10;$ for this reason, in the top panel the red and blue triangles are slightly displayed along the y-direction, for the sake of clarity.

Such differences are better quantified in Figure 5, in terms of percentages. That is, we now display relative ratios among BAO amplitudes in massive neutrino cosmologies confronted with their correspondent values in a massless neutrino scenario, as a function of ${ \mathcal R }$. Colors refer to the same cosmological models considered in Figure 4, but now the various symbols are also used to highlight those models, as reported in the figure. This visualization is quite useful, as it allows one to readily assess the effect of a nonzero neutrino mass on the BAO peaks inferred from the clustering of critical points in configuration space. In particular, as evident, departures from analogous measurements carried out in the baseline M_ν = 0.0 eV framework can reach up to ∼7% at z = 0 when M_ν = 0.1 eV and are much more pronounced for higher neutrino mass values.

From Figures 4 and 5 we can deduce a number of significant features. First, as in Shim et al. (2021), we find an amplification of the BAO peaks with rarity, namely, the amplitude of ξ is higher (i.e., stronger clustering) for lower values of ${ \mathcal R }$, regardless of the critical point type. And BAO features are further amplified by massive neutrinos, as quantified in the two plots. These effects are consistent with—and generalize—the findings of Peloso et al. (2015), who characterized the impact of neutrino masses on the shape and height of the BAO peak of the matter correlation function in real and redshift space, which contains relevant cosmological information; in the nonlinear regime the BAO peak increases with increasing M_ν , and up to 1.2% at z = 0 when M_ν = 0.30 eV, as reported by the same authors. Our approach based on critical points offers a more global multiscale perspective, since critical points carry remarkable topological properties and are good approximations for their corresponding LSS. Critical points are also less sensitive to systematic effects, and this is among the reasons why our autocorrelation BAO peak measurements show that departures from the corresponding massless neutrino scenario are more significant (up to 7%) even when M_ν = 0.1 eV—which is a value closer to the current stringent constraints on the summed neutrino mass reported in the literature. Furthermore, the spatial position of the BAO features in the various autocorrelations is robustly defined (i.e., 102.5 h⁻¹ Mpc) for all of the abundances considered, independently of critical point type and neutrino mass (hence an excellent standard ruler). In addition, generally saddle-point statistics are more advantageous to use for extracting cosmological information because their cosmic evolution is less nonlinear (Gay et al. 2012; Shim et al. 2021); for example, the two-point autocorrelations of walls provide precious information on the characteristic sizes of voids. Finally, we note that all of the autocorrelation functions go to zero at large scales through an inflection point, ¹⁴ an aspect that we address next.

To this end, Figure 6 shows another interesting and remarkable aspect of our analysis: not only the BAO amplitudes of all of the critical point autocorrelations are sensitive to massive neutrinos, but also their correspondent inflection scales—defined as the spatial positions at large r where the correlation functions ξ of a given critical point type computed at different rarities ${ \mathcal R }$ intersect (or are minimally distant), which also coincide with the spatial locations ${r}_{\inf }$ where $\xi ({r}_{\inf })\equiv 0$ (i.e., zero-crossings) within error bars—are altered by a nonzero neutrino mass. Specifically, in Figure 6 we zoom into the 110 h⁻¹ Mpc < r < 140 h⁻¹ Mpc interval and display the various ξ(r) measurements at z = 0 and ${ \mathcal R }=5,10,20$ for minima, wall-type saddles, filament-type saddles, and maxima—from left to right. All of the correlation function measurements are obtained with a "refined zoom-in" technique: in essence, for a given cosmology and within the previously specified spatial interval, we recompute the two-point autocorrelations shown in Figure 3 for the three chosen rarity levels using a finer bin size (1 h⁻¹ Mpc) and average the results over 100 independent QUIJOTE realizations. The inflection points are subsequently determined, and in the panels the vertical gray dashed lines indicate the spatial position of such points, while the shaded horizontal error bars highlight variations of ξ by ±0.1%. For minima, the determination of this scale appears more challenging, but it still falls within the autocorrelation error bars; this is likely due to the relatively small resolution of the simulations used in our study, which impacts underdense regions more severely than other cosmic web components, coupled with the fact that minima depart more significantly from linear theory.

The spatial positions of these inflection points at z = 0 are reported in Figure 7, split by corresponding type. The shaded horizontal error bars indicate the levels where ξ varies by ±0.1%; as evident, all of the inflection point spatial positions fall within this range.

Zoom In Zoom Out Reset image size

The notable results presented in Figures 6 and 7 deserve some further attention. In particular, the interesting fact that such inflection points, independently of rarity, intersect also at the zero-crossing scale where $\xi ({r}_{\inf })\equiv 0$ within error bars suggests that their positions can be well described by standard linear theory, and we will return to this theoretical aspect in a forthcoming publication. Furthermore, these inflection points are also quite sensitive to massive neutrinos, with their spatial position increasing with augmented neutrino mass—within the interval 120–135 h⁻¹ Mpc. As in Shim et al. (2021), it is then intriguing to attempt an analogy, in a multiscale perspective, with the LP of the correlation function (Anselmi et al. 2016), which is likewise subject to massive neutrino effects (Parimbelli et al. 2021). Figure 7 shows indeed a clear sensitivity to a nonzero neutrino mass: through these inflection points it appears feasible to simultaneously quantify the perhaps unique signatures of M_ν on the four key web constituents (i.e., halos, filaments, walls, and voids).

We then move to cross-correlations and carry out a similar analysis to that performed for the autocorrelation case at z = 0 in massless and massive neutrino cosmologies, for the same rarity thresholds previously considered. In general, cross-correlations are helpful in enhancing the signal-to-noise ratio (S/N), if used in combination with autocorrelations, and in mitigating the impact of systematics. In what follows, we compute all of the possible cross-combinations among extrema and saddle points; we show their clustering measurements in Figures 8 and 9.

Zoom In Zoom Out Reset image size

Specifically, Figure 8 displays cross-correlations among overdense and underdense critical points having opposite overdensity sign—i.e., ${ \mathcal P }{ \mathcal V },{ \mathcal P }{ \mathcal W },{ \mathcal V }{ \mathcal F },{ \mathcal W }{ \mathcal F }$, from left to right. The computations are carried out in configuration space at z = 0, as a function of rarity threshold and for the identical models considered in Figure 3, adopting the LS estimator (Equation (4)). Clearly, cross-correlation shapes differ from those of autocorrelations (i.e., Figure 3). In essence, they are mirrored with respect to the x-axis, namely, the clustering is negative (antibiased, or ξ(r) < 0 in the r-interval of interest) until reaching zero at larger spatial separations. Hence, at scales relevant for our analysis, these cross-correlations are always negative, implying that overdense and underdense critical points are anticorrelated. Therefore, BAO features are now "reversed" and appear as dips (rather than peaks) always at r = 102.5 h⁻¹ Mpc. Note that even in this case the spatial positions of the BAO dips are all identical, independently of critical point type pair, neutrino mass, and rarity. And, similarly to the autocorrelation measurements, BAO dips become more pronounced (i.e., showing a steeper negative amplitude) with increasing neutrino mass. As also pointed out by Shim et al. (2021) and Kraljic et al. (2022), anticlustering arises because these critical point pairs are oppositely biased tracers of the underlying DM density field.

Figure 9 shows instead cross-correlations among overdense or underdense critical point pairs (namely, ${ \mathcal P }{ \mathcal F }$ and ${ \mathcal V }{ \mathcal W }$)—characterized by an identical overdensity sign (i.e., similarly biased tracers). Since now critical point pairs are described by the same overdensity sign, the overall shapes of their two-point clustering are comparable to those of autocorrelations (see again Figure 3). Hence, there is no anticlustering at small separations, BAO peaks (local maxima) are detected at r = 102.5 h⁻¹ Mpc, and ξ(r) approaches zero at large spatial separations. Also in this case, BAO peaks are further enhanced by massive neutrinos. Moreover, we note that ${ \mathcal P }{ \mathcal F }$ cross-correlations show higher BAO amplitudes than ${ \mathcal V }{ \mathcal W }$ cross-correlations at a fixed rarity threshold.

Zoom In Zoom Out Reset image size

Figures 8 and 9 highlight another interesting aspect of our analysis, namely that the spatial positions of the BAO dips/peaks detected in cross-correlations are also robustly defined at the scale r = 102.5 h⁻¹ Mpc (as for autocorrelations), independent of critical point type pair, neutrino mass, and rarity. Moreover, cross-correlations of similarly biased tracers exhibit a behavior comparable to those of autocorrelations, while BAO features are instead "reversed" and manifest as dips for oppositely biased tracers. In both scenarios, a nonzero neutrino mass causes an overall enhancement of the amplitudes of the corresponding peaks/dips found in cross-correlations, much more pronounced for higher neutrino mass values. In addition, cross-correlations of minima with other critical point types show less noisy shapes (i.e., smaller error bars), if compared with the autocorrelations of minima alone—which, on the contrary, present the noisiest correlation shapes. In addition, ${ \mathcal W }{ \mathcal F }$ cross-correlations display smooth and less noisy clustering than extrema cross-correlations. Finally, cross-correlations involving maxima and saddle points are also less noisy than maxima autocorrelations alone, emphasizing once more the benefits of cross-correlations.

Next, as previously done for autocorrelation measurements, we focus on two key aspects that can be inferred from the two-point cross-correlation estimations in configuration space (Figures 8 and 9), namely, (1) the amplitudes of the BAO dips/peaks as a function of ${ \mathcal R }$, and (2) the spatial locations of their corresponding inflection points at large r-separations (as defined in the previous section). The first aspect is addressed in Figures 10 and 11, and the second one is quantified via Figures 12 and 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In detail, Figure 10 shows the amplitudes (r² ξ) of the BAO dips or peaks at z = 0 for the various cross-correlations considered along the x-axis, versus their corresponding spatial position marked by horizontal gray dashed lines, for the three cosmologies examined and ${ \mathcal R }=5,10,20$. A clear trend as a function of M_ν can be readily inferred: in essence, BAO dip/peak amplitudes are further enhanced (in absolute value terms) by the presence of massive neutrinos, and the enhancement is more significant with increased neutrino mass. Note also that r² ξ is negative in the first four panels since we are considering antibiased critical points, while r² ξ is positive for similarly biased tracers.

Figure 11 better quantifies such differences in terms of percentages, in analogy to Figure 5 for autocorrelations—namely, we display relative ratios among cross-correlation BAO amplitudes in massive neutrino cosmologies, normalized by the correspondent values in a massless neutrino scenario, as a function of ${ \mathcal R }$. As evident from the plot, departures from analogous cross-correlation measurements carried out in the baseline M_ν = 0.0 eV framework can reach even ∼9% at z = 0 when M_ν = 0.1 eV and are much more distinct for higher M_ν values. Interestingly, cross-correlations involving minima seem to be more effective in distinguishing small neutrino mass signatures (i.e., when M_ν = 0.1 eV) than autocorrelations of minima alone.

Finally, we examine the cross-correlation inflection points at large scales. In this respect, Figure 12 illustrates the remarkable fact that all of the inflection scales of the critical point pair cross-correlations are also sensitive to neutrino mass effects, independently of rarity (note the analogy with Figure 6). Specifically, we consider the 110 h⁻¹ Mpc < r < 140 h⁻¹ Mpc interval, recompute all of the two-point cross-correlations shown in Figures 8 and 9 with a refinement procedure, and average the results over 100 independent QUIJOTE realizations at z = 0. Inflection points are subsequently determined, as indicated in the panels by the vertical gray dashed lines, where from left to right we display ${ \mathcal P }{ \mathcal V },{ \mathcal P }{ \mathcal W },{ \mathcal V }{ \mathcal F },{ \mathcal W }{ \mathcal F },{ \mathcal P }{ \mathcal F },{ \mathcal V }{ \mathcal W }$—while the shaded horizontal error bars highlight the levels where ξ varies by ±0.1%.

The spatial positions of the inflection points determined from all of the possible cross-correlations at z = 0 are also reported in Figure 13, split by corresponding pair type. With the exception of the ${ \mathcal P }{ \mathcal V }$ case when M_ν = 0.1 eV (top panel), we detect a clear trend related to the inflection point positions in neutrino cosmologies—namely, the spatial locations of all of the inflection points are amplified by massive neutrinos, as a function of their mass. The apparent disagreement with this trend represented by the ${ \mathcal P }{ \mathcal V }$ cross-correlation measurement when M_ν = 0.1 eV (top panel) is likely related to the challenges in determining the inflection scale using minima (see for comparison the left panels of Figure 6), considering also their selection procedure (i.e., points below the rarity threshold; see the details in Section 3.4). Remarkably, even in cross-correlations we find an enhanced sensitivity of the inflection scales to massive neutrino effects.

Zoom In Zoom Out Reset image size

Finally, we briefly address redshift evolution effects, in relation to the clustering properties of critical points. To this end, Figures 14 and 15 provide some examples for autocorrelations: the position of the BAO peaks (r = 102.5 h⁻¹ Mpc, independent of redshift) is marked by vertical dashed gray lines.

Zoom In Zoom Out Reset image size

Specifically, Figure 14 displays the clustering statistics of critical points in configuration space at z = 1 (top panels), z = 2 (middle panels), and z = 3 (bottom panels), for a rarity choice corresponding to ${ \mathcal R }=20$ and varying cosmologies. We adopt the same stylistic conventions as in Figure 3. In the various panels, we maintain the same x-y scale for a better visual characterization of the redshift impact in the clustering properties. At any given redshift and for all of the critical point types, a nonzero neutrino mass causes an enhancement of the BAO peak amplitudes: the enhancement is more significant for higher neutrino mass values (see, e.g., the M_ν = 0.4 eV case).

To appreciate possible redshift evolution effects more clearly, in Figure 15 we now show autocorrelation functions for the same rarity threshold previously examined and within a fixed cosmological model—while we vary the redshift in each panel by considering z = 1 (purple), z = 2 (dark green), and z = 3 (light green), respectively, as specified by different colors in the figure inset. Similarly as in Figure 14, from left to right we indicate minima, wall-type saddles, filament-type saddles, and maxima.

Zoom In Zoom Out Reset image size

Around the scale of BAO peaks, all of the autocorrelations show small evolution with redshift. And, as expected, z-evolution (although weak) is more relevant for the autocorrelations involving the most nonlinear critical points (peaks, voids), while those involving less nonlinear critical points (filaments, walls) appear mostly insensitive to redshift evolution effects. In particular, the evolution seems more prominent for minima. Moreover, we find a similar trend in cross-correlation measurements.

Therefore, as also noted in Shim et al. (2021), the overall insensitivity to redshift, especially for saddle statistics, makes the clustering of critical points a topologically robust alternative to more standard clustering methods.