Shear calibration biases in weak-lensing surveys

Abstract

We investigate biases induced by the conversion between the observed image shape and shear distortion in current weak-lensing analysis methods. Such overall calibration biases cannot be detected by the standard tests such as E/B decomposition or calibration with stars. We find that the non-Gaussianity of the point spread function has a significant effect and can lead to an error of up to 15 per cent on the linear amplitude of fluctuations σ₈, depending on the method of analysis. This could explain some of the discrepancies seen in recent amplitude determinations from weak lensing. Using an elliptical Laguerre expansion method we develop a re-Gaussianization method that reduces the error to a calibration error of the order of 1 per cent, even for poorly resolved galaxies. We also discuss a new type of shear selection bias, which results in up to roughly an 8 per cent underestimation of the signal. It is expected to scale with redshift, inducing errors in the growth factor extraction if not properly corrected for. Understanding and correcting for such effects is crucial if weak lensing is to become a high-precision probe of cosmology.

gravitational lensing, methods: data analysis

Introduction

Gravitational lensing is a powerful tool for probing the large-scale matter distribution in the Universe (see, e.g., Bartelmann & Schneider 2001 for a review). As such it offers the possibility of breaking the degeneracies between many cosmological parameters (e.g. Hu & Tegmark 1999) and serves as a powerful complementary method to the cosmic microwave background anisotropies. Its main advantage is that it is sensitive to the dark matter directly and so avoids many of the complications present in other cosmological probes at low redshift such as galaxy clustering or cluster surveys. Recent first detections of shear correlations (Bacon, Refregier & Ellis 2000; Van Waerbeke et al. 2000; Rhodes, Refregier & Groth 2001; Brown et al. 2003; Hoekstra et al. 2002; Jarvis et al. 2003) have already demonstrated the promise of this method. Upcoming surveys such as the NOAO deep survey (www.noao.edu/noao/noaodeep/), the CFHT legacy survey (www.cfht.hawaii.edu/Science/CFHLS/) and the Deep Lens Survey (dls.bell-labs.com/) will go beyond simply detecting cosmic shear and will map it over large areas of the sky. Finally, proposed projects such as DMT/LSST (www.dmtelescope.org/darkhome.html), Pan-STARRS (www.ifa.hawaii.edu/pan-starrs/) or SNAP (snap.lbl.gov) will map a large fraction of the sky with high depth, achieving unprecedented statistical precision.

Advances in statistical precision require parallel progress on the systematic errors if the information from the current and upcoming experiments is to be exploited fully. The effects one is looking for are extremely small, so any small errors in the calibration or data reduction can lead to a large error in the final result. The most commonly used methods of reduction so far are based on work by Kaiser, Squires & Broadhurst (1995) (KSB95). Recently, there have been several new approaches proposed dealing with various aspects of this problem developed by Kaiser (2000), Bernstein & Jarvis (2002) (BJ02) and Refregier & Bacon (2003). What is somewhat unclear from this work is how much difference these new approaches make and what are the remaining errors in them. While most of these papers present at least some tests, these are often performed using simplified models and do not attempt to quantify the error as a function of size of the galaxy or their ellipticity. The original KSB95 method was calibrated using simulations, and other authors have investigated calibration issues in some of the variants and extensions of KSB95 numerically (Kaiser 2000; Bacon et al. 2001; Erben et al. 2001). However, the calibration of more recent methodologies such as BJ02 has not yet been studied in detail, nor have the various methods been compared in the same simulation. [Many of the simulations in the literature are not directly comparable since they have different assumptions concerning the distribution of galaxy shapes and sizes and the point spread function (PSF) and noise properties.] The goal of this paper is to address some of these issues in more detail. We pay particular attention to the problem of point spread function dilution in realistic PSF and galaxy shape models, which only affects the overall calibration of the shear (if seeing is not significantly variable across the lensing field) and does not show up in the systematics tests.

Extracting shear from the images is not the only place where a bias might arise in the weak-lensing analysis. Other possibilities include (but are not limited to) intrinsic correlations, errors in the source redshift distribution and so-called selection bias, where a galaxy is preferentially selected when its shape is aligned with the PSF distortion. However, there is another possible bias caused by the fact that a galaxy may also be preferentially selected when aligned orthogonally to the intrinsic shear, because detection algorithms preferentially detect circular objects. This bias depends on the detection threshold and the relative angular size of the galaxy and PSF; it becomes unimportant for poorly resolved galaxies (i.e. galaxies much smaller than the PSF) because in this case the size and shape of the image of the galaxy are determined by the PSF irrespective of the physical orientation of the galaxy. Since it is directly correlated with the signal it is very difficult to detect it using the standard systematics tests such as the 45° rotation test (E/B decomposition) and must be corrected for using simulations or some other method that explicitly takes it into account.

This paper is organized as follows. In Section 2 we review past approaches to PSF corrections, present the formalism for PSF correction based on the elliptical Laguerre expansion, and derive two approximate methods to correct for non-Gaussianity of the PSF. We compare various methods in the literature and estimate their biases using simulations. In Section 3 we discuss a shear selection bias and how it can affect the amplitude of fluctuations. Conclusions are presented in Section 4 and in the appendices we describe technical details on convolution in the Laguerre expansion formalism, the linear PSF method and the KSB95 method.

Psf Corrections

In actual measurements, we do not measure directly the moments of intrinsic galaxy images, but rather we measure the moments of the galaxy image after convolution with the PSF of the telescope. We will denote the intrinsic galaxy image by ƒ(x), the PSF by g(x) and the measured galaxy image by I(x). Our eventual objective is to take the PSF g(x) (which can be reconstructed by considering the images of stars in the field) and the observed image of the galaxy I(x), and construct from these the ellipticity e_ƒ of the intrinsic image that can (with an appropriate calibration factor) be used as an estimator for the local shear.

Many approaches to this problem appear in the literature. The common theme in most of them is to measure the moments of the image and PSF using some weight function:  

(1)

where x= (x, y) is the position on the sky, w is the weight function and the P_i are polynomials in x and y. These moments M^(I)_i, and their counterparts M^(g)_i for the PSF, are then used to compute the ellipticity estimator. The methods mainly differ in their answers to the following questions.

• (1)

  Is the weight function a circular Gaussian, an elliptical Gaussian or some other shape?

• (2)

  Is the weight function used for the PSF the same as that for the galaxy, or is it matched to the size of the PSF?

• (3)

  Does the method (partially) deconvolve the galaxy image or does it directly correct the galaxy moments for the effect of the PSF? In the latter case, which corrections are performed using perturbation theory, which are performed by comparison with simulations, and which are performed using simple analytical results valid in limiting cases?

Question 1 has typically been decided in favour of a Gaussian weight due to its rapid convergence to zero at large radii (hence the polynomial moments are all convergent) and the absence of singularities, as well as general mathematical convenience. The majority of recent theoretical (Kaiser 2000; Refregier & Bacon 2003) and observational (Bacon et al. 2000; Van Waerbeke et al. 2000; Rhodes et al. 2001; Hoekstra et al. 2002) analyses have followed KSB95 in using a circular Gaussian weight function. The elliptical Gaussian weight function was introduced by BJ02; although it is considerably more complicated to work with, a general ellipse can be a better match to the actual shape of an object than a circle. This property is important for those methods that treat the deviation of an object or PSF from Gaussianity via perturbation theory, including our linear and re-Gaussianization methods.

There is little consensus in the literature on the best answer to Question 2. The advantage of measuring the PSF with its own weight function, independently of the galaxy, is that the PSF weight function can be matched to the size (and, in the case of elliptical moments, the shape) of the PSF. This is necessary for methods (e.g. re-Gaussianization) that treat the non-Gaussianity of the PSF perturbatively. The disadvantage of measuring moments of the PSF with a weight function matched to the PSF is that features of the PSF outside the scale radius of the weight function are suppressed. Features at radii large compared with the characteristic scale of the PSF are present in nearly all real-world PSFs: they are introduced as diffraction rings in space-based observations, and as a power-law (g∝r^−11/3) tail from Kolmogorov turbulence in ground-based observations. If the galaxy is of size comparable to or somewhat larger than the PSF, these features can fall inside the scale radius of the weight function of the galaxy, with potentially disasterous results for attempts to unravel the effects of the PSF on the measured galaxy moments (see Section 2.2). The higher-order moments of the PSF must therefore be considered in order to take into account features at large radius if a PSF weight matched to the size of the PSF itself is used. Numerical tests (Hoekstra et al. 1998) indicate that the KSB95 method supplemented with a pre-seeing shear polarizability calculation (Luppino & Kaiser 1997) works best when the PSF is measured with the same weight function as the galaxy; this is unsurprising, since this method does not consider any moments higher than the fourth order (see Appendix C).

Finally, we arrive at Question 3: how do we convert the measured moments of the galaxy and PSF into an ellipticity of the intrinsic galaxy image? Most methods, including Luppino & Kaiser (1997), Hoekstra et al. (1998) and BJ02, have corrected the galaxy moments for the effect of the PSF using either perturbation theory, simple analytical corrections or some combination thereof. (The original KSB95 work used simulations to compute the diluting effect of the PSF on the galaxy.) Refregier & Bacon (2003) developed a deconvolution method based on the circular Laguerre expansion of the PSF and galaxy (see Section 2.2), writing the PSF convolution as a matrix acting on the vector of galaxy moments and cutting off the hierarchy of Laguerre moments at sixth order to prevent noise amplification in the deconvolution (matrix inversion). BJ02 have investigated matrix inversion methods as well using the raising and lowering operator formalism, with particular attention paid to the effects of the matrix inversion on noise. In this paper, we introduce and investigate two perturbative approaches, expanding around the case of elliptical Gaussian PSFs and galaxies. The ‘linear’ method (Appendix B) is essentially a first-order approximation to the matrix inverse of Refregier & Bacon (2003), generalized to elliptical moments. It improves significantly on the same order method proposed by BJ02 in their Appendix C and should be of particular interest to those using Sloan Digital Sky Survey photo reduction for weak lensing. The ‘re-Gaussianization’ method (Section 2.2) treats the deviations of the PSF from Gaussianity perturbatively in real space, without first measuring moments. In our tests this was the most successful method and the only one that comes close to the precision requirements of the future surveys. The introduction of non-linear calculations such as higher-order perturbation theory or direct matrix inversion can be expected to further improve the accuracy of the method because the full non-linear dependence of measured moments on the PSF is taken into account; however, one must be careful to avoid amplifying noise or creating unstable estimators.

We organize this section as follows. The adaptive moment determination is discussed in Section 2.2. In Section 2.2 (and in Appendix A), we explore the generalization of the Laguerre expansion of Refregier (2003) and BJ02 to the case of an elliptical weight. In Section 2.2, we explore the limitations of BJ02's analytical correction for PSF dilution; the re-Gaussianization method is introduced in Section 2.2. We conclude by comparing the methods in simulations (Section 2.2).

Adaptive moments

A substantial literature exists on approaches using circular Gaussian-weighted moments for g and I. BJ02 introduced ‘adaptive’ moments that use an elliptical weight function, the shape of which is matched to that of the object. The adaptive moments of an image I(x) are determined using the Gaussian with the least-squares deviation from the image, i.e. we minimize  

(2)

over the sextuplet of quantities (A, x₀, M). Defining  

(3)

we note that minimization of E yields the moments satisfying  

(4)

Because of the presence of the w, the adaptive moments are weighted. [It is simpler to work with unweighted moments, but these have divergent noise (KSB95).]

The ellipticity of an object is then defined in terms of the adaptive second moments of I(x):  

(5)

the spin-2 tensor e= (e₊, e_×) is called the ellipticity tensor and the scalar T^(I) is called the trace.

Elliptical Laguerre expansion

This section explores aspects of the Laguerre expansion introduced by BJ02. In particular, we rework BJ02's formalism into a form useful for studying elliptical objects. A Laguerre expansion is constructed around a Gaussian with some symmetric covariance matrix M; we begin by linearly transforming the Cartesian (x, y) coordinate system to another (u, v):  

(6)

where D=M_xxM_yy−M²_xy is the determinant of M and . This provides us with the property that u²+v²=x^TM⁻¹x. The Laguerre basis functions are  

(7)

where  

(8)

where m=p−q is the angular momentum quantum number and the + sign is used in u± iv if m > 0 and the − sign is used if m < 0. (Both signs are equivalent if m= 0.) Equations (7) and (8) can be seen to be equivalent to BJ02's equations (6)–(8) through the use of the explicit expression for the associated Laguerre polynomials. Note that Λ_pq is a polynomial in two variables with complex integer coefficients and Λ_pq is the complex conjugate of Λ_qp. In Table 1 the Λ polynomials are given through fourth order; in the last column we give the result in polar coordinates. We note that the Λ polynomials are related to Refregier (2003)'s polar Hermite polynomials according to Λ_pq(r cos φ, r sin φ) =H_p,q(r) e^imφ.

Table 1.

The α polynomials up to fourth order.

Finally, we consider the Fourier transform  

(9)

The Fourier transform of a Laguerre basis function is given by BJ02's equation (6–62); using our Λ polynomials, it is  

(10)

where N=p+q is the order of the polynomial. The variables k_u and k_v are given by (k_u, k_v) =M^1/2x and satisfy k_uu+k_v=k·x.

Now any image can be expanded in the Laguerre basis functions:  

(11)

(Note that, through a x as an argument.) We note that there is a complex-conjugate relationship among the moments b_pq=b*_qp if the image I is real-valued. Furthermore, by orthogonality of the ψ_pq basis functions, there is an inverse transform for the coefficients:  

(12)

We will find it useful to use the reduced quantities β_pq=b_pq/b₀₀. BJ02 observe that the Laguerre coefficients associated with the adaptively determined covariance M satisfy β₀₁=β₁₁=β₀₂= 0. One can see this directly from the least-squares residual principle by differentiating equation (2) with respect to (A, x₀, M) and setting the result to zero.

PSF corrections

PSF correction is typically performed by correcting the moments of the measured object using the moments of the PSF. There are generally two cases to consider here: either the PSF moments are measured with an elliptical Gaussian weight adapted to the PSF shape rather than the object shape, or they are measured with the same weight function as the object. The former approach is used by BJ02, and we will consider it in detail here; we defer investigation of PSF moments determined with the same weight function as the object.

Another issue that arises is whether to use a rounding-kernel method to symmetrize the PSF. Here we assume that the best-fitting Gaussian to the PSF is circular; this could be achieved using a rounding kernel. An alternative is to apply a shear (to both the image and PSF) to make the best-fitting Gaussian to the PSF circular, and then apply an inverse shear after making the ellipticity measurement. (The transformation of ellipticities is performed according to equation B17.) The rounding kernel method may be better suited to elimination of spurious power from PSF anisotropy, since a rounding kernel can be designed to eliminate arbitrarily many spin-2 moments of the PSF (at least in principle, see, e.g., BJ02); however, rounding kernels increase the effective PSF size which is not desirable. Circularizing the best-fitting Gaussian to the PSF with a shear does not increase PSF size, but the removal of spin-2 moments of the PSF (and hence of the associated spurious power) is incomplete as only the b₀₂ moment is forced to zero by this method.

If the galaxy and PSF are both elliptical Gaussians, the PSF correction is trivial because the covariance matrix of the galaxy and PSF add to yield the covariance of the observed image. In this case, we have M_I=M_ƒ+M_g and hence  

(13)

If we were using unweighted moments, this scheme would be valid for any PSF and galaxy. For weighted moments, equation (13) fails, and a correction for non-Gaussianity of the galaxy and PSF is necessary.

BJ02's prescription (see their Appendix C) for PSF corrections is as follows: begin by applying a rounding kernel; then the ellipticity of the measured object e^(I) is transformed into the ellipticity of the intrinsic object according to e^(ƒ)=e^(I)/R, where  

(14)

This equation is designed to reduce to the correct form in the limit of well-resolved objects with Gaussian PSFs, and in the limit of any Gaussian PSF and Gaussian object (whether or not it is well resolved, see Appendix C of BJ02); however, no rigorous derivation has been provided. Indeed, we can show that for non-Gaussian PSFs, the above equation may be subject to significant error. The measured object has intensity  

(15)

where  

(16)

is the unweighted nth moment of the PSF. (Note that formally q is infinite.) In the limit of a well-resolved object, the effect of the PSF is dominated by the second (n= 2) moment because the series in equation (15) converges rapidly and the first moment vanishes (we assume without loss of generality that the centroid of the PSF lies at the origin). In this limit, BJ02 showed that  

(17)

The relation between q_xx+q_yy and the Laguerre moments can be computed:  

(18)

where ‘⋯‘ represents both higher-order terms in β^(g)₂₂ and terms from the higher-order Laguerre moments b^(g)_pq with p > 2 or q > 2. One would thus expect that for well-resolved galaxies, the correction for non-Gaussianity of the PSF should take the form:  

(19)

where the numerator of the quantity in braces is valid to first order in β^(g)₂₂ and to zeroth order in the higher moments. We note that the BJ02 correction does not exhibit this behaviour; its numerator has the form 1 − 2β^(g)₂₂ instead of 1 + 4β^(g)₂₂. Thus we expect that for PSFs with positive kurtosis β^(g)₂₂ > 0, equation (14) will overestimate the resolution factor and hence underestimate the shear. This expectation is verified in Section 2.2. Note additionally that the above argument does not apply for the less well-resolved galaxies (i.e. those for which T_g is comparable to T_I); in this case, the factor of 1 + 4β^(g)₂₂ is not correct. The correct factor is worked out in Appendix B.

Several methods could be used to improve equation (14). One approach is to compute the correction to the galaxy ellipticity to linear order in the β_pq; this approach is pursued in Appendix B. Unfortunately, realistic PSFs (e.g. turbulence- or diffraction-limited PSFs) have significant power at large radii and hence the convergence of the Laguerre expansion is rather slow. This is manifested in Section 2.2, where significant errors are found from terminating the Laguerre expansion after the 22 moment. Therefore, in the next section we will pursue an alternative approach: we fit a Gaussian to the PSF, compute the residuals and attempt to correct the galaxy image for the effects of these residuals. Then we apply equation (14) to this corrected galaxy image. As shown in the simulations (Section 2.2), this method yields good performance for non-Gaussian PSFs and galaxies, except for the most elongated galaxies.

Re-Gaussianization

Let us suppose that the PSF g can be approximated by a Gaussian G with some covariance M_G:  

(20)

We may then define the residual function ε(x) =g(x) −G(x). It follows that the measured image intensity will satisfy I = g ⊗ ƒ = G ⊗ ƒ +ε⊗ƒ, where ⊗ represents convolution. This can be rearranged to yield  

(21)

This equation thus allows us to determine the Gaussian-convolved intrinsic galaxy image I′, if we know ƒ. At first glance this does not appear helpful, since if we knew ƒ it would be trivial to determine G⊗ƒ. However, ƒ appears in this equation multiplied by (technically, convolved with) a small correction ε, so equation (21) may be reasonably accurate even if we use an approximate form for ƒ. The simplest approach is to approximate ƒ as a Gaussian with covariance:  

(22)

where M_I and M_g are the adaptive covariances of the measured object and PSF, respectively. Then we can define  

(23)

The adaptive moments of I′ can then be computed, and the PSF correction of equation (14) can then be applied to recover the intrinsic ellipticity e^(ƒ).

There are several ‘knobs’ that can be adjusted on the re-Gaussianization method. First, there is the choice of Gaussian approximation to the PSF, equation (20). We have chosen to set M_G to be the adaptive covariance of the PSF, i.e. the covariance of the best-fitting Gaussian. However, rather than using the normalization of the best-fitting Gaussian, we have normalized equation (20) to integrate to unity. While this increases the overall power ∫(ϵ²) of the residual function, it yields ∫ ϵ = 0, which ensures that for well-resolved objects (i.e. objects for which the PSF is essentially a δ-function), the ‘correction’ε⊗ƒ₀ applied by equation (23) does not corrupt the image I. One could also try a more sophisticated function ƒ₀ than a Gaussian, perhaps based on a partial deconvolution of the PSF. We have not explored the possibilities here, but it is clear that an attempt to approximate ƒ will have to make some assumptions regarding the behaviour of ƒ in the regions of Fourier space where the PSF has little power (i.e. where the optical transfer function is close to zero).

Simulation

We may test the PSF correction algorithms via simulation for common categories of galaxy and PSF. The simulation procedure is as follows:

• (i)

  construct an intrinsic image ƒ of the test galaxy;

• (ii)

  construct a sample PSF g;

• (iii)

  convolve the test galaxy and PSF, then pixelize this to yield a ‘measured object’I;

• (iv)

  compute the adaptive moments of the test galaxy ƒ, and measured galaxy I;

• (v)

  apply the PSF correction scheme to the moments of I to yield the ‘recovered’ ellipticity of the test galaxy;

• (vi)

  compare the true and recovered ellipticities e^(ƒ) of the test galaxy; and

• (vii)

  propagate the ellipticity error through to an error in the shear estimate.

Our simulations here will compare five methods of recovering the initial ellipticity. The re-Gaussianization (‘G’) method is that described in Section 2.2. The linear (‘LIN’) method follows our prescription at the end of Appendix B. The ‘BJ’ method uses BJ02's resolution factor from Appendix C, equation (14), instead of ours, equation (B16). [BJ02 discuss several methods of correcting for PSF dilution, including a matrix inversion method which from theoretical expectations should do better than the resolution factor approach. The resolution factor used here is the same one that is used in the cosmic shear study of Jarvis et al. (2003) and can be implemented on Sloan Digital Sky Survey photo reductions.] The no-correction (‘NC’) method computes ellipticities based on equation (13) and hence is correct only for Gaussian galaxies and PSFs. Included for comparison is the ‘KSB’ method, which uses the KSB95 method with pre-seeing shear polarizability estimated by the method of Luppino & Kaiser (1997). Note that there are several versions of the KSB method in the literature, a comparison of some of these is provided by Erben et al. (2001). (For details of our implementation of KSB95, see Appendix C.)

As our objective is to correct for the dilution and smearing of the galaxy image by the PSF, we do not include noise in these simulations; thus noise-rectification biases such as Kaiser (2000)'s ‘noise bias’ (termed ‘centroid bias’ in BJ02) are not addressed here and a separate correction for these must be made. The bias in the PSF-corrected ellipticity e introduced by noise is generally given by the Taylor expansion:  

(24)

where I_i is the flux in pixel i. Here the leading-order term scales with the significance ν as 1/ν², and the higher-order terms scale as 1/ν³, 1/ν⁴, etc. For galaxies detected at high significance, the dominant noise bias contribution is from the leading 1/ν² term (the shape measurement of low-significance galaxies is unstable and not recommended). Present approaches to the evaluation of the noise bias include analytic solutions for Gaussian objects (Kaiser 2000) and analytic scaling relations, the coefficients of which can be determined through simulations (BJ02). There is also a Monte Carlo approach, valid for any object profile and any PSF method: if we add artificial Gaussian noise to the measured image with covariance ζ〈ΔI_iΔI_j〉 (where the parameter ζ≪ 1), and then take an ensemble-averaged change in measured ellipticity over many realizations of the artificial noise, we find the ensemble-averaged change in measured ellipticity . This allows us to determine 〈Δe〉, and thus correct the measured ellipticity for noise bias: . Note that the statistical variance in e resulting from the Monte Carlo method after N realizations of the artificial noise have been computed is 〈Δe²〉/N, whereas the statistical noise in the galaxy shape measurement is 〈Δe²〉; hence the Monte Carlo convergence rapidly becomes a negligible contributor to the overall statistical uncertainty. Thus we can determine and remove the noise bias effect on any PSF correction method by evaluating equation (24), provided only that the significance ν is large enough to treat the noise as a perturbation. Since the O(1/ν²) noise bias is straightforward to remove and in any case the details of its removal process do not depend strongly on the PSF correction scheme, we have presented the simulations here without noise to illustrate the purely PSF-induced effects.

Since the end result of the shape measurements is to produce a shear estimate, we show the calibration error δγ/γ in the shear estimator as a function of object and PSF parameters. For the adaptive moment methods (no correction, linear, BJ02) that yield ellipticity estimates for individual objects, we proceed as follows: consider the (unweighted) shear estimator:  

(25)

where ê is the measured ellipticity. For objects of a fixed ellipticity this has calibration:  

(26)

Now if we use the transformation law for ellipticities under infinitesimal shear, e_α↦e_α+ 2γ_α− 2e_αe_βγ_β (see our equation 41 in the δ → 0 limit, and note that the infinitesimal shear γ differs from that of the distortion shear δ by a factor of 2), we can azimuthally average to yield  

(27)

This is the equation we use to determine the calibration for the no-correction, linear and BJ02 methods. In the plots (Figs 1–3) we show the error .

Figure 1.

The error in reconstructing the shear for the case of a circular Gaussian PSF and elliptical Gaussian galaxy. We have plotted the shear calibration error δγ/γ. Note that the adaptive moments method is exact for this case. The KSB95 method is not exact but gives good performance (calibration errors of the order of 1 per cent for most galaxies) even when the galaxy is not well resolved.

Figure 2.

The error in reconstructing the shear for the case of a circular Gaussian PSF and elliptical exponential galaxy. We have plotted the shear calibration error δγ/γ. Note that for the well-resolved galaxy (a), both BJ02 and linear methods provide good performance; the error in the linear method grows substantially for the poorly resolved galaxies in (b) and (c). The no-correction method underestimates ellipticities in all cases, with the situation worsening for the poorly resolved galaxies.

Figure 3.

The error in reconstructing the shear for the case of a circular PSF generated by atmospheric turbulence and an elliptical Gaussian galaxy. We have plotted the shear calibration error δγ/γ. For well-resolved (a) and moderately resolved (b) galaxies, the linear method substantially reduces errors when compared with the no-correction method. For poorly resolved galaxies, the linear method overcorrects for PSF non-Gaussianity. The PSF non-Gaussianity correction of BJ02 has the incorrect sign in all three cases. Note that the scale is different from that of Fig. 2.

For the KSB95 method, equation (27) is not applicable because the method does not produce an estimator for the adaptive ellipticity. Rather it produces a non-adaptive ellipticity e^K and a pre-seeing shear polarizability tensor P^γ. The shear is then estimated via  

(28)

For an ensemble of galaxies of the same ellipticity but rotated by different angles, 〈P^γ〉 becomes isotropized and is equal to . We can find 〈e^K〉 by an argument similar to that leading to equation (27); the result is  

(29)

This is the calibration equation that we use for the KSB95 method.

We use several different forms for the PSF in our simulations; these are given in Table 3. The first is the Gaussian form, with PSF g(r) and optical transfer function (OTF) ḡ(k) given by  

(30)

[We define the OTF ḡ(k) to be 2π times the Fourier transform g̃(k) of the PSF so that ḡ(0) = 1.] The Gaussian, while simple to analyse, is not a good approximation to the PSF for real experiments. For ground-based experiments, the resolution is generally limited by atmospheric turbulence in the inertial range; the OTF for this case is given by (see Kaiser 2000, Appendix A):  

(31)

where k₀= 2.92/θ_FWHM and the seeing θ_FWHM is the full-width half-maximum width of the PSF; typically this is of the order of 1 arcsec. [The series for g(r) is convergent for all radii, nevertheless the summation is numerically unstable at k₀r≫ 1 and for computational purposes it is better to take the Fourier transform of ḡ(k).] For space-based experiments, the telescope is usually diffraction-limited and hence the OTF is the autocorrelation function of the aperture. If the aperture is circular with diameter D, this is  

(32)

where k₀= 2πD/λ and λ is the wavelength of observation. We do not pursue an analysis of the effects of chromaticity, but note that (particularly for the diffraction-limited PSF, which has a steeper wavelength dependence than the turbulence-limited PSF) the effect may be important.

Simulation results are shown in the figures for three cases: a well-resolved galaxy (σ²_g/σ²_ƒ= 0.2; we remind the reader that subscript g is the PSF and ƒ denotes the galaxy), an intermediate case (σ²_g/σ²_ƒ= 1.0) and a poorly resolved galaxy (σ²_g/σ²_ƒ= 2.0). (Here is proportional to the area covered by the object.) The simple case of a Gaussian galaxy and a Gaussian PSF is shown in Fig. 1. The adaptive moments methods are exact for this case. The KSB95 method, while not exact, has a relatively small calibration error except for the very elongated galaxies.

In Fig. 2, we test the ability of these methods to handle non-Gaussian galaxies by showing simulation results with a circular Gaussian PSF and an elliptical galaxy with exponential (∝ e) radial profile. The exponential profile has a positive kurtosis β^(ƒ)₂₂≈ 0.17. Note that both the BJ02 and linear methods improve upon the no-correction method; BJ02's performance here is superior to that of the linear method, presumably because BJ02 works to all orders in β^(ƒ)₂₂. Once again, the KSB95 method is accurate at the 1 per cent level except for very eccentric galaxies.

In Fig. 3, we display a simulation using an elliptical Gaussian galaxy and a non-Gaussian PSF. The PSF is the Fourier transform of an e function, as appropriate for observations through a turbulent atmosphere with Kolmogorov scaling (Kaiser 2000). This PSF has a kurtosis of β^(g)₂₂≈ 0.046; it is important to note, however, that the PSF drops to zero slowly at large radii and hence it possesses higher moments, e.g. β^(g)₃₃≈ 0.012 (see Table 2) that are not taken into account by either the BJ02 or linear methods. We can demonstrate that the error from the linear method is due to these higher moments by repeating the simulation using a PSF of the ‘quartic-Gaussian’ form:  

(33)

which has the same kurtosis as the turbulent PSF but no higher-order moments. The results using this form for the PSF are shown in Fig. 4. Finally, we have shown simulations where both the galaxy and PSF are non-Gaussian: one with a turbulent PSF and an exponential galaxy (Fig. 5), one with a turbulent PSF and a de Vaucouleurs profile (Fig. 6) and one with a diffraction-limited (Airy) PSF and an exponential galaxy (Fig. 7). Of the methods described, only the re-Gaussianization method retains high accuracy in these cases. This is primarily because of the sixth-order and higher moments of these PSFs, which are ignored by the BJ02 and fourth-order linear methods. For well-resolved objects, the KSB95 method avoids these problems because its non-adaptive weight function allows it to ‘see’ the turbulence-induced tails and/or diffraction rings in the PSF. For the less well-resolved objects, KSB suffers from loss of higher-moment information, although it is apparent from the simulations that the problem is not as serious as for the BJ02 or linear methods.

Table 2.

PSFs used in the simulations. The adaptive (i.e. iterated until β₁₁= 0) radial 2nth moments β_nn characterize the non-Gaussianity of the PSF. Note the slow convergence of the Laguerre expansion for the diffraction-limited PSF.

Figure 4.

The error in reconstructing the shear for the case of a circular quartic-Gaussian PSF (see equation 33) and an elliptical Gaussian galaxy. We have plotted the shear calibration error δγ/γ. For well-resolved (a) and moderately resolved (b) galaxies, the linear method substantially reduces errors when compared with the no-correction method. For poorly resolved galaxies, the linear method overcorrects for PSF non-Gaussianity. The PSF non-Gaussianity correction of BJ02 has the incorrect sign in all three cases.

Figure 5.

The error in reconstructing the shear for the case of a circular turbulence-limited PSF and an elliptical exponential galaxy. We have plotted the shear calibration error δγ/γ. Note the significant loss of accuracy of all the methods except for re-Gaussianization (G) for poorly resolved galaxies.

Figure 6.

The error in reconstructing the shear for the case of a circular turbulence-limited PSF and an elliptical de Vaucouleurs profile galaxy. We have plotted the shear calibration error δγ/γ. Note the significant loss of accuracy of all the methods except for re-Gaussianization (G) and KSB for poorly resolved galaxies.

Figure 7.

The error in reconstructing the shear for the case of a circular diffraction-limited PSF and an elliptical exponential galaxy. Severe loss of accuracy occurs for most of the methods if the galaxy is poorly resolved (c); fortunately, diffraction-limited PSFs occur only for space-based surveys in which the PSF size is usually small compared with the size of the galaxy, hence the comparatively large errors exhibited in (c) are unlikely to be a problem in practice.

The simulations indicate significant (at least several per cent) calibration errors for all of the currently used methods when the PSF is non-Gaussian and has size comparable to that of the galaxy. Both the KSB95 and BJ02 methods yield an underestimate of the shear for turbulence-limited PSFs, with the error for poorly resolved objects being of the order of 5 per cent for KSB95 and 10–15 per cent for BJ02. The errors have the opposite sign for a diffraction-limited PSF, however, the small size (typically of the order of 0.1 arcsec) of diffraction-limited PSFs helps to reduce the error in this case (cf. Figs 7a–c). The re-Gaussianization method introduced here shows promise in reducing these calibration errors to the 1–2 per cent level even for poorly resolved galaxies. It has the best performance of all the methods studied here. While we have not simulated Refregier & Bacon (2003)'s matrix inversion method, we note for comparison that in realistic Monte Carlo simulations (i.e. including noise and pixelization), Refregier & Bacon (2003) found calibration errors of −3 ± 4 per cent in e₊ and 0 ± 4 per cent in e_× using a turbulent PSF with σ²_g/σ²_ƒ= 0.32. This would thus be comparable to our tests of well-resolved galaxies, for which KSB95 and re-Gaussianization methods are accurate to within 2 per cent. It would be interesting to investigate the behaviour of the matrix inversion method at lower resolution where the other methods except for re-Gaussianization show more significant error.

Shear Selection Bias

Even if the PSF correction is perfect, an estimate of the shear can be biased if the object detection algorithm preferentially selects galaxies elongated in one direction over galaxies elongated in another direction. Kaiser (2000) and BJ02 investigate such a selection bias that results from asymmetries of the PSF and offer methods for correcting it. However, there is another selection bias that exists: suppose that we are looking at a part of the sky where the gravitational shear elongates objects in the x-direction (i.e. γ₊ > 0). Then background galaxies intrinsically (i.e. pre-shear) elongated in the x-direction appear more elliptical than otherwise identical background galaxies elongated in the y-direction (although both have the same area, flux and surface brightness). Many object detection algorithms, such as significance cut-offs, will preferentially detect the circular objects in this situation, thus the catalogue of background galaxies that results will preferentially contain galaxies aligned orthogonal to the gravitational shear. This results in underestimation of the lensing signal. We will call this effect ‘shear selection bias’ to distinguish it from the ‘PSF selection bias’ discussed by Kaiser (2000) and BJ02.

Unlike PSF-related biases, shear selection bias is perfectly (anti)correlated with the shear signal. This means that shear selection bias takes the form of a calibration error, with the fractional error depending on the object detection algorithm and on the underlying population of galaxies. As previously mentioned, calibration errors are difficult to detect using, for example, E/B decomposition of the measured shear field, which is insensitive to a position-independent calibration error.

We organize this section as follows. In Section 3.3, we develop a formalism for relating the galaxy selection probability to the selection bias. In Sections 3.2 and 3.3, we apply this formalism to the problem of evaluating shear selection bias for galaxy detection algorithms based on a statistical significance threshold. We show in Section 3.3 that for this class of algorithms, shear selection bias is an effect of the order of several per cent. Suggestions for overcoming the shear selection bias are offered in Section 3.3.

Effect on measured shear

Here we analyse the effect of shear selection bias on measurements of gravitational shear. Following KSB95, we will use lowercase Roman letters for vector indices and Greek letters for spin-2 tensor (i.e. traceless symmetric matrix) indices. To avoid confusion, we call an object with two tensor indices a ‘hypertensor’.

A shear estimator γ̂ will generally take the form  

(34)

where the sum is over the background galaxies in the region of sky where the shear is to be estimated, ε_J is a measure of the shape of a galaxy, w_J are the weights for each galaxy and T is a calibration factor. In order to have the proper first-order behaviour, we must have  

(35)

where the d subscript indicates that the average is a weighted average taken over the detected background galaxies. Because shears observed in weak-lensing experiments have γ≪ 1, there is typically no need for correct O(γ²) and higher behaviour.

In the BJ02 formalism, the shape ε_J=e_J, where e_J is the ellipticity of the galaxy. In the KSB95 formalism, ε_J is a shear estimator computed according to  

(36)

where ě^β is the circular-weighted ellipticity:  

(37)

with p being the unweighted quadrupole moment of the PSF and the hypertensors P^sm and P^γ are the smear and shear polarizabilities, respectively. The weight function W(r²) is usually taken as a Gaussian. The smear polarizability is chosen to eliminate the PSF bias (i.e. the change in measured galaxy ellipticity due to elongation of the PSF), and the shear polarizability is chosen so as to achieve the correct response, i.e. equation (35) with T = 1:  

(38)

Methods for computing P^γ can be found in terms of Laguerre coefficients in our Appendix C and in terms of real-space variables in Hoekstra et al. (1998)'s Appendix A[these methods are based on KSB95 and Luppino & Kaiser (1997), however, KSB95's shear polarizability contains algebra errors that have been corrected by Hoekstra et al. (1998)].

The effect of shear selection bias can now be analysed. If we have an ensemble of background galaxies gravitationally sheared by a small amount γ, then the shear estimated according to equation (34) is  

(39)

where P indicates the probability of a galaxy being detected. (Detection is a stochastic process due to photon counting noise.) Assuming PSF anisotropies and the camera shear are correctly removed, the terms proportional to ϵ^αP vanish. This yields  

(40)

where we have defined the transfer hypertensor Ξ. Note that if the response of the estimator is correctly calibrated in the absence of selection biases (i.e. equation 35 is satisfied), then we have  

(41)

and  

(42)

in which case Ξ represents the calibration error resulting from shear selection bias. Thus, if we can determine Ξ, or even its average value, the shear selection bias is overcome. Below, we examine the transfer hypertensor for several object detection algorithms and find that Ξ is strongly dependent on the details of this algorithm.

Significance threshold algorithms

The simplest object detection algorithm is the significance threshold algorithm. This method convolves the measured image intensity I(x) with some detection kernel K₁ to yield a convolved intensity J:  

(43)

This quantity has an expectation value  

(44)

where ƒ is the intrinsic intensity of the galaxy, g is the PSF and K is the convolution of K₁ and g. The two-point correlation function of the noise is simply the autocorrelation of the kernel:  

(45)

where n is the noise variance per unit area. We will require K₁ to be normalized so that ξ(0) = 1. Then an object is ‘detected’ if at least one point within it has J(x) ≥ν₀, where ν₀ is the cut-off significance. Note that with our normalization, we can think of the point x as being J(x) standard deviations above the background.

Our principal objective here is to determine the probability P of the object being detected. The definition of ‘detected’, and hence the determination of P, is a delicate issue. If we simply assume some intrinsic intensity ƒ(x) peaked at x= 0 and ask whether J(0) ≥ν₀, we will underestimate P because J is a random field and it is unlikely that 0 will be an exact maximum of J. However, the use of ν=J(0) as a detection statistic will serve as a first approximation to the shear selection bias, and due to the complexity of the problem it is unclear whether substantial gains can be made through a more detailed analytical investigation.

The probability of a galaxy being detected is given by the cumulative Gaussian distribution function:  

(46)

where ν̄ = J̄(0) is the expected significance of the galaxy. Defining the function:  

(47)

we can then write the transfer hypertensor as  

(48)

We can compute the partial derivative using integration by parts:  

(49)

where the mixed hypertensor δ^ij_β has components  

(50)

and we have used the tracelessness of δ_β to eliminate the ∂x_i/∂x^j term. In the special case where K is a Gaussian, K(x) ∝ exp(−1/2|x|²/σ²_K), this reduces to  

(51)

For a single detection kernel K₁, equation (51) describes the transfer hypertensor. However, the single-scale detection algorithms considered here are not optimal for detecting objects substantially larger or smaller than the size of the detection kernel. For this reason, most authors have used multiple detection kernels, which we consider next.

Multiscale detection algorithms

The commonly used KSB95 detection algorithm (see, e.g., KSB95, Hoekstra et al. 1998) detects objects by running a significance threshold algorithm on several smoothing scales, i.e. with several functions K₁. We will analyse here only the case of a Gaussian PSF of standard deviation σ_g and a Gaussian smoothing filter of adjustable standard deviation σ_d. Then we have  

(52)

where σ²_K=σ²_d+σ²_g. An object is detected if it is detected with any of these filters. We must therefore evaluate equation (51) at the most significant scale, i.e. where ν̄ is maximized with respect to σ_d. We determine ν̄ by integration:  

(53)

The maximum significance occurs where , that is where  

(54)

(Remember that σ_K depends implicitly on σ_d.) In the special case of a circular Gaussian object , the left-hand side of equation (54) becomes 2/(σ⁻²_K+σ⁻²_ƒ); this gives a solution to equation (54) of  

(55)

(We ignore the imaginary solution σ²_d=−σ²_g.) This means that if the galaxy is a circular Gaussian, the detection algorithm will return maximum significance with a filter matched to the object size.

Whether or not the galaxy is Gaussian, we can use equation (54) to convert the logarithmic derivative expression, equation (51), into  

(56)

since for near-circular Gaussian objects, ϵ ≈ ě. (This follows because the smear polarizability P^γ is approximately the identity matrix in this case, as can be verified through direct computation using the formulae in Appendix C.) The transfer hypertensor is then (see equation 48)  

(57)

How important is shear selection bias?

We can evaluate, at least at the order-of-magnitude level, the importance of the shear selection bias. Since σ_d is matched to the size of the object, the quantity 1 −σ²_g/σ²_d will be near unity for well-resolved galaxies, decreasing to zero for point-like objects. (That Ξ= 0 for point-like objects is unsurprising because a point source remains a point source when sheared.) The T-factor is shown in Fig. 8. The weighted ellipticities ě are half of the adaptive ellipticities (equation 5) for near-circular Gaussian objects with a matched Gaussian weight function. Thus the ě circular ellipticity has a shape variance of approximately V≈ (0.32/2)²≈ 0.026, where 0.32 is the rms ellipticity of galaxies (McKay et al. 2001). (This parameter will depend somewhat on the sample examined, so the value of 0.32 should be thought of as a rough guide rather than a law of nature.) It thus appears that the shear selection bias will be of the order of several per cent for well-resolved galaxies, falling to zero in the limit of a poorly resolved object.

Figure 8.

The T-factor T(ν̄; ν₀) = d 1n P/d 1n ν̄. This factor represents the (logarithmic) increase in probability of detection P if we increase the significance of the galaxy. The shear selection bias increases in proportion to an appropriate average of T (see equation 48). Note that T→ 0 as ν̄ → ∞; this means that there is no shear selection bias for objects far above threshold, as expected.

A more quantitative approach to the problem is possible in the special case where the galaxy significance distribution is a power law and the shear estimator is constructed from an unweighted distribution of galaxies (i.e. w=constant). We suppose that the galaxy counts scale as . The exponent λ is determined by noting that ν̄ is a signal-to-noise ratio and hence is proportional to the total flux of the galaxy, F, divided by the square root of its post-PSF-convolved area A; the number–magnitude relation is roughly N∝F⁻¹. If the galaxies are poorly resolved, A is determined by the PSF and hence is constant, in which case λ= 1. For well-resolved, nearby galaxies (redshift z≪ 1) we can assume a Euclidean universe, in which the usual inverse square laws F∝A∝ 1/r² apply; this gives , since the number of objects scales as r³. This gives λ= 3. Thus we expect that λ will probably lie in the range 1 ≤λ≤ 3, with the distant galaxies of interest to cosmic shear measurements probably having λ closer to 1.

Next, we evaluate the average of equation (57). We define a shear selection bias resolution factor , which is 1 for well-resolved and 0 for poorly resolved objects. We also use the fact that . The average transfer hypertensor is  

(58)

where R_eff is a weighted (by T) average of the shear selection bias resolution factor R over the detected galaxies: . Since T is a rapidly decreasing function of ν̄ (see Fig. 8), R_eff should be evaluated for galaxies near the detection threshold.

At this point we note that the integrals in the numerator and denominator of equation (58) diverge in the ν̄ → 0 limit because P(0) > 0; see equation (46). However, this divergence must be unphysical since real source galaxy catalogues do not contain an infinity of faint objects. (The origin of the problem is that in deriving equation 46, we neglected the fact that objects overlap with each other.) We will thus remove this divergence by supposing that at very small ν̄, the detection probability P vanishes sufficiently rapidly for the integrals in equation (58) to converge; this yields (using integration by parts):  

(59)

Thus for cosmic shear measurements where galaxies are selected via a statistical significance cut, where V= 0.026 and λ≈ 1, the shear selection bias causes an underestimate of the shear by 2.6R_eff per cent. (It vanishes for the poorly resolved galaxies with R_eff≈ 0.) This factor would be increased to roughly 8 per cent in the case of λ= 3, R_eff= 1, i.e. of low-redshift source galaxies for which the Euclidean source counts are more nearly valid and the near-threshold objects are well resolved. Note that this result is independent of the significance threshold, at least when a power-law form for the galaxy count is assumed. Indeed, because it does not depend on the particular form for P, this result is still valid if only those source galaxies are selected that appear in multiple images of the same region of sky (this is the case in most current lensing surveys). It is also important to understand that shear selection bias is small not because the shear is small, but because the background galaxies are so nearly circular. If the galaxies are more elliptical in deeper surveys than in SDSS used here the bias will be larger. Note also that the mean rms ellipticity, λ and R_eff are all likely to be redshift dependent and so will be the bias.

Correcting for shear selection bias

While the above analysis serves to illustrate the order of magnitude of the shear selection bias problem, it may not be adequate to provide an analytical correction due to the large number of assumptions and approximations made. There are also contributions to the selection bias other than significance thresholds; for example, it is generally necessary to remove overlapping objects from the background galaxy catalogue, but the overlap probability is a function of galaxy shape. There are four major approaches that could be taken to overcome shear selection bias.

• (1)

  The bias could simply be ignored.

• (2)

  The galaxy detection algorithm could be designed so as to eliminate, or at least substantially reduce, the bias.

• (3)

  The transfer hypertensor Ξ could be computed and used as a correction to the shear polarizability or resolution factor.

• (4)

  The bias could be determined from simulations.

Of these, Method 1 is appropriate when calibration errors of the order of several per cent are unimportant, for example in typical galaxy–galaxy lensing studies where the statistical error dominates. However, it is clearly inapplicable to future cosmic shear surveys that aim for much higher accuracy. Method 2 requires a detection statistic that does not suffer from shear selection bias, i.e. one that is independent of the shear. Examples would be the magnitude of the galaxy in a given band and the PSF-deconvolved area. Unfortunately, these are not truly independent of the shear when noise is taken into account, and the ‘overlapping object’ effect described above is not addressed. These effects are difficult to compute analytically, hence it is not clear whether we can verify that the residual errors in Method 2 are ‘small enough’ without resorting to simulation.

Method 3 requires that we compute for each galaxy in our survey, then find Ξ^α_β=ε^αE_β, and compute the weighted average 〈Ξ^α_β〉_wd to determine the calibration error due to shear selection bias. The problem is thus one of finding E_β; unfortunately, E_β can only be determined from the intrinsic intensity ƒ of the galaxy, not the observed image. Thus correction of cosmic shear measurements based on the determination of Ξ would require a careful analysis to ensure that noise and PSF effects do not invalidate the calculation. Method 4, the simulation, is straightforward albeit time consuming. Its advantage is that it can accommodate arbitrarily complicated detection algorithms. The main weakness is that the underlying distribution of galaxies in the simulation may not be correct; however, the shear selection bias is only an effect of the order of several per cent anyway, so determination of the bias to many significant figures is unnecessary.

Conclusion

In this paper we have addressed several possible biases present in the data reduction of cosmic shear measurements. We focus on effects that affect the overall normalization only and cannot be detected using the usual methods of finding systematics, such as the tests on stars and E/B decomposition. We pay particular attention to the deviations from the elliptical Gaussian profile of the galaxy or the PSF, which is often the starting assumption of the data reduction methods. The results in Figs 1–7 show that while for a Gaussian PSF most of the methods are successful in reconstructing the shear, the case of a more realistic non-Gaussian PSF is significantly more problematic except for the re-Gaussianization method introduced in this paper. This suggests that in existing surveys the shear calibration may be inaccurate at the 10 per cent level. Some methods, such as the current implementation of BJ02 used in Jarvis et al. (2003) and McKay et al. (2001), are particularly susceptible to this and can underestimate the linear amplitude of fluctuations by 10 per cent, while those based on KSB95 (e.g. Bacon et al. 2000; Van Waerbeke et al. 2000; Rhodes et al. 2001; Brown et al. 2003; Hoekstra et al. 2002) have a negative bias of up to 5 per cent. It is interesting to note that Jarvis et al. (2003) find a significantly (10–20 per cent) lower value for σ₈ than other surveys. Calibrations errors of the type discussed here can explain part of the difference, although other factors such as statistical errors, uncertainties in the redshift distribution of background galaxies, non-linear evolution corrections of the dark matter power spectrum and residual systematics (which are known to be present because their effect on B-type correlation is seen in all the surveys that have looked for it) also play a role. At the moment, residual systematics due to the PSF correction can be as great as 10 per cent in σ₈ and are comparable to the statistical errors in current samples.

There are other biases that one must worry about when an accuracy of the order of 1 per cent is required. One such effect is shear selection bias, which causes the galaxy selection to be correlated with the underlying shear and with the galaxy flux. This error is of the order of a few per cent, but can reach up to roughly 8 per cent in some limits. Since the shear selection bias preferentially selects galaxies with shear orthogonal to cosmic shear the bias translates into an overall underestimate of the amplitude of fluctuations. Moreover, since for a flux-limited survey the smaller and fainter galaxies are at a higher redshift the error in the amplitude may be correlated with redshift and can give an incorrect growth factor if not properly accounted for. Another effect that can create this is if the mean rms ellipticity of the galaxies changes with redshift.

Looking further ahead, weak lensing has the potential to break many of the degeneracies present in the cosmological parameters even after the cosmic microwave background has given us all of its information, especially if one can use the redshift information on background galaxies to extract the information on the growth factor (Hu 1999). It is an especially powerful probe of dark energy and massive neutrinos, since they both affect the growth factor at low redshifts probed by weak lensing. While the statistical precision of some of the proposed experiments such as LSST is impressive, the systematic effects must also be controlled at the 1 per cent level. Current reduction methods are not sufficient for this purpose, but improved methods such as our re-Gaussianization method (which can be further generalized to include higher-order terms) should be up to the task. Similarly, shear selection bias and other biases should be carefully simulated to estimate their importance as a function of galaxy detection significance, PSF, etc. One notes that realistic simulations would allow us to determine both the PSF dilution calibration and the shear selection bias, and indeed shear selection bias is implicitly taken into account in simulated observations such as those of Kaiser (2000), Erben et al. (2001) and Bacon et al. (2001). Future experiments will require that this line of work be extended, for example to include measured PSFs from the telescope and the correct distribution of galaxy shapes and fluxes as a function of redshift.

It is clear that the quality of the data of the upcoming measurements will exceed the current generation of reduction methods. However, this is not a cause for alarm, since none of the problems discussed here is fundamental in nature. We have proposed solutions for some of the problems, while others will require more detailed modelling using realistic simulations of the galaxies and their image processing through the atmosphere (if applicable), telescope, detector and software. Given enough attention to these details weak lensing should fulfil its promise in the era of high-precision cosmology.

Acknowledgments

The authors acknowledge useful discussions with Gary Bernstein, Mike Jarvis, Robert Lupton, Michael Strauss, Joseph Hennawi and Nikhil Padmanabhan. CH is supported through the NASA Graduate Student Researchers Program, grant NASA GSRP-02-OSS-079. This work was supported by NASA NAG5-1993 and NSF CAREER. US is a fellow of the David and Lucille Packard Foundation and the Alfred P. Sloan Foundation.

Appendix

Appendix A: Convolution with Elliptical Laguerre Expansion

Here we examine the problem of convolving two functions ƒ and g given in terms of their Laguerre expansions:  

(A1)

The convolution is  

(A2)

we will consider here the case where M_h=M_ƒ+M_g. This problem is easiest in Fourier space using the convolution theorem ; using equation (10) and writing the result in terms of yields  

(A3)

We can simplify this by defining the matrices A^(ƒ)=M^1/2_ƒM^−1/2_h and similarly A^(g). Then we use the result that any 2 × 2 matrix, whether symmetric or not, can be written as a product of a rotation matrix, a diagonal matrix and another rotation matrix (we will call this the RDR decomposition):  

(A4)

where R_θ represents a rotation by angle θ, and C̄, D̄, and the rotation angles α and χ can found explicitly from the rectangular-to-polar conversions:  

(A5)

This decomposition is unique up to several discrete degeneracies: (i) interchange of C̄ and D̄ with a phase shift of α by +π/2 and χ by −π/2; (ii) sign change of C̄ and D̄ with phase shift of α and χ together by π/2; (iii) phase shift of α and χ together by π; and (iv) phase shift of α and χ independently by multiples of 2π. (In the special case that C̄ = ± D̄, there is an additional continuous degeneracy in that only α±χ, but neither angle is determined individually.) If we decompose A^(ƒ) and A^(g), then the condition M_h=M_ƒ+M_g yields that A^(ƒ)TA^(ƒ)+A^(g)TA^(g) is the identity; this can be expressed in matrix form as  

(A6)

The trace of this equation yields ; the traceless-symmetric part yields

 

(A7)

which implies that either and , or that χ_ƒ and χ_g differ by a multiple of π/2; we may, however, assume the latter case, since in the former case the aforementioned continuous degeneracy allows us to choose χ_ƒ=χ_g. Indeed, the discrete degeneracy allows us to set χ_ƒ=χ_g(≡χ). From this we derive . Finally, the positive determinant of A^(ƒ) and A^(g) implies that C̄_ƒD̄_ƒ and C̄_gD̄_g are positive; the discrete degeneracy then allows us to choose C̄_ƒ, C̄_g, D̄_ƒ and D̄_g positive.

Explicit expressions for , and D̄_g in terms of M_ƒ,g,h may be obtained by directly finding the square roots of these matrices, thereby determining A^(ƒ,g), and then substituting these matrices into equation (A5). In order to construct the As, it is necessary to determine the square root of a symmetric matrix. This is given by  

(A8)

where D=M_xxM_yy−M²_xy.

With this aside completed, we return to the determination of c_pq. Substituting the decomposition (equation A4) into equation (A3), and setting yields  

(A9)

where we have used the result that Λ_pq(R_θk) = e^−imθΛ_pq(k). Then equation (A9) is a linear equation for the c_pq, thus  

(A10)

where the coefficients satisfy  

(A11)

can be calculated by noting that the right-hand side of equation (A11) is a polynomial of the order of p′+q′+p′′+q′′ in the two variables U and V. Since Λ_pq(U, V) with p+q≤N form a basis for the (N+ 1)(N+ 2)/2-dimensional vector space of polynomials in U and V of degree ≤N, the sum on the left-hand side of equation (A11) may be cut-off at p+q≤N. Equating coefficients of U^mVⁿ on both sides then results in a finite system of linear equations for , which is easily solved using standard linear algebra techniques. An alternative route is to exploit the orthonormality relations among the Laguerre basis functions to convert equation (A11) into  

(A12)

this ratio of Gaussian integrals can be computed using Wick's theorem with the Feynman rules 〈U〉=〈V〉=〈UV〉= 0 and 〈U²〉=〈V²〉 = 1/2.

The Gaussian-integral method allows us to construct an explicit formula for in terms of finite sums. To do this, we expand the Λ:  

(A13)

where  

(A14)

and s= (p+q−μ−ν)/2. The sum in equation (A14) is over all values of t such that s, μ−t, t, and q−s−t are all non-negative. The ‘selection rules’ for X are as follows: X_pqμν≠ 0 only if μ+ν is equal to p+q, or less than p+q by an even integer. Having constructed X, we evaluate the integral in equation (A12) as  

(A15)

We have computed several of the lower-order (up to degree six) Z coefficients for the simulations in this paper; it is found that for these moments, it is practical to perform either the linear system solution (equation A11) or the summation method (equation A15) using symbolic manipulation software, and then write ‘hard-wired’ C code to numerically evaluate the resulting expressions. Because the explicit polynomial form of the Λ functions (see equation A14) contains many positive and negative terms for large p and q, it is likely that round-off errors would be an issue for our computational approach if we were to consider moments of sufficiently high order. It is unlikely that moments of very high order would be useful in weak-lensing studies; nevertheless, if this turned out to be the case, or if the high-order Z coefficients were required for another application, it would be necessary to develop a better algorithm (likely involving a recursion relation) for computing them.

Appendix

Appendix B: Linear Psf Correction

We can derive a relation between e^(ƒ) and e^(I), valid to first order in the departure of the galaxy and PSF from Gaussianity as follows. Suppose that the adaptive covariance of the intrinsic galaxy image is M_ƒ. Then the observed image I(x) is the convolution of ƒ with the PSF g, and it can thus be expressed as a Laguerre expansion with covariance M_h=M_ƒ+M_g. Using the convolution results of Appendix A, we find  

(B1)

where  

(B2)

and b^(ƒ)_p′q′ and b^(g)_p″q″ are the Laguerre coefficients of ƒ and g around M_ƒ and M_g, respectively. (The rotation angles χ, α_ƒ, and α_g and diagonal elements C̄_ƒ and D̄_ƒ can be calculated using the approximation M_ƒ≈M_I−M_g since corrections to this result will appear as second-order corrections to the final moments.) Since we are working to first-order in the non-Gaussianity of the galaxy and PSF, we may simplify equation (B2) to give (for pq > 0):  

(B3)

where we have worked to first order in β and used the result that Z^00,00₀₀= 1.

Next, we wish to relate the c_pq to β^(I)_pq. This can be done by use of the generators for point transformations provided by BJ02. To first order in the β and in M_I−M_h, and noting that β₁₁=β₂₂, we obtain (see BJ02 equations 6-45 and 6-46):  

(B4)

where μ is the dilation and η is the shear required to transform from a covariance M_I to M_h. Explicitly,  

(B5)

where ℜ and ℑ denote real and imaginary parts, respectively.

Assuming that we know the PSF, β^(g)_p″q″ can be determined. Furthermore, we know that β^(ƒ)₁₁=β^(ƒ)₀₂=β^(ƒ)₂₀= 0. Now truncate equation (B3) with a finite set L of Laguerre modes, i.e. we sum only over p′q′∈L. We then approximate the higher-order moments c_pq/c₀₀≈β^(I)_pq for in accordance with equation (B4). This gives us a system of linear equations that can be solved for the unknown β^(ƒ)_pq:  

(B6)

for pqϵL̄. The β thus determined can be substituted into equation (B3) to yield c₁₁/c₀₀ and c₀₂/c₀₀, which then yield M_h via equations (B4) and (B5). Finally, we can determine M_ƒ=M_h−M_g and then determine the intrinsic ellipticity from equation (5).

The above procedure is somewhat involved, and we illustrate it with a simple example: that of L̄ = {22}. That is, we will correct for the radial fourth moment or kurtosis of the PSF and intrinsic galaxy image. Obviously, in ignoring the higher-order radial and angular moments we are making an approximation, however, the b₂₂ moment is present for non-Gaussian galaxy profiles such as the exponential or de Vaucouleurs profiles (and is positive in both cases). As we show in Section 2.2, the higher-order moments can be important for some PSFs. The angular moments b_pq for p≠q are not present for purely radial profiles, i.e. profiles for which the isophotes are concentric ellipses that share principal axes and axis ratios. If the angular moments, or the higher-order radial moments b₃₃, b₄₄, etc. are important, these would need to be included in L̄.

Because of the shear and rotation invariance of the adaptive moments method, we can assume without loss of generality that the PSF has e^(g)= 0, and that the measured object has e^(I)_×= 0 and e^(I)=e^(I)₊≥ 0, i.e. that it is elongated along the x-axis. (Alternatively, if a rounding kernel has been applied to circularize the PSF, then we already have e^(g)= 0, and no application of shear invariance is necessary.) We then have  

(B7)

Then to zeroth order we determine  

(B8)

Following the procedure outlined in Appendix A, we compute the RDR decomposition of M^1/2_ƒM^−1/2_I and M^1/2_gM^−1/2_I; this gives α_ƒ=α_g=χ= 0, and  

(B9)

Next, we must construct the system of equations, (B6), which requires knowledge of the Z coefficients. The relevant Z coefficients for the intrinsic image are  

(B10)

and those for the PSF are  

(B11)

From these we construct the linear system, equation (B6); in this case it is only a 1 × 1 system and hence is easy to solve:  

(B12)

Then we can compute the correction M_I ↦ M_h due to non-Gaussianity using equations (B3) and (B4):  

(B13)

(Note that in this case, η is real. This is not necessarily the case if we use m≠ 0 modes in our basis L̄.) Finally, we compute M)_h using equation (B5) and find that  

(B14)

The ellipticity e^(ƒ) of the intrinsic galaxy image are easily computable from M_ƒ using the covariance-ellipticity relations (equation 5):  

(B15)

This is equivalent to BJ02's method, equation (14), except that the resolution factor is now  

(B16)

If a shear has been applied to circularize the PSF, the inverse shear must now be applied to yield the physical e^(ƒ). The transformation of ellipticities e↦e′ under a shear (δ₊, δ_×) (i.e. a shear that maps a circle to an ellipse of ellipticity δ) is given by (see, e.g., BJ02 equation 2–13)  

(B17)

where δ²=δ²₊+δ²_×; the transformation of traces can be determined from conservation of area, .

[Warning: the ‘shear composition’ operator e′=δ⊕e as defined by BJ02 and in equation (B17) does not form a group because it is not associative. This problem arises because the shear e followed by the shear δ does not have the same effect as the shear e′ (they differ by an initial rotation). For this reason, it is probably better to think of ⊕ as representing the transformation law for ellipticities under shear, rather than as a composition operator.]

We can thus list the steps involved in PSF correction using the moments L̄ = {22}. If the PSF has been circularized using a rounding kernel, the steps marked with a ⋆ can be skipped.

• (i)

  Compute the adaptive covariances M_I and M_g for the galaxy and PSF, and the corresponding radial fourth moments β^(I)₂₂ and β^(g)₂₂.

• (ii)

  ⋆ Apply a shear to circularize the PSF, i.e. use equation (B17) with δ=−e^(g). The radial fourth moments are unchanged by this shear.

• (iii)

  Apply a rotation so that the measured galaxy has e^(I)₊ > 0 and e^(I)_×= 0.

• (iv)

  Compute the dilation μ and shear η in accordance with equations (B7)–(B13).

• (v)

  Compute the intrinsic ellipticity of the galaxy in the sheared rotated coordinate system using equation (B15).

• (vi)

  Undo the rotation.

• (vii)

  ⋆ Shear to return the ellipticity to the unsheared unrotated coordinate system. The shear is undone using equation (B17) with δ=e^(g), i.e. the exact opposite shear that was applied initially.

Appendix

Appendix C: The Ksb Method in Laguerre Coefficients

Here we describe our implementation of the KSB95 shear estimator, updated to include the pre-seeing shear polarizability (Luppino & Kaiser 1997). For simplicity, we will assume a Gaussian weight function. We begin by computing the Laguerre expansion (see equation 11) of the measured object I and PSF g to second order with a circular covariance matrix:  

(C1)

where σ is the width of the weight function. Then KSB's ellipticity e^K (which, in general, is distinct from the adaptive ellipticity!) is given by  

(C2)

Then the shear polarizability hypertensor is given by  

(C3)

where  

(C4)

and  

(C5)

The KSB/LK method also requires the smear polarizability hypertensor, given by  

(C6)

where  

(C7)

and  

(C8)

In these equations, the X contributions to the polarizabilities are due to the effect of shear on the quadrupole moments b₂₀ whereas the e_αe^sh,sm_β contributions are due to the effect of the shear on the weighted trace T_w. When the X hypertensors are expressed in terms of the Laguerre moments, the decomposition into spin-0 (scalar) and spin-4 (traceless-symmetric) components is manifest.

Luppino & Kaiser (1997) then defines a pre-seeing shear polarizability hypertensor to take account of the finite PSF size:  

(C9)

Our only remaining issue is the choice of the width σ of the weight function. We use a weight function scale set by the requirement b₁₁= 0, i.e. we use an adaptive circular variance. This is appropriate, at least for the simple case of a circular Gaussian object, for maximum significance detection with a Gaussian filter of variable width (see the discussion in Section 3.3). Following Hoekstra et al. (1998), we use the same σ² for the PSF as for the object.