{'fraction_non_alphanumeric': 0.034285676441760046, 'fraction_numerical': 0.012788418236244669, 'mean_word_length': 4.415184335102568, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Shapes and Shears, Stars and Smears: Optimal Measurements for Weak Lensing

We present the theoretical and analytical bases of optimal techniques to measure weak gravitational shear from images of galaxies. We first characterize the geometric space of shears and ellipticity, then use this geometric interpretation to analyse images. The steps of this analysis include: measurement of object shapes on images, combining measurements of a given galaxy on different images, estimating the underlying shear from an ensemble of galaxy shapes, and compensating for the systematic effects of image distortion, bias from PSF asymmetries, and `"dilution"of the signal by the seeing. These methods minimize the ellipticity measurement noise, provide calculable shear uncertainty estimates, and allow removal of systematic contamination by PSF effects to arbitrary precision. Galaxy images and PSFs are decomposed into a family of orthogonal 2d Gaussian-based functions, making the PSF correction and shape measurement relatively straightforward and computationally efficient. We also discuss sources of noise-induced bias in weak lensing measurements and provide a solution for these and previously identified biases.


Introduction
Gravitational lensing is a powerful tool for studying the distribution of matter in the Universe, because photons are deflected by all forms of matter, regardless of luminosity, composition or dynamical state. Dramatic manifestations of lensing-multiple images, Einstein rings, and giant arcs, so-called strong lensing-provide much information on the highest overdensities in the Universe, namely rich galaxy clusters, cores of individual galaxies, or collapsed objects. To characterize the more typical mass structures, or those without a fortuitously aligned bright background source, we may use weak gravitational lensing, in which we analyze the low-order distortions of the ubiquitous background galaxies in order to infer the mass distribution. Weak gravitational lensing signals are extraordinarily subtle, even by astronomical standards: one seeks a shear (or magnification) of the galaxy images amounting to a few percent at most, more typically 0.2-1% in current studies. Because the undistorted image is not observable, the lensing distortions must be detected as a perturbation to the intrinsic distribution in galaxy shapes (or sizes), which have variation of 30% or more, giving a signal-to-noise ratio (S/N ) of ∼ 1/30 from observation of a single galaxy. Hence a very large number of galaxies must be observed before the weak lensing becomes detectable over this intrinsic shape noise. Weak lensing analyses could not even be attempted until automated means of measuring very large numbers of galaxy shapes became available (Valdes, Tyson, & Jarvis 1983;Tyson et al. 1984). Furthermore, optical and atmospheric distortions in a typical sky image cause coherent shape (and size) distortions that can masquerade as a lensing signal. Such systematic errors are 1-10% in a typical image, up to 50 times larger than the weak lensing signals. A means to remove this contamination is crucial; the necessary analyses can only be conducted with well-calibrated, linear detectors.
Successful detection of a weak lensing signal did not occur until CCD images of sufficient depth and field were available (Tyson, Valdes, & Wenk 1990), and early detections were of the ≈ 10% shears that are found in the inner regions of rich clusters of galaxies (Fahlman et al. 1994;Bonnet, Mellier, & Fort 1994;Smail et al. 1995). In regions of strong shear, the S/N is sufficiently high that a map of the lensing mass can be created (Kaiser, Squires, & Broadhurst 1995). Mellier (1999) includes a review of results from cluster-lensing studies.
With the increase in collecting area of CCD imagers, sufficient background galaxies can be measured to allow convincing detection of smaller shear signals around weaker overdensities: around individual weak clusters  or collections of galaxy groups (Hoekstra et al. 2001); around individual galaxies Smith et al. 2000;Wilson et al. 2001b). Most dramatically, lensing signals on random lines of sight, caused by the background matter fluctuation spectrum, have now been detected and are one tool for "precision cosmology" (Wittman et al. 2000;Van Waerbeke et al. 2000;Bacon et al. 2000;Wilson et al. 2001a). As technology has advanced, weaker and weaker shears have become detectable under the shape noise, sometimes as small as a few tenths of a percent (e.g. Fischer et al. 2000, Jarvis et al., in prep.). As a consequence, the demands for rejection of systematic errors have become more stringent. In many current weak lensing publications, it is clear that the uncorrected systematic effects are only slightly smaller than the signals under study. It is therefore fair to say that, at present, it is the analysis techniques, rather than the ability to collect galaxy images, that bar the way to higher precision in many weak-lensing studies. This paper describes the techniques for extraction of weak-lensing signals from imaging data, which we have developed over the past few years to meet these increasing demands. As described below, our efforts focus on the shear rather than the magnification of the galaxy images by the lens, and hence we are measuring galaxy ellipticities. The desiderata for a weak-lensing methodology include: 1. Shapes of individual galaxy images are determined with the highest possible accuracy in the presence of measurement noise on the image.
2. Each shape measurement should have a known error distribution.
3. Individual galaxy shapes should be combined to yield an estimate of the underlying lens shear with maximal S/N .
4. The shear estimator should have an error level and a calibration that can be derived directly from the data, without recourse to Monte-Carlo simulations.
5. The galaxy shapes should be corrected for the systematic biases due to the point-spread function (PSF), to arbitrary precision.
6. The scheme must allow for a PSF that varies continuously across the image and is different in each exposure.
Given the intrinsic floor on weak lensing accuracy because of shape noise, one might ask why we should expend much effort on goal (1), which is to minimize the effects of measurement noisenormally, we consider that once the ellipticity measurement noise σ e is ≪ σ SN ≡ e 2 /2 1/2 ≈ 0.3, further gains do not increase the shear estimation accuracy-the error σ δ on the lensing distortion δ will just become σ SN / √ N , with N the number of measured galaxies. We note first that the sky density of galaxies scales with apparent magnitude m as 10 αm with 0.3 α 0.4. If we can cut the shape measurement error for a given image noise level, then we can either use fainter galaxies in our lensing measurement (increasing N ), or cut the required exposure time. Second, note that convolution with a PSF suppresses the measured lensing signal and the intrinsic shape noise. Hence the level to which we aim to reduce σ e must, for poorly resolved galaxies, be well below the canonical 0.3. Thirdly, we will see in §4.2 that it may be possible to measure the shear to an accuracy much better than σ SN / √ N , in cases where the distribution of intrinsic galaxy shapes in the ellipticity plane has a cusp or pole at e = 0. In simple cases such as a population of circular disk galaxies, the accuracy to which we can measure the applied shear can increase without limit as the measurement noise is decreased.
The need for traceable uncertainties is also critical as weak lensing is used to measure the power spectrum of mass fluctuations in the Universe. In this application, the measurement uncertainties (including shape noise) contribute to the power spectrum and must be accurately estimated and subtracted to reveal the true cosmic power spectrum. Of course an accurate calibration is also necessary for most applications to precision cosmology; if one must rely on simulated data for the calibration, there is always the danger that the simulations do not properly incorporate some aspect of the real world.
Finally, the need for removal of the systematic PSF ellipticities to arbitrary precision is extremely strong. In the course of this paper we will try to describe other approaches to the problem and compare to our own. The methods in most common current use (e.g. Kaiser, Squires, & Broadhurst (1995), KSB) are formally valid only certain special cases of PSF. While heuristic adjustment and testing has demonstrated that the method works to nominal accuracy in more general cases, the absence of a generally valid method is troubling. A formally exact PSF correction scheme has been put forward by Kaiser (2000)[K00], which is based upon a Fourier-domain calculation of the effects of shear and of PSF convolution. Our approach will be to decompose the image and the PSF into a vector over orthogonal polynomials, and treat the deconvolution as a matrix operation carried to desired order. A very similar approach has been independently put forth by Refregier (2001). This is a longwinded paper, likely to be read in detail only by practitioners of weak lensing. A more casual reading will be beneficial to those who wish to understand the methods and limitations of past and future weak lensing analyses. Some of the techniques we develop may be useful beyond the weak lensing analysis, for example our deconvolution method ( §6.3.5) and the methods for rapid convolution with spatially varying kernels ( §7). As discussed by Refregier (2001), our orthogonalfunction decomposition can be a useful means for compression of galaxy images.
The paper outline is as follows: the following section describes the mathematical space occupied by ellipticities and shears. Understanding the geometry of this space makes it easier to see how our (and other) measurement schemes work. In §3, we describe a scheme which uses our geometrical conceptualization of ellipticity to produce measurements with maximal S/N in images with infinitesimal PSF; a formula for the resultant uncertainty in each ellipticity is also derived. Next, §4 discusses several schemes for combining shape measurements of a given galaxy from different exposures and/or filter bands, to obtain the shape estimate that again offers the best possible S/N and a closed-form error estimate. §5 describes the means to combine shape estimates from different galaxies to form an optimal estimate of the underlying lensing shear. In the absence of measurement noise, this takes a simple closed form; in the presence of measurement noise, some approximations must be made to obtain a closed form for the calibration and error of the shear, and hence we do not fully satisfy goal (4) above. §6 is a very extensive discussion of the effects of the PSF on the image, other approaches to the problem, and our method for optimal extraction of the intrinsic shape in this case. In this section we introduce the Laguerre decomposition technique. §7 uses the Laguerre formalism to construct convolution kernels that can add symmetries to the PSF of an image; this is one means of removing the ellipticity biases due to the PSF. §8 discusses two very important effects that can give rise to biased lensing measurements even when a perfect deconvolution for PSF effects is possible. It is likely that these biases are present in all previously published data.
Finally, §9 puts together all of the methods developed in the paper in a flowchart form describing how raw image data are converted into optimized, calibrated lensing shear data. We reserve for a succeeding paper  the detailed discussion of the code that implements these methods, and a verification of its performance on real and simulated data. In Appendices to this paper, we present the formulae for invoking various transformations on the Laguerre-decomposition representation of an image, and derive some approximate PSF-correction formulae that were used for the analyses of Smith et al. (2000) but which are superseded by the full Laguerre methodology.

Linear Approximation to Lensing
The goal of weak gravitational lensing studies is to infer a distant gravitational potential via the distortions that the potential's deflection of light imparts upon the population of galaxies in the background. The lensing is fully characterized by the map u(x) from the observed angular position x to the source angular position u. The surface brightness observed at x is equal to that which would have been observed at u in the absence of the lens. For an individual background galaxy that is not near a lensing caustic, the map can be accurately approximated by a Taylor expansion (2-1) The displacement u 0 carries no information (unless the source is multiply imaged) because the source plane is unobservable. The 2 × 2 amplification matrix has a unique decomposition of the form dx du = µSR, where R is an orthogonal matrix (rotation); S is a symmetric matrix with unit determinant (shear); and µ is a scalar magnification.
The rotation R is not useful for lensing studies because the unlensed orientation is not known, and the ensemble of background galaxies should be isotropic and hence any collective statistic should be unchanged by rotation. Furthermore the rotation is absent in the limits of single-screen or weak lensing.
The magnification µ increases the angular size by µ and the galaxy flux by µ 2 . While the unlensed quantities are not observable, the magnification is still detectable because the mean flux and size of the population will shift. The magnification also reduces the sky-plane density of sources by µ 2 . The magnification thus modulates the number vs magnitude relations for a given class of background galaxies, in a manner which depends upon the size/magnitude/redshift distribution of the original population.
The shear S has two degrees of freedom, corresponding to the ellipticity and position angle imparted on a circular source galaxy. For weak lensing this shear is undetectable on a single galaxy because the unlensed shape is not necessarily circular and is not observed. The collective distribution of galaxy shapes is assumed to be intrinsically isotropic, and the applied shear breaks this symmetry, rendering it detectable and measurable.
Both the shear and the magnification thus produce measurable effects on the ensemble of galaxies and can in theory be used to quantify the potential. Shear measurements have been used for numerous quantitative studies, but magnification methods still yield at best marginal detections (Dye et al. 2001). There are several factors that favor the shear method: first, the two effects of magnification (increased flux and reduced areal density) push the counts of background sources in opposite directions, weakening the signal. More importantly, the shear is manifested as a variation in the mean orientation of galaxy shapes, and this mean is zero in the absence of lensing; the magnification signal is a modulation of N (m) or some other non-zero quantity. It is always far easier to measure a small change from zero than a small change in a non-zero quantity. For example, exploitation of the magnification effect in the weak-lensing regime would require absolute photometry to much better than 0.01 mag accuracy. Magnification measurements, on the other hand, give a direct measure of the projected mass, whereas mass reconstructions from shear data are degenerate under the addition of a constant-density mass sheet. Hence magnification data are very useful when there is no a priori means of determining the mean mass overdensity in the image.
Henceforth we will ignore the magnification effect and describe how to optimally measure the shear S.

Diagonal Shears
The simplest shear matrix is a small perturbation aligned with the coordinate axes: The effect of this transformation upon a circular source-plane object is to induce an elongation along the x axis, creating an elliptical image with axis ratio q ≡ b/a = 1 − η. We can use this matrix as a generator for the full family of diagonal shear matrices with arbitrary η to obtain The set of diagonal shear matrices forms a group under simple matrix multiplication. The operation is commutative, and clearly corresponds to simple addition of the η parameters: For this reason we will call η the conformal shear and will find it a useful parameterization of shear.
Other common parameterizations of shear include the axis ratio q, the distortion δ (Miralda-Escudé 1991), and the reduced shear g = γ/(1 − κ) (Schneider & Seitz 1995), which are related to η via (2-8) Bonnet & Mellier (1995) define a further set of shear parameterizations, also easily expressed in terms of η: Note that for η ≪ 1, the parameters η, e, ε, τ , and the distortion δ are all equal. Note also that most other author's formulations of shear do not define the matrix S to have unit determinant, so do not form a group.

General Shear
A general (non-diagonal) shear matrix can be decomposed into a diagonal shear and rotations as (2-9) S η,β transforms a circular source to an ellipse with axis ratio q = e −η at positional angle β = θ/2. The shear can be represented as a 2-dimensional vector η ≡ (η + , η × ) ≡ (η cos θ, η sin θ). (2-10) Likewise a shear may be represented as a two-dimensional distortion (δ + , δ × ), etc. A shear (η + , 0) creates ellipses oriented to the x or y axes, while (0, η × ) aligns circular sources to axes at 45 • to the coordinate axes. The shear η is not a vector in the image space, but rather is a vector in a non-Euclidean shear space that we describe below.
The full set of shear matrices do not form a group under matrix multiplication as S η 2 S η 1 may be asymmetric (two-screen lenses can effect a rotation for this reason). But we can form a group with an addition operation for 2-dimensional shears defined as where R is the unique rotation matrix that allows S η 3 to be symmetric. The geometric meaning of the shears is preserved since R will leave a circular source unchanged. The simplest expression of the composition operation in terms of components is cosh η 3 = cosh η 2 cosh η 1 + sinh η 2 sinh η 1 cos(θ 2 − θ 1 ) (2-12) sinh η 3 sin(θ 3 − θ 2 ) = sinh η 1 sin(θ 1 − θ 2 ).
Note that the second equation is not symmetric in the two operands and hence the shear matrix group is non-Abelian. The identity element is η = 0 and the inverse of η = (η + , η × ) is −η = (−η + , −η × ). The addition formula in terms of distortion components is derivable from (2-12), and is given by Miralda-Escudé (1991): We omit the derivations of these equations, which are straightforwardly but tediously executed by composing the transformation matrices. A more elegant derivation follows from noting that the transformation Equation (2-9) transforms the complex plane as 1 (2-14) It will be useful to consider the limit where δ 2 ≪ 1: If we instead make δ 1 ≪ 1, the asymmetry of shear addition is manifested as a change to the azimuthal component formula:

The Shear Manifold
We define a metric distance between two points η 3 and η 1 in shear space as (2-17) The differential form of the metric can be derived by specializing Equations (2-12) to the case ds = η 2 ≪ 1, θ 1 = 0, θ 2 = θ, yielding which means that the metric is Note that the η version of the metric has the normal Euclidean form for the radial component, and the δ parameterization has the Euclidean form for the tangential component of the metric, but neither representation gives a fully Euclidean metric-the shear space is curved. The 2-dimensional shear manifold defined by this metric can be embedded in Euclidean 3-space as illustrated in Figure 1. This geometric depiction of shear is helpful in understanding the transformations of shears. Near η = 0 the surface is tangent to the Euclidean plane, so small shears add with Euclidean component-wise addition. The shear-space surface then curves upwards and as the conformal radius η grows large, the surface approaches a cylinder of unit radius about the z axis. If we project the shear surface onto the z = 0 plane, the radius vector in this plane is equal to the distortion δ. The δ vector is confined within the unit circle.

Definition of Shape
A shear is a transformation of the image plane; we next need a quantity to describe the shape of an arbitrary galaxy image. Let I represent some object whose isophotes are a family of similar ellipses. We can simply parameterize the shape of I by the shear η which produces this object from some object I 0 having circular isophotes, i.e. (2-23) We could thus call η the conformal shape of the object, and can think of a given ellipse as a location on the shear manifold. More commonly the distortion is used to define the shape; an object is said to have ellipticity e if a shear with distortion δ = −e makes it circular. We will use the symbol e since this quantity agrees with the traditional second-moment definition of ellipticities for truly elliptical objects. Equation (2-23) makes it obvious how an ellipse with shape η 1 will be transformed under the action of a shear η 2 : we simply add the shear to the shape using the addition rules of shear space [Equation (2-12)]: Likewise we can also say that a distortion δ maps the ellipticity e → δ ⊕ e. In general, an applied shear may be viewed as a shift of all shapes along the shear manifold. We will use e to represent the shape of an image, whereas δ represents a shear, which is a transformation of the image plane. The shape and shear spaces, however, transform identically under an applied shear.
A real galaxy has some image intensity distribution I(x) which may not have elliptical isophotes; we would like to define a shape for an arbitrary image. By analogy to Equation (2-23), all that is along this surface is the conformal shear η; the radius upon projection onto the xy plane is the distortion δ (or ellipticity e). At small η the manifold is tangent to the δ plane, and at large η approaches the unit cylinder. Two shear vectors s 1 and s 2 of length η = 1 are plotted from the origin, both on the true shear manifold and in the δ plane. The result of adding the two vectors is also plotted; displacements do not commute in this non-Euclidean space. needed is some definition of a "round" image. Let M(I) be any measurement applied to the image which has the simple property that if M(I) = 0, then M (S η I) = 0 for any non-zero shear η. Then the condition M(I) = 0 is our definition of roundness, and we can assign a unique shape η to an arbitrary object I(x) by the condition M(S −η I) = 0.
(2-25) Any shape defined by such a rule clearly transforms under an applied shear just as an ellipse does, namely via Equation (2-24). With any roundness measure M we therefore have a definition of shapes and their mapping under shears that follow the rules of addition in shear space. We do not attempt to prove that the solution to Equation (2-25) exists or is unique.

Shear in Fourier Space
For consideration of the effects of convolution upon sheared images it will be useful to ponder the action of shears in Fourier space. We first note that shearing an image I(x) by η is equivalent to shearing its Fourier transformĨ(k) by −η. For a diagonal shear: (2-27) = I(e η/2 k x , e −η/2 k y ). (2-28) A nondiagonal shear must also satisfy this relation since rotation of the real-space function corresponds to rotation in k-space. We can therefore just as easily define a galaxy's shape by a roundness criterion in k-space as in real space. This is useful when considering finite resolution ( §6) or when analyzing interferometric images with limited Fourier coverage.

Optimal Measurements (without Seeing)
In this section we derive an optimum method of measuring galaxy shapes in the case where the angular resolution and sampling of the instrument are assumed to be perfect. In §6 we will treat the more complex case of seeing-convolved images.

The Ideal Test for Roundness
We have defined the shape of an image in Equation (2-25) by asking what coordinate shear is needed to make the object appear round. We need to choose a measurement M which detects with the highest possible signal-to-noise ratio (S/N ) any small departure of the image from its round state M = 0. It can be shown that, under some sensible simplifying assumptions about M , the solution of Equation (2-25) becomes equivalent to finding the best least-squares fit of an elliptical-isophote model to the galaxy image.
We assume first that the measurement will be a linear function of I(x). Any non-linear method will prove extremely difficult to apply to the case where the image has been convolved with a PSFit will be hard to use the measurements of bright stars to correct the shapes of faint galaxies for convolution. The most general form of M is then where W is some weight function. The weight will be two dimensional, as the measurement must test for departure from roundness in both the η + and η × directions in shear space.
We consider first the weight component to detect a small change in η + . We can decompose the image I(r, θ) into multipole elements I m (r) via We are interested in the change in our measurement upon mapping of the image I(x) →Ĩ(x) = I(Sx) where S is a shear of amplitude η ≪ 1 oriented on the x-axis, as in Equation (2-3). The quantity for which we wish to optimize the S/N can thus be written as with w m an arbitrary radial function for each multipole. A little bit of algebra yields the transformation of multipoles where the primes denote derivatives with respect to r. For an object with truly circular isophotes, we have I m (r) = 0 for m = 0, and the only effect of the shear is to induce a quadrupole term I 2 = −ηrI ′ 0 /4. For objects without perfect circular symmetry, there are terms beyond the monopole. But for an object to be thought of as "round," the monopole term I 0 should dominate the higher multipoles. The monopole is also the only term guaranteed to be positive for all galaxies. Hence the largest effect of the shear will be to alter the m = ±2 quadrupole intensities (which are conjugates of each other as I is real). The optimal sensitivity to small shear should therefore weight only the quadrupole term: r dr I 2 (r)r 2 w(r) = 1 2π r dr dθI(r, θ)w(r)r 2 e −2iθ . (3-6) It is clear that this quadrupole test is the optimal linear measurement for objects with circular symmetry; for more general shapes, the shear has effects on other multipoles that can be measured and used to enhance signal-to-noise (related to the suggestions of Refregier (2001)). But this would require a knowledge of I 4 (r) and the other multipoles to construct the ideal formula; we settle on the simple quadrupole as the best general solution, as we are always guaranteed that I 0 is present and positive for any real source. The measurement of some weighted quadrupole is also the normal definition of ellipticity for weak lensing measurements (Miralda-Escudé 1991).
Combining Equations (3-4), (3-5), and (3-6) we obtain The noise in the measure of M can be derived in two limits: the most common case will be skydominated observations, for which the variance of the flux in area A is nA, where n is the number of sky photons per unit area (it is assumed that I is in units of photons). In this case we have, from Equation (3-6), the variance of each component of δM If we ignore the I 4 terms in Equation (3-7) as being dominated by I 0 terms, then the choice of weight function which optimizes the detectability of the shear η is With this optimal weight, the variance of the measurement M would lead to an error in each component of η equal to (3-10) If the object is much brighter than the night sky, then the noise is no longer uniform and the optimization becomes Var(M ) = 1 4π 2 r dr dθr 4 w 2 (r)I(r) cos 2 2θ (3-11)

Gaussian Objects
An elliptical Gaussian object, when sheared to be circular, will obey I(r, θ) = I 0 (r) = f 2πσ 2 e −r 2 /2σ 2 . (3-12) In the sky-limited case the optimal weight is the same Gaussian: Note that the optimal weight for shape measurement is in this case equal to the optimal filter for detection, i.e. a matched filter. If we define the detection significance ν as the signal-to-noise ratio for detection of the object with the matched filter, we find (3-15) (3-17) We therefore end up with the simple result that the error in each component of the shear is 2 over the detection significance.
The above derivations assumed that the center and the size σ of the Gaussian were known in advance. If there were a sky filled with Gaussian galaxies, we likely would not know in advance the size and location of each. We can determine the centroid in the usual manner by forcing the weighted first moments to vanish: dA wI re iθ = 0. (3-18) The weight for centroiding does not necessarily have to match that used for the shape measurement, but it is convenient to do so. The proper size σ w for the weight can be forced to match the size σ of the object by requiring the significance to be maximized: In the limit of a Gaussian with low background noise, the Equations (3-11) apply and the optimal weight is uniform. The detection significance in this case is just ν = √ f (with the flux f in photons), and we find again that the standard error in η is equal to 2/ν. In practice this situation can never be realized because the weight extends to infinity, and at large radii the sky or read noise will dominate the shot noise from the galaxy, and neighboring objects will impinge upon the integrations. We will henceforth confine our discussion to background-limited observations.

Exponential or Other Profiles
To obtain an optimized measurement of a real galaxy we would have to measure its radial profile and construct a custom optimized weight using Equation (3-9). The majority of galaxies are spirals or dwarfs which are typically described by exponential profiles: ( [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] This weight diverges at the origin, though all the necessary integrals of the weight are convergent. If the galaxy is truly cusped in the center, then the intensity near the center is very sensitive to small shears and is weighted heavily. In practice it is simpler to adopt a weight function that is universal (up to a scale factor), especially at low S/N where attempts to measure each individual profile would be pointless. There are a number of reasons to prefer a Gaussian weight: • The Gaussian drops very quickly at large radii, minimizing interference from neighboring objects. Integrals of all moments are convergent.
• Weights with central divergences or cusps are difficult to use in data with finite sampling, and also amplify the effects of seeing on the galaxy shapes. The Gaussian is flat at r = 0.
• Gaussian weights are analytically convenient, allowing many useful formulae to be rendered in closed form.
• Gaussian weights allow construction of the family of orthonormal basis functions that we will use in later chapters to compensate measured shapes for finite resolution.
• The Gaussian is not far from optimal for most galaxy shapes. For a well-resolved galaxy with an exponential profile, the Gaussian weight measures η with only 7% higher noise than the optimal weight in Equation (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20). In the presence of seeing, the difference between the Gaussian and optimal weight is even smaller. Recall that any weight we choose yields a valid definition of roundness and hence of shape; the Gaussian just incurs a small penalty in noise level.
The procedure for measuring galaxy shapes is therefore as follows: 1. Estimate a shape η for the image I and apply the shear S −η to obtainĨ.
2. Iterate the center and size of the Gaussian weight function until the centroid condition (3-18) and the maximum-significance condition (3-19) are satisfied.
3. Compute the second moments with the Gaussian weight function M (I) = dA w(r)I(r, θ)r 2 e −2iθ , w(r) = e −r 2 /σ 2 . (3-21) 4. If the real and imaginary parts of M are zero, then η is the shape of the object. If not, then we use the measured M to generate another guess for η and return to step 1.
The process is mathematically equivalent to measuring the second moments of I(x) with an elliptical Gaussian weight, and iterating the weight ellipticity, center, and size until they match the measured object shape. It is therefore an adaptive second-moment measurement. The method is also mathematically equivalent to finding the elliptical Gaussian that provides the best least-squares fit to the image.

Uncertainties in Shape Estimates
Once we have settled upon a weight of the form w = e −r 2 /2σ 2 , we can integrate δM from Equation (3-7) by parts and use (3-8) to calculate the variance in η. We will again ignore the I 4 terms; this means we may have a small tendency to over-or under-estimate our shape errors if galaxies tend to be boxy or disky. We first define the weighted flux f w , significance ν, and weighted radial moments Ir m w as The condition (3-19) for optimal significance is  and under this condition the variance in each component of the shape is The quantity a 4 is a form of kurtosis which is zero for a Gaussian image. The terms of order ν −4 arise from errors in the centroid determination, and are discussed further in §8.
The procedure for measuring an object of shape η requires applying a shear −η to the image coordinates x to produce a coordinate system x ′ in which the object appears round. The uncertainties in Equation (3-26) apply in this sheared coordinate system. Because the object is round in this frame, there is no preferred direction in the shear space and the uncertainty region is circular, with an uncertainty of ±σ η on each component (η ′ + , η ′ × ) measured in the sheared frame. We must reapply a shear +η to restore our measurement to the original coordinate system. This process is illustrated in Figure 2. The η coordinate transformations defined in Equations (2-12), in the limit of η 1 = σ η ≪ 1, indicate that the uncertainty region on η will be elliptical, with a shrunken principal axis in the circumferential direction of the shear manifold: If we instead use Equations (2-16), we can find the uncertainty ellipse in the ellipticity plane to be So on the e unit circle, the uncertainty ellipse shrinks radially by 1− e 2 and tangentially by √ 1 − e 2 as we transport the error region from the origin back to the original ellipticity e.
Note that our derivation assumes that the noise characteristics of the image are unchanged when we apply a shear. This is true in the background-limited case because the noise spectrum is flat, and our shear matrices S have unit determinant. For an image that has been smoothed or deconvolved, the power spectrum will have structure and the shear will alter the noise statistics. We will discuss this in §6 in the context of finite image resolution.

Comparison with Other Methods
Galaxy ellipticities for weak lensing were first determined by computing unweighted second moments of the intensity (Tyson, Valdes, & Wenk 1990). If the moment integrals are taken to infinity, then the measured ellipticities transform under shear using the addition rules in Equation (2-13), and furthermore the correction for PSF effects is extremely simple. It is clearly impractical, however, to carry the integrals to infinity, since neighboring objects will interfere, the noise is divergent, and, as noted by K00, many common PSFs have divergent second moments. So the initial methods generally used some sort of isophotal cutoff to the moments. This has the disadvantage of creating moments which are non-linear in the object flux. An alternative would be to use unweighted moments within a fixed circular aperture, but as noted by KSB, the noise properties of unweighted moments are far from optimal.
In the KSB method, the measured ellipticityê KSB is computed from the second moments measured with a circular Gaussian weight with size selected to maximize the detection significance ν.
[A different weight function was suggested by Bonnet & Mellier (1995).] The distinction between the KSB method and ours is that KSB always apply a weight which is circular in the original image plane; in our adaptive method, the weight is circular in the sheared image plane which makes the object round. Or, as viewed in the original image coordinates, the weight is an ellipse with shape iterated to match that of the object. This distinction has two consequences: first, our adaptive method yields lower uncertainties for non-circular objects because the weight is a better match to the image. This effect is minor, though, for objects with e 0.5, and a minority of images are more elliptical then this. triangle marks the true shape of a target galaxy, located in the e plane. The shape is determined by shearing the image until the galaxy appears round. In the ellipticity plane, we are moving the object along the vector to the origin. Panel (b) shows the location of the target (and the original coordinate grid) in the e plane after application of the shear that makes it round. The shaded region represents the uncertainty region for the shape in the sheared view-the error region must be circular because the image is round. Finally in (c) we undo the applied shear, shifting the target and the error region back to the original shape. This mapping, however, shrinks the error region by a factor 1 − e 2 ( √ 1 − e 2 ) in the radial (azimuthal) axes.
The second, more important advantage over the KSB method is that our definition of shape via Equation (2-25) guarantees that our measured ellipticitiesê are transformed by an applied shear via Equations (2-13). The circular-weightê KSB does not have this property-indeed for an object with elliptical isophotes,ê KSB does not equal the true ellipticity. In the KSB method, it is therefore necessary to calculate a "shear polarizability" for each object, describing the response ofê KSB to a small shear δ; this polarizability depends upon the radial profile and m = 4 moments of the object. The "polarizability" of our measured ellipticities is just the δ ≪ 1 limit of Equations (2-13): e → e + δ 1 + δ · e ≈ e + δ − e(δ · e). (3-31) This transformation rule, and the shapes of our uncertainty regions in ellipticity space, arise solely from the geometry of the shear manifold and are independent of the details of the galaxy images. This will simplify the following discussions of methods to derive a shear from an ensemble of measured galaxy shapes.
We implement the adaptive-weighted-moments scheme in the program ellipto, described further in Jarvis et al. (2002).

Combining Exposures
In a typical observing program, a given background galaxy is imaged in a number of different exposures, in one or more bandpasses. This is done to increase exposure time, permit rejection of cosmic rays, and/or gather color information. Multiple exposures can also reduce systematic effects by placing data for a given galaxy on different parts of the detector and in different seeing conditions.
We hence encounter the question of how to combine data on a given galaxy in different images to an optimal single measure of the shape. There are two possible approaches: 1. Measure the shape on each exposure, then create a weighted average of the measurements as the final shape.
2. Register and average the images, then measure the object on the combined image.
We first consider which offers the lowest noise on the final shape. Consider the task of combining N exposures, with the object having significance ν 0 on each exposure. Following Equation (3-26), the uncertainty in the shape of a nearly-round object measured from a single image will be The second term is the uncertainty due to centroiding error, and C and C ′ are constants of order unity. If we average measurements (Method 1), then we decrease σ η by √ N . If we average images (Method 2), then we increase ν by √ N . The net error on the shape in the two cases is then The two methods are equivalent, except for the centroiding noise. If ν 0 is not 5 than averaging images will produce better accuracy on η. Keep in mind that (a) the galaxy will not even be detected on the individual exposures unless ν 0 3, and (b) the galaxy is useless for weak lensing unless σ η 0.5, which requires N ν 2 0 16. When N 3, the centroiding penalty is small for any object that will be useful, so a combined image is extraneous. When N 5, there are many galaxies detectable on the summed image that are not detectable on the individual images, and a summed image has detectability and centroid-noise advantages.
There is a compromise, "Method 1.5," which has the practical advantages (delineated below) of Method 1, while retaining the small S/N edge of Method 2: that is to create a summed image and use it for object detection and centroid determination, so that Var(x 0 ) ≈ 1/N ν 2 0 . Then this centroid is used to measure shapes on individual exposures, and the shape errors are equivalent to Method 2. In practice we will combine deconvolved Laguerre coefficients ( §6.3) rather than measured shapes.
There are several reasons why it may be preferable to average catalogs instead of images: • Correction of shapes for PSF effects is paramount, and only possible if the PSF is constant or slowly varying across the image. If the different exposures in a summed image overlap only partially, then the PSF (and noise level) will jump discontinuously as one crosses the boundaries of component exposures. It is therefore preferable to correct for PSF effects on an exposure-by-exposure basis. If the PSF is very stable (e.g. a space telescope) or if the exposures all have nearly the same pointing (a single deep field), then a summed image will have well-behaved PSF variations.
• For the smallest objects, the exposures with the best seeing will contain nearly all the useful shape information and should be weighted heavily (Kochanski & Fisher 1994). Large objects are, on the other hand, measured equally well in every exposure. Averaging catalogs allows one to adjust the weights of different exposures on an object-by-object basis, whereas this is not possible when combining images.
• If there are exposures in different bands, then the optimal weighting of the exposures is dependent upon the color of the object. This is easily done when averaging cataloged shapes but not easily done by summing images.
• Creation of a summed image requires registration and interpolation of pixels. The latter process smoothes the noise field and causes subtle variations in the PSF, both of which complicate later analyses.
• An especially pernicious hazard to creating a summed image is that slight mis-registration of the component images will cause coherent elongations of the images, which if not corrected will mimic a lensing signal. This is discussed by KSB; in theory such effects are handled by a proper PSF correction scheme. This is a danger for Method 1.5 as well.
Some practical advantages to Method 2, combining images, are that: • The data storage and processing requirements can be lower for a single combined image if N is large.
• In Method 1, outliers (from cosmic rays) are rejected on an object-by-object basis, whereas in Method 2 rejection is pixel-by-pixel. If the galaxies are very oversampled and the cosmic-ray rate is high, Method 2 could salvage the un-contaminated parts of galaxy images that Method 1 discards.
For the simplest circumstances (a single-filter stack of images with common pointing), image averaging is easier and has few drawbacks. For multi-filter or mosaicked data, catalog averaging is needed. The hybrid Method 1.5 is best for such cases, though more work. In the rest of this section we detail procedures for each Method.

Combining Images
There are standard tools for combining exposures into a single image. We remark here upon a few special considerations when doing this for weak lensing observations. First, accurate registration is paramount. Our scheme for image registration is described in Jarvis et al. (2002).
Second, the use of median algorithms is commonplace but dangerous. Proper correction for PSF effects will require that the images of bright stars have precisely the same PSF as do the faint galaxies. But with a median algorithm, the bright, high-S/N stars will be constructed with a PSF which is a median of all the exposures. The images of faint objects, however, will tend toward a PSF that is the mean of all the exposures, because the noise fluctuations will dominate PSF variations. The final PSF will therefore vary with magnitude. A sigma-clipping average is much preferred over the median for the necessary task of cosmic-ray rejection when combining.
Similarly one must be careful about rejecting saturated pixels. There will be many stars which saturate only on the best-seeing exposures; if the saturated pixels are rejected, these stars will have final PSFs which are broader than the PSF for faint objects. One must take care to ignore stars which are saturated in any one of the exposures.

Combining Shape Measurements
Suppose that a given galaxy has been measured to have ellipticity e i in images i ∈ {1, 2, . . . , N }. We desire the e which best estimates the true ellipticity of the object. Using the results of §3.2, we see that in the absence of PSF distortions, the minimum-variance estimate of e will be that which minimizes the χ 2 given by Here, e ⊖ e i is equivalent to e ⊕ (−e i ), where ⊕ corresponds to the addition operator introduced in Equation (2-11), and S i is a covariance matrix which is σ 2 η I in simple cases.
Note that if the e i are measured in different filters, than the galaxy may have no single welldefined ellipticity. By "best estimate," then, we must mean that which offers the best sensitivity to a weak lensing distortion, and the minimum-variance combination of the e i is still the desired quantity.
The ⊖ is a non-linear operator, so we could use a non-linear minimization algorithm to find the value of e at which χ 2 is minimized. However, this is both impractical for time considerations and unnecessary since the values of σ η are usually small. Thus, we can linearize the subtraction operator T can be derived from Equations (2-13) The linearized χ 2 becomes which has a minimum at This is a standard least-squares solution for the mean of the e i given covariances Σ i in a Euclidean e space. In the simple case of S = σ 2 η I, the expression for Σ i simplifies considerably to (4-9) which is equivalent to Equation (3-30). However, when we apply corrections for PSF dilution, we will find that the covariance matrix is more generally an ellipse with axes aligned in the radial and tangential directions. That is, where R θ is a rotation matrix with θ = arctan e × e + . In either case the minimization (4-8) is numerically straightforward, and we are left with an uncertainty ellipse for the meanē. It is wise to implement some outlier-rejection algorithm in this process as well.

Estimating Shear from a Population of Shapes
Now we presume to have measured ellipticities e i for a set of N distinct galaxies, with known measurement uncertainties for each. Our final task is to create a statisticδ from the e i which best estimates the lensing distortion δ that has been applied to this ensemble. There are three main effects which must be considered in constructing the estimator: first, the e i respond differently to an applied distortion δ, as embodied by Equation (3-31) for true ellipticities, or by shear polarizabilities for the KSB estimators, so we need to know the responsivity R ≡ ∂δ/∂δ of our statistic. Second, the variety of ellipticities in the parent (unlensed) galaxy population causes shape noise in the shear estimate. In most weak lensing projects this is the dominant random error, and we wish to minimize its effects. Third, there is measurement error in each ellipticity, which we also wish to minimize in our shear estimator.
Most practitioners have adopted a simple arithmetic mean of e + and e × as estimators for the applied distortion (e.g. . Using the weak-distortion Equations (3-31), it is easy to see that, in the absence of measurement error, this estimator has a responsivity R = 1 − σ 2 SN and a variance Var(δ + ) = σ 2 SN /N , where we have defined the shape noise σ 2 SN ≡ e 2 + (the × component has the same properties and shape noise).
Others have realized, however, that rare, highly elliptical galaxies have too much influence on the arithmetic mean and should be deweighted. Cutoffs on |e| (Bonnet & Mellier 1995) or other weighting functions w(e) (Van Waerbeke et al. 2000) have been applied to the ellipticities and tested with simulations, but without any sort of analytic optimization or justification. Lombardi & Bertin (1998) consider the optimization of a general weighted sum of second moments (rather than ellipticities); this unfortunately couples the ellipticity measurement to the distribution of sizes of the galaxies and leads them to consider only weights which are power-law functions of the moments. Hoekstra et al. (2000) present a weighting scheme which incorporates both measurement error and the shape noise, and K00 gives a detailed discussion of optimal weighting for distortion measurements. Both are similar to our method in many respects, which we comment upon at the end of this section. must have circular symmetry. When each galaxy is sheared by δ + = 0.05, the galaxy distribution shifts to the right as in the middle panel. The right-hand panel shows the change in population under the applied distortion; this is the signal which we wish to detect. Shape noise arises from the Poisson fluctuations in the population, which is proportional to the left-hand panel's P (e). The optimal weight for δ + determination is the ratio of the right-hand to the left-hand panel.

Without Measurement Error
We start with an unlensed background galaxy population with ellipticities distributed within the unit circle according to A fundamental assumption of weak lensing is that the background is isotropic so that the unlensed population can have P depend only upon the amplitude of e ′ , not its orientation. The effect of a distortion δ is to map the background population to a new, anisotropic distribution P δ (e), as illustrated in Figure 3. We are given a sample of N galaxies from the new distribution, and our task is to estimate the δ which gave rise to the distribution from the original P .
One approach is to find the value of δ which maximizes the likelihood of the observed e i . This is true when A similar condition holds for δ × . These equations define the maximum-likelihood δ even for strong distortions-though there is not in general a closed-form solution for strong δ.
For weak lensing (δ ≪ 1), Equation (3-31) and the conservation of number can be used to -25derive P δ (e) to first order in δ.
It can further be shown that the maximum-likelihood estimator for δ takes the form The parenthesized expression is thus a weight function for combining ellipticities into a distortion. We can also show that this weight function is optimal in terms of S/N for weak distortions, as follows. Let us create an estimatorδ which is a general weighted sum of the ellipticities, d 2 e w(e)P (e) . (5-9) The response of this statistic to a small applied shear is where in the last line, we have dropped terms linear in e + or e × which average to zero over an isotropic population. With an isotropic population, the derivative ∂δ × /∂δ × = R as well, and the off-diagonal elements of ∂δ/∂δ are zero.
We may use Equation (5-11) to calculate the response of any weighted estimator by summing over the observed e i , because the small difference between observed and intrinsic distributions does not alter R to first order. In the case where we have some analytic form for P (e), we may replace the sums with integrals over the distribution to obtain R = d 2 e w(e)e + (∂P δ /∂δ + ) d 2 e w(e)P (e) (5-12) = d 2 e w(e)e 2 + P (e) 3 − 1−e 2 e d log P de d 2 e w(e)P (e) (5-13) In the absence of measurement noise, the variance inδ is due to shot noise. Assuming that the background galaxies obey Poisson statistics and their shapes are randomly assigned, we can propagate the Poisson errors through Equation (5-9) to get the expected error In (5-15) it is assumed that the sum is over a sufficiently large ensemble of background galaxies to sample the distribution P (e). Any weak lensing measurement has thousands of background galaxies, so this gives a direct estimate of the error in the shear.
The optimal weight is that which minimizes Var(δ + )/R 2 , which is ( 5-18) where the last line gives the optimized error inδ + /R, which is our calibrated estimate of the distortion. Equation (5-16) reproduces the maximum-likelihood solution in Equation (5-8). This may be compared to the distortion uncertainty for equal weighting w = 1, ( 5-19) We first see that a simple arithmetic mean of the ellipticities is the optimum estimator only if P ∝ (1 − e 2 ) α for some exponent α. For the real galaxy population, there can be a significant gain in accuracy through use of w opt over equal weighting. An extreme case is a population of randomly oriented circular disks, for which (5-21) With w = 1, we would have σ δ = 0.590/ √ N . The optimal weight diverges at e → 0 to take advantage of the extreme sensitivity of P (e) to distortion near e = 0. The integral in Equation (5-18) in fact diverges at e = 0, driving σ δ to zero-which would be a significant improvement over the equal-weighting case! Unfortunately any small measurement error or departures from circularity for the disks will smooth out the central spike in P (e), creating a finite value for σ δ . Figure 4 shows the P (e) measured for well-measured galaxies in the CTIO lensing survey (Jarvis et al., in preparation). These shape histograms are derived from 230,000 galaxies which are well resolved (R > 0.4, under the definitions in Appendix C) and have errors on the intrinsic ellipticity of σ η 0.03-primarily galaxies of magnitude 17 < m R < 22. The shape of P (e) is observed to be highly dependent upon the surface brightness (SB) of the galaxies. The low-SB galaxies show the rise at e > 0.8 expected of a disk population, but there the distribution drops at e > 0.95 instead of diverging-this reflects the finite thickness of the disks. There is also no pole at e = 0 for the low-SB galaxies, showing that the disks are not perfectly circular. The high-SB galaxies are presumably early types since there are very few with e > 0.5. While the value of P (e = 0) increases with surface brightness, it always remains finite, but with dP/de < 0. The ideal weight Equation (5-16) therefore grows as 1/e as e → 0, but the contribution to the S/N does not diverge at zero as for perfect disks. None of the P (e) curves is well fit by a single Gaussian or power law.
The intrinsic ellipticity variance σ SN varies from 0.30 to 0.49 between the highest and lowest SB bins. The optimally weighted distortion S/N per galaxy for high-SB, early types is 2-3 times higher than for the lowest-SB galaxies, indicating the desirability of incorporating some galaxytype discriminant-surface brightness, color, or concentration-into the weighting scheme. The only requirement is of course that the discriminant be independent of ellipticity. It seems likely that P (e) will vary substantially with magnitude.
The heavy histogram in Figure 4 combines all well-measured galaxies with 20 < m R < 21, which we henceforth use as a representative measure of the real galaxy population. The distribution has σ SN = 0.38, which would lead to σ δ = 0.44/ √ N for an unweighted average. The optimal weighting gives σ δ = 0.33/ √ N ; the weighting therefore gives a gain equivalent to a 1.8-fold increase in N . The gain in telescope time is at least as large. This gain is reduced, however, in the presence of measurement noise, which will tend to wash out the sharp feature in P (e), as discussed next.
We reiterate two favorable results of this section: first, the responsivity R and the variance of δ can be expressed exactly as direct sums over the observed population, for arbitrary choice of w; there is no need for calculation of polarizabilities nor recourse to simulated images. Second, we note that the variance inδ can be significantly below the canonical σ 2 SN /N if the ellipticity distribution P (e) has structure that is not washed out by measurement noise.

With Measurement Error
A galaxy image with true ellipticity e will be measured at some ellipticityẽ, with a probability distribution of p(ẽ|e). We consider a population of galaxies all having the same significance ν and resolution parameter R (see §6 and Appendix C), so that they all have a common p(ẽ|e). The The distribution of intrinsic ellipticities for modestly bright galaxies (m R ≈ 20) is plotted; we plot 2πeP (e) rather than P (e) as the nature of the population is more apparent. The distribution is highly dependent upon surface brightness µ R , presumably reflecting the difference between spheroid-and diskdominated galaxies. The dashed line is the distribution for an isotropic population of circular disks. The high-µ galaxies are more useful for distortion measurements. The heavy histogram combines all surface brightnesses in the magnitude range 20 < m R < 21. Though it is difficult to tell from this plot, P (e) is finite and increasing as e → 0. Optimal weighting takes advantage of this structure to reduce the noise in the distortion measurement. measured distribution of ellipticities observed under distortion δ will then bẽ P δ (ẽ) = d 2 e P δ (e)p(ẽ|e), (5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22) where P δ is the distribution of true ellipticities as in the previous section. The symmetry of P and p in the ellipticity plane guarantees that in the absence of distortion, the measured distributioñ P (ẽ) must again depend only upon the magnitude, not the direction, of the measured ellipticity.
Equation (5-22) is not strictly a convolution because the measurement function p(ẽ|e) may depend upon e and not simply upon (ẽ − e)-for example Equation (3-30) describes how the error ellipses contract as e departs from zero, even if the significance of the detection is held fixed. In §6 we will show that the behavior of the measurement error is different when the effects of PSF smearing upon the image are important, so at this point we will consider p(ẽ|e) to be, most generally, some kind of Gaussian whose 1σ ellipse depends only upon the magnitude e. An important point is that the functional form of p(ẽ|e) is unchanged by an applied distortion, since p is determined by ν and R, which are unchanged by a pure shear.
Another fact to keep in mind is that, with finite resolution and noise, it is possible to measurẽ e > 1, if the image noise makes the object appear smaller than the PSF in some dimension. Our formulae should therefore be tractable even forẽ > 1, and we cannot simply discard such measurements without contemplating the consequences.
We proceed as in the previous section, by assuming a distortion estimator of the form Given functional forms for the intrinsic distribution P (e) and the uncertainty function p, we could use (5-28) to derive an optimal weight, and use (5-25) and (5-26) to get the responsivity and noise for the estimator using this or any weight function. In most cases these integrals will not have analytic solutions.
The parenthesized quantity under the integral in Equation (5-13) is a galaxy's responsivity to shear, which depends upon the intrinsic shape. Equation (5-25) is the average of this responsivity for the galaxies with some measured shape. The measurement noise can cause these two quantities to differ; in other words, a naive determination of the responsivity is biased by the measurement noise. KSB-based methods will also suffer a calibration error due to this effect; binning the polarizabilities in parameter space can reduce the noise in the polarizability, but will not remove biases. Precision cosmology will require that such calibration issues be addressed-there are unfortunately no cosmic calibration standards for shear.
The need for P (e) in the above formulae is an unfortunate complication sinceP is the directly observed quantity. Note that the variance of the estimator can be expressed as a closed sum over the observed shapes (5-27), but the responsivity cannot. A precise calibration of the resultant shear/mass maps requires, therefore, that P be estimated either by deconvolving the observedP with the error distribution p, or by recourse to higher-quality images that give P directly.
Derivation of the the optimal weight also requires knowledge of the intrinsic P , but we can explore some generic cases, and make some approximations that give workable methods.

Approximate Form For Responsivity with Errors
We wish to have a form for R as a sum over the observed objects and applied weights, as in Equation (5-11), for the case of finite measurement errors. Toward that end we can take the derivative of Equation (5-23), which is greatly simplified if we assume that the measurement error, i.e.ẽ i − e i , does not have any first-order dependence on δ. While not strictly valid it is a good approximation. In this case where the brackets indicate an average of the true quantity at a given measured value, e.g.
Note that the weight function w may depend uponẽ directly, but also indirectly through some dependence in its covariance matrix Σ. If the measurement error function p(ẽ|e) and the intrinsic distribution P (e) have circular symmetry, then we must be able to write where k 0 and k 1 are functions only of the magnitude, not direction, ofẽ. We must also have e + e × ẽ = k 1 (ẽ)ẽ +ẽ× . Further manipulation, taking advantage of the isotropy of the parent popu-lation, yields This form for R depends only upon the observed quantities and the chosen weight scheme, except through the two functions k 0 and k 1 , which we will approximate below. The resemblance to Equation (5-11) is clear. With this equation and some integration by parts, we may also derive a form for the optimal weighting function: (5-34)

Special Case: Gaussians
In general the functions k 0 (ẽ) and k 1 (ẽ) must be calculated numerically using a presumed underlying P (e) for the background population, but analytic solutions are possible in the case of a Gaussian P (e) with variance σ 2 SN in each component (the shape noise) and a constant measurement error σ 2 e on each component. We find that both k 0 and k 1 are independent ofẽ: (5-35) The quantity f is the fraction of the total ellipticity variance that is attributable to shape noise. When the measurement noise is small, f ≈ 1, the ideal weight is close to (1 −ẽ 2 )/(σ 2 SN + σ 2 e ). This is quite similar to the weight adopted by Hoekstra et al. (2000).
For the Gaussian case, d logP /dẽ is also quite simple, so Equation (5-33) can be used. In our surveys to date (Smith et al. 2001;Jarvis et al. in prep.) we have adopted a weight that results from optimizing the Gaussian case.

Practical, Nearly Ideal Approximation
We obtain an approximation to the correct responsivity R and resultant ideal weight if we adopt the constant k 0 and k 1 functions in Equations (5-35) even for non-Gaussian P (e) distributions. The shape noise σ 2 SN may be defined as the assumed variance of the underlying e + and e × , and may be found by subtracting the measurement noise from the observed ẽ 2 . The measurement noise σ 2 e is known for each galaxy using the methods of this paper; since the covariance matrix forẽ is generally anisotropic, some representative scalar must be selected.
With these guesses for k 0 and k 1 in hand, R may be estimated with a sum over the observed galaxies using Equation (5-33) for any chosen weight function.
For the real-Universe shape distributions measured in the CTIO survey, we find that the following "easy" weight function offers very close to optimal distortion measurements: where σ η is the shape uncertainty that the object would have were it circular [cf. Equation (3-30)].
We can check the accuracy of our approximations numerically for chosen P (e) and p(ẽ|e) functions. We examine the case when P (e) is that shown in Figure 4 for galaxies with 20 < m R < 21, and the measurement error follows Equation (3-30). We find that the weight function given by Equation (5-34) is in fact very close to optimal for all noise levels, even when the simple approximations (5-35) are used for the k functions. The "easy" weight Equation (5-36) also performs nearly optimally, so most applications could use this weight and need not attempt to determine P (e).
A more critical question is whether the approximations (5-35) yield a proper estimate of the calibration factor R when used with Equation (5-33). Figure 5 shows how the simple R estimator compares to the correct value in Equation (5-25) for our choice of underlying distribution and the "easy" weights (5-36). The approximate form yields a responsivity correct to better than 5% for σ η 0.4. It is clear from the Figure that some detailed knowledge of the underlying P (e) distribution will be needed in order to calibrate lensing measurements to the one percent level.
The lower panel of Figure 5 shows the potential advantage of optimal weighting. When the measurement error is 0.2, there is little difference between various weighting schemes. For an unweighted distortion estimator, the accuracy levels out as the measurement noise drops below σ e ≈ 0.2. When optimal weights are used, however, the distortion errors continue to drop as the measurement error is pushed toward zero-the optimal weights take advantage of the e = 0 cusp in the shape distribution. Our "easy" weight scheme recovers nearly all of this potential gain.
To summarize, a practical method of weighting and calibrating the response in the presence of measurement noise is: 1. Determine the underlying σ 2 SN = e 2 + of the the intrinsic distribution. The measurement error σ η (or the full covariance matrix) of each galaxy is known from the formulae in previous sections.
3. Choose a weight function, for example the "easy" form Equation (5-36), or preferably the optimal form Equation (5-34) which can be derived from the observed distributionP (ẽ).
5. The estimator (and variance) must be scaled by the responsivity which is calculated with Equation (5-33). , as a function of the ellipticity measurement noise at e = 0. The weight function is the "easy" Equation (5-36) and the factors k 0 and k 1 adopt the simple heuristic approximation (5-35). The heavy (black) curve is for the population of 20 < m R < 21 galaxies, the upper (green) curve is for a low-SB sample (20 < µ R < 21), and the lower (red) curve for a high-SB sample (17 < µ R < 19), for which the intrinsic P (e) distributions are plotted in Figure 4. The simple formulation yields a calibration accurate to 5% or better in all cases, but 1% accuracy is difficult to achieve. The lower panel shows the uncertainty in the distortion determination when the galaxy shapes are combined with optimum weights (solid lines), the "easy" weights (5-36, dotted), and equal weighting (dashed). The line weights/colors code the galaxy sample, as above. Note that with optimal or easy weighting, the distortion errors continue to shrink even when measurement error is well below the canonical shape noise of 0.3.
If more is known about the intrinsic shape distribution, then more accurate functions for k 0 and k 1 may be derived and used in the sums.

Additional Weighting Considerations
The optimal weight may depend upon parameters other than the observed ellipticityẽ. It must, for example, depend upon the measurement errors as described above. If intrinsic shape distribution P (e) depends upon galaxy type, for example, then it may be advantageous to have different weight functions for each type-as long as galaxy type can be determined independent of e. Ellipticals may be distinguished from spirals, for example, by a concentration ratio such as the b 22 coefficient described below.
The expected shear will depend upon the redshift z of the background galaxy, and hence there may be a z-dependent weight determined from photometric redshift estimates of the background population. More crudely the apparent magnitude may be taken as an indicator of z and used in the weight formulation. The use of such additional weights depends upon the problem being addressed: see Smith et al. (2000) as an example.

Relation to Previous Methods
The optimal weighting scheme of K00 differs from ours in several respects. It assumes a different measure of e which does not follow our geometric transformation relations, so the mean polarizability of galaxies must be calculated in some parameter space, e.g. of flux, size, and e. Binning or smoothing in some such space is common to most of the KSB-based weighting schemes as well (Erben et al. 2001;Hoekstra et al. 2000). The variance of the estimatorsδ i are calculated within each bin, then weights are chosen inversely to the variance of each bin to create a minimumvariance weighted mean estimator. There are three contributors to the variance from each bin: 1. The intrinsic ellipticities of galaxies within the bin are drawn at random from the parent distribution. If |e| is one parameter of the space, then only the direction θ is allowed to vary within the bin. This is of course the shape noise.
2. The measured shape of a given galaxy is drawn from the measurement-error distribution.
3. The number of galaxies within the bin is drawn from a Poisson distribution. If |e| is a parameter, then this Poisson noise includes some elements of the shape noise (1) and measurement noise (2).
All three of these effects are important to optimization; all are included in our formulation and (implicitly) in that of K00, so we expect them to be essentially equivalent in the long run-this is not the case, though, for some of the heuristic or parameter-space weight formulations in the literature. The virtue of our scheme is that the nature of the weight function is apparent given the intrinsic shape distribution P (e) and the measurement errors, and there is no need to choose a parameter space for weight selection and polarizability smoothing. Our formulation tells us when further parameters might be desirable, namely when P (e) changes significantly.

Measurements with Finite Resolution
The preceding sections outline a method for optimal recovery of weak distortion from galaxy images, and rigorous estimation of the uncertainties on these shears, for the case when the detector views the galaxies with perfect resolution. Unfortunately, the finite resolution of real observations has a strong effect upon shape measurements in every weak lensing observation to date, even those using the Hubble Space Telescope. Finite resolution produces two deleterious effects: 1. A PSF which is not circularly symmetric can induce ellipticities on the images, thus breaking the intrinsic isotropy of the background galaxy population and mimicking a lensing distortion. This is a bias induced by asymmetric PSFs. Since present-day weak lensing surveys are seeking distortion signals well below 1%, measured shapes must be compensated for even the slightest asymmetry in the PSF with some smear correction.
2. Convolution by a circularly symmetric PSF will make most galaxies appear rounder, driving e → 0. This is therefore a dilution of the true lensing signal. While this mechanism cannot create a lensing signal where there is none, it can misleadingly modulate the lensing signal or cause a calibration error in the inferred mass distribution.
Most (but not all) approaches to PSF corrections treat the bias and dilution effects in separate steps. To our knowledge, all published weak-lensing observations have incomplete PSF correction, leaving systematic distortion errors of a fraction of a percent or higher. While in most cases these residual systematic effects have not altered the validity of the authors' conclusions, the proper correction of PSF effects is presently the largest barrier to the use of weak lensing for precision cosmology.
In §6.1 we review some existing approaches to these problems, in §6.2 we contemplate how one would ideally wish to approach the problem, and in further sections we develop two means of implementing a nearly-ideal approach: one which treats bias and dilution separately, and another which corrects both problems simultaneously with a limited form of deconvolution.

The Unweighted Ideal
Unweighted second moments of galaxies are ideal measures of ellipticity-not only do the derived ellipticities transform according to the rules of shear space, but correction for PSF effects is in principle quite simple because the unweighted second moment of the image is just the sum of the original moment and the PSF moment. Thus by subtracting the PSF moments from the measured moments we simultaneously correct for bias and dilution, and obtain the image shape. In the case where the PSF is round, the dilution of a true (pre-seeing) image-plane ellipticity e i to an observed (post-seeing) ellipticity e o is described by the exact equation The resolution parameter R is determined by the unweighted second radial moment of the measured image r 2 o relative to that of the PSF r 2 ⋆ . Two things to note: first, the error ellipse on the dilution-corrected, pre-seeing ellipticity e i is magnified by 1/R from the original measurement error Equation (3-30), and is further stretched in the radial direction by the uncertainty in R itself. Thus the error ellipse is no longer simply described by a single σ η . Second, Equation (6-1) can give rise to |e i | > 1, if the noise makes the galaxy look smaller than the PSF about some axis. We cannot arbitrarily discard such measurements without creating a bias in our mean shear. These two phenomena are common to all modes of PSF-dilution correction.
This blissfully simple dilution correction is spoiled by two major problems: first, as discussed above, unweighted second moments have divergent noise properties and for this and other reasons are not practical shape estimators. An equally serious problem noted by K00 is that the second moments themselves are divergent for many realistic PSFs. Further, many galaxies follow deVaucouleurs profiles, for which the second moment converges very slowly.
The simple formulae (6-1) are still valid under the special circumstances that the object and PSF are both Gaussians. The post-PSF object is again a Gaussian, and deconvolution of Gaussians is a simple subtraction of second moments. Hence any shape-measuring algorithm which extracts the proper ellipticity for a Gaussian ellipsoid would allow PSF dilution correction via Equations (6-1) in this limited (and unrealistic) case.
Some early weak lensing measurements (Valdes, Tyson, & Jarvis 1983) adopt second-moment subtraction as a means of PSF correction, despite the fact that this method is not exact when isophotally-bounded or weighted moments are used, and the images are not Gaussian. This would not suffice, however, for the more sensitive measurements being done today.

Heuristic Methods
In the case of unweighted moments or Gaussians, Equation (6-1) would indicate that a regression of the lensing signal against r 2 −1 o would yield a distortion free of PSF effects as r 2 −1 o → 0. Mould et al. (1994) have attempted to measure very weak shears in the presence of PSF effects using such a regression (though against r −1 , in which case a linear relation is not expected). Even with weighted moments, we expect the PSF dilution and bias to decrease as the object becomes well-resolved, so there is some basis to this method, even if it is not exact. Other problems, however, are that the distortion is not likely to be the same for all sizes of galaxy as they likely lie at different distance. Also the regression will lead to substantially higher noise than a more direct dilution correction.
Another approach to the dilution correction is exemplified by FT97, who attempt no analytic correction, instead calibrating the dilution effect by measuring simulated background galaxies which have been subjected to the same distortion, seeing, sampling and noise as the real images. Such a simulation is an essential test of any weak-lensing methodology. The difficulty with sole reliance upon simulated data is that the result is extremely sensitive to one's ability to exactly match the size-magnitude distribution of the true galaxy population, because the dilution correction is a strong function of size (as in Equation [6-1]) in the typical regime of slightly resolved galaxies. Further, as we show below, the dilution correction depends upon detailed higher-order moments of the galaxy images, which would be very difficult to simulate faithfully. One alternative is to use high-resolution, high-S/N images from HST instead of simulated galaxies-but the total sky area imaged to sufficient S/N by HST is a tiny fraction of a square degree, too small for rigorous calibration tests. It would be preferable to have an analytic correction for dilution, and use the simulated data to spot-check the accuracy of the analytic method.

Perturbative Methods
A step beyond the unweighted-mean approximation to the bias correction is taken by KSB and by FT97. Both make the assumption that the anisotropy of the PSF can be described as a small anisotropic convolution applied to a larger, circularly symmetric PSF. In this case, the effect of the tiny asymmetric deconvolution upon the weighted second moments of a given image can be expressed as a fourth-order weighted moment of the image, which KSB christen the smear polarizability. Given the smear polarizability of an image and a measure of the anisotropy of the PSF, the measured second moments are corrected analytically for the PSF bias.
The FT97 method differs in that the correction for PSF anisotropy is applied to the image rather than to the measured moments: a minimal 3 × 3 convolution kernel is created which will "circularize" the PSF. The galaxy shapes are measured after this kernel is applied to the image.
The primary drawback to these methods (K00) is that the approximation upon which they are based often fails: a typical diffraction-limited PSF in no way resembles a small convolution to a round PSF, and even a simple aberration such as coma creates PSFs that violate this condition.
These methods have features, however, that we wish to retain in any improved formulation. They are easily adapted to a PSF that varies across the image-the requisite moments of the PSF are measured wherever a star falls upon the image, and interpolated to the location of each galaxy. The FT97 method, by fixing the image, frees us from having to measure the higher-order moments that comprise the smear polarizability, though at the expense of a slight reduction in image resolution and/or an increase in the image noise.
The perturbative methods correct only the PSF bias, not the dilution, because the dilution cannot be considered a small perturbation in any extant data set. FT97 calibrate the dilution with simulations, as mentioned above. The original KSB work made use of empirical calibration as well, as the shear polarizability (cf. §3.3) measures only the susceptibility of the image to a distortion which might be applied after the PSF convolution. Wilson et al. (1996) suggest an empirical calibration by deconvolving the real images, applying a shear, reconvolving, and remeasuring to determine the response. Luppino & Kaiser (1997) introduce the pre-seeing shear polarizability, which approximates the susceptibility of the KSB weighted moments to a shear applied before the PSF. The pre-seeing shear polarizability is, to its lowest order, similar to the resolution factor R introduced above for unweighted moments. There are, however, fourth-order moments of galaxy and stellar shapes involved, and all are measured with Gaussian weights so the noise does not diverge.
The KSB method updated to use the pre-seeing shear polarizability is exact for Gaussians and in the limit of a compact anisotropy kernel, but is not exact in the general case (K00). There are additional ambiguities regarding the appropriate size of the Gaussian weight to be applied when measuring the PSF moments (Hoekstra et al. 1998). The updated KSB method is in wide use, and several papers have investigated how accurately it performs for simulated galaxies and PSFs that are not Gaussian (Erben et al. 2001;Bacon et al. 2000). While it is clear that in many circumstances KSB is good enough, we would prefer to understand and overcome its limitations.
We will demonstrate below that the KSB and FT97 bias corrections are the first terms in a series expansion of the deconvolution of the galaxy image. Rhodes, Refregier, & Groth (2000) investigate the KSB method in some detail, carrying forth the transformations equations to higher order than do KSB. The PSF corrections, however, still require substantial approximations. Below we construct a hierarchy of all the weighted moments of the image and give general formulae for their transformations under shear, convolution, and other operations. Kuijken (1999) suggests that one determine shear by summing the images of all the galaxies in some cell, then comparing this high-S/N summed galaxy image to the PSF. The summed galaxy image should approach circular symmetry in the source plane; the shape of the summed image can thus be modeled relatively simply as a circular source with arbitrary radial profile, sheared by the lens to have elliptical isophotes, then convolved with the PSF. He adopts a flexible parametric representation of this mean radial profile. A candidate profile with candidate lens distortion can be convolved with the known PSF and compared to the measured mean image. The radial profile and distortion are then adjusted to give a best fit. This method can be viewed as a limited form of deconvolution: the only characteristic of the deconvolved image one cares about is the ellipticity; the multipoles beyond quadrupole are irrelevant to the measurement and are discarded.

Deconvolution Methods
In Kuijken's method, the higher-order multipoles are discarded by averaging over many galaxies. This may be difficult in practical situations, where the PSF and/or the distortion signal are varying across the field too rapidly to gather sufficient background galaxies to sum. Another drawback is that summing galaxies' images may not be the optimal way to combine their information on the shear. Kuijken suggests applying the method to individual galaxies given these potential problems, which amounts to an assumption that the galaxy ellipticity is constant with radius. The accuracy of this method is not discussed in detail. But the method does have the strong advantage of being able to cope with arbitrary PSF behavior, and simultaneously removing both the bias and dilution effects from the measurement.

Commutator Method
K00 introduces a new approach to PSF correction, deriving from the PSF an operator which can be applied to the observed image to effect the transformation of applying a given pre-seeing shear. This operator is derived by considering the commutation of the shear and convolution operators in Fourier space. With the operator in hand, the response to a pre-seeing shear can be determined directly from the post-seeing image, given sufficiently detailed knowledge of the PSF.
The K00 formulation is the first, to our knowledge, to offer an exact correction for PSF bias and dilution. We take a different approach below, and then offer some comparison of the two approaches.

Optimal Methods in the Presence of a Convolution
In §2.5 we noted that a real-space shear is equivalent to an opposing shear of the Fourier-space image, so that we can conduct the roundness test in Fourier space. One virtue of using a Gaussian weight function for our roundness test, as in Equation (3-21), is that this test takes the exact same form in k-space, namely The only difference is that the Gaussian weight has width 1/σ in k-space, i.e. a broader weight in real space is narrower in Fourier space. Furthermore, our assumed uniform white noise in real space transforms to uniform white noise in Fourier space. Therefore the entire derivation of the optimal weight in §3.1 could have been executed in k-space without any change in the result.
The effect of the PSF convolution is to suppress the image by some transfer functionT (k). With perfect knowledge of the PSF, we can deconvolve the observed imageĨ o to retrieve the image-plane transformĨ i =Ĩ o /T (here we mean the image plane of the gravitational lens, before the application of seeing). A deconvolution would remove the PSF bias entirely, but the noise is no longer homogeneous in k-space, having been amplified (perhaps infinitely) by 1/T . If the PSF is close to circularly symmetric, thenT is nearly independent of direction. If we now imagine making our roundness test on the deconvolved k-space image, we adapt Equations (3-7)-(3-9) to give The optimal filter is therefore narrower in k-space than both dĨ/d(k 2 ) andT (k). Hence in real space, on the deconvolved image, the optimal filter is broader than both the object and the PSF. This means that we should, sensibly, restrict our roundness test to the region of k-space that (a) has signal in the true galaxy image, and (b) is not suppressed below the noise by the convolution.
Our strategy might then be to create some kind of deconvolved image, and apply a Gaussian weight to test for roundness in the deconvolved k-space (which is equivalent to using a Gaussian weight in real space). For a Gaussian object with pre-seeing size σ i and Gaussian seeing with size σ * , the optimal weight in deconvolved k-space is a Gaussian with size (2σ 2 * + σ 2 i ) −1/2 . Such a roundness test is equivalent to a roundness test on the observed, real-space image with a Gaussian weight of size σ 2 i + σ 2 * = σ 2 o . So the optimal size of the weight is again matched to the size of the observed image.
Recall that our algorithm for measuring δ requires that we shear the coordinates until the object appears round (M → 0). We want to apply this shear to the deconvolved image. If the object is not round to begin with, then in this sheared coordinate system the transfer functionT (k) will no longer have azimuthal symmetry, which will invalidate the above derivation of the optimal weight. We will still have a valid measurement of the shape of the deconvolved object, but possibly with sub-optimal noise level. The increase in noise is a second-order effect, however, so we will not bother to re-optimize the roundness test for this asymmetric noise spectrum.
When the transfer function is anisotropic, then the noise spectrum in the deconvolved image is also anisotropic. There are subtle second-order effects, described in §8, which bias the orientation of the measured shear in the presence of an anisotropy in the noise. Such a noise anisotropy is present after any deconvolution of an anisotropic PSF. A noise anisotropy is also present in the observed image if the PSF is anisotropic and we are not strictly sky-limited. Furthermore, K00 points out a selection bias that can creep into the shear measurements even with a perfectly unbiased shape measurement algorithm. We will discuss means to defeat these biases in §8.
Applying a Gaussian weight in deconvolved k-space leads, formally, to divergent noise ifT (k) = 0 for some k. This is a real problem, as any finite-sized telescope must produce a transfer function that is identically zero beyond some critical k c . The deconvolved image hence has infinite noise for k > k c , while the Gaussian weight remains finite. As we shear the k-space to make our source appear round, the "wall" of infinite noise moves inwards to e −η/2 k c . In order for our method to remain feasible, our deconvolution algorithm must not attempt to fully deconvolve those portions of k-space at or near k c , so that the noise remains bounded. There is hence a balance to be struck in executing the deconvolution: we want the carry out the deconvolution to sufficiently high order that the effects of the PSF upon the source ellipticity are removed, but we do not want to deconvolve high-order details that will increase the noise. It would seem, intuitively, that this is possible, since the ellipticity we seek, the Gaussian-weighted quadrupole moment, is a low-order characteristic of the deconvolved image.
The method we describe below is based upon an expansion of the image and PSF into hierarchies of Gaussian-weighted moments-essentially an eigenfunction decomposition. Convolution corresponds to a matrix operation on the vector of moments. The moment vector, and hence the convolution matrix, are formally infinite, but we can choose to truncate the description at some order which we believe to contain all the useful information on the image ellipticity. The convolution matrix is then finite and can be inverted, and the deconvolution executed as a matrix operation. The high-order moments are not deconvolved, so the noise in the deconvolved moments remains finite and in fact close to optimally small. This moment-based method reduces the deconvolution (as well as other transformations) to a matrix multiplication, which can be executed on an exposure-by-exposure basis even for objects with very low S/N on a single exposure, and is therefore very practical for the purpose of weak lensing measurements.

Definitions
The simplicity of the formulae for deconvolution in the special case of Gaussian objects, plus the utility of the Gaussian weight for shape measurements, lead us to seek a description of the image in some Gaussian-based expansion. To maintain the simplest form for convolutions, we look for a decomposition of our images into components which are eigenfunction of the Fourier operator. Such functions will also be eigenfunctions of (−∇ 2 + r 2 ), hence we are led to the eigenfunctions of the 2d quantum harmonic oscillator (QHO). 2 In one dimension, the QHO eigenfunctions are each a Gaussian times a Hermite polynomial. The Edgeworth expansion, familiar to many astronomers, is a decomposition into 1d QHO eigenfunctions. The 2d version we call the "Laguerre expansion," using the QHO eigenfunctions q (x) are the Laguerre polynomials, which are defined by the generating function (Abramowitz & Stegun 1965) ( The Laguerre polynomials satisfy many recurrence relations; the following provide a way to calculate them rapidly: A QHO with wavefunction ψ pq has angular momentum m = p − q. We will also make use of the quantum number N = p + q, which is the excitation energy of the state. Any two of {N, m, p, q} suffice to specify the state. The intensity multipole functions I m (r) are composed from the ψ N m for N = |m|, |m| + 2, |m| + 4, . . .. The polynomial in ψ pq has terms up to order N in r, and ψ pq can also be expressed as the Gaussian times a (complex) polynomial of order N in x and y. A few low-order ψ pq are plotted in Figure 6. Refregier (2001) has independently introduced the application of QHO eigenfunctions to galaxy shape analysis. Some of the results and ideas presented in this section are presented therein, and applied in Refregier & Bacon (2001), with different notation. Refregier (2001) also presents useful relations for a Cartesian-based family of 2d QHO eigenfunctions. We make use only of the polar family, which are eigenstates of the angular momentum and hence have strong rotational symmetry.
The Laguerre functions have many properties we will find useful. As eigenfunctions of the harmonic oscillator, they are orthonormal, up to the factor σ which we introduce in order to give where N = p + q is the number of quanta in the state. The ψ 20 component is most important as it responds in first order to shear, and hence its absence is our test for "roundness." the b pq units of flux: Thus b pq is a Gaussian-weighted moment of the intensity. Since I is real, we must have b pq = b qp . When measuring the b pq from an image with uniform white noise, Equation (6-16) yields a covariance matrix for the b pq that is diagonal: (6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17) where n is the number of counts per unit area. The variance in b pq is shared between the real and imaginary parts, except for b pp , which must be real. Non-uniform noise, e.g. shot noise from the galaxy itself, produces a more complicated covariance matrix, as described in Refregier (2001). The significance ν of detection with the Gaussian filter is given by Our algorithm for measuring object shapes requires finding the coordinate system in which the centroid is zero (Equation [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]) and the roundness test (Equation [3-6]) yields zero. When disregarding seeing, we also set the weight size σ by maximizing the significance (Equation [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]). These conditions are succinctly stated by the Laguerre coefficients: The first two of these equations involve both real and complex components, the third is real. We must satisfy these equations by translation, shear, and dilation of the object (or of the underlying coordinate system). These operations can be expressed as transformation matrices acting upon the vector b = {b pq }. The determination of shape is thus efficiently executed by measuring b in the original coordinate frame, converting Equation (6-16) to a sum over pixels. Henceforth we can iterate to a solution of our three conditions by manipulating b, and there is no need to return to the pixel data.
In the case of finite resolution, we wish to satisfy Equations (6-19) and (6-20) for the deconvolved image. We can express the deconvolution as a matrix operation on b as well. So we need to find the matrix equivalents of the translation, dilation, shear, and convolution transformations.

Raising and Lowering Operators, Point Transformations
The raising and lowering operators for the 2d harmonic oscillator wavefunctions have the properties: Note that the commutators of the operators are all zero except for [a p , a † p ] = [a q , a † q ] = 1. These operators can be used to formulate the transformation matrices we need. Consider first a dilation of the image I by factor 1 + µ with µ ≪ 1 to a new image I ′ : An infinitesimal dilation µ thus mixes coefficients (N ± 2, m) into the (N, m) coefficient-as expected, multipole order m is conserved by the transformation. The dilation transformations are thus generated by the matrix d = 1 − a σ p a σ q + a σ † p a σ † q , (6-28) and the transformation matrix D for a finite dilation e µ can be expressed as The finite dilation will preserve m as the generator does, but will mix terms of N + 2j together, with j any integer. The matrix elements for the finite dilation are derivable in closed form and in a rapid recursion form; the latter are derived in Appendix A.
The generator s for the shear transformation is also easily expressed in terms of the raising and lowering operators. For a shear η ≪ 1 along the x axis, = (1 + ηs)I, (6-32)  In (6-34) we have inserted the phase factors which result from a shear at arbitrary position angle β. To leading order, the shear mixes in states with p or q changed by ±2; in the finite case, the shear transformation mixes together b pq for which p and q each change by any even number. Note that this equation generalizes the weighted-moment transformations in Rhodes, Refregier, & Groth (2000) to all orders. Note also that the KSB method identifies the quantity b 20 /b 00 as the ellipticity, so the "shear polarizability" is compactly expressed using Equation (6-34) with p = 2, q = 0.
An infinitesimal translation of the image by (x 0 , y 0 ) gives To leading order the translation changes p or q by ±1, and the finite translation T z mixes all the states together. Appendix A gives the matrix elements for all the finite transformations.
The "smear polarizability" coefficients of KSB could be quickly derived at this point by noting that an infinitesimal convolution along the x direction can be expressed as an average of two images translated by ±x 0 , which will lead to a smear generator that is second order in the raising and lowering operators. Appendix C gives a closely related derivation.

Shape Measurement from Laguerre Decomposition
Our algorithm for shape measurements (with perfect resolution) is to find the translation z, dilation µ, and shear η which must be applied to the image to satisfy the conditions (6-19)-(6-21). If b is the vector of Laguerre coefficients for the image I, we can write  and the elements of b ′ with (pq) = (10), (11), and (20) must vanish. This is in general a complex non-linear equation, but for µ, η, z ≪ 1 the linearized Equations (6-27), (6-34), and (6-41) yield  writing η = η + +iη × . These linear equations may be solved explicitly for the desired transformation coefficients η, µ, and z. The solution appears much simpler if the object is nearly round, such that b pq ≪ b 00 for p = q. In this case the transformation parameters are, to leading order, From these equations it is clear that b 20 , b 10 , and b 11 are the primary carriers of information on shape, centroid, and size, respectively. Since the shape component η + is (η +η)/2, its uncertainty is to leading order Here we have made use of Equation With a more tedious procedure we may derive the solution for η to second order in b pq /b 00 . We take b pp ≪ b 00 for p > 0 this time in order to simplify the expression: The terms in the first line arise from solving for η = 0; the terms in the second and third lines arise when simultaneously solving the centroid and size constraints, respectively. From these terms one can derive the O(ν −4 ) terms in the shape uncertainty if desired. But more importantly, this expression uncovers a very significant source of bias in shape measurements: the presence of the b pq in second order means that it is possible for noise to be rectified in the determination of e. If the noise is anisotropic-as is the case after correction for an anisotropic PSF-then e will be biased. We will explore this in more detail in §8.

Fourier Transforms and Convolution
The observed galaxy intensity I o (x) is the convolution of the image-plane intensity I i (x) and the stellar PSF I ⋆ (x): (6-54) The post-seeing, pre-seeing, and stellar images can be expressed as vectors b o , b i , and b ⋆ , respectively, of coefficients over our eigenfunction sets ψ σo pq , ψ σ i pq , and ψ σ⋆ pq . The convolution can then be expressed as a matrix relation (6-56) We will effect the deconvolution by inverting a truncated version of the matrix C. Its coefficients are determined by the relation The convolution is more easily expressed in k-space. With the usual definitions (denoted as "System 3" by Bracewell (1978))Ĩ 58) the convolution becomes a multiplication, and the matrix coefficients in Equation (6-57) can also be expressed as We now make use of another remarkable property of the Laguerre-exponential eigenfunctions, which is that they are their own Fourier transforms. First we note that ψ σ 00 is a two-dimensional Gaussian, and the transform (6-58) is easily executed to yield a new Gaussian, The functional form ofψ 00 (k) matches that of ψ 00 (x), except that we must send σ → 1/σ. Next we can use the definition of the raising operators (Equation [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22]) to note that Thus the k-space raising operator has the same form as the x-space raising operator, save the σ → 1/σ transformation and an additional factor of −i. The same is true of a † q . Since ψ 00 and the raising operators are each unchanged by the Fourier transform, it must be true that all the ψ pq are their own Fourier transforms as well, with appropriate scaling factors of i and σ. More precisely, we have (cf. [6][7][8]) So the problem of convolving two of our eigenfunctions is reduced to the simpler problem of multiplying the same two eigenfunctions. We may reach several conclusions immediately: • The convolution of ψ σ i N i m i with ψ σ⋆ N⋆m⋆ will produce an observed image with azimuthal order m o = m i + m ⋆ . Furthermore the multipole phase of the observed image must be the sum of the phases for the original and PSF images' coefficients.
• If we choose the scale size σ o for the observed image eigenfunctions such that  then the convolution contains components ψ σo Nomo only for N o ≤ N i + N ⋆ . Recall also that we must have N o ≥ m o = m i + m ⋆ , so it must be true that ψ σ i p i q • ψ σ⋆ p⋆q ∝ ψ σo (p i +p⋆)q for q = 0.
In Appendix B we give a recursion relation to calculate any of the elements C p i q i p⋆q⋆ poqo that we need to calculate the convolution, and thus the deconvolution.

Deconvolution, Noise Amplification, and Truncation
With the formulae for the convolution matrix C in hand, we can investigate the nature of the trade-off between the fidelity of the deconvolution-i.e. the extent to which effects of the PSF upon the shape are removed-and the noise level of the deconvolution. We must truncate the b vectors at some finite order N to implement the deconvolution; higher N will remove more of the PSF effects, but also increase the noise level. We can illustrate this phenomenon by considering the simplest case of convolution by a unit Gaussian PSF, b ⋆ pq = 2 √ πδ p0 δ q0 . Since the PSF has only m ⋆ = 0 terms, we will have m o = m i , so the convolution matrix is block-diagonal and we may deconvolve each m independently. If we choose σ 2 o = σ 2 i + σ 2 ⋆ , then the convolution matrix coefficients are zero for N o > N i , and hence each block matrix is also upper-triangular. This makes the inversion easy.
Using the results of Appendix B, defining a deconvolution parameter D = σ 2 i /σ 2 o = 1 − σ 2 ⋆ /σ 2 o , and using Equation (6-17) for the covariance of the measured moments, the deconvolved value of b i 00 and its variance are In practice we could not measure the b o pp to infinite p and we must truncate the sum (6-64) at some finite p. The more terms we include, the more accurately we describe the deconvolved b i 00 , but each additional term in the deconvolution adds more noise. For D > 1/2, the added noise in successive terms drops as p increases, and we could in principle do a complete deconvolution of b i 00 with finite variance. For D ≤ 1/2, however, the noise in the deconvolved moment diverges for p → ∞, and we are forced to truncate the deconvolution matrix.
For ellipticity measurements we are primarily interested in the b i 20 moment. For this case of a Gaussian PSF, the general b i m0 moment deconvolution is (6-68) We note that these moments also demonstrate the properties that they are noisier than the observed moments, and infinitely so if D < 1/2 and the deconvolution is not truncated. The noise increases as m increases, as we would expect, since higher spatial frequencies must be recovered.
For a purely Gaussian PSF with D > 1/2, it is possible to complete the deconvolution with finite noise. A real PSF must have components beyond b ⋆ 00 , since the transfer function vanishes above a critical k value. When we invert the convolution matrix for such a PSF we will find that the variance series akin to (6-65) would diverge for any D < 1. The best truncation value for the matrix will depend upon the form of the PSF and the accuracy to which we demand correction of PSF effects. We will examine some specific cases in Jarvis et al. (2002).

Analytic or Kernel Correction for PSF Bias?
There are two general means of eliminating the shape bias induced by the PSF. One alternative is to measure the anisotropy of the PSF carefully and apply analytic corrections to the measured objects, as occurs naturally within the deconvolution framework described in the preceding paragraphs. The KSB formalism contains an approximate analytic correction. The other method is to convolve the image with a spatially-varying kernel which removes the anisotropy from the PSF, as first demonstrated by  and further advocated in K00. Removal of the PSF bias is the most critical task for a weak lensing method: the PSF dilution is only a calibration error or signal modulation, but the PSF bias introduces a first-order false signal to the lensing analysis.
Ideally, analytic correction is preferred. The convolution kernel will always degrade the image resolution to some extent. Many PSFs cannot be "rounded out" unless the kernel is of size comparable to the PSF itself (K00). Furthermore the convolution will perturb the noise spectrum of the image, complicating error estimation. The kernel method, however, can for simple kernels be faster than an involved analytic correction, especially if eliminating the bias is more important than calibrating the dilution for the task at hand.
The kernel method may offer peace of mind. The demands for rejection of the PSF bias are very stringent: PSF ellipticities of 10% are not uncommon, yet state-of-the-art lensing surveys require systematic errors below 0.1%, so the method we choose must reduce the PSF bias by over 2 orders of magnitude. If the kernel method can make the PSF truly round, then symmetry principles preclude any artificial coherent ellipticity. We flesh out this statement in §7. We may have more trust in the symmetry principle than we do in an analytic correction, especially when finite sampling is taken into account. Note, however, that the symmetry principle requires that the noise power spectrum is also isotropic; §8 below demonstrates how anisotropic noise can bias shape measurements through centroid bias.
We therefore will describe an alternative procedure: eliminate the bias with a kernel correction; measure the shapes, either isotropizing the noise or making a correction for centroid bias; then correct them for dilution with a heuristic formula such as Equation (6-1). In §7 we describe a rounding-kernel methodology, and in Appendix C we derive a higher-order version of the dilution correction. These methods were used by Fischer et al. (2000), Smith et al. (2000), and the forthcoming CTIO lens survey.

Pixellated Data
Real astronomical data is integrated into pixels and sampled at finite intervals, and the continuum limit that we have assumed in our analysis is not strictly valid. A formal description of the process is that the PSF is altered by convolution with the pixel response function (PRF) to produce the effective PSF (ePSF), which is then sampled at the pixel pitch a-see Lauer (1999) or Bernstein (2001) for further elaboration. If the pixel pitch is coarser than the Nyquist interval of the PSF (λ/2D for a diffraction-limited image), then there is aliasing and we cannot unambiguously reconstruct the original image or the true b pq coefficients.
In the case where the pixel size a is small compared to the PSF size, the b pq may be estimated by converting the integral Equation (6-16) into a sum over pixels. But the formally correct procedure is to fit the vector of sampled pixel fluxes I i with a model galaxy with free b o pq . This is a standard minimization problem for the χ 2 parameter Since the model is linear, the solution for b o pq and its covariance matrix have a closed form, but of course the number of b o pq coefficients allowed in the model must be less than the number of pixels being fit. With the best-fit coefficients and their covariance matrix in hand, we may proceed with the methods outlined above. It is the ePSF rather than the PSF that has convolved the source image, but since we in fact measure the ePSF for stars, the deconvolution procedures are unchanged. Dithered exposures can be handled by extending the pixel sum over multiple exposures, as long as the PSF is unchanging. Poorly sampled images will manifest the aliasing as strong degeneracies in the solution for the b pq .
A potential time-saving procedure would be to fit the pixel data to a Laguerre expansion that is already convolved with the local PSF, minimizing Here C is the Laguerre-coefficient convolution matrix for the local PSF. The intrinsic coefficients b i pq are the unknowns in this linear fit. This is very reminiscent of the method of Kuijken (1999), except that our model for the intrinsic galaxy shape may depart from circular symmetry, and our Laguerre formalism will allow the fit to proceed as an matrix solution, with matrices rapidly built by recursion. Aside from its suitably to poorly-sampled data, this direct-fitting approach has two further advantages over the deconvolution-matrix method of the previous section. First, when the noise is not flat and Equation (6-17) does not apply, the direct-fitting method more directly produces the full covariance matrix for the b i vector. Second, the direct fit may make use of images that are partially contaminated by invalid pixels, e.g. cosmic rays, since they may be excluded from the χ 2 sum. It is this method which we consider most promising for real data.

Rounding-Kernel Methods
In this section we develop a method for producing a convolution kernel that symmetrizes the PSF to some desired degree. As discussed above, this is a potentially efficient means of reducing or eliminating PSF bias from the shear determination. By taking advantage of the Laguerre decomposition, the derivation and application of the spatially-dependent kernel can be efficiently implemented.

The Target Transfer Function
As discussed in §6.2, the effect of seeing is to convolve the initial image, I i , with a PSF transfer function T to produce an observed image, I o . One could theoretically remove all the effects of seeing by convolving the observed image with a kernel K whose Fourier transformK = 1/T to produce a final image, I f .
However, as discussed above, this is impossible to carry out in practice, becauseT (k) = 0 for all k above some critical value. We can avoid this problem by ignoring the high order moments of I o and I f . This is straightforward when one does an eigenfunction expansion of the images using the eigenfunctions introduced in §6.3, since the expansion can be truncated at some value of N = p + q. This captures most of the real information about the PSF without introducing the noise from the higher order moments. After the convolution, we then have a target transfer function, T ′ , which is not exactly a delta function, but which can be made to match a delta function up to some order N . Most generally we have For the ideal case of Thus, if T ′ satisfies Equation (7-5) up to some cutoff order, N , I f will be approximately identical to I i up to the same order. If we can find a kernel K, which produces this target transfer function, T ′ , we will be able to remove all the effects of seeing as well as possible given our ignorance of the high order terms.
How stringently do we need to satisfy Equation (7-5)? Less strict requirements on T ′ make it easier to find an appropriate, compact kernel. Our present goal is to create a transfer function T ′ which does not produce any shear bias. If the original scene I i is unlensed, and we represent the shear measurement process as some operator δ(I), then the isotropy of the Universe guarantees that δ(I i ) = 0 (we consider δ as a complex number δ + + iδ × ). Our demand on the transfer function is that Most generally the results of the shear measurement process can be expanded as a power series in the coefficients of the T ′ PSF: where the a i are some coefficients that depend upon the measurement process, the image characteristics and the size of the PSF. We now examine the consequence of rotating the image and the PSF by some angle β. The measured shear and the T ′ coefficients behave as The individual galaxies in the original scene I i are not invariant under rotation, but any statistical measure of their collective properties must be invariant under rotation, i.e. I i is invariant under rotation in the same sense that the Universe is isotropic. Therefore the a i coefficients of Equation (7-7) are unaffected by the rotation, and to satisfy the conditions (7-8) we must have a p 1 q 1 = 0 for p 1 − q 1 = 2; a p 1 q 1 ,p 2 q 2 = 0 for p 1 − q 1 + p 2 − q 2 = 2; . . . (7-9) Thus only PSF terms with m = p − q = 2 can cause a shear bias, to first order. The primary goal of our kernel, therefore, will be to set b ⋆′ pq = 0 (m = p − q = 2) (7-10) With this condition satisfied, shear bias can only be of order (b ⋆ pq ) 2 where p = q. Higher-order a's are non-zero only if m 1 + m 2 + · · · = 2. Many of these elements are zero as well; for example by considering the invariance of δ(I) under infinitesimal translation, one can demonstrate that a 10,21 vanishes, but a 10,10 must exist. Satisfying (7-10) does not, therefore, guarantee the elimination of shear bias to all orders.
For absolute assurance that a shear bias is absent, we can set b ⋆′ pq = 0 (p − q = 0 [mod 4]), (7-11) i.e. insure that the PSF has fourfold symmetry. In this case the condition m i = 2 can never be satisfied for non-vanishing coefficients of T ′ . Note that it would also suffice to enforce any mfold rotational symmetry on the PSF beyond m = 2, e.g. the diffraction pattern of a triangular secondary-support structure cannot by itself cause shear bias.
In practice we produce a kernel K which enforces condition (7-10) up to some order N . We can set b ⋆′ 10 = 0 by appropriate choice of PSF center. The remaining shear biases must be of order (b 21 /b 00 ) 2 , b 41 b 01 /b 2 00 , etc., which are generally quite small.

The Components of the Kernel
Since we are going to convolve the kernel with a pixellated image, we must construct the kernel as a 2d array instead of a continuous function. The simplest kernel is the identity kernel, composed of an array of 0's with a single 1 in the middle. This kernel conserves flux; we want a flux-conserving kernels, so we only consider kernels which are the identity kernel plus a kernel for which the elements sum to zero.
Next note that taking the derivative of an image is usually approximated by a discrete difference. For example: Another way to write this equation is as a convolution: In fact, all partial derivatives of any order in x and y can be represented as a convolution by a zero-sum kernel. First and second order derivatives can be effected using 3x3 kernels; third and fourth order derivatives require 5x5 kernels, and so on. One can therefore think of the kernel as being made up of a sum of these derivatives (∂/∂x) i (∂/∂y) j including, of course, the identity kernel, (∂/∂x) 0 (∂/∂y) 0 .
The eigenfunction expansion of our images actually suggests a slightly different set of components for the kernel. Namely, These components make it easy to use the raising and lowering operators to determine how a given kernel will act on an image. Note, however, that D ij are complex, which means that to end up with a real image after the convolution, we require that k ji =k ij .
The 3x3 kernel components are: We note at this point that the FT97 method is equivalent to use of just the D 20 and D 02 kernel elements. The 5x5 and larger kernel components can be found by convolving the 3x3 components. For instance, D 30 = D 10 • D 20 , and D 41 = D 11 • D 30 .
The problem with this method is that b ⋆′ will not usually be completely specified. Neither Equation (7-10) nor (7-11) fully constrains b ⋆′ . The easiest way to deal with this is to set k pq = 0 for each unspecified b ⋆′ pq . Then the kernel will be as simple as possible while still satisfying all of the requirements for b ⋆′ . A more sophisticated technique for dealing with this issue is described below in §7.4.

The Kernel with Pixelization
The same general technique applies to the case of finite pixels in that we want to solve Equation (7-25) for k. The only difference is the calculation of D ij b ⋆ . In this case, we must use Equations (7-18) -(7-22) rather than Equation (7-17).
Calculations for the other D ij are similar.

The Kernel in a Distorted Frame of Reference
Most telescopes, especially those with large fields of view, have fairly significant distortion. To deal with this correctly, the shape measurements should be made in the undistorted world coordinates rather than the chip coordinates. The kernel's pixel grid, however, is still is in the distorted coordinate system. The calculation of D ij b ⋆ must therefore take into account the two different coordinate systems.
If the world coordinates are (u, v) and the pixel coordinates are (x, y) then the correct values for z1 and z2 are:

Minimizing the Noise from the Kernel Convolution
The motivation for having a compact kernel is to minimize the noise added by the convolution. A large kernel will use data with significant noise but little signal, adding to the noise in the convolved image.
Therefore, let's consider the noise in the convolved image. Define K(m, n) to be the total convolution mask.
If the noise in the image is dominated by sky noise, then we can define n 2 to be the variance in each pixel of the observed image. The noise in each pixel of the convolved image is In §7.3.1, we set some kernel coefficients to 0 to account for the unspecified components of b ⋆′ . Really, each of these coefficients is arbitrary, so the solution to Equation (7-25) will be k = k 0 + AY  where k 0 is a specific solution, A is a matrix, and Y is an arbitrary vector. The dimension of Y is equal to the number of unspecified components of b ⋆′ .
We can then find the particular Y which minimizes the noise added to the image by the convolution. For this derivation D(m, n) is the vector of each D ij (m, n). So, K(m, n) = D(m, n)k. This is now a matrix equation which can be solved for Y, which then gives you the solution for k which minimizes the noise.
It turns out that if this method is implemented exactly as described, one ends up with a fairly large kernel, which is essentially a smoothing filter. While this solution will work, it is not exactly what we want from the kernel. We would rather have a smaller kernel which minimally changes the image. So, rather than minimize the noise in I f , we minimize the noise in I f − I o . In other words, we leave out the D 00 term in the vector products of Equation . This results in fairly compact kernels which vary smoothly across the image.

Additional Kernel Components
There are 9 free parameters in a general real 3 × 3 kernel. Equations (7-18) -(7-22) define 5 kernels. The identity kernel is another. Thus, there are 3 more independent 3x3 kernels we could construct. Two of these can be made to approximate discrete versions of the derivatives D 10 and D 01 . The other can be made to approximate a discrete version of D 11 . These alternate versions are: AltD 01 = AltD 10 (7-38) We can construct alternate higher order kernels in the same way as we constructed the regular higher order kernels. Namely, AltD 30 = AltD 10 • D 20 , AltD 41 = AltD 11 • D 30 , etc.
These extra kernel components are useful if one is minimizing the noise as described above, since they add extra degrees of freedom for the minimization, and therefore can result in a smaller, less noisy kernel.

Interpolating Across an Image
The above discussion explains how to find the appropriate kernel given a particular PSF. However, in a real image, the PSF varies across the chip. Thus, the kernel will also vary across the chip.
There are two potential ways to deal with this. One can find the appropriate kernel for each star in the image. Then fit the kernel components k ij as functions of (x, y). Alternatively, one can fit the coefficients b ⋆ pq of the measured PSF decomposition as functions of (x, y) and then solve for the appropriate kernel at each point.
We choose the first method in our analysis for two reasons. First, the derived kernel can be directly applied to each star to make sure that it really does make the star round. Occasionally, a star will have significant high order components due to crowding or an uncorrected cosmic ray. When this happens, the derivation above fails, and we reject this kernel from the fit. It is a little cleaner to recognize these outliers with kernel interpolation than with PSF interpolation.
The second reason to prefer fitting the kernel rather than the PSF is computational efficiency. The appropriate kernel must be calculated at pixel. It is significantly faster to evaluate a function than to solve a matrix equation. On a 2K x 4K chip, there are 8 million kernel evaluations. Both methods gain significantly by using a locally linear approximation to the spatial variation, but there is still a significant difference in computation for the two methods.
The kernel scheme described here allows an substantial speed gain over a typical convolution method. Fourier methods are fastest for large convolutions, but are not practical for spatiallyvarying kernels. In our scheme, the output image can be written as (7-40) Instead of calculating the kernel at each output pixel, we can instead produce the images D ij I o (which are constant-kernel convolutions), and then accumulate sums of these images with spatially varying weights k ij (x, y). More importantly, the D ij I o can be produced by recursive application of the 3 × 3 kernels. Hence the convolution by a 7 × 7 kernel can, for example, be reduced to three successive applications of 3 × 3 kernels.
Note that interpolation of the PSF elements across the image is required for the analytic deconvolution described in §6.3.5.

Dilution Correction
Once a kernel has been applied to the image to symmetrize the PSF, the measured galaxy ellipticities need to be scaled by some resolution factor R to account for the PSF dilution. Appendix C describes the scheme we use for estimating R, which closely resembles the KSB "smear polarizability" derivation.

Centroid and Selection Biases
Even if the asymmetries of the PSF have been perfectly removed by a deconvolution or other method, there are two effects which can cause the estimate of the mean shape to be biased in the direction of the original PSF orientation. The first is a selection bias, first noted by K00. The second is a measurement bias which arises in the presence of anisotropic noise or anisotropic PSF, which we term the centroid bias. K00 discusses an effect whereby the errors in centroid appear in second order as biases in the ellipticity, and concludes that such errors are probably negligible. We show below that this bias is significant, and in fact just one of a class of noise-rectification biases that can occur.

The Selection Bias
Kaiser's selection bias operates as follows in the presence of an anisotropic PSF: if the PSF is elongated in the x direction (e ⋆ + > 0), then objects with intrinsic shape e + < 0 cover a larger area after PSF convolution than do objects with intrinsic e + > 0. On the observed image, therefore, such objects have both lower surface brightness and lower significance ν. As most detection algorithms involve some cut in surface brightness or ν, the detected population will be biased toward e + > 0.
The mean e + of the population will be biased, therefore, even if all the detected objects are perfectly corrected for the PSF anisotropy. K00 demonstrates that the bias will scale roughly as e ⋆ σ 2 ⋆ /σ 2 i ν 2 , where ν is the detected significance.
The selection bias may be defeated by careful definition of the sample of target galaxies. The key is to produce some significance statisticν which is independent of the shape of the object. For an image of flux f covering area A in an image with white noise density n, the significance is normally ν ∝ f / √ nA. But the observed area A is shape-dependent in an anisotropic fashion if the PSF is elliptical. We may instead defineν ∝ f / nÃ, whereÃ is the object's area on a version of the image which has been deconvolved, or at least had the PSF rounded by a convolution kernel, as described in the previous section. We then define the target sample by some cutoff limitν min to this "isotropic" significance.
Detection algorithms generally have some statistic which is used as a cutoff. For KSB's imcat, this is ν. For SExtractor and FOCAS, this is the number of pixels N above some isophotal threshold in a filtered version of the image. To eliminate the selection bias, we make a scatter plot of the detection statistic, e.g. N , vs. the "isotropic" significanceν. The selection bias is eliminated if we choose our cutoffν 0 sufficiently high that no members of the selected population approach the cutoff N min of the detection statistic. This is illustrated in Figure 7. In this way we insure that our selected population is free of any anisotropic selection criterion that may have been created by the original detection algorithm. the criteria necessary to avoid the selection bias. To avoid selection bias, we want the galaxy selection to be entirely on the basis of the shape-independent statisticν, and not influenced by the potentially shape-biased SExtractor threshold N min = 5 [some galaxies with N < 5 are present due to object splitting]. We must therefore set our ν min threshold at ≈ 10, which is the lowest value for which the galaxy population does not extend to N < N min . Galaxies in the shaded region are therefore excluded.

Centroid Bias
The centroid bias, discussed in K00, can be qualitatively explained as follows: an error δx in centroid determination along the x axis inflates the measured second moment I xx by some amount ∝ δx 2 , whereas a y-axis error inflates I yy . If the centroid errors are isotropic, the mean effect upon e + is zero. But if the x centroid errors δx 2 exceed those in y, there is a net tendency to measure a positive e + . If the PSF has e ⋆ + > 0, then the PSF-and consequently the mean galaxy image-is more extended in the x direction, and hence centroids are less accurate in x than in y. The noise in centroid estimation therefore pushes the mean measured e in the direction of e ⋆ . The centroid errors scale as ν −1 , and therefore the centroid bias in e scales roughly as e ⋆ ν −2 . Faint galaxies at low ν have a much larger centroid bias than the bright, high-ν stars used as PSF templates, so the effect is not properly removed by the KSB "smear polarizability" corrections.
Deconvolving the image eliminates the PSF anisotropy but does not eliminate the centroid bias, because the deconvolved image will have anisotropic noise. There will be more noise power along k x than k y after deconvolution, hence there will still be anisotropic centroid errors and a bias toward the original e ⋆ .
The same situation arises if we apply a convolution kernel to symmetrize the PSF. This kernel will smooth the image slightly in the direction perpendicular to the PSF, reducing the noise level in that direction somewhat. The convolved image will therefore again have a higher centroid error along the original PSF axis, and a consequent bias in e.
A quantitative understanding of the centroid bias can be gained from Equation (6-53), the second-order expression for η in the simultaneous solution for shape and centroid (and possibly size). Let us presume that each b pq has an true value plus some measurement noise δb pq . The measurement noise has mean value of zero since b pq is a linear function of the intensity. If the object is intrinsically round, centered, and the weight is properly sized, then we have b 00 = b 10 = b 20 = 0, and we would indeed measure η = 0. In the presence of noise, however, the second-order terms in Equation (6-53)-for instance the term ∝ b 2 10 /b 2 00 in the second line-may have non-zero mean values. If b 2 10 = Var(b 10 ) = 0, there will be a non-zero η .
Note that there are a number of second-order terms in the η solution, arising not just from the solution for object center, but also from the solutions for object size and ellipticity. We apply the term "centroid bias" to all of these effects.
According to Equation (6-17), Var(b 10 ) = 0 in the presence of a white noise spectrum P (k) = n (recall that b 10 is a complex number). If we have a PSF of size σ ⋆ and ellipticity e ⋆ along the x axis, we may attempt to round out the kernel by smoothing along the y axis a little bit. This will produce a noise power spectrum P ′ (k) ≈ n(1 − 2e ⋆ σ 2 ⋆ k 2 y ). The k 2 y term in the noise power produces a non-zero value for Var(b 10 ), and also non-zero covariances between all the other second-order elements of Equation (6-53). The leading term in the expression for centroid bias takes the form where K is a constant, e o is the ellipticity measured on the convolved image, σ o is the size observed on this image, and R is the resolution factor of §6. Note that the bias is in the measured shape, before correction for PSF dilution, which would add a factor R in the denominator. The value of K depends upon some higher-order moments of the typical galaxy, and upon the details of the PSF-correction and measurement procedure-in some cases K can be negative. Note also that the functional dependence of the centroid bias is essentially the same as that for the selection bias, hence they are difficult to distinguish. Figure 8 demonstrates the existence of the centroid bias in a very numerical simulation, in which we convolve a circular Gaussian source, with an elliptical Gaussian weight, then measure the ellipticity and centroid with a fixed-size circular Gaussian weight. It is clear that the mean measured ellipticity depends upon the significance (i.e. S/N level) as described by Equation (8-1) with K = −2. If the weight size is iterated to maximize significance, the bias increases to K ≈ −6. For other measurement algorithms or galaxy shapes, K will differ, so we must determine K empirically.
The bias in Equation (8-1) will be equal to the shape-noise uncertainty in the mean of N galaxies when For typical PSF ellipticities and K values, with galaxies somewhat resolved (R ≈ 0.5) and detected at ν = 10, we find the bias equal to the shape noise when N ∼ 10 4 . The effect clearly cannot be ignored in present-day surveys. If the PSF ellipticity varies, then the bias will induce false power into the distortion power spectrum, with a total fluctuation power of ≈ (e ⋆ ) 2 K 2 /ν 4 , for R ≈ 0.5. If we take the RMS value of e ⋆ to be about 3%, take |K| ≈ 5, and ν = 8, we see that the bias power is about 5% of the typical cosmic signal, for which the RMS distortion is ≈ 1%. In unfavorable parts of the power spectrum, the ratio of PSF bias power to cosmic signal will be worse; hence the bias cannot be ignored for cosmic shear studies which aim to move beyond mere detection to true precision measurements.
There are several possible strategies for defeating the centroid bias. The simplest is to find K empirically, then apply an appropriate bias correction to each object using Equation (8-1). Another approach, useful if one has applied a PSF-rounding kernel to the image, is to add noise back to the image to recreate the original isotropic noise spectrum. If both the PSF and the noise spectrum are isotropic, there is no way the mean shape can be biased. Unfortunately, creating the appropriate noise field to add to the image is computationally expensive for all but the simplest convolution kernels.
Symmetrizing the noise is easier if we have used the Laguerre decomposition method to de- convolve the images. If we have properly propagated the original covariance matrix for the b pq , then we know Var(b 10 )-and all the other relevant covariances-in the deconvolved image. We can add noise to the deconvolved b pq elements in order to zero out the asymmetric elements of the covariance matrix.

Procedures for Shape Measurement with PSF Correction
We now have all the necessary tools for a procedure for measuring shapes and shear in the presence of a PSF convolution. Figure 9 is a flowchart for this procedure, and we elaborate on each step below. First note that there are two branch points. The first is the decision whether to measure the shapes on a summed image, or to measure PSF-corrected moment information from each exposure and combine the moments. The former is easier if all exposures cover the same sky area, but the latter is recommended whenever the coverage is interlaced on the sky. The second branch point is deciding whether to use an analytic or a kernel correction for PSF bias. The remainder of this section delineates the overall data reduction procedure for weak lensing measurements.

Exposure Processing
The steps listed in this section must be performed on each exposure. For mosaic detectors, these steps should be done separately for each channel of the detector if the PSF may change abruptly at the boundaries between detectors.
1. Bias Subtraction and Flat-Fielding: These can be done in the usual fashion. Note that field distortion leads to significant variation in pixel area in many CCD cameras. If the flat-field image is, as usual, obtained by exposing a source of uniform brightness, then the flattened pixel data are a properly calibrated map of the surface brightness of the sky rather than a map of the flux from the sky. This is actually what we want, though it means that simple aperture photometry will lead to incorrect flux estimates.
2. Object Detection: Available detection packages such as FOCAS (Valdes 1982), SExtractor (Bertin & Arnouts 1996), or ProFit (P. Fischer, private communication) can be used to identify all objects on the exposure with significance ν 3 and produce a catalog with preliminary position, size, and ellipticity estimation. Each of these packages also produces a useful estimate of total magnitude, which we will preserve. The choice of object detection package is not critical, because the objects that will have useful shape information have sufficiently high significance that any decent detection algorithm will work.
3. Preliminary Measurement: Our shape-measurement algorithm is applied to each detected object in order to obtain an accurate Gaussian-weighted estimate of centroid, size (i.e. the σ  Fig. 9.-The data flow from raw images through a final lensing distortion map or statistic is illustrated here, with each step described in §9. Shaded regions represent steps that may not be used in all circumstances. Implementation details will be presented in Jarvis et al. (2002) which maximizes significance), and shape. Objects with saturated or bad pixels are flagged at this point.

Registration:
The map from pixel coordinates x p to the global celestial coordinate grid x is established by fitting the measured centroid of the bright objects to a collective catalog (e.g. the average of all individual exposures' catalogs) and/or an astrometric catalog. This step is critical as even slight misregistrations will produce a systematic coherent ellipticity to the summed images. The photometric offsets between images are determined at this point as well.
5. Image Summation (optional): The exposures can now be remapped to a common grid and photometric scale and summed to give a deeper image. If the flattened images are surfacebrightness maps as described in Step 1, then a simple interpolation can be done when remapping the images-no Jacobian flux corrections are needed. Recall the caveats about image combination in §4.1. If shapes are to be measured from the summed image, then we must at this point reconstruct the preliminary catalog by executing Steps 2 and 3 on the summed image.
6. Identification of Stars: As many unsaturated stellar images as possible should be identified on the image for use in PSF determination. Stars are typically identified on the size-magnitude plane. This is a difficult task to execute without excessive human intervention; we want our star-finding algorithm to be flexible enough to identify stars successfully in the presence of a spatially varying PSF, but strict enough to insure that galaxies are not counted as stars. Failure on either count will lead to an erroneous estimate of the PSF size and shape on some part of the exposure, which will lead to a false feature in the shear field. Our algorithms for star identification are described in Jarvis (2002). The algorithms are particularly suited to identifying stars in the presence of a spatially varying PSF size.
7. PSF Measurement: The shape-measurement algorithms are now applied to the identified stars. The integrals in Equation (6-16) give b ⋆ at the location of each star. Note that the integral can readily be executed in the global coordinate system x because we know the map from pixel coordinates, and because pixel data are already in terms of surface brightness I. Thus the effects of image distortion are removed at this step.
8. PSF Interpolation: With the PSF b ⋆ measured at the locations of the stars, we need to fit a function b ⋆ (x) which describes the PSF at any location on the image. We use a polynomial to describe the variation of each b ⋆ pq across the image. Note that for complex or undersampled PSFs (such as those on HST), interpolation of the b ⋆ pq components is much easier and more accurate than interpolating a pixel-wise representation. In any case this step is again critical: the interpolation scheme must be flexible enough to follow PSF variation, but must remain well-behaved at the image edges and in regions where PSF stars are sparse. There must also be some form of outlier rejection for PSF stars that have cosmic rays or near neighbors contaminating the measurement, but a true excursion of the PSF in some part of the image must not be rejected.
9. PSF-Rounding Kernel (optional): At this point we derive and apply the kernel to remove the anisotropy from the PSF, if desired. Our method for doing so is described in §7. After the convolution it is necessary to repeat Steps 7 and 8 so as to have an updated map of the PSF.

Object Measurement
We next outline the procedures which yield a shape estimate for every object in the field.
10. Creation of Target Catalog: A list of all targets for shape measurement must be compiled. If there are 5 exposures of each object, then a master target list can be produced by taking the union of the individual exposure catalogs produced in Step 2. One can use a very low detection threshold on the individual exposures, then guard against noise detections by demanding coincidence on 2 or more exposures. If there are 5 exposures per object, then it will be necessary to run a detection algorithm on a summed image in order to find all the potentially useful target galaxies. The target catalog should include an estimate of the centroid (in global coordinates) and some observed size s o of each object. To avoid Kaiser's selection bias ( §8), the criterion for acceptance into the catalog must be some shapeindependent statistic.
The following steps are performed for each object in the target catalog: 11. Measurement of Observed Moments: Given the global coordinates of the object, one can determine all of the exposures on which it should appear. Note that for low-significance objects, this may include exposures for which the object was missed in the preliminary detection of Step 2. We then use Equation (6-16) to measure the Laguerre expansion b o of the image as observed on each exposure. Once again the integrals are performed in the global coordinate system to remove the effects of optical distortions. The integrals require choice of a centroid and a size parameter σ o for the basis functions ψ σ pq . We may use the approximate centroid from the target catalog, and set the weight size σ o equal to the typical observed object size s o from the target catalog. The covariance matrix for the b o vector has the diagonal form given in Equation (6-17) for sky-dominated galaxies. Measurements contaminated by saturated pixels or other defects may be rejected at this step.
12. Correction for PSF Effects: There are two methods available here. The first sequence may be used if the rounding kernel has been applied to remove the PSF bias. In this case, if the significance on the individual exposure is 3, we may proceed as follows: (a) Shape Measurement: The image can be shifted, dilated, and sheared until it passes the centroid, significance, and roundness criteria embodied by Equations (6-19)-(6-21). This yields an optimal measure of the observed ellipticity e o . If the object is of significance ν 0 3 on the individual exposure, then we can solve for centroid and size on each exposure independently and still ignore the higher-order terms in the uncertainty of Equation (4-1).
(b) Centroid Bias Correction: Each measured ellipticity must be corrected for the centroid bias ( §8) using an empirical value of K in Equation (8-1).
(c) Dilution Correction: The PSF is interpolated to the position of this object to determine the resolution parameter R (Equation [6-1]) or the higher-order version derived in Appendix C. The correction for dilution to give the image-plane (pre-seeing) ellipticity is then simply e i = e o /R. The uncertainty ellipse for e (Equation ) is therefore also scaled by 1/R. Note that it is possible to obtain e i > 1 if the noise makes the target appear smaller than the PSF.
(d) Averaging of Exposures: The e i from the collection of exposures are averaged using the weighting procedures of §4.2. Some form of outlier rejection is necessary to remove objects contaminated by cosmic rays or other defects.
If we are to perform the PSF bias correction analytically, or if the significance per exposure is low, then it is best to average the deconvolved moments rather than deriving e i for each exposure: (a) Deconvolution of Moments: The PSF decomposition b ⋆ is interpolated to the position of the object, and the form of the convolution matrix C(b ⋆ ) is calculated using the coefficients in Appendix B. This matrix is truncated at some order p + q ≤ N and inverted; the deconvolved (pre-seeing) Laguerre decomposition is b i = C −1 · b o . Since this is a linear operation on b, the covariance matrix for b i can be propagated from the simple diagonal covariance matrix for b o (Equation [6][7][8][9][10][11][12][13][14][15][16][17]). There is a subtlety involved in the choice of weight scale σ o : in §6.2 we determined that the ideal weight scale for the deconvolved image is σ 2 i = s 2 i + σ 2 ⋆ = s 2 o , where s i and s o are the sizes of the pre-and post-seeing objects, and σ ⋆ is the size of the PSF. The typical value of s o was placed in the target catalog in Step 10. According to Equation (6-63), our deconvolution formulae are simplest if we choose the weight scale by (9-1)

Thus in
Step 11 we in fact want to use a weight scale somewhat larger than the s o which maximizes the significance of detection.
(b) Combination of Moments: From each exposure we have estimated a b i deconvolved moment vector, with known covariance matrix. We average these vectors to obtain a single best estimate of b i . We have a choice of weights to apply in producing the average; an obvious choice is to weight each exposure inversely with Cov(b i 20b i 20 ), since b 20 carries most of the shape information. Again some form of outlier rejection is necessary at this step.
(c) Symmetrization of Noise: Certain elements of the propagated covariance matrix for the b vector must be zero in order to avoid noise anisotropies that produce centroid bias ( §8). This can be achieved by selectively adding Gaussian noise to various elements of the combined moment vector.
(d) Determination of Shape: Given the average of b i from all exposures, we find the translation and shear that must be applied to satisfy the centroid and roundness conditions b 10 = b 20 = 0. The formulae of Appendix A are used for this. The covariance matrix for b i can be propagated through the transformations as well. The uncertainty in the shape η is then the square root of Var(b 20b20 )/(b 2 00 ). Note that we do not want to maximize the significance by dilating to set b 11 = 0. Our optimization criterion is that the shape have minimal error after our transformations. One could dilate the image to satisfy this desire, but our choice of σ i should already have us close to the optimum, according to the arguments of §6.2.
After either of these procedures, we have a measure of the deconvolved shape along with its uncertainty.
13. Combination of Different Wavelengths: If we have imaged the field in a variety of filters, then we will have obtained a shape measurement in each filter. The galaxy's shape and moments may depend upon wavelength so we don't want to average together images or moments measured in different filters. We can, however, use the methods of §4.2 to produce a wavelength-averaged e i which is maximally sensitive to shear. Weighting each filter by the error in its measured shape insures that we obtain the most sensitivity from each galaxy regardless of its color.
After completion of all these steps, we have a catalog of all objects in the field, specifying their location, magnitude, optimal shape measurement, and shape uncertainty.

Determination of Shear
With the shape catalog in hand we are close to the scientific goals. The remaining step is 14. Generation of Shear Data: The target galaxies are binned by position, etc., into subsets for which we wish to determine a shear. The shapes (and their uncertainties) may have to be rotated, e.g. into tangential coordinates about some mass center, depending upon the shear statistic under study. The formulae of §5 take the collection of shapes and uncertainties and allow us to create an optimal shear estimate, as well as to propagate the uncertainties in shape to the shear measurement.

Conclusions
We have attempted to produce, as rigorously as possible, an end-to-end methodology for measuring gravitational distortion that has optimally low noise, with calibration and noise levels derivable entirely from the observations themselves. We have succeeded in many, but not all, aspects of the problem: • The measurement of individual galaxy shapes appears to be optimal and has traceable noise characteristics; indeed there even appears to be a straightforward way to handle undersampled data and retain proper covariance information for the intrinsic Laguerre coefficients b i pq ( §6.4). The Gaussian weights underlying the Laguerre decomposition are nearly optimal for skydominated exponential-profile galaxies. We may wish to reexamine this scheme for the case of when the galaxy is brighter than the sky background.
• The correction of measured moments for the distorting and diluting effects of the PSF can be effected to arbitrary accuracy using the Laguerre-decomposition methods. There will be a tradeoff between elimination of systematic errors-which pushes toward inclusion of higher-order terms in the deconvolution-and the minimization of measurement noise from high-frequency terms. It is not clear if the Laguerre method is optimal with, for example, Airy PSFs with sharp cutoffs in k space; but the method should be better than those yet applied.
• We have identified methods to work around two sources of bias that arise from PSF ellipticities even in the presence of perfect deconvolution.
• Measurements of galaxy shapes from different exposures or filter bands can be optimally combined with standard least-squares techniques, since we know the e uncertainty from each individual exposure.
• The determination of lensing distortion from the ensemble of galaxy shapes has a straightforward, optimal, calibratable solution in the case of no measurement noise. In the presence of measurement noise, however, it is necessary to have some knowledge of the noiseless shape distribution to get the calibration factor exactly correct. Our approximations, however, seem to suffice to get accuracies of 5-10%.
A detailed performance comparison of our methods with other authors' is beyond the scope of the paper. We can guess, however, that the reduction of measurement noise relative to a carefully weighted implementation of KSB will be slight, perhaps a factor of 1.5. But our methods, like those of K00 and Refregier (2001), are formally valid for any reasonable PSF, and hence we expect to have much-reduced systematic errors. Indeed it is likely that the biases of §8 have not been extensively tested because they were lost under the larger errors in PSF correction (as also noted by K00).
Our methods share elements with many previous proposals. An aesthetic difference is that we retain the geometric meaning of ellipticity by, in effect, using adaptively shaped weights instead of fixed circular weights. This makes the "shear polarizability" a purely geometric effect. As a consequence we can examine the P (e) distributions and find ways to exploit surface brightness, b 22 , or color information to separate the spheroid galaxies and weight them more heavily to reduce the shape noise.
It will be of interest to see how our method compares to the commutator method of K00. One would hope that the two independent methods could be applied to the same dataset and yield the same results, bolstering our confidence in these very difficult measurements.
In a succeeding paper  we will present some of the implementation details for the analytical methods here, and test the methods on real and simulated data. Fischer et al. (2000); Smith et al. (2000); Wittman et al. (2000) make use of portions of this methodology, so the systematic-error tests and calibration tests in those papers already serve as demonstrations. Upcoming precision measurements of cosmic shear will make even greater demands upon the systematic-error reduction and accurate calibration that our methods offer. This work was supported by grant AST-9624592 from the National Science Foundation. We would like to thank Phil Fischer, Deano Smith, Tony Tyson, Jordi Miralda-Escudé, and Dave Wittman for many discussions on these methodologies, and help in implemented various imageprocessing algorithms. Thanks also to all our collaborators who have waited several years for a coherent explanation of the methods we are using.
Note that the kurtosis measure a 4 defined in Equation (3-26) is the same as b 22 /b 00 which appears here. We make the ansatz that the correct form for R in the case of finite dilution or finite e is the Gaussian form Equation (C1) with the kurtosis term added: s 2 ≡ 1 − a 4 1 + a 4 σ 2 cosh η.
Note that we apply the kurtosis correction to the PSF size measure s 2 ⋆ in the same way as for the object, to give a well-behaved correction for poorly-resolved objects.