High-precision Astrometry and Photometry with the JWST/MIRI Imager

The first six months after the launch of JWST (Gardner et al. 2023) focused on its Commissioning (Rigby et al. 2023) by performing a series of on-sky and internal source calibration observations and data analyses aimed at understanding the actual performance of the telescope and its instruments. The characterization of each subsystem, together with why and by how much its capabilities differed from pre-launch expectations, was aimed at providing to the community the necessary calibrations to make scientific analyses of JWST data possible from the start. The Commissioning activities covered as many aspects as possible, but due to the limited time allocated to Commissioning, the more detailed characterization of the instrument performance, and development of sophisticated tools for their analysis, was postponed to the Cycle 1.

An example is the lack of tools to obtain high-precision astrometry with JWST's imagers. The experience with the Hubble Space Telescope (HST) has shown that it is possible to achieve exquisite astrometric (and photometric) accuracies and precisions thanks to accurate effective point-spread-function (ePSF) models ¹⁰ and geometric-distortion (GD) corrections. While ePSF models were not released after Commissioning, GD corrections are available in the JWST Calibration Reference Data System (CRDS). However, these GD corrections were designed to meet the mission requirements (GD uncertainty for all detectors to be below 5 mas per coordinate; see Anderson 2016 ¹¹ ), leaving further refinements that can be critical for some scientific investigations to the Cycle-1 Calibration process.

For the cameras onboard of the HST, it took a few years to have the first tools to achieve high-precision astrometry (Anderson & King 2000). For JWST, the community has already provided initial sets of ePSFs and GD corrections for the near-infrared (NIR) detectors of the Near InfraRed Camera (NIRCam, Nardiello et al. 2022; Griggio et al. 2023) and the Near Infrared Imager and Slitless Spectrograph (NIRISS, Libralato et al. 2023). Only the Mid-InfraRed (MIR) capabilities of JWST have not been addressed yet: our paper strives to begin filling this gap in the astrometric cause of JWST.

This manuscript describes the efforts of the Mid-InfraRed Instrument (MIRI, Rieke et al. 2015; Wright et al. 2023) team, and in particular of its imager working group, to provide a set of ePSF models and GD corrections for the MIRI imager, together with the description of how they were made. Our analyses show that high-precision astrometry and photometry are also within the capabilities of the MIRI imager.

The paper is organized as follows: Section 2 presents the data used in this work; Sections 3 and 4 describe how ePSF models and GD corrections were made; Section 5 provides an overview of the tools we release; and finally Section 6 showcase our work with a simple scientific application.

The MIRI imager (MIRIM) uses a Si:As blocked impurity band conduction detector with 1032 × 1024 pixel². Not all pixels are exposed to the incoming light (Ressler et al. 2015; Morrison et al. 2023). What is commonly addressed as the imager is the largest part on the right side of the detector with a Field of View (FoV) of about 74 × 113 arcsec² (nominal pixel scale of 110 mas pixel⁻¹; Bouchet et al. 2015), but also the top-left region, which is designed for Lyot coronography, is also used for imaging, adding extra coverage in each exposure. ¹² MIRI imaging can be performed with nine broadband filters that cover a wavelength range from 5.6 to 25.5 μm (F560W, F770W, F1000W,F1130W, F1280W, F1500W, F1800W, F2100W, and F2550W; for a detailed description of MIRIM and its commissioning, see Bouchet et al. 2015 and D. Dickens et al. 2024, in preparation).

HST and NIR JWST's ePSFs were obtained with specific observations of relatively crowded fields that provide thousands of stars per image for the ePSF modeling (e.g., Anderson & King 2004, 2006; Nardiello et al. 2022; Libralato et al. 2023). At MIR wavelengths, the number of bright unsaturated sources dramatically drops, the background increases and the sensitivity drops at the longest wavelengths due to a decrease in the quantum efficiency of the detectors and significant thermal emission from the telescope. All these characteristics pushed us to try to compensate for the low statistics per exposure by including more images spanning a relatively large temporal baseline. Our MIRI ePSF models were made using publicly available data taken during Commissioning, Cycle-1 GO and calibration programs. ¹³ Only a subset of these images were used to compute the GD corrections, specifically those of the Program ID (PID) 1521. This program targeted a field in the Large Magellanic Cloud (LMC) after Commissioning (when all the major mirror alignment and phasing operations ended). This field contains a relatively high number of sources to use for computing the GD correction. The majority of the data was taken between 2022 March and December, and eight exposures from 2023 observations were included to increase the statistics in some long-wavelength filters. The number of groups, integrations and exposures differ from program to program. Overall, about the 77% of the images were taken with ≤20 groups and ≤2 integrations.

We downloaded the level-1b, uncalibrated (_uncal) products from the Mikulski Archive for Space Telescopes. ¹⁴ These exposures were processed using a development version of the JWST pipeline ¹⁵ (Bushouse et al. 2023) through the stages 1 and 2 to obtain the level-2b (_cal) fits files, fully calibrated and unresampled images (see Morrison et al. 2023, and D. Dickens et al. 2024, in preparation). The data sets were processed at different times, so the versions of the calibration pipeline ¹⁶ and of some reference files ¹⁷ were not the same but did not impact our analyses.

The MIRI ePSFs are Nyquist sampled or better at wavelengths ≳7 μm, leaving only the ePSF at 5.6 μm with F560W to be slightly undersampled. Undersampled PSFs require careful modeling to remove the so-called pixel-phase biases (positions that are systematically measured at a specific location with respect to the pixel boundaries, regardless of where objects truly are), which can be done by combining multiple well-dithered images of the same field (for a detailed description of the modeling of the undersampled ePSF of the JWST/NIRISS detector see Libralato et al. 2023).

Because the F560W MIRI ePSF is slightly undersampled (full width at half maximum of 1.88 pixels), we adopted a different approach from what is usually done for more severe cases: We collected the ePSF sampling from multiple stars in hundreds of images of various fields all at once and took advantage of the fact that dithers, the scene itself and even the variety of stellar fluxes of the sources in the field should distribute stars randomly with respect to the boundaries of the pixel. This solution is similar to how ePSFs are modeled with well-sampled, ground-based data (e.g., Anderson et al. 2006; Libralato et al. 2014, 2015; Kerber et al. 2019; Häberle et al. 2021) rather than with undersampled, space-based images (Anderson & King 2000).

A pixel (i, j) close to the center of a star with position (x_*, y_*) and flux z_* has a value P_i,j that can be defined as:

$$\begin{eqnarray*}&&{P}_{i,j}={z}_{* }\ \cdot \ {\psi }_{{\rm{E}}}(i-{x}_{* },\ j-{y}_{* })\ +\ {\mathrm{sky}}_{* },\end{eqnarray*}$$

where sky_* is the local sky background (measured in an annulus region close to the source). ψ_E(i − x_*, j − y_*) is the fraction of the light of the stars that falls on the pixel (i, j) according to the ePSF model. By inverting the equation above, we can define:

$$\begin{eqnarray*}&&{\psi }_{{\rm{E}}}(i-{x}_{* },j-{y}_{* })=\displaystyle \frac{{P}_{i,j}-{\mathrm{sky}}_{* }}{{z}_{* }}.\end{eqnarray*}$$

This simple relation tells us that every pixel in the vicinity of a star can potentially be used to sample the ePSF. Collecting many samplings from many stars in different images allows us to map the pixel-phase space and build accurate ePSF models.

For every filter, we began by measuring all bright ¹⁸ stars in each image. Positions (x_*, y_*) were initially defined as the center of light of the flux distribution in the core of each star, while fluxes z_* were estimated with aperture photometry. Once ePSF models were available, positions and fluxes were computed via ePSF fit.

Our ePSF models are oversampled by a factor of 4. Each of our ePSFs is centered at a location (151,151), is normalized to have unity flux within a radius of 10 MIRI pixels (40 oversampled pixels) and extends out to 37.5 MIRI pixels (150 oversampled pixels). Smoothness and continuity of our ePSFs (see Anderson & King 2000) were achieved by convolving the ePSFs with low-pass smoothing kernels that removed some high-frequency structures at larger radii, while preserving the shape of real features of the ePSFs. The choice of these kernels was reached with a trial-and-error approach, and tailored for each filter (the ePSF shape changes from 5.6 to 25.5 μm).

We started with a single ePSF model for the entire detector, and progressively let the ePSFs to vary spatially. We found that one ePSF that covers the Lyot region and a 2 × 3 array of ePSFs (Figure 1) for the imager were a good compromise between having enough stars in each region of the detector to model the ePSF and including spatial variations of the ePSF across the FoV. This choice allows us to interpolate the ePSF at any given pixel of the imager region, but forces fitting the same ePSF in all pixels in the Lyot region. Because the number of stars drastically drops from short to long wavelengths, the spatial variability cannot be included in all filters. Also, sometimes this same issue did not let us model the spatial variability in the wings of the ePSFs, and imposed ePSF models with a spatially variable core and constant wings. The details of the spatial variability of our ePSF are summarized as follows:

• 1.  
  F560W and F770W: one ePSF for the Lyot and 2 × 3 ePSFs for the imager. These ePSF models are spatially variable both in the core and in the wings, so to include the change of the location and intensity of the so-called cruciform artifact (e.g., Gáspár et al. 2021);
• 2.  
  F1000W, F1130W ¹⁹ and F1280W: one ePSF for the Lyot and 2 × 3 ePSFs for the imager. These ePSF models are spatially variable in the core out to 7.5 pixels (30 oversampled pixels), but constant in the wings (i.e., all ePSFs share the same profile in the wings);
• 3.  
  F1500W: one ePSF for the Lyot and 2 × 2 ePSFs for the imager (the remaining two ePSFs were linearly interpolated by means of the ePSFs at the corners of the imager). As for the previous filters, these ePSF models are spatially variable in the core out to 7.5 pixels, but constant in the wings;
• 4.  
  F1800W, F2100W, and F2550W: one ePSF model for the entire MIRI detector (Lyot and imager).

The modeling of the MIRI ePSF was made with an iterative process. At each iteration, position and flux of all sources were measured with the available ePSFs, then the ePSF samplings were redefined and the ePSFs updated. Convergence was reached between 10 (one ePSF and no spatial variability) and 20 (all ePSFs and spatial variability) iterations.

Zoom In Zoom Out Reset image size

Figure 2 shows how the ePSF samplings change for the F560W case, the only undersampled ePSF in the MIRI imager. At the first iteration, where positions were defined as the light-center of the flux distribution and fluxes were computed via aperture photometry, the shape of the ePSF was not well constrained, as stars were preferentially measured to be at the corners of the pixels (top panels). Thanks to our careful modeling, the shape of the ePSF becomes smoother and with less scatter than before. Also, stars do not show any preferred positioning in the pixel-phase space (bottom panels).

Zoom In Zoom Out Reset image size

Another way to assess the impact of the pixel-phase errors introduced by non-optimal ePSFs can be obtained by using dithered observations (Anderson & King 2000), because dithering places a source at different locations with respect to the pixel boundaries. By combining positions of stars measured in multiple dithered images, averaging them together on to the same reference-frame system, we can obtain a set of positions free from pixel-phase systematic errors, and these positions can be compared with those measured in each MIRI raw frame (e.g., Libralato et al. 2023). This test was run for the F560W filter.

We measured positions and fluxes of stars in the LMC calibration field and other adjacent fields observed during Commissioning and Cycle-1 Calibration programs. We cross-identified the same stars in multiple catalogs, ²⁰ and averaged their positions and fluxes after they were transformed on to a common reference frame system, hereafter simply master frame, by means of six-parameter linear transformations. The scale and orientation of the master frame were set up by means of the Gaia Early Data Release 3 (EDR3) catalog (Gaia Collaboration et al. 2016, 2021) projected on to a tangent plane centered at (R.A., decl.) = (80.606608, −69.461994) deg, fixing the pixel scale to be the nominal pixel scale of MIRI (110 mas pixel⁻¹). The photometry of our master frame was registered to that of a MIRI catalog.

Finally, we compared the raw positions of bright, well-measured stars in each MIRI image with those predicted by the averaged master frame inverted on to the MIRI raw frame. Figure 3 illustrates the impact of pixel-phase errors in astrometry. In the top panels, the positional residuals show a clear trend as a function of pixel phase when positions were defined as the light-center of the flux distribution. Pixel-phase biases disappear when stellar positions are fit using our ePSF (bottom panels).

Zoom In Zoom Out Reset image size

Figure 4 presents a 3D overview of one ePSF for each filter. The colored profiles refer to the ePSF of each filter, whereas the gray surfaces are the 3D profile of the F560W ePSF for reference. The ePSFs becomes progressively less sharp (from the 19.2% of the total flux in the centermost pixel of the F560W ePSF to the 1.6% of the F2550W ePSF) and broader (from the full-width half maximum of 1.88 pixels of the F560W ePSF to the 7.80 pixels of the F2550W ePSF) from short to long wavelengths.

Zoom In Zoom Out Reset image size

The python package WebbPSF (Perrin et al. 2012, 2014) can provide simulated PSFs that in principle can be used for astrometry and photometry with the JWST's imagers. However, at the moment WebbPSF does not include all detector features that contribute to the effective PSF. Appendix describes a comparison between our ePSFs and WebbPSF PSFs. Our result shows that caution is advised when using WebbPSF models, at least for now, especially for the F560W and F770W filters.

3.1. Quality of the ePSF Fit

To assess the efficacy of our ePSFs, we made use of two parameters: the "Quality of PSF fit," or QFIT, and the "radial-excess" parameter, or RADXS. The QFIT is defined as the absolute fractional error in the ePSF fit of a source (e.g., Anderson et al. 2006; Libralato et al. 2014). The closer the QFIT is to zero, the better is the fit of the ePSF. The RADXS is defined as the the excess/deficiency of flux outside the core (between 1 and 2.5 pixels) of the source compared to the prediction of the ePSF model (Bedin et al. 2008). Stars usually have a RADXS close to 0; cosmic rays or hot pixels have a negative RADXS value (i.e., they appear sharper than the ePSF model); and galaxies exhibit positive RADXS values (they are broader than the ePSF model). Thus, the RADXS is usually a powerful parameter to understand what kind of sources we are dealing with.

We run this test on all images, but discuss the result for the Commissioning program PID 1024, which observed the LMC calibration field in all filters. Figure 5 shows QFIT and RADXS as a function of instrumental magnitude for one image in each filter. The instrumental magnitude is defined as $-2.5\mathrm{log}({\rm{DNs}})$, where DNs are the total accumulated Digital Numbers in the exposure. These values were obtained by converting the pixel values in MJy sr⁻¹ (the physical unit of the level-2b files) in DN s⁻¹ using the conversion factor provided by the header keyword PHOTMJSR, and then multiplying the result by the effective exposure time available in the header keyword EFFEXPTM. Although the actual exposure time of each pixel can differ depending on the actual groups used in the ramp fit, the DN values here computed are a reasonable proxy of the signal of our sources.

Zoom In Zoom Out Reset image size

The left panels highlight the typical trend of the QFIT as a function of instrumental magnitude. The horizontal, gray, dashed line is set at QFIT = 0.05, a value generally used to select well-measured stars. The middle panels present analogous plots but for the RADXS. The horizontal line is instead set at RADXS = 0, the expected value for a point-like source. With the exception of the F2550W filter, for which this specific data set does not contain any bright star, we can measure bright point-like sources well (low QFIT and RADXS ∼ 0) thanks to our ePSF models. At faint magnitudes, the QFIT increases and the RADXS distribution becomes broader because the signal-to-noise ratio of the sources decreases, and so objects becomes progressively poorly measured, but the points are mostly centered at RADXS ∼ 0. Very-bright stars in the short-wavelength filters show an increasing QFIT and a positive trend for the RADXS. This behavior is likely due to a combination of various factors, including saturation, nonlinearity (although it is expected to be small; Morrison et al. 2023) and, most importantly, brighter-fatter (there is more flux just outside the core of the star, meaning the object is "fatter" than what the ePSF predicts; Argyriou et al. 2023, see also Section 3.3). Similar trends have been found for the NIRISS detector (Libralato et al. 2023).

The sensitivity and wavelength range covered by the MIRI imager make it possible to easily find galaxies in almost every exposure. To further help readers to understand the result of ePSF fit with MIRI data, and how it differs for point-like and extended sources, we selected a star and a galaxy which were observed in all our images in this test. The choice fell on two specific objects, i.e., the brightest sources detectable in the F2550W images. In all panels in Figure 5, the star is highlighted with a blue square, while the galaxy is depicted with a green square. The rightmost two panels, using this same color code, are zoomed in around these sources in each image.

For wavelengths ≤10 μm, the star is very bright and its QFIT and RADXS values are larger than what we would expect for a well-measured star for the reasons mentioned above. The galaxy has instead a large QFIT and a positive RADXS value, which place it in a different location with respect to where stars are.

From 11.3 to 21 μm, the star is now well below the line set at QFIT = 0.05, meaning that our ePSF models are adequate. On the other hand, the galaxy progressively displays QFIT and RADXS values closer to those of the stars. At 25 μm, the situation becomes even more extreme, with star and galaxy sharing similar RADXS values. This is likely due the larger FWHMs of the ePSFs at these wavelengths. Specifically, the longer the wavelength, the shallower and broader the ePSF becomes, and so differences in the QFIT and RADXS of point and extended sources are more subtle. For this reason, we advise to visually inspect the targets in the image as an additional cross-check.

3.2. Temporal Stability of the ePSFs

In the following, we refer to the final ePSF models described in the previous Sections as the library ePSFs. Our work made use of images taken over a temporal baseline of a few months. Although limited, this coverage allows us to monitor if, and by how much, the ePSFs change over time. The same limitations discussed for the library-ePSF modeling applies here: the low statistics makes this monitoring feasible only for the short-wavelength filters (F560W, F770W, F1000W) and using images that cover fields in the LMC, where the stellar density is moderate.

We fined-tuned our library ePSF models for every image as described in, e.g., Bellini et al. (2017) or Libralato et al. (2018): We iteratively ePSF-fitted and subtracted well-measured, bright stars in the image, and adjusted the ePSF models to minimize the residuals of the ePSF subtraction. Again, the number of stars at our disposal for this task in each image is low, so we collected the residuals over the entire FoV of MIRI and perturbed all ePSFs in the array by the same quantity (so potential spatial dependencies of the temporal variations are not considered). We find that the ePSF variations are minimal, below 4% in the worst case. These small variations are in agreement with what was found over a similar temporal baseline for NIRCam (Nardiello et al. 2022) and NIRISS (Libralato et al. 2023) imagers.

We made use of WebbPSF to correlate the variations we measured in the centermost pixel of the ePSFs with the in-flight wave front sensing measurements. We retrieved the Optical Path Difference (OPD) files, which include information about the mirrors for the period covered by the observations used to monitor the change of the ePSFs. We looked at the variations over time of the telescope+instrument wave front rms value and of the encircled energy of MIRI, and compared them with the change in the peak of our ePSFs (we chose the centermost ePSF as a benchmark). The result, presented in Figure 6, shows that the mirror configuration did not change significantly in our observations, thus confirming the lack of temporal variations of the ePSFs we measured.

Zoom In Zoom Out Reset image size

As described in Libralato et al. (2023), the perturbation process mainly affects the photometry, while astrometry is not expected to improve significantly. ²¹ Our tests show that the ePSFs are stable, as expected given the overall stability of JWST at these wavelengths. Users should always carefully explore whether the perturbation procedure actually improves the ePSFs, and should use many high signal-to-noise stars to fine-tune the ePSFs.

3.3. A Note on the Brighter-fatter Effect

Argyriou et al. (2023) have first reported evidence of Brighter-Fatter-like (BF) effects in MIRI. When a star is impacted by this effect in the MIRI imager, its profile becomes broader and less peaked than normal (see their Figure 6). We investigated whether, and to what extent, the BF effect impacts our ePSF models and the photometry resulting from their fit.

We followed a similar approach to that of Argyriou et al. (2023). Briefly, we used F560W data of the Commissioning program PID 1464. Each of the 12 images consists of a single integration and 125 groups. Such long ramps are perfect to assess the presence of the BF effect. We processed the data three times modifying the JWST pipeline so to use only the first 10 (i.e., the ramp of very bright stars is still not impacted by the BF), 50 (ramp mildly impacted by the BF) and all 125 (nonlinear deviations toward the end of the ramp due to the BF effect) groups, respectively. We did so by changing the Data Quality group flag of the remaining groups to "saturated" when running the stage-1 pipeline. This strategy allowed us to keep everything except the ramp fit as similar as possible between the three data sets. Finally, we fit our F560W ePSF models to all stars in each sample. Our ePSFs were computed with images taken with various combinations of groups and integrations, but most ramps were relatively short (Section 2). Also, we restricted the sample of stars for the modeling of the core of the ePSFs to objects not too bright, so to avoid dealing with nonlinearity and BF effects. This means that our ePSF models should be a good representation of a point-like source not impacted by these issues.

In the first test, we looked at the QFIT parameters in the three cases. Figure 7 shows the QFIT as a function of magnitude for all images processed using 10 (top panel), 50 (middle panel) and 125 (bottom panel) groups. For simplicity, the instrumental magnitudes were obtained by assuming the same total exposure time in all cases. The main difference between the three plots happens for stars with magnitude brighter than −14: the higher is the number of groups used in the ramp fit, the larger is the QFIT value. This means that the profile of these bright stars progressively deviates from what predicted by our ePSFs. Stars with magnitude >−13 show a different trend when only 10 groups are fit. This is likely due the faintness of the targets (the ramp fit improves when more photons are collected and more groups are fit).

Zoom In Zoom Out Reset image size

The second test evaluated the residuals from the ePSF subtraction of all stars, i.e., the difference between the pixel value in the image (sky subtracted) and the value predicted by the ePSF fit in the same pixel. We collected the centermost 5 × 5 pixels² of each star where the BF effect is expected to show up. Figure 8 presents the residuals for the three cases. Left panels refer to relatively faint stars (<20 DN s⁻¹). The residuals look homogeneously distributed in the pixel space in all three cases, which is what we expect for stars in a brightness regime not affected by the BF. Right panels show the residuals for very bright stars (>40 DN s⁻¹). We can clearly see that at the center of the star the ePSF predicts more flux than what is observed; conversely, the ePSF predicts less flux than what is present in the pixels just outside the core. The discrepancies between observed and predicted values increases the higher the number of groups used in the ramp fit. This means that these bright sources have a profile less peaked and broader than the ePSF models, as expected because of the presence of the BF effect.

Zoom In Zoom Out Reset image size

The last test was a comparison between the photometry obtained from each set of images. For each sample, we made a master frame (similarly to what described in Section 3). We kept only stars present in at least 10 images. We then cross-identified the same stars between the three master frames, and compared the magnitudes of the stars (Figure 9). There is a mild linear relation as a function of magnitude suggesting that the brighter is the star, the fainter it appears in the 50- or 125-group cases with respect to the 10-group case. Although other effects, like e.g., nonlinearity, ²² could contribute to the observed discrepancies, this final piece of evidence seems to point again at the BF effect as possible cause. A correction of the BF effects to apply at the ramp level will be discussed in an ongoing work (D. Gasman et al. 2024, in preparation).

Zoom In Zoom Out Reset image size

To solve for the GD of the MIRI imager, we made use of data from program 1521. As done for other instruments (e.g., Bellini & Bedin 2009; Bellini et al. 2011; Libralato et al. 2014; Häberle et al. 2021; Griggio et al. 2023; Libralato et al. 2023), our GD correction is made up by two parts: a polynomial solution and a look-up table of residuals. We computed the GD corrections for the F560W, F770W, and F1000W filters. We postpone the analysis of the other filters to future releases.

The GD corrections for the MIRI imager were computed using the Gaia EDR3 catalog as a reference. The Gaia positions were propagated at epoch mid 2022 (a representative date for our observations) to remove the contribution of the internal motions of the LMC stars in the GD computation. Then, we projected positions on to a tangent plane centered at (R.A., decl.) = (80.606608, −69.461994) deg. The pixel scale was set to be the nominal pixel scale of MIRI (110 mas pixel⁻¹), and the x and y axes were oriented West and North, respectively.

We start by measuring stellar positions in the MIRI images by fitting our library ePSF models. Only bright, well-measured, unsaturated stars with a QFIT lower than 0.1 were kept in each MIRI catalog. Selections were also applied to the Gaia catalog: we considered in the analysis only stars with G magnitude between 13 (saturation makes Gaia astrometry worse for stars brighter than this threshold) and 19 (to exclude faint stars with poorly measured proper motions) with a positional error (including the contribution for the proper-motion propagation to epoch mid 2022) lower than 0.01 MIRIM pixel.

We cross-identified the same stars in the MIRI and Gaia catalogs, and transformed the Gaia positions on to the raw reference system of the MIRI imager by means of four-parameter linear transformations (rigid shifts in the two coordinates, one rotation, and one change of scale). Positional residuals were defined as the difference between these transformed Gaia and the raw MIRI positions. We collected all residuals and fit them with two fifth-order polynomial functions. The coefficients of the polynomial functions were obtained via least-square fit of all positional residuals. For the polynomial correction, we chose the center of the imager (x_ref, y_ref) = (693.5, 512.5) as the reference pixel with respect to which solve for the GD (Libralato et al. 2023). This pixel is also the reference pixel of the MIRI imager full-frame array used by the Science Instrument Aperture Files of JWST. ²³ The GD correction for each star is defined as:

$$\begin{eqnarray*}\left\{\begin{array}{rcl}\delta x & = & \displaystyle \sum _{i=1,5}\displaystyle \sum _{j=1,5-i}{a}_{{ij}}{\tilde{x}}^{i}{\tilde{y}}^{j}\\ \delta y & = & \displaystyle \sum _{i=1,5}\displaystyle \sum _{j=1,5-i}{b}_{{ij}}{\tilde{x}}^{i}{\tilde{y}}^{j}\end{array}\right.,\end{eqnarray*}$$

where:

$$\begin{eqnarray*}&&\left\{\begin{array}{l}\tilde{x}=\displaystyle \frac{x-{x}_{{\rm{ref}}}}{{x}_{{\rm{ref}}}}\\ \tilde{y}=\displaystyle \frac{y-{y}_{{\rm{ref}}}}{{y}_{{\rm{ref}}}}\end{array}\right.\end{eqnarray*}$$

We set to 0 the coefficients a₀ and a₁ of the polynomial in order to: (i) have the solution at the reference pixel to have its scale equal to that of the chip, and (ii) the corrected y axis aligned to the raw y axis. The terms b₀ and b₁ were instead left to assume whatever values fit best because the detector axes have different scale and not be perpendicular to each other. For the polynomial part, we performed 500 iterations, each time adding only 75% of the coefficient values to the previous estimates to ensure a convergence of the solution. Figure 10 shows the distortion maps and the positional residuals as a function of x/y coordinates before (left) and after (right) applying the polynomial correction in the case of the F560W filter. This first correction improves the positional residuals by a factor ∼50, but we can still notice part of the nonlinear GD is left uncorrected. The result from the polynomial correction provided in this work is similar to what can be obtained using the GD correction in the corresponding CRDS reference file, which, as well, includes a polynomial correction only.

Zoom In Zoom Out Reset image size

The residual GD was empirically corrected using a look-up table of residuals. The positional residuals after the polynomial correction for the F560W, F770W, and F1000W filters had the same shape and amplitude, ²⁴ and therefore, we chose to work with all of them at once to increase the statistics, at the cost of a less accurate distortion than with a purely filter-based solution. We also relaxed the quality selections on both the MIRI and Gaia catalogs, again to have enough stars to tabulate the correction.

We collected all positional residuals and averaged them in two look-up tables, a 3 × 3 for the Lyot region (141 × 92.7 pixel² elements) and a 8 × 12 (127.9 × 128.1 pixel² elements) for the imager. The grid points adjoining to the edges of the Lyot or the imager regions were displaced to edges of the cell, so to allow the correction to be always computed by bi-linearly interpolating between grid elements. Again, the computation of this second correction was iterated 100 times, each time adding 50% of the correction to ensure convergence. The distortion map and positional residuals after all corrections are applied for all three filters is shown in Figure 11. Figure 12 presents the positional residuals as a function of x and y positions for each filter separately. The residual GD systematics are within 0.01 MIRIM pixel and are larger close to the edges of the detector, where we expect our GD solution to be less accurate because of the lack of sources in the region to model the GD.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Astrometric Precision and Short-term Temporal Stability

We evaluated the astrometric precision of our GD correction by combining multiple dithered MIRI observations from three distinct data sets covering fields in the LMC, specifically PIDs 1024, 1040, and 1521. This allowed us to test our GD using images not related to the GD-computation process and, at the same time, to monitor the short-term stability of the GD since these three sets of observations were taken three months apart.

For each program/filter, we made a master frame as described in Section 3 using six-parameter linear transformations. We only kept stars measured in at least three images. The positional rms as a function of the instrumental magnitude for all cases is shown in Figure 13. The median value of the one-dimensional positional residuals is always of the order of 0.01 pixel or better, in agreement with what found for NIRCam and NIRISS (Griggio et al. 2023; Libralato et al. 2023). The results for PID 1024 are a factor two better than in the other cases. However, the PID 1024 images were taken with a small dither and the estimate of the astrometric precision with this data set is blind to systematic residuals (a given star falls in the same region the detector in different images, and the GD usually varies smoothly on small scales). Images of PIDs 1040 and 1521 were instead taken in a 3 × 3 mosaic with large offsets between tiles and can provide a more reliable measurement of the astrometric precision. Finally, we can also conclude that the nonlinear terms of the GD are rather stable over a period of three months.

Zoom In Zoom Out Reset image size

4.2. Pixel Scale

We estimated the absolute pixel scale of the MIRI imager by means of the Gaia DR3 catalog. We cross-identified bright and well-measured stars between each MIRI catalog of PID 1521 and the Gaia catalog, and used six-parameter linear transformations to transform positions from one frame to the other. The relative scale between the two frames is one of the linear parameters solved by these transformations. Given that the pixel scale of our master frame was defined to be 110 mas pixel⁻¹ (see Section 3), we could derive the absolute pixel scale of the MIRI imager from the relative one.

We compute the absolute scale for each filter separately, taking into account for the velocity-aberration correction factor (provided by the header keyword VA_SCALE). The median values of the pixel scale in each filter analyzed for the GD correction is provided in Table 1. We find that the pixel scale is filter dependent.

Table 1. Average Pixel Scale for Three MIRI Filters (corrected for velocity aberration)

Filter	Pixel scale (mas pixel−1)	Error (mas pixel−1)
F560W	110.47645	0.00029
F770W	110.47894	0.00036
F1000W	110.40550	0.00034

Download table as:  ASCIITypeset image

We analyzed the variation of the pixel scale as a function of time by repeating the same process described before using also data from PID 1024 and PID 1040 (Figure 14). We find the pixel scale to be rather stable over three months, with variations at the 5σ level in the worst case. We also noticed that while the pixel scale is almost constant in all images of PID 1024, the values show an apparent trend as a function of time during the observations of PIDs 1040 and 1521 in all filters. The variations are of the order of 4 × 10⁻⁵ in ∼13 hr of observation. This resembles the pattern of the pixel-scale variations seen in the NIRISS detector by Libralato et al. (2023). Interestingly, both the NIRISS and MIRI images showing this modulation in the pixel scale have been taken as part of a 3 × 3 mosaic with large dithers.

Zoom In Zoom Out Reset image size

We release our ePSF models and GD corrections to the community through various channels. We describe here how to use them and provide a few caveats for the users.

We provide two sets of ePSF models. ²⁵ In addition to the 301 × 301 pixel² ePSFs described in Section 3, we also make available a cutout version with size 101 × 101 oversampled pixels², with the ePSF centered at pixel (51,51). Besides the size and the center, the two sets of ePSFs have the same conventions (the ePSFs are normalized to have unity flux within a radius of 10 MIRI pixels—40 oversampled pixels).

These smaller ePSFs are specifically designed to be used with the FORTRAN routine jwst1pass ²⁶ (Libralato et al. 2023; J. Anderson et al. 2024 in preparation), which is the JWST analog of the tool hst1pass Anderson (2022a, 2022b) designed for ePSF fitting with HST images. The larger ePSFs are provided to the community for completeness, but can only be used with other tools, at least for now. Various PSF-fitting software packages like photutils (Bradley et al. 2022) can potentially use our ePSFs as far as these tools follow the ePSF conventions and caveats we described in Section 3.

The ePSF file for each filter contains nine ePSF slots, but only eight are filled. Six (2 × 3) ePSFs are reserved for the imager. They can be linearly interpolated to construct the ePSF for a star centered in any pixel of the imager region. The ePSF made for the Lyot region is copied into the center left and upper left slots of the PSF array, which allows the PSF to be constant within the Lyot region. This formulation allows jwst1pass to use the generic STDPSF-based software developed for all HST and JWST images with MIRI. Specific header keywords are provided in each file with the fiducial detector location for each of the nine PSFs (see Anderson 2022a, 2022b). The ePSFs for the 4QPMs are not provided since this region is not calibrated for direct (as opposed to coronagraphic) imaging.

The tool jwst1pass uses the centermost 5 × 5 pixel² of a star to fit the ePSF model. We have noticed that at long wavelengths (F2100W and F2550W filters) faint stars are not always detected. These sources do not have a well-defined local maximum because of the high pixel-to-pixel noise. Finally, we stress again that the ePSFs at these wavelengths, especially for the F2550W filter, were obtained with few stars, and they are less accurate. For all these reasons, we advise users to carefully interpret the result of the ePSF fit and evaluate it on a case-by-case basis using the various diagnostic parameters and images output by jwst1pass.

The GD corrections are released as data cube FITS file ²⁷ designed to work with the jwst1pass software. We refer to Libralato et al. (2023) and Appendix G of Anderson (2022a, 2022b) for the description of the GD format. We also release ²⁸ a python tool that applies our GD corrections to a list of coordinates. Note that positions must be defined in a 1-index reference frame. Positions in a 0-index python-like frame would need to offset by one pixel in each coordinate before applying the GD correction.

We tested our tools by studying a field in the LMC observed as part of the Commissioning PID 1171 in F560W and F770W filters. We chose this data set because it was employed neither for the ePSF modeling nor for the GD correction, and it thus offers us an independent benchmark.

We fitted all sources in these images using our ePSF models. Positions were then corrected for GD using our corrections. For each filter, we made a master frame as described in Section 3, using the Gaia catalog to setup orientation and scale. Finally, we identified the stars in common between the F560W and F770W star lists. The color–magnitude diagram (CMD) of the stars in this LMC field is shown in the right panel of Figure 15. The left panel presents the CMD obtained using aperture photometry on the LVL-3 mosaic image with the current JWST pipeline. Our photometry was registered on to the VEGA-mag system by computing the zero-point between our photometry and that of the pipeline using bright, well-measured stars in common. PSF and aperture photometry should be comparable in this field because it is not very crowded. However, the MIRI ePSF photometry provides a narrower sequence in the CMD at all magnitudes, not only at the faint end as expected (Libralato et al. 2016a, 2016b).

Zoom In Zoom Out Reset image size

The accuracy of our GD correction was tested by computing proper motions (PMs) for the stars in this field. We cross-identified the stars in our MIRI F560W catalog with those present in the Gaia DR3 catalog and transformed (with global, six-parameter linear transformations) the positions of the stars in the MIRI master frame on to the reference frame of the Gaia catalog. We used only stars belonging to the old LMC population along the red-giant branch in the CMD (left panel of Figure 16) to compute the coefficients of the transformations between frames. Thus, our PMs are initially relative to the bulk motion of these old stars. PMs are defined as the positional displacements divided by the temporal baseline (∼8 yr) and multiplied by the pixel scale of the Gaia catalog (110 mas pixel⁻¹). PM errors are computed as the sum in quadrature of the positional errors, again divided by the temporal baseline and multiplied by the pixel scale. Finally, we linked our relative PMs on to the absolute frame of Gaia by computing the PM zero-point between our and Gaia's PMs. We provide the G versus (G − m_F560W) CMD and the PM error as a function of G in Figure 16. The comparisons between the vector-point diagram of the stars in this LMC field using our (black points) and Gaia's (red points) PMs in two magnitude regimes are also shown. For the brightest stars in common between the two catalogs, the MIRI+Gaia PMs are comparable, with a median PM error of ∼80 μas yr⁻¹. At the faint end, our PMs are more precise because of the long temporal baseline and worse quality of the Gaia PMs. It is worth noticing that these results were obtained with simple techniques, and we expect a more careful analysis (e.g., Libralato et al. 2022) to provide more precise PMs.

Zoom In Zoom Out Reset image size

We presented a careful analysis of the astro-photometric capabilities of the MIRI imager. We described in detail how accurate ePSF models and GD corrections were made, and explain how to use them to obtain PSF-based photometry with JWST data, which is currently not performed by any step of the JWST imager pipeline.

We made ePSF models for all MIRI filters using public Commissioning, Cycle-1 Calibration and GO data. We refer the reader to the discussions throughout Sections 3 and 5 for how to use these models and for the associated caveats. We analyzed the impact of the brighter-fatter effect in MIRI PSF photometry. This, combined with other effects like nonlinearity, can impact the photometry at a few-percent level and result in imperfect PSF subtraction. For the GD, we instead release only the corrections of the three shortest-wavelength filters among those available for MIRI. We show that GD is stable, at least in the short timescale analyzed in our work.

We tested our ePSFs and GD corrections by analyzing a stellar field in the LMC. We showed that PSF photometry outperforms what can be obtained using aperture photometry. We also computed PMs by combining JWST MIRI and Gaia catalogs, which can be a powerful combination to analyze the kinematics of stars in the faint regime of Gaia (e.g., del Pino et al. 2022).

These simple applications highlight the astrometric and photometric potential of the MIRI imager. The complete assessment of the astro-photometric capabilities of this imager is still to be completely understood; only more data spanning a large variety of exposure parameters and the impact of the pipeline and calibration will enable us to learn more about this instrument. The MIRI resolution is worse than that of the NIR JWST detectors but the wavelength coverage of MIRI is unique among JWST's detectors, making the MIRI imager an invaluable tool at disposal of the astronomical community.

M.L. thanks Marshall Perrin, Marcio Melendez Hernandez and Andrea Bellini for the insightful discussions about WebbPSF, ePSFs and GD. M.L. also thanks Jay Anderson for implementing the MIRI module in jwst1pass and for hosting the ePSF and GD-correction files on his web page. I.A. would like to thank the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) for their support in the framework of the PRODEX Programme. P.J.K. acknowledge support from the Science Foundation Ireland/Irish Research Council Pathway program under grant No. 21/PATH-S/9360. Based on observations with the NASA/ESA/CSA JWST, obtained at the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-03127. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of astropy, a community-developed core python package for Astronomy (Astropy Collaboration et al. 2013, 2018).

Facility: JWST - (MIRI), GAIA.

Software: FORTRAN, python, astropy, photutils.

Simulated PSFs for the JWST instruments are available thanks to the python package WebbPSF. In the version used for our test (1.1.2, developer version), WebbPSF PSFs takes into account models of the telescope and optical state of an instrument but does not include detector effects like interpixel capacitance and intrapixel sensitivity variations. The only detector effect included for MIRI is the cruciform artifact, although the simulated cruciform properties are significantly different from what is observed in flight data (see below). This qualitative comparison is meant to highlight the importance of the detector effects.

We simulated one PSF with WebbPSF for each filter and compared it with one of our ePSFs. For filters where the spatial variability was included in the modeling, we considered the centemost ePSF in our array. Figure A1 presents the encircled-flux profiles obtained from the WebbPSF (red lines) and our ePSF (blue lines) models. These profiles are normalized to have unity flux within a radius of 10 MIRIM pixels. Insets in each panel show the centermost part of both models. These profiles and insets suggest that both models are rather similar, although the first minimum of our ePSFs predict more flux than WebbPSF. This is likely due to the smoothing functions used in the ePSF modeling (see Section 3), and to the lack of detector effects in WebbPSF PSFs. The F1130W ePSFs predict more flux than those of WebbPSF at all radii. Our F1130W ePSFs were obtained from data where stars are embedded in dust, which caused a gradient in the wings of the ePSFs in this filter. We choose to include the wings anyway, but we advise caution when using them.

Zoom In Zoom Out Reset image size

The F560W and F770W PSFs include the cruciform artifact. As a benchmark, we focused on the F560W filter and simulated seven PSFs located at the same reference points of our ePSF models. We provide a detailed comparison for the F560W case in Figure A2. The three leftmost panels show, respectively, the WebbPSF PSF, a zoom-in in its core and the difference with respect to the average PSF in the FoV. The rightmost three panels present analogous plots using our ePSFs.

Zoom In Zoom Out Reset image size

The core of both sets of PSFs is rather similar. ²⁹ The main differences are present in the wings of the PSF, where the cruciform artifact is prominent. WebbPSF PSFs seem to predict more flux in the cruciform artifact than our ePSFs that does not vary significantly across the FoV. Instead, our ePSFs show that the cruciform artifact is weaker than in the WebbPSF models and a more evident spatial dependence, especially for the horizontal component. Also, the cruciform seems to be bent (Gáspár et al. 2021).

In light of these pieces of evidence, we advise caution when using WebbPSF models, at least for now. Despite providing an accurate description of JWST and its instruments' optical model, systematic effects might prevent users from achieving high-precision astrometry and photometry with WebbPSF PSFs, because detector effects like interpixel capacitance, intrapixel sensitivity and, for the F560W and F770W filters, the cruciform artifacts is not accurately modeled. Preliminary tests made with the latest WebbPSF version (1.2.1), which includes the effect of the interpixel capacitance of MIRI (M. Engesser et al. 2024, in preparation) and a cruciform model closer to that of flight data, show a better agreement with our ePSF models.