Ray tracing the integrated Sachs-Wolfe effect through the light cones of the dark energy universe simulation-full universe runs

: The late integrated Sachs-Wolfe (ISW) effect correlates the cosmic microwave background (CMB) temperature anisotropies with foreground cosmic large-scale structures. As the correlation depends crucially on the growth history in the era of dark energy, it is a key observational probe for constraining the cosmological model. Here we present a detailed study based on full-sky and deep light cones generated from very large volume numerical N-body simulations, which allow us to avoid the use of standard replica techniques, while capturing the entirety of the late ISW effect on the large scales. We postprocess the light cones using an accurate ray-tracing method and construct full-sky maps of the ISW temperature anisotropy for three different dark energy models. We quantify in detail the extent to which the ISW effect can be used to discriminate between different dark energy scenarios when cross-correlated with the matter distribution or the CMB lensing potential. We also investigate the onset of nonlinearities, the so-called Rees-Sciama effect which provides a complementary probe of the dark sector. We find the signal of the lensing-lensing and ISW-lensing correlation of the three dark energy models to be consistent with measurements from the Planck satellite. Future surveys of the large-scale structures may provide cross-correlation measurements that are suﬀiciently precise to distinguish the signal of these models. Our methodology is very general and can be applied to any dark energy or modified gravity scenario as long as the metric seen by photons can still be characterized by a Weyl potential. The late integrated Sachs-Wolfe (ISW) effect correlates the cosmic microwave background (CMB) temperature anisotropies with foreground cosmic large-scale structures. As the correlation depends crucially on the growth history in the era of dark energy, it is a key observational probe for constraining the cosmological model. Here we present a detailed study based on full-sky and deep light cones generated from very large volume numerical N-body simulations, which allow us to avoid the use of standard replica techniques, while capturing the entirety of the late ISW effect on the large scales. We postprocess the light cones using an accurate ray-tracing method and construct full-sky maps of the ISW temperature anisotropy for three different dark energy models. We quantify in detail the extent to which the ISW effect can be used to discriminate between different dark energy scenarios when cross-correlated with the matter distribution or the CMB lensing potential. We also investigate the onset of nonlinearities, the so-called Rees-Sciama effect which provides a complementary probe of the dark sector. We find the signal of the lensing-lensing and ISW-lensing correlation of the three dark energy models to be consistent with measurements from the Planck satellite. Future surveys of the large-scale structures may provide cross-correlation measurements that are sufficiently precise to distinguish the signal of these models. Our methodology is very general and can be applied to any dark energy or modified gravity scenario as long as the metric seen by photons can still be characterized by a Weyl potential.


I. INTRODUCTION
In a flat Friedmann-Lemaître-Robertson-Walker (FLRW) universe, the presence of dark energy (DE) generates a distinct imprint on the cosmic microwave background radiation through the late integrated Sachs-Wolfe effect [1]. This effect originates in the decay of the gravitational potentials associated with large-scale structure whose growth rate is altered by the increasingly fast expansion driven by the DE component. CMB photons traveling through the structures gain a small energy variation that generates temperature anisotropies at large angular scales. This is not the case in a matter dominated universe where the cosmic expansion exactly compensates the growth rate of structures, rendering the gravitational potentials constant in time. Because of this, the detection of the late ISWeffect in a flat universe is a direct probe of DE. However, this signal needs to be disentangled from that of other effects contributing to the CMB temperature anisotropies. For instance, one can use the cross-correlation with tracers of the distribution of largescale structure (LSS) [2].
The advent of CMB satellite experiments mapping the full sky distribution of temperature anisotropies such as the Wilkinson Microwave Anisotropy Probe (WMAP, [3]) and Planck [4] have made possible the realization of crosscorrelation analyses of the CMB temperature anisotropy maps with the large-scale distribution of structures from galaxy surveys. These studies have resulted in numerous detections of the late ISW signal [5,6] and provided cosmological parameter constraints complementary to those inferred from other cosmic probes (see e.g., [7]).
The next generation of galaxy surveys such as LSST, 1 Euclid 2 and SKA 3 will map the distribution of cosmic structures over volumes of the universe much larger than those currently probed. These will allow for more precise measurements of the ISW and to test the properties of DE [8][9][10][11]. These programs will measure the clustering of galaxies over an unprecedented range of scales, covering large linear modes to small ones, where the nonlinear dynamics of the gravitational collapse of matter induces time-variations of the gravitational potentials. Such nonlinearities produce temperature anisotropies at small scales through the so called Rees-Sciama (RS) effect [12]. Due to the complexity of the nonlinear regime of matter clustering, cosmological model predictions have to rely on N-body simulations. However, this is challenging for studies of the ISW-RS effect, which by spanning from large to small nonlinear scales require large-volume, highresolution simulations.
In the past, several works in the literature have investigated the ISW-RS effect using N-body simulations. These studies differ not only in the characteristics of the N-body simulations used, but also in the adopted ray-tracing methods which are necessary to compute the temperature anisotropies by means of a numerical integration along the photon path [13][14][15][16][17][18]. For instance, [13] used highresolution simulations of a relatively small volume which allowed them to study the RS effect only. In contrast [16] used a simulation of similar resolution with a volume 64 times larger, which allowed them to investigate both the ISW signal on the large scales and the RS contribution at small ones. Both these studies have determined the variation of the gravitational potential along the photon path from multiple patched copies of the simulated volumes. This methodology has the advantage that the variables of interests are continuous at the replica's boundaries due to periodic boundary conditions of the simulations. On the other hand, such replica induce spurious effects since a photon propagating along the direction parallel to any of the principal axes of the simulation box will encounter the same structure several times. Hence, in order to reduce such artifacts, the computation must be limited to oblique photon trajectories crossing the box replica at different points with a maximum radial depth set by the simulation box length. This suppresses the contribution of such spurious effects, but it does not solve the problem entirely. A different approach consists instead in patching repeated boxes that have been randomized by applying a random translation and rotation to the box coordinates. This avoids the repetition of structures along the photon path, but on the other hand it introduces artifacts due to the discontinuity at the replica boundaries. To limit such effects, [18] have used replica in which boxes contributing to the same redshift shell have undergone the same random translation and rotation (see [19]). Hence, each redshift shell around the observer has a different randomization. This reduces the effect of discontinuities at the boundaries of replica, while eliminating any preferred direction in the simulated skymaps. Still, this approach does not account for physical correlations among structures residing in nearby redshift shells. Moreover, in all these works the integration along the photon trajectory is carried out interpolating data from snapshots of the simulation box at given redshift outputs. However, the number of redshift outputs of a simulation remains well below the actual time resolution of the simulations.
Here, we describe a comprehensive methodology which addresses all these limitations. To illustrate the strength of this approach we confront the ISW signal from numerical simulations to the prediction of the linear theory. We also investigate the imprint of nonlinearities on the CMB-LSS correlation and its cosmological dependence on the DE model. To this end we use the full-sky light-cone data from the dark energy universe simulation-full universe runs (DEUS-FUR) [20,21] in combination with a sophisticated ray-tracing technique that solves the photon geodesic equations along the photon trajectory. Our approach takes advantage of the fact that the DEUS-FUR simulations cover the volume of the entire observable universe allowing to generate full-sky light-cone data without the need of recurring to any replica method. Moreover, the light-cone data have been generated during the DEUS-FUR runtime which preserves the full time resolution of the simulations in the computation of the photon trajectories. Finally, the raytracing method allows us to integrate the geodesic equations along the perturbed photon path, rather than the usually adopted Born approximation along the unperturbed path. In the present article we focus mainly on large scales, nevertheless our methodology is very general and can be applied to smaller scales where nonlinearities enhance the differences between predictions of different cosmological models.
The paper is organized as follows. In Sec. II we briefly review the equations underlying the ISW effect. In Sec. III we will describe the numerical simulation datasets and the ray-tracing method. In Sec. IV we will present the results, while in Sec. V we will discuss our conclusions.

II. ISW EFFECT
In order to gain an intuitive understanding of the various contributions to the energy of photons traversing structures in an expanding universe, let us consider a linearly perturbed FLRW metric in longitudinal gauge (given here in natural units), where τ is the conformal time, aðτÞ is the scale factor and ψ is the Newtonian potential (we neglect anisotropic stress, vector and tensor perturbations). In this coordinate system the energy of a photon k 0 evolves as ADAMEK, RASERA, CORASANITI, and ALIMI PHYS. REV. D 101, 023512 (2020) where H ¼ d ln a=dτ and n i is the unit vector pointing into the direction in which the photon is traveling. A first-order integral of this equation yields the following perturbative expression for the observed photon redshift: where v is the peculiar velocity, dχ is the conformal distance element (not the line element which is always null for a photon), the subscript "src" refers to the source and the subscript "obs" refers to the observer. The measured energy of a CMB photon is an observable and hence does not depend on the coordinate system used to describe it. However, the advantage of using longitudinal gauge coordinates is that the terms appearing in the above expression have a clear interpretation. We may notice that apart from the factor related to the overall expansion of the background there are two types of contributions that are local at the source and observer, namely the Doppler shift due to peculiar motion and the gravitational redshift due to time dilation. The last term, which is an integral along the photon path, is precisely what we call the ISW-RS effect. Its origin is also quite intuitive: a photon will receive a net blueshift if it travels through a decaying (∂ψ=∂τ > 0) potential well, since it gains more energy on the infall than it loses climbing out of the increasingly shallow well.
It is worth remarking that N-body simulations are usually not realized in the longitudinal gauge. 4 Hence, in numerical simulation analyses the full expression for the observed redshift should be properly gauge-transformed, which introduces some new local terms that can be of the same order as the gravitational redshift (see [23]). However, the integrated contribution we are interested in here is readily gauge invariant, i.e., up to second-order corrections the integral can be taken directly in the coordinates of the Newtonian simulation.

III. METHODOLOGY
We use numerical data from the DEUS-FUR project taking advantage of the available full-sky and deep light cones (with no replica) for several observers and for three cosmological models.

A. Cosmological models
The cosmological models of the DEUS-FUR simulations consist of a flat ΛCDM model, a quintessence model with Ratra-Peebles potential (RPCDM, [24]) and a phantom dark energy model with constant equation of state w<−1 (wCDM, [25]). The cosmological model parameters have been calibrated to fit the CMB anisotropy power spectra from the WMAP-7 yr data [26] and the luminosity distance measurements from supernova Ia standard candles [27]. These models are known as realistic models [28].I n particular, the values of cosmic matter density Ω m and the normalization of the root-mean-square fluctuations σ 8 of the nonstandard dark energy models have been chosen within the 68% confidence region along the degeneracy line of the Ω m − w and σ 8 − w planes respectively, such as to be statistically indistinguishable from the ΛCDM best-fit model at 1σ (see Fig. 1 in [29]). A summary of the cosmological model parameter values and the simulation characteristics is reported in Table I.
The RPCDM and wCDM models are characterized by a background expansion and the linear growth of density fluctuations which bracket those of the ΛCDM. We can see this in Fig. 1, where we plot the redshift evolution of the Hubble function (left panel), the deceleration parameter (center panel) and the linear growth rate (right panel) for the ΛCDM, RPCDM and wCDM models, respectively. First, we can see that the accelerated expansion starts earlier in RPCDM than ΛCDM (i.e., z acc RPCDM >z acc ΛCDM ), while it starts later in wCDM. Also notice that during the preceding phase, the cosmic expansion in RPCDM is less decelerated than in ΛCDM (i.e., q RPCDM <q ΛCDM ), while is more decelerated in wCDM. As we may notice from the evolution of the linear growth rate, this implies that structures will grow more efficiently in wCDM compared to the ΛCDM case (i.e., f wCDM >f ΛCDM ) and less efficiently in RPCDM during the decelerated phase. These differences are also present during the subsequent accelerated phase of expansion, although exhibiting a different slope at low redshift due to the different rate of cosmic acceleration specific to each model. Given the fact that the variation of the gravitational potentials is proportional to the redshift variation of the growth rate of matter density fluctuations, the trends shown in Fig. 1 entirely characterize the cosmological dependence of ISW signal. While this will be discussed in detail in Sec. IV, for the time being it is informative to present an estimate of the ISW signal for the DEUS-FUR cosmologies assuming the linear theory. In particular, a simple expression of the ISW power spectrum can be obtained through the Limber approximation: where P δδ is the linear power spectrum and the weight functions q i ðχÞ are given by with f being the linear growth rate. Equation (4) provides an explicit relation between the linear growth rate and the ISW power spectrum. Using Eq. (4), we plot in Fig. 2 the increment of the linear ISW power spectrum per unit of comoving distance. We can see that this quantity varies in a cosmological model dependent way over very large distances and reaches a maximum around 2 Gpc=h to about 4 Gpc=h. This implies that the light cones built using N-body numerical simulations and which are used to estimate the ISW effect must be deep enough to include the full time variation of the integrand ∂ψ=∂τ. In particular, light cones must allow to integrate the time variation of the gravitational potential up to comoving distances of more than 7 to 8 Gpc=h.A s we will explain in the following subsections, the light cones generated from the DEUS-FUR simulations have optimal characteristics to perform such an integration.
FIG. 1. Hubble function H, deceleration parameter q and linear growth rate f as a function of redshift z for the ΛCDM (blue solid lines), RPCDM (magenta solid lines) and wCDM (green solid lines) models, respectively. As reference for the latter two we plot as horizontal black lines the constant values of the Einstein-de Sitter model with In what follows, except otherwise stated, we will not rely on the Limber approximation for analytical calculations but rather on the full Bessel integrals as computed by the Boltzmann code CLASS [30].

B. Simulations
DEUS-FUR comprises three N-body simulations of a ð21 Gpc=hÞ 3 volume with 8192 3 particles of a flat ΛCDM model and two DE scenarios with different expansion histories (see Table I). The simulations have been run using the application AMADEUS-a multipurpose application for dark energy universe simulation expressively developed for the realization of the DEUS-FUR project [20,21]. This includes the code generating Gaussian initial conditions for which we use an optimized version of MPGRAFIC [31], the N-body solver for which we use a specifically modified version of the RAMSES code [32] such as to run on a very large number of cores (≈80000) and a parallel friends-offriends halo finder pFoF as described in [33]. RAMSES solves the Vlasov-Poisson equations using an adaptive mesh refinement (AMR) particle-mesh method with the Poisson equation solved with a multigrid technique [34]. We refer the reader to [20] and [21] for a detailed description of the algorithms and optimization schemes adopted in the realization of the DEUS-FUR project.
At coarse level the grid of the DEUS-FUR simulations contains 8192 3 cells, these are allowed to be refined six times reaching a formal spatial resolution of Δx ¼ 40 kpc=h, while the particle mass resolutions is m p ≈ 10 12 M ⊙ =h for the different models. Although these are the first simulations covering the volume of the full observable Universe at such a resolution, this is still relatively low compared to that of other ISW studies in the literature, which have used smaller volumes [16,18]. Hence, our study will focus more on the large scales, where the DEUS-FUR simulations provide cosmic-variance limited errors [35]. We now discuss the characteristics of the light-cone data and the ray-tracing technique.
C. DEUS-FUR light cones DEUS-FUR light cones were built on-the-fly using an onion-shell approach [36,37]. At each coarse time step the shell at the appropriate conformal distance from the observer is recorded (thus moving successively closer and closer to the observer as the simulation advances in time). The time steps multiplied by the speed of light are larger than the resolution of the spatial mesh, consequently each shell has a certain thickness of many coarse cells. The total number of shells is of order ∼400-450 depending on the cosmology. For each of the three cosmologies, five full-sky light cones (fraction of the sky f sky ¼ 1) were stored up to maximum redshift z max ¼ 30. Each of the light cones corresponds to a specific space-time location of the observer given by the position vector (x 1 obs =L box , x 2 obs =L box , x 3 obs =L box , z obs ) where x i obs are the observer Cartesian coordinates in the simulation box going from 0 to L box , L box is the box length and z obs is the redshift of the observer. The locations were chosen as (0.5, 0.5, 0.5, 0.), (0.5, 0.1, 0.1, 0.), (0.1, 0.5, 0.9, 0.), (0.5, 0.5, 0.5, 0.5), and (0.5, 0.5, 0.5, 1.2) so as to maximize the distance between the observers at z ¼ 0.Inthisarticlewewill focus on the two first light cones which do not overlap up to a redshift z ∼ 5.6. For each light cone, two kinds of data are stored: particles and AMR cells. Particle light cones contain particle properties (position, velocity and redshift) while gravity light cones contain mesh properties (position, gravitational field, potential, density, AMR cell's child index). We use the particle light cones to compute maps of the comoving dark matter density on the unperturbed light cone. In contrast, given that our line element Eq. (1) is fully specified by the potential, we perform the ray tracing using the gravity light cones, which allow us to compute ISW-RS maps as well as weak-lensing maps.

D. Ray tracing and map making
In order to compute the ISW temperature anisotropies from the DEUS-FUR light-cone data we use the ray-tracing algorithm part of the fast parallel C++ AMR library Magrathea developed in [38]. For each ray a past-null direction is chosen at the observer, and the geodesic equations, consisting of a coupled set of ordinary differential equations, are numerically integrated backwards in time. Following the true photon path, the gravitational potential and its gradient are obtained from the adaptive mesh by multilinear (cloud-in-cell) interpolation. The set of differential equations is then solved by a Runge-Kutta fourth-order method with adaptive time steps (four time steps per AMR cell). Our methodology is very general as long as the time derivative of the potential can be measured from the gravity light cone. In particular, we do not need to specify any relation between the matter density field and the gravitational potential, which means that a wide range of dark energy and modified gravity models can be studied without changing the analysis pipeline. We consider this approach to be an improvement over previous ones also because the light cones cover the full sky without replica, and the geodesic integration is done at the resolution of the AMR simulation.
DEUS-FUR was not designed to specifically investigate the ISW-RS effect, and because of this the time derivative of the potential was not stored. Nevertheless, it can be recovered at large scales (i.e., larger than the shell size) from the gravitational field (the spatial gradient of the potential computed at a given time) and a total derivative (computed along the light cone). More specifically, we can rewrite the partial time derivative ∂ψ=∂τ as where we have computed the spatial gradient of the potential on the mesh assuming a five-point stencil finite difference approximation. We have generated full-sky maps of the ISW temperature anisotropies of resolution N side by numerically integrating the geodesic equation for all the 12 × N 2 side light rays received from the directions of the HEALPix pixels [39]. Here, we specifically set N side ¼ 512. In addition to the ISW signal we also compute the lensing deflection angle, which allows us to reconstruct the CMB lensing potential that we will discuss in Sec. IV C. To this purpose we first generate a map of the deflection vector, then we extract the lensing potential as the generator of its curlfree part. A small curl part is also present due to higher-order lensing, but we do not study this here.
In Fig. 3 we show full-sky maps of the ISW temperature anisotropies (left panels) and the lensing potential (right panels) generated from the light cones of the three DEUS-FUR cosmological models. In the next section we use the HEALPix package [39] to evaluate angular correlation functions of the full-sky maps.

IV. RESULTS
A. ISW-temperature power spectrum In Fig. 4 we plot the angular power spectrum obtained from the estimator FIG. 3. Full-sky maps of the ISW temperature anisotropy (left panels) and lensing potential (right panels) extracted from the DEUS FUR for (from top to bottom) RPCDM, ΛCDM, and wCDM models for one observer included in the simulation box. As usual, the temperature dipole is not shown since its observation is difficult due to the large kinematic dipole. For a better comparison we therefore also subtract the lensing dipole even though it is taken into account in our analysis.
where a TT lm are the coefficients of the spherical harmonic decomposition of the full-sky maps of the temperature anisotropies as computed using the HEALPix package. For each cosmology, we compute the late ISW-RS signal for four different starting redshifts z Ã , integrating the photon trajectories between z Ã and z ¼ 0. The error bars show statistical errors after combining the estimators from the light cones of two independent observers and binning the data logarithmically in l. From the dependence of the signal on z Ã , we can see that the differences among the various model predictions correlate with the differences in the redshift evolution of the linear growth rate shown in Fig. 1. More specifically, we can see that in the wCDM case (upper right panel) the bulk of the signal is generated at low redshift. As an example, between z Ã ¼ 1 and 3 the ISW signal at l ∼ 10 only increases by 16%, while in the same redshift interval the signal increases by 20% in the ΛCDM case and 40% for RPCDM. This is consistent with the fact that at these redshifts the linear growth rate of the wCDM model remains closer to the Einstein-de Sitter value (f ∼ 1), than ΛCDM and RPCDM, respectively. Moreover, as shown in Fig. 1, the linear growth rate of wCDM (RPCDM) is systematically higher (lower) than the ΛCDM case, thus causing a smaller (larger) ISW signal relative to the ΛCDM prediction as shown in the lower right panel of Fig. 4.
Overall, we find a good agreement with the linear theory for multipoles l ≲ 50 where the linear ISW effect is expected to dominate over the nonlinear RS effect. FIG. 4. Angular power spectra of the ISW signal for three different cosmologies, using logarithmic l-band powers. In order to show where the signal is mostly generated, the temperature anisotropy is calculated by integrating along the line of sight from the observer (at z ¼ 0) to various finite distances z Ã . For the two evolving DE models (wCDM and RPCDM) the lower right panel shows the relative difference with respect to ΛCDM. In all panels, dashed lines indicate the predictions from linear theory as computed by CLASS.
Instead over the range 60 ≲ l ≲ 100 we notice a change of the slope of the power spectrum and a departure from the linear theory that is cosmological model dependent. This departure marks the onset of the nonlinear RS effect. It is worth remarking that while the amplitude of the ISW power spectrum is mostly determined by the linear growth rate through the (1 − f)-factor in Eq. (4), the scale of deviation from the linear theory primarily depends on the amplitude of the matter density power spectrum (i.e., σ 8 D þ , where D þ is the linear growth factor). We can see that at multipoles l ∼ 100 the trends flatten to a plateau. However, at these multipoles our estimation of the temperature fluctuations is dominated by noise that results from the limitations of the numerical computation. In fact, our integration method evaluates the signal from discretized data. In particular, as explained earlier, ∂ψ=∂τ is not directly computed during the simulation run, but estimated from the total derivative and the spatial gradient of the potential over the mesh. The signal is therefore effectively recovered by summing the jumps of the gradient between the time steps of the simulation, or in other words by sampling from a relatively small number of locations. This procedure prevents an exact evaluation of the temperature fluctuations on scales smaller than the sampling rate. Consequently, at the multipoles corresponding to these scales the correlation of the numerical noise with itself (whose spectral properties are consistent with a Poisson process) provides the dominant contribution to the power spectrum, leaving the signal of the RS effect embedded in noise.

B. ISW-matter density correlation
As already mentioned, the ISW signal in the CMB temperature anisotropy autopower spectrum cannot be disentangled from the contribution of other processes which generate temperature fluctuations at the last scattering surface. Instead, it can be detected by cross-correlating CMB temperature maps with the late-time distribution of cosmic structures. Temperature anisotropies generated at early time are largely uncorrelated with the spatial distribution of matter density perturbations at late time. In contrast, ISW anisotropies are sourced by the time varying gravitational potentials associated to late-time matter density perturbations which are traced by cosmic structures.
Here, we compute the cross-correlation between the ISW signal and the matter distribution by constructing full-sky maps of the average density contrast in a given distance interval of the observer's light cone, namely: where WðχÞ is a tophat selection function that specifies the distance interval. Despite the fact that we use the information recorded along the past light cone, it should be noted that the density contrast we are evaluating is not a directly observable quantity. First, we use the distribution of dark matter particles to define δ, while observations of largescale structure have to rely on biased tracers such as galaxies. Second, for simplicity we do not take into account relativistic projection effects such as e.g., weak lensing that would perturb the observed number density of tracers in a given solid angle element. While it may seem that this could be rectified by making use of our ray tracing algorithm, we remind the reader that the particle positions are not provided in the appropriate gauge 5 for this task. This gauge issue is not very important on small scales, but it could contaminate the large angular scales we are mostly interested in [40]. The quantity δðnÞ should therefore be understood merely as an ideal theoretical probe that is easy to compute and that can be used as a proxy for assessing the ISW-matter crosscorrelation in different dark energy models. Fig. 5 shows the ISW-matter cross-correlation power spectrum of the three DEUS-FUR models for three redshift bins covering the range 0 ≤ z<0.5, 0.5 ≤ z<1, and 1 ≤ z<8 respectively. We use a power spectrum estimator similar to that described in the previous section and again applied to the full-sky maps of two independent observers in each case. We also show the theoretical model predictions obtained from linear theory (dashed lines). As expected, for all three cosmologies the cross-correlation signal is larger at lower redshift, and the relative contribution of the different redshift bins follows the trends already discussed in the previous section. Higher multipoles are systematically more correlated at higher redshift which is simply a geometric effect. The lower-right panel of Fig. 5 shows the relative difference of the cross-spectra with respect to the ΛCDM case.
Overall, we can see that at relative low multipoles (l < 100) or equivalently at large angular scales, the linear theory reproduces reasonably well the numerical FIG. 6. The left panel shows the angular cross-power spectra between CMB temperature and matter density contrast (integrated between z ¼ 0 and z ¼ 9) at small scales or large multipoles (l > 80) where nonlinear effects become important. The right panel shows the relative difference of the cross-power spectra to the ΛCDM linear prediction. The different lines correspond to the DEUS-FUR numerical simulation (solid lines), linear prediction (dashed lines) and HALOFIT prediction (dotted lines). 5 As our ray tracer uses longitudinal gauge it requires the number count per coordinate volume in that gauge. This quantity is different in different gauges and only relates to a gaugeinvariant number count per redshift space volume after the relativistic projection effects have been included. cross-spectra from the DEUS-FUR simulations. The agreement is stronger for the higher redshift bins, while we may notice a departure from the linear theory in the lowest redshift bin at l ≈ 100. The exact angular scale of such departure depends on the underlying cosmological model and marks the onset of nonlinearities contributing to the RS effect.
To highlight this point in Fig. 6 we plot the cross spectra (left panel) and their relative differences with respect to the linear ΛCDM prediction (right panel) at high multipoles (l > 100), where we have integrated the signal between z ¼ 0 and z ¼ 9 to have a larger signal-to-noise. Notice that as we explore smaller scales, we approach the spatial resolution of the simulations. This manifests in a drop of the cross-power spectra at large multipoles and in the presence of noise. However, these spurious numerical effects are mostly independent of the underlying cosmology and cancels out when plotting relative differences. This allows us to recover the signature of nonlinearities at small scales which is expected to depend on the cosmological model [28]. We can clearly see this in the right panel of Fig. 6, where the differences of the wCDM and RPCDM cross-spectra with respect to the linear ΛCDM prediction increase at larger multipoles and correlate with their σ 8 D þ value. We also plot an analytical prediction of the nonlinear regime following [41]. At these small scales, the Limber approximation is applicable and we use Eq. (4) with inside the top-hat region of width Δχ and q 2 ðχÞ¼0 outside. P δδ in Eq. (4) becomes the nonlinear power spectrum and f becomes the (scale-dependent) nonlinear growth rate fðk; χÞ¼ 1 2 d ln P=d ln a. As shown in [28], the imprint of nonlinearities on the matter power spectrum depends on the specificities of the underlying dark energy model. However, these are not captured by the standard HALOFIT prescription [42]. Instead, we use the prescription from [43], which accounts for dark energy models with a constant equation of state w. This covers the ΛCDM and wCDM cases, while for the RPCDM model, we approximate its time varying equation of state with its mean value across the cosmic history, corresponding to w ¼ −0.83. Notice though that the imprints of dark energy on nonlinear scales play an important role (at several percent level) only beyond modes k ≳ 1 h Mpc −1 . As we focus on large scales corresponding to k<1 h Mpc −1 , in this regime the prescription by [43] provides accurate enough predictions of the matter power spectra. Moreover, given the level of fluctuations exhibited by the numerical spectra at these scales, nonlinear model predictions accurate at the percent level are not needed in our analysis. In the left panel of Fig. 6 we plot the nonlinear predictions of the ISW-density cross-spectra for the ΛCDM, wCDM and RPCDM models against the numerical simulation results. As we can see the cosmological dependence of analytical predictions and simulations are similar. In the right panel of Fig. 6 we plot the relative differences with respect to the ΛCDM case.
We can see that the onset of the RS effect depends on the model of dark energy, thus confirming previous claims in the literature that such a scale can be used as a probe of dark energy [41,44]. More accurate prescriptions for the computation of the nonlinear matter power spectrum of dynamical dark energy models are required if we had extended the analysis deep in the nonlinear regime. In such a case prescriptions such as those developed in [45,46] are needed. However, this is beyond the scope of this paper.
C. ISW-lensing potential correlation The gravitational potentials of large-scale structure modify the temperature anisotropy pattern through weak gravitational lensing. The lensing potential for the CMB can be constrained from CMB observations directly [47], without the need for external large-scale structure information. Since lensing and ISW effect both originate from the same matter perturbations they are expected to be highly correlated at redshifts where the potentials are decaying.
We construct the lensing potential by taking the curlfree part of the total deflection angle that we compute by integrating the photon geodesic equations for the full sky. In Fig. 7 we plot the angular power spectra of the lensing potentials (top panels) and the cross-correlations with the ISW signal (bottom panels) for our three cosmologies. The right panels show the relative difference with respect to ΛCDM.
Again, for low multipoles (l < 100) we find an excellent agreement with the predictions from the linear theory (dashed lines) both for the lensing potential autocorrelation function and for the lensing-ISW cross-correlation. This validates again our methodology. An analytical calculation of the lensing-ISW cross-power spectrum following [41] with Limber approximation (i.e., using the lensing weight for q 2 ) indicates that the scale for nonlinearity is beyond l ¼ 500 and therefore at scales below the resolution of our data. This is because the CMB-lensing kernel peaks at high redshift pushing the scale of nonlinearities toward larger l.
The cosmological dependence reaches 30% (between RPCDM and wCDM at l ¼ 100) for the autocorrelation of the lensing potential. It appears more pronounced for the cross-correlation between the ISW signal with the lensing FIG. 8. Comparison of DEUS-FUR cross-spectra XðlÞ for κ-κ (left) and T-κ (right) to the Planck measurement (solid black line) from [6] (we follow the convention of that paper). The three models agree with the data. In the future, reducing systematics and the use of multiple probes will allow to disentangle the three "realistic" cosmologies studied in this paper. potential where it reaches 50% and the order of the spectra is reversed. Moreover while wCDM and ΛCDM are indistinguishable from a lensing perspective, they are more than 1-σ away for the ISW-lensing cross-correlation (where σ 2 is the cosmic variance). This illustrates the complementarity of the lensing and ISW probes.
Finally, in Fig. 8 we confront our CMB lensing and ISW spectra to measurements from the Planck collaboration [4]. Given the observational systematics, the measurement errors are not down to cosmic variance. The three models are therefore compatible with the current measurements. It will however be possible in the future to lower the systematics, combine several probes (ISW-density, ISW-lensing) at several redshift, and explore the RS effects. Altogether this should allow the discrimination between various cosmological models compatible with CMB and supernovae data.

V. CONCLUSIONS
We have presented a study of the ISW-RS effect using data from the dark energy universe simulation-full universe runs of three different dark energy models. As these simulations encompass the full observable volume of the simulated cosmologies, we have been able to perform a thorough analysis of the ISW signal using full-sky lightcone data with a photon propagation across the cosmic nonlinear matter distribution up to very high redshifts and for two observers located at two different space-time positions. The usual spurious effects inherent to replica methods are thus completely avoided. Moreover, we used a sophisticated ray-tracing technique that solves the photon geodesic equation along the photon trajectory without recourse to the Born approximation. Our methodology is very general and applicable to both large and small scales.
The cosmological models of the DEUS-FUR simulations are a flat ΛCDM model, a quintessence model with Ratra-Peebles potential (RPCDM) and a phantom dark energy model with constant equation of state w ¼ −1.2 (wCDM). We extracted the angular correlation functions from fullsky maps for such realistic cosmological models. The autocorrelation functions of the ISW signal are very weakly dependent on cosmology for small multipoles (l ≪ 50) and in very good agreement with the linear prediction. For multipoles between l ¼ 50 and l ¼ 100, we can see the onset of the RS effect. At larger multipoles the numerical noise dominates the signal. This is caused by the way our integration method extracts the signal from the discretized data for the time derivative of the gravitational potential which was not specifically stored in our simulations. Therefore we can only observe the beginning of the transition to the RS effect as the nonlinear signal is then drowned in the numerical noise.
The ISW-RS signal in the CMB anisotropy power spectrum is not directly observable due to the contribution of early-time processes which generate temperature anisotropies over the same range of angular scales. Nevertheless, it can be detected through cross-correlation of CMB temperature anisotropy maps with the distribution of large-scale structures. This is because the late-time varying gravitational potentials which sources the ISW-RS effect are traced by cosmic structures. Here, we have performed an analysis of the ISW-matter correlation from the DEUS-FUR light-cone data. We find that the differences among the cross-spectra of the DEUS-FUR models to be enhanced at l ≳ 100. In particular, we clearly find a deviation from the linear prediction as nonlinear effects become dominant at small scales. Such effects are expected to be cosmological model dependent and do account for the observed differences between the ISW-correlation spectra of the DEUS-FUR models. In fact, the linear growth rate of matter density perturbations is larger in wCDM than in the ΛCDM case and even more compared to RPCDM. Consequently the cosmic structure formation is more efficient in wCDM than in ΛCDM and RPCDM, respectively. This results in an attenuation of the ISW-matter correlation signal in the wCDM model with respect to the ΛCDM and RPCDM cases, mainly due to a lower ISW amplitude, and a decrease on the scale of deviation from the linear regime of matter clustering.
We have also investigated the correlation of the ISW effect with the lensing potential from the matter distribution. Again, at low multipoles we find the results to be well reproduced by the linear theory. We find the cross-correlation to be particularly sensitive to the underlying cosmological model with the amplitude of the signal varying more than the cosmic variance error among the DEUS-FUR models. The comparison of our results to ISW and lensing data from Planck shows that these measurements cannot yet disentangle between the three dark energy models studied here. However, in the light of the observational data that will be accessible with future survey programs, the ISW-lensing correlation in combination with the detections of the nonlinear features of the ISW-LSS correlation are likely to provide a powerful cosmological proxy and a probe of the nature of dark energy.