Toward an Optimal Reconstruction of the Shear Field with PDF-folding

Weak lensing provides a direct way of mapping the density distribution in the Universe. To reconstruct the density field from the shear catalog, an important step is to build the shear field from the shear catalog, which can be quite nontrivial due to the inhomogeneity of the background galaxy distribution and the shape noise. We propose the PDF-folding method as a statistically optimal way of reconstructing the shear field. It is an extention of the PDF-SYM method, which was previously designed for optimizing the stacked shear signal as well as the shear-shear correlation for the Fourier_Quad shear estimators. PDF-folding does not require smoothing kernels as in traditional methods, therefore it suffers less information loss on small scales and avoids possible biases due to the spatial variation in the shear on the scale of the kernel. We show with analytic reasoning as well as numerical examples that the new method can reach the optimal signal-to-noise ratio on the reconstructed shear map under general observing conditions, i.e., with inhomogeneous background densities or masks. We also show the performance of the new method on real data around foreground galaxy clusters.

To our knowledge, lensing is the only direct way of mapping out the matter (surface) density field so far, making it a particularly valuable tool in modern cosmology. The generated mass/density map can provide details about the small-scale structure of the Universe and the interaction between galaxies, clusters, and the cosmic web. It holds information about the integrated density fluctuation along the line of sight. Compared with shear two-point correlations, popular applications on the convergence map, such as N-point statistics (Secco et al. 2022), peak statistics (Fan 2007;Dietrich & Hartlap 2010;Liu et al. 2015Liu et al. , 2016Fan et al. 2010;Liu et al. 2023;Zorrilla Matilla et al. 2016;Shan et al. 2017;Chen et al. 2020;Zhang et al. 2022), and Minkowski functionals (Petri et al. 2013;Fang et al. 2017;), can provide complementary information about non-Gaussian density field at late times generated by nonlinear gravitational collapse on small scales. It is usually convenient to implement these methods to the density field directly and obtain the constraint on cosmological parameters and models. The mass maps can also be intrinsically useful. For example, using the DES Science Verification mass map, Clerkin et al. (2017) showed that the one-point distribution of the density field is more consistent with log-normal than with Gaussian. Combining mass maps with the spatial distributions of stellar mass or galaxy clusters enables us to study the relation between the visible baryonic matter and the invisible dark matter. Using mass maps to constrain galaxy bias (Chang et al. 2016), the relation between the distribution of galaxies and matter in turn can aid cosmological probes other than weak lensing. It also enables simple tests for systematic errors in the galaxy-shape catalogs.
The surface density field can be constructed from the background shear field through a linear transformation (Kaiser & Squires 1993). This process is often refined by adopting prior knowledge on the density field through the maximum likelihood method, including Wiener filtering, sparse priors (Leonard et al. 2014;Starck et al. 2015;Li et al. 2021;Price et al. 2021), and deep mass (Jeffrey et al. 2020). However, all these methods focus on deriving the κ map from the shear map, but they focus much less on how to generate a shear map from the shear catalog, which is perhaps an equally important problem. Currently, the shear field is typically made by taking the weighted sum of the shear estimators within a given smoothing kernel. The choice of the kernel size can be quite nontrivial: it should be neither too small so that it keeps enough source galaxies, nor too large so that a reasonably good spatial resolution is reached. A large kernel may also introduce systematic errors in the shear field due to the coupling between the inhomogeneity of the shear field and that of the source density on the kernel scale. In this paper, we aim to solve these problems by proposing a new way to extract the shear map from a shear catalog. Not only do we try to avoid the systematic errors aforementioned, we also hope to reach the optimal statistical uncertainty in the reconstructed shear map. We call this new method PDF-folding (called PF method hereafter) because it is based on symmetrizing the probability distribution function (PDF) of the shear estimators, similar to the PDF-SYM method previously proposed in Zhang et al. (2017).
The rest of the paper is organized as follows. In Section 2 we introduce the PF algorithm. In Section 3 we test the accuracy of this method using the mass map computed from the Illutris simulation (Nelson et al. 2015), and demonstrate the accuracy of the PF method in the presence of an inhomogeneous source distribution. As a check of the actual performance of PF, we apply the method on the real shear catalog around a few massive foreground galaxy clusters in Section 4. We give a brief conclusion in Section 5, and we discuss some outstanding issues in the PF method about which we caution the users.

The PDF-SYM Method for Shear Recovery
The method PDF-SYM is introduced in Zhang et al. (2017) as a new statistical approach to estimating the shear signal from the shear estimators. It aims to achieve the minimum statistical error (the Cramer-Rao bound) without introducing systematic biases. Although it is developed based on the Fourier_Quad shear estimator, the idea is in principle applicable to shear estimators of any form.
To demonstrate the idea, we assume that the shear estimator is simply the galaxy ellipticity e of an ideal form, i.e., the measured e is related to the intrinsic ellipticity e I and the shear signal g via e = e I + g. Note that realistic shear estimators typically require additional corrections, or even take some unconventional forms (Zhang & Komatsu 2011;Sheldon et al. 2017). These changes, however, do not affect our discussion below and can easily be incorporated into the PDF-based algorithms, as shown in Zhang et al. (2017). For simplicity, we have also neglected the subindices of the ellipticity and shear. The basic idea of PDF-SYM is to find the value of a pseudosignalĝ, which can best symmetrize the PDF of the corrected shear estimator = -(ˆ) e e g . The new PDF ofê is related to the PDF of the intrinsic ellipticity e I via = (ˆ)ˆ( ) P e de P e de I I I , and we have = + -(ˆ) (ˆˆ) () P e P e g g . 1 I Because P I is a symmetric function, one can see from Equation (1) that (ˆ) P e is best symmetrized whenĝ is equal to the true shear value g. The symmetry level of the PDF can be quantified by comparing the galaxy number counts within bins that are symmetrically placed on the two sides of zero.
The above algorithm is only well defined for the recovery of a single shear signal. In the following sections, we try to extend it in a way to reconstruct a shear field, which is the main purpose of this work.

Reconstruction of a Shear Field
It is helpful to start with a simple example. We assume that the shear field is one-dimensional, i.e., g(x) = ax + b, in which x is the coordinate in the range of [−1, 1]. For further simplicity, we can assume that the source galaxies are evenly distributed along x. A naive idea of applying PDF-SYM would be to adopt the same form of spatial distribution for the pseudoshear signal, i.e., = + ( )ˆĝ x ax b, and try to find the best values ofâ andb to symmetrize the overall PDF of the galaxies. This is indeed a natural choice for modeling the local distribution of the shear field if the region is small, and this is exactly what we do at the early stage of this project. It turns out that in this case, the parameterb can be constrained by symmetrizing the PDF, but not the parameterâ. The change ofâ has opposite effects on the ellipticity distributions for galaxies located at x < 0 and x > 0, the overall symmetry of the whole PDF is therefore not affected byâ. This is illustrated in Figure 1, in which we show the impact ofâ andb on the PDF. For clarity, in the figure, we split the PDF for the galaxies with x < 0 and x 0, shown in blue and red, respectively. As one can see, the parameterb moves the blue and red populations along the same direction, and it can therefore change the symmetry of the PDF. The parameterâ, on the other hand, moves the two groups of galaxies toward opposite directions, and does thus not affect the overall symmetry of the whole PDF.
A possible remedy for this problem is to invert the sign of the corrected shear estimatorsê for galaxies of x < 0. The symmetry of the resulting new PDF becomes sensitive toâ. Its value can therefore be found by symmetrizing the folded PDF. This is why our new method is called PDF-folding. More generally, for a shear field parameterized by a set of orthogonal functions, we find that one can recover the coefficients one by one, and each time by symmetrizing the PDF that is folded at the places where the corresponding basis function changes its sign.
To realize this idea, we consider a general shear field g(x) in 2D. It can be expanded with a set of orthogonal functions: g(x) = ∑ k a k f k (x). The question is how each mode's magnitude a k can be estimated. To characterize the PDF of the corrected shear estimatorê, we define u i (i = 0, ± 1,K, ± l) as the boundaries of the bins placed symmetrically on the two sides of zero, i.e., u i = − u −i . In total, there are 2(l + 1) bins. Note that the outer boundaries of the two outmost bins are infinite. Assuming we wished to recover a w , the estimated shear field is then written as , where we caution that the summation sign is not present.â w is the presumed value of a w . The number of galaxies N i of the ith bin in our new PDFfolding scheme is defined as For now, we simply assume that the galaxy number densityn is a constant. P is the normalized PDF, which we assume does not Figure 1. The behavior of the ellipticity PDF of galaxies located at x < 0 (in blue) or x 0 (in red) when the parameter a (left panel) or b (right panel) changes. The dotted and solid lines correspond to the PDFs before and after the change of the parameter, respectively. change with position. The form in Equation (2) is quite different from the usual definition of N i , which is simply . The idea of folding is manifested by mixing the galaxies on the two sides of zero in PDF according to the sign of f w (x) (which is known) in Equation (2). The factor |f w (x)| is an additional weight for filtering out the contamination from other orthogonal modes, as we show next.
To find out how the value ofâ w can change the symmetry of the PDF, we calculate Note that integrating without lower and upper bounds means integrating over the whole available range of the variable. Using Equation (1), we can relate the function P with the intrinsic PDF P I (which is symmetric) as . The last step above is from Taylor's expansion to the first order because the shear signal is assumed to be small. Using the results in Equation (4), we obtain Using the result of Equation (5) in Equation (3) and the orthogonality of f k (x), we obtain It is clear that the folded PDF is best symmetrized when = a a w w . This is true for all the bin pairs. Therefore, to estimateâ w , we just need to minimize the following χ 23 : The above procedure is repeated for eachâ w to recover the whole shear field. This is the basic implementation of the PF method.
For an inhomogeneous galaxy distribution, it is not hard to check that a similar conclusion can be reached if we simply replace the weighting of the galaxy number in Equation (2) from |f w (x)| to |f w (x)|/n(x), where n(x) is the galaxy number density. The form of Equation (6) remains the same, but without the factorn. This, however, turns out not to be the end of the story. The formulation can be further improved to achieve two more goals: 1. proper treatment of masks, and 2. minimum statistical noise. This is what we discuss next.

Optimal Reconstruction Method
If we chose the weight to be |f w (x)|/n(x), we can immediately identify a problem: what if there are masked areas? A simple option would be to skip the masked areas. However, in this case, the orthogonality of the functions f k (x) would not lead to Equation (6) because of the incomplete integration domain. Another option would be to use a set of orthogonal functions that are defined in the domain excluding the masked areas. This would be a fine choice, except that the form of f k (x) could be quite complicated.
In this work, we take another route: we fill out the masked area with galaxies of a certain number density, e.g., the average number density of the whole field. Each such galaxy is given a random shear estimator/ellipticity, i.e., the shear field is assumed to be zero in the masked areas. This procedure guarantees the completeness of the domain of integration, and therefore, the validity of the PF method.
More generally, we are interested in finding the optimal way of extracting the shear-field information from the shear catalog in the case of inhomogeneous galaxy distribution. To do so, we find that the following generalized form of series expansion is useful: To recover the coefficient a w , the galaxy weight, and therefore, the definitions of the PDF, should be updated accordingly as In this case, one can show that the difference between the opposite bins is given by Following similar calculations as in Section 2.2, one can show that Equation (6) still holds (without the factorn). Apparently, we now have the exponent α as an additional degree of freedom in our formalism. It turns out that when α = 1/2, the statistical uncertainty of a w reaches its minimum. It is therefore our best choice for the shear-field reconstruction. We show the details of our proof in Appendix A.
The above discussion is all about reconstructing the shear field at a given background redshift. More generally, we should consider the redshift distribution of the source galaxies. For the single thin-lens case (which is what we consider in this work), the shear fields at different redshifts are related through a simple rescaling with the corresponding critical surface densities. The problem is therefore still two-dimensional, and we only need to choose a certain background redshift as a reference. The optimization of the PF method in this case is given in Appendix B.
In the rest of the paper, we demonstrate some advantages of the PF method using numerical examples, and we show some shear/density field reconstruction examples from the real data.

Numerical Examples
In this section, we demonstrate several advantages of the PF method with data from the Illutris-1-Dark simulation (Nelson et al. 2015). The side length of the simulation box is 75 Mpc/h (comoving). We use the snapshot with ID 120 corresponding to z = 0.2. We select a 6 × 6(Mpc/h) 2 (comoving) part containing cluster-like structure and calculate the projected surface density by integrating over the whole box size. The shear field is then deduced from the mass distribution by assuming the source galaxies are all at z = 0.50 and placed on a 120 × 120 grid as shown in Figure 2 for the g 1 component. This shear field is applied to a large number of background galaxies whose intrinsic ellipticities are generated according to Miller et al. (2013). The source galaxies are randomly placed inside the foreground area. To reduce the shape noise, we group every 5 × 5 grid area to form a 24 × 24 shear map.
As our first test, we adopt three strategies to reconstruct the shear map from the shear catalog: 1. by taking the local averages of the shear estimators (called Local_AVE hereafter); 2. by using the PDF-SYM method, which is designed for optimizing the stacked shear signal on the local ensemble of the shear estimators (called Local_PDF hereafter); and 3. by the PF method. In the last method, we expand the shear field with the Fourier series.
In our first test, the shear field is applied onto 10 7 background galaxies whose ellipticities follow the form of disk-dominated galaxies. We set the parameter e max in Miller et al. (2013) to be 1.0. The recovered shear maps (g 1 ) with the methods of Local_AVE, Local_PDF, and PF are shown in the upper panel of Figure 3. The corresponding residuals are accordingly shown in the lower panels. One can see that the PDF-based methods generally generate somewhat less noise than the Local_AVE method. The advantage would be more obvious if the PDF of the ellipticities had a more significant extension to high values.
Note that although in this simple example, the Local_PDF method seems to work similarly well as the PF, its performance is actually more sensitive to the background galaxy number density, i.e., there should be enough galaxies in each grid cell, as we show in the real data examples in Section 4. The PF method avoids this issue by using the background galaxies in the field altogether to form the PDF. It therefore does not require a large smoothing kernel or grid size, and can keep information on smaller scales in principle.

Effect of an Inhomogeneous Galaxy Distribution
The distribution of the background galaxies is often inhomogeneous (as shown in some examples of Section 4). In traditional methods, if both the shear signal and the background density vary significantly within the smoothing kernel, we expect the coupling of these two types of fluctuations to induce systematic errors in the recovered shear field. In the new PF method, we find that this bias can be largely avoided.
To study this effect, we let the background galaxy density vary with the vertical axis y as + ( ) ky const sin , where k = 2π/ 5, whose period equals the grid size of the recovered shear field. To enhance the signal-to-noise ratio of the systematic bias, we use 5 × 10 7 galaxies with ellipticities set in the range of [0.05, 0.1], 90% of which follow the ellipticity distribution of the disk-dominated galaxies, and the other 10% follow that of the bulge-dominated ones. In Figure 4 we show the original shear field in the upper left panel, and the background galaxy distribution in the upper right panel (the unit of which is arbitrary). The lower two panels show the differences between the recovered shear field and the original one for the Local_PDF method (left) and the PF method (right). One can see that the shear residuals of the Local_PDF method are quite significant and correlated with the signal. In contrast, the PF method performs quite well in the presence of the background inhomogeneity. Note that in the current study, the smoothing kernel for the Local_PDF is simply the top-hat function within the grid size (of the recovered shear field). If a larger kernel is chosen, the effect is more obvious. For the PF method, we use the expansion defined in Equation (8) with α = 1/2. The local number density around each galaxy is calculated by averaging the galaxy numbers within its neighboring 3 × 3 finer grid cells (the grid on the original 120 × 120 map).

Effect of a Mask
Masks are just special cases of background inhomogeneity. For similar reasons, the Local_PDF operation can also generate  systematic biases near masks. In this test, we assign a masked area in the central region of the simulated shear map from Illutris, as shown in the upper left panel of Figure 5. It can also be seen from the upper right panel of the same figure, in which we show the distribution of the background galaxy density (the unit is again arbitrary). We apply the shear field to 10 7 background galaxies in total, the ellipticities of which are generated in the same way as those in Section 3.1. The background galaxies are evenly distributed in the whole area except for the masked region. In the lower panels of Figure 5, we again show the performances of the Local_PDF method (left) and the PF method (right) in terms of the differences between their recovered shear fields and the original field. We can see that there are significant residuals near the masked area in the Local_PDF method. As a comparison, the shear field from the PF method shows no such residuals, demonstrating a consistently better performance. We note that to use PF in this case, one needs to add simulated galaxies of random (or zero) ellipticities in the masked area. The number density of the simulated galaxies inside the mask is chosen to be comparable to the mean of the whole field.

Step-by-step Procedure of the PF Method
Here, for simplicity, we still use the ellipticity e to represent the shear estimator and neglect its subindex. The application of the PF method for reconstruction of the shear field follows the following steps: 1. Set up a rectangular grid for the shear field.
2. Count the galaxy number density n(x) in each grid. 3. Fill up the masked area with galaxies of the average number density and random ellipticities/shear estimators.
4. Determine a set of orthogonal and complete functions, e.g., Fourier series f w (x), to parameterize the shear field as = å , where w is the index of each Fourier mode. The rest of the steps are about determining the coefficients a w one by one.
5. Each a w is determined by changing its assumed valueâ w , so that the PDF of the whole background galaxy sample is best symmetrized. The pseudo-shear signal for a galaxy is given by . The resulting shear field can be converted into the κ field through, e.g., the Kennicut-Schmidt (K-S) inversion algorithm.
Note that in the PF method, the choice of the pixel size is rather arbitrary. A large pixel size is fine for those who are only interested in the large-scale behavior of the shear field. On the other hand, if the pixel size is small, for instance, each pixel only contains a handful of source galaxies on average, the method still works fine. In this case, although the recovered shear field in each pixel would become more noisy, the signalto-noise ratio of the large-scale modes would not be affected.

Application on Real Data
In this section, we demonstrate the performance of the PF method with real data. We adopt the shear estimator of the form defined in the Fourier_Quad shear measurement method (Zhang 2010;Zhang et al. 2015), which has been shown to work well with the PDF-SYM method (Zhang et al. 2017). Note that although the Fourier_Quad shear estimator takes an unconventional form, it works equally well as the ideal shear estimators discussed in Section 2, as shown in Zhang et al. (2017). In other words, the whole discussion in Section 2 can be almost trivially replicated for the Fourier_Quad shear estimator.
Our shear catalog comes from the HSC galaxy survey Data Release 2 (Aihara et al. 2019). The images are processed by the Fourier_Quad pipeline in five bands (grizy) individually. The Fourier_Quad pipeline is previously used to process the   imaging data of CFHTLenS ) and DECaLS , and the results of which have all passed the accuracy test using the field-distortion effect. We find that the same pipeline also performs quite well for the HSC data. The details of the measurement are presented in a separate paper (Liu et al. 2023, in preparation). In this work, we use the shear catalog of the rizy four bands. Note that in Fourier_Quad, shear estimators from different exposures or bands are treated as independent even if they are from the same galaxy. The foreground lenses in our examples are three massive galaxy clusters chosen from the Sloan Digital Sky Survey (SDSS) group catalog (Yang et al. 2007(Yang et al. , 2008(Yang et al. , 2009. Their information is shown in Table 1, including their angular positions, redshift (z lens ), and mass (from the catalog). To reconstruct the surface density map, we first build up the background shear field at a reference redshift (z lens + 0.3) following the more general version of the PF method (taking into account the background redshift distribution) in Appendix B, in which a modified version of the step-by-step procedures given in Section 3.3 is also provided. The convergence/density map is then converted from the shear map through the standard procedures given in Kaiser & Squires (1993). One complication is that what we actually measure from the background galaxy shapes are not the shear, but the reduced shear (Zhang 2011). To improve the accuracy of the recovered density map, one therefore needs to modify the shear estimators with the factor 1 − κ, where κ is the convergence just derived from the shear map. Usually, with several iterations, the shear map and κ map can both achieve stable solutions (Liu et al. 2014). To avoid possible contamination from the cluster members, we only use background galaxies with redshifts higher than z lens + 0.1.
For comparison, the Local_PDF method, Local_AVE method and the PF method are implemented. The results are shown in Figures (6), (7), and (8) for clusters A, B, and C, respectively. All the panels in the figures have dimensions of 1°× 1°. In the gray area of each figure, we show the results of Local_PDF, Local_AVE, and PF in the first three rows, with a grid size of ¢´¢ 2.5 2.5. We note that in the Local_AVE method, we use the ensemble average of the local Fourier_Quad shear estimators to evaluate the local shear field g 1 and g 2 , as shown in Equation (28) of Zhang et al. (2017). This approach is usually not encouraged to use for the Fouier_Quad shear estimators because it uses unnormalized quadruple moments of the galaxy, and therefore variations in galaxy luminosity/size can lead to large scatter in shear recovery (this fact, however, does not cause any problem in PDF-based methods such as Local_PDF and PF here). Here, for the purpose of comparison, we still present the results of Local_AVE, but excluding the highest and lowest 8% of the shear estimators to reduce the noise. We would not recommend doing so in practice, as the selection would inevitably introduce the local shear bias.
The left two columns show the recovered shear fields g 1 and g 2 of the two methods. The third and fourth columns show the surface density maps from the E (Σ E ) and B (Σ B ) modes of the convergence fields, respectively. The rightmost panel (out of the gray region) shows the background galaxy image density used in the same area of the sky. It is clear from the figures that  Table 1. The first two columns show the g 1 and g 2 fields for the Local_PDF method (the first row), Local_AVE method (the second row), and the PF method (the third row). The third and fourth columns show the corresponding surface density fields of the E (Σ E ) and B (Σ B ) modes, respectively. The rightmost column shows the distribution of the background galaxy number density in the same area of the sky. The fourth row shows the density fields that were also recovered by the PF method, but with a higher spatial resolution. Figure 7. Same as Figure 6, but for cluster B listed in Table 1. Figure 8. Same as Figure 6, but for cluster C listed in Table 1. the results of the PF method are generally less noisy than those of Local_PDF. Some of the extreme values on the shear maps of Local_PDF arise for two reasons: 1. lack of background galaxies, and 2. outliers in the shear catalog. These two factors fortunately do not affect the performance of the PF method. We consider this feature a significant advantage of PF. It allows us to further increase the spatial resolution of the density map. For example, in the bottom rows of Figures 6, 7, and 8, we show the reconstructed density maps using PF with a much finer grid size of ¢´¢ 0.5 0.5. The final results are smoothed with a Gaussian filter of s = ¢ 1.65 to increase the visibility of the density structure. In this case, one can see that the central overdensities of the clusters stand out more clearly. In contrast, we find that it is not feasible to further increase the spatial resolution in the Local_PDF method.

Conclusion and Discussions
Weak lensing provides a direct way of reconstructing the foreground density distribution. So far, the existing algorithms mostly focus on refining the conversion from the shear field to the density map. The construction of the shear field from the shear catalog, however, receives much less attention. Currently, it is done by locally averaging the shear estimators within some smoothing kernel. The choice of the kernel size can be quite nontrivial: on one hand, it should be large enough to guarantee sufficient background galaxies, and on the other hand, it is desirable to keep the kernel size small for a good spatial resolution. The conventional method may also induce systematic errors in the shear field due to the coupling between its fluctuation and the inhomogeneity of the background source galaxies.
In this work, based on the method of PDF-SYM in Zhang et al. (2017), we propose the PDF-folding method as an alternative way of recovering the shear field from the shear catalog. The details of the new method are given in Section 2. It is significantly different from the conventional weighted-sum methods. The main differences include the following points: 1. Instead of recovering the local shear value, PDF-folding reconstructs the coefficients of the orthogonal modes that are used to decompose the shear field as a whole. It therefore does not require any smoothing kernels.
2. Inhomogeneous distributions of the source galaxies (as well as masks) are automatically taken into account in the method without incurring systematic biases.
3. The statistical errors on the coefficients of the orthogonal modes can reach the theoretical minimum in general observing conditions, i.e., with an inhomogeneous source distribution and realistic shape noise.
The above points are demonstrated with analytical reasoning (see Section 2 and the appendixes) as well as with numerical examples using Illutris-1-Dark simulation data (see Section 3). As real-world examples, we applied in Section 4 the PDFfolding method on the shear catalog of the HSC survey (constructed from the Fourier_Quad pipeline) to reconstruct the surface density maps around three massive galaxy clusters identified in the SDSS. Comparing this to the local reconstruction of the shear field, we show that PDF-folding shows a consistently better performance. This is mainly because the new method is much less susceptible to the local inadequacy of source galaxies or outliers in shape noise.
The main caveat in PDF-folding is in the calculation of N i(>0) −N −i in Section 2, in which we expand the PDF of the shear estimators to the first order in shear. It implies that the shear signal is small enough compared to the dispersion of the shape noise. This condition is not necessarily well satisfied in the neighborhood of massive galaxy clusters, and therefore, there can be a small number of systematic errors. This is not a serious issue, however, as the very central regions of galaxy clusters can be simply masked if necessary. This expansion in PDF may also fail when the PDF itself has singularities. This has been shown as an advantage in the PDF-SYM method discussed in Zhang et al. (2017) because it can in principle lead to an almost infinite signal-to-noise ratio in shear recovery. In PDF-folding, however, we find that these singularities in the PDF lead to significant systematic errors in the recovered shear field because the failure of the expansion invalidates the cancellation of the couplings between different orthogonal modes. Nevertheless, this problem is no serious concern in practice either, as the realistic PDF of the shear estimators rarely contains singularities.
Finally, we did not consider the impact of the shear measurement errors on the recovered shear field. One cause of the error is the uncertainty in the point-spread function, which could depend on the local density of the stars and might therefore induce a variation in the shear field. This topic will be studied in a future work. In the PF method, the value of a w is estimated by minimizing the χ 2 , which is defined as The denominator in the above formula can be worked out in the following way: in which the δ sign refers to the deviation from the statistical mean. In this calculation, we have also neglected the influence of the lensing effect, which yields a slightly nonzero mean for N i(>0) −N −i . According to Equation (9) In the limit of a very small bin size, we can turn the summation in the above formula into an integral form as algorithm, and the κ field can be further converted into the surface density distribution. A point to note is about our assumption that the form of the PDF P I does not depend on the redshift. This is not true in practice. This fact would require us to further change the weighting defined in Equations (B7) and (B8) to take into account the redshift-dependent shape noise for the purpose of achieving the optimal statistical uncertainty. We may study this issue in a future work. Our current formalism of PF is at least close to the optimal form. It is also not hard to show, e.g., with Equation (B4), that even if P I depends on redshift, the PF method would not generate systematic biases at the first order of shear.