Weak-lensing observables in relativistic N-body simulations

We investigate for the first time weak lensing in a high-resolution relativistic N-body simulation, which includes all relevant effects from general relativity. We integrate the photon geodesics backwards from the observer to the emitters. Our analysis is fully relativistic and non-perturbative for the scalar part of the gravitational potential and first-order in the vector part, frame dragging. We solve the Sachs optical equations and study in detail the weak-lensing convergence, ellipticity and rotation. We present the angular power spectra and one-point probability distribution functions for the weak-lensing variables, which we find are broadly in agreement with comparable Newtonian simulations. Our geometric approach, however, is more robust and flexible, and can therefore be applied consistently to non-standard cosmologies and modified theories of gravity.


INTRODUCTION
In the near future galaxy surveys like DESI, LSST, Euclid and others, will map out nearly the entire visible Universesee LSST Dark Energy Science Collaboration (2012); Abell et al. (2009); Aghamousa et al. (2016); Amendola et al. (2018); Santos et al. (2015); Walcher et al. (2019). They will measure redshifts, angular positions but also shapes and sizes of billions of galaxies with unprecedented precision. Because the weak gravitational lensing signal in particular is dominated by nonlinear structures, linear perturbation theory will not be sufficient to interpret these data sets, and numerical simulations are a vital tool which analysis pipelines will have to rely on.
Such simulations are commonly based on Newtonian gravity. However, recently some of us have proposed a new N -body code based on general relativity (Adamek et al. 2016a,b). In this code, all relativistic effects which can be relevant on scales much larger than the biggest black hole horizon are taken into account. In the past, structure formation in ΛCDM and with massive neutrinos have been studied with this code (Adamek et al. 2014(Adamek et al. , 2017, and it has been extended to simulate clustering dark energy (Hassani et al. francesca.lepori@unige.ch † julian.adamek@qmul.ac.uk ‡ ruth.durrer@unige.ch § chris.clarkson@qmul.ac.uk ¶ l.j.c.coates@qmul.ac.uk 2019) as well as f (R) gravity (Reverberi & Daverio 2019).
Recently some of us have also studied fully relativistic photon geodesics in this code to model the distance-redshift relation including clustering . The present work builds upon this development.
Here we want to study weak lensing. We compute the relativistic ellipticity (an observable closely related to the weak-lensing shear), the convergence κ as well as the rotation ω of images due to foreground structure nonperturbatively from a high-resolution simulation. We determine their angular power spectra for different redshifts as well at their one-point probability distribution functions. Because our analysis is non-perturbative, we can study effects like rotation, which is a purely non-linear effect which is absent within linear perturbation theory. A recent study has found (Deshpande et al. 2019) that a linear analysis will also be insufficient for interpreting the weak-lensing shear signal in upcoming surveys.
Following Perlick (2010), we develop a fully relativistic geometric description of lensing which is valid in arbitrary spacetimes, but we also make contact with the more familiar standard treatment in cosmological perturbation theory. This is the first study of weak lensing which uses a relativistic N -body simulation of "production scale"i.e. with multi-billion mass elements. Relativistic simulations of weak lensing have been done in the past for small problem sizes, see Giblin et al. (2017) for an example, but these do not probe the small-scale structure that is responsible for the bulk of the lensing signal. We find that our results agree very well with previous Newtonian and perturbative approaches within the present context of the standard ΛCDM cosmological model.
The remainder of this paper is organised as follows: In the next section we describe our methods. After a brief introduction to the general theory of weak lensing, we present the non-perturbative optical equations in Poisson gauge which are solved numerically in this work. In Section 3 we briefly describe our simulation and Section 4 is devoted to results. We present sky maps of the relevant lensing signals as well as angular power spectra and probability distribution functions. In Section 5 we conclude. Some more technical issues are relegated to four appendices.

Weak-lensing theory
When describing the properties of an infinitesimal bundle of light rays it is useful to introduce a "screen" that, at any given point along the central ray, is defined by two spacelike orthonormal screen vectors e µ 1 , e µ 2 that are normal to the photon four-vector k µ . Each choice of such a screen basis corresponds to a particular reference frame with four-velocity u µ that is normal to the screen vectors. However, one can show that the shape and size of an infinitesimal beam crosssection does not depend on this choice (Sachs 1961). In other words, the proper motion of the screen has no effect on an infinitesimal image that is projected onto it.
One particularly useful choice is then a screen basis that is parallel-transported along the central ray. Such a choice is called a Sachs basis. It draws a straightforward connection between the orientation of an image at an observer location and the aspect of the source in the screen basis. In cosmological spacetimes with a preferred global timelike vector field u µ one sometimes chooses a screen basis that is everywhere orthogonal to u µ , but in this case the screen is in general not parallel-transported. This choice can still be useful, e.g. in perturbative calculations where u µ can be chosen as the matter rest frame, see Pitrou et al. (2013); Marozzi et al. (2018) for examples 1 . However, a unique notion of matter rest frame is not always available in the real Universe, due to the existence of voids and overlapping matter streams. While we will later choose an observer at rest with the CMB (which defines a preferred frame u µ ) we shall discuss lensing purely from a spacetime perspective and work with a Sachs basis.
A small source mapped into a small image can be considered as a linear map from the screen at the source to the screen at the observer. The Jacobi map, D, describes the deformation of such a small image in the Sachs basis. In this treatment we follow Perlick (2010). We first introduce the 1 Note that Marozzi et al. (2018) and some other authors use the term "Sachs basis" for such a screen basis if its projection into the screen space is parallel-transported, which can be arranged for any choice of u µ . This differs from our convention that follows Sachs (1961); Perlick (2010), and the two conventions coincide if and only if u µ itself is parallel-transported. deformation matrix S of a light bundle as Here we use the geometric optics approximation within which the photon wave vector is the gradient of the eikonal so that SAB = e µ A e ν B kµ;ν is symmetric. The deformation matrix is related to spacetime curvature througḣ where a dot denotes the derivative with respect to the affine parameter of the photon geodesic and R is given by contractions of the Ricci and Weyl tensors with the photon fourvector and the so-called screen vectors (the aforementioned Sachs basis), see Perlick (2010) for more details. Written out in components and introducing the complex null shear σ ≡ σ1 + iσ2, one findṡ where Rµν and C αµβν are the Ricci and Weyl curvature tensors, respectively. The Jacobi map D then satisfieṡ In complete generality, we can parametrize D as Here R(α) denotes the rotation by an angle α, With this choice of parametrisation DA = √ det D = √ D+D− is the area distance, where D± = DA exp(±γ) are the eigenvalues of R(ω)D = Σ which is symmetric. The rotation by the angle χ rotates Σ into its eigen directions.
We could also have placed the rotation by the angle ω on the other side: for an arbitrary 2×2 matrix M there exists and angle ω such that R(ω)M = Σ is symmetric. In this case also M R(ω) = Σ is symmetric and Σ = R(−ω)ΣR(ω). The proof of this proposition is very simple and can be done by construction. In our application this would correspond in replacing χ by χ + ω.
It is convenient to introduce the complex ellipticity , Its phase encodes χ and its modulus is the product of the first and second eccentricity of the image of a circle under the Jacobi map. It is therefore closely related to the shear of the image (parametrised by γ1 and γ2) as we will see later.
The differential equation (5) can be written as [see Perlick (2010), eq. (27)] which can be rearranged into following equations for DA, and ω: It is interesting to note that the image rotation ω, which quantifies the net rotation with respect to a paralleltransported screen basis, is in full generality small at second order. This has also been found in Di Dio et al. (2019) where an analysis for scalar, vector and tensor perturbations was presented.
The boundary condition for D is D(s0) = 0 at the observer, which implies that all the quantities DA, σ, and ω can be set to zero in the initial (or rather final) conditions. Note also that χ is not a perturbative parameter. While γ and ω are expected to remain small, χ is simply the rotation into the principal axes which can be a large rotation even if the deformation of the light bundle is arbitrarily small. For this reason it is simpler to work with the complex ellipticity instead of γ and χ. If = 0 the ellipse described by the symmetric matrix Σ = R(ω)D becomes a circle and the angle χ is ill defined. The evolution equation (11) of , however, remains well-behaved at such singular points.
In the literature on weak gravitational lensing the Jacobi map is often written as [see e.g. Durrer (2008), eq. (7.18)] whereDA is the area distance of the background Friedmann metric and A is called the amplification matrix. The trace of A is parametrised by the convergence κ lin , and the shear of the image is given by γ1, γ2. Usually one sets ν = 0. At first order, these quantities are related to the area distance, the complex ellipticity and the rotation via Even though A is the most general ansatz for a 2 × 2 matrix, we prefer to work with (6) since the interpretation of its elements is straightforward: µ =D 2 A /D 2 A is the magnification, γ is the amplitude of the shear and R(χ) determines its principal axes, and finally ω determines the rotation of the image. Note that already at second order in the parametrisation (13) the quantities are mixed, for example the magnification becomes The fully non-perturbative expressions for γ1 + iγ2 and ν are Therefore, at next-to-leading order, both the trace-free symmetric part of A (parametrised by γ1 + iγ2) and the antisymmetric part of A (given by ν) are actually mixtures of image rotation and shear. Since the quantities κ, γ1, γ2 and ν have simple geometrical interpretations only at first order in perturbation theory, we shall work with the fully non-perturbative quantities DA or DA/DA, and ω.
As we discuss in detail in Section 2.2 and 2.3, in our numerical approach we compute DA, and ω without employing any approximations in the scalar sector, and taking into account also the leading-order corrections due to general-relativistic frame dragging.

Optical equations in Poisson gauge
Up to this point our discussion was completely general in the sense that we did not have to specify the spacetime metric. We will now specialise to the case of a perturbed Friedmann-Lemaître-Robertson-Walker (FLRW) metric, where the line element is written in Poisson gauge as Here, τ is conformal time and we assume a spatially flat background with scale factor a(τ ) and comoving Cartesian coordinates x i . In ΛCDM cosmology the metric perturbations are dominated by the two scalar gravitational potentials ψ and φ, which we will treat non-perturbatively, although they do remain small, of the order of 10 −5 at the scales that we probe. The divergence-free vector potential Bi, which captures the frame-dragging effect, is at least two orders of magnitude smaller than the scalar potentials (Lu et al. 2009;Thomas et al. 2015;Adamek et al. 2016a) and will therefore be treated at leading order only. This means we will neglect terms that are quadratic in Bi and also approximate e φ Bi e ψ Bi Bi. The spin-2 field hij, which is transverse and traceless, is even smaller on the scales that we are interested in, and we will therefore neglect it here. It would not pose a conceptual challenge to include it, but it would make our equations more cumbersome. For a classical point particle, let us introduce the canonical peculiar momentum per unit mass as where mp is the mass of the particle and L is the Lagrangian describing the geodesic motion of its comoving coordinate x i p . The particle's four-velocity vector is where q 2 = δ ij qiqj. A local Fermi frame is given by a spatial triad of orthonormal basis vectors s µ (i) that are orthogonal to u µ , i.e. gµν u µ s ν (i) = 0 and gµν s µ (i) s ν (j) = δij. For q = 0 (the cosmological rest frame) we can choose s µ (i) = δ µ i a −1 e φ . For the tangent vector of a photon geodesic, k µ , the null condition implies where n i is normalised as n i n j δij = 1 and denotes the direction of the photon path in the cosmological rest frame, Since we are interested in the properties of null geodesics it is convenient to employ a few conformal transformations in order to arrive at simple equations. First, we definek 0 ≡ a 2 e −2φ k 0 , such that the geodesic equations for the photon can be rewritten in terms of conformal time as where a prime denotes partial derivative with respect to conformal time and ∂i denotes a partial derivative with respect to x i . In order to also cast the geodesic deviation equations, (3) and (10)-(12) into a convenient form, we define the rescaled quantities as in .DA then evolves according to In order to solve the evolution ofσ we first need to construct a screen space, given by a pair of screen vectors e µ A (A = 1, 2) with gµν e µ A e ν B = δAB that are also orthogonal to k µ and are parallel-transported along the ray. We can write them in terms of two spacelike directionsẽ i A with δijn iẽj A = 0 and δijẽ i Aẽ j B = δAB as where theβA are related to the timelike direction to which the screen is orthogonal. For an observer in the cosmological rest frame we haveβA = 0 as a boundary condition, andβA will remain small (of the order of a metric perturbation) along the ray. The parallel transport of the screen basis e µ where we give the second equation here only for completeness as it will not be needed. This is becauseβA would only appear multiplied by terms of order Bi which can then be dropped.
The projection of the Weyl tensor onto the screen basis determines the evolution ofσ as The complex ellipticity and the image rotation ω are then obtained from integrating The set of coupled ordinary differential equations (25), (26), (28), (30), (32) and (33) fully determines our weaklensing observables once we have fixed the boundary conditions of each ray. At the observer location (assumed to be in the cosmological rest frame 2 ) the boundary conditions arẽ The boundary condition fork 0 is arbitrary and corresponds to some reference frequency that may be used to determine the redshift of the observed source. For convenience we set k 0 = ae −2φ−ψ at the observer, in which case the observed redshift of the source simplifies to where qi is the canonical peculiar momentum per unit mass of the source. The proper observed area distance DA is obtained by inverting eq. (27) at the source location. As we will explain in more detail in the next subsection, our ray-tracing method determines the boundary condition for n i for each source through a shooting algorithm. For any given n i a screen basis at the observer is constructed by choosingẽ i A appropriately.

Ray tracing
The weak-lensing observables DA, and ω are obtained by integrating the coupled ordinary differential equations backwards in time from the observation event. The other endpoint of the integration corresponds to an event on the observer's past light cone. A clear physical interpretation arises, for instance, if that event belongs to the world line of a source (e.g. a galaxy). Furthermore, in order to have a notion of observed redshift, one also needs to specify a four-velocity vector at the source location that identifies the source's rest frame. In this work we will consider cold dark matter elements (N -body particles) as the source population. This has the advantage that the bias is trivial, and that there are enough sources available even inside large voids. The issue of bias is very interesting and certainly deserves a separate study.
In the weak-lensing regime, for a given observation event and timelike source world line in the same causal patch of the Universe, there exists exactly one past-directed null ray that connects the two. The situation is different for strong lensing, where there can be multiple images, but we assume that the number of sources for which this happens is negligible. In order to identify that null ray (or one of the possible null rays in the rare case of strong lensing) we use a shooting algorithm similar to the one of Breton et al. (2019) -see in particular their Figure 1 for an illustration.
First, for each source we identify the spacetime event at which its world line would cross the past light cone of the observation event in the absence of metric perturbations. In fact, this is how the light-cone data for particles is constructed "on the fly" during the simulation run. The background FLRW model also provides us with an initial guess for the boundary condition of n at the observer. We then integrate the null geodesic equation in the fully perturbed metric (which is a separate simulation output as explained in the next section) until the ray passes close to the source world line. The source four-velocity allows us to construct a linear segment of the world line, which is sufficient for our purposes. Given the angular distance (estimated asDA), the spatial separation of the two geodesics corresponds to a "deflection angle" by which the shooting direction n must be corrected. The timelike separation corresponds to a "Shapiro delay" by which the source has to be moved along its world line in order to meet the null ray. Due to the non-relativistic peculiar velocities the latter correction is minute but we take it into account anyway. With these corrected boundary conditions we integrate the null geodesic equations anew, obtaining a ray that passes significantly closer to the source. The whole process is then iterated a few times until convergence is achieved.

Lensing potential
We have explained in the last two sections how weak-lensing observables can be constructed non-perturbatively from simulations. In this section we make contact with the wellknown perturbative approach to determine the (linear) convergence κ lin and shear γ1 + iγ2.
In the context of perturbation theory, weak lensing is often discussed in terms of the "lens map" that maps the direction vectors on the observer's sky to points on a distant "source plane" which is in fact a patch of a two-dimensional spacelike sub-manifold. It is important to emphasise the difference between this map and the Jacobi map that maps those direction vectors to points on a screen at the source plane. The two maps are only equivalent once we take into account the coordinate transformation between the coordinate basis used to label points on the source plane and the Sachs basis used to label the points on the screen, see Appendix A.
With sources at a fixed comoving distance r the points on the source plane can be labelled rθ i s , where θ i s is a direction vector on the two-sphere. In the weak-lensing regime the lens map provides a one-to-one correspondence between θ i s and a direction on the observer's sky, θ i o , and we can write where α i is called the deflection angle. Its definition is coordinate dependent because θ i s entirely depends on the chart used at the source plane. The conceptual utility of the deflection angle therefore also hinges on the subtle (often implicit) assumption that the perturbations of the induced metric on the source plane can be neglected. 3 In the literature, the Poisson gauge coordinates are often assumed which then gives a unique notion of α i , and from now on we shall follow this convention.
To first order in the metric perturbations ψ, φ and Bi the deflection angle can be computed by integrating eqs. (23) and (26), If we neglect the frame-dragging effect and furthermore insert the unperturbed trajectory in the integrand, the socalled Born approximation, which means where the index a now refers to coordinates on the twosphere (corresponding to θ i o and hereafter denoted as θ). The fact that αa can be written as a gradient field under this approximation motivates the definition of the lensing potential (Lewis & Challinor 2006) According to eq. (30) the Sachs basis does not rotate (with respect to the coordinate basis) at leading order in the absence of frame dragging. Therefore, the lens map and the Jacobi map are both approximated by and hence The complex shear is a spin-2 quantity which we can decompose into its "electric" component, γE, which has positive parity and its "magnetic" component, γB, which has negative parity, see Bernardeau et al. (2010). It is easy to verify that in the above case which is relevant within linear perturbation theory (in the absence of frame dragging and gravitational waves) γB ≡ 0. Once we include frame dragging or go beyond leading order and take into account e.g. post-Born corrections, the deflection angle is in general no longer a gradient map. We can decompose it into a gradient and a curl component, Here ε b a is the totally anti-symmetric tensor in two dimensions. Note that in two dimensions the curl of a vector is a scalar while the curl of a (pseudo-)scalar as defined above is a vector. Evidently, if Ω = 0 the lens map has an antisymmetric part, but we already know that the image rotation ω vanishes at first order in any spacetime. Indeed, if we look at the first-order calculation including Bi we have to take into account the fact that the Sachs basis rotates with respect to the coordinate basis. It is easy to verify (Dai 2014;Di Dio et al. 2019) that this rotation exactly cancels the one from the lens map at leading order, so that the Jacobi map is again symmetric. This also means that the polarisation (which is parallel-transported in the same way as the Sachs basis) rotates in the Poisson gauge, which has interesting consequences for CMB observables.
It should, however, be evident that the rotation of the Sachs basis does not change the shear signal at all for a statistically isotropic and uncorrelated source population (what makes the effect relevant to the CMB is the fact that the sources have intrinsic correlations). In linear theory Hirata & Seljak (2003) showed that in the "flat-sky approximation" where the C X are angular power spectra as defined in Sec. 4.2 below. Hirata & Seljak (2003) also find that the rotation of the lens map (i.e. its antisymmetric part) has the same power spectrum as γB, but they do not take into account the fact that the Sachs basis rotates. The rotation of the lens map is not directly observable (because it is gaugedependent), as opposed to the image rotation in the Sachs basis which can be observed, e.g. for a polarised source where polarisation and ellipticity are aligned. We note that frame dragging does produce an observable B-mode shear at first order, but no image rotation (in the sense that ω vanishes). However, in standard cosmological perturbation theory the frame-dragging effect itself only appears at second order (Lu et al. 2009), and hence its contribution is typically of the same order as post-Born corrections and other second-order contributions that do allow for image rotation to occur. In Newtonian simulations the lensing potential is often constructed from mass maps [see Hilbert et al. (2019) for a comparison of different state-of-the-art numerical methods].
Here one uses that in the Newtonian sub-horizon limit the two Bardeen potentials are equal and related to the matter over-density δm through the Poisson equation Together with equations (39) and (41) one then finds where the boundary term is usually dropped because it is ∼ 10 −5 and therefore much smaller than the typical value of κ. We note, however, that this approximation would induce spurious effects if one cross-correlates the lensing signal with the gravitational redshift of the sources.
It is worth pointing out that this type of lensing analysis is only applicable if the gravitational fields are indeed generated by matter according to eq. (46), and therefore precludes modifications of gravity or the possibility that in general there can be sources to metric perturbations other than nonrelativistic matter. Our geometric approach, working directly with the metric perturbations, is therefore more robust and flexible.
Furthermore, it has recently been demonstrated (Fidler et al. 2017;Adamek & Fidler 2019) that the weak-lensing analyses of Newtonian N -body simulations would have to include a gauge correction that is commonly neglected. This is due to the fact that the coordinate system of the Newtonian simulations does not coincide with the Poisson gauge for which the lensing calculations have been developed. The correction appears at large angular scales and can be incorporated into a modified lensing potential. This is not an issue in our relativistic simulations as we use Poisson-gauge coordinates consistently throughout.
We want to study how well the first-order lensing potential (neglecting frame dragging) characterises the full weaklensing signal. To this end we construct maps of the lensing potential at fixed comoving distance (corresponding to fixed redshift in the background model) by numerically integrating eq. (39). Due to the Born approximation this can be done very efficiently directly in pixel space, as the gravitational potentials φ and ψ are already conveniently pixelised (this is explained in more detail in the next section). Once the lensing potential Ψ is computed, we generate a map of κ lin by solving eq. (41). This functionality is readily available in the public release of gevolution version 1.2. 4 Appendix D presents a direct validation against other weak-lensing codes within the Newtonian approximation.

SIMULATION
Our numerical results, presented in the next section, are based on a large N -body simulation using the relativistic code gevolution (Adamek et al. 2016a). The simulation has 7680 3 (almost half a trillion) mass elements in a cosmological volume of (2.4 Gpc/h) 3 , which gives a mass resolution of 2.6 × 10 9 M /h. We use a standard cosmology with Ωm = 0.312, Ω b = 0.048, h = 0.67556, and massless neutrinos with N eff = 3.046. We generate the linear transfer functions of baryons and cold dark matter at initial redshift zini = 127 with class (Blas et al. 2011), and choose a primordial amplitude of scalar perturbations As = 2.215 × 10 −9 at the pivot scale k * = 0.05 Mpc −1 and spectral index ns = 0.9619. This corresponds to σ8 = 0.8488.
We choose an observer located in the corner of the box and consider a pencil beam on the past light cone, with a 450 sq. deg. field-of-view centered around (l = 33.4 • , b = 48.6 • ) in "galactic coordinates" aligned with the principal axes.
This specific geometry was chosen such that the conical footprint of the mock survey does not contain any replications up to redshift z 3.25. In other words, the pencil beam is allowed to pass through the periodic domain more than once, but without any self-intersection. In a small region up to a comoving distance of 275 Mpc/h we retain the past light cone on the full sky, in order to have more data around the apex of the pencil beam.
The data on the light cone are recorded as follows. The comoving positions and the peculiar momenta of N -body particles are recorded at the time when their world lines would intersect the past light cone if the metric was unperturbed. Up to a small Shapiro delay this intersection will be very close (in terms of four-dimensional space-time distance) to the intersection with the true light cone of the perturbed metric. The linear segment of the world line that can be constructed from the particle phase space coordinates is therefore sufficient for locating that latter intersection point in the later analysis.
The true past light cone can only be determined from a complete knowledge of the metric, which is only available at the end of the simulation. The metric perturbations are initially sampled on a Cartesian mesh with a comoving spatial resolution of 312.5 kpc/h. However, for the purpose of constructing a light cone it is more convenient -and computationally much more efficient -to work in spherical coordinates. We therefore define a second coordinate mesh, consisting of concentric spherical shells with the observer at the origin. The shells are separated radially by our base resolution of 312.5 kpc/h, and are each pixelised using the HEALPix framework (Górski et al. 2005). The number of pixels is chosen for each shell separately with the requirement that the area of each pixel is always less than (312.5 kpc/h) 2 . Since the number of pixels can only be changed discontinuously this means that our HEALPixmesh is generally denser than the Cartesian mesh.
For each time step of the simulation we first identify the corresponding comoving distance interval on the past light cone in the unperturbed spacetime. We then record the perturbed metric on all the shells within that distance interval by "triangular-shaped particle" interpolation from the Cartesian mesh onto the pixel locations. The same is done for the preceding as well as the following time step, such that we effectively retain a finite portion of the fourdimensional spacetime around the light cone. This allows us to reconstruct the true light cone without any approximations, by tracing null-geodesics backwards in time from the observer.

RESULTS
In this section we discuss the main results of our work. The outline of the procedure is as follows: (1) we perform raytracing over a subset of the particles in our simulation, (2) we compute convergence, shear and rotation angle for each source in our catalogue, (3) we project these fields onto a pixelised sky map, (4) we compute the angular power spectra for the maps of the fields for different redshifts, (5) we compute the one-point probability distribution function (PDF) for the values of the fields.
For our analysis, we run the ray tracer over ∼ 135 mil-  Table 1) are indicated as shaded colour bars. lion particles, randomly selected from the full N -body ensemble. The reason to use a particle catalogue rather than a halo catalogue is twofold: first, the N -body particles are an unbiased tracer of the dark matter and second, the larger sample size offers better statistics compared to a halo-based study.
The normalised redshift distribution of the particle catalogue is show in Fig. 1. The particle catalogue covers a redshift range between z = 0 and z 3.25. At low-redshift, up to z 0.09, the particles are distributed over the whole sky, while for 0.09 z 3.25 the particles are distributed over a sky fraction f sky = 0.01. We focus our analysis on five redshift bins (see Table 1). These bins are also outlined in Fig. 1 as vertical colourbars.
The remainder of this section is organised as follows. In Sec. 4.1 we describe the method to obtain field maps from the optical parameters of the source population, in Sec. 4.2 and Sec. 4.3 we discuss the weak-lensing angular power spectra and PDFs, respectively, as well as their redshift dependence. In Sec. 4.4 we compare the non-linear convergence obtained from the ray tracer to the linear convergence computed within the Born approximation (i.e. along the unperturbed light path), from the lensing potential map.

Weak-lensing maps
The ray tracer numerically computes the area distance, the complex ellipticity and the rotation angle of the image, at the observed redshift and angular position in the sky, for all  the dark matter particles in the catalogue. In our procedure, we work only with observable coordinates, the observed sky position and redshift. Once these observables have been computed, neither the comoving position nor the peculiar velocity of the N -body particles is used in the analysis, as these would be specific to the Poisson gauge and therefore have no invariant meaning.
Convergence and shear, at the observed position of each particle, are estimated from eq. (14). The convergence depends only on the area distance, while the shear coincides with the ellipticity, up to a factor 4. The magnification is a function of the convergence alone and it is computed from eq. (15).
Convergence, ellipticity (respectively its absolute value), magnification and rotation angle are projected onto a HEALPix map in the following way: at fixed redshift z, we select all the particles inside the redshift bin of thickness δz and the value of the field in each pixel is computed as the average value from all the particles inside the pixel volume.
In this section, we adopt a map resolution of N side = 2048 which corresponds to an angular resolution of ∼ 1.71 arcmin. The number of pixels over the whole sky is N pixel = 12 × N 2 side which means that our survey area provides approximately 0.5 megapixels. In Appendix C we discuss the impact of angular resolution on our results.
In Fig. 2 we show the maps for the convergence, the amplitude of the ellipticity, | |, the magnification and the rotation angle at mean redshift z = 1.5. The maps highlight that convergence and magnification show very similar patterns. Ellipticity and rotation angle show a significantly different small-scale behaviour. However, the regions in the map with higher ellipticity and ω (in absolute value) correspond to the regions in the maps with larger convergence and magnification. Despite the fact that ellipticity and convergence have nearly identical spectra, as expected from the linear analysis, their maps look quite different with elongated "filaments" in the ellipticity map which are not present in the convergence map. These patterns are certainly related to similar ones that have been found for the deflection angle (Carbone et al. 2008;Watson et al. 2014).

Angular power spectra
Angular power spectra for our observables can be obtained through spherical harmonic decomposition. Convergence and rotation, being scalar and pseudo scalar quantities, can be expanded in spherical harmonics where κ m and ω m are the coefficients of the spherical harmonic decomposition.
In a similar way, the ellipticity which is a spin-2 object, can be expanded in ±2 spin weighted spherical harmonics (Chon et al. 2004) where E and B denotes the E and B modes of the ellipticity, respectively. The angular power spectra, at fixed redshift 5 , are defined as follows 5 We omit the redshift dependence for brevity. Relative difference between the κ and E /4 power spectra in %.
We plot the spectrum of E /4 and not E since within linear perturbation theory this spectrum agrees with the κ-spectrum.
The dashed line corresponds to 4% relative difference.
The maps we obtain from our simulations and raytracing cover about 1% of the sky. Therefore, in the spherical harmonic decomposition outlined above, masking effects are crucial and need to be corrected for.
In fact, the standard HEALPix routine for the angular power spectrum estimation -anafast -sets to zero the masked pixels before the spherical harmonic decomposition and would introduce a large bias in our analysis.
In order to correct for the masking effects, we estimate the angular power spectra for our maps with the code PolSpice, which implements an unbiased estimator for the power spectra based on the estimation of weighted correlation functions. The details on this estimator can be found in Szapudi et al. (2000); Chon et al. (2004). The estimated power spectra are binned linearly in with bin size ∆ = 100, which corresponds to 1/f sky . In this way, we reduce the cosmic variance errors and spurious oscillations in the spectra.
In Fig. 3 we show the estimated angular power spectra for the convergence and the ellipticity E-modes. In the perturbative approach, the two power spectra are related through 6 Therefore, we expect C E /4 ≈ C κ on large multipoles. On the top panel, we compare the numerical convergence and ellipticity spectra with the expected power spectrum estimated from class, where non-linearities are taken into account through the Halofit corrections to the power spectrum (Takahashi et al. 2012). Note that at 1000 the spectra estimated from our simulations show a significant loss of power. In Appendix C we investigate the cause of this suppression and we find that it is mainly caused by the finite resolution of the simulation mesh. On the bottom panel, we show the relative difference between the convergence and ellipticity E-modes (rescaled by a factor 4). The relative difference is smaller than 5% at all scales and for all the redshift bins.
As discussed in Sec. 2.4, at linear order in perturbation theory shear B-modes and the rotation angle ω are identically zero if frame dragging is neglected (it is usually considered as a second-order effect). However, post-Born corrections do source both these quantities and the curl-mode of the ellipticity also receives contribution from the vector perturbations in the metric.
Since the rotation ω is a pseudo-scalar, extracting its angular power spectrum from a masked map is not particularly problematic. In Fig. 4 we show the angular power spectrum for the rotation, at different redshifts. Data points refer to the results from our simulations, while continuous lines are the theoretical prediction from secondorder post-Born lensing calculations (Pratten & Lewis 2016;Marozzi et al. 2018). The theoretical predictions are computed using the post-Born lensing module implemented in the Code for Anisotropies in the Microwave Background (CAMB) 7 . Similarly to the convergence and ellipticity spectra, non-linearities are modelled with the Halofit prescription (Takahashi et al. 2012).
On large scales we find good agreement between the result of our simulations and the post-Born prediction. On small scales, we notice that ω suffers a power-loss due to finite grid resolution effects, similar to what we observe for the convergence and ellipticity spectra.
The rotation power spectrum has been computed in 7 https://github.com/cmbant/CAMB previous works, based on Newtonian simulations, using a multiple-lens-plane approximation (Becker 2012;Takahashi et al. 2017;Fabbian et al. 2019). Our results qualitatively agree with previous results in the literature, where a fewpercent agreement was found between the spectra extracted from Newtonian simulations and the second-order post-Born prediction. At linear order (and neglecting frame dragging), the ellipticity B-modes and the rotation ω spectra both vanish. Both are generated only at second order and we therefore expect them to be of comparable amplitude (yet not identical -note in particular that the ellipticity receives a contribution that is linear in the frame-dragging potential, while the rotation does not, see also App. A). However, we are not able to extract the B component above the noise from our numerical results for the ellipticity. Considering the four orders of magnitude and more difference in amplitude, we believe that this is due to masking and pixelisation effects that inevitably lead to some leakage between E and B .

Probability distribution functions
In this section, we show the results for the one-point Probability Distribution Function (PDF) for the convergence, the magnification, the ellipticity and the rotation. Our results are complementary to previous work on this topic [see e.g. Takahashi et al. (2012), where a detailed study of the lensing PDFs is presented for high-resolution and small-volume simulations].
The PDFs discussed is this section have been computed in 'pixel-space', i.e. the fields have been estimated on a HEALPix map, as described in Sec. 4.1, and we computed the distribution of the field values in the pixels.
We considered five redshift bins with mean redshifts in the range 1 ≤ z ≤ 3 as detailed in Table 1 and a map resolution N side = 2048.
In Fig. 5 we show the normalised PDF for κ at different redshifts. The shape of the κ-PDF qualitatively agrees with previous results in the literature (Takahashi et al. 2012): the distribution broadens at high z and the peak position shifts toward smaller values of κ. The distribution is strongly skewed towards positive values of κ. This is due to the fact that κ cannot be more negative than its value for an empty beam, κempty given in eq. (58) below. Since modelling the non-Gaussian shape is crucial in order to avoid systematic bias in distance-redshift relation measurements, we have studied the shape of the distribution and its redshift dependence.
We consider the modified lognormal distribution, first proposed in Das & Ostriker (2006) and further validated against simulations in Takahashi et al. (2012): where the convergence of the empty beam, κempty is  Here r(z1, z2) is the comoving distance from z2 to z1 and r(z) = r(0, z). The fitted model parameters are the normalisation Nκ, the parameter σκ, which determines both the width of the distribution of x = κ/|κempty| and the peak position, and the modification to the lognormal distribution Aκ.
The values of the fitted parameters are reported in Table  2. We omit the value of the normalisation. Note that the modification to the log-normal distribution does not depend strongly on redshift, while σκ decreases at high z. This tells us that the width of the distribution as a function of x = κ/|κempty| is decreasing towards higher redshift. However, κempty is is larger at high z so that the width of the κdistribution is increasing (see Fig. 5).
In Fig. 5 we show the result of our fits at z = 1, 1.5, 2, 2.5. The modified log-normal distribution is able to reproduce accurately the shape of the convergence PDF near the position of the peak and for small convergence values, while the fitting function underestimates the PDF for large values of κ. The steep decent of the distribution for κ → κempty is very well captured by (57). These results agree with the analysis in Takahashi et al. (2012).
In Fig. 6 we show the PDF of the magnification. Note that it is clearly related to the convergence PDF since the  two quantities, related by eq. (15), agree up to a shift by one and a factor two within linear perturbation theory.
In Fig. 7 we show the PDF of the real part of the ellipticity, [ ] ≡ ( + * )/2. In linear perturbation theory the convergence and the real part of the ellipticity have the same distribution. However, non-linearities affect the two distributions in quite different ways: while κ is highly skewed (see Fig. 5) both the real and imaginary part of have the same symmetric, albeit non-Gaussian, distribution. Even though the power spectra of κ and /4 are indistinguishable within our error bars, the PDF's of these quantities are clearly different. In Fig. 7 we compare the distribution for [ ] to a Gaussian fit (green dashed lines). Close to the peak,  Table 3. Fitted parameters for the [ ]-PDF and the | |-PDF. σ | | is the standard deviation estimated from the fit of | | to a Rayleigh distribution, while σ [ ] and α are the parameters of the modified Gaussian PDF given in eqs. (59) and (60).
at [ ] = 0, the PDF is well approximated by a normal distribution. However, the distribution has clearly visible, symmetric non-Gaussian tails. The model for the probability distribution dP/d [ ] can be improved by introducing a polynomial correction to the normal distribution of the form where N [ ] is the normalisation, σ [ ] is the standard deviation of the Gaussian, and α is a parameter that models the non-Gaussian tail of the distribution. The fit for the model in eq. (59) is show in Fig. 7 (orange lines). The modified-Gaussian model is a better fit than a simple normal distribution for large values of the ellipticity. The values of the best-fit parameters are reported in Table 3. The parameter α decreases significantly with increasing redshift which is expected as more of the Gaussian signal from higher redshift is accumulated. At the same time, the average lensing signal and therefore σ [ ] increases with redshift. Fig. 8 shows the probability distribution for the absolute value of the ellipticity, | |. The PDF for | | is highly non-Gaussian already at high redshift.
Because of statistical isotropy, the joint distribution of [ ] and [ ] has to be a function of | | = [ ] 2 + [ ] 2 only. We can therefore try to construct such an appropriate joint distribution given the knowledge of the fit for the marginalised distribution for [ ]. Integrating over the phase angle then directly determines the PDF for | |. In the Gaussian case the latter is given by a Rayleigh distribution. These considerations motivate the following form of fitting function, with σ [ ] and α determined through a fit for [ ], while N | | is a normalization parameter that we fit to the data. Using the statistical isotropy argument it is easy to see that this ansatz is consistent with eq. (59), while for α = 0 the standard Rayleigh distribution is recovered. In Fig. 8 we compare the PDF estimated from our simulations to a fit to a Rayleigh distribution (green dashed lines) and to the modified Rayleigh distribution given in eq. (60) (orange solid lines). The unmodified Rayleigh distribution provides a good fit only around the peak of the distribution, while there is a significant discrepancy between data and model both for large value of | | and for | | approaching zero. The Rayleigh distribution underestimates the probability of having large values of the ellipticity amplitude at all redshift. It compensates this somewhat by having a slightly larger width than the modified distribution, σ | | > σ [ ] , but not very successfully. At high redshift the probability for [ ] and [ ] of being both close to zero is significantly reduced compared to the Gaussian case. The modified Rayleigh distribution in eq. (60) is a much better fit to our data for large values of the ellipticity. However, for | | → 0, it reduces to the Rayleigh distribution and therefore still overestimates the probability of having pixels with very small ellipticity.
In Fig. 9 we show the distribution for the rotation angle ω, which is not close to Gaussian at any redshift. In Appendix B we provide a derivation for the expected PDF: for small values of ω it can be well modelled as an exponential, while the tails follow a product-Gumbel distribution [see eq. (B16)] with the respective parameters λmσ 2 lens and s given in Table 4.
Note that here we show the probability distribution functions (PDFs) of the lensing quantities κ, and ω while we can observe these quantities only at the position of a luminous source, where we, e.g. can measure the ellipticity of a galaxy. In our analysis we chose the source population to be drawn randomly from the N -body ensemble (thereby avoiding the additional complication of galaxy bias). What we obtain is therefore closer related to ρ(z) etc., which means that the observables are mass-weighted in each pixel. Nevertheless, since the lensing quantities are integrated effects and  Table 4. Fitted parameters for the ω-PDF models discussed in Appendix B. The parameters s represents the widths of the Gumbel distribution that characterises the tails, while λmσ 2 lens is the scale of the exponential near the maximum. get contributions from density fluctuations along the entire light path, it is most probably a very good approximation to neglect their correlation with density fluctuations at the source redshift z and our PDFs should be representative for the ones which are truly measured.

Convergence from lensing potential
As a consistency test for our results, we finally compare the angular power spectrum and the map for the convergence estimated using the full ray tracing method outlined in the previous sections (from DA) and in the Born approximation, from the lensing potential map (from Ψ). We described the method to obtain the convergence in the Born approximation in Sec. 2.4.
In Fig. 10 we show the angular power spectra of the convergence computed with the two methods. The methods agree within 2%. On scales 1000, where our computation is reliable, the Born approximation systematically underestimates the power by about 1.5%. However, it is not entirely clear whether the difference is physical or a computational artifact. We note in particular that the redshift selection is based on observed redshift for the ray tracer, while the method of the lensing potential assumes a distance-redshift relation. Fig. 11 shows a zoom-in of the κ-maps obtained with the ray tracer (top panel) and from the lensing potential (bottom panel) for the redshift bin at z = 1.5. As visual reference points we also show the observed positions and appar- Figure 10. The angular power spectrum of the convergence obtained with the full ray tracing (i.e. based on the non-perturbative area distance D A ) and within the Born approximation (from the lensing potential Ψ). In the bottom panel we show the relative difference between the two spectra in %, including error bars. ent sizes of dark matter halos with more than 3 × 10 14 M /h and observed redshifts below 1.5. Here, the apparent size is inferred from the proper virial radius and the observed area distance DA. While the agreement between the two methods of computing κ is generally excellent, one may notice at close inspection that the peaks of the κ-map from the lensing potential do not line up perfectly with the positions of the clusters (which is the truly observed one in both panels). This error is due to the Born approximation, and we can estimate it to be of the order of the arcminute for this redshift. Fig. 12 shows the map of the image rotation ω for the same redshift bin. The strongest rotation is not found at the centres of the projected dark matter halos, but typically in scattered places at the periphery of lensing peaks. The pattern can be understood quite well from considering a simplistic model where the rotation is generated by the composition of two linear lens maps, as illustrated in Fig. 13. In the sketched example we place a background lens off-axis and close to the maximum of the shear distortion of a foreground lens. The composition of the two lens maps (done here in the Born approximation) has an anti-symmetric component, i.e. some non-vanishing image rotation. In the shown configuration the rotation map has a parity-odd pattern with strong dipolar and noticeable quadrupolar components. Similar features (though less symmetric) can be seen in Fig. 12.

CONCLUSIONS
In this paper we have presented the first study of weak lensing by cosmological structure using a high-resolution relativistic N -body simulation. To this end we have developed a methodology to solve the Sachs equations along the true photon trajectory. Our ray-tracing code is independent of the underlying cosmology. It only assumes that photons move along geodesics and works for arbitrary forms of the Figure 11. The upper panel shows a small patch of the κ-map constructed from the observable properties of sources in the redshift bin at z = 1.5 as computed with the ray tracer. The superimposed circles show the observed position and apparent size of dark matter halos with masses above 3 × 10 14 M /h in the same field of view. The angular diameter subtended by the largest halo (close to the centre) is ∼ 18 , while the smallest halo (towards top left) is only ∼ 5 across. For comparison, the lower panel shows the map of κ lin constructed from the lensing potential.
scale factor, the Bardeen potentials and the vector perturbations Bi, remaining completely agnostic about how these perturbations are generated. It is nonperturbative in the scalar part and perturbative to first order in the vector part. We have computed the observed positions and redshifts of a source population and the optical scalars DA (or rather κ = 1 − DA/DA), the ellipticity and the rotation ω for Figure 12. For the same patch as shown in Fig. 11 we plot the map of the image rotation ω -the colour scale ranges from −0.002 (blue) to +0.002 (red) radians. Figure 13. Illustration of the generation of image rotation through the composition of two linear lens maps. The left and right panels show the lens maps of two individual lenses, modelled here through rotationally symmetric Gaussian convergence profiles with misaligned centres of symmetry. The composition of the two maps is shown in the centre panel, and the colours indicate the image rotation (red and blue corresponding to opposite-sense rotation).
the respective photon geodesics. We have also studied the angular one-point distributions of these variables and provide the first theoretically motivated fitting functions for the rotation.
Despite the fact that we nowhere invoke a Newtonian approximation, our results overall agree well with previous ones obtained with Newtonian codes. Within the ΛCDM cosmological standard model we therefore find no indication of large non-linear corrections coming from the gravity sector -the non-linearity is completely dominated by the clustering of matter, while the gravitational fields remain weak in the Poisson gauge. The strong clustering of matter however means that the weak-lensing observables deviate significantly from linear perturbation theory in which ω = 0 and both the κ and one-point distributions are Gaussian.
Our resolution has not been sufficient to extract the B-part of the ellipticity from the numerical noise.
While it is included in our results for the first time, we have not systematically studied the effect of the (small) vector perturbations which are induced at second order from the scalar perturbations and are completely absent in a Newtonian treatment. To see them we would need to construct observables that specifically isolate them from other secondorder contaminations. We shall study this problem in future work. Carbon footprint: In this work we re-used existing simulation data. Additional post-processing consumed about 3000 kWh of electrical energy, which has an estimated 8 impact of 600 kg CO2. This project also accounts for one return flight Geneva-London in economy class, causing emissions of approximately 9 180 kg CO2. In this appendix we derive the relation between the Jacobi map and the amplification matrix, which is the Jacobian of the lens map, in full generality, i.e. also for strong deflections. We first introduce the deflection angle α defined via the lens map, see Section 2.4, denoting the image position in our chosen coordinate system by θo = θ(s0) and the source position by θs = θ(s), we set

ACKNOWLEDGEMENTS
Here s is the affine parameter along the photon trajectory, s0 is its value at the observer position and θ is the direction vector tracking the spatial photon position given in some coordinate basis (∂0, ∂i), and α is the deflection angle. Note that the entire definition here is coordinate dependent and cannot be made intrinsic since the quantities α i (s), θ i (s0) and θ i (s) are best understood as coordinate vectors with respect to the coordinates ∂i and not as geometrical elements in a tangent space. In the given coordinates, the amplification matrix is now simply the Jacobian of this map, rescaled 8 Our conversion factor of 0.2 kg CO 2 kWh −1 is taken from Vuarnoz & Jusselme (2018), with the angular diameter distance of the background spacetime, Clearly, this construction makes sense only for lensing due to perturbations on some background, where the photon geodesic at s and the observer position can be covered by one single coordinate patch.
To obtain the Jacobi map from A i j we have to consider a neighbouring geodesic, say which has position θ(s0)+ (s0) at the observer and θs + s = θ(s) + (s) at the source in our chosen (arbitrary) coordinate system. The Jacobi map expresses s to first order, A Sachs basis at the observer ea(s0) , a ∈ {1, 2} is a 2d orthonormal basis of the plane (the "screen") normal to the photon direction θ(s0) and the observer 4-velocity u(s0), and it is defined along the entire photon geodesic by parallel transport. Therefore ea(s) is unique up to a global rotation. Let us denote by e a i(s) the transformation from our coordinate basis ∂i to a fixed Sachs basis ea(s) and by e i a(s) its inverse, so that the basis transformations are given by ea(s) = e i a(s)∂i and ∂i = e a i(s)ea(s). If two (say ∂1 and ∂2) of our arbitrary coordinates ∂i happen to be aligned with the Sachs basis and the deflection is sufficiently weak that we may consider θ n e3 constant (Born approximation), and the peculiar velocity is small, i.e. approximately u ∝ ∂0, we can neglect the evolution of the Sachs basis and 10 e a i(s) δ a i . In this case A i j in the 2d space normal to n agrees with D a b . Interestingly, in Poisson gauge this is the case at first order in perturbation theory if there are only scalar perturbations. The main reason for this is that the spatial part of the metric in Poisson gauge is conformally flat and parallel transport does not rotate the Sachs basis along the photon geodesic [see Di Dio et al. (2019) for details]. The opposite extreme case is geodesic light-cone gauge (Gasperini et al. 2011), where the angular coordinates along a photon geodesic are chosen to be its angular position at the observer so that, by definition α ≡ 0. In these coordinates the amplification matrix is proportional to the identity and the Jacobi map is given by the first term of (A7). In this case the metric is far from conformally flat and parallel transport rotates the Sachs basis along the photpon geodesic so that e a i(s) has a non-trivial s-dependence. Of course for a fixed Sachs basis at the observer, ea(s0) one finds the identical Jacobi map in both coordinate systems.

Our deviation vectors in the Sachs basis are
In general, the Sachs basis varies along the line of sight and D a b is not of the simple form of a gradient map like the Jacobian of eq. (43). This then implies that the power spectrum of the rotation and the B-mode of the Jacobi map do in general not agree. Hirata & Seljak (2003) showed that for a gradient map, A i j (θ) = const + α i ,j (θ) in the flat sky approximation, we always have C κ = C γ E and C ω = C γ B = 0. Since the Jacobi map agrees with the coordinate-dependent amplification matrix only at first order in perturbation theory (Born approximation) in Poisson gauge, these relations are not expected to hold beyond first order perturbation theory.

APPENDIX B: SHEAR AND IMAGE ROTATION AT LEADING ORDER
In this appendix we derive some expressions for the ellipticity and the rotation ω at leading order. For this purpose it is sufficient to work in the Born approximation (n i constant). We can also neglect Shapiro delays and set dτ = −dr, with r denoting the comoving distance from the observer.
Let us begin by inspecting eq. (32) which at lowest order can be approximated as where we drop contributions from Bi which is only generated at second order in ΛCDM cosmology. We introduceψ0 as a shorthand for the complex source term that is generated by Weyl curvature. Sinceσ andψ0 are first-order quantities, we only need the background solution ofDA which is given byDA r or, equivalently, DA r/(1 + z). Thereforẽ We now continue with the same reasoning by approximating eqs. (33) at leading order as so that by inserting the previous result we get an estimate of the complex ellipticity of For the rotation ω we then get which, after some manipulation, can be rearranged to We observe that ω is given by a double integral ofψ0 with a kernel that vanishes in the coincident limit and generates maximum signal for a pair of lenses that are placed at one third and two thirds of the comoving distance between source and observer. More precisely, the kernel takes its maximum value, r 2 /27, at r1 = r/3, r2 = 2r/3, and is parity-odd under r1 ↔ r2.
The last expression can be used, for instance, to estimate the two-point function of ω, which then depends on the first four n-point functions of the gravitational potentials φ and ψ. We can also get some insights into the one-point distribution of ω that we discuss next.
Due to the fact that the kernel vanishes for r1 = r2, large values of ω can only be generated by pairs of gravitational lenses that are well separated along the line of sight. We can therefore assume that these lenses are uncorrelated to a good approximation. Furthermore, if we neglect cosmological evolution, statistical isotropy implies that [ψ0] and [ψ0] are drawn from the same distribution.
Let us now first consider the simplest case of two lens planes, obtained by splitting the line of sight in half and treating the source terms on both legs as uncorrelated random variables. A rough estimate of ω is then given by where Xj, Yj denote the random variables on lens plane j, respectively given by the real and imaginary parts of its "effective" source termψ0. The factor 1/192 comes from integrating the kernel, approximatingψ0 by these piecewise constant "effective" values.
In the case where the Xj, Yj are Gaussian with r.m.s. amplitude σ lens , we can compute the probability density function of ω explicitly, Note that σ lens is the typical amplitude of the effectiveψ0 and therefore has dimensions of inverse area. The above expression would suggest that the r.m.s. value of ω scales as ∼ r 4 , but this is too naive because σ lens would depend on the effective baseline that contributes to each lens plane. We can try to resolve this issue by considering multiple lens planes spaced regularly along the line of sight. Their separation can then be chosen to be commensurate with the correlation length ofψ0, so that the planes can still be considered approximately uncorrelated. In this case σ lens should be approximately constant and can be estimated as the typical scale of the gravitational potential (∼ 10 −5 in cosmological settings) divided by the square of a characteristic length scale, expected to be similar to the correlation length ofψ0. In this approximation of multiple lens planes with random Xj, Yj we get where n is the number of lens planes and the coefficient matrix A jk is computed from eq. (B6) assuming thatψ0 is approximately piecewise constant, In this expression, rj denotes the comoving distance to lens plane j, and ∆r = r/n is the length of the distance interval assigned to each lens plane.
In the Gaussian case the probability density function of ω can be calculated explicitly, assuming for simplicity that the Xj, Y k are all independent and identically distributed 11 with r.m.s. amplitude σ lens , The Gaussian integrals can be solved by defining a 2ndimensional random vector Z ≡ (X1, . . . , Xn, Y1, . . . , Yn) and a 2n × 2n-matrix M as such that Next we make a change of variables to rotate into the orthonormal eigenbasis of M . Due to the symmetries of M , all its eigenvalues have multiplicity 2, and for each eigenvalue λ, −λ is also an eigenvalue due to the antisymmetry of A. The Gaussian integrals can then be solved recursively by working through the list of distinct positive eigenvalues of M . One obtains Therefore, the variance of ω is given by The tail of the distribution is dominated by the largest eigenvalue λmax, and the asymptotic behaviour is p(ω) ∼ exp(−|ω|/2λmaxσ 2 lens ). With A jk given by eq. (B10) we find that λmax ∼ r 4 /n, which means that for ∆r = r/n fixed, the r.m.s. amplitude of ω should scale as ∼ r 3 .
In the real Universe matters are more complicated be-causeψ0 is poorly described by a Gaussian field. In particular, the distribution ofψ0 has a significant non-Gaussian tail of strong lenses. In this situation it may be justified to assume that extreme values of ω are dominated by those lines of sight which encounter two separate extreme values ofψ0. We can then use extreme value theory to make statements about the tails of the probability distribution of ω.
Let us be very simplistic here and neglect the coefficient matrix in eq. (B9) and pretend that instead ω is given by a sum of products of independent and identically distributed pairs of random variables, not necessarily Gaussian. If, as we said, for large ω this sum is dominated by a single term, generated from the product of the maximum value Xmax of all the Xj with the maximum value Ymax of all the Y k , the probability density of ω becomes a product distribution based on the distribution of those maximum values.
The Fisher-Tippett-Gnedenko extreme value theorem states that for a large number of draws Xj from a realistic and possibly non-Gaussian distribution, the distribution of Xmax asymptotically approaches the Gumbel distribution (and likewise for Ymax). Since this distribution is known, we can again estimate the distribution of ω from that extreme pair. We get with two unknown parameters s and µ that depend on the distribution of Xj. The asymptotic behaviour can be studied numerically and the parameters can be determined by fitting to observations. In practice we expect that the above expression is a good description only in the tails of p(ω). To find the best fit, we consider only the bins outside the 95% region centered on ω = 0. For such large values of |ω| the asymptotic distribution (B16) is dominated by the exponential behaviour of the Gumbel distribution for large Xmax, and µ is therefore nearly degenerate with the normalisation. This allows us to set µ 0 and fit only for s and the normalisation. We verified that treating µ as a free parameter does not improve the goodness of fit noticeably.

APPENDIX C: PIXELISATION AND RESOLUTION ERROR
The angular power spectra extracted from our simulations exhibit a suppression of power on small scales, reaching ∼ 50% at = 3000. We identify two sources of error introduced by our method that cause a small-scale power suppression: the pixelisation of the fields and the finite resolution of the simulations. In this section we quantify how much these two effects impact our results.
The pixelisation error is due to the fact that we estimate the value of a field (i.e. κ, or ω) in a pixel by aver- aging over the observed values at the particle positions for all particles whose observed position falls within that pixel. In each pixel we have a different number of particles and the mean of their positions does not correspond to the centre of our pixel. In order to quantify the error introduced by the pixelisation, we generate a high resolution Gaussian map for the convergence (N side = 8192) and we smooth this map into a lower resolution map in the following way: the value of the convergence in each pixel of the smoothed map is computed by averaging over a subset of the high-resolution map values within that pixel. Each value of the convergence within the pixel has a 30% probability of being accepted, which yields a sample rate similar to low-density regions in our source catalogue. We performed this analysis for three different smoothed maps, with N side = 512, 1024, 2048.
In Fig. C1 we show the angular power spectra computed for the three simulated smoothed maps and we compare the results with the Halofit prediction for the convergence spectrum. The pixelisation suppresses power on small scales as expected, and the effect depends on the number of pixels (N pixel = 12 N 2 side ). For N side = 2048, the power suppression is < 5% for 3000. Therefore, we adopt this angular resolution for the analysis presented in the main part of this paper.
Another cause of small-scale power suppression is the finite grid resolution of our simulation. In fact, gevolution is a particle-mesh code and therefore modes smaller than the Nyquist frequency cannot be resolved. The resolution of the grid in our simulation is ∆r = 312.5 kpc/h, which correspond to a Nyquist frequency kNyq ≈ 10 h/Mpc. This resolution allows us to resolve angular scales of around max ∼ r(z)kNyq ∼ 20000 at z = 1. However, the weak-lensing observables are integrated fields. Therefore, the larger resolution error from smaller redshifts propagates into the estimated convergence at all redshifts.
In order to quantify this effect, we model the power suppression in the power spectrum as where Psim is the power spectrum in our particle-mesh simulation, while Ptrue denotes the exact non-linear power spectrum in the continuum limit. A is a smooth function of redshift that can be determined empirically.
In order to quantify the impact of the suppression in eq. (C1) on the convergence and shear spectrum, we compute the convergence spectrum in the Limber approximation (Kaiser & Squires 1993;Scoccimarro et al. 1999;Lemos et al. 2017) where r(z) is the comoving distance at the source redshift (i.e. the mean redshift in the bin), D+(a) is the linear growth factor normalised to 1 at z = 0 and Psim is computed from eq. (C1) assuming the Halofit power spectrum to be the true power spectrum. In eq. (C2) we neglect the scaledependence of the growth for simplicity. However, this approximation does not significantly impact the result of our error analysis where our model for the power suppression in eq. (C1) is mainly valid at large scales.
In Fig. C2 we show the relative difference between the suppressed power spectrum in eq. (C2) and the angular power spectrum for an ideal simulation with kNyq → ∞. The factor A has been set to 1 in this plot.
The factor A in eq. (C1) can be estimated directly from the simulation power spectra. The power spectra computed from gevolution show a drop of power on small scales due to the finite resolution of the grid during the gravitational evolution, which effectively leads to a softening of the gravitational forces. Assuming that the "true" power spectrum is given by the Halofit recipe we can get an empirical fit of A due to modified gravitational evolution.
In addition, since the gravitational potential in the simulation is computed from a smoothed density field, i.e. the density field is estimated through cloud-in-cell (CIC) particle to mesh projection, there is another contribution to the coefficient A due to this projection. In the estimation of the matter power spectrum, the CIC kernel is deconvolved and therefore already accounted for. This correction obscures the fact that the CIC projection still has an effect on the gravitational potential. In other words, the coefficient A has two components, A = Aevo + ACIC, where Aevo describes the empirical suppression (due to force softening) of the matter power spectrum as explained above, and the coefficient ACIC can be estimated analytically from the CIC smoothing kernel as ACIC = π 2 /12.
Taking into account both the finite-resolution loss of power and the cloud-in-cell smoothing, we estimate the coefficient A to be in the range A ∈ [4.5, 6] for redshifts in the range z ∈ [1, 3].
In Fig. C3 we compare the relative difference between the convergence spectrum estimated by our method and the Halofit power spectrum and the model for the power-loss given by eq. (C1). Up to ≈ 2000, the simple model in eqs. (C1) and (C2) explains the small-scale power suppression in the weak-lensing observables obtained from our simulation.

APPENDIX D: VALIDATION OF THE LENSING POTENTIAL CALCULATION
In order to validate our method to compute κ lin directly in pixel space from the lensing potential (see Sec. 2.4), we apply our pipeline to the simulation presented in Hilbert et al. (2019) which they use to compare five different numerical approaches. We re-run the simulation in gevolution, starting from a particle snapshot at z 1.22 and using a mesh of 3072 3 grid points, which yields a spatial resolution of ∼ 166.6 kpc/h. This is not quite as good as the maximum resolution reached by the adaptive codes employed in the original comparison, but it is sufficient for our purpose. We then build the light cone for the Newtonian potential and construct the pixel map of κ lin for a source plane at a comoving distance of 2286 Mpc/h, covering the exact same 10 × 10 degree patch that was used for the code comparison in Hilbert et al. (2019).
After applying a Gaussian smoothing to the scale of one arcminute, we find an excellent agreement of the one-point PDF of κ lin between the different codes, shown in Fig. D1. We therefore conclude that our pixel-based numerical framework is consistent with the current state of the art.