DEMNUni: cross-correlating the nonlinear ISWRS effect with CMB-lensing and galaxies in the presence of massive neutrinos

We present an analytical modelling of the angular cross-correlations between the Integrated Sachs Wolfe-Rees Sciama (ISWRS) effect and large-scale structure tracers in the presence of massive neutrinos. Our method has been validated against large N-body simulations with a massive neutrino particle component, namely the DEMNUni suite. We investigate the impact of different neutrino masses on the cross-correlations between the ISWRS effect and both the galaxy clustering and the lensing of the Cosmic Microwave Background (CMB). We also test the ability of current nonlinear matter power spectrum modellings to reproduce neutrino effects in such cross-correlations. We show that the multipole position of a characteristic sign inversion in the cross-spectra, due to nonlinear effects, is strongly related to the total neutrino mass $M_\nu$ and depends almost linearly on it. While these nonlinear cross-correlation signals may not be able alone to constrain the neutrino mass, our approach paves the way to the detection of such cross-spectra on small scales for their exploitation in combination with main probes from future galaxy surveys and CMB experiments.


Introduction
Since its discovery in 1965 [1], the CMB has provided the most significant evidence for the Big Bang Theory and the study of its anisotropies has given us a unique window onto the primordial Universe.Measurements from WMAP [2][3][4] and Planck [5][6][7] have consistently been in line with the predictions of the ΛCDM model, the so-called Standard Cosmological Model, and have been fundamental in inferring constraints on the cosmological parameters that characterise it.However, there are still many open questions in cosmology, and several tensions exist in the Standard Cosmological Model.One of the most investigated topics today is the nature of Dark Energy (DE), which appears to dominate the energy density of the Universe at late times (see e.g.[8]), counteracting the gravitational collapse and inducing a background accelerated expansion.Another example is the value of the sum of the three neutrino masses, M ν = m ν , whose best bounds come from cosmological data [7,9,10].Among several cosmological probes, both these issues can be addressed also by measuring and analysing the late Integrated Sachs Wolfe (ISW) effect [11], which is the temperature variation that CMB photons experience when they propagate in a time-varying gravitational potential, and the Rees Sciama (RS) effect [12], which is its nonlinear counterpart.DE is mainly responsible for the combined ISWRS signal [13], but it is not the unique source of such effect.The presence of massive neutrinos induces a slow decay of the gravitational potential [14] that generates ISWRS even in the absence of a background expansion.Therefore, a full reconstruction of the ISWRS effect would provide new information on the physics of neutrinos and of DE, allowing to help constraining the DE Equation of State (EoS) and the total neutrino mass.However, in this respect there has not been much progress as this signal is extremely faint compared to CMB primary anisotropies.Moreover, these measurements are limited in the linear regime by cosmic variance and in the nonlinear regime by the lack of data.The only way to reconstruct the ISWRS is to use its cross-correlation with the Large Scale Structure (LSS), in order to exploit the tracers that follow the dark matter field and the variation of the gravitational potential produced by DE and massive neutrinos [15][16][17][18][19][20][21][22].
In this work, we focus on the cross-correlations of the ISWRS with the CMB-lensing (CMBL) potential and the galaxy distribution.As the RS effect anti-correlates with both of them [23][24][25], these cross-spectra are characterised by the presence of a sign inversion that occurs when passing from the linear (ISW) to the nonlinear (RS) regime.Moreover, the position of the sign inversion, which indicates the appearance of nonlinearities, shifts towards smaller cosmological scales due to the suppression of the matter power spectrum resulting from the presence of massive neutrinos [23,25].
We develop an analytical method to produce the ISWRS spectra using the nonlinear modelling of the matter power spectrum implemented in the Boltzmann solver code CAMB1 [26].Given the large redshift extension investigated in this work [27], we consider only two of the revised Halofit [28,29] models implemented: the Takahashi2012+Bird2014 (hereafter TB) model [30], extended to massive neutrino cosmologies by including Bird correction2 [31]; the Mead2020 (hereafter M20) model [32], which is an improvement over the Mead2016 model [33] concerning the treatment of Baryon Acoustic Oscillations (BAO) nonlinear damping.We then validate this analytical approach by comparing it with signals extracted from the "Dark Energy and Massive Neutrino Universe" (DEMNUni) N-body numerical simulations [25].
These studies are particularly relevant in light of ongoing and forthcoming galaxy surveys, such as Euclid [34][35][36] and SKA [37][38][39][40][41][42][43], and CMB surveys, such as CMB-S4 [44,45], which aim to achieve measurements of small scales with high precision and accuracy, such that the possibility of detecting the RS effect could become concrete [46].This paper is organised as follows.In Section 2 we recap the theoretical framework of the cross-correlations of the ISWRS signal with the CMBL potential and the galaxy distribution in the presence of massive neutrinos.In Section 3 we present how we model ISWRS nonlinearities using measurements from the DEMNUni simulations.In Section 4 we test the analytical method for the cross-correlation between the ISWRS signal and the CMBL potential against the DEMNUni simulated cross-spectrum.In Section 5 we validate the cross-correlation between the ISWRS signal and the galaxy distribution against the cross-spectrum measured from the DEMNUni mock maps.Moreover, we model the ISWRS-galaxy cross-spectrum for different redshift ranges.In Section 6, we evaluate the predictability of M ν via the estimation of the sign inversion position in both the analysed cross-correlations.Finally, in Section 7 we draw our conclusions.

The Integrated Sachs Wolfe and the Rees Sciama effect
When CMB photons pass through a time-varying gravitational potential well along their path from the last scattering surface to us, they undergo a variation in energy that translates into a variation in temperature known as the Integrated Sachs Wolfe effect: where n is a unit direction vector on the sphere, T 0 is today (t 0 ) CMB temperature, t ls is the time at the last scattering surface, χ the comoving distance and Φ is the time derivative of the gravitational potential.This effect is mainly due to the presence of DE that induces an accelerated background expansion of the Universe and counteracts the gravitational potential Φ inducing a not vanishing Φ [20,[47][48][49].The nonlinear growth of density perturbations produces additional temperature perturbations, leading to the RS effect.The RS contribution to temperature fluctuations is much more slow over time than the late ISW one, which means that at high redshifts RS predominates over the late ISW [50].Consequently, we are able to observe Φ ̸ = 0 even when the DE effect becomes negligible.Moreover, the presence of massive neutrinos induces a not vanishing derivative of the gravitational potential even during matter domination era.This is because, after becoming non-relativistic, neutrinos free-stream with large thermal velocities that suppress the growth of neutrino density perturbations on scales smaller than the so-called "free-streaming length" [9,21]: where m ν is the mass of the single neutrino species (i.e.ν e , ν µ or ν τ ), H(z) is the Hubble parameter as a function of the redshift z and H 0 ≡ H(z = 0) is the Hubble constant.Moreover, because of the gravitational back reaction effects, the evolution of cold dark matter (CDM) and baryon densities is affected as well by the presence of neutrinos, and the total matter power spectrum is suppressed at scales λ ≪ λ FS [51].Consequently, on small cosmological scales the free-streaming of neutrinos induces a slow decay of the gravitational potential, acting during both the matter and the DE dominated eras, and this effect depends on the total neutrino mass M ν = m ν .The more neutrinos are massive, the more the matter power spectrum will be suppressed with respect to the massless neutrino case, and therefore nonlinearities will appear on smaller scales.
One of the effects of massive neutrinos, which has been poorly investigated so far, is the shift they induce in the sign inversion that characterises the cross-power spectra between the ISWRS signal and the CMBL potential (P ΦΦ ), and between the ISWRS signal and the galaxy distribution (P Φδ ).These sign inversions appear because of the anti-correlation between the RS and both the CMBL potential and the galaxy distribution, that depend on the gravitational potential Φ.While on large scales (i.e. the linear regime) the gravitational potential decays because of the Universe expansion (Φ < 0, Φ < 0), on small scales (i.e. the nonlinear regime) Φ grows because of nonlinear structures formation (Φ < 0, Φ > 0).Consequently, the time derivative of the gravitational potential associated to the ISWRS signal is negative in the linear regime and positive in the nonlinear one.The net result is that ⟨Φ Φ⟩ > 0 in the late ISW regime, producing a positive correlation, and ⟨Φ Φ⟩ < 0 in the RS regime, producing an anti-correlation [24,52].
Following the example of [23,53,54], we report below the dimensionless and scaled (i.e.divided by 3/2[H 0 /(ck)] 2 Ω m ) forms of P ΦΦ and P Φδ in Limber approximation [21,54]: , with where a(z) is the scale factor, P δδ is the matter power spectrum, k is the wavenumber, and with Ω m the matter density parameter and c the speed of light.
We compute each term of Equations (2.3) and (2.4) with CAMB and use both the nonlinear modellings chosen for the matter power spectrum (i.e.TB and M20).The absolute values of the results are shown in the panels of the first and third rows of Figure 1, for both the TB and M20 nonlinear modellings, at different redshift values, for the four total neutrino masses considered in this work.
The panels in the first and third rows of Figure 1 show how the presence of nonlinearities varies with time.As already mentioned, the more we travel back in time, the more the RS effect will dominate over the late ISW.Considering that at large redshifts the RS traces not only nonlinearities at small scales, but even the production of filamentary structures (visible at scales of few hundreds Mpc) [50], we see the shifts of the sign inversion moving towards smaller k (i.e.larger cosmological scales) from z = 0.5 to z = 1.89.Following [54], we verified that, in the redshift range of our interest, the linear growth rate of matter perturbations, f ≡ d ln(D)/d ln(a) (where D is the so-called linear growth factor) is f > 0.5, which is the condition expected to be necessary for the presence of the sign inversions measured in the cross power spectra.
In the panels of the second and fourth rows in Figure 1, we present the trend of k inv as a function of M ν for both the nonlinear models tested.We observe the expected effect: the more massive neutrinos are, the more nonlinearities (and therefore, the sign inversion) move toward smaller cosmological scales (i.e.larger k).The only exception, for both ∆ 2 ΦΦ and ∆ 2 Φδ , is the trend obtained for M ν = 0.17 eV when using the TB model at z = 0.5, where the power spectra show the sign inversion at scales slightly larger than the massless neutrino case.We additionally include the results for M ν values slightly smaller and larger than 0.17 eV.On the one hand, we observe an increase of k inv from M ν = 0 eV to M ν = 0.15, 0.16 eV, on the other hand, we confirm a decrease corresponding to M ν = 0.17 eV, followed by a new increase for M ν = 0.18 eV.This effect is probably due to inaccuracies in the TB nonlinear modelling.
The similarity of the results shown in the panels of the second and fourth rows is due to the fact that the position of k inv as a function of M ν strictly depends on the time derivative of P δδ , that appears and is the same in both Equations (2.3)- (2.4).All the other terms of these two equations are responsible for the differences in the amplitude alone of the analysed cross-spectra, as it is visible in the panels of the first and third rows.
These results highlight the different behaviour of the two models.At small redshifts the TB model predicts an higher level of nonlinearities than the M20 model.Indeed, in the panels of Figure 1 corresponding to z = 0.5, the TB k inv are smaller than the M20 ones.However, it

TB M20
Figure 1: First row : Absolute values of the dimensionless and scaled form of the ISWRS-CMBL potential cross-spectra as a function of the wavenumber k at z = 0.5, 0.9, 1.4, 1.89 for the four neutrino total masses considered (M ν = 0, 0.17, 0.30, 0.53 eV as red, blue, green and purple lines, respectively), obtained via the TB model (solid lines) and the M20 model (dotted lines).Second row : Measured k inv as a function of M ν for the TB (diamonds) and the M20 (squares) models.Third row : Absolute values of the dimensionless and scaled form of the ISWRS-galaxy distribution cross-spectra as a function of the wavenumber k at z = 0.5, 0.9, 1.4, 1.89 for the four neutrino masses considered (M ν = 0, 0.17, 0.30, 0.53 eV as red, blue, green and purple lines, respectively), obtained via the TB model (solid lines) and the M20 model (dotted lines).Fourth row : Measured k inv as a function of M ν for the TB (diamonds) and the M20 (squares) models.
is easy to notice that at z = 0.9 and z = 1.4 there is a switch of trend between TB and M20, and at z = 1.89 the M20 model predicts an higher level of nonlinearities with respect to TB.Such behaviours strictly depend on the approximations adopted for the implementation of the two Halofit models tested, which however do not aim at reproducing correctly the time derivative of the gravitational potential sourcing the ISWRS effect.The presence of massive neutrinos is expected to affect similarly also the angular powerspectra of both these cross-correlations, that actually depend on the time derivative of the matter power spectrum P δδ , and can be computed via [55]: where z ls = 1100, and [21,46,54]: where n(z) and b(z) are the galaxy selection function and the galaxy bias, respectively.The integral limits used in this work are: z min = 0.02 and z max = 1.89.This is because we use Equations (2.5) and (2.6) to compute the theoretical predictions that are compared against the cross-spectra measured from the DEMNUni maps, that cover the redshift range [0.02, 1.89] (see Section 3).
Apart from n(z) and b(z), all the ingredients of Equations (2.5) and (2.6) have been obtained via CAMB.

Nonlinear modelling with N-body simulations
The choice to investigate the nonlinear regime requires the use of large N-body simulations for validation tests.For this work, we use the "Dark Energy and Massive Neutrino Universe" (DEMNUni) [25] N-body cosmological simulations to model at small scales the ISWRS signal and their cross-correlations with the CMBL and the galaxy distributions.

The Dark Energy and Massive Neutrino Universe simulations
The DEMNUni simulations have been produced with the aim of investigating large-scale structures in the presence of massive neutrinos and dynamical dark energy (DDE), and they were conceived for the nonlinear analysis and modelling of different probes, including dark matter, halo, and galaxy clustering [27,[56][57][58][59][60][61][62][63][64][65], weak lensing, CMBL, SZ and late ISW effects [25,[66][67][68], cosmic void statistics [69][70][71][72][73], and cross-correlations among these probes [74].To this end, they combine a good mass resolution with a large volume to include perturbations both at large and small scales.In fact, these simulations follow the evolution of 2048 3 cold dark matter (CDM) and, when present, 2048 3 neutrino particles in a box of side L = 2 h −1 Gpc.The fundamental frequency of the comoving particle snapshot is, therefore, k F ≈ 3 × 10 −3 h Mpc −1 , while the chosen softening length is 20 kpc/h.The simulations are initialised at z ini = 99 with Zel'dovich initial conditions.The initial power spectrum is rescaled to the initial redshift via the rescaling method developed in Ref. [75].Initial conditions are then generated with a modified version of the N-GenIC software, assuming Rayleigh random amplitudes and uniform random phases.The DEMNUni simulations have been performed using the tree particle mesh-smoothed particle hydrodynamics (TreePM-SPH) code GADGET-3, an improved version of the code described in [76], specifically modified in [77] to account for the presence of massive neutrinos.This version of GADGET-3 follows the evolution of CDM and neutrino particles, treating them as two distinct collisionless fluids.Given the large amount of memory required by the simulations, baryon physics is not included.This choice, however, does not affect the results of this work [31,56,78].The DEMNUni suite is composed by a total of sixteen simulations accounting for different combinations of the total neutrino mass M ν and the DE EoS.In this project, we have worked with four of them, characterised by a baseline Planck 2013 cosmology [5], namely a flat ΛCDM model with τ = 0.0925, n s = 0.96 and A s = 2.13 × 10 −9 , generalised to a νΛCDM by varying only the sum of the neutrino masses over the values M ν = 0, 0.17, 0.30, 0.53 eV (whence the corresponding values of Ω ν and Ω cdm , keeping fixed Ω m = 0.32 and Ω b = 0.05), considering a degenerate neutrino mass hierarchy.As aforementioned, in this work we use maps that cover a range of redshift that goes from z min = 0.02 to z max = 1.89.For each simulation, 63 outputs have been produced, logarithmically equispaced in the scale factor a = 1/(1 + z), in the redshift interval z = 0 − 99, 49 of which lay between z = 0 and z = 10.For each of the 63 output times, it has been produced on-the-fly one particle snapshot, composed of both CDM and neutrino particles, one 3D grid of the gravitational potential, Φ, and one 3D grid of its time derivative, Φ, with a mesh of dimension 40963 that covers a comoving volume of (2 h −1 Gpc) 3 .
The map-making procedure has been developed taking into account the standard pixelisation approach introduced by the Hierarchical Equal-Area isoLatitude Pixelization [79] (HEALPix 3 ).

Map-making procedure: ISWRS and CMB-lensing
The procedure followed to build the ISWRS and CMBL potential maps is an adaptation of the one developed in [80].To make the ISWRS maps (left column of Figure 2), CMB photons are ray-traced along the undeflected line of sight through the 3D field Φ.A similar ray-tracing has been applied to the 3D Φ-grids, in order to produce the CMBL potential maps (right column of Figure 2) from the same realisation of the Universe adopted for the ISWRS maps.To this aim (i.e.producing exactly the same realization of the Universe), the simulated Φ and Φ-grids have been stacked around the observer, located at z = 0, applying the replication and randomisation of the DEMNUni comoving outputs along the line of sight, while keeping periodic boundary conditions in the transverse direction, following [80].This particular 3D tessellation scheme is required to avoid both the repetition of the same structures along the line of sight and the generation of artifacts, like ripples, in the simulated deflection-angle field.The latter can be avoided only if the peculiar gravitational potential is continuous transversely to each line of sight.Comparing the left and right columns of Figure 2, it is possible to appreciate at a glance how the ISWRS and the CMBL maps are highly correlated.Comparing the scales of each rows, the slight difference, due to the different M ν values, emerges.

Simulated ISWRS and CMB-lensing auto-spectra
The power spectra extracted from the maps of Figure 2 are reported in Figure 3.For comparison, the orange squares and the orange dashed lines represent theoretical predictions in the linear regime obtained with CAMB for the massless neutrino case, in order to highlight the differences from the nonlinear modelling.The subpanels represent the percentage relative differences of the νΛCDM models with respect to the massless one.There is a good agreement between the simulated ISWRS signals and the ones predicted with CAMB for ℓ < 100 (which roughly corresponds to the transition between the linear and the nonlinear regimes [81,82]).Past the linear limit, CAMB predictions fail, given that the RS contribution is not implemented in the code.The CMBL potential linear prediction for the massless neutrino case, similarly, perfectly agrees with DEMNUni mocks up to the linear limit and then shows a lack of power with respect to the nonlinear ones.Focusing on spectra extracted from maps, it is easy to appreciate at small scales the differences produced in the ISWRS induced temperature anisotropies by free-streaming neutrinos with different total masses.
Similar effects are noticeable also in the CMBL potential power spectra, considering the percentage relative differences with respect to the massless neutrino case (see subpanels of Figure 3).In particular, as expected, the effect of massive neutrinos produces a suppression of power in the CMBL potential which increases with the neutrino mass and decreases with the angular scale (at low ℓ).This effect is reasonably due to the nonlinear excess of suppression of the total matter power spectra with respect to linear expectations in the presence of massive neutrinos, that counteract the nonlinear evolution [25].

C (%)
Figure 3: Left: Angular power spectra of the total ISWRS induced temperature anisotropies computed from DEMNUni maps with M ν = 0, 0.17, 0.30, 0.53 eV (solid red, dashed blue, dot-dashed green, dotted purple lines, respectively) compared with linear prediction (orange squares).Right: Same for angular power spectra of the CMBL potential.In both subpanels, the percentage relative differences with respect to the massless neutrino case is reported.

Map-making procedure: projected galaxy mocks
The procedure followed to construct the galaxy maps is an adaptation of the one developed in [67,83] for convergence maps.The convergence maps are extracted with a postprocessing procedure acting on the N-body particle snapshots to create matter particle full-sky lightcones [84,85].The standard way to build lightcones is to pile up high-resolution comoving snapshots within concentric cells to fill the lightcone to the maximum desired source redshift.With the observer placed at z = 0, the volume of the lightcone is sliced into full-sky spherical shells of the desired thickness in redshift.All the particles distributed within each of these 3D matter shells are then projected onto 2D spherical maps, assigning a specific sky pixel to each particle, via HEALPix pixelisation procedure.Then integration along the undeflected line-of-sight is applied to obtain the convergence maps accounting for the lensing efficiency in the desired redshift range.A very similar method has been applied to the DEMNUni galaxies to construct projected galaxies mocks [27] (see Figure 4).Specifically, in order to reproduce the galaxy distribution in the DEMNUni simulations, the SubHalo Abundance Matching DEMNUni M = 0 eV Galaxies map -1.00 5.89 DEMNUni M = 0.17 eV Galaxies map -1.00 5.44 DEMNUni M = 0.30 eV Galaxies map -1.00 5.44 DEMNUni M = 0.53 eV Galaxies map -1.00 5.44 Figure 4: Full-sky galaxies maps obtained for the four cosmologies used in this work: ΛCDM and νΛCDM with M ν = 0.17, 0.30, 0.53 eV.These maps have a pixel number and resolution given by N side = 2048 and a RING ordering, following the procedure described in Section 3.3.
(SHAM) technique [65] has been used.The SHAM method assumes a one-to-one relation between a physical property of a dark matter halo/subhalo and an observational property of the galaxy that it hosts.In particular, it assumes that the most massive galaxies form in the dark matter haloes characterised by the deepest potential wells.Galaxies are linked to the corresponding dark matter structures using stellar mass or luminosity as galaxy property and a measure of halo mass or the circular velocity, as halo property (i.e. a proxy of the depth of the local potential well).Here we link the stellar mass to the halo mass using the Moster et al. (2010) [86] parameterised stellar-to-halo mass relation, fitted against the DEMNUni simulations.As the SHAM galaxy catalogues are fitted against observations, they guarantee an high-degree of agreement with real data from galaxy surveys [87,88].We do not expect the particular galaxy modelling, and in particular here the SHAM parameters, to significantly affect the sign inversion position or to improve the sensitivity to M ν .In fact, we will show in the following Sections that the cross-correlation between ISWRS and CMBL and the one between ISWRS and galaxy clustering share a very similar trend as far as the l inv is concerned, even if the CMBL does not involve any galaxy bias.Moreover, the SHAM parameters may at most affect the galaxy bias function and its scale-dependence, and we will show in Section 5 that the sign inversion position is not very sensitive to the particular bias modelling but rather to the bias amplitude, i.e. the galaxy population considered and its evolution with z.
Applying to the SHAM galaxy comoving snapshots the same stacking technique as for particle comoving snapshots, we have generated full-sky mock galaxy catalogues from the DEMNUni simulations in different cosmological models.

Galaxy bias and selection function from DEMNUni simulations
The galaxy bias and selection function necessary to compute the analytical predictions of the cross-spectra between ISWRS and galaxies (see Equation (2.6)) have been directly extracted from the DEMNUni simulations.This guarantees a more accurate comparison between mocks and analytical predictions.

Galaxy selection function
The galaxy selection function extracted from the DEMNUni ΛCDM map can be modelled via the following equation: where A = 1.89,B = 1.95,C = 1.41 and z 0 = 0.88.This is the n(z) used to compute all the theoretical predictions to be compared with DEM-NUni data, also for massive neutrino cases.This is because galaxy counts vary mostly because of the variation of the simulated cosmological volume, which depends on H(z) via Ω m -that has been fixed to 0.32 for all the simulations (and consequently for all the theoretical predictions) -and via the DE EoS, whose presence has not been taken into account in this work.As a result, the normalised galaxy counts do not undergo any significant change due to the variation of M ν , and consequently it is possible to exploit the same n(z) to compute the theoretical predictions even in the massive neutrino case.We have verified this by measuring the galaxy counts from the DEMNUni catalogues with different values of M ν , obtained via the SHAM technique applied to COSMOS2020 and SDSS observations [65].2) with z m = 0.9, which is the one assumed for the photometric Euclid survey; the dot-dashed green and the dotted red lines are n(z, z m ) from Equation (3.2) with z m = 2, 3, respectively [21].Right: Scale-independent galaxy bias functions: the dashed black line is the Tutusaus bias [89]; the solid red, dashed blue, dot-dashed green and dotted purple lines are b gg (z, M ν ) in Equation (3.4) for M ν = 0, 0.17, 0.30, 0.53 eV, respectively; finally the dotted pink line is b gg (z) = √ 1 + z used in [35].The shaded grey area highlights the extrapolation of the galaxy bias function above the redshifts tested against the Euclid Flagship 1 simulation [89].
In Figure 5, together with the DEMNUni n(z), we show the galaxy selection function adopted in [21]: where α = 2, β = 1.5, z 0 = z m / √ 2 and z m = 0.9, 2, 3.The case z m = 0.9 corresponds to the galaxy selection function expected for the photometric Euclid survey [34,35].Equation (3.2) (with z m = 0.9, 2, 3) has been used in Section 5.2 to produce forecasts for z ≥ 2, i.e. for redshifts larger than the ones investigated in the DEMNUni maps.

Scale-independent bias
The galaxy scale-independent bias has been computed, for each M ν value, as the linear interpolation in 0 < z < 2 of the bias measurements obtained as the square root of the ratio between the galaxy auto power spectrum, P gg (k, z), and the matter power spectrum, P mm (k, z), from the DEMNUni galaxy and dark matter comoving snaphots, respectively: Here the average has been computed on large linear scales, specifically in the range 0.05 < k < 0.12 Mpc −1 .The linear interpolation in z of Equation (3.3) for M ν = 0, 0.17, 0.30, 0.53 eV, that from now on we call b gg (z, M ν ), are reported in the left panel of Figure 6.We use the measurements in Equation (3.3) to fit, for the different values of M ν , the coefficients in the following functional form adopted to produce the forecasts in Section 5.2 at z > 2: where the coefficients A, B, C and D depend on the neutrino mass, and the specific values A = 1.0,B = 2.5, C = 2.8, D = 1.6 correspond to the galaxy bias model from the Euclid Flagship 1 simulation, reported in [89], and shown in the right panel of Figure 5.We call Equation (3.4), for the different values of M ν , the "a la Tutusaus" bias.When we compare our theoretical predictions with the DEMNUni data, we consider also, for the massless neutrino case alone, the galaxy scale-independent bias computed as the linear interpolation in 0 < z < 2 of the bias measurements obtained as the ratio between the galaxy-matter cross-spectrum, P gm (k, z), and the matter power spectrum, P mm (k, z), from the DEMNUni galaxy and dark matter comoving snaphots: Here the average has been computed again over the range 0.05 < k < 0.12 Mpc −1 and its linear interpolation is reported in the left panel of Figure 6.Since it is affected by a lower level of shot noise, b gm (z, M ν = 0) should allow us to describe better the connection between matter and galaxies, and consequently to improve the cross-correlation reconstruction.However, the results obtained do not differ significantly from those obtained using b gg (z, M ν = 0), as we will show in Section 5.1.2.
Scale-dependent bias A scale-dependent bias fit has been obtained from the DEMNUni massless neutrino mocks using the functional form below [90]: where b 0 is the scale-independent bias, b 2 is a correction due to the peak theory and b 1 does not have a physical meaning, but it is necessary to improve the fit.These parameters have been obtained fitting Equation (3.6) up to k = 1 Mpc −1 against the DEMNUni scale-dependent bias measurements.In Figure 6, we show the evolution of the b i parameters in the left panel, while on the right we show b gg (k, z) as a function of z at some chosen wavenumbers in the range of scales used to measure the scale-averaged b gg (z, M ν = 0).

Simulated projected galaxy auto-spectra
In Figure 7, we report the spectra extracted from the galaxy maps of Figure 4 after subtracting the shot-noise value.
In addition to the DEMNUni maps spectra, we report the theoretical predictions of the C gg ℓ computed in both the linear and nonlinear regimes, for the massless neutrino case, via [91]: where: and for n(z) and b(z, M ν = 0) we use Equation (3.1) and Equation (3.3), respectively.Both linear and nonlinear predictions do not perfectly agree with simulated spectra on large scales.This is due to the limited dimension of the simulation box, which does not allow to represent those scales.Conversely, on small scales there is a good agreement with the nonlinear prediction and the expected disagreement with the linear one.
Despite the presence of shot-noise in the DEMNUni maps, its effect is negligible when cross-correlating these maps with the ISWRS maps (since the noises of the two fields do not correlate).Furthermore, this analysis does not aim to reveal a perfect match between the simulations and the analytical predictions also because of inaccuracies in the adopted nonlinear modelling.C GG (%) Figure 7: Galaxy angular auto-spectra extracted from the DEMNUni maps for M ν = 0, 0.17, 0.30, 0.53 eV (solid red, dashed blue, dot-dashed green, dotted purple lines, respectively) with shot-noise (dotted grey line) subtracted, compared with theoretical predictions in the linear (dashed orange line) and M20 nonlinear (dashed black line) regimes.In the subpanel, the percentage relative difference with respect to the massless neutrino case is reported.
The subpanel in Figure 7 represents the percentage relative differences of the νΛCDM models with respect to the massless neutrino one.Here the larger M ν is, the smaller the difference with respect to the massless neutrino case is.This is because the matter perturbation suppression due to more massive neutrinos is compensated by a stronger clustering [92].

The ISWRS cross-correlation with the CMB-lensing potential
The presence of the sign inversion at large multipoles and its shift due to the variation of M ν in the cross-power spectrum between the ISWRS effect and the CMBL was already verified in [25], where they use the DEMNUni maps obtained via ray-tracing from z min = 0.02 to z max = 20.88.Here, we verify the presence of the same effects in the cross-spectra of DEMNUni maps obtained up to z max = 1.89 and use these results to validate the analytical predictions computed using Equation (2.5).The cross-correlations of the DEMNUni ISWRS and CMBL potential maps for the four neutrino masses are shown in the top panels of Figure 8, compared with the theoretical nonlinear prediction computed for the TB (top left panel of Figure 8) and the M20 (top right panel of Figure 8) models.For completeness, also the theoretical cross-correlation predictions for the massless neutrino case in the linear regime have been reported.
The agreement between the simulation data and both the theoretical predictions obtained with CAMB up to ℓ ∼ 100 is quite impressive.Once exceeded ℓ ∼ 100, the anti-correlation between the RS effect and the CMBL potential takes over, resulting in the characteristic sign inversion that appears in both simulations (as already shown in [25]) and nonlinear predictions.Moreover, both of them validate the expected result, because the more neutrinos are massive, the more they counteract nonlinear evolution (i.e.structure formation) and nonlinear effects appear at smaller cosmological scales.In other words, the more neutrinos are massive, the more the sign inversion shifts towards larger multipoles.
The achieved results disprove the counter-intuitive effect reported in [93], where for larger M ν values the sign inversion position moves towards smaller multipoles.The behaviour found in [93] is presumably due to the use of a very approximated implementation of the neutrino mass dependence in the Halofit models, which leads to different results from ours. Figure 8: ISWRS-CMBL angular cross-spectra for M ν = 0, 0.17, 0.30, 0.53 eV (solid red, blue, green and purple lines, respectively) from the DEMNUni maps compared to the theoretical predictions in the linear (dashed orange line) and nonlinear (dotted lines) regimes, obtained using the TB model (left) and the M20 model (right), respectively.The subpanel shows the behaviour of ℓ inv as a function of M ν for the two nonlinear models, compared with the measurements from the DEMNUni simulations.
In the subpanel of Figure 8, we focus on the sign inversion position as a function of M ν , comparing the two Halofit models to DEMNUni.The 1σ error bars, shown for the DEMNUni measurements, have been estimated as follows: we generate 5000 Gaussian realisations of the DEMNUni maps and estimate the errors as the standard deviation.Further discussion will be provided in Section 6.

The ISWRS cross-correlation with galaxies
The cross-spectrum between the ISWRS signal and the galaxy distribution is the most used tool to investigate the ISWRS, because of its largest signal-to-noise ratio with respect to the spectra analysed in the previous Section.It is therefore the main focus of this work, which consists of verifying the presence of the sign inversion due to the anti-correlation between the RS effect and the galaxy distribution in the cross-correlation angular power spectra, showing that the position of this sign inversion depends on M ν , as in the case of the cross-correlation between the ISWRS signal and the CMBL potential [25].The validation of these results via the DEMNUni maps represents an improvement with respect to previous works: on the one hand, [21] computed this cross-correlation signal in massless and massive neutrino models, focusing only on the effects that neutrinos have in the linear regime; on the other hand, in the nonlinear regime, [54] focused only on the massless neutrino case, finding in the crossspectrum the sign inversion due to the anti-correlation between the RS effect and the galaxy distribution.
The cross-correlations of the DEMNUni ISWRS and the galaxy distribution maps for the four neutrino masses are shown in Figure 9, compared with the theoretical nonlinear predictions, computed with Equation (2.6).For completeness, also the theoretical predictions for the massless neutrino case in the linear regime have been reported.Here, for a fair comparison against the DEMNUni data, the theoretical predictions have been computed using n(z) in Equation (3.1) and the M ν -dependent linear interpolation in z of the measurements in Equation (3.3).
Also in this case, because of the suppression massive neutrinos induce in the matter power spectrum, the larger M ν , the more nonlinearities appear on smaller scales and consequently the sign inversion position that characterises this cross-correlation is shifted towards larger multipoles.
Moreover, the comparison of the simulated spectra with the theoretical predictions for both linear and nonlinear regimes is again quite impressive at large cosmological scales, while for ℓ ≳ 100 nonlinear spectra dramatically differ from the linear one, because of the presence of the sign inversion.This result validates again the analytical method against mocks and makes it an extremely powerful tool to predict cross-spectra of future surveys in both the linear and nonlinear regimes.
In the subpanel of Figure 9, we focus on the sign inversion position, ℓ inv , as a function of M ν , comparing the two Halofit models to DEMNUni.The 1σ error bars, associated to the DEMNUni measurements, have been estimated as described in the previous Section.Further discussion will be provided in Section 6.

Nonlinear ISWRS-galaxy cross-spectra: effects of different bias modellings
As already mentioned, differently from the galaxy selection function, the galaxy bias strongly depends on the neutrino mass.We perform several tests to evaluate how the use of different bias functions affects the comparison between predictions and the DEMNUni measurements.First, we produce theoretical predictions for the massive neutrino cases using both b gg (z, M ν > 0) and b gg (z, M ν = 0), and we verify that the more the neutrinos are massive, the more b gg (z, M ν = 0) fails in predicting the spectra (further details in Section 5.1.1).Second, we Figure 9: ISWRS-galaxy cross-spectra for M ν = 0, 0.17, 0.30, 0.53 eV (solid red, blue, green and purple lines, respectively) compared with theoretical predictions in the linear (dashed orange line) and nonlinear (dotted lines) regimes, obtained using the TB model (left) and the M20 model (right), respectively.The subpanel shows the behaviour of ℓ inv as a function of M ν for the two nonlinear models, compared with the measurements from the DEMNUni simulations.
check, for the massless neutrino case alone, the impact of using a scale-dependent bias function for the theoretical predictions (further details in Section 5.1.2).The results of these tests are shown in Figure 10 and Figure 11, respectively.

Massive neutrino impact on the scale-independent bias
First, we test how the analytical predictions for massive neutrino cosmologies compare with mock data when using the b gg (z, M ν > 0) bias function, or the b gg (z, M ν = 0).Figure 10 shows such a comparison for the cross-correlation signals, for which the results for ℓ ≤ 10 can be neglected because of finite volume effects in the mocks on those scales.
For the M ν = 0.17 eV case, differently from the other cosmologies, there is a small difference in the amplitude of the analytical cross-correlation power spectra, with respect to the DEM-NUni ones, at very small (ℓ < 30) and very large (ℓ > 1000) multipoles.These differences were already visible in Figure 9, and are evident in the top panels of Figure 10.They induce a difference between the theoretical predictions and the DEMNUni cross-spectrum larger than 35% for ℓ < 30 and ℓ > 1000.
Before the sign inversion, for the TB modelling the differences do not overcome 25% for b gg (z, M ν = 0.17) and 15% for b gg (z, M ν = 0).
Figure 10: Sensitivity to the M ν dependence of the galaxy bias.Left: Comparison between the simulated ISWRS-galaxy cross-correlations for M ν = 0.17, 0.30, 0.53 eV (solid black line) and theoretical predictions in the linear (dashed orange line) and nonlinear regimes, using the TB model combined with DEMNUni measurements of b gg in the massless neutrino (dashed red) and massive neutrino cases for M ν = 0.17 eV (top, dot-dashed blue), M ν = 0.30 eV (middle, dot-dashed green) and M ν = 0.53 eV (bottom, dot-dashed purple), respectively.The middle subpanel of each plot shows a zoom-in of the sign inversion region, while the bottom one shows the absolute value of the percentage relative differences of the predictions with respect to the simulated measurements.Right: Same as the left panel for the M20 model.
In the case of the M20 modelling, the difference between the predictions and measurements is slightly smaller (< 20%) when using b gg (z, M ν = 0.17), and ∼ 15% for b gg (z, M ν = 0).The differences in the sign inversion positions are: ∆ℓ DEMNUni−(TB) ∼ −18 and ∆ℓ DEMNUni−M20 ∼ +18.Concerning the signal amplitude, the predictions obtained with the two bias recipes show a difference at a level of ∼ 5%, while the sign inversions occur at the same multipole in the M20 model, and with a ∆ℓ = 1 for the TB one.This result is somehow expected since the two scale-independent biases change only the overall amplitude of the signal but not its shape.Instead, the fact that b gg (z, M ν = 0) in the theoretical predictions leads to smaller differences with respect to the DEMNUni spectra can be explained considering the inaccuracies of the TB and M20 models in reproducing the M ν = 0.17 eV case, as already shown in Section 2.
In the M ν = 0.30 eV case, shown in middle panels of Figure 10, both the nonlinear predictions, obtained with the M20 and TB models, differ from the DEMNUni spectra by no more than 15% when using b gg (z, M ν = 0.30), and no more than 25% when using b gg (z, M ν = 0).The differences in the sign inversion positions are: ∆ℓ DEMNUni−(TB) ∼ −10 and ∆ℓ DEMNUni−M20 ∼ +20.As expected, the difference between the two bias cases increases with the neutrino mass, specifically it is larger by about 10% than in the M ν = 0.17 eV case.Concerning the sign inversions corresponding to the two biases, they occur at the same multipole in the M20 model, and with ∆ℓ = 1 for the TB one.
Finally, as shown in bottom panels of Figure 10, in agreement with the trend already noticed for M ν = 0.17, 0.30 eV, for M ν = 0.53 eV the difference between the predictions using b gg (z, M ν = 0.53) or b gg (z, M ν = 0) increases up to ∼ 20%.The difference with respect to the DEMNUni spectra is ∼ 10% for b gg (z, M ν = 0.53), and ∼ 25% for b gg (z, M ν = 0).Even in this case, the sign inversions occur at the same multipole when using M20, and with ∆ℓ = 2 when using TB.The differences in the sign inversion positions, with respect to the DEMNUni spectra, are: ∆ℓ DEMNUni−(TB) ∼ −10 and ∆ℓ DEMNUni−M20 ∼ +10.

5.1.2
Scale-independent and scale-dependent bias effects in the massless neutrino case In Figure 11 we compare the effects of b mg (z), b gg (k, z) and b gg (z) in the modelling of the theoretical nonlinear predictions.Conversely to the work of [54], we find that the amplitude of the predictions obtained with the scale-dependent bias function does not show any particular difference with respect to the predictions computed with the scale-independent ones.The three reconstruction differ from DEMNUni spectra by no more than ∼ 10%.
With regard to the sign inversion position, we observe in Figure 11 that, when using the scale-independent prescriptions, independently of the TB and M20 models, the sign inversions occur at mostly the same multipoles as for the scale-independent case for which ∆ℓ DEMNUni−TB ∼ −20 and ∆ℓ DEMNUni−M20 ∼ +15, respectively.In fact, the ℓ inv obtained with b gg (k, z) differs by ∆ ℓ = 2 (∆ℓ DEMNUni−TB ∼ −22 ) with respect to the scale-independent prescriptions in the TB case, and by ∆ ℓ = 1 (∆ℓ DEMNUni−M20 ∼ +14) in the M20 one.We stress here that M ν = 0 is the only case in which the M20 model works better than the TB one.
From the results of this test, we deduce that the use of a scale-dependent or scaleindependent bias is almost indifferent for our modelling purposes, and consequently they are both valid prescriptions.Figure 11: Sensitivity to the galaxy bias scale dependence.Left: Comparison between the simulated ISWRS-galaxy cross-correlation for M ν = 0 eV (solid black line) with theoretical predictions in the linear regime (dashed orange line) and nonlinear regime using the TB model.Three different DEMNUni galaxy bias models have been tested: the scaleindependent DEMNUni-b gg bias (dashed red), the scale-independent DEMNUni-b mg bias (dot-dashed brow) and the scale-dependent DEMNUni b gg (k, z) bias (dotted blue).The middle subpanel of the plot shows a zoom in of the sign inversion region, while the bottom one shows the absolute value of the percentage relative differences of the predictions with respect to the simulations.Right: Same as left, but with the M20 model.

Nonlinear ISWRS-galaxy cross-spectra: the case of different galaxy surveys
The last part of this work consists of computing predictions of the cross-spectra that can be observed on higher redshift ranges than the ones analysed so far, forecasting in particular their dependence on the parameter z m , in Equation (3.2), that characterises a galaxy survey.
Here we compute the ISWRS-galaxy cross-spectra for the four M ν values with three n(z, z m ) corresponding to z m = 0.9, 2, 3. To exploit the full redshift range where the n(z, z m ) are not zero, we compute our integrals up to z max = 8, conversely to the previous Sections where z max = 1.89 (i.e. the DEMNUni upper limit), and we use the b T ut gg (z, M ν ) extended up to those redshifts via the fitting formula in Equation (3.4).To investigate the effect of the massive neutrinos presence, we show in Figure 12 the behaviour of ℓ inv as a function of M ν for both the Halofit models.It is evident at first glance that the higher the redshift, the larger is the difference between the four M ν values, for both the Halofit models considered.This is because the late Universe (i.e.small z) is characterised by the presence of DE [94], that affects the appearance of nonlinearities and is expected to have similar effects as neutrinos on this cross-spectrum.Hence, the further we go back in time, the less DE effects will be present, and the shift we measure will be due only to the neutrino mass.
6 Measuring the neutrino mass with nonlinear ISWRS-LSS cross-spectra?
The results shown in the previous Sections demonstrate the validity of the analytical method against the DEMNUni simulations and that both mocks and predictions confirm what we the smaller is the scale (the larger is the ℓ) where nonlinearities appear, because the RS effect dominates over the late ISW one, as we go back in time Focussing on the nonlinear modelling, we find that TB predicts a higher level of nonlinearities with respect to M20 at small redshifts, but this trend is inverted at z ≥ 1.1 (see Figure 1).
The main result of this work is checking and modelling the presence of the same effects in the corresponding angular cross-power spectra, obtained from P ΦΦ and P Φδ , using Equations (2.5)- (2.6).
Concerning the cross-correlation between the ISWRS signal and the CMBL potential, we find that both the DEMNUni maps and the theoretical predictions, computed using the TB and M20 models, show the sign inversion due to the anti-correlation between the RS and the CMBL potential, and the shift of this sign inversion moves towards larger multipoles as M ν increases (see Figure 8).Furthermore, we observe that, while the M20 model tends to overestimate nonlinearities (and the sign inversions appear always on smaller multipoles with respect to the DEMNUni data), the TB model gives different results depending on M ν , but the corresponding sign inversion locations are closer to the DEMNUni ones for all the M ν values, and specifically overlap with the data for M ν = 0.17, 0.30 eV.We conclude that in the case of the cross-correlation between the ISWRS effect and the CMBL potential, the TB model works better than the M20 one in reproducing nonlinear effects and how they are affected by the presence of massive neutrinos.
We have performed a similar analysis for the cross-correlation between the ISWRS effect and the galaxy distribution (see Figure 9).This represents an improvement with respect to [54], where only the massless neutrino case was investigated, and extends the work of [21] properly treating, in the nonlinear regime, the impact of massive neutrinos on this crosscorrelation.In this case, we find that, as for the cross-correlation between the ISWRS effect and the CMBL potential, the M20 model again overestimates nonlinearities.Instead, for the TB model, the trend here is to underestimate nonlinearities for all the M ν values, but again the ℓ inv are closer to DEMNUni ones with respect to the M20 case.
Considering the dependence on the galaxy bias of the cross-correlation between the ISWRS effect and the galaxy distribution, we have performed two further tests to quantify the effects of different bias modellings.First, we test how the analytical predictions fit the simulated signal when considering a bias (measured from the DEMNUni simulations) dependent or not on the neutrino mass.We observe that, for small values of M ν , the differences between the modelling with massless and massive neutrino bias are negligible (Figure 10).Conversely, when considering large values of M ν , the difference cannot be neglected, and a bias model specific to that cosmology is required, as expected.Second, only in the massless neutrino case, we test the effect on the theoretical predictions of a scale-independent or scale-dependent bias models, both calibrated against measurements from the DEMNUni suite.We verify that, for both the TB and the M20 modellings, choosing a scale-dependent bias b gg (k, z) produces analytical cross-spectra comparable to those obtained with a scale-independent bias b gg (z), as shown in Figure 11.Hence, we can conclude that both the bias models are valid for this kind of analysis.
Moreover, we have verified that the ISWRS-galaxy cross-spectra change with redshift because of the different role played by DE at different epochs.Therefore, a straightforward extension of this work would be to include DDE models in the analysis, and study the dependence of the ISWRS-LSS cross-spectra on the DE equation of state.This could help to distinguish between different DE models and provide insights on the properties and the evolution of DE.
Finally, our results highlight the strict connection between the sign inversion location, ℓ inv , and the value of neutrino masses.However, the cosmic variance limit on these measurements makes the detection of ℓ inv highly unlikely.Nevertheless, the accuracy of future surveys might allow to distinguish different M ν values via the estimation of the overall amplitude of the ISWRS signal and its cross-correlation with LSS on very small cosmological scales.
Overall in this work, we have developed an analytical modelling of the nonlinear ISWRS-LSS cross-correlations that exploits the Boltzmann solver code CAMB and has been validated against the DEMNUni simulations [25].From the comparison of the two nonlinear Halofit models, we have observed that the TB model works better for the reconstruction of the two cross-correlations analysed, with respect to the M20 one.Ours is a useful and computationally non-expensive tool to investigate these cross-correlations in different cosmological scenarios without the need to run large N-body simulations for each cosmology.In a future work we will investigate the detectability of the neutrino mass via ISWRS-LSS cross-correlations.We expect future CMB experiments, such as The Simons Observatory [98], LiteBIRD [99], CMB-S4 [44], and galaxy surveys such as Euclid [36], the Roman Space telescope [100], the Vera Rubin Observatory [101], will help us exploiting such nonlinear signals with the aim of measuring, in combination with more traditional probes, the neutrino mass and the DE EoS.

26 k 98 kFigure 2 :
Figure 2: Full-sky maps of ISWRS (left column) and CMBL potential (right column) obtained from the DEMNUni simulations for the four cosmologies used in this work: ΛCDM and νΛCDM with M ν = 0.17, 0.30, 0.53 eV (from top to bottom).These maps have a pixel number and resolution given by N side = 2048 and a RING ordering, following the procedure described in Section 3.2.

Figure 5 :
Figure 5: Left: Galaxy selection functions: the solid blue line is extracted from DEMNUni data; the dashed orange line is n(z, z m ) from Equation (3.2) with z m = 0.9, which is the one assumed for the photometric Euclid survey; the dot-dashed green and the dotted red lines are n(z, z m ) from Equation (3.2) with z m = 2, 3, respectively[21].Right: Scale-independent galaxy bias functions: the dashed black line is the Tutusaus bias[89]; the solid red, dashed blue, dot-dashed green and dotted purple lines are b gg (z, M ν ) in Equation(3.4)  for M ν = 0, 0.17, 0.30, 0.53 eV, respectively; finally the dotted pink line is b gg (z) = √ 1 + z used in[35].The shaded grey area highlights the extrapolation of the galaxy bias function above the redshifts tested against the Euclid Flagship 1 simulation[89].