Euclid: Reconstruction of Weak Lensing mass maps for non-Gaussianity studies

Weak lensing has proven to be an efficient method to constrain models of structure formation and reveal the nature of dark energy. So far, most weak lensing studies have focused on the shear field that can be measured directly from the ellipticity of background galaxies. However, within the context of forthcoming full-sky weak lensing surveys such as Euclid, convergence maps offer an important advantage over shear fields in terms of cosmological exploitation. While carrying the same information, the lensing signal is more compressed in the convergence maps than in the shear field, simplifying otherwise computationally expensive analyses, for instance non-Gaussianity studies. However, the inversion of the non-local shear field requires accurate control of systematic effects due to holes in the data field, field borders, noise and the fact that the shear is not a direct observable (reduced shear). In this paper, we present the two mass inversion methods that are being included in the official Euclid data processing pipeline: the standard Kaiser&Squires method (KS) and a new mass inversion method (KS+) that aims to reduce the information loss during the mass inversion. This new method is based on the KS methodology and includes corrections for mass mapping systematic effects. The results of the KS+ method are compared to the original implementation of the KS method in its simplest form, using the Euclid Flagship mock galaxy catalogue. In particular, we estimate the quality of the reconstruction by comparing the two-point correlation functions, third- and fourth-order moments obtained from shear and convergence maps, and we analyse each systematic effect independently and simultaneously. We show that the KS+ method reduces substantially the errors on the two-point correlation function and moments compared to the KS method. In particular, we show...


Introduction
Gravitational lensing is the process in which light from background galaxies is deflected as it travels towards us, as a result of the gravitation of the intervening mass. The measurement of the deformations in a large sample of galaxies offers a direct probe of the matter distribution in the Universe (including dark matter) and can thus be directly compared to theoretical models of structure formation. The statistical properties of the weak lensing field can be assessed by a statistical analysis either of the shear field or of the convergence field. On the one hand, convergence is a direct tracer of the total matter distribution integrated along the line-of-sight and it has therefore a direct link with the theory. On the other hand the shear (or more exactly the This paper is published on behalf of the Euclid Consortium Email: sandrine.pires@cea.fr reduced shear) is a direct observable and usually preferred for simplicity reasons.
Thus, the most common method to characterise the weak lensing field distribution is the shear two-point correlation function. It is followed very closely by the mass aperture two-point correlation functions, which are the result of the convolution of the shear two-point correlation functions by a compensated filter (Schneider et al. 2002) that is able to separate the E and B-modes of the two-point correlation functions (Crittenden et al. 2002). However, gravitational clustering is a nonlinear process and in particular, at small scales, the mass distribution is highly non-Gaussian. For this reason, several estimators of the three-point correlation functions have been proposed, either in the shear field (Bernardeau et al. 2002b;Benabed & Scoccimarro 2006) or using the mass aperture filter (Kilbinger & Schneider 2005). The three-point correlation functions are the lowest order statistics to quantify non-Gaussianity in the weak lensing field and thus provide additional information on structure formation models.
The convergence field can also be used to measure the two-and three-point correlation functions and other higherorder statistics. Assuming the mass inversion (namely the computation of the convergence map from the measured shear field) is properly conducted, there is the same information in the shear field as in the convergence maps (e.g. Schneider et al. 2002;Shi et al. 2011). While carrying the same information, the lensing signal is more compressed in the convergence maps than in the shear field, which makes it easier to extract and computationally less expensive. The convergence maps becomes a new tool that might bring additional constraints complementary to those that we can obtain from the shear field. However, the weak lensing signal being highly non-Gaussian at small scales, mass inversion methods using smoothing or denoising to regularise the problem are not optimal.
The reconstruction of convergence maps from weak lensing is a difficult task because of noise, irregular sampling, complex survey geometry, and the fact that the shear is not a direct observable. This is an ill-posed inverse problem and requires regularisation to avoid pollution from spurious Bmodes. Several methods have been derived to reconstruct the projected mass distribution from the observed shear field. The first non-parametric mass reconstruction was proposed by Kaiser & Squires (1993) and further improved by Bartelmann (1995), Kaiser (1995), Schneider (1995), and Squires & Kaiser (1996). These linear inversion methods are based on smoothing with a fixed kernel, which acts as a regularisation of the inverse problem. Nonlinear reconstruction methods were also proposed using different sets of priors and noise regularisation techniques (Bridle et al. 1998;Seitz et al. 1998;Marshall et al. 2002;Pires et al. 2009b;Jullo et al. 2014;Lanusse et al. 2016). Convergence mass maps have been built from many surveys including: the COSMOS Survey , the CFHT Lensing Survey CFHTLenS (Van Waerbeke et al. 2013), the CFHT/MegaCam Stripe-82 Survey (Shan et al. 2014), the Dark Energy Survey Science Verification DES SV Vikram et al. 2015;Jeffrey et al. 2018), the Red Cluster Sequence Lensing Survey RCSLenS (Hildebrandt et al. 2016), and the Hyper SuprimeCam Survey (Oguri et al. 2018). With the exception of Jeffrey et al. (2018) who use the nonlinear reconstruction proposed by Lanusse et al. (2016), all these methods are based on the standard Kaiser & Squires methodology.
In the near future, several wide and deep weak lensing surveys are planned: Euclid (Laureijs et al. 2011), LSST (LSST Science Collaboration et al. 2009), and WFIRST (Green et al. 2012). In particular, the Euclid satellite will survey 15 000 deg 2 of the sky to map the geometry of the dark Universe. One of the goals of the Euclid mission is to produce convergence maps for non-Gaussianity studies and constrain cosmological parameters. Therefore, two different mass inversion methods are being included into the official Euclid data processing pipeline. The first method is the standard Kaiser & Squires method (hereafter KS). Although it is well known that the KS method has several shortcomings it is taken as the reference for cross-checking of the results. The second method is a new nonlinear mass inversion method (hereafter KS+) based on the formalism developed in Pires et al. (2009b). The KS+ method aims at performing the mass inversion with minimal information loss. This is done by performing the mass inversion with no other regularisation than binning while controlling systematic effects.
In this paper, the performance of these two mass inversion methods is investigated using the Euclid Flagship mock galaxy catalogue (version 1.3.3,Castander F. et al. in prep) with realistic observational effects (i.e. noise, missing data, the reduced shear observable). The effect of intrinsic alignments is not studied in this paper because of the lack of simulations modelling intrinsic alignments properly. However, intrinsic alignments also need to be considered seriously because they affect second-and higher-order statistics. A contribution of several percent is expected to two-point statistics (see e.g. Joachimi et al. 2013).
In this study, we compare the results obtained with the KS+ method to those obtained with a version of the KS method in which no smoothing step is performed, other than binning. We quantify the quality of the reconstruction using two-point correlation functions and moments of the convergence. Our tests illustrate the efficacy of the different mass inversion methods in preserving the second-order statistics and higher-order moments.
The paper is organised as follows. In Sect. 2, we present the weak lensing mass inversion problem and the standard KS method. Section 3 presents the KS+ method used in this paper to deal with the different systematic effects. In Sect. 4, we explain the methodology to compare these two mass inversion methods. In Sect. 5, we use the Euclid Flagship mock galaxy catalogue with realistic observational effects such as noise, complex survey geometry and considering the reduced shear to investigate the performance of the two mass inversion methods. First, we derive simulations including only one issue at a time to test each systematic effect independently. Then we derive realistic simulations including them all and study the systematic effects simultaneously. We conclude in Sect. 6.

The weak gravitational lensing formalism
In weak lensing surveys, the shear field γ(θ) is derived from the ellipticities of the background galaxies at position θ in the image. The two components of the shear can be written in terms of the lensing potential ψ(θ) as (see e.g. Bartelmann & Schneider 2001) where the partial derivatives ∂ i are with respect to the angular coordinates θ i , i = 1, 2 representing the two dimensions of sky coordinates. The convergence κ(θ) can also be expressed in terms of the lensing potential as For large-scale structure lensing, assuming a spatially flat Universe, the convergence at a sky position θ from sources at comoving distance r is defined by where H 0 is the Hubble constant, Ω m is the matter density, a the scale factor and δ ≡ (ρ −ρ)/ρ the density contrast (where ρ andρ are the 3D density and the mean 3D density respectively). In practice, the expression for κ can be generalised to sources with a distribution in redshift or equivalently in comoving distance f (r), yielding where r H is the comoving distance to the horizon. The convergence map reconstructed over a region on the sky gives us the integrated mass density fluctuation weighted by the lensing weight function p(r ),

The KS mass inversion problem
The weak lensing mass inversion problem consists of reconstructing the convergence κ from the measured shear field γ. We can use complex notation to represent the shear field, γ = γ 1 + iγ 2 , and the convergence field, κ = κ E + iκ B , with κ E corresponding to the curl-free component and κ B to the gradient-free component of the field, called E-and B-modes by analogy with the electromagnetic field. Then, from Eq. (1) and Eq. (2), we can derive the relation between the shear field γ and the convergence field κ. For this purpose, we take the Fourier transform of these equations and obtain where the hat symbol denotes Fourier transforms,P =P 1 + iP 2 , with 2 ≡ 2 1 + 2 2 and i the wave numbers corresponding to the angular coordinates θ i . P is a unitary operator. The inverse operator is its complex conjuguateP * =P 1 −iP 2 as shown by Kaiser & Squires (1993), Note that to recover κ from γ there is a degeneracy when 1 = 2 = 0. Therefore, the mean value of κ cannot be recovered if only shear information is available. This is the socalled mass-sheet degeneracy (see e.g. Bartelmann 1995, for a discussion). In practice, we impose the mean convergence to vanish across the survey by setting the reconstructed = 0 mode to zero. This is a reasonable assumption for large field reconstruction (e.g. Seljak 1998).
We can easily derive an estimator of the E-mode and B-mode convergence in the Fourier domaiñ The weak lensing arising from a scalar potential (the lensing potential ψ), it can be shown that weak lensing only produces E-modes. However, intrinsic alignments, and imperfect corrections of the point spread function (PSF) generally generate both E and B modes. The presence of B-modes can thus be used to test for residual systematic effects in current weak lensing surveys.

The weak lensing missing data problem
The shear is only sampled at the discrete positions of the galaxies where the ellipticity is measured. Therefore, the first step of the mass map inversion method is to bin the observed ellipticities of galaxies on a regular pixel grid to create what we refer to as the observed shear maps γ obs . Some regions remain empty because of various masks applied to the data such as the masking out of bright stars or camera CCD defects. In such cases, the shear is set to zero in the original KS method, with M the binary mask (i.e. M = 1 if we have information at the pixel, M = 0 otherwise) and γ n the noisy shear maps. As the shear at any sky position is non-zero in general, this introduces errors in the reconstructed convergence maps. Some specific methods address this problem by discarding masked pixels at the noise regularisation step (e.g. Van Waerbeke et al. 2013). However, as explained previously, this intrinsic filtering results in subtantial signal loss at small scales. Instead, inpainting techniques are used in the KS+ method to fill the masked regions (see Appendix A).

The weak lensing shape noise
The gravitational shear is derived from the ellipticities of the background galaxies. However, the galaxies are not intrinsically circular, so their measured ellipticity is a combination of their intrinsic ellipticity and the gravitational lensing shear. The shear is also subject to measurement noise and uncertainties in the PSF correction. All these effects can be modelled as an additive noise N = N 1 + iN 2 , The noise terms N 1 and N 2 are assumed to be Gaussian and uncorrelated with zero mean and standard deviation, where N i g is the number of galaxies in pixel i. The root mean square shear dispersion per galaxy, σ , arises both from the measurement uncertainties and the intrinsic shape dispersion of galaxies. The Gaussian assumption is a reasonable assumption and σ is set to 0.3 for each component as is generally found in weak lensing surveys (e.g. Leauthaud et al. 2007;Schrabback et al. 2015Schrabback et al. , 2018. The surface density of usable galaxies is expected to be around n g = 30 gal. arcmin −2 for the Euclid Wide survey (Cropper et al. 2013).
The derived convergence map is also subject to an additive noise, κ n =P * γn =κ +P * N .
In particular, the E-component of the convergence noise is where * denotes the convolution operator, P 1 and P 2 are the inverse Fourier transforms ofP 1 andP 2 . Considering the case where the shear noise terms N 1 and N 2 are Gaussian, uncorrelated and with a standard deviation constant across the field, the convergence noise is also Gaussian and uncorrelated. In practice, the number of galaxies varies slightly across the field. Also, the variances of N 1 and N 2 might be slightly different, which can be modelled by different values of σ for each component. These effects introduce noise correlations in the convergence noise maps, but they were found to remain negligible compared to other effects studied in this paper.
In the KS method a smoothing by a Gaussian filter is frequently applied to the background ellipticities before performing the mass inversion to regularise the solution. Although performed in most applications of the KS method, this noise regularisation step is not mandatory. It was introduced to avoid infinite noise and divergence at very small scales. However, the pixelisation already provide an intrinsic regularisation. Then, there is no need for an additional noise regularisation prior to the inversion. Nonetheless, for specific applications that require denoising in any case, the filtering step can be performed before or after the mass inversion.

Improved Kaiser & Squires method (KS+)
Systematic effects in mass inversion techniques must be fully controlled in order to use convergence maps as cosmological probes for future wide field weak lensing experiment such as Euclid . Here, we introduce the KS+ method based on the formalism developed in Pires et al. (2009b) and Jullo et al. (2014) that integrates the necessary corrections to deal with imperfect and realistic measurements. We summarise its improvements over KS in this section and evaluate its performance in Sect. 5.
In this paper, the mass mapping formalism is developed in the plane. The mass inversion can also be performed on the sphere as proposed in Pichon et al. (2010) and Chang et al. (2018), and the extension of the KS+ method to the curved sky is being investigated. However, the computation time and memory required to process the spherical mass inversion put limitations in terms of convergence maps resolution and/or complexity of the algorithm. Thus, planar mass inversions remain important to reconstruct convergence maps with a good resolution and probe the non-Gaussian features of the weak lensing field (e.g. for peak count studies).

Dealing with missing data
Considering the weak lensing shear field γ sampled on a grid of N × N pixels, we can describe the complex shear and convergence fields by their respective matrices. In the rest of the paper, the notations γ and κ stand for the matrix quantities.
In the standard version of the KS method, the pixels with no galaxies are set to zero. Fig. 1 shows an example of simulated shear maps without noise derived from the Euclid Flagship mock galaxy catalogue (see Sect. 4.3 for more details). The upper panels of Fig. 1 show the two components of the shear with zero values (displayed in black) corresponding to the mask of the missing data. These zero values generate an important leakage during the mass inversion.
With KS+, the problem is reformulated including additional assumptions to regularise the problem. The convergence κ can be analysed using a transformation Φ yielding to a set of coefficients α = Φ T κ (Φ is an orthogonal matrix operator and Φ T represents the transpose matrix of Φ). In the case of the Fourier transformation, Φ T would correspond to the discrete Fourier transform (DFT) matrix and α would be the Fourier coefficients of κ. The KS+ method uses a prior of sparsity, i.e. it assumes that there is a transformation Φ where the convergence κ can be decomposed into a set of coefficients α, where most of its coefficients are close to zero. In this paper, Φ is chosen to be the discrete cosine transform (DCT) following Pires et al. (2009b). The DCT expresses a signal in terms of a sum of cosine functions with different frequencies and amplitudes. It is similar to the DFT but using smoother boundary conditions, providing a sparser representation. Hence the use of the DCT for JPEG compression.
We can rewrite the relation between the observed shear γ obs and the noisy convergence κ n as with M being the mask operator and P the KS mass inversion operator. There is an infinite number of convergence κ n that can fit the observed shear γ obs . With KS+, first, we impose the mean convergence to vanish across the survey, as in the KS method. Then, among all possible solutions, KS+ searches for the sparsest solutionκ n in the DCT Φ (i.e. the convergence κ n that can be represented with the fewest large coefficients). The solution of this mass inversion problem is obtained by solving where ||z|| 0 the pseudo-norm, i.e. the number of nonzero entries in z, ||z|| the classical l 2 norm (i.e. ||z|| = k (z k ) 2 ), and σ stands for the standard deviation of the input shear map measured outside the mask. The solution of such an optimisation task can be obtained through an iterative thresholding algorithm called morphological component analysis (MCA) introduced by Elad et al. (2005) and adapted to the weak lensing problem in Pires et al. (2009b).
In Pires et al. (2009b), an additional constraint is used to force the B-modes to zero. This is optimal when the shear maps have no B-modes. However, any real observation has some residual B-modes due to intrinsic alignments, imperfect PSF correction, etc. The B-mode power is then transferred to the E-modes which degrades the E-mode convergence reconstruction. Here instead, the B-modes are let free and an additional constraint is set on the power spectrum of the convergence map. To this end, we use a wavelet transform to decompose the convergence maps into a set of aperture mass maps using the starlet transform algorithm (Starck et al. 1998;Starck & Murtagh 2002). Then, the constraint consists of renormalizing the standard deviation (or equivalently, the variance) of each aperture mass map inside the mask regions to the values measured in the data, outside the masks, and then reconstructing the convergence via where the pixels with no galaxies are set to zero (displayed in black). The lower panels show the result of the inpainting method that allows us to fill the gaps judiciously.
the inverse wavelet transform. The variance per scale corresponding to the power spectrum at these scales allows us to constrain a broadband power spectrum of the convergence κ inside the gaps.
Adding the power spectrum constraints yields the final sparse optimisation problem where W is the forward wavelet transform and W T its inverse transform and Q is the linear operator used to impose the power spectrum constraint. More details about the KS+ algorithm are given in Appendix A. The KS+ method allows us to reconstruct the inpainted convergence maps and the corresponding inpainted shear maps, where the empty pixels are replaced by non-zero values. These interpolated values preserve the continuity of the signal and reduce the signal leakage during the mass inversion (see lower panels of Fig. 1). The quality of the convergence maps reconstruction with respect to missing data is evaluated in Sect. 5. Additionally, the new constraint allows us to use the residual B-modes of the reconstructed maps to test for the presence of residual systematic effects and possibly validate the shear measurement processing chain.

Dealing with field border effects
The KS and KS+ mass inversion methods relate the convergence and the shear fields in Fourier space. However, the discrete Fourier transform implicitly assumes that the image is periodic along both dimensions. Since there is no reason for opposite borders to be alike, the periodic image generally presents strong discontinuities across the frame border. These discontinuities cause several artefacts at the borders of the reconstructed convergence maps. The field border effects can be addressed by removing the image borders, which throws away a large fraction of the data. Direct finite-field mass inversion methods have also been proposed (e.g. Seitz & Schneider 1996. Although unbiased, convergence maps reconstructed using these methods are noisier than the ones obtained with the KS method. In the KS+ method, the problem of borders is solved by taking a support for the image two times larger and by considering the borders as masked regions to be inpainted. The upper panels of Fig. 2 show the two components of a shear map covering 5 • ×5 • degrees and extended to a field of 10 • ×10 • . The inpainting method is then used to recover the shear at the field boundaries as shown in the lower panels of Fig. 2. Once the mass inversion is performed, the additional borders are removed. This technique reduces the field border effects by pushing the border discontinuities farther away.

Dealing with reduced shear
In Sect. 2.2, we assumed knowledge of the shear, in which case the mass inversion is linear. In practice, the observed galaxy ellipticity is not induced by the shear γ but by the reduced shear g that depends on the convergence κ corresponding to that particular line-of-sight, While the difference between the shear γ and the reduced shear g is small in the regime of cosmic shear (κ 1), its neglect might nevertheless cause a measurable bias at small angular scales (see e.g. White 2005;Shapiro 2009). In the standard version of KS, the Fourier estimators are only valid when the convergence is small (κ 1) and no longer hold near the centre of massive galaxy clusters. The mass inversion problem becomes nonlinear and it is therefore important to properly account for reduced shear.
In the KS+ method, an iterative scheme is used to recover the E-mode convergence map as proposed in Seitz & Schneider (1995). The method consists of solving the linear inverse problem iteratively (see Eq. 9), using at each iteration the previous estimate of the E-mode convergence to correct the reduced shear using Eq. (18). Each iteration then provides a better estimate of the shear. This iterative algorithm is found in Jullo et al. (2014) to converge quickly to the solution (about 3 iterations). The KS+ method uses the same iterative scheme to correct for reduced shear and we find that it is a reasonable assumption in the case of large-scale structure lensing.

Dealing with noise
In the original implementation of KS, the shear maps are first regularised with a smoothing window (i.e. a low-pass filter) to obtain a smoothed version of the shear field. Then, Eq. (9) is applied to derive the convergence maps. In contrast, the KS+ method aims at producing very general convergence maps for many applications. In particular, it produces noisy maps with minimal information loss.
However, for specific applications (e.g galaxy cluster detection and characterisation), it can be useful to add an extra denoising step, using any of the many regularisation techniques that have been proposed (Bridle et al. 1998;Seitz et al. 1998;Marshall et al. 2002;Starck et al. 2006;Lanusse et al. 2016). To compare the results of the KS and KS+ methods on noisy maps, we use a linear Gaussian and the nonlinear MRLens filter (Starck et al. 2006) for noise suppression. Fig. 3 illustrates the effect of noise in the convergence map reconstruction. The upper panels show one E- mode convergence map reconstructed from noise-free (left) and noisy (right) shear data. The convergence map is dominated by the noise. Lower panels show the results of the Gaussian filter (left) and MRLens filter (right). The Gaussian filter gives a smoothed version of the noisy convergence map whose level of smoothness is set by the width of the Gaussian (σ). Thus, the amplitude of the over-densities (in blue) are systematically lowered by the Gaussian filter. In contrast, the MRLens filter uses a prior of sparsity to better recover the amplitude of the structures and uses a parameter namely the false discovery rate (α FDR ) to control the average fraction of false detections (i.e. the number of pixels truly inactive declared positive) made over the total number of detections (Benjamini & Hochberg 1995).
For some other applications (e.g. two or three-point correlation), the integrity of the reconstructed noisy convergence maps might be essential and this denoising step can be avoided.

Comparing second-order statistics
The most common tools for constraining cosmological parameters in weak lensing studies are the shear two-point correlation functions. Following Bartelmann & Schneider (2001), they are defined by considering pairs of positions ϑ and θ+ϑ, and defining the tangential and cross-component of the shear γ t and γ × at position ϑ for this pair as where ϕ is the polar angle of the separation vector θ. Then we define the two independent shear correlation functions where the Bessel function J 0 (J 4 ) corresponds to the "+" ("-") correlation function, P κ ( ) is the power spectrum of the projected matter density, and the Fourier variable on the sky. We can also compute the two-point correlation functions of the convergence (κ = κ E + iκ B ) defined as We can verify that these two quantities agree (Schneider et al. 2002): Article number, page 7 of 17 A&A proofs: manuscript no. MassMapping_Euclid_AA_v5 When the B-modes in the shear field are consistent with zero, the two-point correlation of the shear (ξ + ) is equal to the two-point correlation of the convergence (ξ κE ). Then the differences between the two are due to the errors introduced by the mass inversion to go from shear to convergence. We compute these two-point correlation functions using the tree code athena (Kilbinger et al. 2014). The shear twopoint correlation functions are computed by averaging over pairs of galaxies of the mock galaxy catalogue, whereas the convergence two-point correlation functions are computed by averaging over pairs of pixels in the convergence map. The convergence two-point correlation functions can only be computed for separation vectors θ allowed by the binning of the convergence map.

Comparing higher-order statistics
Two-point statistics cannot fully characterise the weak lensing field at small scales where it becomes non-Gaussian (e.g. Bernardeau et al. 2002a). Since the small-scale features carry important cosmological information, we compute the third-order moment, κ 3 E , and the fourth-order moment, κ 4 E , of the convergence. Computations are performed on the original convergence maps provided by the Flagship simulation, as well as on the convergence maps reconstructed from the shear field with the KS and KS+ methods. We evaluate the moments of convergence at various scales by computing aperture mass maps (Schneider 1996;Schneider et al. 1998). Aperture mass maps are usually obtained by convolving the convergence maps with a filter function of a specific scale (i.e. aperture radii). In this study, this is done by means of a wavelet transform using the starlet transform algorithm (Starck et al. 1998;Starck & Murtagh 2002) that simultaneously produces a set of aperture mass maps on dyadic (powers of two) scales (see Appendix A for more details). In Leonard et al. (2012), it is demonstrated that the aperture mass is formally identical to a wavelet transform at a specific scale and the aperture mass filter corresponding to this transform is derived. The wavelet transform offers significant advantages over the usual aperture mass algorithm in terms of computation time, providing speed-up factors of about 5 to 1200 depending on the scale.

Numerical simulations
For this study, we use the Euclid Flagship mock galaxy catalogue version 1. 3.3 (Castander F. et al.,in prep) derived from N-body cosmological simulation (Potter et al. 2017) with parameters Ω m = 0.319, Ω b = 0.049, Ω Λ = 0.681, σ 8 = 0.83, n s = 0.96, h = 0.67, and the particle mass is m p ∼ 2.398 × 10 9 M h −1 . The galaxy light-cone catalogue contains 2.6 billion galaxies over 5000 deg 2 and it extends up to z = 2.3. It has been built using a hybrid halo occupation distribution and halo abundance matching (HOD+HAM) technique, whose galaxy clustering properties were discussed in detail in Crocce et al. (2015). The lensing properties have been computed using the Born approximation and projected mass density maps (in HEALPix format with N side = 8192) generated from the particle lightcone of the Flagship simulation. More details on the lensing properties of the Flagship mock galaxy catalogue can be found in Fosalba et al. (2015Fosalba et al. ( , 2008.
In order to evaluate the errors introduced by the mass mapping methods, we extract 10 contiguous shear/convergence fields of 10 • × 10 • from the Flagship mock galaxy catalogue, yielding a total area of 1000 deg 2 . The fields correspond to galaxies that lie in the range of 15 • < α < 75 • and 15 • < δ < 35 • where α and δ are the right ascension and declination, respectively. In order to get the density of 30 galaxies per arcmin 2 foreseen for the Euclid Wide survey, we select randomly one quarter of all galaxies in the catalogue. Then projected shear and convergence maps are constructed from the combination of all the redshifts of the selected galaxies. More sophisticated selection methods based on galaxy magnitude would produce slightly different maps. However they would not change the performances of the two methods studied here. The fields are down-sampled to 1024 × 1024 pixels, which corresponds to a pixel size of about 0. 6. Throughout all the paper, the shaded regions stand for the uncertainties on the mean estimated from the total 1000 deg 2 of the 10 fields.

Shear field projection
In this study, we consider fields of 10 • × 10 • . The fields are taken sufficiently small to be approximated by a tangent plane. We use a gnomonic projection to project the points of the celestial sphere onto a tangent plane, following Pires et al. (2012b), who found that it preserves the two-point statistics. We note, however, that higher-order statistics may behave differently under different projections.
The shear field projection is obtained by projecting the galaxy positions from the sphere (α, δ) in the catalogue onto a tangent plane (x, y). The projection of a non-zero spin field such as the shear field requires a projection of both the galaxy positions and their orientations. Projections of the shear do not preserve the spin orientation which can generate substantial B-modes (depending on the declination), if not corrected. There are two problems to consider due to the orientation. First, the projection of the meridians are not parallel, so that the North direction is not the same everywhere in the same projected field of view. Second, the projection of the meridians and great circles are not perpendicular, so that the system is locally non-Cartesian. Since we are dealing properly with the other effects (e.g. noise, missing data or border effects) and considering large fields of view (10 • × 10 • ) possibly at high latitudes, these effects need to be considered. The first effect is dominant and generates substantial B-modes (increasing with latitude) if not corrected. It can be easily corrected by measuring the shear orientation with respect to the local North direction. We find that this correction is sufficient for the residual errors due to projection to become negligible compared to errors due to other effects.

Impact of the systematic effects on the mass map inversion
In this section, we quantify the impact of field borders, missing data, noise and the approximation of shear by reduced shear on the KS and KS+ mass inversion methods.

Dealing with missing data effects
We use the 10 noise-free shear fields of 10 • × 10 • described in Sect. 4.3 and the corresponding noise-free convergence maps. We convert the shear fields into planar convergence maps using the KS and KS+ methods, masking 20% of the data as expected for the Euclid survey. The mask is derived from the Data Challenge 2 catalogues produced by the Euclid collaboration using the code FLASK (Xavier et al. 2016). Fig. 4 compares the results of the KS and KS+ methods when dealing with missing data. It shows the residual maps i.e the pixel difference between the original E-mode convergence map and the reconstructed ones. The amplitude of the residuals is larger with the KS method. Detailed investigation shows that the excess error is essentially localised around the gaps. As the mass inversion operator P is intrinsically non-local, it generates artefacts around the gaps. In order to quantify the average errors, Fig. 5 shows the probability distribution function (PDF) of the residual maps, estimated outside the mask. The standard deviation is 0.0080 with KS and 0.0062 with KS+. The residual errors obtained with KS are then 30% larger than the ones obtained with KS+. and KS+ (in red) convergence maps at the same scales. The KS and KS+ convergence maps are reconstructed from incomplete noise-free shear maps. The estimation of the third-and fourth-order moments is done outside the mask. The shaded area represents the uncertainties on the mean estimated on 1000 deg 2 .
we expect these two quantities to be equal within the precision of the simulations (see Sect. 4.1). The KS method yields a systematic underestimation of the original twopoint correlation function by a factor of about 2 on arcminute scales, but can reach factors of 5 at larger scales. The mass inversion operator P being unitary, the signal energy is conserved by the transformation (i.e. (γ 2 1 + γ 2 2 ) = (κ 2 E +κ 2 B ), where the summation is performed over all the pixels of the maps). We found that about 10% of the total energy leaks into the gaps and about 15% into the B-mode component. In contrast, the errors of the KS+ method are of the order of a few percent at scales smaller than 1 • . At any scale, the KS+ errors are about 5-10 times smaller than the KS errors, remaining in the 1σ uncertainty of the original two-point correlation function. Fig. 7 shows the third-order (upper panel) and fourthorder (lower panel) moments estimated at 6 different wavelet scales ( 2. 34, 4. 68, 9. 37, 18. 75, 37. 5, 75. 0) using the KS (in blue) and KS+ (in red) methods. To that purpose, the pixels inside the mask were set to zero in the reconstructed convergence maps. The aperture mass maps corresponding to each wavelet scale were computed, and the moments calculation was performed outside the masks.
The KS method systematically underestimates the third-and fourth-order moments at all scales. Below 10 , the errors on the moments remain smaller than 50% and they increase with scale up to a factor 3. In comparison, the KS+ errors remain much smaller at all scales, and remain within the 1σ uncertainty. Fig. 8 compares the results of the KS (left) and KS+ (right) methods when dealing with border effects. It shows the residual error maps corresponding to the pixel difference between the original E-mode convergence map and the reconstructed ones. With KS, as expected, the pixel difference shows errors at the border of the field. With KS+, there is also some low-level boundary effects but these errors are considerably reduced and do not show any significant structure at the field border. Again, the PDF of these residuals can be compared to quantify the errors. For the two methods, Fig. 9 shows the residuals PDFs computed at the boundaries (as dotted lines) and in the remaining central part of the image (as solid lines). The border width used to compute the residual PDF is 100 pixels corresponding to about 1 degree. With the KS method, the standard deviation of the residuals in the centre of the field is 0.0062. In the outer regions, the border effect causes errors of 0.0076 (i.e. 25% larger than at the centre). Away from the borders, the KS+ method gives results similar to the KS method (0.0060). However, it performs much better at the border where the error only reaches 0.0061. The small and uniform residuals of the KS+ method show how efficient it is in dealing with borders effects.

Dealing with field border effects
As before, the scale-dependence of the errors can be estimated using the two-point correlation function and higherorder moments computed at different scales. Fig. 10 shows the two-point correlation functions. For both methods, the errors increase with angular scale since the fraction of pairs of pixels including boundaries increase with scale. The loss of amplitude at the image border is responsible for significant errors in the two-point correlation function of the KS convergence maps. In contrast, the errors are about 5 to 10 times smaller with the KS+ method and remain in the 1σ uncertainty range of the original two-point correlation function. Fig. 11 shows the impact of field borders effects on the third-order (upper panel) and fourth-order (lower panel) moments of the convergence maps at different scales. As was observed earlier on the two-point correlation estimation, the KS method introduces errors at large scales on the third-and fourth-order moment estimation. With KS+, the discrepancy is about 1% and within the 1σ uncertainty.

Dealing with reduced shear
In this section, we quantify the errors due to the approximation of shear (γ) by the reduced shear (g). To this end, we used the noise-free shear fields described in Sect. 4.3 and we computed the reduced shear fields using Eq. (18) and the convergence provided by the catalogue. We then derived the reconstructed convergence maps using the KS and KS+ methods.
For both methods, the errors on the convergence maps are dominated by field border effects. We did not find any estimator able to separate these two effects and then identify the reduced shear effect in the convergence maps. But the errors introduced by the reduced shear can be assessed : Field border effects: PDF of the residual errors between the original E-mode convergence map and the convergence maps reconstructed using KS (in blue) and KS+ (in red). The dotted lines correspond to the PDF of the residual errors measured at the boundaries of the field and the solid lines to the PDF of the residual errors measured in the centre of the field. The borders are 100 pixels wide. by comparing the shear and reduced shear two-point correlation functions (see Fig. 12).
While the differences are negligible at large scale, they reach the percent level on arcminute scales (in agreement with White 2005), where they become comparable or larger than the KS+ errors due to border effects.

Dealing with noise
In this section, we study the effect of the noise on convergence maps. We derive noisy shear maps, assuming a Gaussian noise (σ = 0.3). Then, we compare the two mass inversion methods. The pixel difference cannot be used in this case, because the convergence maps are noise dominated (see Fig. 3, upper right panel). However, we can still assess the quality of the convergence maps using two-point    Relative two-point correlation error between the mean two-point correlation functions ξ γ + estimated from the shear fields and corresponding mean two-point correlation function ξ g + estimated from the reduced shear fields without any correction.
Moments of noisy maps are biased and potentially dominated by the noise contribution. For instance, the total variance in the noisy convergence map is expected to be the sum of the variance in the noise-free convergence map and the noise variance. Therefore moments of the noisy KS and KS+ convergence maps cannot be directly compared to moments of the original noise-free convergence maps. Instead, Fig. 14 compares them to the moments of the original convergence maps where noise was added with properties similar to the noise expected on the convergence maps. For this purpose, for each field, we generate noise maps N 1 and N 2 using Eq. (11) and (12) and we derive the noise to be added in the convergence using Eq. (14). Fig. 14 compared to Fig. 11 shows that the third-order moment of the convergence is unaffected by noise. In contrast, the fourth-order moment is biased for scales smaller than 10 . The two methods slightly underestimate the thirdand fourth-order moments at large scale. However, with KS+, the errors are reduced by a factor of 2 and remain roughly within the 1σ uncertainty.

Dealing with all systematic effects simultaneously
In this section, we assess the performance of KS and KS+ for realistic data sets, combining the effects of noise, reduced shear, borders and missing data. Fig. 15 compares the results of the KS method and the KS+ method combined with a filtering step in dealing with all systematic effects in one field. We use a Gaussian filter to reduce the noise in the KS convergence maps, but employ the nonlinear MRLens filter in the case of the KS+ maps. Again, KS+ is better in recovering the over-densities mainly for two reasons. First, KS+ reduces the signal leakage during the mass inversion compared to KS. Second, the nonlinear MRLens filter is particularly suited for detection of isotropic structures (Pires et al. 2009a(Pires et al. , 2012aLin et al. 2016) while the Gaussian filter reduces systematically the amplitude of the over-densities. Fig. 16 shows the two-point correlation computed with the two methods. The masked regions are excluded from the two-point correlation computation, resulting in fewer pairs, and larger noise relative to Fig. 13. Again, the strong leakage due to missing data is clearly observed with the KS method. The results obtained with the KS+ method reduce the errors in the mean convergence two-point correlation function by a factor of about 5 and the errors remain roughly within the 1σ uncertainty.
In Fig. 17, we test the efficacy of the mass inversion methods in preserving higher-order moments of the convergence maps in a realistic setting. As before, a realistic noise has been added to the original convergence maps for comparison. As was observed earlier in the noise-free case, the KS method systematically underestimates the thirdand fourth-order moments at all scales. With KS+, the errors are significantly reduced, by a factor of about 2 on the third-order moment and 10 on the fourth-order moment estimation, at all scales. Although reduced, the errors of the KS+ method on the third-order moment cannot be neglected. These errors might result from noise correlations introduced by the inpainting method in the shear maps. Indeed, inside the gaps, the noise is correlated because it is interpolated from the remaining data. These noise correlations propagate into the convergence maps and can explain the bias in the moments estimation. We note that the two-point correlation functions and higher-order moments are here only used to probe the accuracy of the reconstruction methods. For specific applications, the small residuals of the KS+ method can be reduced even more by using additional treatment such as, for instance, down-weighting the region around the mask when computing the moments (e.g. Van Waerbeke et al. 2013).

Conclusion
This paper is motivated by the use of convergence maps in Euclid to constrain cosmological parameters and to assess other physical constraints. Convergence maps encode the lensing information in a different manner, allowing more optimised computations than shear. However, the mass inversion process is subject to field border effects, missing data, reduced shear, intrinsic alignments, and noise. This requires accurate control of the systematic effects during the mass inversion to reduce as much as possible the information loss. This paper presents and compares the two mass inversion methods being included in the official Euclid data processing pipeline: the standard Kaiser & Squires (KS) method and an improved Kaiser & Squires (KS+) mass inversion technique integrating corrections for the mass mapping systematic effects. The impact of the systematic effects on the reconstructed convergence maps have been studied using the Euclid Flagship mock galaxy catalogue.
In a first step, we have analysed and quantified one by one the impact of systematic effects on reconstructed convergence maps using two-point correlation functions and moments of the convergence. In this manner, we quantify the contribution of each effect to the error budget and better understand the distribution of the errors in the convergence maps. With KS, missing data is the dominant effect at all scales. Field border effects also have an important impact but only affect the border of the maps. These two effects are significantly reduced with KS+. The reduced shear is the smallest effect in terms of contribution and only affects small angular scales. The study also shows that pixellisation provides an intrinsic regularisation and that no additional smoothing step is required to avoid infinite noise in the convergence maps.
In a second step, we have quantified the errors introduced by the KS and KS+ methods in a realistic setting including the aforementioned systematic effects. We have shown that the KS+ method reduces the errors on the twopoint correlation functions and on the moments of the convergence compared to the KS method. The errors introduced by the mass inversion on the two-point correlation of the convergence maps are reduced by a factor of about 5. The errors on the third-order and fourth-order moment estimates are reduced by factors of about 2 and 10 respectively. Some errors remain on the third-order moment that remain within the 2σ uncertainty. They might result from noise correlations introduced by the inpainting method inside the gaps.
Recent studies have shown that combining the shear two-point statistics with higher-order statistics of the convergence such as higher-order moments (Vicinanza et al. 2018), Minkowski functionals (Vicinanza et al. 2019) or peak counts (Liu et al. 2015;Martinet et al. 2018) allows to break common degeneracies. The precision of the KS+ mass inversion makes the E-mode convergence maps a promising tool for such cosmological studies. In future work, it is planned to propagate these errors into cosmological parameter constraints using higher-order moments and peak counts. The upper panels show the original E-mode convergence κ map (left) and the mask that is applied to the shear maps (right). Lower panels show the convergence map reconstructed from an incomplete noisy shear field using the KS method with a linear Gaussian filter of size σ = 3 (left) and using the KS+ method and the nonlinear MRLens filtering with α FDR = 0.05 (right). The field is 10 • × 10 • downsampled to 1024 × 1024 pixels.
of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the Centre National d'Etudes Spatiales, the Deutsches Zentrum für Luft-and Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciênca e a Tecnologia, the Ministerio de Economia y Competitividad, the National Aeronautics and Space Administration, the Netherlandse Onderzoekschool Voor Astronomie, the Norvegian Space Center, the Romanian Space Agency, the State Secretariat for Education, Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. Relative error Fig. 16: All systematic effects: Mean shear two-point correlation function ξ + (in black) and corresponding mean convergence two-point correlation function ξ κE estimated from incomplete noisy shear fields. The convergence maps have been estimated using KS (in blue) and KS+ (in red). The estimation of the convergence two-point correlations is done outside the mask. The shaded area represents the uncertainties on the mean estimated on 1000 deg 2 . The lower panel shows the normalised difference between the two upper curves. This provides a way to estimate a broadband power spectrum of the convergence κ from incomplete data. The power spectrum is then enforced by multiplying each wavelet coefficients by the factor σ out j /σ in j inside the gaps, where σ in j and σ out j are respectively the standard deviation estimated in the wavelet band w j inside and outside the mask. This normalisation can be described by a linear operator Q as used in Eq. (17). The contraint is applied both on the Emode and B-mode components before reconstructing the convergence κ by backward wavelet transform.