A&A 456, 827-838 (2006)
DOI: 10.1051/0004-6361:20054794
K. Kuijken^{1}
Leiden Observatory, PO Box 9513, 2300 RA Leiden, The Netherlands
Received 31 December 2005 / Accepted 15 June 2006
Abstract
Aims. Accurate measurement of gravitational shear from images of distant galaxies is one of the most direct ways of studying the distribution of mass in the universe. We describe a new implementation of a technique for measuring shear that is based on the shapelets formalism.
Methods. The shapelets technique describes PSF and observed images in terms of Gauss-Hermite expansions (Gaussians times polynomials). It allows the various operations that a galaxy image undergoes before being registered in a camera (gravitational shear, PSF convolution, pixelation) to be modeled in a single formalism, so that intrinsic ellipticities can be derived in a single modeling step.
Results. The resulting algorithm, and tests of it on idealized data as well as more realistic simulated images from the STEP project, are described. Results are very promising, with attained calibration accuracy better than four percent (1 percent for round PSFs) and PSF ellipticity correction better than a factor of 20. Residual calibration problems are discussed.
Key words: gravitational lensing - techniques: image processing - cosmology: dark matter
Weak gravitational lensing is recognized as a profitable way to study the dark matter distribution in the universe, and with a series of ever wider-field astronomical cameras coming online, very large weak lensing surveys are being planned and performed. The inevitable source of noise in weak lensing measurements is shape noise, caused by the diversity of projected galaxy shapes on the sky. To beat down the noise, very large numbers of galaxy ellipticities need to be averaged: thus new lensing surveys such as the CFHT Legacy Survey (Hoekstra et al. 2005), the CTIO survey (Jarvis et al. 2003), the VIRMOS-Descart survey (van Waerbeke et al. 2005) or the recently-approved Kilo-Degree Survey (KIDS) on the VLT Survey Telescope, will contain millions of background sources, which in principle should enable very accurate shear measurements.
This averaging down of the shape noise through large number statistics only makes sense if systematic errors can be controlled: the main one is still the blurring of the source images by atmosphere, telescope and detector pixels. The commonly-used Kaiser et al. (1995; henceforth KSB) and Luppino & Kaiser (1997) methods provide recipes for correcting these effects, and are very successful. Nevertheless, they are based on an idealized model of the effect of point spread function (PSF) convolution on ellipticity, and it is possible to construct plausible PSFs that it fails to correct properly (e.g., Hoekstra et al. 1998, Appendix D). Therefore it seems unlikely that the KSB recipes will deliver the factor of 10 to 100 improvement in fidelity that will be required to exploit the new surveys (Erben et al. 2001).
A number of different approaches have been put forward to improve PSF correction (Kuijken 1999; Kaiser 2000; Rhodes et al. 2000; Bernstein & Jarvis 2002; Refregier & Bacon 2003; Mandelbaum et al. 2005). In this paper we present a new technique which combines elements from most of these.
The paper is organized as follows: in Sect. 2 shear and ellipticity are defined, and the effect of one on the other. Section 3 summarizes the shapelets formalism, and describes how a shapelet description of a source and its PSF can be used to generate an ellipticity estimate that is useful for shear estimation. Section 4 substantiates the approach with idealized tests of the algorithm. In Sect. 5 a software pipeline is presented that implements the full processing chain from an astronomical image to a shear estimate. Results of applying the pipeline to test data from the STEP project are shown in Sect. 6. In Sect. 7 we compare with other techniques, and Sect. 8 gives the conclusions.
Following the usual practice, we parameterize the effect of weak
gravitational lensing on a distant source in terms of a shear
and a convergence :
the distorted image
is derived from the original I(x,y) via the
transformation
,
where
Without knowledge of the intrinsic source size, only the distortion , which affects the shape of the source, can be measured (Schneider & Seitz 1995).
We define the ellipticity of an object's image I as follows:
Let be the model image with constant-ellipticity isophotes that best approximates I. Then the ellipticity (e_{1},e_{2}) of I is defined such that a distortion of (-e_{1},-e_{2}) makes circular.
The major axis position angle
and the axis ratio q of an
elliptical source are simply related to (e_{1},e_{2}):
(2) |
As discussed by BJ02, expressing the shapes of objects in terms of
distortions e has a practical advantage: in this formulation it is
simple to calculate the response of object shapes to small
distortions. An elliptical source with ellipticity (e_{1},e_{2}) that is
sheared by a small amount (g_{1},g_{2}) can be viewed as a circular
source that is sheared twice, first by e_{i} and then by g_{i}, giving
a combined distortion matrix
(3) |
(4) |
K=e_{1}g_{1}+e_{2}g_{2} | (5) | ||
R=e_{1}g_{2}-e_{2}g_{1} | (6) | ||
(7) | |||
(8) |
Assuming now that we have an ensemble of elliptical sources, of random
orientations, so that before distortion
and
,
the average
ellipticity of the population after a distortion (g_{1},g_{2}) is simply
The shapelets basis is described in Refregier
(2003). It consists of the two-dimensional Cartesian
Gauss-Hermite functions, famous as the energy eigenstates of the 2D
quantum harmonic oscillator:
(10) |
Shapelets are a convenient basis set for describing astronomical images because of the compact way in which various operators (translation, magnification, rotation, shear) can be expressed as matrices that act on the shapelet coefficients. Shapelets have a free scale radius (the size of the Gaussian core of the functions), and R03 shows how the coefficients transform under change of , and how to convolve objects with different scale radii.
To avoid introducing a preferred direction, the expansion should be truncated in combined order N=a+b, not in a or b separately. (A basis set truncated in N is complete under rotation. Effectively, such a truncation describes an image as a product of a circular Gaussian with an inhomogeneous polynomial in x and y of order N. Rotation of such an image will mix the x^{i}y^{j} terms at constant i+j.)
The reason for choosing shapelets as a formalism for weak lensing analysis is its ability to describe the main operations that a galaxy image undergoes before it is registered at a telescope focal plane: in reverse order, pixelation, convolution with a PSF, and distortion.
In this paper we concentrate on well-sampled (PSF FWHM at least 3-4 pixels), ground-based seeing-limited images. It remains to be seen to what extent diffraction-limited PSFs, and undersampling, can be handled with this formalism.
An image that is registered on a CCD is pixelated: the flux on the
surface of the detector is read out in binned form. Mathematically,
the flux is first boxcar smoothed (i.e., convolved with a pixel), and
then sampled at a spatial frequency of once per pixel.
Therefore, if we fit a shapelet expansion to the binned image
I(k,l) as a linear superposition of the basis functions B^{ab}:
Convolution with the PSF can be expressed as multiplication with a PSF
matrix
:
if
the shapelet coefficients of the PSF are p_{ab}, and those of a
model source are m_{ab}, then their convolution is
(12) |
Given a PSF and a PSF-convolved source, both expressed as shapelet series, we determine the ellipticity of the source by modeling it as a PSF-convolved, distorted circular source of arbitrary radial brightness profile. This approach is similar to the one described in Kuijken (1999), but it is more effective when expressed in terms of shapelets.
A circular source of arbitrary radial brightness profile can be
written as a series of circular shapelets C^{n} of the form
,
where the c_{n} are free
coefficients. The C^{n} (see Appendix) are normalized to have unit
integral over x,y, so c_{n} gives the total counts in each
component. After distortion of such a source, it becomes an elliptical
source which, to leading order in ellipticity e, can be written as
This machinery allows us to write the model for the observed source as
To improve the accuracy, we make two modifications: we add centroid
error parameters, and we only fit the model to a subset of the
shapelet expansions. The centroid error parameters are included to
allow for a mismatch between the center of the object and its shapelet
expansion. If the PSF and/or the galaxy have some lopsidedness to
them, the center of their best-fit shapelet expansion may not be at
the flux-weighted center of the source (since our centering technique
simply requires the 01 and 10 components to be exactly zero). Hence
the centroid may move under convolution, which would spoil the fit. To
guard against this, instead of fitting Eq. (15) we fit a
model of the form
The second modification is made to contain truncation errors. While the shapelet basis is complete, and hence can describe any source given enough terms, in practice the fact that the source is only sampled in a finite set of pixels means that the expansion needs to be truncated. Hence, except in very special cases, the PSF and galaxy are not described perfectly by a truncated shapelet series. The missing information propagates through the analysis, and is a source of systematic error in the PSF convolution (as some PSF terms may be missing) and in the calculation of the action of shear (since shearing high-order shapelets generates also lower-order terms).
Truncation effects can be seen most clearly if we re-express the shapelets in polar coordinates (they become Gaussians times Laguerre polynomials of r - see BJ02, R03, Massey & Refregier 2005). Polar shapelets are combinations of Cartesian shapelets B^{ab} of the same order N=a+b, whose angular dependence is a pure sine or cosine of , for angular order m. The order of the Laguerre polynomial is N, with and N+m even.
The translation and shear operators, when applied to a polar shapelet of order (N,m), generate terms at order and , respectively. If we truncate the shapelet series of the best-fit circular model for the pre-seeing, pre-shear galaxy at order , then to be consistent the m=1 and m=2 series must be truncated at order and , respectively.
Figure 1: The information used in the ellipticity determination, for shapelet expansion to order N=8. On the left, the Cartesian shapelet coefficients that are fitted to describe a source and its associated PSF. On the right, the same information has been rearranged into a polar shapelet expansion (the two may be transformed into one another by appropriate mixing of the terms at order ). The heavy line indicates the polar shapelets from which the ellipticity is estimated by a fit to Eq. (16). Note the safety margin at order N-1 and N, and the fact that only azimuthal orders between -2 and +2 are fitted. | |
Open with DEXTER |
Note that we are never completely safe from truncation effects: in particular, complex PSFs can in principle mix coefficients of all orders. The problem is minimized, but not completely eliminated, by adopting suitable scale radii so that the amount of information carried in the high-order coefficients is small. We further include a "safety margin'' by setting , in case the highest-order shapelet coefficients are affected by PSF structure at even higher order. The scheme is illustrated for the typical case of N=8 in Fig. 1. The highest-order polar shapelet coefficients that should be included in the fit are , and .
In summary, from the shapelet series for each source, up to Cartesian order , we form a truncated polar shapelet series including terms of order , and . The effect of the shear, translation and PSF operators on circular basis functions up to order is then calculated up to order , and the result likewise converted to polar shapelets up to order .
Least-squares fitting the model to each source yields an ellipticity estimate (e_{1},e_{2}), expressed as the shear that needs to be applied to a circular source to fit the object optimally.
Performing the least-squares fit is straightforward to do numerically. is a fourth-order polynomial in the fit parameters , and the minimum can typically be found in a few Levenberg-Marquardt iterations (Press et al. 1986).
The errors on the shapelet coefficients for each source can be derived from the photon noise, and these can be propagated through in the function. The second partial derivatives of at the best fit give the inverse covariance matrix, which can be inverted to show the variance and covariances between the fit parameters. This results in proper error estimates on all parameters, in particular on e_{1} and e_{2}. In practice the errors on e_{i} are only weakly correlated with those on d_{i} and c_{n}.
We now describe some elementary tests of this approach.
For small ellipticity e the linear shear operator provides a good
description of the action of the shear, but for larger ehigher-order corrections come into play, and these corrections depend
on the radial brightness profile. We have calibrated these corrections
empirically using a set of model sources that follow a Sersic
(1968) distribution of brightness, of index 1-4. Each
source was sheared by varying amounts and encoded into shapelets using
a range of different scale radii. The first-order ellipticity estimate
derived by fitting a model of
Eq. (14) was then compared to the true ellipticity
(see Fig. 2). For small e the correct ellipticity
is recovered, but at larger e the discrepancy grows. Fitting the
residuals versus the radial profile shape parameters c_{0}, c_{2},
c_{4}, it turns out that the true ellipticity
can be
derived from the fitted 1st-order estimate by applying the correction
factor
Figure 2: Empirical calibration of the post-linear correction to the measurement of e. Top left: fractional error on derived 1st-order ellipticity e for a range of different scale parameters and Sersic indices. Top right: fractional error of the corrected ellipticity. Bottom left: the coverage of the c_{2}/c_{0}, c_{4}/c_{0}plane by the models. The four groups of points correspond, from top to bottom, to Sersic index 4, 3, 2 and 1. The horizontal spread is mostly a consequence of using different scale radii . Bottom Right: residuals on 1st-order ellipticity vs. the linear combination of c'/c_{0}=(0.085c_{2}+0.63c_{4})/c_{0} for input ellipticity 0.4, showing that this parameter drives the scatter. | |
Open with DEXTER |
Below we will, in fact, NOT apply this correction, because the errors on c_{0}, c_{2} and c_{4} are typically so high, and (in the case of small, barely resolved sources) correlated, that the correction factor cannot be determined accurately. Fortunately most galaxies in the sky have ellipticity below 0.3, where the correction is below 3%. If necessary, the accuracy could be further increased by evaluating the effect of shear on the shapelet basis functions to higher order in Eq. (16).
A truncated shapelet expansion can only describe deviations from a Gaussian over a particular range of spatial scales (which widens with order N). For ellipticity determinations the outer parts of galaxy and PSF images are most important (e.g., the classic second moments depend on the 3rd moment of the radial profile), so there is some advantage to taking as large a scale radius as possible. On the other hand, this radius should not be so large that the inner structure of the source cannot be resolved.
We have found that, for shapelet expansions up to order N=8 or higher, taking a scale radius which is 1.3 times the dispersion of the best-fit Gaussian works well for a range of model PSFs. Figure 3 shows an example for a Moffat PSF with index 2: both the core and the wings can be fitted adequately with this choice of .
Figure 3: Different shapelet fits to a Moffat PSF with FWHM of 12 pixels. Top panel: the solid line shows a cross-section of an 8th-order fit, using the dispersion of the best-fit Gaussian as scale radius. The dotted line shows the corresponding fit for a scale radius that is a factor of 1.3 times larger. The lower panel shows the PSF multiplied by radius^{3}, in order to accentuate the residuals in the wings. Dots are the actual (monte-carlo realized) pixelated PSFs used in the simulations of Sect. 4. Note how the increased scale radius makes for a much better fit in the outer regions, without a serious degradation in the core of the PSF. | |
Open with DEXTER |
The impact of the choice of scale radius on ellipticity measurements is shown in Fig. 4, for typical images. The factor 1.3 represents a good compromise, though its exact value is not critical.
Figure 4: Illustration of the effect of the choice of scale radius on the ellipticity measurement. Each curve shows, for a different galaxy profile, how the derived e depends on the choice of scale radius (expressed as a multiple of the dispersion of the best-fitting round Gaussian). Each model galaxy had an effective radius of 4 pixels and ellipticity 0.3, and was convolved with a PSF of Moffat index 2 and FWHM 8 pixels. The for convolved source and PSF are both scaled by the same factor. | |
Open with DEXTER |
For the tests we generate simulated, high signal-to-noise (S/N) sources. To be sure that pixelation effects are taken into account properly, and that the PSF convolutions are done accurately, a brute-force technique is used: each source is built of individual "photons'' that are drawn from a 2D Sersic distribution, sheared if required, then have a "PSF'' displacement added to them, and finally are added to the pixel in which they fall. We use 10 million photons per source, which gives effectively noise-free images (S/N 1000).
We ran three sets of tests. In each case we explore galaxies with Sersic laws with indices (Gaussian), 1 (exponential), 2, 3 and 4 (de Vaucouleur), and use PSFs with a Moffat profile of index , 3 and 9 (nearly Gaussian). All galaxies are scaled to an intrinsic effective radius of 4 pixels, and the PSF FWHM's range from 4 to 12 pixels. We use shapelet expansions to order N=8 and 12, and scale radii of 1.3 times the dispersion of the best-fit Gaussian.
First, to test the "shear calibration'' factor, we check how well we can recover the shear of a galaxy that is sheared by g_{1}=0.1 and convolved with a round PSF. Figure 5 shows the result: any calibration error is at the sub-percent level; the worst results are obtained for the most non-Gaussian PSFs ( ). The noise in the curves suggests that we are also limited by the accuracy of the Monte-Carlo simulations of the galaxies, and by the numerics of the software implementation of the method.
The second test shows to what extent PSF ellipticity pollutes the shear estimate. The input PSFs were given an elipticity e_{1}=0.1(axis ratio 0.82), and convolved with round galaxies. The recovered e is shown in Fig. 6. The residual effect is at most half a percent (worst case), which represents a correction of the PSF ellipticity by a factor of better than 20. The best results are obtained for lowish Sersic indices (below 2) and for PSFs with not too extended wings (Moffat index 3 or higher). Perhaps surprisingly, the larger the PSF (for the same galaxy size) the better the correction: presumably this is a small sampling effect.
Finally, we introduced a lopsided PSF (by giving 1/4 of the photons an extra offset of half the FWHM) and repeated the analysis. As can be seen from Fig. 7, this combination of dipole and quadrupole PSF distortion can also be handled.
In all cases, taking the expansion to N=12 increases accuracy, though not spectacularly. We conclude that the algorithm works: shapelets provide a promising technique for measuring galaxy ellipticities, and for correcting ellipticities for smearing by the PSF.
Figure 5: Fractional error in recovering a 10% shear, using round PSFs. Each panel represents a different Moffat PSF; the rightmost panels are very nearly Gaussian. The simulated galaxies have effective radii of 4 pixels. Top row: 8th-order shapelets; bottom row: 12th-order shapelets. | |
Open with DEXTER |
Figure 6: Residual shear after correction for an elliptical PSF (the same PSFs as Fig. 5, sheared by 10%). | |
Open with DEXTER |
Figure 7: Residual shear after correction for a lopsided PSF (the same PSFs as Fig. 5, but a third of their flux is displaced by (2/3) FWHM). | |
Open with DEXTER |
Figure 8: The result of Monte-Carlo simulations, in which many noise realizations of Sersic profile galaxies were run through the ellipticity-fitting procedure described in this Paper. Shapelet order N=8 was used throughout. Each plotted dot represents the average ellipticity of 2500 different noise realizations of the same galaxy image. The same data are plotted in both panels, but coded by different model parameters: the Sersic index on the left, and the PSF size on the right. The vertical axis shows the fractional scatter of the measured fluxes, , of the sources, as determined by integrating their shapelet series. No trend of the mean shear with S/N is seen: noise leads to scatter but no bias. | |
Open with DEXTER |
A final series of tests was used to check the behaviour of the algorithm on noisy images. We added Poisson noise to simulations such as those just described, and compared the result of many realizations. The PSF FWHM ranged from 4 to 12 pixels, the galaxies' effective radii were set to 4 pixels. The same ranges of Sersic (0.5-4) and Moffat (2, 3, and 9) indices were used as above, and the input shear was 0.1. Different levels of noise were added, roughly to span the S/N range between 10 and 100.
The results are summarized in Figs. 8 and 9, which show that the noise causes a scatter on the ellipticity estimates, but that this does not lead to a bias; and that the propagated error estimates on the ellipticities are a good measure for the rms scatter among the different realizations.
We have implemented the above ideas into a "shear pipeline''. It starts from a fully reduced image, detects sources, determines the PSF and its variation across the image, decomposes detected sources into shapelets, and obtains a shear estimate for each object. The pipeline consists of a number of stand-alone programmes that run in sequence. The modules operate on source catalogues and FITS images, and generate new catalogues and diagnostic plots. The pipeline can run fully automatically.
The first step in the reduction process is the detection of sources. For this the SExtractor software (Bertin & Arnouts 1996) is used. A few basic parameters are measured during extraction: position, the flux and its error, the FWHM^{}, major and minor axis length and position angle, and source quality flags. SExtractor is fast and effective, particularly for the relatively high S/N sources that are used here.
The PSF is determined from the stars in the image itself. We assume that the pixel values are linearly related to intensity.
First, the stellar locus is determined from the plot of magnitude vs. FWHM in the standard way (e.g., KSB). We fit a circular Gaussian to each star, and adopt the median dispersion from these fits, multiplied by 1.3 (see Sect. 4.2) as the scale radius .
For all objects near the stellar locus a shapelet expansion is fitted using as scale radius. The SExtractor centroid is chosen as initial centroid for the expansion, but after the expansion is completed the centroid is adjusted (using the 1st-order translation operators of R03) until the coefficients of B^{01} and B^{10} are exactly zero. Finding the required offsets involves solving a simple linear equation. Each shapelet is finally normalized to unit integral, by analytically integrating the counts in each shapelet term and dividing by the total.
Once this is completed, we have a shapelet description of (candidate) PSF objects scattered over the image. A map of the PSF variation across the image can now be made by interpolating the shapelets. We have found that a straightforward polynomial fit, coefficient by coefficient, works well, though complex PSF variations may require more sophisticated schemes such as weighted nearest neighbour averages (Christen 2006), Padé interpolants (Hoekstra 2004), or even physically-motivated model fits (Jarvis & Jain 2004; Jain, et al. 2005). During the fitting of the variation of each coefficient over the image, deviant points can be rejected, leaving a cleaned sample of PSF objects.
The result of this step is a recipe for the shapelet coefficients of the PSF at any point in the image.
Once the PSF is determined, all other detected sources are expressed as shapelets as well. As for the PSF objects, a shapelet expansion centered on the SExtractor coordinates, is fitted directly to the observed pixel values. The statistical errors on the pixel values are propagated through the least-squares fitting, leading to errors (and if desired, covariances) on the shapelet coefficients. In the case of well-resolved shapelets and uniform noise level across the source, the shapelet normalization is such that the rms error on each coefficient is the same, and the correlation between errors on different coefficients small.
Each source is encoded into shapelets with scale radius derived as described in Sect. 4.2. All pixel values within a radius of from the SExtractor centroid (at least 10 pixels) are used in the fit. For efficiency reasons in the shear estimation step, the allowed values are quantized: allowed values are , After fitting, the center of expansion for each object is shifted in the image plane by means of the shapelet translation operators until the 01 and 10 coefficients are zero, as before.
This procedure describes the source as seen in the image plane, i.e.,
after it has been convolved with the PSF and pixelated. An alternative
approach also make sense: to convolve all basis functions with the
pixelated PSF, and fit the observed source image as a combination of
those. This yields a shapelet description for the intrinsic,
pre-seeing, object shapes (Massey & Refregier 2005).
As long as the same procedure is followed for the sources and PSF, the
end result should be the same: the deconvolved image will be free of
pixelation and PSF. We prefer our approach because it leaves the
covariance between the fitted shapelet coefficients small.
All sources are encoded to the same shapelet order as the PSF
(typically 8 or 12), in order to avoid signal-to-noise dependent
smoothing effects. For faint sources, the higher-order coefficients
will therefore be very noisy, but still unbiased.
With a description of the shape of each source and the corresponding
PSF, the next task is to determine the intrinsic shape parameters that
are needed for a weak lensing analysis.
As explained in Sect. 3.3, we derive the ellipticity
as the shear that needs to be applied to a suitable round source in
order to fit the observed image. The fit is applied completely in
shapelet space. We use a scale radius equal to
for the intrinsic,
deconvolved, circular model galaxy, and normalize all sources to unit
flux before fitting so that es with
are
used. The C_{nml} convolutthe c_{n} can be used as radial profile
shape parameters. Only sourcion coefficients need to be evaluated only
once per value of .
The resulting catalogue of sources with shape estimators needs to be cleaned in order to remove sources which are affected by neighbours, edge effects, poor fits, etc. We apply various cuts:
Once individual estimates have been obtained for each source these need to be combined in some way to generate a shear estimate.
Conceptually the simplest methods are (i) to average the e_{i} and divide by (Eq. (9)), and (ii) to identify the mode of the ellipticity distribution (provided it is centrally peaked), which identifies the intrinsically round galaxies.
A better technique is to form a weighted mean, where the weight is driven by the amount of information about the shear field each source provides. The scatter in the ellipticity measurements of sources is due to the intrinsic dispersion in shapes , and to measurement errors. The latter can be estimated by propagating the noise in each image through the fitting procedure, and the former can be estimated as the excess variance of e in the source population. We therefore adopt individual weights when forming the mean of all measured ellipticities.
The same weighting is then applied to determine the value of
in Eq. (9). So the shear estimate
is determined as
(20) |
Depending on the form of the intrinsic shape distribution of galaxies, different weightings are optimal: for example, for a very peaked distribution of ellipticities higher weight can be given to nearly round sources, whereas for a top-hat distribution the sources with large ellipticity carry more information - for a discussion see BJ02. A problem with weight factors that depend on e is that the centroid needs to be found first as it is the intrinsic, pre-shear ellipticity that counts. The weight adopted above is optimal for a Gaussian distribution of intrinsic ellipticities.
Figure 9: Comparison between the scatter in ellipticity measurements from sets of 2500 random noise realizations, and the error predicted by propagating the pixel noise through the calculations. | |
Open with DEXTER |
In the tests below we will compare the weighting scheme just described, and a simple unweighted median. To the extent that the median identifies the center of the distribution of ellipticities of the source population, i.e., the intrinsically round sources, no correction needs to be applied to the median.
A rather realistic test of the whole procedure is provided by the "Shear Testing Programme'' (STEP, Heymans et al. 2006). Phase 1 of STEP has produced a large set of realistic simulated images across which a constant shear (in g_{1}) and PSF smearing has been applied. The PSFs mimic realistic optical aberrations. These images provide an important test, since unlike the tests presented in Sect. 4, the pre-seeing sources do not have a constant ellipticity with radius. Furthermore, the STEP images include photon noise, and the sources may overlap. The PSFs are also somewhat smaller (in pixels) than in the tests of Sect. 3.3: scale radii used range from 2.2 to 3.2.
Table 1: Summary of the results of the STEP1 simulations. For each PSF, the slope m_{1} and intercepts c_{i} of the best-fit linear relation between input and recovered shear are shown. The c_{i} have been multiplied by 100 for clarity. As no input g_{2} distortion was applied in the STEP1 simulations, m_{2} cannot be measured.
An analysis of the STEP1 data with an earlier version of this software was reported in Heymans et al. (2006). The main improvement in the implementation since then has been the use of polar shapelets in the shear determinations, which allows truncation effects to be curtailed properly, and the use of larger scale radii as discussed in Sect. 4.2.
Figure 10: Results from the STEP1 simulations, which are based on modeling of the optics of the CFHT. Each plotted point represents the average ellipticity of about 2200 sources in one STEP1 image. Shapelets to order were used. PSFs 0 to 5 are, respectively, round, with coma, with astigmatism, with defocus, and with 3rd and 4th-order astigmatism, and results for images with applied shears of 0, 0.05 and 0.1 are shown as three clusters of points in each panel. The correction factor is about 0.93 in all cases. Results from cluster to cluster are not statistically independent, but within each cluster of points they are. | |
Open with DEXTER |
Figure 11: As Fig. 10, but now the median ellipticity is used as shear estimator. The STEP1 galaxy ellipticity distribution is very peaked, which makes the median a very efficient estimator of the center of the distribution. | |
Open with DEXTER |
Results of the use of the present pipeline on the STEP1 data are shown in Fig. 10. For each of three g_{1} shear values, and six different PSFs, 64 separate simulated images were run through the pipeline, each image yielding shear estimates based on about 2200 galaxies. Table 1 summarizes the results per PSF in terms of a multiplicative correction factor m, and an additive offset c. It can be seen that in general the method suffers from very little bias. For the non-elliptical PSFs (0, 3, 4, 5), the recovery is perfect within the noise (m=1, c=0), which indicates that the correction for dilution by PSF smearing works. On the other hand, for the comatic (1) and elliptical PSF (2) there is a residual additive term, of nearly 0.003 and 0.005 in shear respectively. This result is consistent with that of the simulations in Sect. 4.3. In addition the elliptical PSF suffers from a multiplicative bias of nearly 4%. The origin of this discrepancy is not clear at the moment.
The STEP1 analysis revealed that many methods show a small systematic trend between the error in the derived shear and the magnitudes and sizes of the objects (Heymans et al. 2006). The technique presented in this Paper is no different, as illustrated in Fig. 12. These trends are still a puzzle, but it is clearly important to trace their origin and further improve the accuracy of the shear measurement. Possible causes include ellipticity-dependent incompleteness in the catalogues, problems estimating the intrinsic ellipticity dispersion , magnitude- or size-dependent neighbour contamination, or residual systematic issues in the method itself. Further work addressing these issues is on-going.
Figure 12: The recovered shear for the STEP1 simulations for PSF0 and input shears of 0, 0.05 and 0.1, split up by source brightness ( top) and size ( bottom). In both cases there is a systematic, so far unexplained trend. The upper panel in each plot shows the ellipticity dispersion correction factor derived for and applied to each bin. | |
Open with DEXTER |
The method described here has several advantages. It goes beyond the KSB description of PSF anisotropy as a convolution with a very compact PSF, and is in principle applicable to all PSF shapes. The correction for the PSF is performed in a single step, which avoids the need to separate the effect of the PSF into an anisotropic part that shifts the ellipticities, and a round part that causes a dilution of the source ellipticities. The forward-fitting approach of a PSF-convolved model for the intrinsic galaxy shapes to the observations allows error propagation. The fact that the ellipticities derived are "geometric'' in the sense of BJ02 (i.e., they represent the shear to apply to a round object in order to fit the source) means that there is no need to derive a higher-order "shear polarizability'', but instead the response of an ensemble of sources to shear can be predicted simply from the dispersion in ellipticities.
The use of shapelets in this method is not essential, but does help to speed up the calculations, and gives a natural framework for isolating the m=2 components that carry the ellipticity information about a galaxy. It also allows the source lists to be filtered efficiently based on objective shape criteria, and gives a robust way of interpolating the PSF shape across an image. In cases where the PSF cannot we well-described by a truncated shapelet expansion (for example, poorly-sampled space-based observations) it is possible to extend the technique, by performing the PSF convolution of Eq. (16) in pixel space instead of in shapelet space (Kuijken 1999).
As was shown by the STEP1 project (Heymans et al. 2006) a variety of techniques can be used to derive shears from ground-based image data, with residual errors around the 1% level. While the focus of STEP1 was on different variations of the KSB method, we have shown here that the shapelets technique can do as well. As KSB is expected to hit a fundamental level of systematic error (because of its good but imperfect description of the PSF effects), which may well be inadequate for the next generation of weak shear surveys, it is worthwhile to look to higher-order methods.
The approach we have taken here differs in several ways from that of Refregier & Bacon (2003) and Massey & Refregier (2005), which is also shapelets-based. We perform the shapelet expansions to a fixed order, and prefer not to introduce signal-to-noise dependent thresholds that may lead to biases in the derived shears (S/N dependent truncation and averaging do not commute). Our non-iterative procedure is also much faster. Our shapelet expansions describe the observed, post-seeing images, which means that the coefficients are statistically almost independent (whereas shapelet coefficients of the deconvolvedshapes are necessarily correlated), ensuring that they are virtually unaffected by truncation of the series. We explicitly allow for errors in the centroiding of the sources in our ellipticity estimates, and forward-fit the observed images to account for PSF effects rather than deconvolving them using a truncated PSF shapelet expansion. Finally, rather than modeling the response of a (model or empirical) population of images to a shear, we derive "geometric'' ellipticities (BJ02) for each individual source, and average these only at the end.
Compared to the Bernstein & Jarvis (BJ02) technique, the two main differences are the way ellipticity is measured, and the way PSF effects are handled. BJ02 derive geometric ellipticities by fitting an elliptical Gaussian form to the observed images. This is similar to, but not quite the same as, a low-order shapelet fit as used in this paper. For sources with constant-ellipticity isophotes the BJ02 method is exact, whereas ours is accurate only to o(e^{2}) (see Fig. 2). The benefit is a much faster numerical convergence. BJ02 correct the ellipticities for PSF convolution by means of a description of the PSF as a Gauss-Laguerre expansion (equivalent to polar shapelets). Optionally, the images can be convolved with a rounding kernel before the ellipticities and PSF are measured, in order to improve the accuracy of the PSF model. Thus the main difference with the approach here is that BJ02 first derive a post-seeing ellipticity, which is then corrected for the PSF; here we forward-fit the intrinsic ellipticity in one step. BJ02 separate the effects of PSF anisotropy (ellipticity bias) and circularly-symmetric smearing (ellipticity dilution). They correct for the latter by assuming an intrinsic light profile of the source that is a perturbation about a Gaussian, expanded up to the kurtosis. The fact that here we model the full intrinsic radial profile of the source as higher-order shapelets should provide higher accuracy.
The Kaiser (2000) technique is a more sophisticated kernel convolution technique, in which a convolution kernel is constructed which turns an image into one for which the effect of pre-seeing shear on all sources is known exactly. This allows one to find the shear that makes the source ellipticity distribution isotropic - this is then the opposite of the amount by which the population was sheared on its way to the telescope. The method is in principle exact, and operates in pixel space. It appears to have received relatively little use thus far (Wilson et al. 2001a,b; Dahle et al. 2002).
Shapelets provide a neat framework in which to describe the transformations that an astronomical source image undergoes until it is registered on a detector. Gravitational shear, convolution with a PSF, and pixelation can all be modeled within the shapelets formalism. All these elements can be combined into an efficient algorithm for extracting image ellipticities that can be used for accurate gravitational lensing shear measurements.
The implementation of these techniques into a working pipeline is presented in this paper. Tests show that the pipeline is able to recover input gravitational shears with very small calibration error (of the order of a percent) and PSF residual (better than a factor of 30 in PSF ellipticity).
It remains to be seen to what extent this approach can be applied successfully to diffraction-limited PSFs, which cannot be described easily with a shapelet expansion. A different set of basis functions for the expansion might be the answer. Further possible improvements are also under investigation.
Acknowledgements
Ludo van Waerbeke provided the simulated "STEP'' images used in Sect. 6. The Leids Kerkhoven Bosscha Fonds is thanked for travel support.
(A.1) |