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

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\'e 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.


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 et al. 1996) and of the twopoint correlations of shear (Van Waerbeke et al. 2000;Bacon et al. 2000; Wittman et al. 2000;Kaiser 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  8 to precision better than a few percent (Hikage et al. 2019;Hamana et al. 2020;Asgari et al. 2021;Amon et al. 2022;Secco 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 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  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 (Laureĳs et al. 2011); and the Nancy Grace Roman Space Telescope, whose Wide Field Instrument will survey the sky at high angular resolution in 4 NIR bands spanning 0.9-2.0m (in the Reference Survey design 1 ; 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  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 /, where  is the wavelength of observation and  is the diameter of the entrance pupil.If the pixel scale  is larger than /(2), then the image is undersampled in the Nyquist sense: there are Fourier modes present in the image with more than 1 2 cycle per sample, with the consequence that the images cannot be unambiguously interpolated, and thus the image intensity  () 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. 2 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 et al. 2021;Finner et al. 2023).One particular step is calibration of the shear estimator -that is, determining how the measured ellipticity of a galaxy   in a catalog responds to an applied shear   . 3While one might hope for an ellipticity measurement algo- 1 The Reference Survey is an example survey used during the design and construction phases to show that the Roman mission meets its requirements.The actual survey conducted will be designed through a community process and could be different. 2If significant, read noise can impose an additional penalty for a wellsampled survey due to the the low sky flux on each pixel. 3We include an index since  and  are 2-component quantities.
rithm that has unit shear response ⟨  ⟩/  =    , 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  max , and the shear operation moves modes across the  =  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  cut <  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 et al. 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). 4Further 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 percent-level biases (Kannawadi et al. 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. 2020Sheldon et al. , 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--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.
resolution images of a small sample of the galaxies in the same wavelength range to determine the statistical properties of the galaxy images at  >  max , with the application of using Hubble Space Telescope images to calibrate ground-based shapes.
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 et al. 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 et al. 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.

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  in and in general have their own rotations, distortions, PSFs, and masks; and produce a well-sampled output image at an output pixel scale  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 point spread function (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  ×  matrix T, where  is the number of output pixels and  is the number of input pixels.These statements apply to many of the common algorithms that have been applied to groundbased weak lensing data sets (either for shapes, source selection, or photo-s): for example, linear predictive codes for interpolating bad pixels (Bosch et al. 2018, §4.5); Lanczos-3 (Bertin et al. 2002;Bosch et al. 2018, §3.3) or polynomial (Huff et al. 2014, §4.4) interpolation for resampling; and pixel-domain (Bernstein & Jarvis 2002, §7) or Fourier-domain (Huff et al. 2014, §4.1) rounding kernels, and PSF Gaussianization (Hildebrandt et al. 2012, §3;Kuĳken et al. 2015, §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 wellsampled 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 et al. (2011) searches for a matrix T that accomplishes all 4 steps with minimum noise and target PSF error (as measured with quadratic metrics). 5It 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   is determined by the overlap of output pixel  with a "shrunken" version of input pixel .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 et al. (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.

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/  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.

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 -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 et al. 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 et al. 2011;Spergel et al. 2015): 4 in J129, and 3 in each of H158 and F184.The Y106 band strategy was not required to achieve full sampling and so used 3 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 Sec.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 Sec.5.4.
The "simple" input model also comes with a PSF model (corresponding to a flat spectral energy distribution or SED in   ) 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).

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 (Troxel et al. 2023, Appendix B).Out of the 18 chips initially selected for flight, an average of 3.01% 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 et al. 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 −4 per pixel, which is the product of the 10 −6 cm 2 pixel area, the 140 s exposure time, and the expected rate of 5.5 events cm −2 s −1 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 (Rigby et al. 2023, §6.6), Roman has smaller pixels (so potentially a higher leakage into neighbors 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.

Injected sources
We have also created two layers that are grids of injected stars.The stars have unit flux and are injected on a HEALPix6 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, 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.

Simulated noise fields
We construct two types of simulated noise fields: white noise and 1/  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/  noise input is constructed as follows.For each of the 32 readout channels, we construct a 1-dimensional Fourier domain signal of length  = 2 20 : , where the real and imaginary components of S are Gaussians with mean 0 and variance 1/(2|   |) (except for the constant mode S0 = 0), where   = / is the normalized frequency.This is then discrete Fourier transformed,   =  /2−1 =−  /2 S e 2 i  / , 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  min and  max is ≈ ln(  max /  min ).We take a length 2 19 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/  noise that we can use to assess the impact on weak lensing with coadded images.An example of a 1/  noise field generated in this way is shown in Fig. 3.
Both noise inputs use the numpy.random.default_rngrandom 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.

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  t × 4096 × 4224 cubes (where  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 (Troxel et al. 2023, Appendix B) 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% of the pixels.Lab noise field analysis and results will be presented in detail in the forthcoming paper (Laliotis et al., in prep).

IMAGE COADDITION
The image coaddition process is an updated version of the Imcom framework for coaddition of postage stamps (Rowe et al. 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.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 −4 ).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 Sec.4.1, and indicate where there have been algorithm updates.We describe the block and mosaic construction in Sec.4.2.The choice of output PSFs -essentially determining the resolution of the coadded images -is discussed in Sec.4.3.
A diagram of the block coaddition script is shown in Fig. 5.

Postage stamp coaddition
The Imcom formalism is developed in Rowe et al. (2011) 7 , 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  vector {  } −1 =0 . 8Masked pixels can simply be removed from this vector with a corresponding decrease in .These are to be combined into an output image with  pixels:   = −1 =0     , where T is an  ×  coaddition matrix.We denote the positions of each pixel center by {  } −1 =0 (input images) or {  } −1 =0 (output image).If the input images have effective PSF   , then the output pixel  has an effective PSF at offset  given by 9 PSF out,  () (1) 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   that make PSF out,  uniform and close to Γ. Quantitatively, we wish to minimize the leakage function 10 , where the square norm is written in Fourier space with weighting Υ().The numerator   is the square norm of the difference between 7 See also Bolton & Schlegel (2010), who develop a related formalism for combining the pixels in 2D fiber spectra into a 1D spectrum. 8Since the new implementation is in Python+C, in this paper we follow the Python+C convention of indices starting at 0. 9 There is a subtlety in defining the "PSF" of a pixel when the PSF does not act exactly as a convolution, i.e., Eq. ( 5) of Mandelbaum et al. (2023) is not satisfied.In general a pixel at position   responds to a point source at some position  according to a PSF  (  −  ).The term "point spread function" literally refers to this function at fixed  and varying   .Equation (1) fixes the output pixel   and varies , which is a more natural choice when pixels are discrete and positions on the sky are continuous.One might think of this object as a "reverse" of the usual PSF. 10 We intended to take Υ() = 1.As noted in Section 4.4, the production runs were run with an alternative weighting in the objective function,  ) is defined by a center, a map projection, and a number of blocks (BLOCK×BLOCK).Each block (panel [b]) is itself composed of postage stamps; we make an  1 ×  1 array, with padding of PAD postage stampps around the rim so that the blocks overlap.The postage stamps (panel [c]) are composed of an  2 ×  2 grid of output pixels, with a transition region around the edge that is merged at the block processing level before writing to a FITS file.The postage stamp is built from all un-masked input pixels in all input images within a given acceptance radius of the stamp.the as-realized output PSF and the target, while the denominator  is the square norm of the target (thus forming a dimensionless ratio).
We will refer to   / itself as the "leakage" (smaller is better), and to the quantity −10 log 10 (  /) as the "fidelity" (larger is better: one can think of this as a suppression of PSF error in decibels).If the  ×  noise covariance matrix of the input is N, then the  ×  output covariance is  = TNT T .Imcom also tries to minimize the output covariance Σ  for input white noise (N = I × ). 11aving two objective functions, leakage   / 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   +   Σ  , 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   given   is described in Rowe et al. (2011); the algorithm performs a bisection search in log 10   , 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 et al. (2011) requires the construction of the correlations [  •  ] () 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 [  •  ] () 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, 1 12 " kernel here.
• Rather than simply specifying a target leakage   / or a target noise Σ  , we specify a maximum noise (Σ  ) max as our first priority and a maximum leakage (  /) 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.

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[, ,,] contains the pixel in layer , observation ID/SCA , and pixel (, ).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 centers 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 Tmatrix ( §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 2 pixels; we define a sequence {  } 2 =1 .If we have a "left" postage stamp  left (, ) (whose output region, excluding the transition pixels, is  <  cut ) and a "right" postage stamp  right (, ) (whose output region, again excluding the transition pixels, is  ≥  cut ), then in the transition region  cut −  ≤  <  cut + , we may make a merged image, (3) This approach does not affect the output PSF as long as the transition sequence satisfies   +  2+1− = 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 Eq. ( 3) is trivially extendable to the  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 center, and the Roman footprint can be efficiently tiled with square regions 12 , although we may revisit this choice in the future.At an angle  from the projection center, the stereographic projection introduces a lowest-order plate scale error of 1 4  2 , or 2.4 × 10 −5 at the corner of a 0.8 degree square ( = 0.8/ √ 2 × /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 (BLOCK − 1, 0) is in the lower-right (southwest corner).

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 (T-matrix), rather 12 The Mercator projection satisfies the first two criteria, and would be appropriate to a great circle strip, e.g., Huff et al. (2014); but since strips converge on a sphere, the unique area must be tapered.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  = / is the diffraction scale at the central wavelength of that filter,  1 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/ 3 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   / 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.)An example of the resolution versus output fidelity trade-off in the Roman Y106 band is shown in Fig. 7.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.
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 et al. (2021) to use a larger weighting scale for galaxy ellipticities 13 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.

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 th input image overlapping that postage stamp, the cosmic ray mask was read from the th 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 center of every 2 × 2 group of postage stamps) went off 13 As an example of how these ideas are related, let us consider the quadrupole moment of an image  at radius , Let us think about what happens once full sampling is recovered so that we can really think of this operation as an integral.If a new image  sm is constructed as the convolution of  with a Gaussian of scale  0 , then   [  sm ] = ( / tot ) 6  tot [𝐼 ] where  2 tot =  2 +  2 0 .That is, the second moment of the smeared image is (aside from a scaling factor) the second moment of the original image at a larger scale.
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 2 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, Υ() = 1 in the language of Rowe et al. (2011).An experimental feature to place more weight on the longerwavelength modes, setting Υ() 2 , was accidentally left on for the production runs with  Υ = 1 and  Υ = 3 2 × 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.

RESULTS
We present some general results, examples, and statistics of the output PSF leakage   / 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.

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 PSFthis 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 −4 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.

Output PSF fidelity
We describe the output PSF fidelity in terms of the quantity −10 log 10 (  /).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% 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 5 (Y106), 4 (J129), 3 (H158), and 2 (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 Y10614 and so we accepted fewer dither positions.The trouble comes

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  and (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 ( 1 2 pixel offset), with a wavelength of /.This behavior is shown in Fig. 12 for the simple case of 2 input images in 1 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 ∼ 1 4 SCA dither in the direction is ∼ 1.2%, 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 patternthey 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  > 2 dithered images since then  ( − 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, 1 4 pixels of each other can be expressed analytically15 , and is  12.An illustration of a Moiré pattern, with 2 input images in 1 dimension.The top row (blue + signs) represents the pixel positions in image #1, and the middle row (red × signs) represents the pixel positions in image #2.The interlaced pixel positions (from both #1 and #2) are shown in the bottom row, and one sees alternating regions of "degenerate" combination (in the linear algebra sense: the two input images give information at the same locations) and "ideal" combination (half pixel offsets, so effectively achieving twice the sampling rate).The wavelength of the pattern is the harmonic difference of the two input wavelengths, 1/[1/ − 1/ (1 + ) ] ≈ /, where  ≪ 1 is the fractional difference in plate scale.3 16 = 18.75% for  = 3; 1 16 = 6.25% for  = 4; and 5 256 ≈ 1.95% for  = 5.If one goes up to 2 dimensions but with no rolls, then this argument applies in both  and  directions: for example, for  = 3, one would expect 1 − (1 − 3 16 ) 2 ≈ 34% of the area to have samples clustered within 1 4  in , , or both (this drops to 12% for  = 4 and 4% for  = 5).For the case with arbitrary rolls, we are unaware of any practical analytical results on these the "left" side of the cluster, and then a probability 1  4 for each later dither to be between 0 and 1 4 pixel to its right.
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; Refregier et al. 2012;Melchior & Viola 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.

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/, 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 2 ) with a coverage of only 3 exposures.This generally consists of 2 exposures at one roll angle and 1 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.35pixels (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.

CONCLUSION
We have presented an implementation of the Imcom image combination algorithm (Rowe et al. 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 "onepoint" 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. ) Once again, we can see the features in the coverage map (Fig. 1).The regions of large cosmetic defects in the SCAs can be seen in the output noise maps as strings of higher-noise splotches corresponding to the dithering sequence.

Noise standard deviation (as fraction of input images)
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 et al. 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.and Since  is even, S is symmetric.
The most obvious choice is to weight all Fourier modes in the band limit equally: () = 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  is constant for || <  and 0 otherwise).
Since the matrix S −1 does not depend on , but only on the parameters  and , 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 , but we want the tightest control for small .This is the case, for example, of a diffraction limited image whose power spectrum only truly drops to zero at  =  but is small at  ≳ /2.We might then modify the objective function to be () = 1 − ||/.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.,  =1−   () 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 (Ω) 2-point interpolation scheme that has zero error for the constant mode, i.e., the  = 0 mode.This problem can be solved by adding a formally infinite -function to (), i.e., we construct a new objective function: where e is a length 2 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 −1 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  depends only on the choice of interpolation scheme and can be pre-tabulated.Note that by construction, the entries in  sum to unity, guaranteeing that  =1−  ′  = 1.The use of any vector  that sums to unity would fix the background conservation issue (for example, the rescaling in Bernstein & Gruen 2014 corresponds to   /    ), but Eq. (A18) does so at the least cost to the weighted mean square error over the specified range in .
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  and small  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 −8 ), 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 Eq. (A8), we may define an error metric that is not an integral over Fourier modes, but rather a discrete sum: where {  }  =1 are the abscissae (rescaled so that 0 <   < 1); {  }  =1 are the weights, assumed to be positive; and we have enforced symmetry of positive and negative frequencies.We require  ≥  so that the matrix S has full rank.Then where we defined   = 2  .The advantage of this is that the matrix S can be factorized as S = M T M, where M is the 2 × 2 matrix (These results can be verified by direct multiplication, as well as the identity cos( − ) = cos  cos  + sin  sin .) Next we perform the singular value decomposition of M, so that M = UDV, where U is a 2 × 2 orthogonal matrix; D is a 2 × 2 matrix with all zeros except for the diagonal {   } 2 =1 ; and V is a 2 × 2 orthogonal matrix.In this case, Eq. (A10) simplifies to This can be expressed as where the 2 × 2 matrix H can be written as where in the prefactor we define  + =   .Note that H depends on , , and the weights   and   , but not on ; thus once a scheme is chosen, it can be pre-computed.The expression contains a division by a singular value   , but since S = M T M, the condition ratio of M can be better than that of S (if  = , 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   and   corresponding to 2-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  > , we have not found significant gains that offset the increased computational cost (we must compute 2 trigonometric functions for each 1D interpolation).Thus we have retained  = .
The discrete window interpolators with  =  are in principle exact at the spatial frequencies   , since then M is invertible so S −1 b = M −1 c and one can show that an input function  () that is a cosine or sine wave leads to a vector f T = (  (1 − ),  (2 − ) ...  ()) that is a constant times a row of M. Then the interpolated value, f () = f T M −1 c = (M T−1 f )c, collapses.However, in finiteprecision 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, sinc  = sin()/().Also the oversampling relative to Nyquist is 1/(2): thus  = 1 4 is 2× Nyquist,  = 1 8 is 4× Nyquist, etc.

A4.1 6-point interpolation
A well-studied choice in weak lensing image processing is 4× Nyquist sampling (e.g., by zero-padding an FFT) and 6-point ( = 3) interpolation.Bernstein & Gruen (2014) found that multiplicative errors and ghosting could be reduced to the < 0.1% level with quintic polynomial interpolation (P3).The square window LSE method gives the interpolation weights for the "S3, 1 8 " scheme: where the  (S) coefficients are given in Table A2.Similarly, for the triangle window LSE method or "T3, 1 8 " scheme, we get: The  vectors are also shown in the table, so one may use Eq.(A17) to construct the background-conserving kernels   ( |S ′ 3, 1 8 ) and   ( |T ′ 3, 1 8 ).A comparison of the errors for the various 6-point interpolation schemes is shown in the upper panel of Fig. A1.The quintic interpolation scheme performs best at very low spatial frequencies ( < 1 16 ).The Lanczos-3 kernel has large errors, several tenths of a percent, across all of the low spatial frequencies, but with this sacrifice it gains improved behavior at intermediate spatial frequencies ( ∼ 0.3).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.memmapfor the input data typically reduced the memory usage by ∼ 30%.
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 this appendix, 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 , which is the number of pixels corresponding to one cycle at the maximum spatial frequency;  < 2 corresponds to undersampling by the usual Nyquist crtierion, and smaller  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 <  < 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 near infrared bands on Euclid have "strong undersampling" ( < 1), so the mathematical framework of standard Metacalibration does not apply to these cases.The alternative Deep Metacalibration algorithm (Zhang et al. 2023) in principle could address these problems when  < 1, but many Fourier modes crisscrossing the (, )-plane must be excluded from the re-convolution PSF.The analytic shear response method of Zhang et al. (2023) faces similar issues: in the weakly undersampled case, convolution in preprocessing 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 ( in = 1).

C1 Fourier transforms and operators
We perform continuous Fourier transforms in the "waves per pixel" convention: where (C6) When we discuss the magnitudes of the wave vectors, it is helpful to keep in mind that in the weak lensing regime || + || < 1, the maximum singular value of M −1 T ,,  is that is, for (, ) in a circle of radius , the farthest that M −1 T  , (, ) can be from the origin is Λ.
A discretely sampled function  , , measured only at integer pixel positions (, ) ∈ Z 2 , instead has a Fourier transform: where is the first Brillouin zone or the unit square centered on the origin (shaded square in Fig. C1).The two types of Fourier transforms are related to each other via the X ("Shah" or Dirac comb) function: which is its own Fourier transform ( X = X).Ifwediscretelysample a function  at the pixel positions,  , =  (, ) for integer , , then One sees that the Fourier transform of the discretely sampled data contains a term associated with the original field (Δ = Δ = 0), as well as other terms offset by an integer numbers of cycles per pixel.We define the fractional part operator that maps  to the range − 1 2 to + 1 2 .Then (F, F) is the single mode in the first Brillouin zone that aliases to (, ).

C2 Astronomical scenes and available information
Let us define an image which is a convolution of intrinsic projected scene  and the effective point spread function  (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 Eq. ( C13) is infinite, in practice  has a band limit or maximum spatial frequency present: G (, ) = 0 if √  2 +  2 ≥ 1/, where  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  = /( in ), where  is the wavelength of light,  is the diameter of the entrance pupil, and  in is the pixel scale.For a ground-based observatory, atmospheric smearing usually suppresses the highest spatial frequencies.While turbulent contributions to G usually decline smoothly rather than going to zero at a "hard" cutoff in spatial frequency, there is still some radius in the (, )-plane beyond which G (, ) is negligible.
The sampling parameter  restricts the number N (, ) of nonzero terms in Eq. (C13), with the number of terms growing as  decreases (more undersampled).To determine N (, ), we draw a square lattice with unit spacing centered at (, ), and count the number of lattice points in a circle of radius 1/ centered at the origin (Fig. C1).The results are visualized in Figure C2.Each panel shows the region B in (, )-space covered by a discrete Fourier transform centered on (0,0) and this goes from − 1 2 to 1 2 cycles per pixel on both axes.The color scale shows the number of terms that contribute in Eq. (C13), i.e., the number of Fourier modes in the original image that alias to (, ) 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 √  2 +  2 ).It can be seen that more Fourier modes are aliased as the sampling factor  decreases (i.e., increasing the radius 1/ of the circle).Intuitively, as the circle begins to overflow outside the region, neighboring circles centered at (, ) ∈ Z 2 begin to overlap and the frequency in the overlapped region is aliased.

√
2 +  2 < 1/} of radius 1/.If 1 <  < 2 (cases b and c), then the image is weakly undersampled: while some modes are aliased, others are not and N −1 (1) is non-empty; the largest disc contained in it is D 1−1/ .If  < 1 (cases d and beyond), then the image is strongly undersampled and N −1 (1) is the null set: there are no Fourier modes S(, ) that can be recovered.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  exp images with dither offsets  1 ...  exp , then the th image can be thought of as measuring the displaced scene T   .In this case, the discrete Fourier transform of the th image is f  (, ) = e 2 i(   +   ) ∑︁ (Δ,Δ) ∈Z 2 e 2 i(   Δ+   Δ)

C3 Recovering information with translational dithers
× G ( + Δ,  + Δ) S( + Δ,  + Δ) (C14) (compare to Eq. C13).This is a system of  exp equations with N (, ) unknowns, and it can be solved if the  exp × N (, ) matrix L with coefficients where  represents one of the  exp images and  represents one of the (Δ, Δ) terms in the sum, has rank N (, ).Note that it is necessary but not sufficient that  exp ≥ N (, ).
Cases with more small dither steps for other sampling cases could be discussed: for example, the  exp = 5 "knight's move" pattern   = ( 2 5 , 1 5 ) works for case (d), and Gonzaga et al. (2012) discuss the  exp = 8 point dither (for cases e-g) and  exp = 9 point dither (for case h).[The  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 (Δ, Δ) = (cos , sin ), and the "network" of coupled modes (Eq.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  exp A where A is the survey area (again in units where  in = 1), whereas the number of "unknowns" (output Fourier modes) is A/ 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  exp ≥ / 2 . (C18) The Imcom formalism (Rowe et al. 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 E the region where S(, ) can be reconstructed, with the understanding that for  exp = 1 dither, E = 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 E, and actually contains a "buffer" region before one reaches the edge of E; see, e.g., Bernstein (2010, §4).We will see that this behavior is generic.
We start with the "standard" Metacalibration (Huff & Mandelbaum 2017;Sheldon & Huff 2017), and then investigate "deep" Metacalibration (Zhang et al. 2023), with a particular emphasis on the operations that are potentially problematic on undersampled data (e.g.Hoekstra et al. 2021).Finally, we consider the use of pixel responses as in Li & Mandelbaum (2023).

(C20)
We suppose that the re-convolution PSF has a Fourier transform Gr ( + Δ,  + Δ) with support in some circular region out to radius  r .
For non-zero shear, the Fourier modes M −1 T 0,,0 ( + Δ,  + Δ) do not alias to each other.In this case, we requires knowledge of S at all the points that contribute to the sum.That means that for every point ( ′ ,  ′ ) in the support of Gr , M −1 T 0,,0 ( ′ ,  ′ ) must be in the reconstructible region E.This leads to the following cases: • For oversampled data ( > 2), E = D 1/ .This means that the radius of support of Gr must have  r ≤ (1 − ||)/ (see Eq. C7).Since the original PSF G has support out to a radius 1/, 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 .
• For weakly undersampled data (1 <  < 2) without dithering, the largest disc contained within E has radius 1− −1 .Thus we can carry out the standard metacalibration operation if the re-convolution PSF satisfies  r ≤ (1 − ||) (1 −  −1 ).Note in this case that as the native PSF gets smaller at fixed pixel scale ( gets smaller), the Fourierspace cutoff  r must get smaller and hence the re-convolution PSF has to get bigger in real space.
• For strongly undersampled data ( < 1) without dithering, 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 is again limited by (1 − ||)/.
For single-epoch Roman images, "standard" Metacalibration does not apply in Y106.It theoretically can be applied in J129 ( = 1.021), however 1 −  −1 is so small that it would not be practical (the reconvolution 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 et al. 2023).Here there is both a deep scene  deep and a wide scene  wide that have to be rendered with a reconvolution PSF  r .The shear operation, Eq. (C20), is applied only to the deep data, whereas no shear is applied to the wide data.In the context of both ground-and spacebased 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  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 (, )-plane E deep .Because of the large number of dithers, this region may be larger than E wide , but it is still bounded by the radius 1/.

(C22)
This means that for a given (, ), if there is a (Δ, Δ) ∈ Z 2 satisfying Gwide (+Δ, +Δ) ≠ 0 and M −1 T 0,,0 (+Δ, +Δ) ∉ E deep , (C23) then the linear combination appearing in Eq. (C22) must simply be zero, and Gr must be zero for all modes that alias to (, ).And in particular, if there is a lattice point ( +Δ,  +Δ) in the troublesome annulus: then there is an orientation for the shear such that Eq. ( 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 Gr (, ) = 0 in these regions.In practice, since one uses Deep Metacalibration with small values of ||, the troublesome regions of (, )-space are those near the circle √︁ ( + Δ) 2 + ( + Δ) 2 = 1/, 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  = 0.834 (Y106 band) and || = 0.01, then Gr must go to zero at a radius of 0.187 cycles pixel −1 , 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 groundbased 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 (, ) ∈ Z 2 can be written in terms of the sky scene  via is the pixel basis function.[Li & Mandelbaum (2023) define  , (, ), which is our Φ , (, )/(2) 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

Figure 1 .
Figure 1.The coverage (number of exposures) in each of the 4 bands in the 48 × 48 arcmin region considered in this paper.Each sub-panel shows one of the filters.The 18-chip "pawprint" feature of the Roman focal plane is easily visible, as is the presence of two roll angles from the two passes in each filter.

Figure 2 .
Figure 2.An example of a simulated mask.This is the 512 × 512 lower-right corner of SCA 11 in a H158-band observation (ID 8836).Reference pixels are shown in dark yellow; permanently masked pixels are shown in black; and pixels rejected in this observation only due to cosmic rays are shown in red-orange.

Figure 4 .
Figure 4.The hierarchical structure of mosaic coadds in this paper.The mosaic (panel [a]) is defined by a center, a map projection, and a number of blocks (BLOCK×BLOCK).Each block(panel [b]) is itself composed of postage stamps; we make an  1 ×  1 array, with padding of PAD postage stampps around the rim so that the blocks overlap.The postage stamps(panel [c]) are composed of an  2 ×  2 grid of output pixels, with a transition region around the edge that is merged at the block processing level before writing to a FITS file.The postage stamp is built from all un-masked input pixels in all input images within a given acceptance radius of the stamp.

Figure 5 .
Figure 5.The workflow for coaddition of a block in this paper.There are two repositories: the postage stamp coaddition (furry-parakeet) and the mosaic driver with interfaces to the simulations (fluffy-garbanzo).

Figure 6 .
Figure6.The input and output modulation transfer functions (MTFs: absolute value of the Fourier transform of the PSF).The upper panel shows the MTFs as computed in the Exposure Time Calculator(Hirata et al. 2013) based on the requirement aberrations.Two curves are shown for each filter, either for wave vector aligned with the pixel grid or at a 45 • angle; for F184, the difference is not visible on the plot.The Nyquist frequency at the native (input) pixel scale is marked.The lower panel shows the output PSFs, with fractions of the Nyquist frequency at the output pixel scale  out marked.Note that all of the input images are undersampled (the MTF is nonzero for spatial frequencies extending past  Ny,in ) but the output images are well sampled, with the MTF dropping to negligible levels (< 10 −4 ) by ∼ 0.36  Ny,out in all filters.

Figure 7 .
Figure 7. Examples of how output fidelity depends on the output PSF and on the number of exposures covering a given region.The first three rows show the fidelity maps, −10 log 10 (  / ), for a block of the Roman Y106 band centered at right ascension 53.5142 • and declination −40.3898 • .On this scale, 60 corresponds to the specified leakage   / = 10 −6; lower values correspond to worse matching of the output PSF to the target.We show three resolutions of the target PSF ( FWHM indicated), with sharpest resolution at the top.The bottom row shows the coverage map (number of exposures).The first column shows the full 1.08 × 1.08 arcmin block; the other three columns show zoom-ins of particular regions (the tics are in units of output pixels).Note the particularly poor reconstruction when  exp = 3 (second column), corresponding to an intersection of chip gaps.Even when  exp = 5 (third column), if the target PSF is too small (top), one sees imprinting of the 0.11 arcsec (4.4 output pixel) input pixel grid where there is difficulty interpolating from the small number of samples to certain sub-pixel positions.We also see losses where a cosmetic defect reduces the effective number of exposures to 4. Similar behavior can be seen in the right column, which overlaps with a 4-exposure region.

Figure 8 .
Figure 8.An assortment of objects in the simulated coadded images.Each panel is a Y106 (blue) + J129 (green) + F184 (red) composite, with a field of 18 arcsec on a side, with a scale stretched from −8 to +1200 / 2 in /exposure.(a) An elliptical galaxy at  = 0.18.(b) A deep field; the objects labeled are galaxies at  = 0.81 (object 1),  = 0.79 (2), and  = 0.93 (3).(c) A bright star (visible magnitude  550 nm = 16.4);note the square pattern imprint of the postage stamp boundaries and the diffraction spikes.(d) A disc galaxy at  = 0.22.(e) An assortment of galaxies; objects 1 and 2 are at  = 0.34, and object 3 at  = 2.81 appears red since the Balmer + 4000Å feature is redshifted into the H158 band.(f) An assortment of stars (1 and 2) and galaxies (3-5 are in a group at  = 0.28).

Figure 9 .
Figure 9.An example of a simulated coadded star in the Y106 band.The left part of the figure shows the 6 images of the star, at the as-observed roll angles in the two passes, in the "simple" DC2 simulation.The right panel shows the coadded map.The Imcom coadd (this project) is shown at top; the drizzled coadd is shown below.A logarithmic color stretch is used.The sky level (61 e/p) is present in the inputs but has been removed in the outputs.The black squares in the inputs are masked pixels (for the Imcom run).

FidelityFigure 10 .
Figure10.The fidelity map over the 48 × 48 arcmin region simulated, in each of the 4 bands.The lowest values are seen in Y106, since it has the most undersampled PSF and the number of dither positions is lower than in J129.Note that the same features in the coverage map (Fig.1) also appear here: the output PSF is not as well matched to the target in regions of intersecting chip gaps.

Figure 13 .Figure 14 .
Figure 13.An example of Moiré patterns in the simulated coadd in the Y106 filter.The left panel is a coverage map in a region with intersecting chip gaps in the two passes.The middle panel shows the RMS of the "whitenoise1" noise field, as measured in 1.25 × 1.25 arcsec postage stamps.The wave-like features are Moiré patterns; note that these patterns usually start and end at a chip gap, and that a range of wave vectors are present.The right panel shows a zoom-in of a 16.5 × 16.5 arcsec (660 × 660 output pixel) region of the coadded "whitenoise1" noise field.Increased variance in the form of ripples with specific wave vectors can be seen in the degeneracy bands.

Figure 15 .
Figure 15.The noise map over the 48 × 48 arcmin region simulated, in each of the 4 bands, measured as the RMS of the "whitenoise1" noise field, as measured in 1.25 × 1.25 arcsec postage stamps.(The Y106 panel is a zoom-out of the center panel in Fig.13.)Once again, we can see the features in the coverage map (Fig.1).The regions of large cosmetic defects in the SCAs can be seen in the output noise maps as strings of higher-noise splotches corresponding to the dithering sequence.

Figure C2 .
Figure C2.The number N (, ) of Fourier modes contributing to each part of the (, )-plane for discretely sampled data.A sequence of cases is shown depending on the sampling parameter , ranging from oversampled (a) to increasingly undersampled cases (b,c,d,...).The numbers in square brackets in the caption are the minimum and maximum number of Fourier modes that contribute.

Table 1 .
The layers of input data used in this paper.Data may be in "all" of the bands (Y106, J129, H158, and F184), or in the indicated bands only (the decision to add the last two layers was made after processing of Y106 and F184 had started).All layers are generated as 4088 × 4088 SCA images prior to being fed into the coaddition.

Table 2 .
The permanent mask flags based on the BADPIX HDU of the SCA files.The last line shows the union of these flags; since some pixels have multiple flags, the sum of the fractions is greater than ALL.

Table 3 .
Sizes and dimensions used in this work.

Table 4 .
Kannawadi et al. (2021) and output PSFs.The PSF effective area is defined in the square-norm sense Ω psf = 1/ ∫ [ (s) ] 2  2 s and expressed in input pixels.The sampling factors are in the convention of / in ; this is the same convention used inRowe et al. (2011), but is a factor of 2 larger than the convention ofKannawadi et al. (2021).The output smearing is listed as the full width at half maximum of the Gaussian in units of native pixels.
The cumulative distribution of the output PSF fidelity − log 10 (  / ) for each of the Roman filters.The vertical axis shows the fraction of the output 0.8 × 0.8 degree mosaic with PSF fidelity worse than the indicated value.

Table A1 .
The histogram of the output image fidelity in the Y106 filter in the 3-exposure region of block(44, 34).Results are shown for 5 values of the charge diffusion length: 0 (no charge diffusion), 0.1, 0.2, 0.3, and 0.4 pixels rms per axis.Note the improvement in PSF fidelity as the charge diffusion is increased since the input images are better sampled.Interpolation methods considered in this appendix.The allowed parameters are  (kernel size: based on 2 nearest pixels),  (number of abscissae, for the discretely optimized kernels), and  (band limit for which the kernel is optimized).

Table A2 .
Interpolation weight coefficients  Eqs.A26 and A27).The  coefficients are also given so that one may build the S ′ 3, 1 8 and T ′ 3, 1 8 schemes.The data in these tables are included as Tables_A2A3.txt in the furry-parakeet GitHub repository.

Table B1 .
Computing resources for the image coaddition runs in this paper.
)(This is the description of shear in terms of the polar decomposition of the 2 × 2 matrix M  ,,  .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 T   (, ) = e 2 i(  +  ) f (, )