Differencing and Coadding JWST Images with Matched Point-spread Function

We present an algorithm to derive difference images for data taken with JWST with matched point-spread functions (PSFs). It is based on the saccadic fast Fourier transform method but with revisions to accommodate the rotations and spatial variations of the PSFs. It allows for spatially varying kernels in B-spline form with separately controlled photometric scaling and Tikhonov kernel regularization for harnessing the ultimate fitting flexibility. We present this method using the JWST/NIRCam images of galaxy cluster Abell 2744 acquired in JWST Cycle 1 as the test data. The algorithm can be useful for time-domain source detection and differential photometry with JWST. It can also coadd images of multiple exposures taken at different field orientations. The coadded images preserve the sharpness of the central cores of the PSFs, and the positions and shapes of the objects are matched precisely with B-splines across the field.


INTRODUCTION
The James Webb Space Telescope (JWST) provides a unique opportunity for time-domain astronomy.The superb image quality enables the detection of faint transients out to the explosions of the first generation of massive stars and white dwarfs in the Universe (Riess & Livio 2006;Wang et al. 2017;Regős & Vinkó 2019;Lu et al. 2022).These diverse transients as different types of supernovae (SNe) can be direct probes to trace the cosmic star-formation history in the early Universe and expand our understanding substantially about the physics of the events at the epoch of the cosmic dawn.Recent observations have demonstrated that very faint transients (mag ∼ 29) are abundant in the near-infrared images taken by the JWST (Yan et al. 2023a,b;Chen et al. 2022a,b;Hu et al. 2022c;Yan et al. 2023c;De-Coursey et al. 2023;Chen et al. 2023).
Corresponding author: Lifan Wang lifan@tamu.edu,leihu@andrew.cmu.eduImage difference analysis is a major enabling technique in time-domain astronomy (Alard & Lupton 1998;Alard 2000;Bramich 2008;Miller et al. 2008;Becker 2015;Hu et al. 2022b).However, for the JWST, the highly structured point spread function (PSF) and the uncertainties in precisely matching the astrometries of images taken at different epochs can pose significant challenges in the identification of transients in the vicinity of bright sources like the central regions of galaxies.These include, for example, strongly gravitationally lensed transients in the vicinity of foreground lensing galaxies (Sheu et al. 2023) and the various nuclear transients including SNe as well as tidal disruption events (e.g., Zhu et al. 2021;Grishin et al. 2021;van Velzen et al. 2011;Regős et al. 2021).Calculating the difference images accurately is also important in highprecision differential photometry, which is employed in the detections of micro-lensing events (Mao & Paczynski 1991;Sumi et al. 2003Sumi et al. , 2006Sumi et al. , 2013) ) and exoplanet transits (e.g., Oelkers & Stassun 2018;Montalto et al. 2020).Recent works also reveal that the highly structured JWST/NIRCam PSF exhibits prominent spatial variation across the field of view up to 20% and shows
This study presents a method based on the saccadic fast Fourier transform (SFFT) algorithm (Hu et al. 2022b, hereafter H22) to accommodate such complicated PSFs as those of the JWST for image differences.Unlike previously published algorithms, the SFFT can model the spatial variations of the PSF using a B-Spline function and allows for more accurate image matching across the field in terms of PSF homogenization and compensation of astrometric misalignments.By design, our algorithm presents image subtraction in Fourier space, thereby achieving exceptional computational efficiency.Moreover, it can be used for both sparse and crowded stellar fields.The SFFT has been applied and extensively examined in several ongoing time-domain surveys and some transient analyses (e.g., Zhang et al. 2020;Palmese et al. 2022;Sun et al. 2022;Hu et al. 2017;Yang et al. 2022;Wang et al. 2022;Sheu et al. 2023).Recently, it also enabled new transient discoveries in JWST multi-epochs imaging observations (Hu et al. 2022c,a;Hu & Wang 2023a,b).The code of this study is built upon the SFFT algorithm proposed in H22 with significant improvements for the JWST.The improved version of SFFT is publicly available on Github1 , and a tutorial demonstrating how to perform and evaluate SFFT subtraction on JWST/NIRCam is also provided2 .The package is easily adaptable for data from other telescopes, such as with the Nancy Grace Roman Space Telescope (Roman; Spergel et al. 2015) the Legacy Survey of Space and Time (LSST; Ivezić et al. 2019).

TEST DATA
In this work, we demonstrate our method using the public JWST/NIRCam images of the Hubble Frontier Field (HFF) galaxy cluster Abell 2744 (Castellano et al. 2016;Merlin et al. 2016).In Cycle 1, JWST has carried out imaging observations of the well-studied lensing cluster at multiple epochs by several JWST programs, which has built a great data set for testing image difference methods that can be used for JWST time-domain analyses.
JWST revisited the cluster in November 2022 and conducted ultra-deep NIRCam observations with 4-6 hour exposures (∼ 29-30 AB magnitudes) in seven filters (F115W, F150W, F200W, F277W, F356W, F410M, and F444W) as a part of another early JWST program that targets Abell 2744: the Ultradeep NIRSpec and NIRCam ObserVations before the Epoch of Reionization (UNCOVER) Treasury survey (JWST- GO-2561, PI Labbé & Bezanson;Bezanson et al. 2022).We take the NIRCam observations collected by UNCOVER program on November 2, 2022, as the science images for subtractions.Note that the two JWST visits with different pointing and orientation overlap partially in a sky area covering ∼ 1.4 arcmin 2 centered at R.A. = 3  • 3673410 (see the JWST footprint of Abell 2744 in JWST Cycle 1 in Weaver et al. 2023).All the JWST data used in this paper can be found in MAST: 10.17909/7qx3-zt80.

IMAGE REDUCTION AND MOSAICS
We process the raw NIRCam data using the official STScI JWST Calibration Pipeline version 1.9.0 3 (Bushouse et al. 2022) in the context of jwst 1027.pmap 4 that includes in-flight reference calibration files released on 2023 February 20.To enhance the reduction quality, we incorporate a few augmentations into the official pipeline largely following the data processing prescriptions of the Cosmic Evolution Early Release Science Survey (CEERS; Finkelstein et al. 2022) program described in Bagley et al. (2023) We reduce the uncalibrated images through Stage 1 (Detector1Pipeline) of the JWST Calibration Pipeline that performs detector-level corrections and converts ramps to count-rate (slope) images.Given that the presence of stray light reflected off a secondary mirror support bar can introduce contamination to observations, giving rise to the characteristic "wisp" features on images, we address this problem by undertaking the subtraction of wisp patterns on count-rate images in F150W and F200W (most prominent filters) using the available wisp templates released on 2022 August 26 5 .Next, our processing turns to deal with the 1/f noise that is introduced during detectors read out (Schlawin et al. 2020) and manifest itself in random horizontal and vertical striping patterns.We identify and reduce the 1/f noise on count-rate images following the approach of amp-row and column subtraction proposed in Bagley et al. (2023).Note that the stripes are measured after the sources in the field have been well masked, unlike Bagley et al. (2023), here we identify the source mask using NoiseChisel (Akhlaghi & Ichikawa 2015), a noise-based method tailored for the detection of very extended and diffuse objects.However, we note that the correction of snowballs (Rigby et al. 2023), the circular defects on NIRCam images caused by cosmic-ray events, is yet to be included in our processing.We perform additional instrumental corrections and calibrations (e.g., flat-fielding) by Stage 2 (Image2Pipeline) of the JWST Calibration Pipeline.In this step, the count-rate images are converted to units of MJy Sr −1 .
We employ the JWST Stage 3 (Image3Pipeline) routine to create a mosaic image for each filter at the reference (science) epoch that combines all detectors and dithers with drizzling.This final stage consists of reduction steps including astrometric alignment (TweakReg), background matching (SkyMatch), outlier detection (OutlierDetection) and resampling (Resample).As the SkyMatch step in Stage 3 may have difficulties in background matching for the cases with small dither (Bagley et al. 2023), we skip this step in our processing but remove a single sky value for each individual image prior to Stage 3 using SkyMatch function 6 .In addition, we opt to deactivate the TweakReg step in Stage 3, which is used to calculate the coordinate transformations for aligning individual images to an absolute World Coordinate System (WCS) frame.Instead, we first combine all detector images of each exposure to a single exposure image using the Resample function.We then select the exposure that has the maximal overlapping area with other exposures as the agent of mosaic creation to provide a reference WCS frame.We use SourceExtractor (Bertin & Arnouts 1996, hereafter, SExtractor) to create a source catalog for each exposure and perform relative astrometry with respect to the sky coordinates measured on the agent exposure using SCAMP (Bertin 2006).This step harmonizes the WCS information across all exposures involved in the mosaic creation without invoking any absolute WCS reference.Upon running the OutlierDetection step to identify the outliers in the data, the exposures are drizzled to a 6 we adopt the mode with the parameter skymethod = "local" single mosaic with a drizzling parameter pixfrac = 1 using JWST Stage 3 routine7 .
Finally, we make a custom sky subtraction on the mosaic by NoiseChisel to eliminate the residual background.Recall that the relative astrometry above is separately performed for each mosaic, thereby not guaranteeing WCS consistency across the mosaics.We then undertake additional relative astrometric calibrations at the mosaic level with respect to an agent mosaic in F200W.The preprocessing ends with the final image resampling of all mosaics aligned to the agent mosaic using SWarp (Bertin 2010).Throughout the paper, unless explicitly stated otherwise, the term "mosaic" will refer to the astronomically aligned version after resampling.

IMAGE SUBTRACTION
Our image subtraction method developed for optimal difference imaging of JWST data is based on the SFFT algorithm proposed in H22.Here, we briefly recap the algorithm framework and introduce the improvements in Section 4.1.The specific subtraction scheme for JWST imaging observations is described in Section 4.2.

Improved SFFT Method
SFFT is an algorithm for astronomical image difference that presents the least-squares problem of image subtraction in Fourier space, bringing about a remarkable advancement in computational performance.It employs a δ-function basis to allow for ultimately flexible image matching with "shape-free" convolution kernels.Furthermore, SFFT can accommodate spatial variations of the matching kernel across the image field modeled by polynomials or B-splines.
For JWST observations, a number of modifications have been implemented within the SFFT framework.We summarize here: • The B-spline form spatial variation of convolution kernel, originally proposed in H22, is now well integrated into our software and extensively tested on JWST imaging data.
• The photometric scaling factor through convolution, that is, the sum of the matching kernel, can be separately controlled, though it was fully entangled with the kernel pixels in H22.For instance, the improved SFFT allows for a B-spline form matching kernel while imposing user-defined constraints on the kernel sum, such as being less flexible polynomials or constant across the field.
• The improved SFFT enables Tikhonov regularization (Press et al. 2007) to suppress the overfitting trend of matching kernel due to the ultimate flexibility of δ-function basis.In particular, we provide an option to adjust our regularization to create a trivial penalty when the optimal matching kernel has a profile close to a δ function, e.g., for subtractions between images that already have similar PSFs.

The SFFT: a recap
The problem of image subtraction can be written as the minimization of the difference image D (Alard & Lupton 1998;Alard 2000), in the form where R and S are input images with same dimension (N 0 , N 1 ), x and y are image coordinate indices in the ranges of [0, N 0 ) and [0, N 1 ), respectively.The spatially varying convolution of image R is denoted by R ⊛ K. K x,y is the matching kernel attached to image coordinate (x, y) with shape (2w 0 + 1, 2w 1 + 1), while u and v are kernel coordinate indices in the range of [−w 0 , w 0 +1) and [−w 1 , w 1 + 1), respectively.Note that the matching kernel K x,y , as a function of image coordinate, can vary across the image field to adapt to the ubiquitous spatial variations of PSF, photometric scaling and astrometry misalignment.Furthermore, the additional offset map B is used to account for the background difference between images R and S. Following Bramich (2008);Miller et al. (2008), SFFT decomposes the matching kernel K x,y into a complete δ-function basis K as follows: where (α, β) is the kernel coordinate of a non-center kernel pixel, i.e., (α, β) ̸ = ⃗ 08 .The basis K is defined as, and where ⃗ δ is a binary function such that ⃗ δ(ρ, ϵ) = 1 if ρ = ϵ = 0 and zero otherwise with ρ and ϵ being any integers.Note that all basis vectors other than K 00 have a zero-sum.As a result, the photometric scaling encapsulated in the convolution is uniquely determined by the coefficient Åxy00 .
The kernel spatial variation is completely encoded in the coefficients Åxy00 (particular) and Åxyαβ (general), which can be modeled by a two-dimensional smooth surface of either polynomials or B-splines across the image field: Åxy00 = rs års00 U rs (x, y) and The base functions U rs (V ij ) can have a k-order polynomial from: where ρ and ϵ are polynomial power indices in the range of [0, k] and [0, k − ρ], respectively; Alternatively, they may follow a more flexible B-spline form: where B ρ;k,tx (B ϵ;k,ty ) are the one-dimensional B-spline basis functions of given degree k and knots t x (t y ) along x (y) axis, the indices ρ and ϵ are in range of [0, k + N tx ] and [0, k+N ty ], respectively.We note that the improved SFFT, unlike in H22, formulates the coefficients Åxy00 and Åxyαβ independently.It signifies that the photometric scaling factor can be separately controlled, a feature that has been validated as beneficial in Bramich et al. (2013).
The differential background B(x, y) is also fitted by a polynomial/B-spline form function: where the base functions W pq are P pq or B pq .With an approximation (see Appendix A of H22) based on the fact that the scale of the spatial variations under consideration is significantly larger than that of convolution kernel, the Equation (1) can be rewritten as, where the notation • indicates circular convolution and U rs = U rs R and V ij = V ij R for abbreviations.The Equation ( 10) can be also derived using a different approach without the approximation (see the alternative perspective in Appendix A).
It is characteristic of the SFFT method to deliver the image subtraction problem forward into Fourier space.In the Fourier domain, we obtain where the symbols with a hat denote the Fourier transform of the images, a rs00 = N 0 N 1 års00 and a ijαβ = N 0 N 1 åijαβ .By Parseval's theorem, the least-squares optimization of difference image D is equivalent to minimizing its power spectrum in Fourier space.Let G = D D * ( * stands for complex conjugate) that represents the power spectrum of the difference image, and we can define the loss of the minimization in Fourier space as follows, where θ are the free parameters (a rs00 , a ijαβ and b pq ) of Equation 11that characterize the image subtraction.Thus, a least-square solution can be attained by optimizing its gradients to satisfy the condition ∇L 0 = 0.

Kernel regularization
The minimization of image subtraction with δfunction basis is prone to overfitting problems (Becker et al. 2012;Bramich et al. 2016;Masci et al. 2017).It is often characterized by irregularities (excessively noisy) resulting from matching kernels showing undesired adaptions to the noise of input data.The improved SFFT leverages Tikhonov regularization to address this overfitting issue.
Following the prescription outlined in Becker et al. (2012); Bramich et al. (2016), we regularize the shape of fitted matching kernels to have minimal local second derivatives by applying an additional penalty in the loss function.Since SFFT has invariant matching kernels across the field, we must implement the regularization for the convolutional kernels realized at different positions.With this consideration, we modify Equation 12as follows, where Kg is the flatten version of matching kernel K at image coordinate (x g , y g ).λ is an empirically tuned parameter to adjust the overall strength of regularization.w g is the specific weight of regularization at the coordinate (x g , y g ).The specific weights accommodate the situations when some subregions may be more susceptible to overfitting so that a higher local suppression is required accordingly.For simplicity, we only use a uniform weighting scheme, i.e., w g = 1, in this work.L is the symmetric (N K , N K ) Laplacian matrix (also see Bramich et al. (2016)) that represents the connectivity graph of kernel pixels, or equivalently, the standard Cartesian δ-function basis, with elements (14) Note that L K is an array of approximations to the local second derivative at each kernel pixel of K.It is locally calculated using at most five kernel pixels, generally fol- ) except for those kernel pixels adjacent to the boundary.
However, an important concern regarding the aforementioned definition lies in the fact that the regularization penalty remains non-trivial for a trivial matching kernel, i.e., a δ-function kernel with its only non-zero element at the kernel center.For example, the optimal matching kernels can be close to such a δ-function when the image subtractions are performed on images with PSFs already broadly aligned to each other.Under such circumstances, it is no longer reasonable to allow the regularization to hinder the image subtraction from finding a δ-function-like kernel as its optimal solution.Hence the improved SFFT offers users an option to remove the "barrier" by further tweaking the central rows of the Laplacian matrix as follows: where µ + represents the central five kernel pixels of shape "+" at (−1, 0), (1, 0), (0, −1), (0, 1) and (0, 0).The modification is equivalent to dropping the contributions of the local second derivatives at these central kernel pixels from the summation of Equation 13.As a result, the penalty is no longer a function of the pixel value of the kernel center at (0, 0).In this work, we zero out the central rows of L following Equation 15.

Hu & Wang
Invoking Equations 2-6 and using abbreviations V g ij = V ij (x g , y g ) and U g rs = U rs (x g , y g ), we flatten the matching kernel using the following equations and where η is the flatten index of non-center kernel pixel (α, β), i.e., η = (w 0 + α)(2w 1 + 1) + (w 1 + β); ϕ is the flatten index of the center kernel pixel (0, 0), i.e., ϕ = w 0 (2w 1 + 1) + w1;.With taking regularization into account and abbreviation of L = L T L, the loss function described in Equation 13 can be rewritten as It remains to minimize of the loss function L(θ) by optimizing its gradient with ∇L = 0, which still yields a linear system: One can solve the linear equations following Appendix B. Finally, taking the solution of Equation 19, we arrive at the difference image D by applying the solution to Equation 11.

Image Subtraction Scheme for JWST
This section presents the custom procedures of image subtractions we developed for JWST imaging data.To demonstrate our method, we perform the image subtraction between the NIRCam mosaic at the reference epoch (hereafter, reference mosaic) and the NIRCam mosaic at science epoch (hereafter, science mosaic) in each band (F115W, F150W, F200W, F277W, F356W, and F444W) created in Section 3.

Cross convolution
Firstly, we convolve the reference (science) mosaic with the PSF model of science (reference) mosaic, i.e., the so-called cross-convolution.More specifically, the cross-convolved mosaics R * and S * are generated following where R (S) is the reference (science) mosaic, and P R (P S ) is a corresponding PSF model at the image center retrieved from WebbPSF tool9 .The image subtraction approach with crossconvolution involved to homogenize PSFs was initially proposed by Gal-Yam et al. (2008) and Yuan & Akerlof (2008) and used in Zackay et al. (2016).Given that the JWST observations with a varying position angle lead to the rotation of the highly-structured PSF in the images taken at different epochs, a cross-convolution before a more sophisticated subtraction can broadly align the PSFs of input images to each other in a numerically stable way.As a result, it can effectively avoid deconvolution in the subsequent image subtraction that amplifies the noise and results in pathological correlation-induced patterns on the difference image (Zackay et al. 2016).
Note that we do not construct the PSFs from the observed frames; the WebbPSF models are used instead as approximations to the PSFs of the observations.We do this because it is not always possible to construct PSFs from the observed images.The cross-convolution brings the images to a level of approximately matching PSFs.Further processing by SFFT will match the PSFs with the matching kernel given in Equation 2.

SFFT preprocessing
As described in H22, it is essential for the SFFT method to solve the image subtraction on a masked version of input images rather than the original ones.A proper image mask serves to eliminate the affection of "bad" pixels in real observations, e.g., saturation, cosmic rays, and variable sources.It thus directs the SFFT minimization towards an unbiased construction of convolution kernels.
In this work, we create a binary mask shared for the given science and reference mosaics, largely following the preprocessing routine of sparse-flavor SFFT described in H22 with some modifications: 1.We identify the "bad" sources for SFFT subtraction using SExtractor.A field source is seen as "bad" if its SExtractor catalog value FLAG is zero and satisfies one of the following conditions: (1) it shows significant brightness change larger than one magnitude between reference and science mosaics measured by SExtractor catalog value MAG AUTO; (2) it is only detectable on one of the mosaic images.One may notice that the source exclusion here is a simplified version of the sparse-flavor SFFT in H22.We have skipped the identification of point-like sources as the deficiency of stars in JWST observations renders it difficult to detect a line feature regarding point sources using the Hough transformation.
2. We create a binary mask initialized by ones, and then zeros out all pixels associated with the identified "bad" objects or background regions, according to SExtractor segmentation maps of reference and science mosaics.
3. We tweak the binary mask by further clipping (converting mask value 1 to 0) all pixels that fall below three times of background standard deviation on any of the cross-convolved mosaics.The mask refinement is to make a better rejection of the noise-dominated pixels.Furthermore, our software allows for a final custom adjustment on the binary mask by clipping specific pixel regions defined in a SAOImage DS9 region file using polygon format.Although this function is not activated in our subtraction tests, it can be useful to exclude the saturated stars or other defects in the images when necessary.
Note that the masked regions (mask value 1) are considered reliable for solving the image subtraction.To create the masked version of cross-convolved mosaics, we replace the pixels of cross-convolved mosaics located in unmasked regions (mask value 0) with trivial zeros.

SFFT subtraction
We adopt the improved SFFT algorithm described in Section 4.1 to perform image subtraction between the cross-convolved mosaics produced in Section 4.2.1.This step aims to achieve better image matching to compensate the simple cross-convolution with limited matching accuracy.
We configure the tunable parameters of SFFT subtraction for JWST/NIRCam as follows: • The half-width of matching kernels is tuned to 11 and 5 pixels for short-wavelength and longwavelength channels, respectively.
• The general spatial variation of matching kernel (i.e., Equation 6) is modeled by a flexible B-spline function following Equation 8 with degree two and knots uniformly distributed at regular intervals of 300 (150) pixels along both x and y axes for shortwavelength (long-wavelength) channels; • The particular spatial variation of matching kernel sum (i.e., Equation 5) is fitted by a quadratic function following Equation 7. Note that a relatively low flexibility of kernel sum can preclude the photometric scaling from any undesired local structures.Otherwise, for instance, a variable object not successfully masked in Section 4.2.2 could incur a local underestimate or overestimate of photometric scaling at the position of the variable source; • Given that the background component of the mosaic images has been removed in Section 3, we assume a trivial flat differential background by a zero-order polynomial function (see Equation 9).
• The regularization strength parameter λ (see Equation 13) is finely tuned to be 3×10 −5 .We use 512 positions randomly sampled within the image frame to regular the solution.
We solve the image matching on the masked version of cross-convolved mosaics created in Section 4.2.2, and subsequently carry out the image subtraction between the (unmasked) cross-convolved mosaics by applying the matching solution10 .That is, where K is the spatially variant matching kernel solved from the masked mosaics.An initial difference image D * is thus obtained.We hereafter refer to it as undecorrelated difference.

Noise Decorrelation
Given that the undecorrelated difference must possess highly correlated background noise introduced in the cross-convolution and SFFT subtraction, we whiten the difference image through a noise decorrelation procedure following H22 (also see Zackay et al. 2016).
The formulation of noise decorrelation in H22 (see Appendix C of H22) does not take kernel spatial variations into account, while the SFFT matching kernels are spatially variant across the image.We adopt a straightforward strategy by dividing the image frame into a grid of small tiles so that we can perform the noise decorrelation separately for each tile, where the matching kernel is approximately constant.Here we use a grid with tile size as small as the matching kernel size.The noise decorrelation is performed by convolving undecorrelated difference with a spatially varying decorrelation kernel Q that is constant for each tile in the grid: and where σ R (σ S ) is the background standard deviation of reference mosaic R (science mosaic S), and K is the mathcing kernel realized at each tile center.Z is a factor that normalizes the decorrelation kernel Q to have a unit kernel sum for preserving the flux zero-point.As the noise decorrelation is derived in Fourier space, the decorrelation kernel Q is obtained after an inverse discrete Fourier transform (IDFT).To avoid confusion, we refer to the final difference image D as decorrelated difference.

Differential SNR
A key objective of optimal image subtraction is to yield a difference image where the flux residues of nonvariable sources and background are only dominated by their intrinsic statistical fluctuations, such as photon noise.The least-squares minimization of SFFT subtraction serves this purpose by generally suppressing the flux residues as far as possible.However, the effort may also be overkill in this direction: significant overfitting can adapt to the noise in data and artificially drive the residues toward zeros.The concern motivates us to modulate the loss function of SFFT subtraction with Tikhonov regularization.To evaluate whether the differential residues exactly stand at the optimal state, it is useful to derive the expected statistical noise of the difference image to provide a fiducial level of optimal subtraction.
Using a simple Monte Carlo sampling, we generate a propagated error map for decorrelated difference.JWST Stage 3 (Image3Pipeline) has produced an associated error map for each mosaic image.Firstly, we resample the JWST error map of the unaligned reference (science) mosaic to the target frame using SWarp.Subsequently, we calibrate the resampled error map using a constant scaling factor so that the background errors can coincide with the actual flux distribution measured on the reference (science) mosaic.Here, we have assumed that each pixel of reference (science) mosaic approximately follows an independent Gaussian distribution.
We randomly sample a zero-mean noise image 1024 times following the calibrated error map of reference (science) mosaic to trace the noise propagation through image subtraction.Next, we apply the same pixel operations involved in generating the decorrelated difference (cross convolution, SFFT image subtraction, and decorrelation) to each randomly sampled noise image pair.This step outputs 1024 propagated noise images at decorrelated difference stage, and we calculate their standard deviation at each pixel to construct the final propagated error map.
Finally, the signal-to-noise ratio (SNR) of decorrelated difference is calculated as the decorrelated difference D divided by the propagated error map N : hereafter referred to as differential SNR map.
The differential SNR map is a convenient check image to evaluate the quality of image subtraction.An optimal image subtraction should yield a differential SNR map that broadly adheres to an independent standard Gaussian distribution across the entire field.However, this criterion is overly idealistic as the assumption of independent noise distribution of input reference (science) mosaic is only an approximation.Instead, we opt for a more pragmatic standard by considering the measured distribution of differential SNRs on the background as the fiducial level.One can judge the quality of image subtraction at sources by comparing its differential SNRs to the benchmark background level.

PERFORMANCE AND COMPARISONS
In this section, we demonstrate our image differencing method described in Section 4.2 using NIRCam F200W mosaics of the Abell 2744 cluster created in Section 3. We present the subtraction performance of our method in Section 5.1.We compare our subtraction results with other methods in Section 5.2.

Subtraction Performance
Figure 1 shows the subtraction performance of the SFFT method over a section of Abell 2744 covering 25 ′′ × 25 ′′ centered at R.A. = 3 • • 5270483, Decl.= −30 • • 3645990.Our SFFT method can subtract the JWST data with few conspicuous subtraction-induced artifacts in the decorrelated difference, indicating an excellent image matching.The behavior of background noise correlation throughout the image subtraction process is depicted in Figure 1.Initially, the reference (science) mosaic exhibits only weak local correlation.After the subsequent cross-convolution, there is a notable surge in noise correlation, which remains pronounced in the undecorrelated difference.However, our noise decorrelation efficiently reduces the prominent correlation to a level comparable to the original reference (science) mosaic.
Figure 2 zooms in on a specific region of Abell 2744 centered at R.A. = 3 • • 5295458, Decl.= −30 • • 3648780, where our subtraction method unveiled a transient candidate, AT 2022acew (Hu et al. 2022c).AT 2022acew appears at a sky position between two galaxies.One of the galaxies is a typical elliptical galaxy (SExtractor detection G1 in Figure 2), while the other might be a (lensed) galaxy that exhibits more intricate structures (SExtractor detections G2, G3, and G4 in Figure 2).Notably, a bright foreground star is near the galaxies (SExtractor detection P1 in Figure 2).This selected region showcases how our subtraction method performs on different sources with diverse morphology.The differential SNR map in Figure 2 reveals a distinct and prominent detection of AT 2022acew, despite its faint nature.Meanwhile, most pixels associated with the galaxies exhibit desirable low differential SNRs comparable to the background fiducial level.As shown in Figure 2, the subtraction quality on the central regions of bright objects G1, G3, and G4 have reached a (nearly) optimal level.For the object G2, our subtraction reveals structured residues at its core.Aperture photometry of G2 on the SFFT difference shows a significant positive net flux, which suggests that the residues are not likely a stand-alone subtraction-induced artifact, i.e., they may originate from the true flux variability of an AGN or a nuclear transient candidate at G2.At the bright star P1, a circular subtraction artifact appears, but its contamination is well confined to a relatively small area with a radius of 0 ′′ • 22 (equivalently, ∼3.5 times of FWHM) to the star's centroid.We conduct a statistical analysis on the differential SNRs across the field presented in Figure 2. It also confirms that the differential SNRs in signal-dominated regions are broadly consistent with the fiducial background level.Only a moderate performance deterioration appears at the brightest pixels (SNR > 100).We emphasize that the subtraction quality on these regions is critical for the search and accu- rate photometry of the transients close to galaxy nuclei.Figure 2 also traces the PSF changes through the processes of image subtraction.One may notice that the noise decorrelation not only whitens the noise but also narrows down the PSF size of the difference image, rendering it comparable to original un-convolved mosaics.

Comparisons with Other Subtraction Methods
To compare the performance of different subtraction methods, we conduct the following tests on NIRCam F200W mosaics of the Abell 2744 cluster using various approaches: (a) Direct HOTPANTS: The subtraction is directly performed between science mosaic and reference mosaic using software HOTPANTS 11 (Becker 2015), 11 HOTPANTS is configured with the following parameters: -c=t, -r=11, -ko=2, -bgo=0, -rss=30, -nsx=NX/100 and -nsy=NY/100.ferential SNR map is produced based on the simplified noise propagation.
(c) Direct B-spline SFFT: The subtraction is directly performed between science mosaic and reference mosaic using B-spline form SFFT. The same binary mask is used as in the direct polynomial SFFT test, and SFFT configurations completely follow Section 4.2.3.Like the direct HOTPANTS test, a differential SNR map is generated using the simplified noise propagation.
(d) Cross-convolved arithmetic subtraction: The subtraction is performed between cross-convolved science mosaic and cross-convolved reference mosaic using a straightforward pixel-wise arithmetic operation of subtraction.Subsequently, we whiten the correlated background noise of the resulting difference image introduced in cross-convolution fol-lowing the noise decorrelation method outlined in Section 4.2.4.Finally, a differential SNR map is created using the Monte Carlo sampling described in Section 4.2.5.
(e) Cross-convolved polynomial SFFT: The subtraction is conducted between cross-convolved science mosaic and cross-convolved reference mosaic using polynomial form SFFT, with the same configurations as the direct polynomial SFFT test.We proceed to decorrelate the background noise of the resulting difference image introduced in crossconvolution and SFFT subtraction following the strategy in Section 4.2.4.Again, we create a differential SNR map with the Monte Carlo sampling.
(f) Cross-convolved B-spline SFFT: The subtraction test introduced in Section 4.2. Figure 3 shows the subtraction performances on NIR-Cam mosaics of the Abell 2744 cluster using the abovementioned methods.The statistics of differential SNR indicate that the two B-spline SFFT tests, with or without cross-convolution, have the best subtraction quality among the methods.Their differential SNRs in signal-dominated regions can broadly reach the fiducial background levels.It suggests that the introduction of flexible B-spline form spatial variations for convolutional kernels, coupled with the δ-function kernel basis of SFFT, contributes to achieving a nearly optimal level of image subtraction for JWST/NIRCam data.We note that the less optimal direct polynomial SFFT subtraction already demonstrates a considerable improvement over the direct HOTPANTS (see H a and H b in Figure 3).Despite both approaches adopting polynomial modeling for the PSF variation, SFFT's utilization of a δ-function basis proves more effective for handling sophisticated PSF homogenization and compensation of astrometric misalignments for JWST, compared to HOTPANTS with Gaussian basis functions.A detailed evaluation of polynomial form SFFT against other existing implementations (including HOTPANTS) is presented in H22 with an emphasis on ground-based observation tests.
Considering the distributions of differential SNR alone, one may conclude that the direct B-spline SFFT appears to be as good as the cross-convolved B-spline SFFT.However, it is important to emphasize that the inclusion of cross-convolution, in addition to its numerical stability, improves the subtraction quality for point sources.The advantage is evident for the bright star P1: the direct B-spline SFFT gives rise to a more extended subtraction artifact, as shown by the residues surrounding the square-like pattern in Figure 3. Due to the limited number of stars in our test data, this advantage is not adequately reflected in the statistics.
To visualize the flexibility of the B-spline kernel, we have extracted the matching kernels at different image positions for the cross-convolved polynomial and Bspline SFFTs. Figure 4 shows the spatial variations of five central kernel pixels and the kernel sum for these two subtraction tests.Note that the photometric scaling for cross-convolved B-spline SFFT has been constrained to a low-degree quadratic function in the same form as cross-convolved polynomial SFFT.The two kernel sum surfaces depicted in Figure 4 are roughly consistent with each other and close to a flat unit plane.As anticipated, central kernel pixels for the B-spline form exhibit more structures across the field that can better adapt to the PSF spatial variation and compensate for local astrometrical misalignments.

THE EFFECT OF KERNEL REGULARIZATION
Overfitting is an important concern regarding any image subtraction method, especially with δ-function basis.In this section, we demonstrate how the regularization technique assists the SFFT subtraction in properly harnessing the least-squares minimization process to mitigate overfitting.
We reconduct our subtraction test on the NIRCam F200W mosaics six times, each with a different regularization strength λ ranging from 0 to 10 −3 .To evaluate the quality of undecorrelated differences of these subtraction tests, we calculate a propagated error map for each undecorrelated difference using the Monte Carlo ) of the undecorrelated differential SNR maps for SFFT subtraction tests with increasing regularization strength λ being 0, 10 −6 , 10 −5 , 3 10 −5 , 10 −4 and 10 −3 , respectively.In each postage stamp, the cyan circle highlights the central region of the examined galaxy with a radius of 0 ′′ • 1763, twice the FWHM of the cross-convolved PSF model.In contrast, the yellow circle marks a region selected from the ambient background with a radius of 15 pixels (0 ′′ • 4650).The standard deviation σ of the pixel values enclosed by the cyan (yellow) circle is accordingly labeled at the left corner.(K1) -(K6): From left to right, the panels show the matching kernels realized at the position of the examined galaxy for SFFT subtraction tests with the increasing λ as labeled.(H1) -(H6): From left to right, the histograms show the probability distributions of undecorrelated differential SNRs, enclosed by different pixel regions shown as different histogram colors, for SFFT subtraction tests with the increasing λ.The pixel regions are identified following Figure 2 but defined on the cross-convolved image pair.Each distribution's standard deviation σ is labeled at the right border in the corresponding color.
sampling described in Section 4.2.5 and generate a corresponding undecorrelated differential SNR map.
Figure 5 shows the impact of Tikhonov regularization on SFFT subtraction tests at an examined galaxy at R.A. = 3 • • 5216322, Decl.= −30 • • 3627870 in Abell 2744.When no regularization is applied, one can notice that the overfitting manifests itself in the highly flattened undecorrelated differential SNRs at the examined galaxy compared with the neighboring background.Accordingly, the matching kernel at the examined galaxy is excessively noisy, signifying an undesirable adaptation to the galaxy's Poisson noise.Statistically, the overall distribution of undecorrelated differential SNRs at signal-dominated regions is also pathological: it is even better than the fiducial background level.As shown in Figure 5, adjusting the regularization strength can adequately address the overfitting problem.Increasing the regularization parameter λ results in a smoother matching kernel and broadens the distribution of undecorrelated differential SNRs in signal-dominated regions to be more reasonable.It is worth mentioning that the regularization strength we used in Section 4.2, λ = 3 × 10 −5 , is a moderate value that effectively curbs overfitting without significantly compromising the quality of subtraction.
Despite its effectiveness in regularizing the noise levels of difference images, the Tikhonov regularization may not be the ultimate solution to address the overfitting problem in image subtraction.The penalty term of Tikhonov regularization only modulates the shape of the matching kernels but does not guarantee or properly quantify the optimal solution.It is not mathematically well-defined, and users must fine-tune the regularization strength as a hyper-parameter.For optimal subtraction, the ultimate goal is to eliminate any structured residues in the difference image for all the objects without true variability.It is equivalent to minimizing the information content or maximizing the entropy of the residual images.In a future study, we will aim to incorporate a term to the loss function that can steer the fitting towards maximum entropy.

IMAGE COADD
Naturally, the precise image matching accomplished by the B-spline form of SFFT can also be employed in the image co-addition to construct deeper mosaics.We illustrate this capability by performing a co-addition of the two NIRCam F200W mosaics of Abell 2744 that were used in Section 5.
The co-addition scheme mostly inherits the procedures described in Section 4.2.We use the same cross-convolution and B-spline form SFFT to align the input mosaics.The matched mosaics are subjected to weighted co-addition (instead of subtraction) followed by a noise decorrelation.More specifically, we co-add the reference mosaic R and science mosaic S following and where simple inverse variance weights w S = 1/σ 2 S and w R = 1/σ 2 R are adopted.Again, Z c is a normalization factor, and the noise decorrelation is performed over the same grid of tiles.The example image co-add is shown in Figure 6.We note the sharp NIRCam PSF is preserved through the image co-addition.As expected, the co-add image bears reduced background noise and increased signal levels.With the detection parameters set to be identical 12 , the image co-add lead to detection limits that is ∼ 0.25 mag fainter, consistent with the expectation from Poisson statistics of the noise.The advantage of using SFFT for image co-adding is that it automatically aligns and matches the positions and profiles of the objects in the image field and creates co-add images that preserve the sharpness of the central cores of the JWST PSFs.

SUMMARY AND CONCLUSIONS
We have introduced an image differencing pipeline adapted to improve difference image analyses of JWST/NIRCam observations.We briefly summarize here the major steps of the pipeline: 1. Starting from uncalibrated NIRCam observations, we utilize an augmented version of the official STScI JWST Calibration Pipeline to create the reference and science mosaics taken at different epochs (see Section 3).
2. Next, we perform cross-convolution by convolving the reference mosaic with the PSF model of the science mosaic, and vice versa, using the PSF models provided by WebbPSF (see Section 4.2.1).This process can broadly align the PSFs of the two mosaics and minimize the numerical instability during the subsequent image subtraction.
3. We introduce a B-spline form of kernel variations in the SFFT method and modulate it by the Tikhonov regularization to perform the image subtraction between the cross-convolved mosaics (see Section 4.2.2 and Section 4.2.3).This step aims to achieve an accurate image matching with PSF homogenization and corrections of astrometrical misalignments.This version of the SFFT method is characterized by several new features, outlined as follows: (1) It allows for flexible B-spline functions to depict the spatial variation of matching kernels, and in parallel, it can independently control the photometric scaling (kernel sum), typically by modeling with polynomial functions with lower degrees of freedom.
(2) Tikhonov regularization has been incorporated to suppress the undesired noise adaptions (i.e., overfitting problem) caused by the high level of flexibility in the fitting process.The new method is detailed in Section 4.
4. We apply the prescription of noise decorrelation outlined in H22 to the difference image obtained from the SFFT subtraction (see Section 4.2.4).This step effectively removes the convolutioninduced correlations of background noise, and the FWHM of the central cores of the PSF becomes comparable to that in the original mosaics.
5. Finally, the pipeline provides a differential SNR map as a check image to evaluate the quality of subtraction and diagnose possible overfitting (see Section 4.2.5).
This paper demonstrates the performance of the pipeline using JWST/NIRCam imaging data of the Abell 2744 cluster acquired in JWST Cycle 1 by the GLASS and UNCOVER programs.We exemplify that our method can achieve high subtraction quality, for which the residues on signal-dominated regions statistically harmonize with those at the fiducial background level.Moreover, we show the regularization technique can properly suppress the overfitting trend stemming from the high degree of freedom in SFFT subtraction.We also make a comparison of subtraction performance using different techniques.Among them, a regularized B-spline form SFFT coupled with cross-convolution can achieve the best quality of image subtraction for the JWST/NIRCam data.The method can also be used for accurate co-adding of JWST images.The algorithm is potentially useful in studying variable stars/transients in nearby galaxies, in searching for exoplanets through microlensing, and in finding SNe, especially those that are gravitationally lensed by nearby, relatively bright galaxies (e.g., Yuan et al. 2022;Riess et al. 2023;Mayker Chen et al. 2023;Penny et al. 2019;Chen et al. 2022c).
Figure 7 presents more examples of the subtraction performance on Abell 2744 cluster in different filter bands, zoomed in on several galaxies.The data reduction and image subtraction follow the steps described in Section 3 and Section 4.2, respectively.The PSF mismatch can lead to spurious features extending to ≳ 0 ′′ .5 from the center of a bright point source if the difference images are calculated using the original mosaics (See Figures 3 and 7 for examples).While the total area affected by the PSF mismatch is usually insignificant compared to the entire image field, the central regions of galaxies require more careful treatment and can not be ignored for many important studies.For transient searches around diffuse galaxies, the effect of the PSF structure may not matter much.Our method is most important in searching for transients around galaxies with bright central cores.Some examples are the SNe and tidal disruption events close to the central regions of galaxies, AGN variabilities, and SNe lensed by foreground galaxies (whose Einstein ring is of the size ∼ 1 ′′ ).
As shown in Figure 7, the difference images of these galaxies derived from our algorithm are clean and reveal no variabilities, whereas residual patterns are conspicuous when using simple arithmetic subtraction.Our algorithm enables robust discoveries of variabilities and transients down to the nuclei of galaxies.The 8th row of Figure 7 shows such an example (see the white cross).It shows a transient phenomenon that is significantly de-

Figure 1 .
Figure 1.Image subtraction performance of the SFFT method on NIRCam (F200W) mosaics of the Abell 2744 cluster.Top panels: from left to right, reference mosaic (R), science mosaic (S), and images decorrelated difference (D).Bottom panels: from left to right, the cross-convolved reference mosaic (R * ), cross-convolved science mosaic (S * ), and undecorrelated difference (D * ).The inset panel of each subplot shows the local noise correlation measured on the samples of the background pixels enclosed in the blue square.More specifically, it represents the covariance matrix of a multivariate random vector, consisting of 25 adjacent background pixels (see definition in Appendix C of H22).The black dashed square indicates the position of transient candidate AT 2022acew.The arrows between different panels represent the corresponding operations as labeled on the right side.

Figure 2 .
Figure 2. Image subtraction performance of the SFFT method on NIRCam (F200W) mosaics in a close-up view around the transient candidate AT 2022acew (white cross) in the Abell 2744 cluster.The grayscale images presented in the first three panel columns are the same as Figure 1 but zoomed in a 3 ′′ × 3 ′′ section around AT 2022acew.However, the inset panels show the corresponding PSF models with measured FWHM values (in units of arcsec) labeled at the top.Note that the PSF models of R and S are retrieved from WebbPSF, while those of R * and D are generated by applying the corresponding convolutions to original WebbPSF models.In the reference mosaic R, the overlaid contour in red (purple) marks the detection profiles for pixels with flux over 10 (100) times the sky background standard deviation σs on both reference mosaic R and science mosaic S. The top rightmost panel shows the differential SNR map D, while the bottom rightmost panel presents the probability distributions (histograms in light colors) of differential SNRs within different regions across the image as labeled.The cyan crosses mark the centroid positions of the five most prominent detections in this field measured by SExtractor on S. The arrows between different image panels represent the corresponding operations as labeled on the right side.A Gaussian profile is fitted for each distribution and illustrated by a dark-colored curve.The skewness of each distribution is marked in the corresponding color.The vertical black short lines at the right side of the x-axis indicate the differential SNR values of AT 2022acew higher than five.

Figure 3 .
Figure 3. Image subtraction performance on NIRCam (F200W) mosaics of the Abell 2744 cluster using different methods.(S): A 25 ′′ × 25 ′′ image section of science mosaic as shown in Figure 1, and its inset panel shows a 3 ′′ × 3 ′′ zoomed-in postage stamp around AT 2022acew as shown Figure 2. ( Da) -( Df ): Postage stamps (with the same view as the inset panel of S) of differential SNR maps from the six subtraction tests in Section 5.2: direct HOTPANTS, direct SFFT, direct B-spline SFFT, cross-convolved arithmetic subtraction, cross-convolved polynomial SFFT, and cross-convolved B-spline SFFT, respectively.In each postage stamp, the cyan circle highlights the central region of the elliptical galaxy at the right corner, centered at R.A. = 3 • • 5293333, Decl.= −30 • • 3648317, with a radius of 0 ′′ • 3098, five times of FWHM size of the PSF model of science mosaic; The yellow circle marks a region selected from the ambient background with the same radius.The standard deviation σ of the pixel values enclosed by the cyan (yellow) circle is accordingly labeled at the left corner.The five prominent detections shown in Figure 2 are also marked with cyan crosses in the inset panel of S and the panel Df .(Ha) -(H f ): The histograms, from left to right, show the probability distributions of differential SNRs for different pixel regions for the six subtraction tests, respectively.The pixel regions are identified following Figure 2.Each distribution's standard deviation σ is labeled at the right border using the corresponding color.
Figure 4. Comparisons of spatial variations of matching kernels between the polynomial form SFFT and B-spline form SFFT. The leftmost panel shows the 25 ′′ × 25 ′′ image section of science mosaic same as Figure 1 and Figure 2. The top-right panels that are labeled by K P (α, β) show the spatial variations across the image field of the matching kernel pixel at (α, β) for cross-convolved polynomial SFFT test.In contrast, the top rightmost panel shows the spatial variations of the kernel sum.The bottom-right panels have the same style but for the cross-convolved B-spline SFFT test.In these panels, the dashed white lines indicate the positions of the inner knots of B-splines.The black dashed square indicates the position of AT 2022acew, same as Figure 3.

Figure 5 .
Figure 5.The effect of Tikhonov regularizations in SFFT subtraction.(S): A ′′ × 25 ′′ pixel section of science mosaic, centered at R.A. = 3 • • 5234850, Decl.= −30 • • 3615950.(S * ): A 4 ′′ × 4 ′′ postage stamp of the cross-convolved science mosaic centered at the examined galaxy.( D * 1 ) -( D * 6 ): Postage stamps (with a view same as S *) of the undecorrelated differential SNR maps for SFFT subtraction tests with increasing regularization strength λ being 0, 10 −6 , 10 −5 , 3 10 −5 , 10 −4 and 10 −3 , respectively.In each postage stamp, the cyan circle highlights the central region of the examined galaxy with a radius of 0 ′′• 1763, twice the FWHM of the cross-convolved PSF model.In contrast, the yellow circle marks a region selected from the ambient background with a radius of 15 pixels (0 ′′• 4650).The standard deviation σ of the pixel values enclosed by the cyan (yellow) circle is accordingly labeled at the left corner.(K1) -(K6): From left to right, the panels show the matching kernels realized at the position of the examined galaxy for SFFT subtraction tests with the increasing λ as labeled.(H1) -(H6): From left to right, the histograms show the probability distributions of undecorrelated differential SNRs, enclosed by different pixel regions shown as different histogram colors, for SFFT subtraction tests with the increasing λ.The pixel regions are identified following Figure2but defined on the cross-convolved image pair.Each distribution's standard deviation σ is labeled at the right border in the corresponding color.

Figure 6 .
Figure 6.Image co-add of NICam (F200W) mosaics of Abell 2744 cluster using the SFFT method.From left to right, the mosaic image was taken on 2022 Jun.28-29, 2022 Nov. 2, and the SFFT co-add mosaic.The inset panels show the corresponding PSF models with measured FWHM values (in units of arcsec) labeled at the top.Note the PSF models of R and S are retrieved from WebbPSF, while the PSF models of co-added image C are generated by applying the corresponding convolutions to the WebbPSF models.In each panel, the blue circles are the sources detected by SExtractor on each panel image.The red squares represent the sources detected only in the co-added image.The SExtractor parameters were kept identical in detecting these sources.The black dashed circles indicate the pixel regions contaminated by the uncorrected snowballs and sources inside them are not shown.