Simulating image coaddition with the Nancy Grace Roman Space Telescope – I. Simulation methodology and general results

ABSTRACT

The upcoming Nancy Grace Roman Space Telescope will carry out a wide-area survey in the near-infrared. A key science objective is the measurement of cosmic structure via weak gravitational lensing. Roman data will be undersampled, which introduces new challenges in the measurement of source galaxy shapes; a potential solution is to use linear algebra-based coaddition techniques such as imcom that combine multiple undersampled images to produce a single oversampled output mosaic with a desired ‘target’ point spread function (PSF). We present here an initial application of imcom to 0.64 square degrees of simulated Roman data, based on the Roman branch of the Legacy Survey of Space and Time (LSST) Dark Energy Science Collaboration (DESC) Data Challenge 2 (DC2) simulation. We show that imcom runs successfully on simulated data that includes features such as plate scale distortions, chip gaps, detector defects, and cosmic ray masks. We simultaneously propagate grids of injected sources and simulated noise fields as well as the full simulation. We quantify the residual deviations of the PSF from the target (the ‘leakage’), as well as noise properties of the output images; we discuss how the overall tiling pattern as well as Moiré patterns appear in the final leakage and noise maps. We include appendices on interpolation algorithms and the interaction of undersampling with image processing operations that may be of broader applicability. The companion paper (‘Paper II’) explores the implications for weak lensing analyses.

1 INTRODUCTION

Weak gravitational lensing – the distortion of the shapes of distant galaxies as their light passes through the gravitational potential of foreground structures – has emerged as one of the powerful tools for probing the growth of structure in the Universe [see Hoekstra & Jain (2008), Weinberg et al. (2013), and Mandelbaum (2018) for recent reviews]. It has now been more than two decades since the early detections of lensing by galaxies (Brainerd, Blandford & Smail 1996) and of the two-point correlations of shear (Bacon, Refregier & Ellis 2000; Kaiser, Wilson & Luppino 2000; Van Waerbeke et al. 2000; Wittman et al. 2000). In that time, both the size of cosmological surveys and the understanding of the instruments and data processing necessary to extract the weak lensing signal has advanced considerably. The current generation of weak lensing surveys – the Kilo Square Degree survey (KiDS), the Dark Energy Survey (DES), and the Hyper Suprime Cam (HSC) – have measured the amplitude of cosmic structure S₈ to precision better than a few per cent (Hikage et al. 2019; Hamana et al. 2020; Asgari et al. 2021; Amon et al. 2022; Secco, Samuroff et al. 2022). With these surveys, and the detailed picture of the initial conditions for the growth of structure provided by cosmic microwave background observations (Planck Collaboration 2020), it is now possible to do detailed comparisons of low redshift structures to the predictions of general relativity in the ΛCDM (cosmological constant + cold dark matter) model and various alternatives (e.g. Lemos et al. 2021).

In the 2020s, several major surveys are expected to come online that will lead to large advances in both the statistical constraining power and control of systematic uncertainties in weak lensing measurements. These include the ground-based Vera Rubin Observatory, which will carry out a ‘Wide, Fast, Deep’ survey of the optical sky in the ugrizy bands (LSST Dark Energy Science Collaboration 2012; Ivezić et al. 2019); the Euclid satellite, which features a wide-field visible channel for galaxy shape measurements and a near-infrared (NIR) instrument to improve photometric redshifts (Laureijs et al. 2011); and the Nancy Grace Roman Space Telescope, whose Wide Field Instrument will survey the sky at high angular resolution in four NIR bands spanning 0.9–2.0 μm (in the Reference Survey design;¹ Spergel et al. 2015; Akeson et al. 2019). All of these observatories bring advantages relative to current surveys in terms of systematic error control for galaxy shapes: Rubin will obtain hundreds of images of each field, thus enabling much better internal constraints of instrument systematics, while the space missions provide the high angular resolution and stability possible above the Earth’s atmosphere. Furthermore, the coverage from u band in Rubin through ∼2 μm for the space missions (with Roman data reaching NIR depths comparable to Rubin sensitivity, and Euclid data being shallower but covering a wider footprint) provides an enormous wavelength baseline for photometric redshifts (Hemmati et al. 2019; Newman & Gruen 2022).

Although space is in many ways an ideal location for a weak lensing experiment, the high angular resolution brings some challenges related to the pixel scale. A diffraction-limited telescope produces a point spread function (PSF) with a characteristic angular width of λ/D, where λ is the wavelength of observation and D is the diameter of the entrance pupil. If the pixel scale P is larger than λ/(2D), then the image is undersampled in the Nyquist sense: there are Fourier modes present in the image with more than |$\frac{1}{2}$| cycle per sample, with the consequence that the images cannot be unambiguously interpolated, and thus the image intensity |$I({\boldsymbol r})$| cannot be treated as a continuous field. It also produces biases in the moments of images, including the first moment (centroid, relevant to astrometry) and the second moments (sizes and shapes, relevant to weak lensing), which have been studied in many contexts (e.g. Lauer 1999b; Anderson & King 2000; High et al. 2007; Samsing & Kim 2011). One possible solution to the undersampling problem is to simply use small pixels – but if one has a fixed pixel count in the focal plane (often limited by available resources or technical considerations), then shrinking the pixels leads to a smaller field of view and a slower survey if we fix the survey depth and hence the required exposure time per pointing.² An alternative, adopted for Euclid and Roman, is to accept undersampling, and use multiple dithered exposures of each field. This increases survey speed, but requires the development of algorithms for each application that take multiple undersampled images as input.

Undersampling can affect several stages of a weak lensing analysis (e.g. Kannawadi, Rosenberg & Hoekstra 2021; Finner et al. 2023). One particular step is calibration of the shear estimator – that is, determining how the measured ellipticity of a galaxy e_i in a catalog responds to an applied shear γ_j.³ While one might hope for an ellipticity measurement algorithm that has unit shear response ∂〈e_i〉/∂γ_j = δ_ij, any stable ellipticity estimator has a response that depends on the galaxy population (Massey et al. 2007; Zhang & Komatsu 2011), and so modern lensing analyses contain a ‘shear calibration’ step that determines this response given the ensemble of galaxy morphologies present in that survey (at that depth, resolution, and observed wavelength). Whether the data are well-sampled or not, shear calibration must work with the fact that Fourier modes in the image are only well-measured up through some k_max, and the shear operation moves modes across the k = k_max boundary (e.g. Bernstein 2010), so one cannot take an observed image and infer what the sheared image would look like at the same resolution. In the past decade, several principled shear calibration approaches have been introduced that apply some re-smoothing (or in Fourier space, a cut k_cut < k_max) before measuring galaxy shapes (or moments) and have been successful at mitigating the galaxy population-dependent shear calibration biases in simulations. These include Metacalibration, which numerically applies a shear to each object to compute the ensemble response (Huff & Mandelbaum 2017; Sheldon & Huff 2017; Zhang, Sheldon & Becker 2023); the Bayesian Fourier Domain technique, which builds the probability distribution of Fourier-space moments (Bernstein & Armstrong 2014; Bernstein et al. 2016); and approaches that analytically build shear responses (Li & Mandelbaum 2023).⁴ Further development of these methods is anticipated in support of final analyses of the Stage III ground-based survey data and the upcoming Vera Rubin Observatory. However, all of these techniques rely fundamentally on operations such as cuts in Fourier space that cannot be performed on undersampled data. Indeed, the first attempts to simulate Metacalibration for undersampled images resulted in per cent-level biases (Kannawadi, Rosenberg & Hoekstra 2021; Yamamoto et al. 2023) that exceed requirements for the upcoming surveys. This motivates us to develop image processing techniques to recover full sampling, as a step in the lensing analysis that precedes shape measurement (and even source selection), so that these shear calibration techniques can be applied to Roman data.

The shear response is a property not just of the shape measurement algorithm, but of the sample selection as well, since applying a shear to a galaxy could cause it to cross a selection threshold (Kaiser 2000; Hirata & Seljak 2003). The same tools that have been developed to measure the shear response can be extended to incorporate the response from source selection (e.g. Sheldon et al. 2020, 2023). Furthermore, in a tomographic weak lensing analysis, the assignment of a galaxy to a particular tomographic bin is itself a form of selection, and thus one needs to derive the shear response of the photometry in each filter used for constructing the tomographic bins (e.g. Troxel et al. 2018; Gatti et al. 2021; Myles et al. 2021). The next-generation weak lensing surveys have specified ‘shape measurement’ filters (J129+H158 + F184 in the case of Roman) as well as additional filters that are used for photometric redshifts (Y106 for Roman). This means that although the source galaxies do not need to be resolved in the photo-z-only filters, one still needs to either recover full sampling in these filters (and use the aforementioned frameworks) or develop an alternative framework for determining the statistical shear response. (This issue was not fully understood at the time the Roman Reference Survey was designed.) Therefore, this paper has investigated recovery of a fully sampled mosaic in Roman Y106 band as well.

The recent joint Roman + Rubin simulations (Troxel et al. 2023), built on top of the Legacy Survey of Space and Time (LSST) Dark Energy Science Collaboration (DESC) Data Challenge 2 (DC2) simulations (Korytov et al. 2019; LSST Dark Energy Science Collaboration. 2021a,b; Kovacs et al. 2022), provide an opportunity to test out algorithms to recover full sampling as part of an integrated simulation and processing pipeline suite, and explore the implications for Roman weak lensing analyses. This is the first in a series of papers that develops and tests an implementation of the imcom algorithm (Rowe, Hirata & Rhodes 2011) on a simulation of part of the Roman Reference Survey. In this paper (‘Paper I’), we focus on the characteristics of the simulation, relevant mathematical background, image combination machinery, and basic properties of the outputs. The companion paper (‘Paper II’) covers statistical analyses of the output images, including the ellipticities of simulated stars, correlation functions, noise power spectra, and noise-induced biases. Both papers use a 48 × 48 arcmin region from the simulations, large enough to contain ∼2 Roman fields of view and a representative portion of the tiling pattern (see Fig. 1 ).

This paper is organized as follows. The problem of combining images to create fully sampled output with uniform PSF, and the goals of this simulation effort, are described in more detail in Section 2. The input data, including ancillary information such as masks, are described in Section 3. The image coaddition simulation methodology is described in Section 4. Some of the basic results on output PSF and noise properties are described in Section 5, and we conclude in Section 6. Appendix  A presents some useful results on optimal interpolation methods for functions with known oversampling rates that we derived for this work, and that may be of more general interest, while a description of the computing resources for the project can be found in Appendix  B. A summary of how undersampling interacts with the operations used in modern shear calibration methods is provided in Appendix  C.

2 STATEMENT OF THE PROBLEM

A variety of methods have been explored in the literature for combining multiple images of the sky into a single ‘coadded’ image. If the coadded image is to be used as the starting point for a standard weak lensing analysis, one wants it to be both oversampled and have a well-defined PSF in the sense of Mandelbaum et al. (2023): the output image should be the true astronomical scene convolved with a PSF. We will go further here and ask for an algorithm that produces a specific desired output PSF Γ that is uniform and circular – thus we choose the output PSF and design a coaddition scheme to achieve it, rather than running a coaddition code and accepting the measured PSF at the end. This is advantageous for survey uniformity and mitigation of additive biases, but as we will see this is only possible under certain circumstances.

Thus in the context of this paper, the goal of the coaddition pipeline is to take in several input images of the sky, which are at some native pixel scale s_in and in general have their own rotations, distortions, PSFs, and masks; and produce a well-sampled output image at an output pixel scale s_out, and with a uniform, round output PSF. Implicit in this statement of charge is that we are trying to accomplish at once several tasks that are sometimes distinct steps in an image processing pipeline:

1. Interpolation over masked pixels (e.g. cosmic ray hits, bad columns, or hot, dead, or unstable pixels).

2. Resampling onto a common grid.

3. Rounding and homogenization of the output PSF (for some weak lensing pipelines; this step could also be performed last or not at all).

4. Averaging of the intensities from each input image to yield a single output image.

For oversampled data, it may be reasonable to treat these as separate, since operations such as resampling, convolution, and (sometimes) filling in a single missing sample, can be carried out on an oversampled function without introducing biases. However, they may be viewed in a unified framework if all of the operations used are linear. In this case, each step is a matrix operation, leading to an output image that is a linear combination of input images, with the mapping described by an m × n matrix T, where m is the number of output pixels and n is the number of input pixels. These statements apply to many of the common algorithms that have been applied to ground-based weak lensing data sets (either for shapes, source selection, or photo-zs): for example, linear predictive codes for interpolating bad pixels (section 4.5 Bosch et al. 2018); Lanczos-3 (Bertin et al. 2002; Bosch et al. 2018, section 3.3) or polynomial (section 4.4 Huff et al. 2014) interpolation for resampling; and pixel-domain (section 7 of Bernstein & Jarvis 2002) or Fourier-domain (section 4.1 of Huff et al. 2014) rounding kernels, and PSF Gaussianization (Hildebrandt et al. 2012, section 3; Kuijken et al. 2015, section 4.2).

None of steps 1–3 are possible individually on undersampled data without introducing biases or making additional assumptions about the astronomical scene. However, the end goal of constructing a well-sampled output image with uniform PSF with some matrix T may still be possible, even if T cannot be factored into the individual steps. The Fourier-domain algorithm of Lauer (1999a) and the iterative algorithm of Fruchter (2011) are examples that combine some of these steps for some types of input data. The imcom technique of Rowe, Hirata & Rhodes (2011) searches for a matrix T that accomplishes all four steps with minimum noise and target PSF error (as measured with quadratic metrics).⁵ It is computationally expensive but can work with rolls, geometric distortions, varying input PSFs, and complex masks, and thus is a promising choice for a space-based weak lensing survey with Roman. imcom returns a residual estimate for the output PSF; this allows us to identify cases where the desired output PSF is impossible to build (e.g. due to aliasing with insufficient dithers, or contains Fourier modes not represented in the input image).

We note that the Drizzle algorithm (Fruchter & Hook 2002) that is commonly used to combine undersampled space-based images is itself a linear operation that can be described using a coaddition matrix T. In the case of Drizzle, T is sparse; the entry T_αi is determined by the overlap of output pixel α with a ‘shrunken’ version of input pixel i. The Drizzle matrix T is within the search space for imcom, and therefore by imcom’s target metrics of sum-of-squares error in the PSF and noise, imcom will always perform at least as well as Drizzle (usually much better). But by choosing a particular sparse T, the Drizzle algorithm has a lower memory footprint and much faster run time, and therefore is likely to remain useful in the weak lensing analysis as a ‘quick look’ tool (this is also how it was used in the DC2 simulation; Troxel et al. 2023) and for flagging purposes.

While Rowe, Hirata & Rhodes (2011) demonstrated their method on some test problems with ∼6 × 6 arcsec postage stamps, and the method has also been demonstrated on laboratory data (Shapiro et al. 2013), the method has not yet been applied to Roman simulations over an area large enough to measure the statistical properties of the coadded images. The availability of the Rubin Data Challenge 2 (DC2) + Roman simulation suite, and updated knowledge of Roman properties including characterization of the flight detectors, makes this an excellent time to embark on such a simulation. The qualitative and quantitative objectives for this simulation are as follows:

1. Wrap the algorithm in a driver that can tile the sky and pipe the appropriate (simulated) observations to the linear algebra kernel and re-assemble the output into a community standard format such as FITS.

2. Find an example (not necessarily final) set of parameters for the imcom algorithm that avoid basic problems such as noise amplification or ghosting across postage stamp boundaries (or reduce them to acceptable levels).

3. Characterize how chip gaps and cosmetic defects map into the assembled mosaics.

4. Use injected sources to test the output normalization, astrometry, size, shape, and higher moments of the final coadd PSF after propagation through imcom.

5. Measure the correlation function of the ellipticities at scales overlapping the range likely to be included in the Roman weak lensing analysis.

6. Measure the noise properties of the output image for both uncorrelated and correlated input noise.

7. Compare the noise measured on the output images to the predictions from the Exposure Time Calculator (Hirata et al. 2013) and analytical descriptions in the literature (e.g. Bernstein 2002).

8. Determine the as-realized computing time requirements for the coadd on a modern computing cluster (in order to inform both resource estimates and priorities for further optimization).

9. Compare the moments (i.e. shape and size) of bright unsaturated stars (suitable for PSF characterization) in the output images with the measurements performed on Drizzled coadd from DC2 + Roman simulations.

Most of these objectives can only be achieved by simulating >1 field of view. This paper presents the initial set of coaddition simulations and first look results; follow-on papers will go into more detail on some of the individual objectives.

3 INPUT DATA

We generate several types of input data: full and noiseless image simulations, noise fields, and grids of injected sources. All of these are common types of inputs to image processing pipelines in weak lensing data analysis. Since the coaddition process is linear, we may use the distributive property to superpose outputs and generate, e.g. an output sky image with additional 1/f noise, or an injected object in the real survey (e.g. Suchyta et al. 2016). The processing of all layers at the same time allows the T matrix to be computed only once, a major advantage since it is both the most computationally demanding step and is so large that only a few examples can be stored to disk rather than the T for the whole survey. The layers used in this analysis are shown in Table 1.

3.1 The Roman + Rubin simulations

The principal input data to this study is the Roman arm of the joint Roman + Rubin DC2 simulation, described in Troxel et al. (2023). The simulation begins with a cosmological N-body simulation (Heitmann et al. 2019) populated with galaxies (Benson 2012; Hearin et al. 2020). This forms the basis for the cosmoDC2 simulated extragalactic sky (Korytov et al. 2019; Kovacs et al. 2022), which includes a superset of the region simulated in this paper. This is integrated with a local Universe simulation and a telescope + instrument simulation for the Rubin Observatory (LSST Dark Energy Science Collaboration 2021a,b). The Roman arm of the image simulation observes the same sky, but with the current version of the Roman telescope + instrument simulation framework (an update of Troxel et al. 2021).

The simulation includes a portion of the Roman Reference Survey tiling strategy, which observes the sky in four bands: Y106, J129, H158, and F184. Each band has two passes over the sky at different roll angles, with small-step diagonal dithers to cover the chip gaps. The number of dither positions at each roll was based on preliminary estimates of the number of dithers required to mitigate sampling issues (Rowe, Hirata & Rhodes 2011; Spergel et al. 2015): four in J129, and three in each of H158 and F184. The Y106 band strategy was not required to achieve full sampling and so used three positions. The resulting coverage pattern is shown in Fig. 1. A more in-depth description of the sky tiling can be found in appendix A of Troxel et al. (2023).

The main input layer used for this simulation is the ‘simple’ model sky images from Troxel et al. (2023). These include the convolution of sky objects with the PSF and Poisson noise. They do not include the full detector physics model (this version was chosen so that we can test downstream processing steps, such as image combination, before corrections for the detector effects are ready). The two most important parts of the detector physics missing for our purpose are (i) that there is no pixel mask in the ‘simple’ simulation (we introduce our own in Section 3.2); and (ii) the charge diffusion is not included. Since charge diffusion smears out the effective PSF before pixelization, it actually improves sampling, so we expect that ignoring it is conservative from the perspective of a code that attempts to reconstruct a fully sampled image. This expectation is explicitly tested for a small subsample of the data in Section 5.4.

The ‘simple’ input model also comes with a PSF model (corresponding to a flat spectral energy distribution or SED in F_λ) and an astrometric solution (the World Coordinate System or WCS in the FITS header). The WCS solution in the simulation is ‘perfect’ in the sense that we are working with the same WCS used to draw the image; propagation of systematic errors in the WCS solutions will be considered in a future paper.

We have also included as additional layers the ‘truth’ (noiseless) images from the simulation, and – for J129 and H158 – the ‘err’ HDU (the realization of background Poisson noise used in the simulation).

3.2 Masks

The ‘simple’ DC2 + Roman simulations produce output for all of the pixels. However, in the real mission there will be masked pixels. These can be divided into two types: pixels that are defective (e.g. disconnected, hot, unstable) and are flagged by the calibration pipeline; and pixels affected by a cosmic ray in that particular exposure. These two classes of masks affect image combination differently, since defective pixels tend to be spatially clustered and are the same pixels when we dither, whereas the cosmic ray impacts randomly knock out small groups of pixels.

For this simulation, we have taken the ‘permanent mask’ based on the SCA files constructed from acceptance testing (appendix B Troxel et al. 2023). Out of the 18 chips initially selected for flight, an average of 3.01 per cent of the pixels were flagged (BADPIX HDU) by at least one of the steps in the pipeline. The sources of these flags in order are shown in Table 2. This should be considered only a first approximation to the likely mask used in flight, since the hot pixel populations change as the detector ages or is exposed to radiation (e.g. Sunnquist, Brammer & Baggett 2019).

The cosmic ray mask is based on a random number generator, with the seeds chosen based on the observation ID and SCA so that the same mask is generated for each observation even if it is called in a different block. The density of cosmic ray strikes is taken to be 7.7 × 10⁻⁴ per pixel, which is the product of the |$10^{-6}\,$| cm² pixel area, the 140 s exposure time, and the expected rate of 5.5 events cm⁻² s⁻¹ expected for Roman (see Kruk et al. 2016, and note that at L2 we do not have the trapped electron population). While experience with similar detectors on the James Webb Space Telescope is that usually 1–2 pixels are affected (section 6.6 of Rigby et al. 2023), Roman has smaller pixels (so potentially a higher leakage into neighbours via charge diffusion) and we are aiming for higher precision. Therefore, for this simulation, we have masked a 3 × 3 pixel region surrounding each cosmic ray. The less common larger events such as ‘snowballs’ (e.g. Rieke et al. 2023, fig. 15) are not included in the current simulation; we plan to add them in the future once we understand better how they should be implemented. However for the purposes of testing how imcom responds to masked pixels, it is the more frequent events that are likely to have the greatest impact.

An example of a mask is shown in Fig. 2.

3.3 Injected sources

We have also created two layers that are grids of injected stars. The stars have unit flux and are injected on a healpix⁶ grid (Górski et al. 2005) of resolution 14 (nside|$\, =2^{14}=16384$|⁠) in the equatorial coordinate system. The selection of grid points that land on each SCA was performed with healpy routines (Zonca et al. 2019). Two layers are generated – one with the external GalSim package, and one with the internal interpolation machinery in our pipeline.

We constructed the grid of injected sources with galsim as follows. For each PSF that is characterized by the observation and SCA that overlaps with the output region, we interpolated the PSF made by the DC2 + Roman simulation, which is oversampled by the factor of 8, using the interpolant, lanczos50. The interpolated PSF image is then convolved with the delta function of unity, and is drawn at each grid point.

The internal simulation was generated from the same PSF, but interpolated using the ‘D5, 5, |$\tfrac{1}{12}$|’ interpolation kernel. In principle one should get the same results from the two methods aside from the choice of interpolation kernel or treatment of edge effects; however, having both serves as an important cross-check on the implementation.

3.4 Simulated noise fields

We construct two types of simulated noise fields: white noise and 1/f noise. Since the image combination is a linear process, the proper normalization of these noise fields, or the choice to include both, can be chosen in post-processing. We therefore make the simplest choices to normalize the inputs.

The white noise input is a 4088 × 4088 Gaussian random field with mean 0 and variance 1.

The 1/f noise input is constructed as follows. For each of the 32 readout channels, we construct a 1D Fourier domain signal of length N = 2²⁰: |$\lbrace \tilde{S}_k\rbrace _{k=-N/2}^{N/2-1}$|⁠, where the real and imaginary components of |$\tilde{S}_k$| are Gaussians with mean 0 and variance 1/(2|f_k|) (except for the constant mode |$\tilde{S}_0=0$|⁠), where f_k = k/N is the normalized frequency. This is then discrete Fourier transformed, |$S_j = \sum _{k=-N/2}^{N/2-1} \tilde{S}_k {\rm e}^{2\pi {\rm i}jk/N}$|⁠, and we take the real part. This gives a signal whose variance per logarithmic range in frequency is 1, i.e. the variance coming from all modes between f_min and f_max is ≈ln (f_max/f_min). We take a length 2¹⁹ portion of this signal (to avoid wrapping effects), reformat it into a 4096 × 128 array, and broadcast it into the appropriate portion of a 4096 × 4096 image. The even channels are flipped left-to-right (see e.g. fig. 2 of Freudenburg et al. 2020 for the ordering of the readout pattern). This does not include guide window or row overheads, but it does provide a noise field with the characteristic horizontal banding pattern of 1/f noise that we can use to assess the impact on weak lensing with coadded images. An example of a 1/f noise field generated in this way is shown in Fig. 3.

Both noise inputs use the numpy.random.default_rng random number generator, with a seed chosen based on the observation ID, the SCA, and an integer specified by the user. In order to reduce storage requirements, the noise fields are generated on the fly at the same time the science data is read. The seed construction ensures that when we make a mosaic coadd, and a given SCA image contributes to more than one tile of the mosaic, that the same noise realization is generated each time.

3.5 Laboratory noise fields

In addition to simulated noise fields, we construct a layer containing noise from laboratory testing of the SCA detectors. By including detector read noise in our analysis, we are able to test in a new way how the detectors themselves might impact galaxy shape measurement. The lab detector test noise fields are intended for a separate work, and thus a more complete and detailed analysis will be forthcoming in a companion paper to this. Here we will present the basic principles of the laboratory noise fields, and refer the interested reader to this future paper for further details.

The lab noise data, taken at the Detector Characterization Laboratory, are N_t × 4096 × 4224 cubes (where N_t is the number of time slices), including dark frames and low and high level flat frames for each SCA. Data is split into 32 readout channels, plus one reference output channel and four rows of reference pixels on each side of the grid. For construction of the noise frames, lab data is averaged into six effective time bins (‘Multi-Accum’ processing; see e.g. section 7.7 of Dressel 2022), slope fitted, and reference pixel-corrected, and converted to electrons using the gain determined using solid-waffle (Freudenburg et al. 2020). We impose each SCA’s lab-tested pixel mask (appendix B Troxel et al. 2023) on each frame so that the mask is the same in all input layers. This ensures that we are still able to process all the layers at the same time and thus only compute the T matrix once. For the total focal plane, this additional masking amounted to 0.095 per cent of the pixels. Lab noise field analysis and results will be presented in detail in the forthcoming paper (Laliotis et al., in preparation).

4 IMAGE COADDITION

The image coaddition process is an updated version of the imcom framework for coaddition of postage stamps (Rowe, Hirata & Rhodes 2011). Major changes include replacing the original Fortran 90 code with a python interface and C back end; and the new driver script to produce mosaic coadds.

The mosaic process is organized hierarchically, controlled by keywords in the configuration file (see Fig. 4). The top level is a mosaic, consisting of a single world coordinate system over a region of sky small enough to be treated as ‘flat’ to within weak lensing requirements (i.e. the coordinate system shear and plate scale variation should be <10⁻⁴). The mosaic is divided into a square grid of blocks; the processing of each block in each filter is a single run of the python script run_coadd.py, and the resulting coadded images and metadata are stored in a single FITS file. The blocks are themselves divided into a square grid of postage stamps, and are extended by some number of postage stamps (PAD) so that there is overlap among the blocks. Each postage stamp has a grid of output pixels, and each output pixel is constructed as a linear combination of input pixels. The algorithm allows any input pixels that are un-masked and within a specified radius (INPAD arc seconds) of the square stamp. The output pixels in a postage stamp are surrounded by a ring of transition pixels that allow the weights to vary smoothly when going to the next postage stamp.

The sizes of each hierarchical layer used in this paper are shown in Table 3. These are of course subject to refinement as we prepare for the eventual Roman science analysis.

We briefly review the postage stamp coaddition problem in Section 4.1, and indicate where there have been algorithm updates. We describe the block and mosaic construction in Section 4.2. The choice of output PSFs – essentially determining the resolution of the coadded images – is discussed in Section 4.3.

A diagram of the block coaddition script is shown in Fig. 5.

4.1 Postage stamp coaddition

The imcom formalism is developed in Rowe, Hirata & Rhodes (2011),⁷ and here we summarize only the key points that are needed here. imcom produces a coadded postage stamp starting from a set of input postage stamps: these images can be flattened and concatenated into a length n vector |$\lbrace I_i\rbrace _{i=0}^{n-1}$|⁠.⁸ Masked pixels can simply be removed from this vector with a corresponding decrease in n. These are to be combined into an output image with m pixels: |$H_\alpha = \sum _{i=0}^{n-1} T_{\alpha i}I_i$|⁠, where T is an m × n coaddition matrix. We denote the positions of each pixel centre by |$\lbrace {\boldsymbol r}_i\rbrace _{i=0}^{n-1}$| (input images) or |$\lbrace {\boldsymbol R}_\alpha \rbrace _{\alpha =0}^{m-1}$| (output image). If the input images have effective PSF G_i, then the output pixel α has an effective PSF at offset |${\boldsymbol s}$| given by⁹

It is a common difficulty with image combination algorithms that PSF_{out, α} varies from pixel to pixel in the coadded image. imcom takes as input a ‘target PSF’ Γ, and attempts to find coaddition weights T_αi that make PSF_{out, α} uniform and close to Γ. Quantitatively, we wish to minimize the leakage function,¹⁰

where the square norm is written in Fourier space with weighting |$\tilde{\Upsilon }({\boldsymbol u})$|⁠. The numerator U_α is the square norm of the difference between the as-realized output PSF and the target, while the denominator C is the square norm of the target (thus forming a dimensionless ratio). We will refer to U_α/C itself as the ‘leakage’ (smaller is better), and to the quantity −10log₁₀(U_α/C) as the ‘fidelity’ (larger is better: one can think of this as a suppression of PSF error in decibels). If the n × n noise covariance matrix of the input is N, then the m × m output covariance is |$\mathbf{\Sigma } = \mathbf{T}\mathbf{N}\mathbf{T}^{\rm {T}}$|⁠. imcom also tries to minimize the output covariance Σ_αα for input white noise (⁠|$\mathbf{N} = {\mathbb {I}}_{n\times n}$|⁠).¹¹

Having two objective functions, leakage U_α/C and noise Σ_αα, implies a trade-off: usually one can be improved at the expense of the other. We thus optimize row α of the T-matrix with respect to the objective function U_α + κ_αΣ_αα, where κ_α is a Lagrange multiplier that controls the trade-off (decreasing κ_α decreases leakage and increases noise; increasing κ_α increases leakage and decreases noise). The method of solving for T_αi given κ_α is described in Rowe, Hirata & Rhodes (2011); the algorithm performs a bisection search in log₁₀κ_α, with the user specifying the criterion for whether to increase or decrease κ_α.

Our postage stamp coaddition pipeline implements almost the same algorithm as in the original imcom, but with a python interface and using numpy routines for the linear algebra and Fast Fourier Transform (FFT) steps; the for loop over κ_α and the interpolation are coded in C and wrapped using the numpy C Application Programming Interface. Aside from changes in language, the substantive differences in the relative to the original imcom are

• The construction of T in Rowe, Hirata & Rhodes (2011) requires the construction of the correlations |$[G_i\circ G_j]({\boldsymbol s})$| for every pair of input PSFs. These are now computed and saved in the PSF_Overlap class, because these correlations only need to be re-computed on the scale of the PSF variation, not for every single postage stamp. This saves time when constructing a large-area mosaic.

• The interpolation of the correlations |$[G_i\circ G_j]({\boldsymbol s})$| in the original imcom was performed using bivariate polynomial interpolation. In our case, since the PSFs are band-limited, we know that their correlation is band-limited and given a grid size we know the oversampling factor (i.e. sample spacing relative to the Nyquist spacing). We have derived optimal interpolation kernels for this case in Appendix  A; we used the ‘D5, 5, |$\tfrac{1}{12}$|’ kernel here.

• Rather than simply specifying a target leakage U_α/C or a target noise Σ_αα, we specify a maximum noise (Σ_αα)_max as our first priority and a maximum leakage (U_α/C)_max as our second priority. That is, if both the leakage and noise targets can be met, the Lagrange multiplier search looks for the minimum noise that meets the leakage requirement; but if they cannot both be met, it treats the noise target as a hard constraint and finds the minimum leakage subject to that constraint. Of course in the latter case we would consider masking that pixel in a weak lensing analysis, but it avoids spectacular noise amplifications that could confuse some downstream analysis packages.

The postage stamp coaddition tools are in a separate GitHub repository (furry-parakeet) from the block driver.

4.2 Blocks and mosaics

The coaddition of a block (Fig. 5) begins with the selection of observation ID/SCA pairs that overlap with the block. Once we have these observations, a 4D numpy array is constructed to include all of the input data: in_data[i, j, y, x] contains the pixel in layer i, observation ID/SCA j, and pixel (x, y). Each layer consists of simulated images, noise fields, or grids of injected sources (Table 1). The input data are stored as 32-bit floating point numbers.

The algorithm then executes a loop over the postage stamps in that block. To save computing time, the PSFs are extracted and a PSF_Overlap object created every 4th postage stamp (it is treated as uniform in a block of 2 × 2 postage stamps, based on the PSF for the position at the common corner of the stamps). Then a suite of input postage stamps is created by mapping the centres of the output postage stamp back to each input image; these input stamps are deliberately oversized, since they can then be reduced by combining the input image mask with the logical test for whether a pixel is within the acceptance region (right panel of Fig. 4). Then we build the T-matrix (Section 4.1), and multiply to get an output cube (postage stamp for each layer).

A complication arises when the postage stamps are tiled to make a block: if the postage stamps are simply placed next to each other, there are noise discontinuities in the output image since one suddenly jumps to a different set of input pixels from which one can construct the output pixels. For some applications, these discontinuities may present no difficulties. However, some common applications such as peak finders that run at the resolution of the output image may be confused by these effects (e.g. Bertin & Arnouts 1996). Therefore, we use the transition pixels to smoothly transition from one postage stamp to the next. The overlap of the output postage stamps (including transition pixels) is 2k pixels; we define a sequence |$\lbrace a_m\rbrace _{m=1}^{2k}$|⁠. If we have a ‘left’ postage stamp I_left(x, y) (whose output region, excluding the transition pixels, is x < x_cut) and a ‘right’ postage stamp I_right(x, y) (whose output region, again excluding the transition pixels, is x ≥ x_cut), then in the transition region x_cut − k ≤ x < x_cut + k, we may make a merged image,

This approach does not affect the output PSF as long as the transition sequence satisfies a_m + a_{2k + 1 − m} = 1. We use the truncated sine function (Li & Mandelbaum 2023),

which avoids the discontinuities of the first derivative of the noise that would occur without the second term. The merging technique in equation (3) is trivially extendable to the y direction as well.

The mosaics are built using a single map projection with an array of blocks. In principle, one could select from any of the projections in the FITS WCS standard (Calabretta & Greisen 2002). We have used the stereographic projection here since it introduces no shear distortion, has zero plate scale gradient at the projection centre, and the Roman footprint can be efficiently tiled with square regions,¹² although we may revisit this choice in the future. At an angle θ from the projection centre, the stereographic projection introduces a lowest-order plate scale error of |$\frac{1}{4}\theta ^2$|⁠, or 2.4 × 10⁻⁵ at the corner of a 0.8 degree square (⁠|$\theta =0.8/\sqrt{2}\times \pi /180$|⁠). Thus for this test region we only need one mosaic.

We refer to a block within a mosaic by its array position: (0,0) is in the lower-left (southeast) corner, and |$({\tt BLOCK}-1,0)$| is in the lower-right (southwest corner).

4.3 Output PSFs

The imcom formalism is different from many other coaddition algorithms in that the target output PSF Γ is specified by the user and this is used to determine the input-to-output mapping (ssT-matrix), rather than the other way around. It is therefore advantageous to choose an output PSF that is circularly symmetric, so that the image coaddition and ‘rounding kernel’ are accomplished in a single step.

A straightforward choice is the convolution of an obstructed Airy disc (e.g. Rivolta 1986) with a Gaussian (whose width determines the resolution of the coadded image):

where ξ = λ/D is the diffraction scale at the central wavelength of that filter, J₁ is the Bessel function of the first kind, and ε is the linear obstruction fraction. This profile has an exact band limit at the diffraction scale, just like all of the input PSFs, and it has the usual ∝ 1/s³ diffraction wings with the same amplitude as for the input PSFs. Being circularly symmetric, it does not contain the diffraction spikes.

The choice of σ involves a trade-off. If σ is larger, then the output PSF is larger, which makes it more difficult to measure the shapes of small galaxies and increases blending. If σ is smaller, then the algorithm will attempt to reconstruct the higher spatial frequency Fourier modes in the image, which usually suffer more aliasing. This means that the output PSF leakage U_α/C is larger (worse). This effect is more significant with small numbers of dithers where not all combinations of input Fourier modes can be de-aliased. The implication is that for smaller σ (hence smaller target PSF), we may need to mask a larger portion of the final output image to control shear measurement systematics. (A possible way around this, which we have not explored in this paper, is to build several mosaics at different output resolutions, with the lowest resolution being available everywhere and the highest resolution available in select regions with more dithers.)

The smoothing scales chosen for each of the Roman filters in this simulation are shown in Table 4, and the modulation transfer functions are shown in Fig. 6. An example of the resolution versus output fidelity trade-off in the Roman Y106 band is shown in Fig. 7.

It may seem strange to additionally smooth the image to mitigate sampling effects; after all, the high angular resolution is one of the reasons to build a space mission in the first place. However, it is actually quite similar to the suggestion in Kannawadi, Rosenberg & Hoekstra (2021) to use a larger weighting scale for galaxy ellipticities¹³ and the convolution suggested by Shen et al. (2022). And we note that the output PSFs in Table 4 are still much smaller than those obtained in ground-based surveys.

4.4 Known issues

There are several issues that arose during the production runs for this project. These were not considered significant enough to warrant re-running the simulations, nor will they affect any of our primary results. However, we do plan to address them in a future version of the simulation:

• There was an indexing error that led to the wrong cosmic ray map being read in the Y106 and F184 simulations (this was fixed prior to the J129 and H158 runs). When a postage stamp was being created using the ith input image overlapping that postage stamp, the cosmic ray mask was read from the ith input image for the overall block. Since the cosmic ray maps are realizations of the same random process, this almost always has no effect. But it does mean that in a transition region between two postage stamps at the edge of a chip, there could be cases where the two postage stamps have inconsistent cosmic ray masks.

• The lookup algorithm that chooses input images to be used for a block used a search radius that did not account for plate scale variations. This means that there are a few instances where the corners of a chip overlap the block but that input image was not used. In the few cases where this happened, the sense of the effect is that the chip gaps in the coaddition simulation are larger than what we will have with the real data. In addition, sometimes the PSF computation points (at the centre of every 2 × 2 group of postage stamps) went off the edge of the SCA and hence that SCA was rejected from all 4 of the surrounding postage stamps; this also has the effect of increasing the effective chip gaps.

• The focal plane layout in the image simulations has some chip spacings that are different from the as-built Roman focal plane. The simulation and the WCS used to coadd the images are however internally consistent.

• The output WCS is shifted two postage stamps (2.5 arcsec) from the intended output region due to an error in writing the WCS. This has no effect on the validity of the results since all downstream steps consistently used this WCS; we simply coadded a slightly different region than intended.

• The weighting of Fourier modes in the leakage function was intended to be uniform, |$\tilde{\Upsilon }({\boldsymbol u})=1$| in the language of Rowe, Hirata & Rhodes (2011). An experimental feature to place more weight on the longer-wavelength modes, setting |$\tilde{\Upsilon }({\boldsymbol u}) = [1 + A_\Upsilon {\rm e}^{-2\pi ^2 {\sigma} _\Upsilon ^2(u^2+v^2)}]^2$|⁠, was accidentally left on for the production runs with |$A_\Upsilon =1$| and |$\sigma _\Upsilon = \frac{3}{2}\times 0.11$| arcsec. This placed more weight in the optimization than intended on the lowest spatial frequencies (relative to the highest spatial frequencies) by a factor of ≈4; since the resulting T is still a valid coaddition matrix, the results on imcom performance and applicability to weak lensing shape measurement remain valid.

5 RESULTS

We present some general results, examples, and statistics of the output PSF leakage U_α/C and features in the noise here. Statistical analyses of the results in the context of weak lensing are presented in the companion Paper II.

Some examples of regions in the coadded images are shown in Fig. 8.

5.1 Stars

An example of a star propagating through the coadd pipeline is shown in Fig. 9. We show both the imcom coaddition as well as the same star coadded in with the commonly used ‘Drizzle’ algorithm (Fruchter & Hook 2002) as described in appendix C of (Troxel et al. 2023). The Drizzle output has pixel scale 0.0575 arcsec and pixfrac = 0.7. Note that the imcom algorithm attempts to produce a round output PSF – this means, in particular, that the diffraction spikes are suppressed to the extent possible by adjusting the T-matrix. Some imperfections at the level of a few× 10⁻⁴ of the central surface brightness remain visible. Drizzle is a more local operation, and thus the asymmetry of the original PSFs as well as the diffraction spikes are more prominent.

5.2 Output PSF fidelity

We describe the output PSF fidelity in terms of the quantity −10log₁₀(U_α/C). This is a measure of how well the imcom algorithm thinks it has done on matching the PSF of the output image to the target PSF. The requested fidelity on this scale is 60, which corresponds to a 0.1 per cent error (in a root-sum-square sense) of all of the moments of the PSF.

The fidelity maps over the full simulation region for this paper are shown in Fig. 10. The chip gaps are easily visible as regions of lower fidelity: in these cases, there were insufficient dither positions to break the degeneracies among the various Fourier modes in the image. Equation (C18) gives a mode-counting argument for the minimum number of dithers required to disambiguate all Fourier modes in the astronomical scene (note that this is a necessary but not sufficient condition): this argument leads us to a number of dithers of five (Y106), four (J129), three (H158), and two (F184). However, these numbers could prove insufficient due to accidental degeneracies in the dithering pattern (e.g. two dithers at the same roll have an integer pixel offset) or masked pixels; and the required number may also be lower if the high spatial frequency modes have low enough transfer function that they can be ignored even if they are formally present in the image. The imcom simulations in this paper provide a more robust determination of the needed number of dithers in each filter. The worst performance is in Y106 band, since the Reference observing strategy was not originally designed for shape measurement using Y106¹⁴ and so we accepted fewer dither positions. The trouble comes mainly from regions with four or fewer dithers in Y106; given this result, and the importance of the Y106 filter in photometric redshifts (e.g. Hemmati et al. 2019), in a future version of the tiling strategy we might consider adding a dither position in Y106 even if we have to use shorter exposure times or accept one less dither in F184. The statistics of the PSF fidelity maps are shown in Fig. 11.

5.3 Output noise and Moiré patterns

Another aspect of working with multiple undersampled input images is the existence of Moiré patterns. If one combines two input images with no roll, but with a fractional difference ε in plate scale [so the plate scales are P and P(1 + ε)], then the images interlace to produce alternating regions where the samples land on top of each other, and where the samples are ideally interlaced (⁠|$\tfrac{1}{2}$| pixel offset), with a wavelength of P/ε. This behaviour is shown in Fig. 12 for the simple case of two input images in one dimension. One can see that in some parts of the Moiré pattern, the image is effectively sampled at the native pixel scale, whereas in other regions the sampling rate has been improved by a factor of 2; and an obvious concern for a survey is the resulting heterogeneity of the survey data. In particular, in the ‘degenerate’ regions, reconstructing a fully sampled image at output pixel positions in between the samples may prove impossible, or at least result in a large amplification of the noise in the input images (e.g. Lauer 1999a). For Roman, the plate scale is 0.11 arcsec, and the fractional variation in plate scale over a |$\sim \frac{1}{4}$| SCA dither in the y-direction is |$\sim 1.2{{\ \rm per\ cent}}$|⁠, leading to a predicted Moiré wavelength of ∼9 arcsec; this should be thought of as representative since the distortion gradient itself is spatially variable. An example of such a feature in the Roman image combination simulations is shown in Fig. 13. Such features are more common in regions with multiple intersecting chip gaps. While the Moiré patterns are intrinsic to the survey pattern – they represent spatial variations in what degeneracies are present in the data – the specific way they appear in higher-level products such as coadded images or ellipticity biases depends on the algorithm.

The full situation with the Moiré patterns is more complicated with N > 2 dithered images since then N(N − 1)/2 Moiré wavelengths are present and the pattern is not strictly periodic. In some regions, the dither positions will be closer to the degenerate case in Fig. 12; the fraction of the area where of all the dither positions land within, say, |$\frac{1}{4}$| pixels of each other can be expressed analytically,¹⁵ and is |$\frac{3}{16}=18.75{{\ \rm per\ cent}}$| for N = 3; |$\frac{1}{16}=6.25{{\ \rm per\ cent}}$| for N = 4; and |$\frac{5}{256}\approx 1.95{{\ \rm per\ cent}}$| for N = 5. If one goes up to 2 dimensions but with no rolls, then this argument applies in both x and y directions: for example, for N = 3, one would expect |$1-(1-\frac{3}{16})^2 \approx 34{{\ \rm per\ cent}}$| of the area to have samples clustered within |$\frac{1}{4}P$| in x, y, or both (this drops to 12 per cent for N = 4 and 4 per cent for N = 5). For the case with arbitrary rolls, we are unaware of any practical analytical results on these degeneracies, so numerical simulation remains the best approach to assessing the sampling impact of survey strategy choices for Roman. The histograms of the output noise properties from this simulation are shown in Fig. 14.

The banded features in Fig. 13 are visually striking and are a potential concern because coadded noise with an anisotropic power spectrum could introduce an additive noise bias (e.g. the centroid bias, Kaiser 2000; Bernstein & Jarvis 2002, although this is just one of a hierarchy of biases related to noise; Melchior & Viola 2012; Refregier et al. 2012). While the separation of the bands themselves is at the Moiré scale of order 10 arcsec, the band structures are coherent over scales of >1 arcmin (see the middle panel of Fig. 13) and one may wonder whether they introduce power at scales of interest for the weak lensing analysis. A quantitative investigation of the additive noise biases is carried out in the companion Paper II.

A fully zoomed-out output noise map is shown in Fig. 15.

5.4 Impact of charge diffusion

The ‘simple’ image simulations did not include a charge diffusion, and so we have not included charge diffusion in the main results of this paper. We intuitively expect that this is a ‘stress test’ of the application of imcom, since charge diffusion occurs before pixelization (Mosby et al. 2020) and thus a PSF including charge diffusion is better sampled than one without. However, it is important to check this expectation. We tested this by re-running one block in the most difficult filter (Y106) with several values of the charge diffusion, ranging from 0.0 to 0.4 pixels rms. As an intuitive guide to the importance of charge diffusion, the least-squares fit Gaussian to an Airy disc has width σ_Airy = 0.41λ/D, so by root-sum-square addition we would expect that inclusion of 0.24 pix rms charge diffusion effectively smears a diffraction-limited Y106 PSF to the resolution of J129, and 0.38 pix rms smears it to the resolution of H158; however, since the Airy disc and features induced by aberrations are not Gaussian, these should be viewed as rough estimates and the full output of the imcom should be used for planning purposes. We chose block location (44,34) because it includes some intersecting chip gaps (the right panel of Fig. 13 is in this block).

Block (44,34) contains a total of 1180 postage stamps (0.51 arcmin²) with a coverage of only three exposures. This generally consists of two exposures at one roll angle and one exposure at the other. For these postage stamps, we show the histogram of output PSF fidelity in Fig. 16. The fidelity improves with increasing charge diffusion. The measured charge diffusion for the Roman detectors is 0.27–0.35 pixels (Givans et al. 2022). We caution that those measurements were performed at an illumination wavelength of 500 nm, and it is possible that the charge diffusion length could be shorter at the longer wavelengths.

6 CONCLUSION

We have presented an implementation of the imcom image combination algorithm (Rowe, Hirata & Rhodes 2011) and its initial application to the Roman simulations of Troxel et al. (2023). The algorithm attempts to create a coadded image with a specified target PSF (generally taken to be circular and constant). We find a set of parameters and choice of target PSF where the algorithm is successful on the simulated Roman wide-area imaging survey data with the reference survey dithering pattern and the sampling in the four reference survey filters (Y106, J129, H158, and F184), including allowing for chip gaps and cosmetic defects. The quality of coadded PSF is reported in an output ‘fidelity’ map. The output noise maps show Moiré patterns characteristic of combining undersampled images with plate scale variations, especially in regions with small numbers of exposures. We have found that this effect is reduced if charge diffusion is included, which motivates more detailed characterization of the charge diffusion in Roman detectors before the final tiling strategy is selected. While the algorithm is computationally expensive, a mitigating factor is that simulated noise realizations and grids of injected sources can be processed at the same time at little additional cost.

Weak lensing applications of imcom coadds require not just ‘one-point’ statistics of the PSF and noise properties, but the correlation functions of the PSF residuals and the noise power spectrum and its spatial variation. These properties are investigated in the companion Paper II.

ACKNOWLEDGEMENTS

During the preparation of this work, CH and EM were supported by the National Aeronautics and Space Administration (NASA), contract 15-WFIRST15-0008, Simons Foundation grant 256298, and David & Lucile Packard Foundation grant 2021-72096. MY and MT were also supported by NASA award 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.’ Work at Argonne National Laboratory was supported under the U.S. Department of Energy (DOE) contract DE-AC02-06CH11357. This research used resources of the Argonne Leadership Computing Facility, which is supported by DOE Office of Science (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 Paul Martini, Peter Taylor, Chun-Hao To, 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 Center (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: INTERPOLATION KERNELS

The computation of the PSF overlap matrices for the image combination algorithm requires interpolation of a band-limited function from a 2D grid, |${\mathbb {Z}}^2$|⁠. A standard linear algorithm¹⁶ for interpolating a point based on (2K)² neighbours, separable in x and y, is

where α and γ are integers, 0 ≤ ξ_α, ξ_γ ≤ 1 are the fractional parts of the points to which we interpolate, and the functions w_μ are interpolation weights. We use |$\hat{f}$| to denote the interpolated function and f to denote the original (known only on the grid). The functions w_μ(ξ) are the interpolation weights and differ by interpolation scheme. Appendix  A is dedicated to the choice of weights. We summarize the major options in Table A1.

A1 Common choices

Choices common in the image processing literature include polynomial interpolation of order 2K − 1,

(this is a consequence of the polynomial interpolation formula of Waring 1779; see e.g. equation (A7) of Hirata et al. 2004, for this form) and Lanczos-K interpolation,

The Lanczos kernel does not preserve a constant background, i.e. it has ∑_μw_μ(ξ) ≠ 1, so one may also construct the background-conserving version

All of these interpolation schemes are exact in the limit that K → ∞ and the image is Nyquist-sampled (i.e. has Fourier wavevectors u with |$-\frac{1}{2}\lt u_x,u_y\lt \frac{1}{2}$|⁠). However for moderate values of K, they suffer from multiplicative errors as well as introduction of aliased Fourier modes (Bernstein & Gruen 2014). The fundamental issue is that these kernels are optimized for specific purposes: the polynomial kernel is best for functions with extremely low spatial frequencies, and the Lanczos kernel is a well-tested choice when moderate accuracy is required for functions that contain frequencies within a factor of ∼2 of the Nyquist limit (Turkowski 1990). Neither of these is especially well suited to the needs of weak lensing or other precision photometry programs, where one works with functions that are a few times Nyquist sampled (e.g. 4 × is common) but where four or more significant digits of accuracy in the interpolated function is required.

A2 Schemes optimized using a least-square error metric

To define an optimal interpolation scheme, we return to the 1D interpolation problem (since we have enforced separability) and define the various errors that may occur when a function is interpolated.¹⁷ We suppose that we are interpolating a complex exponential function of spatial frequency u, f_u(x) = e^2πiux (any function can be constructed as a superposition of these). Then there is an interpolated function |$\hat{f}_u(x)$|⁠. The original spatial frequency is present with some multiplicative error m(u),

There is also a leakage into the other modes (or ‘ghosting’ in the terminology of Bernstein & Gruen 2014),

The total mean squared error in the reconstruction is the sum of these:

That is, ε(u) is the root-sum-square of the multiplicative bias and all of the ghost amplitudes. We are now ready to define some conditions for an optimized kernel. We have tested several choices; these are summarized in Table A1.

If we are working with a band-limited function with spatial frequencies |u| < R, then we might choose to minimize

where ρ is a non-negative even window function. Setting the functional derivative of Ω to zero gives

This is a linear equation for w_ν(ξ). It has a solution of the form

where the system matrix is

and

Since C is even, S is symmetric.

The most obvious choice is to weight all Fourier modes in the band limit equally: ρ(u) = 1. This leads to what we call the ‘square window least-squares error (LSE)’ interpolation method (‘square’ because the weighting of different input Fourier modes u is constant for |u| < R and 0 otherwise),

where we use the numpy sinc convention, |${\rm sinc\, }z=\sin (\pi z)/(\pi z)$|⁠. Since the matrix S⁻¹ does not depend on ξ, but only on the parameters K and R, it can be computed once for an interpolation scheme and saved.

In some cases, we might want some control over the errors for input Fourier modes out to the band limit R, but we want the tightest control for small u. This is the case, for example, of a diffraction limited image whose power spectrum only truly drops to zero at u = R but is small at u ≳ R/2. We might then modify the objective function to be ρ(u) = 1 − |u|/R. This leads to

This leads to the ‘triangle window LSE.’

A2.1 Background conservation

The aforementioned interpolation algorithms do not conserve a constant background, i.e. |$\sum _{ {\mu} =1-K}^K w_ {\mu} (\xi)$| is not exactly equal to unity. This is not surprising since we optimized the weights so that the integral of the LSE over a range of frequencies is minimized. For interpolation of images with a large background, one might want the constant mode to be represented perfectly. Thus we are led to considering the lowest cost (Ω) 2K-point interpolation scheme that has zero error for the constant mode, i.e. the u = 0 mode.

This problem can be solved by adding a formally infinite δ-function to ρ(u), i.e. we construct a new objective function,

where e is a length 2K column vector of all ones, and taking the limit of λ → ∞. This leads to a new weight that can be obtained from the Sherman–Morrison formula, in the form of Bartlett (1951),

Now at this stage, we see that S⁻¹b is simply w, the weight vector that is obtained without enforcing background conservation. Then we have a new weight vector given by

where

The vector |${\boldsymbol \eta }$| depends only on the choice of interpolation scheme and can be pre-tabulated. Note that by construction, the entries in |${\boldsymbol \eta }$| sum to unity, guaranteeing that |$\sum _{ {\mu} =1-K}^K w^{\prime }_ {\mu} =1$|⁠. The use of any vector |${\boldsymbol \eta }$| that sums to unity would fix the background conservation issue (for example, the rescaling in Bernstein & Gruen 2014 corresponds to w_μ/∑_σw_σ), but equation (A18) does so at the least cost to the weighted mean square error over the specified range in u.

Interpolation kernels that are ‘improved’ using this approach are denoted with a ′ (e.g. S′ or T′).

A2.2 Rounding error issues

A disadvantage of the LSE methods for large K and small R is that the S matrix becomes nearly singular, resulting in amplification of roundoff error when one takes the inverse; eventually, these errors can exceed the formal analytical error of the interpolation method. Tools such as direct system solution without the inverse (e.g. numpy.linalg.solve) reduce the problem but do not eliminate it. The cases we have shown in the figures do not suffer from these problems when floating point computations are performed using IEEE 754 standard 64-bit arithmetic. But for the most demanding applications (where we aim for ε ≲ 10⁻⁸), we aim for alternatives to direct computation of S.

A3 Discretized error metric

We recall that an integral is often evaluated by taking a linear combination of values of the functions at specified abscissae. Applying this idea to equation (A8), we may define an error metric that is not an integral over Fourier modes, but rather a discrete sum,

where |$\lbrace z_j\rbrace _{j=1}^L$| are the abscissae (rescaled so that 0 < z_j < 1); |$\lbrace \alpha _j\rbrace _{j=1}^L$| are the weights, assumed to be positive; and we have enforced symmetry of positive and negative frequencies. We require L ≥ K so that the matrix S has full rank. Then

where we defined ζ_j = 2πRz_j. The advantage of this is that the matrix S can be factorized as S = M^TM, where M is the 2L × 2K matrix

and we have defined |$\beta _j \equiv \sqrt{2\alpha _j}$| and |$s_j \equiv -K-\frac{1}{2}+j$|⁠. The length 2K vector it b can be expressed as it b = M^Tc where

(These results can be verified by direct multiplication, as well as the identity cos (x − y) = cos xcos y + sin xsin y.) Next we perform the singular value decomposition of M, so that |$\mathbf{M}=\mathbf{UDV}$|⁠, where U is a 2L × 2L orthogonal matrix; D is a 2L × 2K matrix with all zeros except for the diagonal |$\lbrace D_{jj}\rbrace _{j=1}^{2K}$|⁠; and V is a 2K × 2K orthogonal matrix. In this case, equation (A10) simplifies to

This can be expressed as

where the 2K × 2L matrix H can be written as

where in the prefactor we define β_{l + L} = β_l. Note that H depends on K, L, and the weights z_j and α_j, but not on ξ; thus once a scheme is chosen, it can be pre-computed. The expression contains a division by a singular value D_kk, but since S = M^TM, the condition ratio of M can be better than that of S (if L = K, it it the square root). This means that roundoff errors are correspondingly less pernicious.

The choice remains of the weights. The simplest choice, which we adopt here, is to use z_j and α_j corresponding to 2L-point Gauss-Legendre quadrature; we expect it to be the best discrete approximation to the ‘square wave LSE’ scheme. Furthermore, while the method works for L > K, we have not found significant gains that offset the increased computational cost (we must compute 2L trigonometric functions for each 1D interpolation). Thus we have retained L = K.

The discrete window interpolators with L = K are in principle exact at the spatial frequencies u_j, since then M is invertible so S⁻¹b = M⁻¹c and one can show that an input function f(x) that is a cosine or sine wave leads to a vector |${\boldsymbol f}^{\rm T} = (f(1-K),\, f(2-K)\, ...f(K))$| that is a constant times a row of M. Then the interpolated value, |$\hat{f}(\xi)=\mathit {f}^{\rm {T}}\mathbf{M}^{-1}{\mathit {c}} = (\mathbf{M}^{\rm {T}-1}{\mathit {f}})\mathit {c}$|⁠, collapses. However, in finite-precision arithmetic, this property may not hold.

A4 Examples

We now present some specific choices of kernel that are optimized for various choices of oversampling factor and precision. Recall again that in what follows, ‘sinc’ follows the numpy convention, |${\rm sinc\, }z = \sin (\pi z)/(\pi z)$|⁠. Also the oversampling relative to Nyquist is 1/(2R): thus |$R=\frac{1}{4}$| is 2 × Nyquist, |$R=\frac{1}{8}$| is 4 × Nyquist, etc.

A4.1 Six-point interpolation

A well-studied choice in weak lensing image processing is 4 × Nyquist sampling (e.g. by zero-padding an FFT) and six-point (K = 3) interpolation. Bernstein & Gruen (2014) found that multiplicative errors and ghosting could be reduced to the |$\lt 0.1{{\ \rm per\ cent}}$| level with quintic polynomial interpolation (P3). The square window LSE method gives the interpolation weights for the ‘S3, |$\tfrac{1}{8}$|’ scheme,

where the |$q^{\rm (S)}_{ {\mu} {\sigma} }$| coefficients are given in Table A2. Similarly, for the triangle window LSE method or ‘T3, |$\tfrac{1}{8}$|’ scheme, we get,

The |${\boldsymbol \eta }$| vectors are also shown in the table, so one may use equation (A17) to construct the background-conserving kernels |$w_{ {\mu} }(\xi |{\rm S}^{\prime }3,\tfrac{1}{8})$| and |$w_{ {\mu} }(\xi |{\rm T}^{\prime }3,\tfrac{1}{8})$|⁠.

A comparison of the errors for the various six-point interpolation schemes is shown in the upper panel of Fig. A1. The quintic interpolation scheme performs best at very low spatial frequencies (⁠|$u\lt \tfrac{1}{16}$|⁠). The Lanczos-3 kernel has large errors, several tenths of a per cent, across all of the low spatial frequencies, but with this sacrifice it gains improved behaviour at intermediate spatial frequencies (u ∼ 0.3). The ‘S3, |$\tfrac{1}{8}$|’ kernel is best at showing small (<10⁻⁴) errors over the widest range. In particular, all the modes in a 4 × Nyquist sampled image (⁠|$u\lt \tfrac{1}{8}$|⁠) are reconstructed with errors <10⁻⁴, whereas the quintic has an error of |$\varepsilon (\tfrac{1}{8})= 7.4\times 10^{-4}$|⁠. The ‘T3, |$\tfrac{1}{8}$|’ kernel behaves similarly, although as expected the lowest spatial frequencies behave better (error <10⁻⁵) at the expense of slightly worse performance right at the band limit (reaching error of 1.4 × 10⁻⁴ at |$u=\tfrac{1}{8}$|⁠).

A4.2 10-point interpolation

For very high precision applications – for example, if interpolation is to be used as a step in constructing the A matrix in Rowe, Hirata & Rhodes (2011), where small errors can be magnified in subsequent linear algebra steps – even the ‘S3, |$\tfrac{1}{8}$|’ and ‘T3, |$\tfrac{1}{8}$|’ schemes may not be sufficient. We have therefore implemented the 10-point ‘D5, 5, |$\tfrac{1}{12}$|’ scheme, which is to be applied to images that are at least 6 × Nyquist sampled.

The coefficients w_μ(ξ) are a sum of 10 trigonometric terms, as given by equation (A24). The coefficients are given in Table A3. The performance is shown in the lower panel of Fig. A1; note that at |$u\lt \tfrac{1}{12}$|⁠, the error is very small, ε(u) < 1.5 × 10⁻⁹. The nonic polynomial interpolator (P5), by contrast, ‘only’ achieves |$\varepsilon (\tfrac{1}{12}) \approx 2.4\times 10^{-7}$|⁠; but at the very smallest values of u (⁠|$\lt \tfrac{1}{24}$|⁠) it has the best performance. The Lanczos interpolator (L5), as expected, has significant errors (>10⁻³) at low frequencies, but has good performance (<1.63 × 10⁻²) all the way out to u = 0.38 (i.e. modes that are only 1.3 × Nyquist sampled).

The ‘D5, 5, |$\tfrac{1}{12}$|’ interpolator should be ‘exact’ at the five frequencies u_j that are |$\frac{1}{12}$| times the Gauss-Legendre abscissae, i.e. the roots of the P₁₀ Legendre polynomial. One can see these as ‘dips’ in the bottom panel of Fig. A1. However, in IEEE 754 64-bit (double precision) arithmetic, one never reaches exactly zero; the ε(u) curves have minima at ∼1.1 × 10⁻¹⁰ (inspection of the data file shows these are actual minima, not simply the closest point to u_j that was plotted). This is unsurprising since the amplitudes in Table A3 reach ∼10⁶ and given the limitations of double precision representations of numbers (53 bits of precision: 2⁻⁵³ ∼ 10⁻¹⁶). Further improvement would probably require the construction of alternate basis sets that are more orthogonal than cosines and sines; this presents no fundamental challenge but we suspect on many architectures the computation time would be increased. Since errors of ‘a few × 10⁻¹⁰’ are well below our tolerance, and are likely good enough for the vast majority of image processing needs in astronomy, we stop here and recommend 6 × Nyquist sampling and the ‘D5, 5 , |$\tfrac{1}{12}$|’ interpolation scheme for the most demanding applications.

APPENDIX B: COMPUTING RESOURCES

Most of the computing for the image combination runs was carried out at the Pitzer cluster (Ohio Supercomputer Center 2018) at the Ohio Supercomputer Center (OSC), and the Duke Compute Cluster (DCC). A small number of test runs were carried out on the Amazon Web Services (AWS) Elastic Compute Cloud (EC2), or on a MacBook Pro laptop (‘Mac’). The statistics for the computing costs are in Table B1. Note that several experiments were carried out near the beginning of the runs, including running the code on multiple cores using numpy multithreading; this turned out not to be efficient for the current code structure since only some operations benefited from the multithreading, and we dropped it for the remainder of this project. For most of the runs, memory (≈10 GB) rather than cores was the limiting factor; therefore the rate at which the runs could be completed (and, depending on platform, billing rates) were set by the memory usage rather than core-hours.

The central processing units (CPUs) used were the Intel Xeon Gold 6148 (2.4 GHz) or Platinum 8268 (2.9 GHz) on OSC; the Intel Xeon E5-2680 v3 (2.5 GHz) on DCC; the Intel Xeon Platinum 8124M (3.0 GHz) on AWS; and the Apple M1 Pro18,3 on the Mac laptop.

A change partway through the J129-band runs was to introduce memory mapping (via the numpy.memmap function) for the input data cube to reduce the memory usage with many input layers. The input data cube is used least frequently, and on some nodes memory rather than core count limited the rate of processing. We found that using numpy.memmap for the input data typically reduced the memory usage by ∼30 per cent.

Note that the number of operations can vary by band: J129 has the largest number of dither positions in the Reference Survey, which means that the matrices used to construct T are larger. Thus one expects J129 to require more core-hours than the other bands, even for the same code and computing architecture.

APPENDIX C: THE EFFECT OF PIXEL-LEVEL OPERATIONS ON UNDERSAMPLED IMAGES

In Appendix  C, we review the mathematical properties of undersampled data, and how these properties affect the operations commonly used to calibrate shear estimators. Sampling is parameterized by Q, which is the number of pixels corresponding to one cycle at the maximum spatial frequency; Q < 2 corresponds to undersampling by the usual Nyquist crtierion, and smaller Q represents more severe undersampling. We show that the usual implementation of Metacalibration (Huff & Mandelbaum 2017; Sheldon & Huff 2017) suffers from aliasing when acting on undersampled images. If the undersampling is not too severe (‘weak undersampling’ or 1 < Q < 2), then Metacalibration can be implemented by re-convolving the image with a PSF that excludes the problematic Fourier modes. However, the Roman Y106 band and the NIR bands on Euclid have ‘strong undersampling’ (Q < 1), so the mathematical framework of standard Metacalibration does not apply to these cases. The alternative Deep Metacalibration algorithm (Zhang, Sheldon & Becker 2023) in principle could address these problems when Q < 1, but many Fourier modes criss-crossing the (u, v)-plane must be excluded from the re-convolution PSF. The analytic shear response method of Zhang, Sheldon & Becker (2023) faces similar issues: in the weakly undersampled case, convolution in pre-processing can eliminate the aliased Fourier modes, but in the strong undersampling case the formalism as currently implemented cannot be applied.

We leave open the possibility that some new mathematical framework might be developed in the future to deal with these issues in the general undersampled case. But given these results, we are motivated to focus our attention on image processing algorithms that can recover a fully sampled image.

To avoid clutter, we work in units of the input pixel scale (s_in = 1).

C1 Fourier transforms and operators

We perform continuous Fourier transforms in the ‘waves per pixel’ convention,

Convolution in one domain (real or Fourier) is equivalent to multiplication in the other; we denote the convolution of two functions by * and multiplication by ·.

Shear measurement and shear calibration algorithms often use the finite translation operator |${\mathcal {T}}_{\boldsymbol{s}}$| with displacement |${\boldsymbol{s}}$|⁠,

and the distortion operator |${\mathcal {Y}}_{\kappa ,{\boldsymbol{\gamma }},\varphi }$| with magnification κ, shear |${\boldsymbol{\gamma }}$|⁠, and rotation φ,

where

(This is the description of shear in terms of the polar decomposition of the 2 × 2 matrix |$\mathbf{M}_{\kappa ,{\mathbf{\gamma }},\varphi }$|⁠. The polar decomposition is convenient because it allows us to rotate a galaxy image by angle φ, then shear and magnify it, and then insert it into an image. It is thus convenient for problems where we want to know what a galaxy looks like at a random orientation φ, or if we want to insert both a galaxy and its 90°-rotated image to cancel shape noise in a simulation.) These operate on Fourier transforms in accordance with

and

When we discuss the magnitudes of the wave vectors, it is helpful to keep in mind that in the weak lensing regime |$|\kappa |+|{\boldsymbol{\gamma }}|\lt 1$|⁠, the maximum singular value of |$\mathbf{M}_{\kappa ,{\mathbf{\gamma }},\varphi }^{-1\, \rm T}$| is

that is, for (u, v) in a circle of radius r, the farthest that |$\mathbf{M}_{\kappa ,{\mathbf{\gamma }}}^{-1\, \rm T}(u,v)$| can be from the origin is Λr.

A discretely sampled function f_{x, y}, measured only at integer pixel positions |$(x,y)\in {\mathbb {Z}}^2$|⁠, instead has a Fourier transform,

where |${\mathcal {B}} = \lbrace (u,v): -\frac{1}{2}\le u\lt \frac{1}{2}, -\frac{1}{2}\le v\lt \frac{1}{2} \rbrace$| is the first Brillouin zone or the unit square centred on the origin (shaded square in Fig. C1). The two types of Fourier transforms are related to each other via the III (‘Shah’ or Dirac comb) function,

which is its own Fourier transform (⁠|$\widetilde{\rm III}={\rm III}$|⁠). If we discretely sample a function f at the pixel positions, f_{x, y} = f(x, y) for integer x, y, then

One sees that the Fourier transform of the discretely sampled data contains a term associated with the original field (Δu = Δv = 0), as well as other terms offset by an integer numbers of cycles per pixel.

We define the fractional part operator

that maps u to the range |$-\frac{1}{2}$| to |$+\frac{1}{2}$|⁠. Then |$({\mathbb {F}}u,{\mathbb {F}}v)$| is the single mode in the first Brillouin zone that aliases to (u, v).

C2 Astronomical scenes and available information

Let us define an image which is a convolution of intrinsic projected scene S and the effective PSF G (which includes the pixel tophat), which is then sampled at integer positions,

The Fourier transform of the discretely sampled data is then related to the Fourier transform of the sky scene by

While in principle the sum in equation (C13) is infinite, in practice G has a band limit or maximum spatial frequency present: |$\tilde{G}(u,v)=0$| if |$\sqrt{u^2+v^2}\ge 1/Q$|⁠, where Q is a sampling parameter (equal to the number of pixels across a single cycle at the maximum spatial frequency). For a space-based observatory, the optics usually operate at the diffraction limit: thus Q = λ/(Ds_in), where λ is the wavelength of light, D is the diameter of the entrance pupil, and s_in is the pixel scale. For a ground-based observatory, atmospheric smearing usually suppresses the highest spatial frequencies. While turbulent contributions to |$\tilde{G}$| usually decline smoothly rather than going to zero at a ‘hard’ cutoff in spatial frequency, there is still some radius in the (u, v)-plane beyond which |$\tilde{G}(u,v)$| is negligible.

The sampling parameter Q restricts the number |${\mathcal {N}}(u,v)$| of non-zero terms in equation (C13), with the number of terms growing as Q decreases (more undersampled). To determine |${\mathcal {N}}(u,v)$|⁠, we draw a square lattice with unit spacing centred at (u, v), and count the number of lattice points in a circle of radius 1/Q centred at the origin (Fig. C1). The results are visualized in Fig. C2. Each panel shows the region |${\mathcal {B}}$| in (u, v)-space covered by a discrete Fourier transform centred on (0,0) and this goes from |$-\frac{1}{2}$| to |$\frac{1}{2}$| cycles per pixel on both axes. The colour scale shows the number of terms that contribute in equation (C13), i.e. the number of Fourier modes in the original image that alias to (u, v) in that region. The top left panel shows the case where the image is Nyquist-sampled and each Fourier mode in the image is unique (the maximum allowed frequency can be seen as the radius of the circle |$\sqrt{u^2 + v^2}$|⁠). It can be seen that more Fourier modes are aliased as the sampling factor Q decreases (i.e. increasing the radius 1/Q of the circle). Intuitively, as the circle begins to overflow outside the region, neighbouring circles centred at |$(u, v) \in \mathbb {Z}^2$| begin to overlap and the frequency in the overlapped region is aliased.

If |${\mathcal {N}}(u,v)=1$|⁠, then only (Δu, Δv) = (0, 0) contributes, which means that equation (C13) becomes |$\check{f}(u, v) = \tilde{G}(u, v) \tilde{S}(u, v)$|⁠, and |$\tilde{S}(u,v)$| can be reconstructed by division in Fourier space. However, when |${\mathcal {N}}(u,v)\ge 2$|⁠, there are more than one |$\tilde{S}(u+\Delta u, v+\Delta v)$| to solve but only one constraint |$\check{f}(u, v)$|⁠, so |$\tilde{S}(u,v)$| is not recoverable. Thus the region that can be reconstructed unambiguously is the region |${\mathcal {N}}^{-1}(1) = \lbrace (u,v): {\mathcal {N}}(u,v)=1\rbrace$| where exactly one mode contributes. If Q > 2 (case a in Fig. C2), then we say that the image is oversampled and |${\mathcal {N}}^{-1}(1)$| is the disc |${\mathcal {D}}_{1/Q} = \lbrace (u,v): \sqrt{u^2+v^2}\lt 1/Q\rbrace$| of radius 1/Q. If 1 < Q < 2 (cases b and c), then the image is weakly undersampled: while some modes are aliased, others are not and |${\mathcal {N}}^{-1}(1)$| is non-empty; the largest disc contained in it is |${\mathcal {D}}_{1-1/Q}$|⁠. If Q < 1 (cases d and beyond), then the image is strongly undersampled and |${\mathcal {N}}^{-1}(1)$| is the null set: there are no Fourier modes |$\tilde{S}(u,v)$| that can be recovered.

C3 Recovering information with translational dithers

Lauer (1999a) discusses the linear solution to de-alias these modes for the case where multiple dithered images are taken with no rolls: if there are N_exp images with dither offsets |${\boldsymbol{s}}_1 ... {\boldsymbol{s}}_{N_{\rm exp}}$|⁠, then the jth image can be thought of as measuring the displaced scene |${\mathcal {T}}_{{\boldsymbol{s}}_j}S$|⁠. In this case, the discrete Fourier transform of the jth image is

(compare to equation C13). This is a system of N_exp equations with |${\mathcal {N}}(u,v)$| unknowns, and it can be solved if the |$N_{\rm exp}\times {\mathcal {N}}(u,v)$| matrix L with coefficients

where j represents one of the N_exp images and m represents one of the (Δu, Δv) terms in the sum, has rank |${\mathcal {N}}(u,v)$|⁠. Note that it is necessary but not sufficient that |$N_{\rm exp}\ge {\mathcal {N}}(u,v)$|⁠.

Particular dithering patterns have been proposed for the various sampling cases, For example, a diagonal dither (fig. 2.2 of Gonzaga et al. 2012) with N_exp = 2 positions at |${\boldsymbol{s}} = (0,0), (\frac{1}{2},\frac{1}{2})$| works for case (b), because in the |${\mathcal {N}}(u,v)=2$| regions (⁠|${\mathcal {N}}=2$| in Fig. C2b) the lattice points contributing to the sum are (Δu, Δv) = (0, 0) and one of (1,0), (−1, 0), (0,1), or (0, −1). Then the matrix L is

which has rank 2. Similarly, a ‘2 × 2’ dither with N_exp = 4 positions (fig. 2.3 of Gonzaga et al. 2012) at |${\boldsymbol{s}} = (0,0), (\frac{1}{2},0), (0,\frac{1}{2}), (\frac{1}{2},\frac{1}{2})$| works for case (c). For example, in the lower-left corner region of Fig. C2(c), we have |${\mathcal {N}}(u,v)=4$| and the lattice points contributing to the sum are (Δu, Δv) = (0, 0), (1,0), (0,1), and (1,1). Then the matrix L is

which has rank 4.

Cases with more small dither steps for other sampling cases could be discussed: for example, the N_exp = 5 ‘knight’s move’ pattern |${\boldsymbol{s}}_j = (\frac{2}{5}j,\frac{1}{5}j)$| works for case (d), and Gonzaga et al. (2012) discuss the N_exp = 8 point dither (for cases e–g) and N_exp = 9 point dither (for case h). [The N_exp = 6 pattern considered in fig. 2.4 of Gonzaga et al. (2012) does not lead to an L matrix of full rank for cases (e) and (f).] However such large numbers of small deterministic dithers – which would have to be supplemented by large dithers over chip gaps – do not fit naturally into a fast wide-angle survey and thus were not options for the Roman HLIS.

This approach does not apply to cases with rolls, since in an exposure rolled by angle α there are aliased Fourier modes with (Δu, Δv) = (cos α, sin α), and the ‘network’ of coupled modes (equation C14) does not close with a finite number of constraints and unknowns. The one straightforward analytic result in this case comes from mode counting: the number of ‘measurements’ (input pixels) is |$N_{\rm exp}{\cal A}$| where |${\cal A}$| is the survey area (again in units where s_in = 1), whereas the number of ‘unknowns’ (output Fourier modes) is |$\pi {\cal A}/Q^2$| (recall the number of modes is area in real space times area in Fourier space). Therefore to disambiguate all the modes, it is necessary (but not sufficient) that

The imcom formalism (Rowe, Hirata & Rhodes 2011) was developed to handle this additional case numerically, which is needed for the proposed Roman strategy since it includes multiple roll angles.

In general, we will denote by |${\mathcal {E}}$| the region where |$\tilde{S}(u,v)$| can be reconstructed, with the understanding that for N_exp = 1 dither, |${\mathcal {E}} = {\mathcal {N}}^{-1}(1)$|⁠.

C4 The operations in shear calibration techniques

We now recall the operations used in the various shear calibration techniques. There is a class of techniques that works directly with Fourier-domain quantities (Bernstein & Armstrong 2014; Bernstein et al. 2016), and in this case one wants to work in a region of the Fourier plane that is contained within |${\mathcal {E}}$|⁠, and actually contains a ‘buffer’ region before one reaches the edge of |${\mathcal {E}}$|⁠; see e.g. Bernstein (2010, section 4). We will see that this behaviour is generic.

We start with the ‘standard’ Metacalibration (Huff & Mandelbaum 2017; Sheldon & Huff 2017), and then investigate ‘deep’ Metacalibration (Zhang, Sheldon & Becker 2023), with a particular emphasis on the operations that are potentially problematic on undersampled data (e.g. Hoekstra, Kannawadi & Kitching 2021). Finally, we consider the use of pixel responses as in Li & Mandelbaum (2023).

C4.1 Standard Metacalibration

The basic operation in standard Metacalibration is to take an image; infer the sky scene S (e.g. by de-convolving the PSF); apply a magnification κ and shear |${\boldsymbol{\gamma }}$|⁠; and re-convolve it with a new PSF G_r. The resulting image is

so that

We suppose that the re-convolution PSF has a Fourier transform |$\tilde{G}_{\rm r}(u+\Delta u,v+\Delta v)$| with support in some circular region out to radius R_r.

For non-zero shear, the Fourier modes |$\mathbf{M}_{0,{\mathbf{\gamma }},0}^{-1\, \rm T}(u+\Delta u,v+\Delta v)$| do not alias to each other. In this case, we requires knowledge of |$\tilde{S}$| at all the points that contribute to the sum. That means that for every point (u′, v′) in the support of |$\tilde{G}_{\rm r}$|⁠, |$\mathbf{M}_{0,{\mathbf{\gamma }},0}^{-1\, \rm T}(u^{\prime },v^{\prime })$| must be in the reconstructible region |${\mathcal {E}}$|⁠. This leads to the following cases:

• For oversampled data (Q > 2), |${\mathcal {E}} = {\mathcal {D}}_{1/Q}$|⁠. This means that the radius of support of |$\tilde{G}_{\rm r}$| must have |$R_{\rm r}\le (1-|{\boldsymbol{\gamma }}|)/Q$| (see equation C7). Since the original PSF |$\tilde{G}$| has support out to a radius 1/Q, this encapsulates the usual notion that the re-convolution PSF must be at least a little bit bigger than the original PSF, with the definition of ‘a little bit’ being determined by the applied shear |${\boldsymbol{\gamma }}$|⁠.

• For weakly undersampled data (1 < Q < 2) without dithering, the largest disc contained within |${\mathcal {E}}$| has radius 1 − Q⁻¹. Thus we can carry out the standard Metacalibration operation if the re-convolution PSF satisfies |$R_{\rm r}\le (1-|{\boldsymbol{\gamma }}|)(1-Q^{-1})$|⁠. Note in this case that as the native PSF gets smaller at fixed pixel scale (Q gets smaller), the Fourier-space cutoff R_r must get smaller and hence the re-convolution PSF has to get bigger in real space.

• For strongly undersampled data (Q < 1) without dithering, |${\mathcal {E}}$| is the null set and standard Metacalibration cannot be implemented.

• For data with enough dithers to de-alias all of the Fourier modes, the radius R_r is again limited by |$(1-|{\boldsymbol{\gamma }}|)/Q$|⁠.

For single-epoch Roman images, ‘standard’ Metacalibration does not apply in Y106. It theoretically can be applied in J129 (Q = 1.021); however, 1 − Q⁻¹ is so small that it would not be practical (the re-convolution PSF would have to be more compact in Fourier space, and thus larger in real space, than even a typical ground-based PSF). The algorithm could be of interest in the redder filters.

C4.2 Deep Metacalibration

We now turn to ‘deep’ Metacalibration (Zhang, Sheldon & Becker 2023). Here there is both a deep scene S_deep and a wide scene S_wide that have to be rendered with a reconvolution PSF G_r. The shear operation, equation (C20), is applied only to the deep data, whereas no shear is applied to the wide data. In the context of both ground- and space-based surveys, deep fields usually have a much larger number of observations than the wide survey; but for a diffraction-limited space telescope the band limit and sampling parameter Q are the same in the deep and wide surveys. (This may or may not be the case on the ground.) The deep survey recovers information over some region of the (u, v)-plane |${\mathcal {E}}_{\rm deep}$|⁠. Because of the large number of dithers, this region may be larger than |${\mathcal {E}}_{\rm wide}$|⁠, but it is still bounded by the radius 1/Q.

In order to recover all of the modes in the deep data that we need for equation (C20), we must have

However, if we are working with an undersampled single-epoch image in the wide survey (with PSF G_wide), we do not have access to all Fourier modes independently, but only the linear combinations,

This means that for a given (u, v), if there is a |$(\Delta u,\Delta v)\in {\mathbb {Z}}^2$| satisfying

then the linear combination appearing in equation (C22) must simply be zero, and |$\tilde{G}_{\rm r}$| must be zero for all modes that alias to (u, v). And in particular, if there is a lattice point (u + Δu, v + Δv) in the troublesome annulus,

then there is an orientation for the shear such that equation (C23) is satisfied, and that Fourier mode must not be present in the reconvolved image. So if we were to implement Deep Metacalibration on an undersampled image, we must avoid all modes that alias to this annulus, i.e. set |$\tilde{G}_{\rm r}(u,v)=0$| in these regions. In practice, since one uses Deep Metacalibration with small values of |$|{\boldsymbol{\gamma }}|$|⁠, the troublesome regions of (u, v)-space are those near the circle |$\sqrt{(u+\Delta u)^2 + (v+\Delta v)^2} = 1/Q$|⁠, i.e. the ‘edges’ seen in Fig. C2.

The aforementioned procedure is in principle possible, and the conditions required for it to work are slightly less restrictive than for standard Metacalibration – indeed, these conditions could be met even in the strongly undersampled case. However, in practice the re-convolution PSF must be very large – for example, if Q = 0.834 (Y106 band) and |$|{\boldsymbol{\gamma }}|=0.01$|⁠, then |$\tilde{G}_{\rm r}$| must go to zero at a radius of 0.187 cycles per pixel, which is equivalent to an Airy disc with a full width at half maximum of 0.61 arcsec. This means that the re-convolution PSF effectively degrades the sky image to ground-based resolution, so even though the mathematics works this is not an attractive option for analyzing Roman data.

C4.3 Pixel basis function technique

Li & Mandelbaum (2023) present another approach to shear calibration, using many of the same ideas as Metacalibration but using analytic differentiation to compute the shear response. The idea is that the flux in a given pixel |$(x,y)\in {\mathbb {Z}}^2$| can be written in terms of the sky scene S via

where

is the pixel basis function. [Li & Mandelbaum (2023) define ϕ_{x, y}(u, v), which is our Φ_{x, y}(u, v)/(2π)²; the difference results from the choice of cycles per pixel versus radians per pixel as the unit of spatial frequency.] One may then consider a partial derivative of the pixel with respect to a shear component (a = 1, 2),

where

and

are the shear responses of the pixel basis functions. (These are equation 19 of Li & Mandelbaum 2023, but written for a non-Gaussian PSF.)

We are interested in the conditions under which f_{x, y; a} can be determined from the data, i.e. when the shear response functions can be written in terms of a linear combination of the pixel response functions themselves,

We write the Fourier transform of the weights on the (Δx, Δy) indices,

Substituting this and the expressions for Φ_{x, y} and Φ_{x, y; a} into equation (C30), we find that in the summation over Δx and Δy, only the cases where (u, v) and (u′, v′) alias to each other survive. Then if we cancel the common factors of e^{−2πi(ux + vy)} in equation (C30), we arrive at

and

For oversampled data, this is only non-trivial for modes (u, v) that are within the band limit and hence in |${\mathcal {B}}$|⁠, and we can solve the problem by dividing by |$\tilde{G}^\ast (u,v)$|⁠. A subtlety occurs at the band limit itself, |$\sqrt{u^2+v^2} =1/Q$|⁠, where both G(u, v) and its derivatives approach zero. Li & Mandelbaum (2023) considered the case where the Fourier transform of the PSF declines smoothly toward zero as (u, v) increases, in which case beyond some point both |$\tilde{G}$| and its derivatives can be treated as negligible to the accuracy required for a given project. With oversampled data, one may always accomplish this by a convolution in pre-processing even if the native PSF of the telescope has a sharper band limit.

For undersampled data, however, equations (C32) and (C33) are overdetermined: we have |${\mathcal {N}}(u,v)$| constraints (each u, v contributes a constraint) for the same 1 unknown |$\check{W}^{(a)}_{x,y}({\mathbb {F}}u,{\mathbb {F}}v)$|⁠. With the exception of some special unrealistic cases, there is no solution. If the data are weakly undersampled, then a filter in pre-processing can remove all of the aliased modes and we can work in the region where |${\mathcal {N}}(u,v)=1$|⁠. In the strongly undersampled case, however, there is no such simple fix for the problem.