© C. Gouin et al. 2019

1. Introduction

Gravitational lensing has become a versatile tool for probing the cosmological model and scenarios of galaxy evolution. From the coherent distortions, generated by the intervening matter along the line of sight, of the last scattering surface (e.g. Planck Collaboration VIII 2018) or intermediate-redshift galaxies (Bartelmann & Schneider 2001; Kilbinger 2015), to the inner parts of massive galaxies (Treu 2010), lensing directly measures the fractional energy density in matter of the Universe. Because it does not rely on assumptions about the relative distribution between the galaxies and the underlying dark matter (DM), which drives the dynamical evolution of cosmological structures, weak lensing plays a key role in recent, ongoing, or upcoming ground-based imaging surveys, such as the Canada France Hawaii Lensing Survey (Heymans et al. 2012), the Dark Energy Survey (Dark Energy Survey Collaboration 2005; Abbott et al. 2016), the Kilo-Degree Survey (KiDS: Kuijken et al. 2015), the Hyper Suprime-Cam Subaru Strategic Survey (Miyazaki et al. 2012), and the Large Synoptic Survey Telescope (LSST Science Collaborations 2009). It is also at the centre of the planned Euclid and WFIRST satellites (Laureijs et al. 2012; Spergel et al. 2015).

The statistical power of these experiments dramatically increases and drives on its way enormous efforts for the control of systematic effects. One of them concerns the accuracy to which theoretical predictions on the statistical properties of the matter distribution when it has evolved into the non-linear regime can be made on a small scale. Arguably, cosmological N-body numerical simulations have been playing a key role in solving the complex dynamical evolution of DM on scales smaller than a few Megaparsec (e.g. Springel et al. 2006). The upcoming Euclid or LSST missions require an extreme accuracy on the matter density power spectrum and the associated covariances that may enter a likelihood analysis of these data. The effort is currently culminating with the Flagship simulation, for instance (Potter et al. 2017). It also motivated earlier very large simulations such as Horizon-4π (Teyssier et al. 2009; Pichon et al. 2010), DEUS (Rasera et al. 2010), or MICE (Fosalba et al. 2015a), however. It has early been envisioned to propagate light rays through such DM simulations in order to reproduce the deflection and distortions of light bundles in a lumpy universe. The motivation is to derive lensing observables such as convergence maps and one-point probability distribution functions of this field or its topological properties (peaks, voids, etc.) or two-point shear correlation functions (Jain et al. 2000; Pichon et al. 2010; Hamana & Mellier 2001; Vale & White 2003; Hennawi & Spergel 2005; Hilbert et al. 2007, 2009; Sato et al. 2009). Much progress has since been made on large and mildly non-linear scales with the production of full-sky maps with an angular resolution of a few arcminutes (e.g. Fosalba et al. 2015b; Giocoli et al. 2016; Takahashi et al. 2017).

In order to make the most of the upcoming surveys, the matter distribution for Fourier modes as large as k ∼ 10 h Mpc⁻¹ must be predicted to the percent accuracy, which today still represents a challenge (Schneider et al. 2016). Furthermore, at these scales, the physics of baryons can differ from the dynamics of DM, and even though it amounts to ∼17% of the total cosmological matter budget, it has to be taken into account (van Daalen et al. 2011). For weak-lensing statistics, Semboloni et al. (2011) showed that the modelling of the two-point shear correlation function can be significantly biased when the baryons are simply treated like the collision-less DM. Even the number of convergence peaks itself is altered by baryons, but to a lesser extent than the power spectrum (Yang et al. 2013).

Recently, significant progress has been made on hydrodynamical simulations, which are now able to reproduce a morphological mix of galaxies in a cosmological context by considering baryonic physics such as radiative cooling, star formation, and feedback from supernovae and active galactic nuclei (AGN). Despite the balancing act that is required to be achieved between the high-resolution needs for properly describing the galaxies that formed at the centre of DM halos and the necessity of simulating sizeable cosmological volumes, recent simulations such as Horizon-AGN (Dubois et al. 2014), Illustris/Illustris-TNG (Vogelsberger et al. 2014; Pillepich et al. 2018), or EAGLE (Schaye et al. 2015) have now reached volumes of about 100 Mpc on a side and a resolution of about 1 kpc. This opens the possibility to quantify the effect of baryons (experiencing adiabatic pressure support, dissipative cooling, star formation, feedback, etc.) on the total matter distribution and its effect on lensing cosmological observables (see e.g. van Daalen et al. 2011; Tenneti et al. 2015; Hellwing et al. 2016; Springel et al. 2018; Chisari et al. 2018). Prescriptions to account for this effect (e.g. Semboloni et al. 2013; Schneider & Teyssier 2015; Mead et al. 2015; Rabold & Teyssier 2017) have been explored, and some start to be incorporated into cosmic shear studies (Hildebrandt et al. 2017).

In this paper, we further investigate the effect of baryons on lensing observables in the Horizon-AGN simulation. By taking advantage of the lightcone that is generated during the simulation run, we are able to fully account for projection effects (mixing physical scales) and small-scale non-linearities occurring in the propagation of light rays (e.g., Born approximation, lens-lens coupling, corrections for shear – reduced shear) that may be boosted by the steepening of the gravitational potential wells that are caused by cooled gas that sinks to the bottom of DM halos. This extends the analysis of Chisari et al. (2018), who mostly focused on the effect of baryons on the three-dimensional matter power spectrum and compared the Horizon-AGN results with those of Illustris, OWLS, EAGLE, and Illustris-TNG and found a broad qualitative agreement. The common picture is that hot baryons that are prevented from sinking into halos like DM induce a deficit of power inside halos (in a proportion of about Ω_b/Ω_M), and at still smaller scales (k ≳ 30 h Mpc⁻¹), baryons in the form of stars (and to a lesser extent, cooled gas) dramatically boost the amplitude of density fluctuations. However, even though these results seem to converge from one simulation to another, they substantially depend on the assumptions about sub-grid physics, and in particular, about AGN feedback.

In addition to these encouraging successes at quantifying the nuisance of baryons on cosmological studies, hydrodynamical simulations entail a wealth of information on the relation between galaxies or galaxy properties and the halo they live in. It is therefore a way to understand the large-scale biasing of these galaxies with respect to the overall total matter density field. We also explore the small-scale relation between galaxies and their surrounding gravitational potential that sources the lensing deflection field. In particular, the correlation between galaxies and the tangential distortion of background sources (so-called galaxy-galaxy lensing signal, GGL) has proven to be a way to constrain the galaxy-mass correlation function (e.g. Brainerd et al. 1996; Guzik & Seljak 2001; Mandelbaum et al. 2006, 2013; Leauthaud et al. 2012; Velander et al. 2014; Hudson et al. 2015; Coupon et al. 2015). In this vein, Velliscig et al. (2017) recently showed that the GGL around z ∼ 0.18 galaxies in the EAGLE simulation is consistent with the GGL measured around the Galaxies And Mass Assembly (GAMA) groups using KiDS data (Dvornik et al. 2017).

Finally, subtle observational effects entering GGL by high-redshift deflectors (z ≳ 0.8) are investigated from the lensing information over the full past lightcone of the Horizon-AGN simulation. The magnification bias affecting the selection of deflectors (Ziour & Hui 2008) complicates the interpretation of GGL substantially. Currently, no such high-z lens sample has been studied because even higher faintly lensed sources that carry the shear signal are scarce, but the situation may change with Euclid. Its slit-less grism spectroscopy will provide a large sample of Hα emitters in the 0.9 ≤ z ≤ 1.8 redshift range. A thorough understanding of the clustering properties of this sample may be achieved with the GGL measurement of this sample by using the high-z tail of the shape catalogue obtained with the visible imager (VIS) for Euclid. Some raytracing through cosmological simulations (Hilbert et al. 2009; Fosalba et al. 2015b) has briefly mentioned some aspects of the problem of magnification bias that was raised by Ziour & Hui (2008). The Horizon-AGN lightcone is a good opportunity to quantify these effects in order to correctly interpret upcoming GGLs. In this paper, cosmic shear or GGL quantities are directly measured from the lensing quantities obtained by ray-tracing methods. They are not inferred from the shape of galaxies, as is done in observations. A forthcoming paper will present the generation of mock wide-field images including lensing distortions from the full view of Horizon-AGN lightcone and the light emission predicted for the simulated stars, taking us one step closer to a full end-to-end generation of mock lensing observations.

The paper is organised as follows. Section 2 presents the Horizon-AGN hydrodynamical simulation, the structure of its lightcone, and some properties of the galaxy population therein. Section 3 describes the methods we implemented to generate the deflection field on thin lens planes and to propagate light rays through them. Section 4 describes the one- and two-point statistics of the resulting convergence and (reduced-)shear fields. The validity of the ray-tracing method is quantified by comparing our results with independent methods. Section 5 measures the GGL around the galaxies in the Horizon-AGN simulation. A comparison with observations is made for low-redshift deflectors. The problem of magnification bias is investigated for future observations of high-z GGL. Section 6 summarises our results.

2. Horizon-AGN simulation lightcone

2.1. Characteristics

The Horizon-AGN simulation is a cosmological hydrodynamical simulation performed with RAMSES (Teyssier 2002). The details of the simulations can be found in Dubois et al. (2014). We first briefly summarise the main characteristics. Horizon-AGN contains 1024³ DM particles with a mass resolution of 8 × 10⁷ h⁻¹ M_⊙ in a box of comoving size L_box = 100 h⁻¹ Mpc on a side. The gravity and hydrodynamics are treated in RAMSES with a multiscale approach with adaptive mesh refinement (AMR): starting from a uniform 1024³ grid, cells are then adaptively refined when the mass inside the cell exceeds eight times the initial mass resolution. Cells are recursively refined (or de-refined according to the refinement criterion) down to a minimum cell size of almost constant 1 proper kpc (an additional level is triggered at each expansion factor a = 0.1, 0.2, 0.4, 0.8). The underlying cosmology is a standard ΛCDM model consistent with the WMAP7 data (Komatsu et al. 2011), with total matter density Ω_m = 0.272, dark energy density Ω_Λ = 0.728, amplitude of the matter power spectrum σ₈ = 0.81, baryon density Ω_b = 0.045, Hubble constant of H₀ = 70.4 km s⁻¹ Mpc⁻¹, and scalar spectral index n_s = 0.967.

The evolution of the gas is solved on the RAMSES grid using a Godunov method with the approximate Harten-Lax-van Leer-Contact Riemann solver on the interpolated conservative hydrodynamical quantities, which are linearly interpolated at cell boundaries from their cell-centred values using a MinMod total variation diminishing scheme. In addition, accurate models of unresolved sub-grid physics have been implemented. The gas heating comes from a uniform UV background that started at the re-ionisation z_reion = 10 (Haardt & Madau 1996). The cooling function of the gas follows Sutherland & Dopita (1993), from H and He collision and from the contribution of other metals. Star formation is modelled following the Schmidt law (Kennicutt 1998), with a constant star formation efficiency of 2% per free fall time. It occurs when the density of the gas exceeds the threshold 0.1 H cm⁻³. The temperature at gas densities higher than 0.1 H cm⁻³ is modified by a polytropic equation of state with polytropic index of 4/3 and scaling temperature of 10⁴ K (Springel & Hernquist 2003). Stellar evolution is performed assuming a Salpeter (1955) initial stellar mass function. The sub-grid physics also includes stellar winds and supernova feedback in the form of heating, metal enrichment of the gas, and kinetic energy transfer to the ambient gas (see Kaviraj et al. 2017, for more details). Finally, black holes (BH) are created when the gas density exceeds 0.1 H cm⁻³, and when no other BH lies in the close environment. They grow by direct accretion of gas following an Eddington-limited Bondi-Hoyle-Littleton accretion rate, and by merger when BH binaries are sufficiently close. The AGN feedback is treated by either an isotropic injection of thermal energy, or by a jet as a bipolar outflow, depending on the ratio between the Bondi and the Eddington accretion rates (see Dubois et al. 2012; Volonteri et al. 2016, for details).

The past lightcone of the simulation was created on-the-fly as the simulation was running. Its geometry is sketched in Fig. 1. The opening angle of the cone is 2.25 deg out to redshift z = 1 and 1 deg all the way to z = 8. These two values correspond to the angular size of the full simulation box at these redshifts. We can therefore safely work in the flat sky (or infinitely remote observer) approximation. Up to z = 1, the volume of the cone is filled with about 7 replicates of the box. Between z = 0 and z = 4, the narrow cone contains about 14 replicates of the box, and the union of the two cones contains about 19 copies. This should be kept in mind when the statistical robustness of our results is quantified.

Fig. 1. 2D sketch of the past lightcone around redshift z = 1 (orange vertical line). Each mesh is a replicate of the Horizon-AGN simulation box (bounded with cyan lines). The tiling is performed all the way up to redshift z ∼ 8.

In order to limit projection effects, a non-canonical direction was chosen for the past lightcone, but in order to preserve periodic boundary conditions between replicates, no random rotation was applied. Projection effects will still be present and induce characteristic spectral distortions on large scales that must be taken into account. Particles and AMR cells were extracted on-the-fly at each coarse simulation time step (when all levels were synchronized in time because a factor 2 of subcycling is used between levels) of the simulation according to their proper distance to a fiducial observer located at the origin of the simulation box. The lightcone of the simulation thus consists of 22 000 portions of concentric shells. Each of them contains stellar BH DM particles (with their position and velocity, mass, and age), along with AMR Eulerian cells storing the gas properties (position, density, velocity, temperature, chemical composition, and cell size) and the total gravitational acceleration vector.

2.2. Properties of galaxies and host halos

The ADAPTAHOP halo finder (Aubert et al. 2004) was run on the lightcone to identify galaxies from the stellar particle distribution. The local stellar particle density was computed from the 20 nearest neighbours, and structures were selected with a density threshold equal to 178 times the average matter density at that redshift. Galaxies resulting in fewer than 50 particles (≃10⁸ M_⊙) were not included in the catalogue. Because the identification technique is redshift dependent, ADAPTAHOP was run iteratively on thin lightcone slices. Slices overlapped to avoid edge effects (i.e. cutting galaxies in the extraction) and duplicates were removed. In a second step, DM haloes were extracted independently from the DM particle distribution, with a density threshold of 80 times the average matter density, and keeping only haloes with more than 100 particles. The centre of the halo was temporarily defined as the densest particle in the halo, where the density was computed from the 20 nearest neighbours. In a subsequent step, a sphere of the size of the virial radius was drawn around it and a shrinking sphere method (Power et al. 2003) was implemented to recursively find the centre of mass of the halo. In each iteration, the radius of the halo was reduced by 10%. The search was stopped when a sphere three times larger than our spatial resolution was reached. Each galaxy was matched with its closest halo.

The simulation contains about 116 000 galaxies and halos in the simulation box at z = 0, with a limit of about M_* ≳ 2 × 10⁹ M_⊙. These yields have been extensively studied in previous papers of the Horizon-AGN series. For instance, Kaviraj et al. (2017) compared the statistical properties of the produced galaxies, showing a reasonable agreement with observed stellar mass functions all the way to z ∼ 6. The colour and star formation histories are also well recovered, and so are the BH – bulge relations and duty-cycles of AGNs (Volonteri et al. 2016).

Following up on an earlier work (Dubois et al. 2013) that focused on a handful of zoomed galaxy simulations with RAMSES, Dubois et al. (2016) confirmed with a much greater statistical significance in Horizon-AGN that the morphological diversity of galaxies is well reproduced (fraction of rotation- versus dispersion-supported objects, and how this dichotomy maps into the star-forming versus quiescent dichotomy). Taking advantage of a parallel simulation run with the same initial conditions and in which the AGN feedback is turned off (Horizon-noAGN), the key role of the latter in shaping the galaxy morphology was emphasised. Furthermore, Peirani et al. (2017) studied the effect of AGN feedback on the innermost density profiles (stars, gas, DM, and total) and found a good agreement of the density profile, size-mass relation, and DM fraction inside the effective radius of galaxies with observations. In particular, Peirani et al. (2019) showed that the innermost parts of Horizon-AGN galaxies are consistent with strong-lensing observations of Sonnenfeld et al. (2013) and Newman et al. (2013, 2015).

Populating the lightcone yields a volume-limited sample of 1.73 × 10⁶ galaxies in the narrow 1 deg cone. However, a large portion of the low-mass high-redshift galaxies would not be of much practical use in a flux-limited survey, as shown in Fig. 2, which plots the redshift-dependent limit in stellar mass that is attained with several i-band apparent limiting magnitudes. This was obtained using the COSMOS2015 photometric catalogue of Laigle et al. (2016).

Fig. 2. Distribution in the redshift – stellar mass plane of the 1.7 million galaxies in the Horizon-AGN lightcone. For guidance, the stellar mass limit for completeness is shown as well as fiducial cuts in mass that would be obtained with a flux-limited survey of various i-band limiting magnitudes.

3. Ray-tracing through the lightcone

After briefly describing the basics of the propagation of light rays in a clumpy universe and the numerical transcription of this formalism, we now describe the ray-tracing computation in the Horizon-AGN lightcone. Our implementation of the multiple lens plane (but also the Born approximation) builds on similar past efforts (Hilbert et al. 2008; Metcalf & Petkova 2014; Petkova et al. 2014; Barreira et al. 2016). It has been tailored for the post-treatment of the Horizon-AGN past lightcone, but provided the flat-sky approximation holds, our implementation could readily be applied to any other RAMSES lightcone output (Teyssier et al. 2009).

As detailed below, two methods are investigated to infer deflection angles from either the distribution of various particle-like matter components or the total gravitational acceleration stored by RAMSES. The light rays are then propagated plane by plane (both within and beyond the Born approximation) for these two different estimates of the deflection field.

3.1. Thin lens plane

We define β the (un-perturbed and unobservable) source plane angular position and θ the observed angular position of a light ray. Considering a unique thin lens plane, the relation between the angular position of the source β, the deflection angle α, and the image θ is

$$\begin{array}{c}\hfill \mathit{\beta }=\mathit{\theta }-\frac{D_{\mathrm{ls}}}{D_{\mathrm{s}}}\mathit{\alpha }(\mathit{\theta }),\end{array}$$(1)

where D_ls and D_s are the angular diameter distance between the source and the lens, and between the observer and the source, respectively. The deflection angle α(θ) is obtained by integrating the gravitational potential Φ(r) along the line of sight (here, radial proper coordinate x₃),

$$\begin{array}{c}\hfill \mathit{\alpha }(\mathit{\theta })=\frac{2}{c^2}{\displaystyle \int {\mathbf{\nabla }}_{\mathbf{\perp }}\mathrm{\Phi }(\mathit{\theta },x_3)}\hspace{0.17em}\mathrm{d}{x}_3.\end{array}$$(2)

Hence, across a thin lens plane, the lensing potential ϕ(θ) is related to the deflection field by the Poisson equation:

$$\begin{array}{c}\hfill \mathrm{\Delta }\varphi =\mathbf{\nabla }·\mathit{\alpha }\equiv 2\kappa ,\end{array}$$(3)

where the convergence κ is the projected surface mass density Σ(θ) in the lens plane expressed in units of the critical density Σ_crit

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit}}\kappa (\mathit{\theta })=\mathrm{\Sigma }(\mathit{\theta })\equiv {\displaystyle \int \rho (\mathit{\theta },z)}\mathrm{d}z.\end{array}$$(4)

The critical density reads

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit}}=\frac{c^2}{4\pi G}\frac{D_{\mathrm{s}}}{D_{\mathrm{l}}{D}_{\mathrm{ls}}},\end{array}$$(5)

with D_l, the angular diameter distance between the observer and the lens. In the above equations, all distances and transverse gradients are expressed in physical (proper) coordinates.

A Taylor expansion of the so-called lens Eq. (1) yields the Jacobian of the θ → β mapping, which defines the magnification tensor (e.g. Bartelmann & Schneider 2001)

$$\begin{array}{c}\hfill a_{\mathit{ij}}(\mathit{\theta })=\frac{\partial \mathit{\beta }}{\partial \mathit{\theta }}=({\delta }_{\mathit{ij}}-{\varphi }_{,ij})\equiv \left(\begin{array}{cc}1-\kappa -{\gamma }_1& -{\gamma }_2\\ -{\gamma }_2& 1-\kappa +{\gamma }_1\end{array}\right),\end{array}$$(6)

where δ_ij is the Kronecker symbol, and the two components γ_1/2 of the complex spin-2 shear have been introduced. We note that subscripts following a comma denote partial derivatives along that coordinate. Both shear and convergence are first derivatives of the deflection field α (or second derivatives of the lensing potential)

$$\begin{array}{cc}& \kappa =\frac{1}{2}({\alpha }_{1,1}+{\alpha }_{2,2}),\hfill \end{array}$$(7)

$$\begin{array}{cc}& {\gamma }_1=\frac{1}{2}({\alpha }_{1,1}-{\alpha }_{2,2}),\hfill \end{array}$$(8)

$$\begin{array}{cc}& {\gamma }_2={\alpha }_{1,2}={\alpha }_{2,1}.\hfill \end{array}$$(9)

Therefore, starting from pixelised maps of the deflection field α_1/2(i, j) in a thin slice of the lightcone, we can easily derive γ_1/2(i, j) and κ(i, j) with finite differences or fast Fourier transforms (FFTs), even if α is only known on a finite aperture, without periodic boundary conditions. Conversely, starting from a convergence map κ(i, j), it is impossible to integrate Eq. (3) with FFTs to obtain α (and then differentiate again to obtain γ) without introducing edge effects if the periodic boundary conditions are not satisfied. Additionally, we also introduce the scalar magnification μ, which is the inverse determinant of the magnification tensor a_i, j of Eq. (6).

3.2. Propagation of rays in a continuous lumpy universe

On cosmological scales, light rays cross many over- or under-dense extended regions at different locations. Therefore, the thin lens approximation does not hold. The transverse deflection induced by an infinitely thin lens plane is still given by the above equations, but the trajectory of rays along their path needs to be fully integrated. For a given source plane at comoving distance χ_s, the source plane position of a ray, initially observed at position θ, is therefore given by the continuous implicit (Voltera) integral equation (Jain & Seljak 1997)

$$\begin{array}{c}\hfill \mathit{\beta }(\mathit{\theta },{\chi }_{\mathrm{s}})=\mathit{\theta }-\frac{2}{c^2}{\displaystyle {\int }_0^{{\chi }_{\mathrm{s}}}\mathrm{d}\chi \frac{{\chi }_{\mathrm{s}}-\chi }{{\chi }_{\mathrm{s}}\hspace{0.17em}\chi }{\mathbf{\nabla }}_{\beta }\varphi (\mathit{\beta }(\mathit{\theta },\chi ),\chi )}.\end{array}$$(10)

To first order, the gravitational potential along an unperturbed path can be evaluated, so that

$$\begin{array}{c}\hfill \mathit{\beta }(\mathit{\theta },{\chi }_{\mathrm{s}})=\mathit{\theta }-\frac{2}{c^2}{\displaystyle {\int }_0^{{\chi }_{\mathrm{s}}}\mathrm{d}\chi \frac{{\chi }_{\mathrm{s}}-\chi }{{\chi }_{\mathrm{s}}\hspace{0.17em}\chi }{\mathbf{\nabla }}_{\theta }\varphi (\mathit{\theta },\chi )}.\end{array}$$(11)

This is known as the Born approximation, which is common in many diffusion problems of physics. An interesting property of the Born approximation is that the relation between β and α can be reduced to an effective thin lens identical to Eq. (1) allowing the definition of an effective convergence, which is the divergence of the effective (curl-free) deflection field: 2κ_eff = ∇ ⋅ α_eff.

When the approximation does not hold, the relation between β and α can no longer be reduced to an effective potential and some curl-component may be generated, implying that the magnification tensor is no longer symmetric but requires the addition of a rotation term ω and so-called B-modes in the shear field. In this more general framework, the magnification tensor should be rewritten

$$\begin{array}{c}\hfill a_{\mathit{ij}}(\mathit{\theta })=\left(\begin{array}{cc}1-\kappa -{\gamma }_1& -{\gamma }_2-\omega \\ -{\gamma }_2+\omega & 1-\kappa +{\gamma }_1\end{array}\right).\end{array}$$(12)

3.3. Multiple lens planes approximation

The numerical transcription of Eq. (10) in the Horizon-AGN past lightcone requires the slicing of the latter into a series of parallel transverse planes, which could simply be the 22 000 slabs dumped by RAMSES at runtime every coarse time step. These are too numerous and can safely be stacked into thicker planes by packing together 40 consecutive slabs¹. Here 500 slices of varying co-moving thickness were produced all the way to redshift z = 7 to compute either the deflection field or the projected surface mass density as described below.

The discrete version of the equation of ray propagation (10) for a fiducial source plane corresponding to the distance of the plane j + 1 reads

$$\begin{array}{c}\hfill {\mathit{\beta }}^{j+1}=\mathit{\theta }-{\displaystyle \sum _{i=1}^j\frac{D_{i;j+1}}{D_{j+1}}{\mathit{\alpha }}^i({\mathit{\beta }}^i)},\end{array}$$(13)

where αⁱ is the deflection field in the lens plane i, D_j + 1 is the angular diameter distance between the observer and the plane j + 1, and D_{i; j + 1} is the angular diameter distance between planes i and j + 1. Therefore, as sketched in Fig. 3, rays are recursively deflected one plane after the other, starting from unperturbed positions on a regular grid θ ≡ β¹.

Fig. 3. Schematic view of the propagation of a light ray through a lightcone sliced into multiple discrete lens planes. The ray (red line) is deflected at each intersection with a thin lens plane. The deflection field is defined for each plane depending of the angular position on this plane αj(βj).

The practical implementation of the recursion in Eq. (13) is computationally cumbersome and demanding in terms of memory because the computation of the source plane positions β^j + 1 requires holding all the j previously computed source plane positions. Instead, this paper follows the approach of Hilbert et al. (2009), who showed that Eq. (13) can be rewritten as a recursion over only three consecutive planes²

$$\begin{array}{c}\hfill {\mathit{\beta }}^{j+1}=(1-\frac{D_j}{D_{j+1}}\frac{D_{j-1;j+1}}{D_{j-1;j}}){\mathit{\beta }}^{j-1}+\frac{D_j}{D_{j+1}}\frac{D_{j-1;j+1}}{D_{j-1;j}}{\mathit{\beta }}^j-\frac{D_{j;j+1}}{D_j}{\mathit{\alpha }}^j({\mathit{\beta }}^j).\end{array}$$(14)

In addition to this thorough propagation of light rays, the source plane positions and associated quantities (convergence κ, shear γ, and rotation ω) were additionally computed using the Born approximation, following the discrete version of Eq. (11):

$$\begin{array}{c}\hfill {\mathit{\beta }}^{j+1}=\mathit{\theta }-{\displaystyle \sum _{i=1}^j\frac{D_{i;j+1}}{D_{j+1}}{\mathit{\alpha }}^i(\mathit{\theta })}.\end{array}$$(15)

The deflection maps in each lens plane were computed on a very fine grid of pixels of constant angular size. In order to preserve the spatial resolution of about 1 kpc that is allowed by the simulation at high redshift, 36 000 × 36 000 deflection maps were built in the narrow 1 deg lightcone. The deflection maps in the low-redshift 2.25 sq deg wide cone that reach z = 1 were computed on a coarser 20 000 × 20 000 pixel grid because the actual physical resolution of the simulation at low redshift does justify the 0.1 arcsec resolution of the narrow 1 sq deg field of view. Even though the image plane positions θ = β¹ are placed on the regular pixel grid, the deflections they experience must be interpolated in between the nodes of the regular deflection map as they progress backward to a given source plane. This was done with a simple bilinear interpolation scheme.

3.4. Total deflections from the RAMSES accelerations

We now describe how we obtained α that we used in Eqs. (14) and (15). The first method uses the gravitational acceleration field, which is registered on each (possibly refined) grid location inside the lightcone. The very same gravitational field that was used to move particles and evolve Eulerian quantities in RAMSES was interpolated at every cell position and was therefore used to consistently derive the deflection field. The merits of the complex three-dimensional multi-resolution Poisson solver are therefore preserved and the transverse components of the acceleration fields can readily be used to infer the deflection field. By integrating the transverse component of the acceleration along the light of sight, we can compute the deflection field according to Eq. (2).

To do so, gas cells that intersect the ray were considered for each light ray, and the intersection length along the line-of-sight l_i was computed. Knowing the cell size δ_i, and its orientation with respect to the line of sight, we deduced l_i with an oriented-box-boundary (OBB) algorithm (e.g. Akenine-Möller et al. 2008) that assumes that all cells share the same orientation (flat-sky approximation), and we factorised out expensive dot products between normals to cell edges and the line of sight,

$$\begin{array}{c}\hfill \mathit{\alpha }(\mathit{\theta })=\frac{2}{c^2}{\displaystyle \sum _{i\in \mathcal{V}(\mathit{\theta })}{\mathbf{\nabla }}_{\mathbf{\perp }}{\varphi }_i(\mathit{\theta })\hspace{0.17em}{l}_i},\end{array}$$(16)

where 𝒱(θ) denotes the projected vicinity of a sky position θ. As shown in Fig. 4, a fiducial light ray is drawn: at each lens plane, the deviation of the light is calculated as the direct sum of the transverse acceleration components recorded on the cells i, weighted by the intersection length l_i. Here, the field of view is small and we can safely assume that light rays share the same orientation (flat-sky approximation) and are parallel to the line of sight.

Fig. 4. Schematic view of the problem induced by cells at the boundary of slabs j and j + 1, which become refined between time t and t + dt. Missing cells (devoid of dots) or cells in excess (overlapping dotted cells of different colour) can end up as lightcone particles. A fiducial light ray is drawn to illustrate the intersection length li between the ray and RAMSES cells.

This method has the main advantage that it preserves the gravitational force that was used when the simulation was evolved. In particular, the way shot noise is smoothed out in the simulation to recover the acceleration field from a mixture of Lagrangian particles and Eulerian gas cells is faithfully respected in the ray-tracing. In other word, the force felt by photons is very similar to the force felt by particles in the simulation. Acceleration is also local, in the sense that the deflection experienced by a light ray (and related derivatives leading to e.g. shear and convergence) depends only on the acceleration of cells that this ray crosses. The mass distribution outside the lightcone is therefore consistently taken into account through the acceleration field.

However, this method is sensitive to small artefacts that are present at the lightcone generation stage (i.e. simulation runtime) and that could not be corrected without a prohibitive post-processing of the lightcone outputs. When the simulation dumps two given neighbouring slabs at two consecutive time steps, problems can arise if cells on the boundary between the two slabs have been (de-)refined in the mean time. As illustrated in Fig. 4, such cells can be counted twice or can be missing if they are refined (or de-refined) at the next time step. These bumps and dips in the deflection map translate into saw-tooth patterns in the convergence maps. They are quite scarce and of very modest amplitude, however.

A 100 arcsec wide zoom into the convergence map obtained with this method is shown in the left panel of Fig. 5. The source redshift is z_s = 0.8. A few subdominant artefacts due to missing acceleration cells are spotted. They induce small correlations on scales smaller than a few arcseconds and are otherwise completely negligible for our cosmological applications.

Fig. 5. Comparison of zs = 0.8 convergence maps obtained with the OBB method (integration of transverse accelerations in cells, left) and with the SPL method (projection of particles onto convergence planes after adaptative Gaussian smoothing, right). The latter method applies a more aggressive smoothing that better erases shot noise. Inaccuracies of long-range deflections in the SPL method that are due to edge effects translate into a global shift for some galaxies, as compared to OBB. With this method, some missing acceleration cells occasionally produce modest artefacts on a small scale.

3.5. Projection of smoothed particle density

The second method of computing the deflection maps in thin lens planes is more classical: it relies on the projection of particles onto surface density maps that are then turned into deflection maps. If the line-of-sight integration is performed under the Born approximation, the Fourier inversion going from the projected density to the deflection is just done once starting from the effective convergence. Otherwise, with the full propagation, many FFT inversions on projected density maps that do not fulfil the periodic boundary condition criterion imply an accumulation of the inaccuracies in the Fourier inversion.

First of all, this method allows us to separate the contribution of each matter component from the total deflection field. We can therefore compute the contribution of stars or gas to the overall lensing near a given deflector. This is not possible with the acceleration method because only the total acceleration is computed by the simulation.

In addition, particles can be projected with an efficient and adaptive smoothing scheme. Instead of a standard nearest grid point or cloud-in-cell projection, a Gaussian filter (truncated at four times the standard deviation σ) is used in which the width of the smoothing filter σ is tuned to the local density, hence following the smooth particle lensing (SPL) method of Aubert et al. (2007). Because the AMR grid of RAMSES is adaptive, the resolution level around a given particle position from the neighbouring gas cells can be recovered. This thus bypasses the time-consuming step of building a tree in the distribution of particles, which is at the heart of the SPL method.

To illustrate the merits of this method and for comparison with the previous one, we show the same region of simulated convergence fields for a source redshift z_s = 0.8 in the right panel of Fig. 5. This adaptive Gaussian smoothing (referred to as SPL method below) seems more efficient at smoothing the particle noise out. Between the two methods, we note small displacements of some galaxies of a few arcseconds. They are due to the long-range inaccuracies generated by the Fourier inversions.

3.6. Lensing of galaxy and halo catalogues

In order to correlate galaxies (or halos) in the lightcone with the convergence or shear field around them and, hence, measure their GGL, their catalogue positions β (which are intrinsic source plane coordinates) need to be shifted and their observed lensed image plane positions θ need to be inferred. These are related by the thorough lens Eq. (10), or its numerical translation (13). However, this equation is explicit only for the θ → β mapping. The inverse relation, which can be multi-valued when strong lensing occurs, has to be solved numerically by testing for every image plane mesh θ_ij whether it surrounds the coordinates β^gal of the deflected galaxy when cast into the source plane β_ij (e.g. Schneider et al. 1992; Keeton 2001; Bartelmann 2003). Because the method should work in the strong-lensing regime, regular rectangular meshes may no longer remain convex in the source plane and it is therefore preferable to split each mesh into two triangles. These triangles will map into triangles in the source plane and we can safely test whether β^gal is inside them. In order to speed up the test on our large pixel grids, the image plane was partitioned into a quad-tree structure that recursively explores finer and finer meshes. The method is very fast and yields all the image plane antecedents of a given galaxy position β^gal. This provides us the updated catalogues of halos and galaxies. Obviously, when the GGL signal is measured in the Born approximation, catalogue entries do not need to be deflected and therefore source plane and image plane coordinates are identical.

3.7. Summary of the generated deflection maps

Table 1 summarises the main advantages and drawbacks of the OBB and SPL methods. Altogether, 2 × 2 (OBB/SPL and Born approximation/full propagation) deflection maps were generated for each of the 246 source planes all the way to z = 1 in the wide opening angle field. Likewise, we obtained 2 × 2 maps for each of the 500 source planes all the way to z = 7 in the narrow opening angle field.

4. Cosmic shear

This section assesses the validity of our ray-tracing methods by measuring one- and two-point statistics of the lensing quantities such as convergence and (reduced-)shear. It also compares these finding with other methods.

The focus is on the effect of baryons on small scales for multipoles ℓ ≳ 2000 to check whether the baryonic component is connected to other non-linear effects like the shear – reduced shear correction and beyond-Born treatments.

4.1. Convergence of one-point statistics

The most basic quantity that can be derived from the convergence field shown in the right panel of Fig. 6 is the probability distribution function (PDF) of the convergence. The Fig. 6 shows this quantity, which is extremely non-Gaussian at the ∼1″ resolution of the map. The skewness of the field is visible, with a prominent high-end tail and a sharp decrease in negative convergence values.

Fig. 6. Left panel: convergence map generated with a $0\stackrel{″}{.}1$ pixel grid over a 2.25 × 2.25 sq. deg. field of view for a fiducial source plane at zs ∼ 1. Right panel: convergence map with a field of view of 1 sq. deg. at zs ∼ 2, and its corresponding convergence PDF showing the characteristic skewed distribution.

4.2. Convergence power spectrum

In Fourier space, the statistical properties of the convergence field are commonly characterised by its angular power spectrum P_κ(l),

$$\begin{array}{c}\hfill ⟨\widehat{\kappa }(\mathit{\ell }){\widehat{\kappa }}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }})⟩={(2\pi )}^2{\delta }_D(\mathit{\ell }-{\mathit{\ell }}^{\mathbf{\prime }})P_{\kappa }(\ell ),\end{array}$$(17)

where δ_D(ℓ is the Dirac delta function. For two fiducial source redshifts (z_s = 0.5 and z_s = 1), Fig. 7 shows the angular power spectrum of the convergence obtained with the two ray-tracing techniques: the OBB and SPL methods (solid magenta and solid cyan curves, respectively). The low-redshift methods are based on the 2.25 deg wide lightcone. They are thus more accurate on larger scales ℓ ≲ 10³, even though the large sample variance will not permit quantitative statements. On small scales (ℓ ∼ 2 × 10⁵), the additional amount of smoothing implied by the SPL projection of particles onto the lens planes induces a deficit of power with respect to the less aggressive softening of the OBB method in which shot noise has not been entirely suppressed (see Fig. 5).

Fig. 7. Upper panel: convergence power spectra for source redshift zs = 1 (top) and zs = 0.5 (bottom) derived with the OBB (magenta) and the SPL (green) methods. The more aggressive smoothing of this latter method translates into a faster high-ℓ fall-off. The cyan curves (DM) only account for the DM component (rescaled by 1 + Ωb/ΩM). The red curve corresponds to the direct integration of the three-dimensional total matter power spectrum (Limber approximation) in the Horizon-AGN simulation (Hz-AGN). The blue curves are the direct integration of the Horizon-DM (DM-only) matter power spectrum (Hz-DM). Dashes reflect regimes where the three-dimensional spectra of Chisari et al. (2018) were extrapolated by a simple power law (extrapolation). The yellow lines show the particle shot-noise contribution at two different redshifts. Middle panel: ratio of the zs = 0.5 convergence power spectra obtained with the Born approximation and the proper multiple lens plane integration showing only very small changes up to ℓ ∼ 105. Bottom panel: ratio of the DM-only to total convergence power spectra at zs = 0.5,  1.0, and 1.5 for the SPL method.

The middle panel of Fig. 7 shows the difference between power spectra inferred using the Born approximation or with the full multiple lens plane approach for the OBB method. For angular scales ℓ ≲ 8 × 10⁴, we find differences between the two propagation methods that are smaller than about 0.5%, which is totally negligible given possible numerical errors and sampling variance limitations. At lower angular scales ℓ ≳ 10⁵, departures rise above the few percent level. We note that this scale also corresponds to the scale where shot noise (from DM particles) and convergence power spectral are of equal amplitude (yellow shaded area). Below these very small scales, close to the strong-lens regime, the Born approximation may start to break down (Schäfer et al. 2012).

Under the Limber and Born approximations, the convergence power spectrum can be expressed as an integral of the three-dimensional non-linear matter power spectrum P_δ (Limber 1953; Blandford et al. 1991; Miralda-Escudé 1991; Kaiser 1992) from the observer to the source plane redshift or corresponding co-moving distance χ_s:

$$\begin{array}{c}\hfill P_{\kappa }(\ell )={\left(\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2}\right)}^2{\displaystyle {\int }_0^{{\chi }_{\mathrm{s}}}\mathrm{d}\chi }{\left(\frac{\chi ({\chi }_{\mathrm{s}}-\chi )}{{\chi }_{\mathrm{s}}a(\chi )}\right)}^2{P}_{\delta }(\frac{\ell }{\chi },\chi ),\end{array}$$(18)

where a is the scale factor and where no spatial curvature of the universe was assumed for conciseness and because the cosmological model in Horizon-AGN is flat. As a validation test of our light-deflection recipes, the lensing power spectrum derived from the actual ray-tracing was compared to an integration of the three-dimensional matter power spectrum measured by Chisari et al. (2018) in the Horizon-AGN simulation box. The red curve is the direct integration of P_δ(k) power spectra, and the dashed parts of the lines corresponds to a power-law extrapolation of the P_δ(k) down to smaller scales. In the range 3000 ≲ ℓ ≲ 3 × 10⁵, an excellent agreement is found between the red curve and the spectra inferred with our two ray-tracing techniques. On larger scales, the cosmic variance (which is different in the full simulation box and the intercept of the box with the lightcone) prevents any further agreement. This is also the case for ℓ ≳ 3 × 10⁵, where some possibly left-over shot noise in the ray-tracing maps and the hazardous high-ℓ extrapolation of the three-dimensional power spectra complicate the comparison. In addition, the low-ℓ oscillations of the spectrum is likely to originate from the replicates of the simulation box throughout the past lightcone.

Chisari et al. (2018) also measured matter power spectra in the Horizon-DM simulation at various redshifts. This simulation is identical to Horizon-AGN in terms of initial conditions, but has been run without any baryonic physics in it after the mass of DM particles was rescalded to conserve the same total matter density (Peirani et al. 2017; Chisari et al. 2018). The integration of this DM-only power spectrum allows us to obtain a sense on the effect of baryons in the DM-distribution itself. In the same way as the red curve showed the result of the Limber integral in Eq. (18) for Horizon-AGN, the dark blue curve shows the same integral for Horizon-DM. The latter has much less power for ℓ ≳ 2 × 10⁴ than either the integration of the full physics Horizon-AGN matter power spectrum (red) or that derived directly from ray-tracing (purple or green). The boost of spectral amplitude is due to cool baryons in the form of stars at the centre of halos. Moreover, we note a deficit of power on scales 2 × 10³ ≲ ℓ ≲ 2 × 10⁴ for the full physics simulation. As pointed out by Semboloni et al. (2011), the pressure acting on baryons prevents them from falling onto halos as efficiently as DM particles, hence reducing the depth of the potential wells when compared to a DM-only run. This effect has previously been investigated with more sensitivity on the three-dimensional matter power spectrum in the Horizon-AGN simulation (Chisari et al. 2018), and a clear dip in the matter density power spectrum of the full physics simulation is observed on scales 1 ≲ k ≲ 10 h Mpc⁻¹. Here, the projection somewhat smears out this dip over a larger range of scales, but a ∼15% decrease in amplitude is typically observed for ℓ = 10⁴ at z_s = 0.5. In order to show the changes due to the inclusion of the baryonic component more clearly, we traced rays through the lightcone by considering only the DM particles of the Horizon-AGN run with the SPL method. For this particular integration of rays trajectories, we multiplied the mass of the DM particules by a factor 1 + Ω_b/Ω_DM (where Ω_DM = Ω_m − Ω_b) to obtain the same overall cosmic mean matter density. The cyan curve in the upper panel shows the resulting convergence power spectrum. The ratios between the total full physics convergence power spectrum and the rescaled DM contribution of this power spectrum at z_s = 0.5, 1.0 and 1.5 are shown in the bottom panel and further illustrate the two different effects of baryons on intermediate and small scales.

By considering two raytracing methods to derive the convergence power spectrum, and by asserting that consistent results are obtained by integrating the three-dimensional matter power spectrum, we now search for small scale effects that involve the possible coupling between the baryonic component and shear – reduced shear corrections.

4.3. Shear – reduced shear corrections to two-point functions

In practical situations, rather than the convergence power spectrum, which is not directly observable, wide field surveys give access to the angular correlation of pairs of galaxy ellipticities. The complex ellipticity³ε is directly related to the shear γ. The relation between the ensemble mean ellipticity and the shear is

$$\begin{array}{c}\hfill ⟨\epsilon ⟩=g\equiv \frac{\gamma }{1-\kappa }\simeq \gamma ,\end{array}$$(19)

with g the so-called reduced shear. Therefore, the two point correlations of ellipticities and shear only match when the convergence κ is small. Because the regions of large convergence are typically the centres of halos where the contribution of cooled baryons is highest, we might expect a coupling between the inclusion of baryons and the shear – reduced-shear corrections needed to properly interpret the cosmological signal carried by the two-point statistics (e.g. White 2005; Kilbinger 2010)

Owing to the spin-2 nature of ellipticity, we can define the angular correlation functions ξ_±

$$\begin{array}{cc}\hfill {\xi }_{\pm }(\theta )& ={⟨{\gamma }_{+}(\vartheta +\theta ){\gamma }_{+}(\vartheta )⟩}_{\vartheta }\pm {⟨{\gamma }_{\times }(\vartheta +\theta ){\gamma }_{\times }(\vartheta )⟩}_{\vartheta },\hfill \\ \hfill & =2\pi {\displaystyle \int \mathrm{d}\ell \ell J_{0/4}(\theta \ell )P_{\kappa }(\ell )},\hfill \end{array}$$(20)

where γ₊ and γ_× are defined with respect to the separation vector between two galaxies or, here, any two image plane positions at separation θ. J₀ and J₄ are zeroth- and fourth-order Bessel functions.

Instead of the shear, observers can only measure associated ellipticities ϵ, which should thus replace γ in Eq. (20) in practical measurements. The reduced shear maps were computed together with shear and convergence maps, so as to measure the modified ξ₊ and ξ₋ angular correlations to compare them with the actual correlation functions. For efficiency, the Athena code⁴ was used to compute correlation functions.

The results are shown in Fig. 8 for a fiducial source redshift z_s = 0.5. Here ${\xi }_{+}^g$ and ${\xi }_{+}^{\gamma }$ only depart from one another at the ∼2−3% level on angular separations ∼1′. The effect is slightly stronger for ξ₋, which is known to be more sensitive to smaller non-linear scales than ξ₊, but is also more difficult to measure in the data because of its lower amplitude. On 1′ scales, ${\xi }_{-}^g/{\xi }_{-}^{\gamma }-1\simeq 7-8\%$. Like for the power spectra in the previous subsection, the cyan curves represent the correlations ${\xi }_{\pm }^{\gamma }$ for the rescaled DM contribution. The bottom panel shows the ratio of rescaled DM over full physics reduced shear correlation functions, further illustrating the effect of baryons on small scales. Again, ξ₋ responds more substantially to the inclusion of baryons. The deficit of correlation amplitude when baryons are taken into account peaks at 3−4′ and is of about 10%. Below 1′, the effect starts to increase, but those scales are never used in practical cosmic shear applications. We show in the next section that these scales remain perfectly relevant for galaxy evolution studies by means of the galaxy-galaxy weak-lensing signal.

Fig. 8. Upper panel: two-point shear correlation functions ξ+ (solid lines) and ξ− (dotted lines) for a fiducial source redshift zs = 0.5. We either correlate actual shear (red) or reduced shear (green) in the calculation to highlight the small-scale effect of baryons on this non-linear correction. Middle panel: ratio of shear correlation functions for the two cases. Bottom panel: ratio of shear correlation functions for a ray-tracing that only includes rescaled DM particles or all the components.

5. Galaxy-galaxy lensing

Focussing further on DM halos, we now investigate the yields of the simulation in terms of the galaxy-galaxy weak-lensing signal. The tangential alignment of background galaxies around foreground deflectors is substantially altered by the aforementioned baryonic physics, and we also expect a strong signature in this particular lensing regime.

For a circularly symmetric mass distribution Σ(R), we can relate shear, convergence, and the mean convergence enclosed inside a radius R centred on a foreground galaxy or halo as

$$\begin{array}{c}\hfill \overline{\kappa }(<R)=\frac{2}{R^2}{\displaystyle {\int }_0^R\kappa (R\prime )R\prime \mathrm{d}R\prime =\kappa (R)+\gamma (R)}.\end{array}$$(21)

Using the definition of the critical density given in Eq. (5), we can define the excess density

$$\begin{array}{cc}\hfill \mathrm{\Delta }\mathrm{\Sigma }(R)& =\frac{M(<R)}{\pi R^2}-\mathrm{\Sigma }(R),\hfill \\ \hfill & ={\mathrm{\Sigma }}_{\mathrm{crit}}\gamma (R).\hfill \end{array}$$(22)

The previous section has shown that the lensing convergence or shear maps have adequate statistical properties, and in Sect. 3.6 we showed how to use the associated deflection maps to map our lightcone galaxy catalogue into the image plane. In addition, galaxies are also expected to become magnified when lensed. Future extensions of this work will include the realistic photometry of the Horizon-AGN galaxies. We can easily account for the magnification bias by multiplying stellar masses by the magnification μ, however, as if luminosity or flux were a direct proxy for stellar mass. In the following, we refer to M_* for the intrinsic and μM_* for the magnified mass proxy.

For any given source redshift, we average the tangential shear around galaxies of any given stellar mass M_* (or more realistically magnified stellar mass μM_*), in order to estimate the GGL around Horizon-AGN galaxies. This is done around deflected galaxy positions.

5.1. Comparison with CMASS galaxies

We first compared the GGL around Horizon-AGN galaxies with the GGL excess mass profiles obtained by Leauthaud et al. (2017), who analysed the spectroscopic CMASS sample of massive galaxies in the footprint of the CFHTLS and CS82 imaging surveys, which covered ∼250 deg². These authors paid particular attention to quantifying the stellar mass of the CMASS galaxies that are centred around lens redshift z ∼ 0.55. The CMASS sample is not a simple mass selection, and includes a set of colour cuts, which makes this just a broad comparison. These results are somewhat sensitive to the detailed distribution in stellar mass above that threshold. The sample mean mass only slightly changes with redshift, but remains close to 3 × 10¹¹ M_⊙.

In order to match this lens sample, we extracted from the wide low-redshift lightcone the galaxies in the redshift range 0.4 ≤ z ≤ 0.70, and with a stellar mass above a threshold that was chosen to match the CMASS mean stellar mass. Even though these galaxies centred around lens redshift z ∼ 0.52 were treated as lens galaxies, they experience a modest amount of magnification (they behave like sources behind the mass distribution at yet lower redshift, see Sect. 5.2). We therefore chose galaxies that satisfy μM_* >  1.7 × 10¹¹ M_⊙. At this stage, selecting for M_* or μM_* does not make any significant difference (≲4%) because of the relatively low redshift of the lens sample. By doing so, we obtained the same sample mean stellar mass as the CMASS sample.

We then measured the mean tangential shear around these galaxies for a fiducial unimportant source redshift z_s = 1 and converted shear into excess density ΔΣ. The result is shown in Fig. 9. A good agreement between our predictions (OBB method, green with lighter envelope) and the observations of Leauthaud et al. (2017) (blue dots) is found, further suggesting that Horizon-AGN galaxies live in the correct massive halos (M_h ≃ 10¹³ M_⊙), or at the very least, produce the same shear profile as CMASS galaxies around them. We note that we split the 2.25 deg field of view into four quadrants and used the dispersion in these areas to compute a rough estimate of the model uncertainties.

Fig. 9. Comparison of the GGL tangential shear signal around z = 0.55 Horizon-AGN galaxies (green curve surrounded by a light green ribbon) and the GGL observations of Leauthaud et al. (2017) (blue dots with error bars). Units are all physical (and not comoving!). Model uncertainties in the simulation past lightcone are roughly estimated by splitting the 2.25 deg wide field of view into four quadrants. They may be underestimated beyond 1 h−1 Mpc. Cuts in stellar mass are expressed in units of 1011 M⊙. Black, blue, and red curves show the GGL shear signal predicted with the SPL method for the total, DM, and baryonic mass distributions, respectively. For clarity, uncertainties are omitted. They are similar to those in the case of the OBB method (green).

On scales R ≲ 0.2 h⁻¹ Mpc, the shear profile is 10−15% above the observations. Answering whether the discrepancy is due to faulty subgrid baryonic physics, a missing cosmological ingredient (or not perfectly adequate cosmological parameters), or leftover systematics in the data will certainly require more GGL observations, possibly combined with yet smaller scale strong-lensing and kinematical data (e.g. Sonnenfeld et al. 2018). Small-scale GGL clearly is a unique tool for addressing these issues (e.g. Velliscig et al. 2017), and asserting that the galaxy-halo connection is correctly reproduced by the simulations all the way to z ≳ 1 is arguably one of the foremost goals of galaxy formation models.

Figure 9 also shows our GGL results for the same population of lenses at the same redshift, but as inferred from the SPL method (solid black), which allows us to split the total lensing signal into its DM (blue and baryonic components (red). First of all, the agreement between the two methods for the total lensing signal is remarkable, except on scales ≳2 Mpc ∼ 5′, where differences begin to exceed the percent level. As we described in the previous section, this is due to inaccuracies of the Fourier transforms performed with the SPL method. We can use this latter technic to compare the contribution of DM and baryons (stars+gas), however. Clearly, the total and DM profiles look very similar beyond ∼0.2 Mpc up to a ∼17% renormalisation of the matter density. Only below these scales begin cooled baryons (stars) to contribute substantially. We predict an equal contribution of DM and stars to the total shear signal near a radius ∼15 kpc. We refer to Peirani et al. (2017) for further details about the innermost density profiles around Horizon-AGN galaxies in the context of the cusp-core problem.

5.2. High-redshift magnification bias

For z_l ≳ 0.6, the lens population begins to be lensed by yet nearer structures. This can lead to a magnification bias, which was studied by Ziour & Hui (2008).

The spatial density of a lensed population of background sources can be enhanced or decreased by magnification as light rays travel through over- or under-dense sight-lines (e.g. Moessner & Jain 1998; Moessner et al. 1998; Ménard & Bartelmann 2002; Scranton et al. 2005). Furthermore, the fraction of sources that are positively or negatively magnified depends on the slope of the luminosity function of the population. If it is very steep (typically the bright end of a population), we can observe a dramatic increase in the number of bright lensed objects. These deflectors appear brighter than they actually are Fig. 10 shows the mean magnification experienced by Horizon-AGN lightcone galaxies above a given stellar mass threshold (mimicking a more realistic flux limit) as a function of redshift and minimum mass. The upper panel does not take into account the effect of magnification bias, whereas the lower panel does. Those that are consistently magnified and pass a given threshold (bottom panel) are slightly magnified on average, whereas the top panel only shows a tiny constant μ ∼ 1−3% systematic residual magnification. This residual excess does not depend on wether the SPL or OBB method are used, or whether we properly integrate rays or use the Born approximation. This is likely because the replicates of the simulation box fill up the lightcone, which slightly increases the probability of rays leaving an over-dense region to cross other over-dense regions on their way to the observer. This residual magnification is however tiny for sight-lines that are populated by galaxies, and it completely vanishes for rays coming for random positions.

Fig. 10. Average magnification experienced by presumably foreground deflectors including (bottom) or excluding (top) for magnification bias effect that mostly affects the rapidly declining high end of the stellar mass function. Without magnification bias, a flat nearly unity mean magnification at all redshifts is recovered to within ∼1%. When the magnification bias is turned on, as expected in actual observations, no rapid rise is found (∼10% at z ∼ 1 for the most massive or luminous galaxies). Cuts in stellar mass are expressed in units of 1011 M⊙.

The massive end of the galaxy stellar mass function appears tobe significantly magnification biased. A ∼8% effect for galaxies at 0.6 ≤ z ≤ 1.2 and M_* ≳ 2 × 10¹¹ M_⊙ is typical. It can be as high at ∼20−50% at 1.5 ≤ z ≤ 2 for μM_* ≳ 3 × 10¹¹ M_⊙. A thorough investigation of the effect of this magnification bias when we try to place constraints on the high end of the z ≳ 2 luminosity function from observations is left for a forthcoming paper.

Taking the magnification bias into account, we now explore three fiducial populations of massive deflectors to highlight the changes induced in the projected excess density profiles. The first population consists of the CMASS galaxies at z = 0.54 and μM_* ≥ 1.7 × 10¹¹ M_⊙, the second case corresponds to the same lower limit on the mass, but pushed to z = 0.74. In both cases, the excess density is measured for source redshift z_s = 0.8. The last lens sample corresponds to the population of Hα emitters in the 0.9 ≤ z ≤ 1.8 redshift range that will be detected by the Euclid slit-less grism spectrograph above a line flux of ∼2 × 10⁻¹⁶ erg s⁻¹ cm⁻². About 2000 such sources per square degree are expected; therefore the 2000 most massive Horizon-AGN lightcone sources are picked in that redshift interval to crudely mimic an Hα line flux selection. To account for magnification bias, the selection was also made on μM_*, and the source redshift for this populations was set to z_s = 2. Results for these three populations are shown in the top panel of Fig. 11, where we distinguish the excess density profiles that include (dotted) or exclude (solid) for magnification. As anticipated, no significant change is obtained for the z = 0.54 CMASS-like sample (green), but differences are more noticeable as lens redshift increases, and on large scales (R ≳ 1 Mpc), we observe a 20−50% increase in ΔΣ, consistent with the large-scale linear scale-invariance bias model used by Ziour & Hui (2008). Between z = 0.54 and z = 0.74, galaxies of the same mass seem to live in halos of the same mass (very little evolution of the M_* − M_h relation), leading to no evolution of ΔΣ below ∼200 kpc. The only difference occurs farther out where the two-halo term starts to be important in this galaxy-mass correlation function. There, galaxies of the same mass at z = 0.54 and z = 0.74 live in rarer excursions of the initial density field, and are thus more highly biased, which leads to an increase of ΔΣ on the large scale. For the Euclid-like distant lens population, the trend is similar, and the amplitude of the magnification bias effect would suggest a bias of the lens population about 30% higher than it really is.

Fig. 11. Upper panel: effect of magnification bias on GGL for several high-z fiducial lens samples showing an increase in excess density ΔΣ (or tangential shear) for R ≳ 1 Mpc. Solid curves ignore the magnification, whereas dotted lines account for it. Lower panel: dependence of this effect on the source redshift. In both panels, cuts in stellar mass are expressed in units of 1011 M⊙.

The lower panel of Fig. 11 shows the evolution of the magnification-bias-induced excess density profile with source redshift for massive deflectors at z = 0.74. In principle, according to Eq. (22), the excess density should not depend on source redshift. However, the magnification bias favours over-densities in front of deflectors. The response of distance sources carrying shear to these over-densities will depend on the source redshift in a way that is not absorbed by Eq. (22). Hence, a scale-dependent distortion of the profiles is observed. The closer the source redshift from the deflector, the smaller the scale it kicks in. As already stressed by Ziour & Hui (2008), this hampers a direct application of shear-ratio tests with high-redshift deflectors (e.g. Jain & Taylor 2003).

6. Summary and prospects

Using two complementary methods for projecting the density or gravitational acceleration field from the Horizon-AGN lightcone, we propagated light rays and derived various gravitational lensing observables in the simulated field of view. The simulated area was 2.25 deg² out to z = 1 and 1 deg² all the way to z = 7. The effect of baryons on the convergence angular power spectrum P_κ(ℓ) was quantified, together with the two-point shear correlations ξ_±(θ) and the galaxy-galaxy lensing profile around massive simulated galaxies.

For cosmic shear, the inclusion of baryons induces a deficit of power in the convergence power spectrum of about 10% for 10³ <  ℓ <  10⁴ at z_s = 0.5. The amplitude of the distortion is about the same at z_s = 1, but is slightly shifted to roughly twice as high ℓ multipole values. On yet higher multipoles, the cooled baryons, essentially in the form of stars, produce a dramatic boost of power, nearly a factor 2 for ℓ ∼ 10⁵. As emphasised in Chisari et al. (2018), it is worth stressing that detailed quantitative statements on such small angular scales may still depend on the numerical implementation of baryonic processes.

For galaxy-galaxy lensing, the projected excess density profiles for a sample of simulated galaxies consistent with the CMASS sample at z ∼ 0.52 (analysed by Leauthaud et al. 2017) were found to be in excellent agreement. To properly analyse this signal around high-redshift deflectors, the magnification bias affecting the bright end of a population of distant galaxies was carefully taken into account, showing a large-scale increase of the signal as high as 30% beyond 1 Mpc for lenses at z ≳ 1. This type of effect is particularly pronounced for future samples of distant deflectors, such as the spectroscopic Euclid sources that are detected based on their Hα line intensity.

Peirani et al. (2019) have shown that the innermost parts of Horizon-AGN galaxies are consistent with strong-lensing observations of Sonnenfeld et al. (2013) and Newman et al. (2013, 2015) at z_lens ≲ 0.3. We intend to make more predictions on the optical depth for strong lensing in the Horizon-AGN lightcone with our implemented ray-tracing machinery. Likewise, in a forthcoming paper we will present the results of the deflection field applied to simulated images derived from the light that is emitted by the stars that are produced in the simulation, hence enabling the possibility of measuring lensing quantities (shear, magnification, etc.) in the very same way as in observations: shape measurement in the presence of noise, PSF, pixel sampling, photometric redshift determinations, realistic galaxy biasing, and more generally, directly predicted galaxy-mass relation, and also the intrinsic alignment of galaxies and their surrounding halos (Codis et al. 2015; Chisari et al. 2015, 2016).

Acknowledgments

The authors would like to thank D. Aubert for making his SPL code available to us. We acknowledge fruitful discussions with K. Benabed, S. Colombi, M. Kilbinger, and S. Prunet in early phases of the project. We also thank G. Lavaux, Y. Rasera, and M-A Breton for stimulating interactions around this project. We are also grateful to A. Leauthaud for constructive comments about the comparison with her GGL lensing results. CL is supported by a Beecroft Fellowship. This work was supported by the Agence Nationale de la Recherche (ANR) as part of the SPIN(E) ANR-13-BS05-0005 (http://cosmicorigin.org) ERC grant 670193, and AMALGAM ANR-12-JS05-0006 projects, and by the Centre National des Etudes Spatiales (CNES). NEC is supported by a RAS research fellowship. This research is also funded by the European Research Council (ERC) under the Horizon 2020 research and innovation programme grant agreement of the European Union: ERC-2015-AdG 695561 (ByoPiC, https://byopic.eu). This work has made use of the Horizon Cluster hosted by the Institut d’Astrophysique de Paris. We thank S. Rouberol for running the cluster smoothly for us.

References

All Tables

All Figures

	Fig. 1. 2D sketch of the past lightcone around redshift z = 1 (orange vertical line). Each mesh is a replicate of the Horizon-AGN simulation box (bounded with cyan lines). The tiling is performed all the way up to redshift z ∼ 8.
In the text

	Fig. 2. Distribution in the redshift – stellar mass plane of the 1.7 million galaxies in the Horizon-AGN lightcone. For guidance, the stellar mass limit for completeness is shown as well as fiducial cuts in mass that would be obtained with a flux-limited survey of various i-band limiting magnitudes.
In the text

	Fig. 3. Schematic view of the propagation of a light ray through a lightcone sliced into multiple discrete lens planes. The ray (red line) is deflected at each intersection with a thin lens plane. The deflection field is defined for each plane depending of the angular position on this plane αj(βj).
In the text

	Fig. 4. Schematic view of the problem induced by cells at the boundary of slabs j and j + 1, which become refined between time t and t + dt. Missing cells (devoid of dots) or cells in excess (overlapping dotted cells of different colour) can end up as lightcone particles. A fiducial light ray is drawn to illustrate the intersection length li between the ray and RAMSES cells.
In the text

Fig. 5. Comparison of zs = 0.8 convergence maps obtained with the OBB method (integration of transverse accelerations in cells, left) and with the SPL method (projection of particles onto convergence planes after adaptative Gaussian smoothing, right). The latter method applies a more aggressive smoothing that better erases shot noise. Inaccuracies of long-range deflections in the SPL method that are due to edge effects translate into a global shift for some galaxies, as compared to OBB. With this method, some missing acceleration cells occasionally produce modest artefacts on a small scale.

In the text

	Fig. 6. Left panel: convergence map generated with a $0\stackrel{″}{.}1$ pixel grid over a 2.25 × 2.25 sq. deg. field of view for a fiducial source plane at zs ∼ 1. Right panel: convergence map with a field of view of 1 sq. deg. at zs ∼ 2, and its corresponding convergence PDF showing the characteristic skewed distribution.
In the text

Fig. 7. Upper panel: convergence power spectra for source redshift zs = 1 (top) and zs = 0.5 (bottom) derived with the OBB (magenta) and the SPL (green) methods. The more aggressive smoothing of this latter method translates into a faster high-ℓ fall-off. The cyan curves (DM) only account for the DM component (rescaled by 1 + Ωb/ΩM). The red curve corresponds to the direct integration of the three-dimensional total matter power spectrum (Limber approximation) in the Horizon-AGN simulation (Hz-AGN). The blue curves are the direct integration of the Horizon-DM (DM-only) matter power spectrum (Hz-DM). Dashes reflect regimes where the three-dimensional spectra of Chisari et al. (2018) were extrapolated by a simple power law (extrapolation). The yellow lines show the particle shot-noise contribution at two different redshifts. Middle panel: ratio of the zs = 0.5 convergence power spectra obtained with the Born approximation and the proper multiple lens plane integration showing only very small changes up to ℓ ∼ 105. Bottom panel: ratio of the DM-only to total convergence power spectra at zs = 0.5,  1.0, and 1.5 for the SPL method.

In the text

Fig. 8. Upper panel: two-point shear correlation functions ξ+ (solid lines) and ξ− (dotted lines) for a fiducial source redshift zs = 0.5. We either correlate actual shear (red) or reduced shear (green) in the calculation to highlight the small-scale effect of baryons on this non-linear correction. Middle panel: ratio of shear correlation functions for the two cases. Bottom panel: ratio of shear correlation functions for a ray-tracing that only includes rescaled DM particles or all the components.

In the text

Fig. 9. Comparison of the GGL tangential shear signal around z = 0.55 Horizon-AGN galaxies (green curve surrounded by a light green ribbon) and the GGL observations of Leauthaud et al. (2017) (blue dots with error bars). Units are all physical (and not comoving!). Model uncertainties in the simulation past lightcone are roughly estimated by splitting the 2.25 deg wide field of view into four quadrants. They may be underestimated beyond 1 h−1 Mpc. Cuts in stellar mass are expressed in units of 1011 M⊙. Black, blue, and red curves show the GGL shear signal predicted with the SPL method for the total, DM, and baryonic mass distributions, respectively. For clarity, uncertainties are omitted. They are similar to those in the case of the OBB method (green).

In the text

Fig. 10. Average magnification experienced by presumably foreground deflectors including (bottom) or excluding (top) for magnification bias effect that mostly affects the rapidly declining high end of the stellar mass function. Without magnification bias, a flat nearly unity mean magnification at all redshifts is recovered to within ∼1%. When the magnification bias is turned on, as expected in actual observations, no rapid rise is found (∼10% at z ∼ 1 for the most massive or luminous galaxies). Cuts in stellar mass are expressed in units of 1011 M⊙.

In the text

	Fig. 11. Upper panel: effect of magnification bias on GGL for several high-z fiducial lens samples showing an increase in excess density ΔΣ (or tangential shear) for R ≳ 1 Mpc. Solid curves ignore the magnification, whereas dotted lines account for it. Lower panel: dependence of this effect on the source redshift. In both panels, cuts in stellar mass are expressed in units of 1011 M⊙.
In the text