Estimating small angular scale CMB anisotropy with high resolution N-body simulations: weak lensing M.J. Fullana1, J.V. Arnau2, R.J. Thacker3, H.M.P. Couchman4, and D. S´aez5 0 1 0 [email protected] 2 n a J Received ; accepted 7 2 ] O C . h p - o r t s a [ 1 v 1 9 9 4 . 1 0 0 1Institut de Matem`atica Multidisciplina`ria, Universitat Polit`ecnica de Val`encia, 46022 1 : v Val`encia, Spain i X 2Departamento de Matem´atica Aplicada, Universidad de Valencia, 46100 Burjassot, Va- r a lencia, Spain 3Department of Astronomy and Physics, Saint Mary’s University, Halifax, Nova Scotia, B3H 3C3 Canada 4Department of Physics and Astronomy, McMaster University, 1280 Main St. West, Hamilton, Ontario, L8S 4M1, Canada 5Departamento de Astronom´ıa y Astrof´ısica, Universidad de Valencia, 46100 Burjassot, Valencia, Spain – 2 – ABSTRACT We estimate the impact of weak lensing by strongly nonlinear cosmological structures on the cosmic microwave background. Accurate calculation of large ℓ multipoles requires N-body simulations and ray-tracing schemes with both high spatial and temporal resolution. To this end we have developed a new code that combines a gravitational Adaptive Particle-Particle, Particle-Mesh (AP3M) solver with a weak lensing evaluation routine. The lensing deviations are eval- uated while structure evolves during the simulation so that all evolution steps– rather than just a few outputs–are used in the lensing computations. The new code also includes a ray-tracing procedure that avoids periodicity effects in a universe that is modeled as a 3-D torus in the standard way. Results from our new simulations are compared with previous ones based on Particle-Mesh simu- lations. We also systematically investigate the impact of box volume, resolution, and ray-tracing directions on the variance of the computed power spectra. We find that a box size of 512h−1 Mpc is sufficient to provide a robust estimate of the weak lensing angular power spectrum in the ℓ-interval (2,000–7,000). For a reaslistic cosmological model the power [ℓ(ℓ + 1)C /2π]1/2 takes on values of a ℓ few µK in this interval, which suggests that a future detection is feasible and may explain the excess power at high ℓ in the BIMA and CBI observations. Subject headings: methods:numerical — cosmic background radiation — cosmol- ogy:theory — large-scale structure of the universe – 3 – 1. Introduction Photons from the last scattering surface of the cosmic microwave background (CMB) detected at redshift zero, are inevitably lensed by cosmological structures. This lensing produces a number of modifications to the CMB angular power spectrum, deviations from Gaussianity, and B-polarization. While these effects have been extensively studied, and a review of their impact can be found in Lewis & Challinor (2006), we are here concerned with the lensing due to strongly nonlinear objects such as galaxy clusters and smaller-scale structures. In particular, we focus our attention on the estimation of the small angular scale (high ℓ) correlations induced by this lensing. This marks a somewhat different direction to the extensive recent efforts to extend N-body simulation techniques to all-sky modelling (e.g. Fosalba et al. 2008; Das & Bode 2008; Carbone et al. 2008; Teyssier et al. 2009; Carbone et al. 2009). N-body simulations are necessary to estimate some CMB anisotropies. This is the case for the gravitational anisotropies (weak lensing and Rees-Sciama effects) produced by strongly nonlinear cosmological structures. While simulations with sufficient resolution to model the formation of galaxies within an entire Hubble volume are currently out of reach, researchers are getting close to this goal (Kim et al. 2008; Teyssier et al. 2009). However, until such simulations are possible in order to estimate gravitational anisotropies in the CMB it is necessary to create an artificial periodic universe where copies of the simulation box at various redshifts are replicated back to a given redshift. While numerous alternatives exist for how this is achieved all the methods propagate CMB photons along appropriate directions within these model universes. When conducting this ray-tracing, box sizes, spatial scales and ray directions must be chosen in such a way that the periodicity of the artificial universe has no appreciable effects on the resulting anisotropy. The ray-tracing procedure used here achieves these goals and was described and applied in earlier papers, – 4 – where it was also compared with previous methods (Jain et al. 2000; White & Hu 2000). Two of these earlier papers dealt with weak lensing (Cerd´a-Dur´an et al. 2004; Ant´on et al. 2005), and two others with the Rees-Sciama effect (S´aez et al. 2006; Puchades et al. 2006). The simulations discussed in these papers were conducted with the Particle-mesh (PM, e.g. Hockney & Eastwood 1988) N-body simulation algorithm. For the large simulation boxes required in our study (see below) the spatial resolution of PM simulations is moderate, at best. For this reason it was suggested in these papers that better N-body simulations might be necessary to estimate good CMB angular power spectra for large multipole indexes (ℓ values). This suggestion is confirmed in the present study, where AP3M (Adaptive Particle-Particle, Particle-Mesh; Couchman 1991) simulations with high resolution are performed to build up an appropriate periodic universe, in which the CMB photons are then moved according to the aforementioned ray-tracing procedure. A CMB angular power spectrum has been estimated from ℓ ∼ 200 to ℓ ∼ 10000, and the accuracy of this spectrum in various ℓ-intervals is discussed. Motivated by a desire to match a number of cosmological constraints, including the five-year WMAP (Wilkinson microwave anisotropy probe) observations, high-redshift supernovae measurements and galaxy surveys, (e.g. Spergel et al.2007; Vianna & Liddle 1996; Riess et al.1998; Perlmutter et al.1999) we consider a cosmological model in which the universe is flat and the initial fluctuation spectrum has an inflationary origin. It is assumed that only scalar perturbations are present and that weak lensing from gravitational waves is negligible. The resulting distribution of perturbations is then Gaussian and the perturbations themselves are purely adiabatic. Post inflation, the power spectrum of the perturbations is almost exactly of the Harrison-Zel’dovich (Harrison 1970; Zel’dovich 1970) form. This power spectrum is then modified during the radiation-dominated phase, as acoustic oscillations in the photon-baryon fluid damp growth, an effect accounted for by our use of a CMBFAST calculation of the transfer function. The following cosmological – 5 – parameters are used for consistency with the standard model (Komatsu et al. 2008; Hinshaw et al. 2008): (1) a reduced Hubble constant h = 10−2H = 0.7 (where H is the 0 0 Hubble constant in units of km s−1 Mpc−1); (2) density parameters Ω = 0.046, Ω = 0.233, b d and Ω = 0.721 for the baryonic, dark matter, and vacuum energy, respectively (the matter Λ density parameter is then Ω = Ω +Ω = 0.279); (3) an optical depth τ = 0.084 which m b d characterizes the reionization; and (4) a parameter σ = 0.817 normalizing the power 8 spectrum of the energy perturbations. The layout of this paper is as follows: in Section 2 we outline salient issues in our map-making technique. We follow this in Section 3 with an outline of N-body methods and our ray-tracing technique. Results from our new simulations are presented in Section 4 along with a critical comparison to earlier work. We then compare our results to current observations and conclude with a brief discussion. 2. Map Construction We begin by describing the procedure for constructing lensed maps of the CMB from small unlensed maps at the last scattering surface. Note that in this section and those following, we choose units such that c = 8πG = 1, where c is the speed of light and G the gravitational constant. For a quantity A, A and A denote the value of A at, respectively, e 0 the time of emission (last scattering surface) and the present. The scale factor is a(t), where t is the cosmological time, and its present value, a , is assumed to be unity, which is always 0 possible in flat universes. The unit vector ~n defines the observation direction (line of sight). Small, unlensed maps of CMB temperature contrasts (∆ = δT/T) must be constructed to be subsequently deformed by lensing. These maps have been obtained with a method based on the fast Fourier transform (S´aez et al. 1996) and lead to small squared-Gaussian – 6 – maps of the contrast ∆. These maps are uniformly pixelised. The code designed to build up the unlensed maps (map making procedure), requires the CMB angular power spectrum, which has been obtained by running the CMBFAST code (Seljak & Zaldarriaga 1996) for the model described above. In order to deform the unlensed maps, the lens deviations corresponding to a set of directions, covering an appropriate region of the sky, must be calculated. These deviations are the quantities (Seljak 1996): λ0 ~δ = −2 W(λ)∇~ φ dλ , (1) ⊥ Z λe where ∇~ φ = −~n∧~n∧∇~φ is the transverse gradient of the peculiar gravitational potential ⊥ φ, and W(λ) = (λ −λ)/λ . The variable λ is e e 1 db −1 λ(a) = H . (2) 0 Z (Ω b+Ω b4)1/2 a m0 Λ In the case of weak lensing, the integral in Eq. (1) is usually integrated along straight paths, ignoring the small deflections associated with the lensing effect. This approximation, commonly known as the Born approximation, is equivalent to keeping first order terms in a positional expansion of the transverse potential as a function of the normal ray and its lensing offset. While this approximation is known to be comparatively inaccurate at small ℓ (e.g. Van Waerbeke et al. 2001), at 1000 . ℓ . 10000 detailed calculations using higher order perturbation theory suggest that corrections to the first order assumption are approximately two orders of magnitude lower than the first order term (Shapiro & Cooray 2006). Recent N-body simulation work examining the validity of the Born approximation (Hilbert et al. 2009) has shown corrections only begin to become significant at the 5% level for ℓ > 20000. Since we are interested in calculating 1000 < ℓ < 10000 values we cautiously accept the errors inherent in the first order approach. Once the deviations have been calculated, they can be easily used to get the lensed – 7 – maps from the unlensed ones. This is achieved using the relation ~ ∆ (~n) = ∆ (~n+δ) , (3) L U where ∆ and ∆ are the temperature contrasts of the lensed and unlensed maps, L U respectively. Given the unlensed map ∆ , and the map ∆ obtained from it after deformation by U L lensing (the lensed map), the chosen power spectrum estimator can be used to get the quantities C (U) and C (L), whose differences C (LU) = C (L)−C (U) can be considered ℓ ℓ ℓ ℓ ℓ as an appropriate measure of the weak lensing effect on the CMB. Moreover, a map of deformations ∆ = ∆ −∆ can be obtained and these can be analyzed to get another D L U angular power spectrum C (D). Since the maps L and U are not statistically independent, ℓ the spectra C (D) and C (LU) appear to be very different (see Ant´on et al. (2005) for ℓ ℓ ~ details). Eight hundred unlensed maps are lensed (deformed) by using the same δ field and the average C (D) and C (LU) spectra are calculated and analyzed. Both spectra are ℓ ℓ very distinct measures of the weak lensing effect under consideration. Many times, only the customary oscillating C (LU) spectrum is shown; however, in a few appropriate cases, ℓ the C (D) spectra are displayed. Our power spectrum estimator was described in detail in ℓ Arnau et al. (2002) and Burigana & S´aez (2003). Results obtained with this estimator were compared with those of the code ANAFAST of the HEALPix (Go´rski et al. 1999) package in Arnau et al. (2002) and also in Puchades et al. (2006). These comparisons showed that our estimator is a very good one in the case of regularly pixelised squared maps as analyzed here. – 8 – 3. N-body simulations and ray-tracing procedure 3.1. N-body simulation technique The primary simulations presented in this work were run using a parallel OpenMP- based implementation of the “HYDRA” code (Thacker & Couchman 2006). This code uses the AP3M algorithm to calculate gravitational forces within a simulation containing N particles. In the AP3M algorithm a cubic “base” mesh of size N cells per side is p c supplemented by a series of refined-mesh P3M calculations to provide sub-mesh resolution. Gravitational softening is implemented using the S2 softening kernel (Hockney & Eastwood 1988) which is remarkably similar in shape to the cubic spline softening kernel used in many treecodes (e.g. Hernquist & Katz 1989). The S2 softening used in the kernel is 2.34×S p where S is an equivalent Plummer softening length which we quote throughout the paper p to enable a simple comparison to other work. The softening length is held constant in physical coordinates subject to the resolution not falling below 0.6 of the interparticle spacing at high redshift. This technique is widely applied (e.g. Springel et al. 2005) and is a compromise between assuring that the potential energy of clusters does not evolve significantly at low redshift, while still ensuring structures and linear perturbations at high redshift are followed with reasonable accuracy. Initial conditions were calculated using the standard Zel’dovich approximation technique (Efstathiou et al. 1985), and all simulations were started at a redshift of z = 50, which is sufficiently early to place modes in the linear regime. To account for the impact of varying box sizes, L , we considered simulations of size 256h−1 Mpc, 512h−1 Mpc, and box 1024h−1 Mpc and a full list of N-body simulation and ray-tracing parameters is given in Table 1. At z = 6 (the beginning of our ray-tracing epoch in the simulation), the box sizes of 256h−1 Mpc, 512h−1 Mpc, and 1024h−1 Mpc correspond to square sky patches with angular sizes, Φ , of 2.48◦, 4.96◦, and 9.92◦, respectively. These are then the sizes of our map – 9 – Table 1. Lensing Simulations Name Algorithm Lbox Modeset Np Mp Nc Sp Ndir zin ∆ps ∆ang Direction /h−1Mpc 1010M⊙ /h−1kpc /h−1kpc ’ RLSAA AP3M 512 A 5123 11 1024 12 512 6 25 0.59 D512A RLSAB AP3M 512 A 5123 11 1024 12 512 6 25 0.59 D512B RLSBA AP3M 512 B 5123 11 1024 12 512 6 25 0.59 D512A RLSCA AP3M 512 C 5123 11 1024 12 512 6 25 0.59 D512A LSLDA AP3M 1024 D 5123 88 1024 24 1024 6 50 0.59 D1024A LSMAB AP3M 512 A 2563 88 512 24 512 6 50 0.59 D512B LS∆S1AA AP3M 512 A 5123 11 1024 24 512 6 50 0.59 D512A LS∆S2AA AP3M 512 A 5123 11 1024 36 512 6 75 0.59 D512A LS∆1AA AP3M 512 A 5123 11 1024 12 512 6 12 0.59 D512A LS∆1AB AP3M 512 A 5123 11 1024 12 512 6 12 0.59 D512B LS∆1BA AP3M 512 B 5123 11 1024 12 512 6 12 0.59 D512A LS∆1BB AP3M 512 B 5123 11 1024 12 512 6 12 0.59 D512B LS∆2AA AP3M 512 A 5123 11 1024 12 512 6 40 0.59 D512A LS∆2AB AP3M 512 A 5123 11 1024 12 512 6 40 0.59 D512B LS∆2BA AP3M 512 B 5123 11 1024 12 512 6 40 0.59 D512A LS∆2BB AP3M 512 B 5123 11 1024 12 512 6 40 0.59 D512B LSHEA AP3M 256 E 5123 1.4 1024 6 512 6 15 0.29 D256A LSMGA AP3M 512 G 2563 88 512 24 512 6 60 0.59 D512A LSLHA AP3M 512 H 1283 704 256 48 512 6 120 0.59 D512A PMLS PM 256 F 5123 1.4 512 1000 256 6 500 0.59 D256A Note. —Listoftheparametersusedinthelensingsimulations. ColumnsthreetoeightgiveinformationassociatedwiththeN-bodysimulations,namely theboxsize,Lbox,thesetofinitialmodes(determinedbytheboxsizeandasinglerandomseed),thenumberofparticlesused,Np,themassoftheindividual particles,Mp,thenumberofFouriercellsalongeachsideofthesimulation,Nc,andthePlummersofteninglengthused,Sp. AllN-bodysimulationswere startedataredshiftofz =50. Thefollowingfivecolumnsgiveparametersassociated withthelensingcalculation,namelythenumberofraysalongthe edgeofthemap,Ndir,theinitialredshiftatwhichraytracingbegins,zin,thedistancebetweenevaluationsonthegeodesic∆ps,theangularresolutionof theCMBlensedmaps∆ang,andthepreferreddirectionused(seeTable2). – 10 – constructed CMB maps for each simulation. 3.2. Ray-tracing technique 3.2.1. Lensing regimes As in our previous work (Ant´on et al. 2005) we divide the total lensing effect into three parts: • AWL (A weak lensing), namely the effect due to scales k > 2π/L (where max L = 42h−1 Mpc) at redshifts z < 6. This signal is dominated by strongly nonlinear max scales • BWL, the lensing signal due to scales k < 2π/L which corresponds to modes that max are always in the linear regime down to z = 0 • CWL, the lensing signal due to scales k ≥ 2π/L but at redshifs z > 6 max The goal of this paper is to calculate the AWL signal, using ray-tracing through N-body simulations, to ℓ values in the range 1000−10000. Since the modes associated with the BWL component are specifically chosen to correspond to scales that are linear down to z = 0, which is set by wavenumbers k . 0.15h Mpc−1 (Smith et al. 2003), the BWL component can be calculated with the linear approach implemented in CMBFAST. The CWL component involves modes in the mildly nonlinear regime which can, nonetheless, be evolved via approximation schemes (Zel’dovich 1970; Sandarin & Zel’dovich 1989; Moutarde et al. 1991). Standard semi-analytical methods designed to study weak lensing in the nonlinear regime should also apply in this case, hence, the CWL effect can be calculated without resorting to N-body techniques. This calculation can be performed using the nonlinear version of CMBFAST which is based on semi-analytical approaches