All frames were debiased and flat-fielded in the standard manner. Sky-flats were used to divide out the overall illumination of the CCD in each filter and dome flats used to derive the pixel-to-pixel variations.

Each frame was cleaned individually of bad pixels including cosmic rays. Such pixels were identified by fitting a 5x5 pixel 2D Gaussian of width less than the seeing on each frame to each pixel in the image. The residuals of the difference between the data and the model effectively discriminate between neighboring pixels that are consistent with the seeing and those that are associated with chip defects and cosmic rays. Once such pixels were identified simple linear interpolation using good neighboring pixels was used to replace a bad pixel. The advantage of cleaning individual frames in this way is to allow a weighted combination of frames taken under variable seeing to be made which is clean of defects, as compared with a more robust combination such as a median which can degrade the final PSF of a combination of frames with very different seeing.

Registration of all frames to common coordinates was made with $\sim 50$ bright, unsaturated stars on each frame whose centers were determined using the DAOPHOT PHOT routine and the IRAF routines GEOMAP and GEOTRAN. A 5th order polynomial interpolant was used in the geometric transformation and flux was conserved. Final registration of frames was accurate to within 0.1 pixels. The PSFs of stars on all frames were well sampled with minimum FWHM's of $2.1$ pixels, thus resampling errors were small and PSFs were not significantly degraded in this registration process. At this point frames were combined for different timescales: for the analysis of sub-night timescales sequences typically five images were combined giving images with exposure times of 12.5$^m$ in each filter, and on nightly timescales all frames taken in each filter were combined.

\subsection{PSF Matching with a Fourier Algorithm}

In some test frames taken in similar seeing C92 was able to show that the residuals in the difference image were comparable to the photon noise. As stressed in {\S 1}, however, the residuals in the difference image are expected to be primarily influenced by differences in seeing between a pair of frames to be differenced. Since a typical PSF approximates a Gaussian we have found that simply convolving one frame with a Gaussian kernel to match the FWHM of the PSF of the better seeing frame is very effective. In general, systematic residuals in the difference frame can be minimized by matching the PSFs between frames as closely as possible.

In order to monitor partially resolved globular clusters in M31 for nova eruptions over a number of years through differential photometry Ciardullo, Tamblyn \& Phillips (1990) developed a Fourier algorithm to match seeing variations of different epochs frames. In essence each good seeing frame was degraded to a common seeing value. With frames of identical seeing the relative flux of all non-varying point-like sources in a frame is the same for any given finite photometric aperture and thus meaningful differential photometry can be performed between frames. Phillips \& Davis (1995) have developed this algorithm further and we are very grateful to Andrew Phillips for providing us with his software, which we have employed as the basis of our PSF matching technique. To apply it to its full extent to a typical wide field imaging telescope is not straightforward. In {\S 3.3} we discuss the details and the algorithm we have developed to solve these problems and apply this to the KPNO data. In {\S 3.6} we show how the PSF matching can be simplified with an optical corrector for a wide field imager and demonstrate its application to the VATT data.

First, ignoring noise, consider a frame $r$ with a narrower PSF than a frame $i$,

$$ i = r * k, \eqno (1) $$ \noindent

where $k$ is a convolution kernel that describes the difference in the seeing and guiding between the two frames. The Convolution Theorem states that the Fourier transform (FT) of these three variables has the form,

$$ FT (i) = FT (r) \times FT (k), \eqno (2) $$ \noindent

then,

$$ k = FT [ FT (i) / FT (r) ], \eqno (3) $$ \noindent

Thus $k$ can be determined empirically with a high S/N isolated star on a frame pair. Convolving the good seeing frame with this kernel will in principle provide a match to the PSF of the poorer seeing frame. In practice the determination of $k$ is not straightforward since the high frequency components of the FT become dominated by the noise in the wings of the PSF where the signal is weakest. An effective method of dealing with this problem was determined by Ciardullo, Tamblyn \& Phillips (1990). Since the FT of a typical PSF is roughly Gaussian the convolution kernel will also be approximately Gaussian. By modeling the high S/N low-frequency components of the PSF FT with an elliptical Gaussian these noise-contaminated components can be replaced with the model fit yielding a convolution kernel close to the ideal.

Figure 1 shows an application of this method to some KPNO data. We show a 128x128 pixel subimage near the center of an image, which clearly shows a background of unresolved stars. A suitable bright star near this region has been used to determine $k$. Typically the kernel is determined over a region extending up to five times the FWHM of the PSF. After photometric scaling of the frame pair ({\S 3.5}) and the degradation of the better seeing frame with its convolution with $k$, the difference image shows the clear removal of structure in the background, including a good subtraction of all but the brightest star on the frame which is saturated. However, also shown in Figure 1 is a subimage of the difference frame located about 500 pixels away from the location where $k$ was determined. It is clear that this region is plagued by large systematic residuals associated with a poorly matched PSF in this region. Critically, this is not simply a problem affecting the bright stars on the frame, but also the background of unresolved stars where we are most concerned with achieving the best subtraction. This demonstrates that there is no unique solution of $k$ applicable to the entire frame.

Application of the Phillips \& Davis (1995) algorithm to the KPNO data. The left side shows a $128\times 128$ pixel subimage close to the center of the original frame (upper left panel) with the mean galaxy background subtracted from the image, and its difference image (lower left panel). A suitable star close to this region has been used to empirically determine the PSF matching convolution kernel to match the image pair being differenced ({\S 3.2}). All structure in the unsubtracted data has been effectively removed; the residuals around the brightest star are due to it being saturated on the CCD. However, applying the same convolution kernel to a region located 500 pixels away (upper right panel) shows large systematic residuals in its difference frame on the scale of the PSF (lower right panel). This indicates a poor match of the PSF in this region and shows that there is no unique solution to the matching convolution kernel applicable to the entire frame. Effective subtraction of the full image can only be done by modeling the spatially varying PSF kernel using the limited number of appropriate PSF matching stars on the frame {\S 3.4}. Once applied in this case the quality of the subtraction for the entire frame becomes comparable to the lower left panel. The inset image in the lower right corner is the new difference image located in the region of the box in the upper right panel. A clear detection of a point source is now made, which was almost completely hidden in the systematic residuals in the first attempt at matching the PSF. All differences are displayed in the same intensity range; the intensity range of the upper panels is five times larger.

\subsection{Understanding the Spatial Dependence of the Matching PSF Function}

Figure 2 illustrates why $k$ becomes spatially dependent in a frame pair to be differenced. We plot the FWHM of bright, unsaturated stars (we ignore partially resolved M31 globular clusters) as a function of distance from some radial symmetry point of the PSF variation close to the center of the CCD ({\S 3.5} shows how this was located) for two frames taken close together in time and in the same filter. The radial variation of the PSF is clear and is quite dramatic at furthest radial distance which is located near one corner of the CCD. However, also apparent is the different form this radial variability takes in the two separate frames. This form is affected by small differences in the focus of the telescope between the two exposures. It is this time variability of the spatial functional form of the PSF that, independent of seeing and guiding variations which are themselves not spatially dependent across the camera, causes the PSF matching convolution kernel $k$ to become a function of position on the frame.

The FWHM of point sources on two separate R band KPNO frames taken two hours apart (shown as open and filled circle points) as a function of distance from the PSF radial-variation symmetry point on each frame. The functional form of the radial variability of the FWHM is different in the two images due to a small change in the telescope focus between the two exposures.

In principle the solution to this problem is simply to match frames in a piecemeal fashion by PSF-matching subregions around bright, isolated and unsaturated stars which can be used to determine the correct convolution kernel for the local region. In practice, however, the extent over which good subtraction can be made local to a PSF matching star is frequently only as large as $50\times 50$ pixels. Such a method therefore requires a star suitable for determining a convolution kernel in over 1500 subframes to adequately match the entire $2048^2$ pixel frame. The requirements for the choice of PSF matching stars are rather stringent: the star must (i) have high S/N, (ii) must not be saturated in any part of the PSF or corrupted by any defects such as cosmic rays or bad pixels, (iii) must have an amplitude that significantly exceeds the surface brightness fluctuations associated with the unresolved stars in the galaxy, and (iv) must be isolated from any bright neighbors over the extent to which $k$ is being determined. In these data the number of point-sources on each frame that satisfy these criteria is typically $\sim 200$. Furthermore, few of these sources reside in the bulge region where the underlying bulge light requires correspondingly brighter stars for good PSF matching, leaving large and important areas of the frame without good kernel determinations. In the next section we describe how this problem was addressed by deriving a model for the spatial variability of the matching convolution kernel over the entire frame from every available location where local convolution kernels could be measured.

\subsection{Derivation of Spatially Dependent PSF-matching Convolution Kernel}

The search for variables in various frames proceeded by differencing each frame with a high S/N template reference image with good seeing. For this purpose the KPNO reference images for both filters were generated from the combined image stacks on the final night since this was the best quality data. The analysis of each frame began with measures of the raw convolution kernels for the individual frame/reference frame pair, using the method just described, at the location of all isolated (no companions within a radius of 15 pixels), high S/N and unsaturated stars - a list comprising 220 stars. Each kernel is 15x15 pixels in size, which is significantly larger than the typical PSF with FWHM in the range 2-4 pixels. Since the PSF's to be matched have essentially the same functional form, the largest value in this array is almost exclusively located in the central pixel.

To determine the behavior of the overall spatial variability of the kernels a polynomial surface fit was made to the value of this pixel at the $x,y$ positions of each kernel determination on the frame. Figure 3 shows a contour plot of the result of a typical high order polynomial fit for a frame pair. In this example the fit excludes the smaller region marked by the inner contour with a value of unity where the relative size of the PSF on each frame reverses. The unit contour represents the location where the PSFs are close to identical. In Figure 2 this would be located where the FWHMs radial forms cross at a radial distance of $\sim 450$ pixels.

The spatial variability of the PSF matching convolution kernel for a pair of images in the KPNO 4-m data. The empirical PSF matching convolution kernel has been measured at the location of points marked on the plot. The $\times$ marks are distinguished from $+$ marks as locations where the ratio of the size of the PSF on the two frames either exceeds or is less than unity. A 6th order polynomial surface fit (as a function of CCD coordinates) has been made to the central pixel intensity of the matching convolution kernels at the $+$ mark locations and is shown as a contour plot. The inner contour marks the region where the PSFs are close to the same size on the original frames and within this region a similar determination of the PSF matching convolution kernel is made but with the images reversed.

Since the normal variation of the PSF within these frames generally exceeds the variation in seeing between frames the situation in shown in Figure 3, with an inner elliptical or near circular region within which the relative size of the PSF's of the frame pair reverses, is quite typical. The overall radial variability of the convolution kernel in both regions is always concentric about the same point, however the position of this point was found to vary up to 100 pixels between frames from the symmetry point apparent in Figure 3. The location of this point most likely represents the principal axis of the telescope which is slightly displaced from the center of the CCD and telescope flexure may account for the movement of this point from one frame to another %[will get a private communication from someone at KPNO to confirm this.]

Spatially dependent kernel models were derived from ($r,\cos\theta$) polynomial fits to each pixel, $i,j$, in the 15x15 array of the 220 raw kernels derived for each frame. The radial variation of the PSF is by far the dominant term, but models were generally better fit with some azimuthal component, invariably with symmetry about an axis $\theta_o$. Typical models were parameterized by locating the radial point, ($x_o,y_o$), and an azimuthal symmetry axis from contour plots such as Figure 3. The general fit is given by:

$$k_{i,j} (x(r,\cos\theta),y(r,\cos\theta)) = \sum_{n,m} {a_{i,j}(n,m) ~r^{n} \cos^{m}\theta}, \eqno(4) $$ \noindent

where $r$ is the radial distance from the PSF symmetry point and $\theta$ measured with respect to the azimuthal symmetry axis. The best fits were cubic ($n=3$) in $r$ and linear ($m=1$) in $\cos\theta$. Prior to the overall fit an interactive examination of the radial plots of the central kernel pixel was performed in a small number of azimuth sections in order to delete any outlying raw kernels before the final fit. In the situation illustrated in Figures 2 and 3, with two regions defined by the ratio of the sizes of the PSF on each frame, two separate model fits were made.

\subsection{Image Subtraction}

Before subtraction both the image and reference image were processed further. Images were ``unsharp masked'' - the underlying smoothed galaxy and sky background were removed from both images by subtracting the large-scale median smoothed image of both frames. This leaves only the data to reference frame photometric scaling factor to be determined. The measures of the 220 raw kernels were used to do this. First, a 10x10 pixel box around each of the PSF stars was matched, then the intensity scaling factor was derived from a linear fit to the data versus reference pixel intensity match within this box (using software provided by A. Phillips). The scaling factor for the entire data image was the median value of the independent scaling factors determined for the 220 PSF stars, rendering it insensitive to any stellar variability, and was accurate to 1.5\%.

As a convenient method of performing later photometry, a well spaced grid of 16 of the 220 PSF stars were removed from the data frame by zeroing the pixel values in a 16x16 pixel box centered on each of these stars. Since this step is performed prior to overall PSF-matching of the data and reference image and the subsequent subtraction, this has the advantage of ensuring properly scaled, high S/N and PSF-matched reference stars from the reference image on the final subtracted image for the purposes of differential photometry.

For computational efficiency the models derived from equation (4) were used to compute kernels for a grid of $(r,\cos\theta)$ positions. In matching a frame to its reference image, the frame was divided up into an $(x,y)$ grid of subimages and the nearest kernel to the center of each subimage from the model grid was used to perform the PSF-matching convolution for that subimage. The resolution of both the model kernel grid in $r$ and $\cos\theta$ and the image processing $x,y$ grid were chosen to be at least at the point where no significant improvement in the subtraction quality could be obtained with increasing resolution. For the kernel model this was typically around $\Delta r$ of 50 pixels and $\Delta\cos\theta$ of 0.2. The final image-reference differenced subimage comprises either a data frame PSF-matched to its reference frame or a reference frame PSF-matched to the data frame, depending on the ratio of the sizes of the PSF's of the data and reference image within the subimage. In {\S 4} we assess the quality of the subtraction we have achieved with this algorithm.