ICECOLA: fast simulations for weak lensing observables
Abstract
Approximate methods to full Nbody simulations provide a fast and accurate solution to the development of mock catalogues for the modeling of galaxy clustering observables. In this paper we extend ICECOLA, based on an optimized implementation of the approximate COLA method, to produce weak lensing maps and halo catalogues in the light cone using an integrated and self consistent approach. We show that despite the approximate dynamics, the catalogues thus produced enable an accurate modeling of weak lensing observables one decade beyond the characteristic scale where the growth becomes nonlinear. In particular, we compare ICECOLA to the MICEGC Nbody simulation for some fiducial cases representative of upcoming surveys and find that, for sources at redshift , their convergence power spectra agree to within one percent up to high multipoles (i.e., of order ). The corresponding shear two point functions, and , yield similar accuracy down to and arcmin respectively, while tangential shear around a lens sample is accurate down to arcmin. We show that such accuracy is stable against an increased angular resolution of the weak lensing maps. Hence, this opens the possibility of using approximate methods for the joint modeling of galaxy clustering and weak lensing observables and their covariance in ongoing and future galaxy surveys.
keywords:
methods: numerical – dark matter – largescale structure of Universe1 Introduction
The images of distant sources are distorted as their photons travel across the gravitational field of darkmatter structures along the line of sight of observers today. The strength of the phenomenon depends on the expansion history and the growth of structure (for reviews see Bartelmann & Schneider 2001; Kilbinger 2015). Precise measurements of this effect can tightly constrain cosmological models in a complementary way to other traditional techniques (see e.g. Weinberg et al. 2013), and elucidate the properties of the accelerated expansion and the largescale structure formation process, which are directly related to the properties of dark energy and dark matter, respectively. Present and planned galaxy surveys such as the Dark Energy Survey (DES)^{1}^{1}1http://www.darkenergysurvey.org/, the KiloDegree Survey (KiDS)^{2}^{2}2http://kids.strw.leidenuniv.nl, the Hyper SuprimeCam Subaru Strategic Program (HSCSSP)^{3}^{3}3http://hsc.mtk.nao.ac.jp, Euclid^{4}^{4}4http://www.euclidec.org/, and the Wide Field InfraRed Survey Telescope (WFIRST)^{5}^{5}5https://wfirst.gsfc.nasa.gov, the Large Synoptic Survey Telescope (LSST)^{6}^{6}6http://www.lsst.org, represent a large increase in the volume sampled by weak gravitational lensing experiments, which results in a tremendous narrowing down of the statistical uncertainties associated with the measurements. This has to go along with an equivalent reduction of systematic errors, otherwise the constraining power on cosmological models will not scale accordingly.
The process of converting the data acquired by experiments into catalogues ready for scientific usage involves many calibration and analysis steps. It is a nonlinear process in which many systematic effects may be introduced, from instrumental noise and limitations to atmospheric variations (Leistedt et al., 2015) or astrophysical sources of errors (e.g., the confusion of objects in the same region of the sky, see Hartlap et al. 2011). The science analysis is also prone to introduce biases if the physical properties of the sample observed (such as the signal measured and its covariance matrix) are not well understood and modeled accurately (Dodelson & Schneider, 2013; Taylor et al., 2013; Blot et al., 2015; Sellentin & Heavens, 2016, 2017). In recent years, mock galaxy surveys have become an essential tool to model and mitigate systematics, as well as optimizing analysis pipelines (e.g., in the Baryon Acoustic Oscillations reconstruction, see Takahashi et al. 2009; Manera et al. 2013; Kazin et al. 2014; Ross et al. 2015).
Synthetic data has been traditionally produced by means of body simulations, which require huge computational resources in order to meet the resolution and volumes sampled by galaxy surveys (Kim et al. 2011; Angulo et al. 2012; Alimi et al. 2012; Skillman et al. 2014; Heitmann et al. 2014; Fosalba et al. 2015b; for a review see Kuhlen et al. 2012). In this context, weak lensing can be accurately modeled in simulations via raytracing techniques, which follow the ray propagation from the source to the observer along the perturbed path (see e.g. Blandford & Narayan 1986; Jain et al. 2000; Das & Bode 2008; Teyssier et al. 2009; Li et al. 2011). This involves intensive computations because the deflection angle needs to be constantly updated as the ray travels in order to determine the geometry at each encounter of the multiplelens system. However, under the Born approximation, which we use in this work, integrations are much faster since photons propagate along the unperturbed path across each of the (uncorrelated) lens planes (see White & Hu 2000; Fosalba et al. 2008; Kiessling et al. 2011; Borzyszkowski et al. 2017). This removes higherorder contributions and the coupling of lenses at different distances (Krause & Hirata, 2010), and amount to subpercent corrections for the relevant scales and redshifts sampled by planned surveys and thus can be safely neglected in their modeling (see e.g, Petri et al. 2016).
Large ensembles of realizations are necessary for a precise determination of covariance matrices and the optimization of some pipelines and this becomes computationally prohibitive using conventional numerical simulations. This problem has been overcome thanks to the emergence of techniques known as approximate methods, that provide an optimal solution in the compromise between accuracy and speed. The latter comes from avoiding solving the highly nonlinear dynamics in virialized regions, that are expensive to resolve. Using an approximate gravity solver to evolve the dark matter density field, galaxies or haloes are modeled either using a biasing prescription (Coles & Jones 1991; White et al. 2014; Chuang et al. 2015a; Kitaura et al. 2015) or by identifying collapsed regions in the density field (Monaco et al., 2002; Scoccimarro & Sheth, 2002; Tassev et al., 2013; Manera et al., 2015; Munari et al., 2016). For a comparison of the performance of some of these methodologies see Chuang et al. (2015b) and for a recent review see Monaco (2016).
Among the suite of approximate methods recently developed, COLA (Tassev et al., 2013, 2015) stands out as a simple yet powerful alternative to exact Nbody simulations. In particular, it accurately solves mildly nonlinear scales by integrating numerically the equations of motion of particles in a frame comoving with respect to Lagrangian Perturbation Theory (LPT) trajectories. By using broad time steps and solving the gravitational forces using a fine ParticleMesh (PM) algorithm, the spatial distribution of darkmatter halos can be faithfully captured whereas the computational cost is reduced by 23 orders of magnitude with respect to standard body methods. In Izard et al. (2016) an optimal configuration of the method, named ICECOLA, was presented for the production of accurate mock halo catalogues adapting the parallel version of COLA developed in Koda et al. (2015). In this paper we extend the ICECOLA code to produce allsky light cone catalogues onthefly^{7}^{7}7At run time of the simulation., what allows for the fast and accurate modeling of weak lensing observables. Other recent efforts involving COLA have been developed to extend the capabilities of the original code, including Modified Gravity solvers and massive neutrinos (Howlett et al., 2015; Valogiannis & Bean, 2016; Winther et al., 2017; Wright et al., 2017) or relativistic corrections (Borzyszkowski et al., 2017). Plain PM simulations have also been useful to produce mock catalogs (Merz et al., 2005; White et al., 2010), as well as variations of the algorithm using few time steps (White et al., 2014; Feng et al., 2016).
Fast modeling of weak lensing is a specially challenging task, since the gravitational evolution of the relevant scales is well within the nonlinear regime of structure formation. At the same time, it entails probing large volumes extending far in the observer past light cone, up to the most distant sources (e.g, the last scattering surface for CMB lensing). In this paper, we discuss the capabilities and limitations of the ICECOLA method for light cone observables, and show that it can meet these two basic requirements. To do so, we implement the production of matter and halo light cone catalogues onthefly. The dark matter density field is then projected into twodimensional sky maps, from which we can construct weak lensing maps in the Born approximation, following the procedure presented in Fosalba et al. (2008) and Fosalba et al. (2015a). Halo catalogues in the light cone are stored and can then be used to produce galaxy mocks based on the Halo Occupation Distribution approach (Carretero et al., 2015; Crocce et al., 2015; Fosalba et al., 2015a). Weak lensing properties are then imprinted into the halos (or galaxies) interpolating them from the weak lensing maps.
Other fast methods to produce weak lensing maps have been recently proposed in the literature. Giocoli et al. (2017) is based on a halo model formalism (see references therein too) and Yu et al. (2016) use a gaussianization technique. These works differ to what is presented here in that they are not self consistently built upon an Nbody solver of the gravitational dynamics.
This paper is organized as follows: Section 2 presents the light cone construction, its outputs and postprocessing pipeline. Section 3 describes the ICECOLA and MICE Nbody simulations used for the analysis. The methodology for computing allsky weak lensing observables and their analysis is discussed in section 4, whereas the corresponding results for the light cone halo catalogues is the subject of section 5. We conclude summarizing our main results in section 6.
2 Allsky light cone construction
Sources in the Universe are seen by observers today as they were at the time photons from the source were emitted. The time elapsed corresponds to what it takes photons to travel from the source to the observer. This “light cone” geometry can be implemented in a simulation by computing the light cone crossing time of each particle, determined by the intersection of its trajectory with the surface of a sphere that is centered at the observer and its physical radius shrinks at the speed of light (reaching zero at the present time).
In an body simulation, particle positions are sampled in a collection of discrete time steps. In the socalled leapfrog integration, each time step contains a pair of drift and kick operators that integrate the equations of motions and update the positions and velocities of particles respectively. The two operators are interleaved, that is, they are separated in time by half the size of the time step (see Quinn et al. 1997 and Dehnen & Read 2011 for a more detailed description of the leapfrog integration). We build the particle light cone output as the set of interpolated particle coordinates determined at their light cone crossing time. In such way, the time evolves smoothly in the radial coordinate and there are no jumps in the catalogues, as may occur in light cones built from stitching snapshots. In ICECOLA, the light cone routine is called after each drift and kick operators, producing a shell with the geometry sketched in Fig. 1. The algorithm takes the following steps:

The time of the last drift and kick operators to evolve particle positions and velocities ( and respectively) determine the radius of the shell being constructed at a given time of the simulation^{8}^{8}8Particle coordinates at intermediate time values between the last drift and kick operator can be interpolated using an extra pair of drift and kick operators..

The shell is enlarged with a velocity buffer zone around to consider moving particles that may enter into the shell during the interpolation process^{9}^{9}9There is no buffer necessary around since at that region the crossing time is close to the time of the drift operator and therefore the displacements during the interpolation tend to zero..

Particles that do not belong to the volume of the shell, including the buffer zone, are discarded.

Given the position and radial velocity of each particle, its crossing time is computed as:
(1) where is the time when the light cone boundary is at the distance of the particle position at the last drift operator.

Particle coordinates are interpolated at crossing time and checked whether they are within the volume of the shell (without the buffer zones now).

The particles satisfying the previous conditions are copied to the light cone catalogue. The process may be repeated for different box replicas, in which the coordinates are shifted by the box side length around the observer, as we discuss later in this section.
The light cone routine is called after each pair of kick and drift operators. To make it more efficient, we precompute and tabulate all the time dependent functions and do an interpolation when these have to be evaluated. This is the case for conversions between cosmological distances and times and the timevarying quantities in the temporal operators that interpolate particle coordinates to their crossing times. We have estimated that with as few as 10 time samples within each timestep the accuracy achieved is sufficient, since they are smoothlyvarying quantities. In particular, we have checked that the errors in particle coordinates arising from these approximations are much lower than the spatial resolution of the simulation.
Given that the simulated comoving volume has periodic boundary conditions, it makes adjacent box replicas used to construct the full light cone around the observer have a continuous density field distribution across the entire light cone volume. We take advantage of this to produce light cones sampling larger volumes (allowing to produce catalogues that extend to higher redshifts and cover more sky area). The disadvantage of this approach is that repeated structures along the line of sight could potentially appear in the catalogue. However, they will be separated a distance given by the box size and they have typically crossed the light cone boundary at a different cosmological time, since their distance to the observer will likely be different. Therefore, no obvious repetitions of structures will be observed in the resulting light cone volume (see also Fosalba et al. 2008, 2015a).
Writing light cone catalogs containing the particle information results in huge data volumes, that demand not only large storage quotas but also require large computing time to be processed because of the limited input/output bandwidth. In the scope of a fast method, it is desirable to generate highlevel processed catalogs that can be postprocessed more easily. We describe next the two types of catalogs we generated onthefly.
2.1 Onthefly projected matter density fields
We implement and adapt to the particularities of COLA the “onion universe” method originally developed by Fosalba et al. (2008) to generate compressed outputs of the matter density field in the light cone that are used to derive weak lensing maps in postprocessing. The particle distribution in the light cone is interpolated into a threedimensional grid using spherical coordinates around the observer. The radial coordinate is binned into several thin shells and the sky area is discretized into twodimensional pixels using the Healpix software (Górski et al., 2005). In that way, the volume is subdivided into several slices of Healpix maps, whose values give the darkmatter particle count in a given pixel.
For the ease of comparison with MICEGC reference Nbody simulation, we follow Fosalba et al. (2015b) and use 265 shells in the redshift range having a width of at low (high) redshift, corresponding to steps of in lookback time. It is left for future work the study of how many radial bins are enough to recover accurately the matter field within the accuracy of COLA. The Healpix maps have an angular resolution given by the parameter , which produces million pixels of length. The data volume for storing all the maps is 50 GB, which represents a compression factor of two orders of magnitude with respect to the whole particle information generated by the light cone output.
2.2 Onthefly halo catalogues
The COLA version of Koda et al. (2015) includes a parallel FriendsofFriends (FoF) algorithm (Davis et al., 1985) that runs onthefly on particle snapshots, and it is based on the publicly available serial code from the body shop of the University of Washington^{10}^{10}10http://wwwhpcc.astro.washington.edu/tools/fof.html. In ICECOLA we adapted it to produce FoF catalogues in the light cone.
A halo finder explores the environment around particles and links those that belong to the same halo. In a comoving snaphsot, the environment of regions near the edges of the simulated box can be completed with the particles close to the opposite side of the box, thanks to the periodic boundary conditions. However, this is no longer possible when the algorithm is run in the light cone particle distribution, since the crossing time varies across the volume and thus differs for each box replica. Therefore, at the edges of the box it is needed the information of the neighboring replica, which has to be computed using the adequate crossing times. That is, the environment of a particle close to the box edge has to be completed with the volume corresponding to the neighboring replica. But since only one replica is processed at a time, it is necessary to add buffer zones around box edges and provide in that way the required spatial environment to the algorithm.
Besides, each call to the light cone routine builds a spherical shell of the light cone volume. Therefore, buffer zones are also needed around the spherical shell boundaries. In summary, when FoF catalogues are requested, the light cone algorithm builds a subvolume containing the shell and the buffer zones^{11}^{11}11Note that it shall add an additional buffer for the velocities, as described before.. Then the FoF code is run in the corresponding subsample of particles and finally only the halos whose center of mass lay within a given spherical shell and the local box replica are added to the catalogue.
An incorrect implementation of buffer zones could lead to an incorrect detection of haloes, especially with the largest ones, that would be fragmented into small pieces. We set the width of the buffer zones to . We try to check for the fragmentation effect by comparing the mass function in the light cone and in snapshots at redshift 0.5 and 1. We make a subsample of halos in a spherical shell centered at these redshifts and having a volume equal to the volume of the box (this can be achieved by choosing the appropriate width of the shell). In this way we expect that the measurements have similar error bars. In Fig. 2 we show the ratio of the mass functions measured in the light cone and in the snapshot. Differences are well within the percent level at low masses and compatible with error bars at high masses. We also check that varying the buffer size does not modify our results. We thus conclude from this analysis that there are no signs of systematic biases in the algorithm for generating FoF in the light cone.
2.3 Postprocessing pipeline
The darkmatter light cone catalogues presented in this section are generated onthefly and ease considerably the postprocessing steps needed to produce science ready catalogues. However, an additional processing step is required to produce weak lensing observables. One has to transform darkmatter counts maps into weak lensing ones and combine them with the light cone halo catalogues. We do this in a postprocessing pipeline that we briefly describe below^{12}^{12}12For these computations we employ the Healpy software tools, https://healpy.readthedocs.io/en/latest/ and that follows the approach presented in Fosalba et al. (2008) (see also Fosalba et al. 2015a). For further details and the assumptions involved we refer the reader to Sec. 4 and 5, where results are also discussed. The basic steps of the weaklensing postprocessing pipeline are:

The convergence field is determined in the Born approximation by integrating along the lineofsight the 2D dark matter counts maps weighted by the lensing kernel (eq. 5).

We transform back to angular space to obtain the shear map. Note that these transformations require the weak lensing maps to be allsky to avoid boundary effects. All these steps are repeated for each of the thin shells that sample the whole redshift range , having at the end a suite of 2D maps (the “onion universe” as introduced in Fosalba et al. 2008) of the weak lensing effect based on the Healpix discretization.

Finally, halos are assigned convergence and shear values. To do so, we determine the pixel in the 2D map, for a given redshift shell, that corresponds to the halo position and assign its weak lensing quantities to the halo.
3 Simulations
In this section we describe the numerical simulations used in this paper to validate the performance of ICECOLA for the modeling of weaklensing observables.
3.1 IceCola simulations
In our previous work Izard et al. (2016) we investigated the impact of various parameters of the COLA method on the matter power spectrum, the halo abundance and clustering. We showed that, in order to capture halo formation accurately at a given output time, more than 10 simulation timesteps need to be computed from the initial conditions to that given output time. We suggested a configuration using 40 time steps linearly distributed with the scale factor between an initial redshift of and and using a cell size for the computation of forces that is 3 times smaller than the mean interparticle distance separation. Simulations using this configuration predict the mass function with a accuracy and the matter power spectrum within the percent level at scales . In this paper we use this optimal configuration for runs with particles and a box size of , which results in a particle mass of . Both the particle mass, the cosmological model and the input linear power spectrum are identical to the reference MICEGC Nbody simulation (see Sec. 3.2).
For this paper we make use of the capabilities of ICECOLA to generate an allsky light cone onthefly. The simulated box is replicated twice in each cartesian direction^{13}^{13}13This generates 64 box replicas in total., allowing to build a light cone that samples the redshift interval . From the full particle information, two derived catalogues are produced and stored: twodimensional matter density maps projected on the sky (see Sec. 2.1) and FriendsofFriends (FoF) halo catalogue (Davis et al., 1985) with a linking length of 0.2 (see Sec. 2.2). The allsky FoF catalogue include objects with 50 or more particles, generating 15 Gb of data. This data volume is to be added to the 50 Gb of the particle darkmatter light cone outputs in Healpix maps. The simulation was run on 1024 cores in 100 minutes, where roughly 50 of the time was spent running the FoF algorithm in the light cone volume covering the 64 box replicas.
Our default Healpix resolution for the weak lensing maps is . However we produced a higher resolution version of them to test in Sec. 5 whether the accuracy of the method is limited by the approximated dynamics of COLA or by the pixel scale employed. These maps have , which correspond to pixels, and cover only one octant of the sky to reduce the total data volume. We generated them at postprocessing from a light cone catalogue containing threedimensional particle information, which we projected at postprocessing into twodimensional maps as explained in Sec. 2.1.
3.2 MICEGC: the benchmark body run
Our benchmark body simulation is the MICE Grand Challenge (MICEGC)^{14}^{14}14More information is available at http://www.ice.cat/mice.. It evolved particles in a volume of using the gadget2 code (Springel, 2005) assuming a flat CDM cosmology with , , , , and . This results in a particle mass of . The initial conditions were generated at using the Zel’dovich approximation and a linear power spectrum generated with camb^{15}^{15}15http://camb.info. Haloes were identified using a the FoF algorithm (Davis et al., 1985) with a linking length of 0.2. Note that MICEGC and the ICECOLA runs have different box sizes and therefore the comparison between them will be affected by sample variance.
This simulation and its products have been extensively validated. Of special interest in this work is the modeling of weak gravitational lensing, presented in Fosalba et al. (2015a). The darkmatter and halo outputs are described in Fosalba et al. (2015b) and Crocce et al. (2015). The Halo Occupation Distribution implementation used to produce galaxy mocks is detailed in Carretero et al. (2015) while Hoffmann et al. (2015) focuses on the higherorder clustering. In this work we make use of the weak lensing maps and halo catalogues in its past light cone.
The weak lensing maps in MICEGC employ a finer pixel size of (i.e, Healpix ), meaning that they have a better angular resolution than the catalogues of ICECOLA, but we have checked that degrading to the same resolution than ICECOLA does not change results down to . As for the light cone outputs, the convergence and matter density maps are allsky, whereas the corresponding MICEGC halo catalogues are produced only for one octant of the sky.
4 Modeling weak lensing
As mentioned in Sec. 2.1, we make use of the twodimensional projected matter density field generated onthefly to compute weak lensing maps in the Born approximation. The resulting convergence maps receive contributions from the intervening matter fluctuations weighted by the socalled lensing efficiency, which peaks around half the distance to the source. Throughout this paper we choose a particular configuration to study weak lensing, in which the sources are centered at redshift 1 and the lenses at redshift 0.5.. It is thus informative to investigate the accuracy in modeling the darkmatter clustering at the latter distance to assess the accuracy of ICECOLA to model weak lensing observables. To illustrate this, in Fig. 17 we compare the angular matter power spectrum at for ICECOLA (circles) against MICEGC (squares), where we have subtracted the shotnoise contribution to the power spectrum according to
(2) 
where is the number of objects in the sample and is the fractional area of the sky of the catalogue (1 if it is allsky). At small scales it represents a 1% correction for COLA. We display Gaussian error bars, calculated by considering the number of modes sampled at each multipole bin
(3) 
where is the multipole bin width. In the lower panel displaying the ratio, we add the respective errorbars in quadrature.
We find that the agreement between ICECOLA and MICEGC is very good up to and within at . At and in the flatsky limit (which should be very accurate for the relevant scales) these multipoles project onto and respectively, where we have used that for the cosmology assumed. We note that, at these scales, most of the contribution to the power spectrum comes from the internal halo structure (van Daalen & Schaye, 2015), and therefore we expect a decline in the power for COLA since dynamics within haloes are not fully resolved (see also the matter power spectrum of the same simulation in Fig. 8 of Izard et al. 2016).
4.1 Convergence maps
The convergence field characterizes the fractional change in the image size. In the Born approximation, the contribution from all the lens planes between the observer and the sources at a distance can be calculated by an integral:
(4) 
where is the Hubble constant, is the matter density parameter, is the speed of light, is the scale factor, and is the 3D matter density contrast at the angular position and radial distance . Following Fosalba et al. (2008), we compute the convergence maps by discretizing the equation in several radial bins and angular pixels as
(5) 
where is the width of the radial bins. The color map in Fig. 4 shows the convergence field in a patch of , which features the characteristic high density peaks of the cosmic web surrounded by large void regions. The white ticks display the shear field, which we discuss in more detail below (see Sec. 4.2).
It is useful to derive an expression for the convergence angular power spectrum starting from Eq. 4. Using the Limber approximation (Limber, 1953) to compute the twopoint statistics we can write (Kaiser, 1992, 1998):
(6) 
where is the 3D matter power spectrum which can be derived from a theoretical model or measured from a numerical simulation. In the latter case, we can replace the integral by a sum over the redshifts where measurements for the matter power spectrum (produced onthefly after each time step) are available. The shotnoise contribution to the convergence angular power spectrum, , can be computed using Eq. 6 by integrating the 3D shotnoise power spectrum, , where is the mean number density of particles, weighted by the weaklensing kernel.
We show the measurements of the convergence angular power spectrum at a source redshift of in Fig. 5. The curves display different theoretical predictions obtained with Eq. 6. In particular, the dotted curve uses the nonlinear 3D matter power spectrum from ICECOLA. This prediction is in close agreement with the actual measurements from the maps, although it underestimates the signal at small scales^{18}^{18}18Note that at large scales, the deviations of the prediction using the ICECOLA nonlinear matter power spectrum are due to sample variance.. This is because the measured 3D power spectrum used in equation 6 is available up to a maximum wavenumber. As a consequence, the prediction doesn’t account for the contributions above that wavenumber, which are relevant in the limit of short distances and high multipoles , that correspond to high . When compared to our fiducial MICEGC reference measurements, COLA recovers the power accurately up to and declines thereafter. The difference is at . Note that theses are highly nonlinear scales, one decade above the transition scale to the nonlinear regime (i.e., where the nonlinear growth starts being significant, at ), and roughly correspond to at the distance where the lensing efficiency peaks.
The convergence field is largely correlated with the intervening mass distribution due precisely to the lensing effect^{19}^{19}19This effect is closely related to what is known as galaxygalaxy lensing. This is shown in Fig. 6 for the crosscorrelation between the convergence at and a lens sample of the matter density field in a lens plane at with width . In this case, a good agreement is observed between ICECOLA and MICEGC up to and the deviations are within 10% at most. Thus the crosscorrelation with the matter field is better reproduced than the autopower spectrum of the convergence. While the former measures the contribution to the convergence by the matter density field at the lens position, the latter has contributions also from high modes of structures at low redshift, which are more nonlinear. Therefore, at a given angular scale, the convergence autopower spectrum has more nonlinear contributions, hard to model in COLA, explaining why differences are larger in Fig 5 than in Fig. 6 and the latter resembles the behaviour seen in Fig. 17.
4.2 Shear maps
The shear field, , is a spin2 tensor that in the complex notation can be expressed in terms of the polarization Stokes parameters (): . In the allsky formalism, the decomposition in harmonic space is written in terms of the spin2 spherical harmonics (Hu, 2000)
(7) 
The shear can be decomposed in the socalled  and modes as In the weak gravitational lensing limit and in the Born approximation, the mode shear is zero, , due to the parity symmetry of the weak lensing field (see e.g. Cooray & Hu 2002). Then the shear is related to the convergence in harmonic space as (see Hu 2000):
(8) 
Note that in the small angle limit () the coefficient is just 1 and the power spectra of both quantities coincide.
Finally, the Stokes parameters of the shear field are given by (in the sign convention adopted by Healpix and in the absence of modes):
(9)  
(10) 
with defined as: .
In summary, starting from an allsky convergence map we compute the shear by {enumerate*}
transforming the map to harmonic space,
computing the mode shear using equation (8),
determining the two components in angular space of the shear by the inverse harmonic transforms given by equations (9) and (10) ^{20}^{20}20The inverse harmonic transforms are performed using Healpy, https://healpy.readthedocs.io/en/latest/.
Fig. 4 displays the resulting shear map (white ticks), overlaid on the convergence field (scalar field, color coded). The shear is tangentially oriented to the convergence peaks, as expected for a cosmological weak lensing field (i.e, no mode, which would generate swirl or rotated patterns). An extensive validation of the shear is presented in the next section, where we also discuss twopoint clustering statistics for the light cone halo catalogues.
5 Halo catalogues for weak lensing and clustering
Although we concentrate here on halo catalogues, we can think of them as a simple galaxy catalogues populated only with central galaxies that live in the center of darkmatter halos. In the last step of the postprocessing pipeline, we combine the two types of light cone outputs of the simulations: the halo catalogues and the weak lensing maps derived from the projected matter density field. In particular, the position of each halo determines the pixel in the 3D volume discretizaton of the light cone (as explained in Sec. 2.1) that is used to assign the convergence and shear values. For illustrative purposes, we generate two halo samples at well separated redshifts: the source sample at and the lens sample at . Both catalogues are allsky, have a radial width of and include objects with 100 or more particles (i.e., halos of mass larger than ). This results in 6M and 11M objects for the lens and source samples, respectively, which correspond to number densities of 0.3 and . The same criteria are used to draw the reference body samples from MICEGC, except the latter were limited to only one octant of the sky.
In this section we use an additional ICECOLA simulation (with the same exact parameters and initial conditions as the default one) with the aim to distinguish the effects of the approximate dynamics of the simulation from the limited angular resolution employed in the Healpix sky maps. For that, we stored a light cone catalogue with particle information^{21}^{21}21The overhead due to writting particles in the catalogues was a 10% of the running time for storing one octant of the light cone. that we used to construct at postprocessing a higher resolution version of the 2D projected matter density field^{22}^{22}22The light cone and therefore the catalogues generated are restricted to one octant for this simulation. However, this is sufficient to appropriately sample all the relevant scales discussed in this paper.. We recall that the maps produced onthefly use the Healpix parameter while for this test we use . In turn, MICEGC catalogues were built with , but have checked that degrading them to makes no difference on scales , which is well within the range of scales discussed in this paper. Thus for all practical purposes there is no difference between the highresolution maps of ICECOLA and MICEGC. The reference we quote the corresponding pixel size for these resolutions, 1.72, 0.85 and 0.43 arcmin, for and respectively. Lastly, using these higher resolution Healpix maps we run the postprocessing pipeline and produce the corresponding halo catalogues with weak lensing properties.
5.1 Convergence autocorrelation
The weak lensing maps describe the image distortions that an object would experiment at a given location of the sky. In practice, however, the maps are only sampled at those positions where there are sources (halos or galaxies). Nevertheless, this does not affect the 2point weak lensing statistics of the halo lens sample, that are assigned from the weak lensing maps built from the underlying darkmatter field, and the measurement is also independent of the source sample galaxy bias.
Fig. 7 shows the convergence correlation function in angular space for the source sample at . In this work we measure correlation functions using the public code Athena^{23}^{23}23Athena is based on a twodimensional tree, http://www.cosmostat.org/software/athena/.. We find that ICECOLA reproduces the reference MICEGC measurements at large scales and down to separations of arcmin. At smaller scales, ICECOLA shows a lack of power regardless of the angular resolution of the weak lensing maps, although this effect is smaller for the higher resolution maps (e.g., at 2 or 1 arcmin for or , respectively). Similarly as for the convergence power spectrum (see Fig. 5), it is remarkable that the signal is well reproduced to separations one decade smaller than the scale where nonlinearities arise. These results are in agreement with Heitmann et al. (2010) (see their Fig. 1 and also sec. 3.3 in Fosalba et al. 2015a), where it is investigated how a systematic effect in the matter power spectrum translates into a suppression of the shear correlation function at small scales. In the simulations used in this work, COLA recovers 50% of the power in the matter power spectrum at , which according to Heitmann et al. (2010) should translate into a percent level accuracy down to scales of arcmin in the convergence angular correlation function, as we find.
The accuracy in the convergence twopoint statistics is in rough agreement with the corresponding one from the angular power spectrum. Since angular scales of arcmin are approximately associated with multipoles in the Limber limit, we should expect a similar accuracy in both statistics using this 1to1 relation in scales. However we find a better accuracy for the correlation function () than for the ’s (). A plausible explanation is that the approximate dynamics in ICECOLA introduce an error that resembles a scattering of the position of particles with a characteristic scale. Clustering measurements in configuration space will be correct at separations equal or larger than that scale. However, when transforming to harmonic space, small and large scales are mixed and errors are propagated to a wider range of scales (see e.g., Monaco et al. 2013).
5.2 Shear correlations
The two components of the shear, that in the polar coordinates read , can be expressed in a new and rotated base that is aligned with the separation vector between two points (where one point is associated with the shear map and the other with the positions of the objects in the lens sample, see e.g. Kilbinger 2015). The tangential and the cross components of the shear are then defined as:
(11)  
(12) 
where is the polar angle of the separation vector. The nonvanishing shear correlation functions are usually expressed in the two following combinations of the tangential and cross components (see e.g. Bartelmann & Schneider 2001)
(13) 
In the absence of shear mode and under the Limber approximation, the shear correlation functions can be expressed as (see Bartelmann & Schneider 2001):
(14) 
where is the Bessel function of the first kind.
Fig. 8 shows the shear correlation functions for the source sample at . On large angles, ICECOLA is in good agreement both with the MICEGC simulation and nonlinear theory expectations. The only caveat being the excess power in for both simulations with respect to theory. One possible explanation is that we are using the smallangle limit for the prediction. Another is that the volume is not enough and both measurements are slightly biased (note that the points are very correlated). Nonetheless in this paper we are concerned about the agreement between ICECOLA and MICEGC, and that is well achieved in Fig. 8. On small scales, ICECOLA is percent accurate, as compared to MICEGC, down to and arcmin for the and components, respectively. Such different “characteristic accuracy scale” (i.e, the smallest angular scale at which ICECOLA matches Nbody results) comes from the fact that the minus component probes smaller scales of the density field, which are more nonlinear, than the plus component. We can seen it in eq. 14, since the Bessel function peaks at larger values of its argument than . In both cases ICECOLA resolves accurately scales down to one order of magnitude smaller than the nonlinear scale (i.e., the scale where nonlinearities become relevant), as we also found for the convergence correlation function. Also clear from Fig. 8 The angular resolution of the lensing maps used does not alter this characteristic scale. But the coarser the angular resolution of the maps, the larger is the power loss at that characteristic scale, as shown in the lower panels of Fig. 8.
Another important observable for weak lensing is the tangential shear around galaxies, also known as the galaxygalaxy lensing. The correlation function of the tangential shear in the flatsky limit, which is accurate enough for separations angles below few degrees (see de Putter & Takada 2010), it is given by
(15) 
where is the crosspower spectrum between the source convergence and the lens number density. Fig. 9 shows the measurements of the tangential shear correlation function in ICECOLA for the source and lens samples at and respectively for halos with . There is a good agreement with MICEGC for scales larger than 4 arcmin, regardless of the angular resolution employed for the weak lensing maps. Again, as for the previous lensing observables discussed, this matches the “characteristic accuracy scale” of ICECOLA, i.e, one decade below the nonlinear scale. In the theoretical curves, the matter power spectrum has been multiplied by the linear bias factor that matches the clustering amplitude. The disagreement found on small scales with respect to theory predictions (which assume linear halo bias) may be attributed to a (positive) nonlinear halo bias, that increases the clustering amplitude at small scales. The tangential shear depends not only on the weak lensing properties of the background sample but also on the clustering of the foreground one. To discriminate the impact of both on the tangential shear, we show in appendix A the angular correlation function of the lens halo sample. We conclude that at small scales, the errors from the halo bias seem to be subdominant in Fig 9, although halo autocorrelations start to depart from those in MICEGC also at about arcmin.
6 Conclusions
In the advent of next generation of galaxy surveys, mock galaxy catalogues are key in the optimization of survey strategies and analysis pipelines. The increasingly large range of scales probed by weak lensing experiments put stringent requirements on the simulations needed to generate mock catalogues. We presented ICECOLA, a methodology to produce efficiently weak lensing maps and halo catalogues, both in the light cone geometry and generated onthefly. The speedup thus gained with respect to standard body simulations is of 23 orders of magnitude and it requires a minimal postprocessing of the catalogues compared to other techniques based on concatenation of snapshots to produce light cones. In this way, the amount of data to be stored is also greatly reduced.
We demonstrated the accuracy of the method by analysing the weak lensing properties of allsky halo samples and compared them to the fiducial MICEGC simulation. Running ICECOLA with the optimal parameter configuration described in Izard et al. (2016), we show that we can accurately model the twopoint clustering of weak lensing observables one decade beyond the characteristic scale where the growth becomes nonlinear. In particular, the convergence angular power spectrum is consistent with MICEGC to multipoles of 1000, or 23 arcmin for the correlation function. For a fiducial case with sources at and lenses at , we show that the shear correlation functions are accurate down to 2 and 20 arcmin for the and components respectively, whereas the tangential shear is accurately modeled down to 4 arcmin. These “characteristic accuracy scale” where ICECOLA breaks down is due to the approximate dynamics. We tested that our choice for the angular resolution of the weak lensing maps is balanced since improving it by a factor of two does not affect the scale where ICECOLA deviates, although it makes the power loss with respect to exact Nbody simulations smaller.
Our implementation of ICECOLA opens the possibility of using approximate methods for the joint modeling of galaxy clustering and weak lensing observables and their covariance in the optimal exploitation of future galaxy surveys.
Acknowledgments
We thank Linda Blot for her great help in some of the ICECOLA runs used in this paper and for helpful comments along the completion of the project. We are also grateful to the DEUS consortium (JeanMichel Alimi, PierStefano Corasaniti, and Yann Rasera) and Alina Kiessling for their help and infrastructure support. AI was supported in part by Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. AI was also supported in part by NASA ROSES 13ATP130019, NASA ROSES 14MIROPROs0064, NASA ROSES 12EUCLID120004, and acknowledges support from the JAE program grant from the Spanish National Science Council (CSIC). MC has been funded by AYA201344327 and AYA201571825P. MC acknowledges support from the Ramon y Cajal MICINN program. PF is funded by MINECO, project ESP201348274C31P, ESP201458384C31P, and ESP201566861C31R. Simulations in this paper were run at the MareNostrum supercomputer  Barcelona Supercomputing Center (BSCCNS, www.bsc.es), through grants AECT201510013, 201520011, 201620010, and 201630015. Funding for this project was partially provided by the Spanish Ministerio de Ciencia e Innovacion (MICINN). Copyright 2017. All rights reserved.
References
 Alimi et al. (2012) Alimi J.M., et al., 2012, preprint, (arXiv:1206.2838)
 Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
 Bartelmann & Schneider (2001) Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
 Blandford & Narayan (1986) Blandford R., Narayan R., 1986, ApJ, 310, 568
 Blot et al. (2015) Blot L., Corasaniti P. S., Alimi J.M., Reverdy V., Rasera Y., 2015, MNRAS, 446, 1756
 Borzyszkowski et al. (2017) Borzyszkowski M., Bertacca D., Porciani C., 2017, preprint, (arXiv:1703.03407)
 Carretero et al. (2015) Carretero J., Castander F. J., Gaztañaga E., Crocce M., Fosalba P., 2015, MNRAS, 447, 646
 Chuang et al. (2015a) Chuang C.H., Kitaura F.S., Prada F., Zhao C., Yepes G., 2015a, MNRAS, 446, 2621
 Chuang et al. (2015b) Chuang C.H., et al., 2015b, MNRAS, 452, 686
 Coles & Jones (1991) Coles P., Jones B., 1991, MNRAS, 248, 1
 Cooray & Hu (2002) Cooray A., Hu W., 2002, ApJ, 574, 19
 Crocce et al. (2015) Crocce M., Castander F. J., Gaztañaga E., Fosalba P., Carretero J., 2015, MNRAS, 453, 1513
 Das & Bode (2008) Das S., Bode P., 2008, ApJ, 682, 1
 Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
 Dehnen & Read (2011) Dehnen W., Read J. I., 2011, European Physical Journal Plus, 126, 55
 Dodelson & Schneider (2013) Dodelson S., Schneider M. D., 2013, Phys. Rev. D, 88, 063537
 Feng et al. (2016) Feng Y., Chu M.Y., Seljak U., 2016, preprint, (arXiv:1603.00476)
 Fosalba et al. (2008) Fosalba P., Gaztañaga E., Castander F. J., Manera M., 2008, MNRAS, 391, 435
 Fosalba et al. (2015a) Fosalba P., Gaztañaga E., Castander F. J., Crocce M., 2015a, MNRAS, 447, 1319
 Fosalba et al. (2015b) Fosalba P., Crocce M., Gaztañaga E., Castander F. J., 2015b, MNRAS, 448, 2987
 Giocoli et al. (2017) Giocoli C., Di Meo S., Meneghetti M., Jullo E., de la Torre S., Moscardini L., Baldi M., Mazzotta P., 2017, preprint, (arXiv:1701.02739)
 Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759
 Hartlap et al. (2011) Hartlap J., Hilbert S., Schneider P., Hildebrandt H., 2011, A&A, 528, A51
 Heitmann et al. (2010) Heitmann K., White M., Wagner C., Habib S., Higdon D., 2010, ApJ, 715, 104
 Heitmann et al. (2014) Heitmann K., Lawrence E., Kwan J., Habib S., Higdon D., 2014, ApJ, 780, 111
 Hoffmann et al. (2015) Hoffmann K., Bel J., Gaztañaga E., Crocce M., Fosalba P., Castander F. J., 2015, MNRAS, 447, 1724
 Howlett et al. (2015) Howlett C., Manera M., Percival W. J., 2015, Astronomy and Computing, 12, 109
 Hu (2000) Hu W., 2000, Phys. Rev. D, 62, 043007
 Izard et al. (2016) Izard A., Crocce M., Fosalba P., 2016, MNRAS, 459, 2327
 Jain et al. (2000) Jain B., Seljak U., White S., 2000, ApJ, 530, 547
 Kaiser (1992) Kaiser N., 1992, ApJ, 388, 272
 Kaiser (1998) Kaiser N., 1998, ApJ, 498, 26
 Kazin et al. (2014) Kazin E. A., Koda J., Blake C., Padmanabhan N., et al. 2014, MNRAS, 441, 3524
 Kiessling et al. (2011) Kiessling A., Heavens A. F., Taylor A. N., Joachimi B., 2011, MNRAS, 414, 2235
 Kilbinger (2015) Kilbinger M., 2015, Reports on Progress in Physics, 78, 086901
 Kim et al. (2011) Kim J., Park C., Rossi G., Lee S. M., Gott III J. R., 2011, Journal of Korean Astronomical Society, 44, 217
 Kitaura et al. (2015) Kitaura F.S., GilMarin H., Scoccola C. G., Chuang C.H., Muller V., Yepes G., Prada F., 2015, MNRAS, 450, 1836
 Koda et al. (2015) Koda J., Blake C., Beutler F., Kazin E., Marin F., 2015, preprint, (arXiv:1507.05329)
 Krause & Hirata (2010) Krause E., Hirata C. M., 2010, A&A, 523, A28
 Kuhlen et al. (2012) Kuhlen M., Vogelsberger M., Angulo R., 2012, Physics of the Dark Universe, 1, 50
 Leistedt et al. (2015) Leistedt B., et al., 2015, preprint, (arXiv:1507.05647)
 Li et al. (2011) Li B., King L. J., Zhao G.B., Zhao H., 2011, MNRAS, 415, 881
 Limber (1953) Limber D. N., 1953, ApJ, 117, 134
 Manera et al. (2013) Manera M., et al., 2013, MNRAS, 428, 1036
 Manera et al. (2015) Manera M., et al., 2015, MNRAS, 447, 437
 Merz et al. (2005) Merz H., Pen U.L., Trac H., 2005, New Astron., 10, 393
 Monaco (2016) Monaco P., 2016, preprint, (arXiv:1605.07752)
 Monaco et al. (2002) Monaco P., Theuns T., Taffoni G., 2002, MNRAS, 331, 587
 Monaco et al. (2013) Monaco P., Sefusatti E., Borgani S., Crocce M., Fosalba P., Sheth R. K., Theuns T., 2013, MNRAS, 433, 2389
 Munari et al. (2016) Munari E., Monaco P., Sefusatti E., Castorina E., Mohammad F. G., Anselmi S., Borgani S., 2016, preprint, (arXiv:1605.04788)
 Petri et al. (2016) Petri A., Haiman Z., May M., 2016, preprint, (arXiv:1612.00852)
 Quinn et al. (1997) Quinn T., Katz N., Stadel J., Lake G., 1997, preprint, (arXiv:astroph/9710043)
 Ross et al. (2015) Ross A. J., Samushia L., Howlett C., Percival W. J., Burden A., Manera M., 2015, MNRAS, 449, 835
 Scoccimarro & Sheth (2002) Scoccimarro R., Sheth R. K., 2002, MNRAS, 329, 629
 Sellentin & Heavens (2016) Sellentin E., Heavens A. F., 2016, MNRAS, 456, L132
 Sellentin & Heavens (2017) Sellentin E., Heavens A. F., 2017, MNRAS, 464, 4658
 Skillman et al. (2014) Skillman S. W., Warren M. S., Turk M. J., Wechsler R. H., Holz D. E., Sutter P. M., 2014, preprint, (arXiv:1407.2600)
 Springel (2005) Springel V., 2005, MNRAS, 364, 1105
 Takahashi et al. (2009) Takahashi R., et al., 2009, ApJ, 700, 479
 Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152
 Tassev et al. (2013) Tassev S., Zaldarriaga M., Eisenstein D. J., 2013, J. Cosmology Astropart. Phys., 6, 36
 Tassev et al. (2015) Tassev S., Eisenstein D. J., Wandelt B. D., Zaldarriaga M., 2015, preprint, (arXiv:1502.07751)
 Taylor et al. (2013) Taylor A., Joachimi B., Kitching T., 2013, MNRAS, 432, 1928
 Teyssier et al. (2009) Teyssier R., et al., 2009, A&A, 497, 335
 Valogiannis & Bean (2016) Valogiannis G., Bean R., 2016, preprint, (arXiv:1612.06469)
 Weinberg et al. (2013) Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87
 White & Hu (2000) White M., Hu W., 2000, ApJ, 537, 1
 White et al. (2010) White M., Pope A., Carlson J., Heitmann K., Habib S., Fasel P., Daniel D., Lukic Z., 2010, ApJ, 713, 383
 White et al. (2014) White M., Tinker J. L., McBride C. K., 2014, MNRAS, 437, 2594
 Winther et al. (2017) Winther H. A., Koyama K., Manera M., Wright B. S., Zhao G.B., 2017, preprint, (arXiv:1703.00879)
 Wright et al. (2017) Wright B. S., Winther H. A., Koyama K., 2017, preprint, (arXiv:1705.08165)
 Yu et al. (2016) Yu Y., Zhang P., Jing Y., 2016, Phys. Rev. D, 94, 083520
 de Putter & Takada (2010) de Putter R., Takada M., 2010, Phys. Rev. D, 82, 103522
 van Daalen & Schaye (2015) van Daalen M. P., Schaye J., 2015, preprint, (arXiv:1501.05950)
Appendix A Halo correlation function
In this appendix, we discuss an observable that affects the tangential shear or galaxygalaxy lensing: the halo clustering of a given lens population. In particular, the tangential shear contains information from both the weak lensing of the background sample and the bias of the foreground one. Therefore, for completeness, it is important to assess the accuracy on reproducing the clustering properties of the foreground sample. Fig. 10 displays the halo angular correlation function. Simulation measurements deviate with respect to the nonlinear prediction possibly due to both nonlinear bias and exclusion effects, that reduce the power at scales close to the halo size. In fact, for the sample used, with halos of mass or corresponding size of , at they subtend an angle of arcmin, in agreement with the scale at which we find possible exclusion effects in the measurements. The clustering amplitude at linear scales is lower in ICECOLA, which is consistent with the underestimation of the halo linear bias (in turn due to an underestimation of halo masses, see Izard et al. 2016).