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

Critical points represent a subset of special points tracing cosmological structures, carrying remarkable topological properties. They thus offer a richer high-level description of the multiscale cosmic web, being more robust to systematic effects. For the first time, we characterize here their clustering statistics in massive neutrino cosmologies, including cross-correlations, and quantify their simultaneous imprints on the corresponding web constituents - i.e., halos, filaments, walls, and voids - for a series of rarity levels. Our first analysis is centered on a density-threshold-based approach in configuration space. In particular, we show that the presence of massive neutrinos does affect the baryon acoustic oscillation peak amplitudes of all of the critical point correlation functions above/below the rarity threshold, as well as the positions of their correspondent inflection points at large scales: departures from analogous measurements carried out in the baseline massless neutrino scenario can reach up to ~7% in autocorrelations and ~9% in cross-correlations at z=0 when M_nu=0.1 eV, and are more pronounced for higher neutrino mass values. In turn, these combined multiscale effects can be used as a novel technique to set upper limits on the summed neutrino mass and infer the type of hierarchy. Our study is particularly relevant for ongoing and future large-volume redshift surveys such as the Dark Energy Spectroscopic Instrument and the Rubin Observatory Legacy Survey of Space and Time, which will provide unique datasets suitable for establishing competitive neutrino mass constraints.


INTRODUCTION
Neutrinos represent a clear indication of physics beyond the standard model, since the confirmation via oscillation experiments that they are massive particles -see, e.g., Lesgourgues & Pastor (2006) and Gonzalez-Garcia & Maltoni (2008) for seminal reviews on the phenomenology of massive neutrinos. In this respect, the baseline ΛCDM cosmological framework characterized by massless neutrinos (or at best by a minimal non-zero neutrino mass), a spatially flat cosmology dominated by collisionless cold dark matter (CDM), and a dark energy (DE) component in the form of a cosmological constant (Λ), should be extended accordingly.
Therefore, it comes with no surprise that determining the neutrino mass scale and type of hierarchy are among the major enterprises of all of the ongoing and upcoming large-volume astronomical surveys, and one of the primary targets of future space missions. This is for example the case of the Dark Energy Survey ( Takada et al. 2014), the Nancy Grace Roman Space Telescope (Roman; Spergel et al. 2015), and the Euclid Consortium (Laureĳs et al. 2011), as well as of next-generation cosmic microwave background (CMB) experiments such as CMB Stage 4 (CMB-S4; Abazajian et al. 2016Abazajian et al. , 2019Abitbol et al. 2017)  To this end, cosmology is mainly sensitive to the summed neutrino mass -denoted as ≡ i i throughout this paper, with i ( = 1, 2, 3) representing the three individual masses of active neutrinos -and provides competitive upper bounds on . On the other hand, particle physics experiments set stringent lower bounds (i.e., ≥ 0.05 eV). Moreover, flavor oscillations bring no information on the absolute neutrino mass scale and are insensitive to individual neutrino masses, as they allow only for a measurement of squared differences. Namely, Δ 2 21 = 2 2 − 2 1 +7.50 × 10 −5 eV 2 can be obtained from solar neutrino oscillations, and Δ 2 31 = 2 3 − 2 1 Δ 2 32 = 2 3 − 2 2 ±2.45 × 10 −3 eV 2 can be inferred from atmospheric neutrino oscillations (Gonzalez-Garcia & Maltoni 2008). Since we only know experimentally that Δ 21 > 0, as well as the value of |Δ 31 |, this leads to the possibility of either normal hierarchy (NH; 1 < 2 3 ) or inverted hierarchy (IH; 3 1 2 ) for the mass ordering. In the NH configuration, the minimal sum of neutrino masses is = 0.057 eV, while in the IH configuration the minimal summed mass is = 0.097 eV. The latest neutrino mass upper bounds reported by the Planck collaboration are < 0.12 eV (Planck Collaboration et al. 2020a), while the most constraining combination of data including Data Release 16 (DR16) from the extended Baryon Oscillation Spectroscopic Survey (eBOSS; Dawson et al. 2016), part of the fourth generation of the Sloan Digital Sky Survey (SDSS-IV; York et al. 2000;Blanton et al. 2017), gives the upper limit on the sum of neutrino masses at < 0.115 eV (eBOSS Collaboration et al. 2021). Clearly, such stringent bounds put considerable pressure on the IH scenario as a plausible possibility for the neutrino mass ordering.
A direct neutrino mass detection, or at least more compet-itive upper limits on are expected in the next few years from Stage-IV cosmological experiments. For example, DESI will be able to measure with an uncertainty of 0.020 eV for max < 0.2ℎMpc −1 , sufficient to make the first direct detection of at more than 3 significance and rule out the IH at the 99% confidence level (CL) if the hierarchy is normal and the masses are minimal (DESI Collaboration et al. 2016). Moreover, future CMB-S4 measurements combined with late-time measurements of galaxy clustering and cosmic shear from the Rubin Observatory LSST would allow one to achieve a 3 detection of the minimal mass sum min = 0.06 eV (Mishra-Sharma et al. 2018). In synergy with neutrino decay experiments sensitive to the electron neutrino mass (i.e., KATRIN;KATRIN Collaboration 2001), and neutrinoless double decay probes sensitive to the effective Majorana mass (i.e., Kam-LAND, GERDA, Cuore; Gando et al. 2013;Ackermann et al. 2013;CUORE Collaboration et al. 2015), it should be possible to finally close on the neutrino mass scale and type of hierarchy within this decade. For extensive details to this end, see, i.e., Gariazzo et al. (2018).
Constraints on massive neutrinos from cosmology can be obtained exploiting a variety of tracers and probes, spanning different scales. The traditional and perhaps most straightforward route is via the CMB, typically involving CMB gravitational lensing or using the early integrated Sachs-Wolfe (ISW) effect in polarization maps -see, e.g., Battye & Moss (2014), Brinckmann et al. (2019), Planck Collaboration et al. (2020b, and references therein. Regarding baryonic tracers of the large-scale structure (LSS), a plethora of techniques have been proposed in the literature to date. Among them, a (largely incomplete) recent list includes three-dimensional (3D) galaxy clustering via the matter power spectrum (Angulo et al. 2021;Bose et al. 2021), bispectrum Hahn & Villaescusa-Navarro 2021), and marked power spectrum , cosmic shear, peaks, and bispectrum through weak lensing (Coulton et al. 2019;Ajani et al. 2020), galaxy clusters via the Sunyaev-Zel'dovich (SZ) effect (Roncarelli et al. 2017;Bolliet et al. 2020), peculiar velocities (Whitford et al. 2022), voids statistics (Kreisch et al. 2019;Zhang et al. 2020), the Lyman-(Ly ) forest flux power spectrum (Seljak et al. 2005(Seljak et al. , 2006Viel et al. 2010;Rossi 2017Rossi , 2020, and much more. The common denominator among all these different methodologies is to rely on a single observable, and/or on a specific scale (or a limited scale-range). Only recently, there have been some interesting attempts to use instead combinations of probes -see for example the work of Bayer et al. (2021).
Furthermore, gaining a deeper theoretical understanding of the impact of massive neutrinos on structure formation (particularly at small scales) and on the primary observables typically used to characterize neutrino mass effects are necessary, to obtain robust constraints free from systematic biases. Consistent progress has been made in this direction over the last few years. Relevant examples include the detailed assessment of neutrino effects on Ly forest observables (Rossi 2017), the quantification of the consequences of non-zero neutrino masses and asymmetries on dark matter halo assembly bias (Lazeyras et al. 2021;Wong & Chu 2022), the implications of massive neutrinos in the realm of the Hubble tension (Di Valentino & Melchiorri 2022), the challenges related to parameter degeneracies when combining CMB and LSS probes in massive neutrino cosmologies (Archidiacono et al. 2017), and the estimation of the impact of massive neutrinos on the baryon acoustic oscillation (BAO) peak (Peloso et al. 2015) and on the linear point (LP) of the spatial correlation function (Parimbelli et al. 2021).
On an apparently unrelated subject, recently an interesting analysis of the cosmic web in terms of critical points (i.e., a subset of special points in position-smoothing space, tracing cosmological structures) has been proposed by Cadiou et al. (2020) , and the corresponding clustering properties of those points have been characterized in Shim et al. (2021) and in Kraljic et al. (2022) as a function of rarity threshold in the baseline ΛCDM model. Critical points (extrema and saddles, where the spatial gradient of the density field vanishes) carry remarkable topological properties, and provide a more fundamental view of the cosmic web as a whole -in a multiscale perspective. This is primarily because the topology of the initial density field, at a fixed smoothing scale, is encoded in the positions and heights of such points: hence, in principle, it it possible to predict the evolution of cosmic web structures at later times from their clustering properties jointly with the power spectrum of the underlying initial Gaussian field. Moreover, the special scales at which two critical points coalesce produce merging effects corresponding to the mergers of halos, filaments, walls, and voids (Cadiou et al. 2020). In addition, the drift of critical points with smoothing defines the skeleton tree (Sousbie et al. 2008;Pogosyan et al. 2009;Gay et al. 2012), capturing topological variations with increasing smoothing scale -that acts effectively as a time variable. Besides being more robust to systematic effects, critical points are thus useful because they represent a meaningful and efficient compression of information of the 3D density field, capturing its most salient features -as they carry significance in terms of cosmology and/or galaxy formation. In fact, with some care (see our discussion later on in Section 4), critical points can be associated to corresponding physical structures of the same type, although their typical size will depend on the smoothing scale of the density field. In particular, and of direct relevance to this work, Shim et al. (2021) found an amplification of the BAO features with increasing density threshold in the autocorrelations of critical points (reversed for cross-correlations) within the ΛCDM, and reported of an 'inflection scale' around ∼ 133ℎ −1 Mpc which appears to be common to all of the autoand cross-correlation functions.
Here, we combine all of the previous (seemingly unrelated) topics in a coherent framework, and characterize, for the first time, the clustering statistics of critical points in massive neutrino cosmologies -addressing their sensitivity to small neutrino masses, with the goal of identifying multiscale signatures. Specifically, we compute the auto-and cross-correlation functions of critical points in configuration space, for a series of rarity thresholds, and also quantify redshift evolution effects. Our main analysis is centered on BAO scales, and it is focused on two key aspects: (1) the multiscale effects of massive neutrinos on the BAO peak amplitudes of all of the critical point correlation functions above/below rarity threshold; and, (2) the multiscale effects of massive neutrinos on the spatial positions of their correspondent correlation function inflection points at large scales. The first aspect is inferred from the fact that massive neutrinos do impact the BAO peak (Peloso et al. 2015), while the second one is inspired by a sort Note though that their primary focus is on critical events, not to be confused with critical points. See Section 4.2 for the definition of such inflection points at large spatial separations. C P M N 3 of resemblance (although in a multiscale perspective) with the LP of the spatial correlation function, which is also affected by a non-zero neutrino mass (Parimbelli et al. 2021). We carry out our measurements using a subset of realizations from the QUIJOTE suite , as explained in Section 2. In particular, we utilize full snapshots at redshifts = 0, 1, 2, 3 for a choice of representative neutrino masses, in order to assess redshift evolution effects. Our first analysis is centered on a 'density thresholdbased' approach: the methodology to construct density fields from the output of -body simulations, extract and classify critical points, and perform clustering measurements in configuration space is thoroughly explained in Section 3. Noticeably, in this study we consider three cuts in rarity R, namely R ∈ {5, 10, 20} -where the abundances are expressed in percentages.
Our main results are detailed in Section 4. Here, we show that the presence of massive neutrinos do affect the BAO peak amplitudes of all of the critical point correlation functions above/below rarity threshold, as well as the positions of their correspondent inflection points at large scales. Departures from analogous measurements carried out in the baseline massless neutrino scenario can reach up to ∼ 7% in autocorrelations and ∼ 9% in cross-correlations at = 0 when = 0.1 eV, and are more pronounced for higher neutrino mass values. In a companion publication, we show how these combined multiscale effects derived from the clustering of critical points can be used as a novel powerful technique (complementary to more traditional methods) to set upper limits on and infer the type of hierarchy. This work is the first of a series of investigations that aim at exploring the sensitivity of critical points and critical events to massive neutrinos, and more generally in relation to the dark sector. The layout of the paper is organized as follows. Section 2 briefly describes the simulations considered in our study. Section 3 presents the methodology at the core of the 'density threshold-based' approach. The main results of our analysis are provided in Section 4, including cross-correlation measurements as well as the characterization of redshift evolution effects. We conclude in Section 5, where we summarize the various results and highlight novel aspects -along with ongoing applications and future avenues. We leave in Appendix A a number of technicalities on the analysis performed, related in particular to our choices of bin size and smoothing scale of the density field, and provide some complementary tables in Appendix B.

SIMULATIONS
In this work, we use a subset of realizations from the QUI-JOTE suite. The QUIJOTE simulations (Villaescusa-Navarro et al. 2020) are a set of 44,100 publicly available -body runs spanning over 7000 cosmological models, following the gravitational evolution of CDM CDM plus neutrino particles (if present) with different resolutions (i.e., 256 3 , 512 3 , 1024 3 particles per species). Initial conditions (ICs) are set at redshift in = 127 either adopting the Zel'dovich approximation (Zel'dovich 1970) or via second order Lagrangian Perturbation Theory (2LPT). Moreover, the QUIJOTE suite contains standard, fixed, and paired fixed simulations, the difference being in the way ICs are generated. The input matter power spectra and transfer functions are obtained via the Boltzmann code CAMB (Lewis et al. 2000). All of the realizations are produced https://quijote-simulations.readthedocs.io/en/latest/ with the Three Particle Mesh (TreePM) code GADGET-III (Springel 2005), over a periodic box size of 1ℎ −1 Gpc, for a total volume of 44, 100 (ℎ −1 Gpc) 3 . Snapshots are saved at = 0, 0.5, 1, 2, 3, respectively. The gravitational softening is set to 1/40 of the mean inter-particle distance for all of the particle species, and the ICs random seeds are the same for an identical realization in different models, but vary from realization to realization within the same model. The baseline cosmological parameters are closer to those reported by the Planck Collaboration in 2018 (Planck Collaboration et al. 2020a), representing a flat (i.e., Ω Λ + Ω m = 1) massless neutrino ΛCDM model with the total matter density Ω m = 0.3175, the DE density parameter Ω Λ = 0.6825, the baryon density Ω b = 0.0490, the primordial scalar spectrum power-law index s = 0.9624, the Hubble parameter ℎ = 0.6711, the DE equation of state (EoS) parameter = −1, and the primordial power spectrum amplitude s = 2.13 × 10 −9 at the pivot scale pivot ≡ 0 = 0.05 Mpc −1 . For convenience, these reference parameters are reported in the upper part of Table 1. Note that the normalization choice for s in the baseline massless neutrino cosmology implies 8 = 0.834 at = 0, with 8 the linear theory root mean square (rms) matter fluctuation in 8ℎ −1 Mpc spheres; by construction, all of the other realizations are tuned to match an identical 8 value at the present time.
In massive neutrino cosmologies, the rescaling method developed by Zennaro et al. (2017) is used to determine initial conditions. In addition, all of the massive neutrino realizations assume degenerate masses: the summed neutrino mass is indicated with in Table 1. Regarding neutrino implementations within the -body framework, all of the the neutrino runs are performed employing a popular particle-based approach, treated as a collisionless and pressureless fluid (in the same fashion as CDM). Moreover, the forces at small scales for the neutrino species are also properly computed. The particlebased approach, despite shot noise challenges, is automatically able to capture the full nonlinear neutrino clustering and it accurately reproduces the nonlinear evolution at small scales -a crucial aspect for obtaining reliable and robust cosmological constraints. See Rossi (2017Rossi ( , 2020 for additional details on this implementation. In our study, we use 100 independent realizations (i.e., having different initial random seeds) from the QUIJOTE suite for any of the following cosmologies: fiducial framework (massless neutrino scenario or 'best guess'), massive neutrino model with = 0.1 eV, and massive neutrino model with = 0.4 eV, respectively. Specifically, we exploit the runs denoted as standard, characterized by Zel'dovich ICs at = 127. The resolutions considered are relatively low, 512 3 particles per species, on a periodic box size of 1ℎ −1 Gpc. We utilize full snapshots at = 0, 1, 2, 3, in order to also assess redshift evolution effects. As previously noted, 8 ( = 0) = 0.834 in all of the runs; hence, s varies with different neutrino masses. This normalization convention is what has been termed 'NORM' in Rossi (2020), and is well-motivated observationally, since 8 at the present epoch is effectively dictated by observational constraints. The main characteristics of the simulations adopted here are reported in the bottom part of Table 1, and more technical details can also be found in Villaescusa-Navarro et al. (2020).
As a visual example, Figure 1 shows a spatial projection of the dark matter (DM) overdensity field along the − plane with a depth of 50ℎ −1 Mpc across the -axis, obtained from three selected QUIJOTE simulations at = 0. From left to right, the plots refer to the fiducial massless neutrino cosmology run ( = 0.0 eV), the massive neutrino model with = 0.1 eV, and the massive neutrino sce-nario with = 0.4 eV, respectively. Top panels display 100 × 100 [ℎ −1 Mpc] 2 slices of the simulation cubes, while bottom panels are enlargements of a small 10 × 10 [ℎ −1 Mpc] 2 inset with identical spatial location and the same thickness as in the top figures. The normalized density fields are rendered using a rainbow palette, with the corresponding values indicated by the right side color bar. While rather small, morphological and topological differences induced by a non-zero neutrino mass are visually perceptible: they are quantifiable via the clustering statistics of critical points in the complex structure of the evolving cosmic web network, as we show in Section 4.

METHODOLOGY
In this section, we briefly describe our technique to construct density fields from the output of -body simulations, extract and classify critical points, and subsequently perform clustering measurements in configuration space. We generically refer to this methodology as the 'density threshold-based' approach.

Density Threshold-Based Approach: Overview
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 -body simulations. In essence, in the first step DM particles are turned into cells of given density, that are properly smoothed up to a desired scale. Critical points are C P M N 5 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 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.

Construction of Density Fields
We begin with the production of density fields, which are constructed from the DM particle output of selected QUIJOTE -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 -body particles using a triangular-shape cloud (TSC) massassignment scheme to deposit particle masses into the grid; in this process, no weights are associated to each particle. Density contrast fields, defined as (x) = (x)/¯− 1 with the corresponding average background density, are readily inferred from the global knowledge of (x). We adopt grid = 1024 voxels per dimension on a regular grid spacing -corresponding to 2 p , with 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 grid , and found negligible dependence on our results.

Extraction and Classification of Critical Points
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 G over G = 3ℎ −1 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 on 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 cut-off, erasing small-scale structures: at a given -interval, structures of different sizes and masses are present, and therefore selecting an explicit G would correspond to averaging over specific density regions. A suitable smoothing choice is then 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: See https://github.com/franciscovillaescusa/Pylians See https://github.com/cphyc/py_extrema 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 (|x i − x c |) < 1 pixel -with = 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 (V or minima); 1 for walls (W or 'saddle 1' or walltype saddles); 2 for filaments (F or 'saddle 2' or filament-type saddles); 3 for peaks (P or maxima).

Clustering Measurements in Configuration Space
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 carries 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's first introduce the notation ∈ {0, 1, 2, 3} ≡ {V, W, F , P} to specify the critical point type. We will also use the indices ℓ and in place of , whenever necessary (i.e., ≡ ℓ, ). In essence, for peak-and filament-type (respectively, void-and wall-type) critical points, we select ensambles with higher (respectively, lower) than a given threshold , R , where 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: with ∈ {P, F } ≡ {3, 2} for peaks and filaments, respectively, and: with ∈ {W, V} ≡ {1, 0} for walls and voids, respectively. In the previous expressions, C, indicates the entire ensemble of critical points of type , while C, ( ≥ , R ) or C, ( ≤ , R ) represent the subset containing critical points of type above or below the selected threshold , R . For ease of notation, in what follows we define ≡ , R , and also indicate C, ( ≥ ) or C, ( ≤ ) generically with C, . In our Throughout the paper, we will refer to a specific critical point via its type, although the reader should keep in mind the important remark detailed in Section 4.1: the correspondence between density-based critical points and cosmological structures depends on the smoothing scale choice.
Note that throughout the paper we will refer indistinctly to either ,R or R to indicate the rarity level/threshold, with the latter quantity defined as in Equations 2 and 3 and expressed in percentage terms. J. M , G. R , H. Y (2022) study we consider three cuts in rarity, namely R ∈ {5, 10, 20} -where the abundances are expressed in percentages. We then compute the spatial clustering statistics of critical points above (below) threshold -i.e., real space two-point auto-and cross-correlations ( ) as a function of the separation -using the H 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 of critical points above (below) a selected density threshold having total size C, , and assuming a corresponding random catalog of total size R, characterized by a random uniform probability distribution of points within the same volume, the LS estimator for the two-point correlation function ,LS kℓ above (below) reads : where the normalization factors are specified by Note that all of the correlation measurements obtained via Equation (4) are reported as a function of the spatial separation , conventionally expressed in units of ℎ −1 Mpc. Unless specified otherwise, we always use a bin size of 5ℎ −1 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ℎ −1 Mpc in selected regions near such points.

DENSITY THRESHOLD-BASED APPROACH: RESULTS
This section contains the main results of our first analysis. After some basic considerations related to the 'density threshold-based' approach, we present here auto-and crosscorrelation measurements of the clustering of critical points in massive neutrino cosmologies for a series of rarity levels, confronted with analogous computations in the baseline massless neutrino model. In particular, we show how massive neutrinos affect the BAO peak amplitudes of all of the critical point correlation functions above/below rarity threshold, as well as See https://github.com/astropy/halotools We have also extensively tested the impact of increasing the size of random catalogs in this process, and found negligible effects on our results.
Throughout the rest of the paper, for ease of simplicity, we omit to indicate the threshold levels ≡ ,R in the notation of the two-point correlation functions. The adopted rarity choices will be readily distinguishable in the various plots presented, via the usage of different line styles, point shapes, or contrasting colors. the positions of their corresponding inflection points at large scales. Finally, we also address redshift evolution effects.
In what follows, unless specified otherwise, we always express rarity thresholds R in percentage and consider two massive neutrino models having = 0.1 eV and = 0.4 eV, respectively, besides the baseline massless neutrino scenario. Moreover, all of the measurements represent averages over 100 independent QUIJOTE realizations at fixed cosmology, and the associated errorbars are the corresponding 1 variations.

Abundance and Visualization of Critical Points in
Massive Neutrino Cosmologies 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 = 0, color-coded by type, for three different rarity thresholds. The squared patches in the figure are characterized by a side of 100ℎ −1 Mpc in the − plane, and a 50ℎ −1 Mpc depth along the -coordinate. The underlying normalized density fields are the results of sampling the density fields with Pylians assuming grid = 1024. Top panels refer to the baseline cosmology (massless neutrino scenario), middle panels are for a massive neutrino cosmology with = 0.1 eV, while bottom panels display results for = 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.
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 cut-off. To this end, see also the discussion in Shim et al. (2021), who adopted a smoothing scale of G = 6ℎ −1 Mpc corresponding to ∼ 10 15 M average density regions. As noted by the same authors, with such a choice of G only large voids can be resolved, and only 5% of all density peaks represent virialized galaxy clusters at = 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 R manifest into smaller dispersions/scatter in ( ), 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 ≠ 0.0 eV translate into morphological and topological differences that carry rich cosmological information.
In this regards, Table 2 reports the overall abundance of the entire set of critical points at = 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 number of peaks (P) and voids (V) is similar, as well as the number C P M N 7 Figure 2. Visualization of the spatial distribution of critical points (color-coded by their corresponding type) in massless and massive neutrino cosmologies at = 0, from small QUIJOTE-simulated cubic density patches of 100ℎ −1 Mpc side and 50ℎ −1 Mpc depth, as a function of rarity -with R increasing from left to right, as indicated in the various panels. Red, green, cyan, and blue are used respectively to display maxima, filament-saddles ('saddle 2'), wall-saddles ('saddle 1'), and minima. A 'density threshold-based' approach is adopted for the identification and classification of critical points, as reported in Section 3. See the main text for more details. of filaments (F ) versus walls (W). Moreover, in Gaussian fields, the ratio between filaments and peaks (F /P) or walls and voids (W/V) is estimated to be ∼ 3.05, and in addition the total number of extrema and saddle points is preserved at the first non-Gaussian perturbative order. Also, for sufficiently large volumes, the ratio between the number of peaks and walls over voids and filaments ([P/W]/[V/F ]) in a Gaussian field is very close to unity throughout the entire -evolution. This is because the genus topology, equal to the alternate sum of critical points, should be preserved, namely: with (i) cp the mean number density of peaks, filaments, walls, and voids -respectively.
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 non-zero 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 -   (2017), their Section 4.3.

Auto-Correlations of Critical Points in Massive Neutrino
Cosmologies Next, we move to auto-correlation measurements. Figure  3 shows the clustering statistics of critical points in configuration space at = 0 for three different rarity thresholds (R = 5, 10, 20), computed with the LS estimator via Equation (4). From left to right, minima, wall-type saddles, filamenttype saddles, and maxima are displayed, respectively, for the three cosmologies considered in our analysis. BAO peaks are clearly visible in all of the panels, with the vertical dashed grey lines marking their exact scale (expressed with in ℎ −1 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ℎ −1 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.
Several interesting features can be inferred from Figure 3. First, it is evident that a non-zero neutrino mass generally corresponds to higher BAO amplitudes, regardless of the specific critical point type, and those amplitudes increase as 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 (respectively, maxima and saddles-2) can be explained by considering how critical points below (above) 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 smoother correlation function shapes when compared to saddle-point clustering: notice in fact that their corresponding errorbars 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 for example the top panels of the figure, where R = 5). This finding can be simply interpreted via linear bias, meaning that the more the tracer gets biased, the stronger the clustering -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, and also 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): namely, (1) the amplitudes of the various BAO peaks as a function of 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.
Specifically, Figure 4 shows the measured BAO amplitudes ( 2 ) at = 0 along the -axis, versus their corresponding spatial position marked by horizontal grey 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 = 0.1 eV are rather small, especially for the case of minima -for which their values almost coincide when R = 10; for this reason, in the top panel the red and blue triangles are slightly displayed along the -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 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 non-zero neutrino mass on the BAO peaks inferred from the clustering of crit-C P M N 9    Figure 3, indicating a clear sensitivity to massive neutrinos. In the various panels, measurements of ( ) at = 0 are performed in the 110ℎ −1 Mpc < < 140ℎ −1 Mpc interval, using a refined binsize of 1ℎ −1 Mpc. Top panels are for the baseline massless neutrino model, middle panels refer to the scenario with = 0.1 eV, and bottom panels show the case of = 0.4 eV. Results are averaged over 100 independent QUIJOTE realizations per given cosmology. The three rarity cuts are represented by different colors as indicated in the panels, and minima, wall-type saddles, filament-type saddles, and maxima are reported from left to right. Shaded horizontal errorbars indicate ±0.1% variations in . Remarkably, the inflection scales (highlighted by the vertical grey dashed lines), independently of rarity, are altered by a non-zero neutrino mass.
cation of the BAO peaks with rarity: namely, the amplitude of is higher (i.e., stronger clustering) for lower values of 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 , and up to 1.2% at = 0 when = 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 auto-correlation BAO peak measurements show that departures from the corresponding massless neutrino scenario are more significant (up to 7%) even when = 0.1 eVwhich 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 auto-correlations is robustly defined (i.e., 102.5ℎ −1 Mpc) for all of the abundances considered, independently of critical point type and neutrino mass (hence an excellent standard ruler). Also, generally saddle-point statistics are more advan-tageous 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 auto-correlations of walls provide precious information on the characteristic sizes of voids. Finally, we note that all of the auto-correlation 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 auto-correlations are sensitive to massive neutrinos, but also their correspondent inflection scales -defined as the spatial positions at large where the correlation functions 's of a given critical point type computed at different rarities R intersect (or are minimally distant), which also coincide with the spatial locations inf where ( inf ) ≡ 0 (i.e., zero-crossings) within errorbars -are altered by a nonzero neutrino mass. Specifically, in Figure 6 we zoom into the 110ℎ −1 Mpc < < 140ℎ −1 Mpc interval and display the various ( )'s measurements at = 0 and R = 5, 10, 20 for minima, wall-type saddles, filament-type saddles, and maximafrom left to right. All of the correlation function measurements are obtained with a 'refined zoom-in' technique: in essence, We use here the same terminology of Shim et al. (2021) to denote such point, namely the spatial location inf where ( inf ) ≡ 0 at larger (i.e., zerocrossing of ), although our definition of inflection scale is slightly different, as reported in the main text. for a given cosmology and within the previously specified spatial interval, we recompute the two-point auto-correlations shown in Figure 3 for the three chosen rarity levels using a finer bin size (1ℎ −1 Mpc), and average the results over 100 independent QUIJOTE realizations. The inflection points are subsequently determined, and in the panels the vertical grey dashed lines indicate the spatial position of such points, while the shaded horizontal errorbars highlight variations of by ±0.1%. For minima, the determination of this scale appears more challenging, but it still falls within the auto-correlation errorbars: this is likely due to the relatively small resolution of the simulations used in our study, that 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 = 0 are reported in Figure 7, split by corresponding type. The shaded horizontal errorbars indicate the levels where varies by ±0.1%; as evident, all of the inflection point spatial positions fall within this range.
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 ( inf ) ≡ 0 within errorbars suggests that their positions can be well-described by standard linear theory, and we will return on 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]ℎ −1 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 non-zero neutrino mass: through these inflection points it appears feasible to simultaneously quantify the perhaps unique signatures of on the four key web constituents (i.e., halos, filaments, walls, and voids).

Cross-Correlations of Critical Points in Massive
Neutrino Cosmologies We then move to cross-correlations, and carry out a similar analysis as performed for the auto-correlation case at = 0 in massless and massive neutrino cosmologies, for the same rarity thresholds previously considered. In general, crosscorrelations are helpful in enhancing the signal-to-noise ratio (SNR) if used in combination with auto-correlations, 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.
Specifically, Figure 8 displays cross-correlations among overdense and underdense critical points having opposite overdensity sign -i.e., PV, PW, VF , WF , from left to right, respectively. The computations are carried out in configuration space at = 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 auto-correlations (i.e., Figure 3). In essence, they are mirrored with respect to the -axis: namely, the clustering is negative (anti-biased, or ( ) < 0 in the -interval of interest) until reaching zero at larger spatial separations. Hence, at scales relevant for our analysis, these crosscorrelations are always negative, implying that overdense and underdense critical points are anti-correlated. Therefore, BAO features are now 'reversed' and appear as dips (rather than peaks) always at = 102.5ℎ −1 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 auto-correlation 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), anti-clustering 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, PF and VW)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 twopoint clustering are comparable to those of auto-correlations (see again Figure 3). Hence, there is no anti-clustering at small separations, BAO peaks (local maxima) are detected at = 102.5ℎ −1 Mpc, and ( ) approaches zero at large spatial separations. Also in this case, BAO peaks are further enhanced by massive neutrinos. Moreover, we note that PF crosscorrelations show higher BAO amplitudes than VW crosscorrelations at a fixed rarity threshold. 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 = 102.5ℎ −1 Mpc (as for auto-correlations), independent of critical point type pair, neutrino mass, and rarity. Moreover, cross-correlations of similarly biased tracers exhibit a behavior comparable to those of auto-correlations, while BAO features are instead 'reversed' and manifest as dips for oppositely biased tracers. In both scenarios, a non-zero 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 errorbars), if compared with the auto-correlations of minima alone -which, on the contrary, present the noisiest correlation shapes. Also, WF cross-correlations display smooth and less noisy clustering than extrema cross-correlations. Finally, cross-correlations involving maxima and saddle points are also less noisier than maxima auto-correlations alone, emphasizing once more the benefits of cross-correlations. Next, as previously done for auto-correlation 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 R, and (2) the spatial locations of their corresponding inflection points at large -separations (as defined in the previous section). The first aspect is addressed in Figures 10 and 11, the second one is quantified via Figures  12 and 13.
In detail, Figure 10 shows the amplitudes ( 2 ) of the BAO dips or peaks at = 0 for the various cross-correlations considered along the -axis, versus their corresponding spatial position marked by horizontal grey dashed lines, for the three cosmologies examined and R = 5, 10, 20. A clear trend as a function of 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 2 is negative in the first four panels since we are considering anti-biased critical points, while 2 is positive for similarly biased tracers. Figure 11 better quantifies such differences in terms of percentages, in analogy to Figure 5 for auto-correlations: 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 R. As evident from the plot, departures from analogous cross-correlation measurements carried out in the baseline = 0.0 eV framework can reach even ∼ 9% at = 0 when = 0.1 eV, and are much more distinct for higher values. Interestingly, cross-correlations involving minima seem to be more effective in distinguishing small neutrino mass signatures (i.e., when = 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ℎ −1 Mpc < < 140ℎ −1 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 = 0. Inflection points are C P M N 13 Figure 9. Configuration-space cross-correlations of critical points characterized by an identical overdensity sign, as a function of rarity threshold, at = 0. Colors, line styles, rarity thresholds, and cosmologies same as those adopted in Figure 8. Left panels display cross-correlations of overdense critical point pairs (P F), while right panels are for underdense critical point cross-correlations (V W). Note that the shapes of the two-point crosscorrelation clustering are akin to those of auto-correlation measurements (i.e., Figure 3), since here we consider cross-correlations of similarly biased tracers. BAO features appear as peaks at = 102.5ℎ −1 Mpc, with the vertical dashed grey lines marking their exact scale.
subsequently determined, as indicated in the panels by the vertical grey dashed lines, where from left to right we display PV, PW, VF , WF , PF , VW, respectively -while the shaded horizontal errorbars highlight the levels where varies by ±0.1%. The spatial positions of the inflection points determined from all of the possible cross-correlations at = 0 are also reported in Figure 13, split by corresponding pair type. With the exception of the PV case when = 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 PV cross-correlation measurement when = 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 rarity threshold, see the details in Section 3.4). Remarkably, even in crosscorrelations we find an enhanced sensitivity of the inflection scales to massive neutrino effects.

Clustering of Critical Points in Massive Neutrino
Cosmologies: Redshift Evolution Effects Finally, we briefly address redshift evolution effects, in relation to the clustering properties of critical points. To Specifically, Figure 14 displays the clustering statistics of critical points in configuration space at = 1 (top panels), = 2 (middle panels), and = 3 (bottom panels), for a rarity choice corresponding to R = 20 and varying cosmologies. We adopt the same stylistic conventions as in Figure 3. In the various panels, we maintain the same − 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 non-zero neutrino mass causes an enhancement of the BAO peak amplitudes: the enhancement is more significant for higher neutrino mass values (see for example the = 0.4 eV case).
To appreciate possible redshift evolution effects more clearly, in Figure 15 we now show auto-correlation functions for the same rarity threshold previously examined and within a fixed cosmological model -while we vary the redshift in each panel by considering = 1 (purple), = 2 (dark green), and = 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.
Around the scale of BAO peaks, all of the auto-correlations show small evolution with redshift. And, as expected,evolution (although weak) is more relevant for the autocorrelations involving the most non-linear critical points Figure 11. BAO dip/peak amplitudes in massive neutrino cosmologies obtained from cross-correlation measurements, normalized by their correspondent values in a massless neutrino scenario and expressed in percentages as a function of rarity threshold at = 0 (with R = 5, 10, 20, respectively). Colors and symbols same as in Figure 5. In cross-correlations, departures from analogous measurements carried out in the baseline = 0.0 eV framework can be as high as ∼ 9% at = 0 when = 0.1 eV. (peaks, voids), while those involving less non-linear 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.

SUMMARY AND OUTLOOK
Determining the neutrino mass scale and type of hierarchy are primary targets for all of the ongoing and upcoming large-volume astronomical surveys, as well as a major goal of future space missions. Current limits from cosmology are already putting considerable pressure on the IH scenario as a plausible possibility for the neutrino mass ordering (Planck Collaboration et al. 2020a;eBOSS Collaboration et al. 2021), and a direct neutrino mass detection, or at least more competitive upper bounds on , are expected in the next few years from Stage-IV cosmological experiments. Moreover, gaining a deeper theoretical understanding of the impact of massive neutrinos on structure formation particularly at small scales and on the major observables typically used to characterize neutrino mass effects are necessary, to obtain robust constraints free from systematic biases.
Traditional methodologies for detecting and constraining neutrino mass effects from cosmology generally rely on a single observable, or on a limited scale-range. Examples include CMB gravitational lensing, the ISW effect in polarization maps, as well as a large number of LSS tracers such as 3D galaxy clustering via the matter power spectrum, bispectrum, cosmic shear, galaxy clusters, voids statistics, and much more.
In this work, we pursued a novel and different multiscale route, exploiting the remarkable topological properties of critical points -whose positions and heights at fixed smoothing scale retain precious cosmological information (Sousbie et al. 2008;Pogosyan et al. 2009;Gay et al. 2012;Cadiou et al. 2020;Shim et al. 2021;Kraljic et al. 2022). Besides being more robust to systematic effects, and in particular to non-linear evolution, clustering bias, and redshift space distortions, such special points are useful because they represent a meaningful and efficient compression of information of the 3D density field capturing its most significant features; thus, they provide a more fundamental and global view of the cosmic web as a whole. The overall weak sensitivity to systematics ultimately stems from their close relationship with the underlying topology, as described for example by the genus or Betti numbers. And, noticeably, all of the topological changes of a space occur only at critical points, and critical points of the density field are responsible for both the formation and destruction of a given feature.
For the first time, we characterized their clustering statistics in massive neutrino cosmologies, addressing the critical point sensitivity to small neutrino masses with the goal of identifying (possibly unique) neutrino multiscale signatures on the corresponding web constituents -i.e., halos, filaments, walls, and voids. We carried out our measurements on a subset of realizations from the QUIJOTE suite (Villaescusa-Navarro et al. 2020), using full snapshots at = 0, 1, 2, 3 for a choice of representative neutrino masses, as explained in Section 2. Exploiting a 'density threshold-based' methodology, described in Section 3, we computed critical point auto-and cross-correlation functions in configuration space for a series of rarity thresholds, and also characterized redshift evolution effects. Our major results, presented in Section 4 and summarized through Figures 3-15, are focused on BAO scales and are centered on two key aspects: (1) the multiscale effects of massive neutrinos on the BAO peak amplitudes of all of the critical point correlation functions above/below rarity threshold, and (2) the multiscale effects of massive neutrinos on the spatial positions of their correspondent correlation function inflection points at large scales.
The main findings and pivotal outcomes of this first work can be summarized as follows: • The BAO scale determined from critical point autoand cross-correlation measurements above/below rarity threshold is always robustly detected at 102.5ℎ −1 Mpc. The scale coincides with the BAO expected location in the reference cosmology, and it is independent of critical point type, neutrino mass, and rarity. 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 corresponding structures.
• All of the BAO peaks in auto-correlations are amplified with decreasing rarity: the amplitude of is higher for lower values of R, regardless of the critical point type.
Similarly, such amplifications of the BAO peaks are also found in critical point cross-correlations characterized by an identical overdensity sign (i.e., PF and VW, similarly biased tracers), and in the BAO dips of the critical point cross-correlations of overdense and underdense regions having opposite overdensity sign (i.e., PV, PW, VF , WF ). Cross-correlations of similarly biased tracers exhibit a behavior comparable to those of auto-correlations, while BAO features are 'reversed' and manifest as dips for oppositely biased tracers. Note also that the effect of a rarity threshold is also similar to that of smoothing.
• Massive neutrinos affect the BAO peak amplitudes of all of the critical point auto-and cross-correlation functions above/below rarity threshold. Differences at = 0 between a massless neutrino cosmology and a scenario with = 0.1 eV can reach up to ∼ 7% in autocorrelations, and ∼ 9% in cross-correlations. BAO peaks/dips become more pronounced (in absolute value terms) with increasing neutrino mass.
• Inflection scales are detected around 125ℎ −1 Mpc, both from auto-and cross-correlation measurements. Their exact values differ according to the critical point type. And, remarkably, inflection scales are altered by a nonzero neutrino mass -independently of rarity -with their spatial position increasing with augmented neutrino mass.  • Generally, extrema show noisier and less smoother autoand cross-correlation function shapes when compared to saddle-point statistics, as the time evolution of extrema is more nonlinear than that of saddles. Hence, saddle point statistics may be more advantageous for extracting cosmological information, and in relation to neutrino mass constraints.
• Around the BAO scale, all of the auto-correlations barely show evolution with redshift. Also, as expected, -evolution (although weak) is more relevant for autocorrelations involving the most nonlinear critical points (peaks, voids), while those involving less nonlinear critical points (filaments, walls) seem mostly insensitive to redshift evolution effects. And, at any given redshift that we considered and for all of the critical point types, a non-zero neutrino mass causes an enhancement of the BAO peak amplitudes, which appears more significant for higher neutrino mass values.
Our novel approach to neutrino mass effects, based on the statistics of critical points, offers a multiscale perspective, since such points carry remarkable topological properties and are faithful representation of their corresponding LSS. Moreover, critical points are also less sensitive to systematics, and this is among the reasons why our auto-and cross-correlation measurements appear more sensitive to neutrino signatures than analogous measurements performed on the total matter correlation function. In this view, part of our results can be seen as a multiscale generalization of the Peloso et al. (2015) findings about the impact of massive neutrinos on BAOs, and perhaps also as the multiscale extension of the LP of the spatial correlation function -which is also subject to non-zero neutrino mass effects (Parimbelli et al. 2021). In addition, our technique is complementary to more traditional methods, and can be combined with such methods to enhance the SNR: we will address this aspect in upcoming publications. Hence, our study is particularly relevant for ongoing and future largevolume redshift surveys such as DESI and the Rubin Observatory LSST, that will provide unique datasets suitable for establishing competitive neutrino mass constraints, as well as for future high-redshift probes like the WHT Extended Aperture Velocity Explorer (WEAVE; Dalton et al. 2012).
To this end, much work still remains to be performed in order to bring the methodologies presented here to real data applications: the current study is just the first of a series of investigations that aim at exploring the sensitivity of critical points and critical events to massive neutrinos, and more generally in relation to the dark sector. Ongoing efforts are along two major directions: towards galaxy clustering at relatively low redshift, and towards the high-universe as mapped for example by the Ly forest. In forthcoming publications, we will focus on these aspects, and also characterize the relation of our methodology with a persistence-based approach. Moreover, we will show how the combined multiscale effects presented here -inferred from the clustering of critical points -can be used to set upper limits on the summed neutrino mass and infer the type of hierarchy. This work is supported by the National Research Foundation of Korea (NRF) through grant No. 2020R1A2C1005655 funded by the Korean Ministry of Education, Science and Technology (MoEST) and by the faculty research fund of Sejong University. We acknowledge extensive usage of our com-puting resources (Xeon Silver 4114 master node, Xeon Gold 6126 computing nodes, and Lustre file system architecture) at Sejong University. G.R. would like to thank Simon White and Volker Springel for kind hospitality at the Max-Planck Institute of Astrophysics (MPA Garching) in Summer 2022, where this manuscript was finalized. We are also grateful to the referee for insightful feedback and suggestions. APPENDIX A. ANALYSIS TECHNICALITIES: ON BIN SIZE AND SMOOTHING In this section, we provide more details on the bin size used in the computations of the auto-and cross-correlation functions presented in the main text, and on the smoothing scale adopted. Specifically, the left panel of Figure 16 shows the effect of a varying bin size on the critical point auto-correlation function. As an example, we consider only maxima at = 0 for a rarity threshold R = 5, when grid = 1024 and the smoothing scale is fixed to G = 3ℎ −1 Mpc. Contrasting colors are used to indicate different bin sizes ( bin ), with bin = 1, 2, 3, 5, 10, respectively, in units of ℎ −1 Mpc. As evident from the plot, altering the bin size does not affect the overall shape of the auto-correlation function, including the location of the BAO peak and the zero-crossing scale. However, clearly smaller bin sizes manifest noisier shapes (i.e., bigger errorbars). Our choice of 5ℎ −1 Mpc, represented by the green line with errorbars in the panel, guarantees a sufficiently smooth : this is what is typically assumed in surveys like eBOSS for characterizing two-point correlation functions. Note also that we adopt a smaller bin size of 1ℎ −1 Mpc in selected regions near inflection points, employing a 'refinement' technique, in order to accurately determine their spatial locations at large scales. The right panel of Figure 16 shows instead the effect of smoothing on correlation function measurements. In detail, at a fixed grid = 1024 and for a bin size of 5ℎ −1 Mpc, we compute the two-point auto-correlations of maxima at = 0 when R = 5, for different choices of G . Contrasting colors represent alternative smoothing scale choices, as specified in the panel, with G = 1.5, 3.0, 4.5 and G expressed in ℎ −1 Mpc units. As pointed out in the main text, smoothing acts similarly to the effect of a rarity threshold choice -see also Kraljic et al. (2022). Namely, an increase in smoothing causes a decrement in the number of volume elements along with an increase in bias, and consequently the spatial correlation function becomes noisier because of enlarged statistical uncertainties. In essence, it is the combination of smoothing and rarity which is relevant: in fact, smoothing erases small-scale structures, and too rare events may be too noisy while less rare events will manifest less enhanced characteristic features. Generally, decreasing rarity or augmenting smoothing provides a more significant signal, and increases the statistical uncertainty. Hence, selecting an appropriate smoothing scale along with a suitable rarity is relevant for the 'density-threshold C P M N 19  based' approach, and the choice of rarity is eventually a compromise targeted to the specific problem addressed. In our specific case, considering three rarities R = 5, 10, 20 and setting G = 3.0ℎ −1 Mpc appears to be ideal, since the positions and heights of maxima are better constrained with smaller smoothing.

B. CRITICAL POINT ABUNDANCE: USEFUL TABLES
Complementing the information at = 0 reported in Tables 2 and 3 that appear in the main text (Section 4.1), we provide here additional details related to redshift evolution effects. Specifically, Table 4 contains the overall abundance of the entire set of critical points at = 1, 2, 3, classified by type, in the three cosmological frameworks examined in this work (i.e., baseline massless neutrino model, cosmology with = 0.1 eV, and scenario with = 0.4 eV, respectively). The associated errorbars are the corresponding 1 variations. Table 3 reports instead various ratios between the number of peaks (P), voids (V), filaments (F ), and walls (W) at the same redshifts considered in Table 4, showing clear departures from Gaussian expectations due to nonlinear evolution -in addition to the presence of a non-zero neutrino mass. Note also that the number of critical points decreases with time (or redshift), although the ratio of extrema over saddles remains constant.