Simulating image coaddition with the Nancy Grace Roman Space Telescope – II. Analysis of the simulated images and implications for weak lensing Open Access

ABSTRACT

1 INTRODUCTION

Weak gravitational lensing has long been recognized as a powerful probe to infer matter distribution and matter content in the Universe (e.g. Hoekstra, Yee & Gladders 2002; Albrecht et al. 2006; Weinberg et al. 2013). It has received great attention as a test of dark energy and modified gravity models, and is a leading tool in testing the consistency of Λ-CDM model between early-time (e.g. cosmic microwave background; Planck Collaboration 2020) and late-time (e.g. weak lensing) observations (Lemos et al. 2021). Further improvement when the next-generation observatories Vera C. Rubin Observatory (LSST Science Collaboration 2009; Ivezić et al. 2019), Euclid (Laureijs et al. 2011), and the Nancy Grace Roman Space Telescope (Spergel et al. 2015) start operating relies on our capability to accurately process raw images and analyse the processed images. Recent efforts to minimize the weak-lensing shear systematics were mainly proposed and applied in large optical imaging survey collaborations such as the Dark Energy Survey¹, Hyper-Suprime Cam² and Kilo-Degree Survey.³ These include a new image coaddition technique (e.g. Becker et al. in prep.) to avoid PSF discontinuity, accurate PSF modelling (e.g. Jarvis et al. 2021), and shear calibration and characterization (e.g. Sheldon et al. 2020; Li & Mandelbaum 2023). These are mainly developed to mitigate various biases related to cosmic shear measurement using ground-based telescopes. Images from space-based telescopes, however, are often undersampled at the native plate scale, which introduces additional complications into precision measurements such as galaxy shapes.

Since it is not possible to recover the original astronomical scene from undersampled images through operations such as deconvolution, object shapes measured in undersampled images can suffer biases (see e.g. Finner et al. 2023 for a study with the Hubble Space Telescope). A study of the fundamental mathematical issues can be found in appendix C of our companion paper (Hirata et al. in prep.; hereafter ‘Paper I’). Kannawadi, Rosenberg & Hoekstra (2021) and Yamamoto et al. (2023) measured galaxy shapes with the Metacalibration technique (Huff & Mandelbaum 2017; Sheldon & Huff 2017) on simulated Euclid or Roman images and presented a few per cent shear calibration bias. Even though these studies did not particularly identify the sources of other biases except aliasing bias, in order to meet the systematic error budget on shear calibration defined by the Science Requirements Document (SRD).⁴ We need to mitigate this aliasing bias at the image level by reconstructing well-sampled images with survey-specific dithering strategy (e.g. telescope rotation and pointing) and realistic simulated PSFs.

As discussed in Mandelbaum et al. (2023), in order for a reconstructed image to have a well-defined PSF, the image coaddition algorithm needs to be linear and the weights on each exposure need to be independent of signals of sky scene (e.g. inverse-variance weights, including source Poisson noise break the linearity assumption). To produce a well-defined coadded PSF, Rowe, Hirata & Rhodes (2011) presented the linear image combination algorithm imcom. imcom finds an optimal matrix T mapping from input (native) to output (coadd) pixels by minimizing the cost function that consists of the ‘leakage’ (squared L² norm of the difference between the output PSF compared to a user-specified ‘target’ PSF) and the output noise variance. The relative weight of the leakage and noise in the objective is controlled by a Lagrange multiplier κ in Imcom; the logic in the iterative solver for κ can be configured for different outcomes (e.g. to minimize noise subject to leakage being less than a specified value).

Imcom is fundamentally different from most other image reconstruction techniques in the sense that users are able to design their desired PSF for a given area and noise correlation level. Qualitatively, it is a process that transforms the input images with PSFs into a combined image with a common PSF in every pixel for a given area. If one would like a reconstructed image to be well-sampled and to have an isotropic and homogeneous PSF (indeed the case that is explored in Paper I to avoid PSFs with diffraction spikes), the size of the output PSF needs to be larger than the original PSF; hence the reconstructed image would be lower resolution than a ‘typical’ coadd image. Despite this trade-off, the size of the output PSF is still much smaller than that of a ground-based PSF (see table 4 of Paper I), meaning that this method still takes advantage of space-based imaging. Readers may refer to Rowe et al. (2011) and section 2 of Paper I for the full details of the algorithm and the motivation in the current context of Roman, respectively.

With the recent development of realistic image simulations (Troxel et al. 2021, 2023) for Roman following the specifications of the telescope, Paper I re-implemented Imcom with a python interface and c back-end so that the pipeline is compatible with the image simulations. Paper I describes the implementation of Imcom and updates from the original version (e.g. how to compute input PSF correlations. see section 4 of Paper I). It then applied Imcom on various input layers of 48 × 48 arcmin footprint (including realistic Roman single-exposure images produced in Troxel et al. (2023); the description of different input layers is in section 3 and table 1 of Paper I). The target PSF specified in Paper I is the Airy disc convolved with a Gaussian kernel whose width depends on the bandpass (Paper I table 4), and Imcom is configured to produce (if possible) an output PSF within 0.1  per cent of the target PSF in an L² norm sense. This implies that the root-sum-square of the error in the moments is 0.1  per cent if the moments are defined in an orthonormal basis such as shapelets (Refregier 2003). Paper I demonstrates several aspects of Imcom on Roman simulations:

• Reconstruction of Nyquist-sampled images;

• Homogenization and isotropization of well-defined coadd PSF; and

• Minimization of noise covariances in coadded images.

This paper goes further and examines the moments of simulated sources and the correlation functions and power spectra of the output noise and residual systematics in the coadded images. These studies are needed to support error budgeting for the Roman weak lensing analysis. Most of the analyses are carried out on stars, since this is expected to be a ‘stress test’ maximizing the impact of undersampling.

This paper is organized as follows. In Section 2, we briefly review the input data. Section 3 computes the noise power spectra of the output images for both uncorrelated (white) and correlated (1/f) input noise. A series of the quantitative analyses of the coadded simulated input images is presented in Section 4: we

1. Use injected stars to test the output normalization, astrometry, size, and shape of the final coadd PSF after propagation through Imcom (Section 4.1);

2. Measure the properties (i.e. shape and size) of the stars based on the truth positions of the input stars (Section 4.2);

3. Measure correlation functions of the ellipticities of both injected and simulated stars at scales overlapping the range likely to be included in the Roman weak lensing analysis (Section 4.3);

4. Measure correlations of the fourth moments of the injected stars (Section 4.4); and

5. Estimate the noise-induced additive bias on star ellipticities using the injected stars and the Monte Carlo noise realizations, and compare this to analytic expectations (Section 4.5).

Section 5 presents a synthetic wide band image (averaging all 4 bands, smoothed to a common PSF) and the ellipticity correlations of simulated stars. We consider directions for future work in Section 6. The appendices cover analytic models for the 2D noise power spectra (Appendix  A) and the analytic theory of additive noise bias (Appendix  B).

2 SIMULATED DATA

This study utilizes the simulated data products from Paper I, which covers 0.64 square degrees in four filters: Y106, J129, H158, and F184 from shortest to longest wavelength. The sky inputs are from the Large Synoptic Survey Telescope Dark Energy Science Collaboration Data Challenge 2 (LSST DESC DC2; Korytov et al. 2019; LSST Dark Energy Science Collaboration (LSST DESC) et al. 2021a, b; Kovacs et al. 2022). Table 1 of Paper I illustrates various input images that have been fed into the Imcom algorithm, which can be classified as follows.

• Sky image (layer names: SCI, truth): The images were generated through the Roman image simulations pipeline (Troxel et al. 2021, 2023) utilizing the truth catalogs of stars and galaxies from the LSST DESC CosmoDC2 simulations (Korytov et al. 2019; Kovacs et al. 2022). The layer labelled SCI includes saturation and simple detector physics models such as read noise and dark current, simulated in GalSim, while the layer labelled truth refers to noiseless images; star and galaxy profiles are drawn at designated locations.

• Point sources (layer names: gsstar14, cstar14): These are grids of sources located at healpix resolution 14 (nside = 2¹⁴) pixels (Górski et al. 2005). Each source is a δ-function convolved with a simulated Roman PSF; the PSFs are simulated in Troxel et al. (2023) through GalSim using flat spectral energy distribution (SED). gsstar14 refers to images where the sources were drawn using GalSim, while cstar14 refers to images where they were drawn with internal routines in the coaddition code. We call them ‘injected stars’ throughout this paper.

• Noise images (layer names: whitenoise1, 1fnoise2): We generate realizations of white noise (whitenoise1) and 1/f noise (1fnoise2). The white noise is uncorrelated between pixels and has unit variance. The 1/f noise is scale-invariant in the time stream and is then re-formatted into a 2D image according to the pixel readout order (Mosby et al. 2020). It has unit variance per logarithmic range in frequency.

Examples of the injected stars and noise realizations are shown in Fig. 1. And examples of a variety of objects identified in the coadded sky images are shown in fig. 8 of Paper I.

Figure 1.

Coadded injected GalSim stars (left), white noise (centre), and 1/f noise (right) realizations, displayed as three-colour F184 (red)/J129 (green)/Y106 (blue) combinations. Each image shows a 700 × 700 output pixel (17.5 × 17.5 arcsec) region of the coadded images, from the gsstar14, whitenoise1, and 1fnoise2 layers, respectively. The colour scale is a fourth-root stretch (0–0.2 input flux per input pixel) in the injected star image (left); and for the noise realizations it is linear, spanning ±1.25 (centre) or ±4 (right) in input units. Note that the input white noise layer leads to output noise correlated on the scale of the input pixels, whereas the 1/f noise layer shows the characteristic striping at each of the two input rolls. These rolls are at different angles in each filter, hence the colour pattern. The region shown is 270 ≤ x < 970, 301 ≤ y < 1001 of block (2,30).

Open in new tabDownload slide

Throughout this paper, we use the quantity called ‘fidelity’ as the basic quality indicator of the Imcom output, and many statistics are reported as a function of the fidelity. This defines how well the output PSF recovered the target PSF we specified. Mathematically, fidelity depends on the difference between the output and target PSF,

where the index α indicates an output pixel; i indicates an input pixel; G_i is the PSF in the ith input pixel; and |${\boldsymbol R}_\alpha$| and |${\boldsymbol r}_i$| are the output and input pixel positions. The fidelity is defined as the square norm of the output PSF residual scaled to the square norm of the target PSF and re-written as an inverse logarithmic measure:

The fidelity is usually – but not always–better (larger) if there are more exposures. The fidelity map in our simulation footprint for each bandpass we simulated can be found in fig. 10 of Paper I.

3 MEASUREMENT OF NOISE CORRELATIONS

We first investigate the noise fields, giving an overview of noise effects (Section 3.1), and then describing the 2D (Section 3.2) and 1D azimuthally averaged (Section 3.3) noise power spectra.

3.1 Overview

Noise is a significant source of error for the precision necessary to observe weak gravitational lensing of galaxies. The impact of noise on the measurement of a shear signal can be parametrized into two types of systematic observational errors: additive and multiplicative biases (e.g. Heymans et al. 2006). If the true signal is given by γ_true, noise (or other) biases give a γ_meas of

The additive bias (also called ‘spurious shear’) is given by the constant c, and manifests as a shear signal that is present even when the true population of galaxies is unlensed. The factor of m gives the multiplicative bias (also called ‘calibration error’). Multiplicative bias occurs when a real lensing signal is detected, but the measurement is larger or smaller than the true signal by some factor (Bacon et al. 2001; Erben et al. 2001; Hirata & Seljak 2003). The next-generation weak lensing surveys have stringent requirements for both types of errors (Paulin-Henriksson et al. 2008; Cropper et al. 2013; Massey et al. 2013; The LSST Dark Energy Science Collaboration 2018; Euclid Collaboration 2020; Troxel et al. 2021).

To quantify the impact that noise might have on observations and determine whether these requirements will be met, we measure and analyse noise correlations in the coadded images from two types of input noise fields: white noise and 1/f noise (see section 3.4 of Paper I). Real noise has both of these components, as well as some smaller additional terms (Rauscher 2015). We focus on the noise power spectrum in the output images, since both the lowest-order noise-induced bias and the variance of the measured shapes are proportional to second moments of the noise [i.e. depend on (S/N)⁻²; Bernstein & Jarvis 2002] and are therefore captured by the power spectrum. We then compare this power spectrum to the expected output noise power spectrum if we ignored sampling issues. Appendix  B describes in detail the impact of anisotropic noise (such as 1/f noise and non-ideal white noise) on ellipticity measurements.

Different properties of a galaxy are affected by noise at different wave numbers. If a large scale uniform (zero wavenumber) noise offset is present in an image, the error in this flux would impact the estimated photometric flux from the galaxy, but not the position or shape. An astrometric measurement will, however, be affected by an overall tilt (gradient of the noise, or in Fourier space the noise weighted by one factor of wavenumber); this would bias the centroid towards one direction or another. The shape of a galaxy will only be biased if the second derivative of the noise field is biased towards a preferred direction. In Fourier space, the ellipticity error thus comes from the noise weighted by two powers of the wave number (the impact of noise on shapes is discussed more in Section 4.5). Photometry, astrometry, and shape measurements are thus dependent on higher and higher spatial frequencies respective to one another. This is demonstrated quantitatively in Appendix  B, including the factor of wave vector squared (u² − v² for g₁ and 2uv for g₂) appearing in the noise-induced error in the shapes.

3.2 2D noise power spectra

We take the convention that the power spectrum P(u, v) of a field is given by the 2D Fourier Transform of the correlation function ξ(Δx, Δy):

where ξ(Δx, Δy) = 〈S(x, y)S(x′, y′)〉 is the correlation function of the noise field S with pixel coordinates given by (x, y) in Cartesian space and (u, v) in Fourier space. In the case of discrete data, we replace the integral with a sum over pixels; we preserve the convention that P has units of |$S^2\times \,$|area (e.g. for a field sampled at the output pixel scale s_out as considered here, we sum over pixels and multiply by |$s_{\rm out}^2$|⁠).

We accomplish this by making use of the Numpy Fast Fourier Transform (FFT) function (Harris et al. 2020). Each noise field in each N × N pixel block (here N = 2600) is FFTed and normalized by the size of the blocks and the size of the output pixels:

where u and v are sampled at integer multiples of 1/(Ns_out). The 2D power spectra that result from this are further binned into 8 × 8 bins, resulting in power spectra that are 325 × 325 pixels covering the range of u and v from −20 to +20 cycles arcsec⁻¹ (the Nyquist limit for s_out = 0.025 arcsec pixel⁻¹) at spacing

For each observation filter, we then averaged together the power spectra from each of the 2304 blocks (the total 48 × 48 arcmin² region). The resulting average power spectrum for the full mosaic in each observing filter can be seen in Fig. 2. Note that this procedure included the padding regions that appear in more than one block, so these regions are overrepresented in the averaged power spectrum. Since the padding regions are not special relative to the detectors or tiling strategy, we do not think this is a major issue.

Figure 2.

2D averaged power spectrum of the white noise field of each band, plotted on a logarithmic colour scale. The horizontal and vertical axes show wave vector components (u and v respectively) ranging from −20 to +20 cycles arcsec⁻¹. The colour scale shows the power P(u, v) in units of arcsec² [equation (4)]. The minimum and maximum values of each power spectrum are as follows: 6.6 × 10⁻⁸ to 6.7 × 10⁻³ for Y106; 3.4 × 10⁻⁸ to 3.5 × 10⁻³ for J129; 5.9 × 10⁻⁸ to 3.2 × 10⁻³ for H158; and 2.0 × 10⁻⁸ to 4.8 × 10⁻³ for F184.

Open in new tabDownload slide

In the absence of sampling issues and if the input and output PSFs were all identical, input white noise should lead to output white noise and a constant power spectrum for the 2D image. Since we have an output PSF that is larger than the input PSF in real space (narrower in Fourier space), we expect that the high-wave number Fourier modes will have reduced power (proportional to the square of the Fourier transform of the output PSF). The examples in Fig. 2 clearly show this behaviour. Note that the spectra are plotted on a log scale, so that although several features are visible, any features outside of the central maxima have very little power. The circular region in the centre of the 2D power spectra comes from the output PSF: the Airy disc convolved with a Gaussian. In the Y106 and J129 bands, where the Gaussian dominates the cutoff of the PSF, we see a smooth fade into uncorrelated wavenumber space. In the H158 band, especially the F184 band, the Airy disc defines the limit of the PSF, so beyond this limit the image contains only noise – this causes Imcom to set the background to zero, causing the hard edge of the circle.

We can confirm this by comparing measurements of the image with expectations. The band limit (maximum number of cycles per unit angle) for a diffraction-limited telescope is given by D/λ, where D is the telescope diameter and λ the wavelength of light. For the F184 band,

This is indeed the radius of the circle on the F184 band image, confirming that this feature is caused by the Fourier modes outside the circular regions being beyond the target PSF band limits.

The large ‘+ signs’ extended throughout the images correspond to the directions of the postage stamp boundaries. (Note that the input images are at various roll angles, so these features must instead be associated with the output.) The step function-like feature of the postage stamp boundaries becomes a continuous line in Fourier space, so we see these extended + signs in the power spectra.

The spots in the F184 band are a small print-through of the initial pixel positions in the input images into the final output images. The first ring of eight points is located at a radius of ∼9 cycles arcsec⁻¹ from the centre, i.e. roughly the inverse of the input pixel scale |$s_{\rm in}=0.11\,$|arcsec. As expected, the directions of the eight points corresponds to the four-grid directions of the input exposures.

For the 1/f noise power spectra in Fig. 3, we have combined two sets of images with horizontal banding at different roll angles, so the output noise image will have at least two distinct preferred directions (see the right-hand panel of Fig. 1). These features appear strongly in each panel of Fig. 3 as a distinct X through the centres of the images. The roll angles are slightly different in each band in order to maximize coverage, leading to different orientations of the X’s in each image.

Figure 3.

The 2D averaged power spectrum of the coadded 1/f noise field of each band, plotted on a logarithmic colour scale. The horizontal and vertical axes show wave vector components (u and v, respectively) ranging from −20 to +20 cycles arcsec⁻¹. The colour scale shows the power P(u, v) in units of arcsec² [equation (4)]. The minimum and maximum values of each power spectrum are as follows: 1.8 × 10⁻⁶ to 5.4 × 10¹ for Y106; 8.5 × 10⁻⁷ to 4.6 × 10¹ for J129; 6.9 × 10⁻⁷ to 4.8 × 10¹ for H158; and 2.1 × 10⁻⁷ to 5.6 × 10¹ for F184. The X-shape, + sign, spots (in H158 and F184), and vertical fringes (in J129) are discussed in the main text.

Open in new tabDownload slide

As in the white noise spectra, the boundary between postage stamps contained in the images creates step-like features in the noise fields, which manifest as two distinct perpendicular Fourier modes in the power spectra. However, this feature is clearly more prominent in the 1/f noise than in the white noise. In 1/f noise, we have striping across entire channels within the detector, which then crosses over postage stamp boundaries. The correlations between postage stamp boundary-caused features extend over much larger scales in real space, corresponding to sharper features in Fourier space. Thus we see the postage stamp boundary feature (the large + sign) more distinctly in these images.

The 1/f bands also display the same reduction of power at large wave number as in the white noise, and the sharp circular boundary in H158 and especially F184 caused by the band limit of the PSFs.

In addition to features coinciding with behaviours seen in the white noise images, we see one additional pattern appear in the 1/f noise that is not present in white noise data. Fig. 3 shows this feature most clearly in the J129 band: alternating bright and dark vertical fringes across the centre of the 2D power spectrum image. We measure the spacing between fringes to be 0.8 cycles arcsec⁻¹, or one cycle per 1.25 arcsec postage stamp. We therefore believe that the fringes result from steps at the postage stamp edges, with a sense that is coherent across multiple postage stamps (due to the large correlation length of the 1/f noise).

While these rough analytical arguments allow us to identify specific features in the output 2D power spectra with properties of the input images and Imcom algorithm, the quantitative details – how much power appears in each feature, and how this varies as a function of sampling going from Y106 (bluest/worst sampled) to F184 (reddest/best sampled) – must be determined via simulations.

3.3 Azimuthally averaged power spectra

In addition to the 2D power spectra, we generate for each filter a set of 1D azimuthally averaged power spectra. First, the 2D power spectra of each block are saved from the previous step of analysis. We then calculate the mean exposure coverage in each block, and group the power spectra into one of five bins in mean coverage. Since the root-mean-square of noise in the output noise images shows the mean exposure coverage dependence, it is good to check how noise correlations depend on the coverage as well (see section 5.3 of Paper I for more details). These power spectra are also separated by observing filter, as exposure coverage varies significantly depending on the band in which observations are taken. We average together the 2D power spectra for blocks within a given mean coverage bin into a single 2D power spectrum, which is then azimuthally averaged over thin annuli, using the method from Casey et al. (2023). For this work we used 162 radial bins, effectively taking rings of width one pixel in the power spectrum map [equation (6)] and averaging the values of the power at each circular aperture of wavenumbers.

Fig. 4 shows the mean coverage-binned 1D power spectra for the output white and 1/f noise images in each observation filter. Rather than just P(υ), we plot 2πυP(υ) so that the area under the curve corresponds to the variance in the output pixels. We additionally include in each figure a purely analytical model of the 1D power spectrum, based on N_in = 5 input exposures and ignoring pixelization/sampling issues. Derivations of the analytical models for the output noise power spectra can be found in Appendix  A.

Figure 4.

Top row: 1D power spectra for the output white noise fields in each filter. Bottom row: 1D power spectra for the output 1/f noise fields in each filter. Each filter’s spectra are divided into five even-width bins of mean coverage (‘mc’ in short), and plotted against the analytical expectation for noise power spectra for combining 5 exposures in the absence of sampling issues (see Appendix  A for derivations).

Open in new tabDownload slide

For all bands and both types of noise, the peak power and width of the power spectrum depends consistently on the mean coverage in the image. Exposures with small mean coverage values have the highest peaks and the slowest dropoffs down to zero power at large wave number. Our analytic expectation represents an idealized case, and the behaviour as compared with the spectra for the mean coverage bins shows that with decreasing mean coverage, we move farther above the ideal power spectrum. However, we note that in all of these cases the noise power spectrum is above the analytic expectation for most wave numbers.

In the white noise 1D power spectra, several features are evident in all bands. Each of the observation bands shows that the contributions to the power spectrum (weighted by 2πυ) go up, rise to some peak, and come back down to zero, qualitatively consistent with the analytic expectation. Noise increase due to aliasing often gets worse as one moves to larger wave number (the zero-mode corresponds to total flux, which is conserved in the absence of intrapixel sensitivity variation). This can result in 2πυP(υ) having a steeper dependence than ∝ υ at small υ, as seen prominently for Y106 (the most undersampled filter) in Fig. 4. Fig. 5 shows the central feature in the 2D power spectrum from input white noise for each filter with the colour scale significantly stretched to show this behaviour.

Figure 5.

Zoomed in image of the central features in the input white noise power spectra. By stretching the colour scale we can see more clearly that the zero-wavenumber modes do not contribute the most to power in the output noise, particularly in Y106.

Open in new tabDownload slide

In contrast, the 1/f noise spectra have a distinct peak in the centre of the image, caused by the overlaps in the Fourier modes from the roll angles. The power spectrum as a function of frequency on one dimension goes like 1/f, so for small wave numbers the power should reach very large values – a result that is clearly visible in the spectra for all bands. The analytic estimation for the 1/f noise spectra in the absence of pixelization/sampling effects (Appendix  A) is qualitatively correct for the redder filters, but the normalization is slightly low, and it misses the ‘bump’ caused by aliasing in the bluer filters (especially Y106).

Analysing these power spectra allows us to understand the correlations between features caused by noise in the simulated images being operated on by Imcom. In Section 4.5, we provide a more detailed discussion of the quantitative links between these power spectra and weak lensing measurement biases.

4 MOMENTS ANALYSIS OF THE COADDED IMAGES

In this section, we present the quantitative analysis of the simulated coadded images of multiple layers. This includes object centroid, ellipticity, and size measured by computing image moments of objects located in grids and simulated fields with and without noise.

4.1 Moments of the injected sources

Our first set of tests is conducted on the grids of simulated stars. Our main objective is to show the properties of the output PSFs and demonstrate that these properties have an expected dependency on how divergent the output PSF is from the target. The expected properties are exactly known since we chose the target PSF to be the Airy disc convolved with the Gaussian kernel.

Our simulation footprint contains 54 597 unique simulated stars, each at the centre of a healpix⁵ (Górski et al. 2005) resolution 14 pixel. For each injected star, we cut out a 99 × 99 output pixel (2.475 × 2.475 arcsec) postage stamp from the output block containing that star (for the ‘block’ definition refer to fig. 4 of Paper I). We consider the first the GalSim noiseless star layer (gstar14).⁶ We measure the first and second-order moments of the simulated sources in the output images using the galsim.hsm module (Hirata & Seljak 2003; Mandelbaum et al. 2005), which implements the fitting of an adaptive elliptical Gaussian to the image (e.g. Bernstein & Jarvis 2002).⁷ The first moments (centroids) can be compared with the expected position based on the World Coordinate System of the coadded image; the result is reported as the astrometric error. The second moments (covariance of the Gaussian: a 2 × 2 matrix M) are reported as three real numbers: the shear-invariant width

and the shape components

that are zero for a perfectly circular object and generally satisfy |$g_1^2+g_2^2\lt 1$| for a positive definite second moment matrix. (The use of ‘g’ for the latter indicates consistency with the conventions of Schneider & Seitz 1995.) We also report the mean fidelity 〈 − log₁₀(U_α/C)〉 for the pixels in the central 20 × 20 output pixel (i.e. 0.5 × 0.5 arcsec) region surrounding the star. This is a measure of how well Imcom estimates it has done at building an output image with the desired PSF in that region of the survey.

Fig. 6 shows measured ellipticity, astrometric error, and fractional size difference of the sources drawn by GalSim. It is expected that the anisotropy of shapes decreases as the fidelity of the output image is maximized. Most of the postage stamps (for the ‘stamp’ definition refer to fig. 4 of Paper I) achieved high fidelity (above 50) in all bandpasses except Y106. This is likely due to the bandpasses being the most undersampled and that the algorithm is unable to find the transformation matrix that the leakage and noise correlations are minimal (this is also shown in fig. 11 of Paper I). In all bandpasses, however, the total RMS ellipticity of the coadded injected stars in each shape component is ≲5.7 × 10⁻⁴ which is the PSF ellipticity requirement on angular multipole scales (32 < ℓ < 3200) determined by the SRD. Here, we also demonstrate that the astrometric error of the sources induced by Imcom is a small contribution to the relative error budget in the astrometric calibration (<1.3 mas) defined by the SRD. Finally, the bottom panel of Fig. 6 shows the fractional size error of the injected stars relative to the target PSF. We cannot directly compare the RMS size error in each bandpass to the required PSF size error in the SRD (≲3.6 × 10⁻⁴)⁸ because our RMS errors are computed at individual points, whereas the size error requirement is computed at scales ℓ < 3200 (or smoothed at a scale of ∼4 arcmin, i.e. where there are 3200² pixels on the full celestial sphere). We see that the size errors are largest in Y106 and F184. If we bin the stars into pixels, we see that the size error declines to (1.58, 0.59, 1.01, and 2.09) × 10⁻⁴ in Y106/J129/H158/F184 respectively in 1 × 1 arcmin pixels; and then it declines to (0.46, 0.39, 0.70, and 1.46) × 10⁻⁴ in 4 × 4 arcmin pixels. Thus we see that the image combination step is using only a small fraction |$(1.46/3.6)^2 = 16~{{\ \rm per\ cent}}$| of the overall PSF size error budget in a root-sum-square sense.

Figure 6.

The distributions (in an arcsinh scale) of moments of the injected stars drawn in the input images, coadded in Imcom, and then measured by the galsim.hsm module in each of the filters, with no noise. Top: the ellipticity |$g=\sqrt{g_1^2+g_2^2}$| of the injected GalSim stars. Middle: the astrometric displacement from the coordinates where the star was injected. Bottom: the relative size error of the GalSim stars and target PSF. The definition of size in this figure is equation (8). The horizontal axis is the output PSF fidelity (larger means a better match to the target PSF; see Section 2). The solid line shows the RMS ellipticity (top), RMS astrometric displacement (middle), and RMS fractional size error (bottom) in each fidelity bin. The dashed line is the same thing but is cumulative (i.e. the RMS for all stars in regions with that fidelity or better). This would be applicable to cases where we impose a mask based on the fidelity.

Open in new tabDownload slide

The spatial distribution of the second moment errors, averaged in 1 arcmin² pixels, is shown in Fig. 7. Most of the area has ellipticity errors at the ≲ 10⁻⁴ level; one sees a few outliers, with corresponding errors in the size. The large outliers seen in Y106 and F184 correspond to low-coverage regions (see fig. 1 of Paper I). What matters most for weak lensing purposes are the two-point statistical features of these maps; the correlation function will be investigated in Section 4.3.

Figure 7.

A spatial map of the size error (colour scale) shape (whiskers) of the injected stars. Each panel shows the 48 × 48 arcmin footprint of the simulation in one of the four filters. Stars are averaged into pixels (1 arcmin² each). The colour scale shows size error (σ_star/σ_target − 1) on an arcsinh scale, and the whiskers are scaled so that a length of 1 pixel corresponds to g = 5 × 10⁻⁴. This is presented on a square grid, however we have checked that the main features are also present if gridded in healpix and so they are not gridding artefacts.

Open in new tabDownload slide

We note that these quantities are also measured on the sources drawn by our internal ‘croutines’, and we find consistent results as GalSim sources.

4.2 Measurement of stars in the simulated images

For precise measurements of galaxy properties, of particular interest is the propagation of the PSF through coaddition algorithms. In principle, Imcom tries to remove any anisotropy or inhomogeneity from the coadded PSF and coadd stars should result in being round. However, since the leakage (square norm of the difference of output and target PSF) is not exactly zero there is some residual ellipticity of the coadded stars even in the absence of noise. Although stars in the coadded images are not utilized directly for modelling the PSF (with Imcom, the PSF determination step must be performed on single-epoch images prior to coaddition⁹), the coadded stars are an important test of PSF propagation and any ellipticity that we see at the end will have implications for galaxy measurement.

Stars are identified in the coadded sky images based on the truth locations of these stars. We use the simple simulated sky images for Roman simulated in Troxel et al. (2023), which utilized the LSST DESC DC2 simulations (LSST DESC et al. 2021b). Its stellar catalog was based on Galfast (Jurić et al. 2008; Jurić 2018) simulations. We cut out postage stamps around the truth locations and measure the star properties using the galsim.hsm module (which allows the flux, centroid, and second moments to float). The truth catalog contains locations, magnitudes estimated in Roman bandpasses from object flux based on objects’ SEDs, and whether the star is a candidate for PSF modelling. The criteria for being a candidate star defined in Troxel et al. (2023) are that:

• No pixel in the simulated star stamp is saturated in at least one exposure of the star (saturation was chosen at 10⁵ e in the images with simple detector models); and

• The star has a detection signal-to-noise (S/N)_det above 50, as defined by (S/N)_det  = |$0.015 \times (\rm total\ flux)$|⁠, where 0.015 was the typical background inverse noise level. (The observed signal-to-noise ratio including source Poisson noise is slightly lower.)

The S/N level was chosen here to be permissive to allow more restrictive selections later. We used this criteria for selecting clean stellar samples to conduct our measurement. We additionally select stars whose magnitudes are fainter than 18 in all bandpasses since we found that some single-exposure images contained saturated stars. (In the real survey, some information might be obtained from saturated stars because the detectors are read non-destructively and the early, non-saturated reads can be used; however, it is possible that the PSF would be different since one does not average over telescope motion during the full exposure, and in any case the present simulations do not include the early reads.)

With these selections, we consider 2623 (Y106), 2674 (J129), 2588 (H158), and 2647 (F184) stars in our simulation footprint. For each star in each of the Imcom images (SCI and truth), we cut out 99 × 99 output pixels postage stamps resulting in a stamp 2.475 arcsec on each side, whereas we extract 43 × 43 pixels around the truth centroid of a star in Drizzle images resulting in a stamp 2.473 arcsec on each side. Although the cutouts from Drizzle coadded ‘SCI’ images (Drizzle: SCI) are directly comparable with those of Imcom coadded ‘SCI’ images (Imcom: SCI), the measurements of stars on ‘truth’ images are not. While we cut out postage stamps from Imcom coadded ‘truth’ images (Imcom: truth), Drizzle coadded images of the ‘truth’ layer were not simulated by Troxel et al. (2023) at the time of this project. Examples of Imcom and Drizzle images are shown in Fig. 8.

Figure 8.

Coadded images of SCI and truth layers for a selected region for Y106, J129, H158 and F184 bandpasses. Left: Drizzle coadded SCI image, Middle: Imcom coadded SCI image. Right: Imcom coadded truth image. The flux in Drizzle images are scaled by (0.11/0.0575)² to account for the difference in the output pixel size and normalized to show the same colour range. The area that is shown here is 32.7 × 19.2 arcsec centred around RA|$=53.006\,$|deg and Dec=−40.027. Note that this is near the centre of our simulated output region.

Open in new tabDownload slide

Compared to the measurement of injected stars (Section 4.1), the measurement on the cutout of stars is expected to contain complications that could directly affect the shape and size of an object. These complications include background and simulated-detector noise, blending, and chromaticity (since the injected sources are drawn with the same SED passed to Imcom, but the simulated stars are drawn with their ‘true’ SEDs). The simulated stars therefore test more of the sources of systematic biases in weak lensing surveys. Some of these effects can be seen in Fig. 9 as outliers from the stellar locus in each bandpass. Fig. 9 additionally displays stellar locus of the same stars in Drizzle images, which is located at a smaller size than Imcom stars. This is expected because Drizzle does not smear the PSF. However, there is a large dispersion in sizes especially in the bluer, more undersampled bands. The Drizzle algorithm is designed to preserve total flux in its ‘raining down’ procedure (Fruchter & Hook 2002), but the spatial spread in the output image will depend on the specific locations of the pixels. Even though the output PSF from Imcom is larger than Drizzle (especially in Y106), it is designed to be uniform: it must be large enough so that that output resolution could be achieved even in a part of the image where the pixels interlace each other differently (see Paper I, fig. 7). The Imcom output PSF is still much smaller than ground-based PSF and the narrowness of the locus will help with star-galaxy separation in the real mission.

Figure 9.

Size-magnitude diagram for the selected stars (mag>18) in our simulated Imcom and Drizzle images. The size of the stars is measured with the adaptive moments method (Bernstein & Jarvis 2002). Note that Imcom produces a larger output PSF, but it is much more uniform (see main text).

Open in new tabDownload slide

Fig. 10 shows the measured ellipticity of stars found in simulated images as a function of object magnitude. The mean ellipticity components g₁ and g₂ in Imcom are at the level of a ≲few× 10⁻⁴ for all bandpasses. By comparison, the measured shapes on Drizzle images are consistently different from zero by an amount of order ∼10⁻² (Drizzle: SCI). It shows that the Drizzle process is not able to average out the shape of PSFs, regardless of the fact that output pixel size is larger in Drizzle images (0.0575 arcsec pixel⁻¹) than Imcom (0.025 arcsec pixel⁻¹) and pixel size does have an effect on the measured PSF ellipticity (fig. E1 of Troxel et al. 2021).

Figure 10.

The mean ellipticity (g₁, g₂) of PSF candidate stars as a function of star magnitude of corresponding bandpass is shown. Four coloured data points in a panel show the values from the stars found in Y106, J129, H158, and F184 images. These stars were cut out from mosaics of coadded images using the coordinates from truth star catalog. The magnitude is also taken from the truth catalog. The ellipicity is measured with adaptive moments method in GalSim and the error bar is a standard error of the mean. Top: the ellipticiity (g₁) from the cutouts of three sets of images (Imcom and Drizzle coadded SCI images, Imcom coadded truth images. Bottom: g₂ for the same set of stars. We use the Drizzle images produced in Troxel et al. 2023.

Open in new tabDownload slide

Fig. 11 shows the residual ellipticity measured on the cutouts of stars from ‘SCI’ and ‘truth’ images, and the ellipticity compared to the target PSF for Imcom (which in this case is only different from zero due to numerical precision). The difference between ‘SCI’ and ‘truth’ shows the residual shape introduced by the simple detector physics noise model (which includes dark current, saturation, and read and Poisson noise). Here, the average residual ellipticity is ≲2 × 10⁻⁴ for bright stars and ≲8 × 10⁻⁴ for faint stars, verifying that Imcom shows a tolerance to simple noise models on star shapes. We also investigated whether achieving the target fidelity can affect star shapes of different magnitudes. The right column of Fig. 11 that shows the cumulative effect of noise and fidelity indicates that for g₁ the residual behaviour is consistent with the simple noise model case (top left-hand panel of the same figure) whereas Δg₂ shows a non-zero residual ellipticity over a range of magnitude. Although this is likely caused by the output PSF that did not exactly match the target, the leakage into star shapes is less than the PSF ellipticity error requirement.

Figure 11.

The residual ellipticity (g₁: Top, g₂: Bottom) of the stars extracted from several layers of simulated images is presented. Left: the difference in the shapes from Imcom coadded SCI and truth images. This shows the effect of simple noise models on the ellipticity as SCI images include them and truth images do not. Right: the difference in the shapes from the Imcom coadded SCI images and target PSF we chose in our Imcom simulations. This shows the cumulative effect of noise in sky images and output PSF fidelity.

Open in new tabDownload slide

Additionally, we have estimated the magnitudes of the stars from the total image intensity for best-fitting elliptical Gaussian in GalSim. These magnitudes are estimated on the stars of both Imcom and Drizzle ‘SCI’ images. Fig. 12 shows the relative error of the measured magnitude compared to the truth for each filter. It is estimated from the photon flux in the following way. The measured number of photons (N_star) can be used to measure the magnitudes of the stars, using

where N_zero point is the number of photons that Roman would observe for a 0 AB magnitude source (Oke & Gunn 1983) and m_correct is an aperture correction term to account for the fluxes in the wings of stars, because the diffraction wings of the PSF are not captured by the best-fitting Gaussian; thus the HSM fluxes are less than the ‘total’ (integrated to ∞) fluxes. We compute the zero point based on the model throughput curve of Roman that was used in the simulation:

where A_eff is the effective collecting area of the telescope for a given bandpass, |$t_{\rm obs}=139.8\,$|s is the exposure time, and hν is the energy of photons (where h is Planck’s constant). The effective collecting area for a given bandpass can be integrated over the bandpass¹⁰, and ∫(A_eff/λ)dλ  = 0.5915, 0.6051, 0.5978, and 0.3929 m² for Y106, J129, H158, and F184 respectively. We then correct for m_correct by measuring the flux and magnitude on the target PSF for each bandpass; the corrections are 0.0977, 0.1471, 0.2180, and 0.3018 mag for Y106, J129, H158, and F184 respectively.

Figure 12.

The deviation of measured magnitudes of the stars in Drizzle and Imcom images from the magnitudes from the truth catalog for all the bandpasses is shown. The magnitude of the stars is estimated from the flux measurement using galsim.hsm. Since the flux units in GalSim are photons cm⁻² s⁻¹, we multiply the flux measured on Imcom star cutouts by (0.025/0.11)² to obtain flux from surface brightness. There is no need for this scaling for Drizzle cutouts because Drizzle operates on fluxes rather than surface brightnesses. The number of photons is then converted to the AB magnitude through equation (10). The texts in each panel display the interquartile range to indicate the spread that is unaffected by outliers, and the standard deviation of the sample.

Open in new tabDownload slide

Since we specify the target PSF in Imcom to be the Airy disc convolved with the Gaussian kernel, the adaptive moments method captures the fluxes from these stars quite well. Their RMS residual magnitudes are 28, 32, 41, and 51 mmag in Y/J/H/F band whereas the SRD states relative photometric calibration in each filter to be better than 10 mmag. We should note that this measurement includes neighbouring fluxes from nearby sources (which are not included in the SRD requirement) and stars with low fidelity (outliers in Fig. 9). A visual inspection of outliers from the size-magnitude diagram in Y106 showed that most are due to leaking flux from neighbouring objects (not necessarily fully blended). For the adaptive moment magnitudes measured on Drizzle images, on the other hand, as the dispersion in the difference suggests, the coadded PSF obtained through Drizzle is not uniform across our footprint and the magnitude measurement is not reliable in the strongly undersampled case (Y106). (We note that Drizzle preserves total flux by construction; but the fraction of the flux captured by the adaptive moment method depends on the PSF.) Even for the weakly undersampled cases, the mean of the residual magnitude is one order-of-magnitude worse than the Imcom case.

In this section, we have only explored the moments on single-band images; the moments measured on multiband (synthetic wide-band) images will be analysed in Section 5.

4.3 Correlation functions of stars in various outputs

We have so far looked at one-point statistics of observed properties of stars in the coadded images of various input layers. In the end, however, it is crucial to verify that systematic biases related to undersampling and the image reconstruction process do not contaminate the cosmological signal (e.g. fig. 12 of Troxel et al. (2023) for their estimated ξ_± over 20 deg² simulated sky). We approach this by measuring two-point statistics, especially shape–shape correlations of these stars. The calculations of two-point correlation functions were performed using the publicly available TreeCorr package (Jarvis, Bernstein & Jain 2004).

In Fig. 13, we measure shape-shape correlations (ξ₊) of observed stars in Imcom coadded injected stars and ‘SCI’ layers, and in Drizzle coadded ‘SCI’ image. We also show a simple comparison of these star shape correlations to Roman PSF mitigation requirements. Here we compute the approximate requirement on ξ₊ from the requirements on additive errors in SRD in each angular multipole moment bin. We rewrite the Hankel transform of angular power spectra C_ℓ as a Riemann sum (Givans et al. 2022):

where J₀ is the Bessel function of first kind. It is clear here that the shape correlations of Imcom coadded stars are at a statistically acceptable level compared to the requirement. The most idealized case is shown as the correlations of injected stars, drawn with exactly the same input PSF given to Imcom: here any deviations of the measured shapes from isotropy and homogeneity can be attributed to the algorithm itself. We confirm that the level of correlations is two orders of magnitude smaller than the estimated systematic requirement. As can be seen in Fig. 13, disparities exist between ξ₊ on Imcom: SCI and Imcom: injected stars since the stars drawn in the noiseless injected source layer are isolated and the PSFs are made with a flat SED, while the stars drawn in Troxel et al. (2023) are fully chromatic. We should note that we do not expect the measurement using Drizzle stars to be consistent with zero since Drizzle does not attempt to smooth the output PSF to be isotropic.

Figure 13.

Top: Autocorrelation (ξ₊) of observed star ellipticity from various coadded images. This includes Drizzle and Imcom images with simple detector noise models (see Section 2 for detailed explanation), and GalSim-drawn injected stars from Imcom images. The ellipticity was measured with adaptive moments. Dark grey shaded region shows the required ξ₊ signal approximated with the requirements on additive shear errors from SRD. Each panel displays the signals corresponding to each bandpass, and ‘diamond’ marker is the positive signal while ‘square’ is the negative signal but taken absolute value. Bottom: Cross-correlation of the sky coordinates and the observed ellipticity of the same sources. Readers may refer to Fig. 7 for the residual ellipticity around the positions where there are fewer dither positions.

Open in new tabDownload slide

It is worth mentioning that the ξ₊ presented here is the contamination solely from star or PSF shapes – error contributions from residual PSF shape/size errors and shape measurement are not included. We consider the shapes of the injected stars as e_PSF and those of stars in sky image as e_star. Typically, we characterize the additive shear systematics as the following in terms of observed shear, PSF modelling errors and noise:

here |$\delta e^{\rm sys}_{\rm PSF}$| (additive shear systematic biases due to PSF modelling) can be decomposed into

following Paulin-Henriksson, Refregier & Amara (2009) and Jarvis et al. (2016). Our test in this section only quantifies the contribution to e_PSF from residuals in the image combination process, but not the coefficients α, β, and η. This will be investigated in a future paper. We will, however, quantify the contribution to the observed shear due to additive noise biases (δe_noise) in Section 4.5.

Fig. 13 additionally demonstrates the tangential shear around the positions of simulated (in the ‘SCI’ image) stars and the injected stars. While we expect these signals to be zero given that we have not put large scale structure in the position catalog and they seem to be below the additive shear error requirement (only in 2.0 < log₁₀ℓ < 2.5 due to the relevant scale of galaxy-galaxy lensing), γ_t signal around injected stars might be underestimated because the angular distance between sources is discrete.

4.4 Correlations of fourth moments

In addition to the PSF second moments, Zhang et al. (2023c) demonstrated that the correlation function of the spin-2 components of the PSF higher moments (e.g. fourth moment) can cause shear additive bias to cosmic shear correlation function, contingent on the shear estimation method used in the analysis. We measure the standardized higher moments defined in Zhang et al. (2023b),

Here the (u, v) are transformed coordinates from the image coordinate (x, y), such that the second moment shapes in equation (9) vanish, and the second moment size in equation (8) is normalized to 1; p and q are the integer moment indices. For fourth moments, p + q = 4. Here |$\omega (x,y) = {\rm e}^{-{\boldsymbol x}\cdot {\mathbf M}^{-1}{\boldsymbol x}/2}$| is the adaptive weight function given by HSM (Hirata & Seljak 2003), and I(x, y) is the image of the PSF. The complex spin-2 fourth moment is defined as

In Fig. 14, we show the correlation function ξ_± of the PSF fourth moments, measured on the coadded image of the injected stars drawn from the GalSim routine. This is presented with an approximated requirement for the fourth moments correlation function from the SRD, by assuming the fourth moment leakage factor α⁽²⁾ is six times larger than the second moment leakage factor (α in equation (14)). This decision is supported by the results of Zhang et al. (2023c). The fourth moment correlation function of the injected stars are safely within the requirement for the Y106, J129, and H158 band, showing that the higher moments of the anisotropy of the Imcom PSF will not cause significant additive bias in these bands. However, the correlation functions in F184 are only slightly below the requirement for the most angular range and are above the requirement (though not by a statistically significant amount) for the largest scales. This indicates that the PSF fourth moments could potentially contaminate the real cosmological signal in F184. In fact, the fourth moments correlation function is ∼3 order of magnitude larger than the second moment correlation in terms of signal-to-requirement fraction, mainly because Imcom PSFs have extremely low large-scale coherence in the second moment shapes, which is not necessarily the same for the fourth moments.

Figure 14.

Autocorrelation (ξ_±) of the fourth moment spin-2 components (⁠|$M^{\rm (4)}_{\rm PSF}$|⁠), measured by the injected stars drawn from the GalSim routine. Each panel displays the signals corresponding to each bandpass, and ‘diamond’ marker is the positive signal while ‘square’ is the negative signal but taken absolute value. The grey-shaded region shows the requirement on the fourth moment ξ₊, computed from the SRD and the assumption that the fourth moment leakage will be 6-times larger than that of second moments. The PSF fourth moments correlation functions are safely within the requirement for Y106, J129, and H158 band, and are on par with the requirement for the F184 band.

Open in new tabDownload slide

This result suggests that the PSF fourth moments should be actively monitored in the Roman HLIS cosmological analysis. However, further work is required to investigate the leakage of PSF fourth moments into the galaxy shape for the Roman HLIS, by cross-correlating the galaxy shape with the star moments [characterization of the α, β, and η coefficients in equation (14)]. A larger-area simulation will also be valuable here for better understanding the leakage at scales larger than 60 arcmin.

4.5 Estimation of additive noise bias

A concern with measuring shapes on a coadded image is that noise correlations can result in a bias in the shape measurement even if the PSF has been made round and is exactly known (Kaiser 2000; Bernstein & Jarvis 2002; Hirata et al. 2004; Melchior & Viola 2012; Refregier et al. 2012; Mandelbaum et al. 2015; Gurvich & Mandelbaum 2016; Herbonnet, Buddendiek & Kuijken 2017; Okura & Futamase 2018). For example, if there are noise correlations in the x-direction, then there are different centroid uncertainties in the x and y directions; this propagates into different biases in object moments 〈x²〉 and 〈y²〉, and hence a bias in the ellipticity g₁. There are in fact a number of such terms, each resulting from the fact that the ellipticity is a non-linear function of the data, and leading overall to a bias proportional to ν⁻² where ν is the detection significance (see e.g. the detailed discussion in Bernstein & Jarvis 2002, section 8.2). This results in the additive noise bias; the aforementioned leading terms result from the second derivative of the measured shape with respect to the input image, and so we refer to it as the ‘second-order’ additive noise bias. (There are third-order and higher terms as well.) A more detailed treatment will be presented in a future paper; here our aim is a first look at the second-order noise bias of the ellipticities of injected stars.

One may construct a Monte Carlo estimate of the additive bias as follows:

where I is the image of the object; N is a noise image; A and ϵ are scaling factors; and the 〈〉 denotes an average over the Monte Carlo noise realizations. If the shape measurement of the object is performed on a coadded image, then I should be the coadded image and N should be a coadded noise realization (and thus includes any correlations imprinted by the coaddition procedure).

If the object image is normalized to have unit flux, and N is normalized to be a noise image with unit input variance per input pixel, then the proper scaling is

where ν_SE is the single-epoch signal-to-noise ratio and |$\Omega _{\rm psf}/s_{\rm in}^2$| is the input PSF effective area in pixels (shown in table 4 of Paper I). Then equation (17) returns an estimate of the second-order additive noise bias. In principle, the use of multiple Monte Carlo realizations of the noise should allow us to reduce the uncertainty in equation (17). But even a single realization is unbiased and can be used in a correlation function code such as TreeCorr to estimate the additive noise bias correlation function.

There may also be noise bias from the input 1/f noise. Typically the 1/f noise is specified by a ‘knee frequency’ f_knee (with the pixels time-ordered, so units of Hz) where the 1/f noise power spectrum is equal to the white noise component power spectrum. For readout, the 4096 × 4096 pixel arrays used by Roman split into 32 channels of 4096 rows and 128 columns each. The pixels in each row are read out with a cadence of 5 μs; after all the pixels in one row are read out, we move to the next row (so the time to read each row is |$128\times 5 = 640\, \mu$|s plus overheads). This pattern is shown in, e.g. fig. 2 of Freudenburg et al. (2020), and maps 1/f noise into a ‘banding’ pattern (fig. 3 of Paper I). Since our input 1/f noise fields are normalized to unit variance per logarithmic range in f, the appropriate normalization for the 1/f noise fields is

where Δf_band is the bandwidth for the white noise (equal to half the sampling rate, so 100 kHz for Roman); and |$\sigma ^2_{\rm read}/\sigma ^2_{\rm tot}$| is the fraction of the noise variance coming from read noise (as opposed to Poisson noise).

For both types of noise, we first attempted to mesure at single-epoch signal-to-noise ratio ν_SE = 10; if the shape measurement did not converge (see Table 1 for statistics on how often this happened) then we computed Aϵ² using ν_SE = 10^1.5 = 31.6.

We report the noise-induced additive bias (c₁, c₂) for different noise models applied on the noise-less shapes of the GalSim-made injected stars in Table 1. As can be seen, the magnitude of additive biases surpasses the total additive systematic budget of 2.7 × 10⁻⁴ according to the SRD. Furthermore, we correlated Δg (which is the effect of noise in object ellipticity), and as can be seen in Fig. 15 noise-induced additive bias does show correlations over many scales. Hence, in the future, there will be a need to correct for these noise biases with a more comprehensive set of noise images if Imcom is used for image processing for Roman. Possible correction schemes range from numerical approaches based on shearing a noise field (as done in Metacalibration; Sheldon & Huff 2017); subtraction based on second derivatives of the weighted shape estimator (Li & Mandelbaum 2023); or using analytic models such as Appendix  B (presumably allowing the overall normalization to be empirically determined, as suggested in Bernstein & Jarvis 2002). It may be best in a cosmology analysis to implement more than one of these approaches as an internal check. It is clear that for the current set of Imcom settings, the noise-induced biases are larger than residual PSF biases, suggesting that we might benefit from tuning the balance between leakage and noise in the future.

Figure 15.

Autocorrelations of the noise bias (Δg) estimated with white noise (top) and 1/f (bottom) noise models for all the four bandpasses are shown. ‘diamond’ marker is the positive signal while ‘square’ is the negative signal but taken absolute value. These are normalized to single-epoch signal-to-noise ratio ν_SE = 10 and (for the 1/f noise case) a knee frequency of f_knee = 1 kHz. Note that white noise biases generally scale as |$\propto 1/\nu _{\rm SE}^2$| and 1/f noise biases scale as |$\propto f_{\rm knee}/\nu _{\rm SE}^2$|⁠; this means that their correlation functions scale as |$\propto 1/\nu _{\rm SE}^4$| and |$\propto f_{\rm knee}^2/\nu _{\rm SE}^4$|⁠, respectively. These are the ‘raw’ biases, with no attempt at mitigation.

Open in new tabDownload slide

We also show in Table 1 the additive noise biases predicted by second-order biasing theory for Gaussian model objects equation (B26), derived in Appendix  B. (Note that all radii and power spectra must be scaled to the output pixel scale s_out in order to use that result.) The predicted biases are in the correct direction where significant and have approximately the right magnitude (the two largest biases in the table are greater than the prediction by 20 per cent and 37 per cent, although only the latter is statistically significant). The remaining discrepancy may be due to the assumption of a Gaussian profile used in the analytic bias prediction.

5 SYNTHETIC WIDE BAND IMAGES

Synthetic wide band images, created by stacking several standard-width filters, are often used as an intermediate product in weak lensing analyses (e.g. r + i + z in DES (Dark Energy Survey Collaboration; Abbott et al. 2018, 2021). They can be used for making a deep detection image, or as a reference for forced photometry in individual filters for photometric redshifts. (This could include filters on other observatories, e.g. Vera Rubin Observatory data if one were to make a Y106+J129+H158 + F184 image with Roman.) One could also make a shape catalog from the combined Y106+J129+H158 + F184 imaging; such a catalog would have the advantages of providing a deeper catalog and reducing noise biases,¹¹ but the disadvantage that with only one catalog using all the data one can only measure a shape autocorrelation (we would not have the option of cross-correlating two different versions of the shape catalog from different data). It is also possible that one would use both the single-band and the synthetic wide band shape catalogs to test for different systematics–for example, the individual bands allow for cross-correlations to test for PSF systematics associated with the tiling pattern (which is different in each filter; Spergel et al. 2015), but the wide-band image has lower noise bias.

We have built a set of wide band images in post-processing as follows. First, the output images in each filter a are smoothed to a common output PSF Γ_out by convolving with a kernel K_a. Ideally, we want the coadd PSF Γ_a convolved with this kernel to be the output PSF, Γ_a⊗K_a ≈ Γ_out; in practice, we do this by writing

Without the ϵ term, this kernel would exactly transform the PSF Γ_a into Γ_out; but we set ϵ = 10⁻⁴ in the denominator to avoid division by a near-zero quantity. The kernel is clipped to a 201 × 201 output pixel region. We choose the output PSF Γ_out to be an Airy disc convolved with a Gaussian as in Paper I; the Airy disc parameters are λ/D = 0.112 arcsec (the size for J129), obscuration 0.31, and the Gaussian has a scale length of σ = 0.1051 arcsec (the largest of the Gaussian widths used in Paper I). The full width at half-maximum of this output PSF is 0.287 arcsec.

Next, the images are added with some weights:

where H_a is the input image in band a, and H_out is the output image. The normalization factors C_a convert the input units (electrons per |$s_{\rm in}^2$| per exposure) into surface brightness units (μJy |$s_{\rm out}^{-2}$|⁠, or μJy per output pixel) using the effective area curves provided by the Roman project.¹² The weights (summing to ∑_aw_a = 1) are proportional to the inverse square depths given by Akeson et al. (2019): 0.294, 0.323, 0.294, and 0.089 for the Y106, J129, H158, and F184 bands, respectively. Examples of the synthetic wide images are shown in Fig. 16.

Figure 16.

Two examples of fields in the synthetic wide band images (Section 5), from block (0.14). The left column shows the ‘simple’ detector model, the right column shows ‘truth’ (the noiseless image). The top row contains a bright star that saturates on the colour scale (reaching 0.90 μJy |$s_{\rm out}^{-2}$| at its centre), with a log stretch to show the wings of the output PSF. Both images show a 450 × 400 output pixel (11.25 × 10 arcsec) region.

Open in new tabDownload slide

These images are then used again to extract stars to measure their properties, especially the two-point correlation functions. We measure the centroids and shapes of the stars using adaptive moments, and Fig. 17 shows their shape-shape and shape-position correlation functions. Although these signals suggest that the star ellipticity correlations are within the requirement defined in SRD, we do not see an improvement in removing the residual contamination compared to the cases with individual bandpasses shown in Fig. 13.

Figure 17.

Shape-shape (left) and shape-position (right) correlation functions of stars observed in the synthetic wide band images. The grey region shows the estimated requirement on the signals from the PSF ellipticity additive error requirement in the SRD. ‘diamond’ marker is the positive signal while ‘square’ is the negative signal but taken absolute value.

Open in new tabDownload slide

6 SUMMARY AND DISCUSSION

In this paper, we have analysed the coadded images that are produced in Paper I. The layers that were simulated in Paper I include the simulated sky image produced in Troxel et al. (2023), the injected point sources (δ-function convolved with Roman PSF), and two types of noise fields (white noise and 1/f noise). We have measured the following statistics on these images that are relevant to better characterize the weak lensing shear systematics:

• The power spectra of coadded white noise and 1/f noise, as measured in the output images;

• One point statistics of magnitude, astrometric error, second moments (shape and size), and fourth order moments of the injected stars (an idealized case with no blending or chromaticity effects) and stars in the coadded sky image; and

• Two-point statistics (shape-shape ξ_± and position-shape γ_t/x) of shapes and/or positions of the same sources.

If Roman images are to be processed using Imcom, it is essential to understand the implications of noise biasing on measurements of galaxy shapes. Noise power spectra are directly involved in the calculation of shape measurement bias (Appendix  A), so we first investigated the power spectra of the coadded white and 1/f noise fields. In both cases, the noise power spectra of the coadded images decline to zero beyond ∼5–6 cycles arcsec⁻¹, depending on the band. Specific features in the 2D and azimuthally averaged 1D power spectra indicate correlations in the noise fields that will impact shape measurements. We find that the most prominent features in the noise power spectra come from: band limits on the output PSFs, input image pixel positions, and the postage stamp boundaries and angular size. By comparing with analytic models of the expected 1D power spectra for well-sampled images (following derivations found in Appendix  A), we recover some of the expected qualitative features, but the noise power spectra are generally above the simple expectation. This highlights the importance of using full simulations of the image-processing steps to predict the output noise and associated biases, rather than using simple formulae appropriate for well-sampled data. Each filter performs slightly differently, with the increase in power spectrum relative to the simplistic expectation being worst in Y106. This is unsurprising since it has the worst coverage and sampling of all the filters (see Paper I for further discussion of coverage regions).

We have additionally presented the average astrometric offset from input catalog positions and the shape and size deviation from target PSF measured on injected stars (coadded PSF), and compared them to the PSF requirements documented in the Roman SRD. We have shown (Fig. 6) that they are well within the requirements in a root-mean-square sense in the four bandpasses in the Reference survey design (Y106, J129, H158, F184). The Imcom algorithm, although much more computationally expensive, has yielded large gains in output PSF quality relative to the current ‘industry standard’ for combining space images (Drizzle).

Since what matters to weak lensing studies in the end is the level of contamination each systematic effect has in shape correlations, we have measured two-point correlation functions of shapes (second moments: Fig. 13) of the injected and simulated (sky image) stars and fourth order moments (Fig. 14) of injected stars. While the PSF shape-shape correlations are well below the mission requirement for all the bandpasses, the fourth moment-fourth moment correlations in F184 is on the verge of falling short (other fourth moment correlations are not of concern).

We finally investigated the additive bias imprinted on the shapes of the injected stars by both white and 1/f noise in the input images, by adding the appropriate noise fields in and re-measuring the shapes. We find biases that exceed Roman requirements when scaled to signal-to-noise ratio per observation ν_SE = 10, indicating that noise bias will have to be corrected in Roman analyses (as is already planned for Roman and other weak lensing projects). Appendix  B develops an analytic model for the noise bias, which we find is in agreement with the observed trends.

While this set of simulations represents a major step toward use of the Imcom algorithm to the Roman weak lensing program, there are several more studies that should be carried out prior to implementation in an operational pipeline. These include:

• Computational efficiency: The existing implementation is computationally expensive, and would require ∼10⁹ CPU-hours to run on the full baseline HLIS survey if there were no speed-ups, whereas Drizzle would require ∼10⁵ CPU hours (extrapolated from the implementation in Troxel et al. 2023). On the algorithm side, one could investigate optimizing postage stamp and padding parameters (see Paper I), or re-arranging the linear algebra operations in Imcom if we can search over a limited range in the Lagrange multiplier κ.¹³ On the hardware/firmware side, one could make use of ongoing advances in graphics processing unit-accelerated linear algebra.

• Extended source injection: The tests with injected stars in this paper used point sources since these are a ‘stress test’ for undersampling effects. However, we also want to implement grids of extended sources so that we can test metacalibration (Huff & Mandelbaum 2017; Sheldon & Huff 2017) or analytic differentiation (Li & Mandelbaum 2023) techniques wrapped around an algorithm that operates on an imcom coadd.

• Error propagation: We would like to study propagation of astrometric errors, relative flux calibration between images, and PSF model errors through the image coaddition pipeline. This will involve the insertion of specialized layers with these types of errors injected.

• Laboratory noise fields: The analysis of the laboratory noise fields from this project, and an assessment of their implications for additive noise bias, is ongoing. Additional noise fields from the focal plane level and Wide-Field Instrument level tests will be incorporated as they become available.

• Poisson noise bias corrections: While imcom coaddition is a linear operation on the input pixels, shape measurement algorithms are non-linear operations on the data and hence are subject to noise biases. The current implementation allows one to coadd simulated noise relatizations, and therefore one can generate simulated output noise, as needed by tools such as Deep Metacalibration (Zhang, Sheldon & Becker 2023a).¹⁴ Such realizations would also allow a Monte Carlo evaluation of the noise bias terms in twice-differentiable shape measurement algorithms (e.g. Li, Li & Massey 2022; Li & Mandelbaum 2023). The performance of the method has not been tested on simulations with source Poisson noise. But further work will be needed to address biases arising from self-Poisson noise of the source galaxies (e.g. Appendix D of Li & Mandelbaum 2023).

• Chromatic effects: When running image combination in mosaic mode, one has to choose a ‘reference’ intraband SED in order to have a well-defined PSF, but of course the sources in the field have a range of SEDs. Chromatic terms in the PSF will arise from diffraction, dispersion in the filter substrate, chromatic wavefront from mirror, filter, and detector coatings, and depth-dependent absorption effects in the detectors (Mosby et al. 2020). Additionally, there is a field-dependent filter bandpass due to angle of incidence effects, which results in objects redder than the reference SED having a normalization that is larger for input exposures near the field centre and smaller for input exposures near the corners of the field. We plan to test how these chromatic effects propagate through Imcom and assess how they should be corrected.

• Deep fields: Like other weak lensing surveys, the Roman HLIS requires deep fields for photo-z, selection, and noise bias calibration. Some shear calibration algorithms also require deep fields (Zhang et al. 2023a). The Imcom algorithm will run on deep fields, but since the current version has |${\mathcal {O}}(n^3)$| (where n is the number of input pixels), a deep survey with 10 as many exposures but only 1 per cent of the area of the wide survey actually requires 10³ × 0.01 = 10 times as many operations to run Imcom as the wide survey if done by brute force. One solution to mitigate this is to coadd subsets of the deep-field exposures to obtain full sampling at a common PSF, and then do a pixel-by-pixel coadd, but other options should be investigated for feasibility and performance. Another opportunity for the deep fields, since the observations include large dithers, would be to mitigate field-dependent bandpass effects. This would involve extending the linear algebra formalism in Rowe et al. (2011) to simultaneously reconstruct both an image through the ‘mean’ bandpass in each filter and a ‘first principal component’ (PC1) bandpass.¹⁵

• Other survey strategies: The Reference survey strategy was developed to support hardware trades for the Roman mission, and the actual survey design may be different (see e.g. Eifler et al. 2021 for one proposal). Simulations similar to those in this series of papers should be carried out for the possible alternative survey designs to ensure the data will be usable.

ACKNOWLEDGEMENTS

During the preparation of this work, MY and MT were supported by NASA contract 15-WFIRST15-0008 as part of the Roman Cosmology with the High-Latitude Survey Science Investigation Team and by NASA under JPL Contract Task 70–711320, ‘Maximizing Science Exploitation of Simulated Cosmological Survey Data Across Surveys’. EM and CH were also supported by NASA contract 15-WFIRST15-0008, Simons Foundation grant 256298, and David & Lucile Packard Foundation grant 2021–72096. Work at Argonne National Laboratory was supported under the US Department of Energy contract DE-AC02-06CH11357. This research used resources of the Argonne Leadership Computing Facility, which is supported by DOE/SC under contract DE-AC02-06CH11357. RM was supported by a grant from the Simons Foundation (Simons Investigator in Astrophysics, Award ID 620789).

We thank Kirsten Casey for assistance with the azimuthally averaged power spectra, and Eric Huff and Chaz Shapiro for comments on the derivation in Appendix  B. We thank Paul Martini, Chun-Hao To, Peter Taylor, and David Weinberg for feedback on the methodology and analysis choices for this project.

The detector data files used in the Roman image simulations are based on data acquired in the Detector Characterization Laboratory (DCL) at the NASA Goddard Space Flight Center. We thank the personnel at the DCL for making the data available for this project. Computations for this project used the Pitzer cluster at the Ohio Supercomputer Centre (Ohio Supercomputer Center 2018) and the Duke Compute Cluster.

This project made use of the numpy (Harris et al. 2020) and astropy (Astropy Collaboration 2013, 2018, 2022) packages. Some of the figures were made using ds9 (Joye & Mandel 2003) and matplotlib (Hunter 2007). Some of the results in this paper have been derived using the healpy and healpix package (Górski et al. 2005; Zonca et al. 2019).

DATA AVAILABILITY

The codes for this project, along with sample configuration files and setup instructions, are available in the two GitHub repositories:

• https://github.com/hirata10/furry-parakeet.git (postage stamp coaddition)

• https://github.com/hirata10/fluffy-garbanzo.git (mosaic driver)

The version of the code used for this project is in tag 0.1.

Footnotes

References

APPENDIX A: MODELS FOR THE OUTPUT NOISE POWER SPECTRA IN THE ABSENCE OF SAMPLING EFFECTS

This appendix derives the analytic predictions for output noise power spectra in the absence of sampling effects, as used in Section 3. There are two cases for the input noise: white noise, and 1/f noise.

A1 Input white noise

To calculate the expected output power spectrum for white noise, we start by writing the input and output images ({|$I_{a}\rbrace _{a=0}^{N_{\rm in}-1}$| and H_out) as convolutions of the input and output PSFs G_in and Γ_out with the sky scene |$f({\boldsymbol r})$|⁠:

where the input PSF G_in is the Airy disc convolved with a top-hat function, a is the index of the input PSF, and the target PSF Γ_out is the Airy disc convolved with a Gaussian. Taking the Fourier transform and combining these equations for each input image j gives:

where |${\boldsymbol{\upsilon }}=(u,v)$| is the wave vector. For independent input images, the power spectrum of the output is

Here, the Airy disc components of the PSFs cancel each other out when we take the ratio |$\tilde{\Gamma }_{\rm out}/\tilde{G}_{\rm in}$|⁠, so this ratio is simply the ratio of the Gaussian to the sinc (Fourier transform of the top-hat). Following numpy convention and Paper I, we use the definition of the sinc function sinc(ξ) = sin πξ/πξ = j₀(πξ), where j₀ is the zeroth-order spherical Bessel function. Then

In this case, we are interested in input white noise with unit variance in pixels of area |$s_{\rm in}^2$|⁠, so the input power spectrum is |$P_{\rm in}({\boldsymbol{\upsilon }})=s_{\rm in}^2$|⁠.

We now taking an angular average over the direction of the Fourier mode φ ∈ [0, 2π) at fixed magnitude υ. In this case, u = υcos φ and v = υsin φ. After a lengthy but straightforward calculation, we get

where the ‘..’. denote higher terms in the Taylor expansion of υ. Note that the power spectrum has units of area since it is of an input field that is dimensionless by construction. This is the analytic expectation for the white noise power spectrum used in the analysis in Section 3 in this work.

A2 Input 1/f noise

We now consider the case of input 1/f noise. To avoid clutter, we work initially in units of the input pixel scale [up through equation (A20)], and then transform units at the end.

The input 1/f noise fields are generated as 1D time-ordered data (TOD), and re-formatted into a 2D image. Thus for the prediction of an output noise power spectrum for input 1/f noise, we start with the power spectrum as a function of frequency on the 1D pixel axis,

where |$\frac{1}{2}$| is the Nyquist frequency (in units where the time to sample each pixel is unity) and 1/N² is the cutoff for a sequence of N² pixels. Any two pixels separated by some time t have a correlation function

where we have substituted P_TOD(f), used the fact that the f < 0 range contributes the conjugate of the f > 0 range, and set f_max = 1/2 (the Nyquist frequency).

Pixels separated by t in the time series are separated by some (Δx, Δy) in the input data. The 1D (TOD) and 2D correlation functions are thus of the form ξ_2D(Δx, Δy) and ξ_TOD(Δx + NΔy). Summing over (Δx, Δy) gives us a 2D power spectrum (in Fourier space):

In the last line, we have factored out an exponential so that we may substitute, in analogy with the first line,

This gives the 2D power spectrum of the image as a function of the 1D power spectrum of the TOD:

which, using the exponential form definition of the -function (the ‘Shah’ or Dirac Comb function), becomes:

Equation (A11) is general in the case that one considers the full plane, extending to infinity. However, the presence of the -function (and hence δ-function) presents a challenge to numerical computation. The solution is to ‘smear out’ the δ-function according to the uncertainty principle given that our detector array is finite in size. In a finite region, the measured power spectrum |${P}_{\rm meas}({\boldsymbol{\upsilon }})$| is related to the true power spectrum |${P}({\boldsymbol{\upsilon }})$| by

where A is the area of the space and |$\tilde{W}$| is the window function (Fourier transform of the region measured; see discussion in section 13.4 of Press et al. 1992).

We consider a rectangle of width N (in the x-direction) and N_y (in the y-direction) gives A = NN_y and a window function

Applying this to our 2D power spectrum for 1/f noise gives

where we have defined u′ = u − Δu and v′ = v − Δv. Performing the integral over Δv thus takes v′ → v, resulting in

.

The -function can equivalently be written as , which will allow us to trivially perform the integral over Δu as well:

Finally, we substitute in our time-ordered power spectrum given for 1/f noise, equation (A6), and find that the measured power spectrum for 1/f noise should be

To include this result in our 1D representations in Fig. 4, we convert to polar coordinates |$\upsilon = \sqrt{u^2+v^2}$| and |$\varphi =\, {\rm arctan2}(v,u)$| and take the angular average:

We can restrict the range of j since the TOD power spectrum has a Nyquist frequency, so we include only

The integral can be turned into a sum using |$\varphi =2\pi \, K_\varphi /N_\varphi$|⁠, resulting in the final form we use to compute the expectations for the 1D 1/f noise power spectra:

We note that this is in a single input image, in units where s_in = 1. To compare to the power spectra in Section 3.3, we must convert back to general units, include the factor of 1/N_in for the input exposures, and include the square of ratios of modulation transfer functions to convert to the output images [as in equation (A5)]. The result is

We use N = 128 as the width of the output channels (Mosby et al. 2020).

APPENDIX B: ADDITIVE BIAS FROM ANISOTROPIC NOISE

While we have used injected stars with and without noise added to estimate the additive noise bias in the main text of the paper (Section 4.5), it is useful to have an analytic model in order to understand the orders of magnitude in the problem and the scaling behaviours. These analytic models, once calibrated with simulations, are also useful for setting requirements on knowledge of correlated noise—in particular, they were used to determine how many dark images needed to be taken at each calibration epoch to measure read noise correlations in the Roman Design Reference Mission¹⁶ The main objective of this appendix is to predict the additive bias c for ellipticity of stars or of galaxies in terms of the noise power spectrum P_N(u, v). For simplicity, this appendix does so for the adaptive moments method (Bernstein & Jarvis 2002). A related previous study on noise biases in fitting elliptical Gaussians can be found in Condon (1997).

We follow the description of biases in appendix C of Hirata et al. (2004) and appendix A of Refregier et al. (2012). This description applies to the biases in a least-squares fitting algorithm where the galaxy image |$I({\boldsymbol x})$| is fit by a model |$J({\boldsymbol x}|{\boldsymbol p})$|⁠, where |${\boldsymbol p}$| is a point in a parameter space with parameters {p^α} (with N parameters). The metric for the least-squares fit is assumed to be an inner product, denoted by 〈, 〉: that is, the ‘energy functional’ that is minimized is

This definition guarantees the functional derivative

We assume we are working with a well-sampled image (in the context of this paper, that means the output image generated by Imcom) so that summations over pixels instead of integrals, and hence partial instead of functional derivatives, may also be used. The case considered here is more difficult than that in Hirata et al. (2004), because in the Hirata et al. (2004) calculation the least-squares fit is inverse-noise-weighted with respect to the true noise covariance matrix. In contrast, here we are interested in the case where the true noise covariance is not the weight used in the inner product of equation (B1). Indeed, the true noise covariance may not be known exactly, and one of our purposes is to understand what impact errors in the noise model have on output shape measurements.

We further use subscripts on J to indicate partial derivatives with respect to p^α, and drop the arguments unless explicitly required for clarity: thus |$J_\alpha = \partial J({\boldsymbol x}|{\boldsymbol p})/\partial p^\alpha$|⁠.

We first consider the general calculation of mean shifts and biases, and then proceed to the case of galaxy ellipticity.

B1 General results

We note that the minimization of E with respect to a parameter p^α leads to the minimization equation:

This is a set of N nonlinear equations for the N parameters |$\lbrace p^\alpha \rbrace _{\alpha =0}^{N-1}$| in terms of the image I. Clearly if |$I = J({\boldsymbol P})$| for some point |${\boldsymbol P}$| in parameter space, then the solution to equation (B3) is that p^α = P^α.

This solution can be Taylor-expanded. We must first take the functional derivative of equation (B3) with respect to the image:

and then one more derivative:

If we focus next on expanding around the point where I = J, then – defining H^αβ to be the N × N matrix inverse of the symmetric matrix 〈J_β, J_α〉 – we find

and (after some simplification)

Now we suppose that there is some noise covariance matrix |${\mathbf N}({\boldsymbol r},{\boldsymbol s}) = \langle \Delta I({\boldsymbol r}) \Delta I({\boldsymbol s})\rangle$|⁠. Then let us define the symmetric matrix

where |$P_{\rm N}({\boldsymbol{\upsilon }})$| is the noise power spectrum at wave vector |${\boldsymbol{\upsilon }}=(u,v)$|⁠.

The covariance of the model parameters, from equation (B6), is then

and the bias is

After simplifying using the idenitity (provable using the product rule)

this becomes

where from equation (B11):

Note the appearance of the bias tensor, as discussed in Cox & Snell (1968); this is very much like a Christoffel symbol if H is interpreted as a contravariant metric on parameter space. It appears in other contexts where we consider noise biases in photometry (e.g. Portillo, Speagle & Finkbeiner 2020). The first term, involving the derivative of H^γδQ_γδ, does not appear in the usual formula for the bias of a maximum-likelihood estimator; it results essentially from the fact that the weight in the inner product 〈, 〉 used to define the energy functional E is not an inverse-noise weighting with the true noise. In particular, one can show that for white noise with |$P_{\rm N}({\boldsymbol{\upsilon }})=A$| that Q = AH⁻¹, and hence |$\partial (H^{\gamma \delta }Q_{\gamma \delta })/\partial p^\alpha =0$|⁠. But it needs to be taken into account in the present context, where we allow an arbitrary noise power spectrum but then fit objects with uniform weighting. Equation (B12) is additive in the noise correlation function, as is always the case for second-order biases.

B2 Application to ellipticities

We now consider the case where the particular model of interest is the fitting of an elliptical Gaussian:

where the parameters are the flux F, the centroid (x_c, y_c), and the three moments of the 2 × 2 symmetric matrix M, usually written in the form

This is the Bernstein & Jarvis (2002) convention for ellipticities; the scale parameter R was chosen to make M simple, and is defined in terms of the principal axes by R² = (a² + b²)/2. The Fourier transforms are

(it is easiest to compute the inner products in Fourier space).

The derivatives of J are

These lead to the inverse matrix of inner products (where |$e^2=e_1^2+e_2^2$| and ξ = 1 − e², and we have used many applications of Wick’s theorem):

From equation (B8), the noise correlations projected into the data space can be written as an integral over the noise power spectrum:

where the column vector |${\boldsymbol{\tau }}$| is a function of |${{\boldsymbol{\upsilon }}}$|⁠:

We also need to use equation (B13) to evaluate |$\langle J_{\gamma \delta }, J_{e_1} \rangle$| at zero ellipticity:

After a tedious calculation, we find that the total ellipticity bias is

An analogous equation holds for Δe₂:

The combined equation can be written in the form:

where |$\varphi = \arctan (v/u)$| is the position angle of mode (u, v), |$P_{\rm N}(u,v)\, {\rm d}u\, {\rm d}v$| is the differential contribution to the noise intensity variance |$\sigma _I^2$|⁠, and the function

is a dimensionless weight (Fig. B1). The weight function goes to zero both at z → 0 and z → ∞: this is because correlated noise at zero spatial frequency is equivalent to an overall ‘pedestal’ in the background at each galaxy, which may affect the flux measurement but does not bias the ellipticity. Noise at high-spatial frequency υ ≳ 0.6/R is outside the band limit of the shape measurement, and hence has no effect.

Note that the above considerations apply to the ‘as-observed’ galaxy, e.g. ‘R’ is the scale radius after convolution with the PSF, and (Δe₁, Δe₂) is the noise-induced ellipticity bias before PSF correction. In for example a shape measurement method where the observed ellipticity is ‘corrected’ by dividing by a responsivity factor |${\cal R}_2$| (usually between ∼0.3 and 1), the bias on the reported ellipticity would also be divided by this factor. Also note that one must divide by 2 to get an ellipticity in the g = (a − b)/(a + b) convention (where a and b are major and minor axes) instead of the e = (a² − b²)/(a² + b²) convention:

The analytic formula is useful since it is applicable to any shape of the noise power spectrum: weakly anisotropic noise, ‘striping’ noise (correlated in the row direction), and at any spatial scale. This allows us to rapidly estimate the order of magnitude of the bias introduced by a noise source and understand its scaling behaviour. However, it is also critical to compare the analytic model to simulations.

B3 Comparison to toy simulations

It is useful to compare the analytical bias estimate to a ‘toy’ simulation designed to mimic the assumptions in the analytic model. This is a check of whether the analytical calculation is performed correctly and also allows simple tests of its range of validity.

The first set of such simulations is performed on a circular galaxy with flux F = 200 counts, scale radius R = 4 pixels, and centred at pixel (32,32) of a 64 × 64 simulation. A Gaussian random field was generated for the noise, with power spectrum:

where A is an amplitude parameter (with units of counts per square pixel), and the default values of the constants are α = 4/π and β = 16/π. This results in a noise field correlated in the y direction; see Fig. B2.

Figure B1.

The weight function W(z) in equation (B24) as a function of wave number. Note that an anisotropic contribution to the noise is most damaging at a scale of |$\upsilon =\sqrt{u^2+v^2} \approx 0.3/R$| (or a wavelength of υ⁻¹ ≈ 3R).

Open in new tabDownload slide

Figure B2.

An example circular galaxy with flux F = 200 counts and scale radius R = 4 pixels, superposed on the noise field of equation (B27) with A = 1, α = 4/π, and β = 16/π. That is, the correlation length of the noise is twice as large in the y-direction as the x-direction. Noise of this form leads to a bias of e₁ toward negative values, i.e. toward inferring a vertical orientation for the galaxy.

Open in new tabDownload slide

We may now investigate the biases in ellipticity resulting from this noise field. We generated 10⁶ random realizations of the noise, and measured the galaxy shape each time using the least-squares elliptical Gaussian method. We compare the ellipticity bias with the analytical prediction of equation (B22) for a range of galaxy parameters. For the noise spectrum of equation (B27), the predicted bias can be obtained via integration over the Gaussian:

and this can be compared to the numerical results.

The outcome of these tests is shown in Fig. B3. We see that, as expected, the noise bias scales as the noise variance ∝ A² (panel a) and as the inverse square of the galaxy flux, ∝ 1/F² (panel b). The radius dependence is more complex (panel c) – indeed, from equation (B28), one can see that the bias approaches a constant at large R, i.e. when the size of the galaxy is large compared to the correlation length of the noise. This is because for fixed flux, the reduction in signal-to-noise as the light is spread over more pixels cancels against the decreasing importance of the noise correlations (since they are over a small fraction of a galaxy scale length). The noise anisotropy (β/α, or aspect ratio squared) is varied in panel (d), with fixed α. Good agreement with the analytical model is seen in all of these cases.

Figure B3.

The shape measurement biases measured from simulations (red points with 1σ Monte Carlo errors) and calculated from the analytic approximation, equation (B22) (blue lines). Each panel is based on the reference case (circular galaxy, F = 200 counts, R = 4 pixels, noise at A = 1 with the noise spectrum of equation (B27) except for the variations indicated. In panels (a–d), the agreement is in general excellent, but for small signal-to-noise ratio, statistically significant deviations can be seen. Panels (e) and (f) show cases that violate the assumptions of the analytic model, and hence show larger errors. The number of simulations per point is 10⁶, except for the 0 < η < 1 points in panel (f), which are more computationally expensive; in these cases 10⁵ simulations were run.

Open in new tabDownload slide

Panels (e) and (f) of Fig. B3 show effects not treated in the analytical model. In panel (e), we have varied the ellipticity of the galaxies—that is, instead of inserting a circular Gaussian galaxy, we have given it the indicated ellipticity |$e=\sqrt{e_1^2+e_2^2}$| and assigned a random position angle |$\chi = \frac{1}{2}\arctan (e_2/e_1)$|⁠, and kept the same flux and scale radius R. It can be seen that for more elliptical galaxies, the bias becomes ‘smaller’ (in absolute value, i.e. closer to zero).

Panel (f) shows the variation with the galaxy profile. For each value of η, a galaxy was generated via the convolution of an exponential profile and a Gaussian profile:

and with a ratio of radii r_e: r_G = η: 1 − η. This interpolates between the pure Gaussian limit (η = 0) and pure exponential (η = 1). This is intended to be a crude representation of a centrally peaked galaxy with some smearing. The normalization (amplitude) and total radius are set so that the least-squares Gaussian fit has the default values of F = 200 and R = 4.