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

Astrometry is one of the main pillars of astronomy, and one of its oldest branches. Over the years, an increasing number of astrometric works by means of Hubble Space Telescope (HST) data have revolutionized our understanding of various phenomena. With the launch of JWST, it becomes almost instinctive to want to replicate or improve these results with data taken with the newest, state-of-the-art, space-based telescope. In this regard, the initial focus of the community has been on the Near-Infrared detectors on board of JWST because of their high spatial resolution. This paper begins the effort to capture and apply what has been learned from HST to the Mid-InfraRed Instrument (MIRI) of JWST by developing the tools to obtain high-precision astrometry and photometry with its imager. We describe in detail how to create accurate effective point-spread-function (ePSF) models and geometric-distortion corrections, analyze their temporal stability, and test their quality to the extent of what is currently possible with the available data in the JWST MAST archive. We show that careful data reduction provides deep insight on the performance and intricacies of the MIRI imager, and of JWST in general. In an effort to help the community devise new observing programs, we make our ePSF models and geometric-distortion corrections publicly available.


INTRODUCTION
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 Corresponding author: Mattia Libralato libra@stsci.edutelescope 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 de-velopment of sophisticated tools for their analysis, was postponed to the Cycle 1.
An example is the lack of tools to obtain highprecision 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) models1 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 2 ), 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 InfaRed 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 highprecision astrometry and photometry are also within the capabilities of the MIRI imager.
The paper is organized as follows: Sect. 2 presents the data used in this work; Sects.3 and 4 describe how ePSF models and GD corrections were made; Sect. 5 provides an overview of the tools we release; and finally Sect.6 showcase our work with a simple scientific application.

INSTRUMENT DESCRIPTION AND DATA SETS
The MIRI imager (MIRIM) uses a Si:As blocked impurity band conduction detector with 1032×1024 pixel 2 .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 2 (nominal pixel scale of 110 mas pixel −1 ; 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 exposure3 .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 Dicken et al., 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 programs4 .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 March and December 2022, 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 (MAST) 5 .These exposures were processed using a de-Figure 1. MIRI ePSF layout using the F560W filter data as an example.The image in the background is a cal fits file of PID 1521.The position of each ePSF is indicative of the region within which the stars for the given ePSF modeling were selected.Slots (1,1) and (1,2) are not shown because the bottom-left region of the imager is reserved for the 4-quadrant phase mask coronagraphs, thus making it not available for direct imaging.
velopment version of the JWST pipeline 6 (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 Dicken et al., in preparation).The data sets were processed at different times, so the versions of the calibration pipeline 7 and of some reference files 8 were not the same but did not impact our analyses.

EPSF MODELING
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 6 https://github.com/spacetelescope/jwst. 7The MIRI data was processed using either one of the following pipeline versions: 1.9.5.dev29+g40f282e, 1.9.5.dev26+g4285c4e, 1.10.2.dev11+g3f269f5.the pixel boundaries, regardless of where objects truly are), which can be done by combining multiple welldithered 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 modelled with well-sampled, ground-based data (e.g., Anderson et al. 2006;Libralato et al. 2014Libralato et al. , 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: where sky * is the local sky background (measured in an annulus region close to the source).
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: 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 bright9 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 (Fig. 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: • 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. 2021a); • F1000W, F1130W 10 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); • 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; • F1800W, F2100W and F2550W: one ePSF model for the entire MIRI detector (Lyot and imager).
10 Most of the stars used to model the wings of the F1130W ePSFs are embedded in dust, which caused a gradient in the wings of the ePSFs in this filter.We advise caution when using the wings of these particular ePSFs.
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.
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).
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 pixelphase 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 catalogs11 , and averaged their positions and fluxes after they were transformed on to a common reference frame system, hereafter simply master frame, by means of sixparameter 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.,Dec.)= (80.606608,−69.461994) deg, fixing the pixel scale to be the nominal pixel scale of MIRI (110 mas pixel −1 ).The photometry of our master frame was registered to that of a MIRI catalog.
Figure 2. ePSF samplings for the F560W filter at the beginning (top row) and at the end (bottom row) of our iterative ePSF-modeling process.The first two columns from the left show the ePSF samplings (ψE) as a function of distance from the center of the ePSF along the x and y axes, respectively.Only 1% of the points in a strip within ±0.1 pixel from the center of the stars are plotted.The density plots in the third column (color-coded as in the colorbar on the right) highlight where the measured centers of the stars are with respect to the boundaries of the pixels.At the beginning of the iterative process, the ePSF samplings are preferentially measured at the corners of the pixels, and the ePSF model is too sharp and with many outliers.Once accurate ePSF models are used to find the center of the stars, the ePSF samplings are more homogeneously distributed with respect to the pixel boundaries, and the shape of the ePSF is smooth and well constrained.Finally, we compared the raw positions of bright, wellmeasured 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).
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.
The python package WebbPSF (Perrin et al. 2012(Perrin et al. , 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 A describes a compar- ison 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.

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 instrumen-tal magnitude for one image in each filter.The instrumental magnitude is defined as −2.5 log(DNs), where DNs are the total accumulated Digital Numbers in the exposure.These values were obtained by converting the pixel values in MJy sr −1 (the physical unit of the level-2b files) in DN s −1 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.
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 pointlike 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, non-linearity (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 Sect.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 Fig. 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.

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 shortwavelength 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 wavefront 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 wavefront rms value and of the encircled energy (EE) 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 Fig. 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.
As described in Libralato et al. (2023), the perturbation process mainly affects the photometry, while as-trometry is not expected to improve significantly 12 .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 12 Some changes in the reference files and pipeline used in the processing of the data can also require a perturbation of the library ePSFs.However, turning on/off some pipeline steps like the interpixel-capacitance correction can significantly impact how the flux of point sources is distributed across pixels and could require new ePSF models.should use many high signal-to-noise stars to fine-tune the ePSFs.

A note on the brighter-fatter effect
Argyriou et al. ( 2023) have firstly 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 Fig. 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 (non-linear deviations towards the end of the ramp due to the BF effect) groups, respectively.We did so by changing the Data Quality (DQ) 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 (Sect.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 non-linearity 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).
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 2 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 −1 ).The residuals look homogeneously distributed in the pixel space in all three cases, The left panels shows the residuals from the ePSF subtraction in the centermost 5×5 pixel 2 for all stars in the faint regime, whereas the right panels focus on very bright sources that are likely impacted by the BF effect.From top to bottom, each row corresponds to the result obtained using images processed using only the first 10, 50 or all 125 groups in the ramp fit, respectively.
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 −1 ).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.
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 Sect.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 (Fig. 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.non-linearity 13 , 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 (Gasman et al., in preparation).

GEOMETRIC-DISTORTION CORRECTION
To solve for the GD of the MIRI imager, we made use of data from program 1521.As done for other instruments (e.g., Libralato et al. 2014;Bellini & Bedin 2009;Bellini et al. 2011;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 13 The uncertainty on the linearization of the ramps is larger in the presence of the BF effect, resulting in a larger uncertainty in the flux measured by each pixel 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.,Dec.)= (80.606608,−69.461994) deg.The pixel scale was set to be the nominal pixel scale of MIRI (110 mas pixel −1 ), 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 (SIAF) of JWST 14 .The GD correction for each star is defined as:  We set to 0 the coefficients a 0 and a 1 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 0 and b 1 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 non-linear 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.
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 15 , 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 2 elements) and a 8×12 (127.9×128.1 pixel 2 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 Fig. 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 15 It is worth noticing that the ∆y as a function of y plot in the rightmost panel of Fig. 10 shows a periodic pattern with period of about 374 MIRIM pixels.The pattern is present regardless of the filter, as we would expect for a detector-related feature (e.g.Libralato et al. 2014Libralato et al. , 2023)).
Figure 11.As in Fig. 10, but after the look-up table of residuals correction is also applied.Vectors are now magnified by a factor 5 000.
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.

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 Sect. 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 Fig. 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 dif- ferent 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 non-linear terms of the GD are rather stable over a period of three months.

Pixel scale
We estimated the absolute pixel scale of the MIRI imager by means of the Gaia DR3 catalog.We crossidentified 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 −1 (see Sect. 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.
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 (Fig. 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 −5 in ∼13 hours of observation.This resembles the pattern of the pixelscale 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.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.

PUBLIC RELEASE OF THE MIRI EPSF MODELS AND GD CORRECTIONS
We provide two sets of ePSF models 16 .In addition to the 301 × 301 pixel 2 ePSFs described in Sect.3, we also make available a cutout version with size 101 × 101 oversampled pixels 2 , 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 17 (Libralato et al. 2023, Anderson et al., in preparation), which is the JWST analog of the tool hst1pass Ander- 16 The MIRI ePSF models are available at STDPSFs. 17 The FORTRAN code is available at CODE. son (2022a,b) 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 Sect.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,b).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 2 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 pixelto-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 file18 designed to work with the jwst1pass software.We refer to Libralato et al. (2023) and Appendix G of Anderson (2022a,b) for the description of the GD format.We also release19 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.

EXAMPLE OF SCIENTIFIC APPLICATION
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 Sect.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 Fig. 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,b).
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 Fig. 16) to compute the coefficients of the transfor- mations 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 years) and multiplied by the pixel scale of the Gaia catalog (110 mas pixel −1 ).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 Fig. 16.The comparisons between the vector-point diagram (VPD) 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 −1 .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.

CONCLUSIONS
We presented a careful analysis of the astrophotometric 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 Sects.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 non-linearity, 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.

ACKNOWLEDGMENTS APPENDIX
A. COMPARISON WITH WEBBPSF MODELS 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 Sect. 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.
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 Fig. 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 field of view.The rightmost three panels present analogous plots using our ePSFs.
The core of both sets of PSFs is rather similar 20 .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 field of view.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. 2021b).
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 modelled.Preliminary tests made with the latest WebbPSF version (1.2.1), which includes the effect of the interpixel capacitance of MIRI (Engesser et al., in preparation) and a cruciform model closer to that of flight data, show a better agreement with our ePSF models.

Figure 3 .
Figure 3. Pixel-phase errors for the ∆x and ∆y positional residuals.The top row refers to the first iteration of our ePSF modeling when the center of the stars was defined as the light-center of the flux distribution.The bottom row shows the results for the last iteration when the ePSF models were fit to compute the position of the stars.The same stars are shown in all panels.

Figure 4 .
Figure 4. 3D view of our ePSF models.In each panel, the colored surface represent the ePSF in the filter reported on top, while the gray surface is the profile of the F560W ePSF as a reference.At the top of each panel, we provide the percentage of flux of the star within the centermost pixel and the full width at half maximum (FWHM) in pixel.

Figure 5 .
Figure 5. Quality of the ePSF fit in each filter for a representative image of PID 1024.The left panels show the QFIT as a function of instrumental magnitude.The horizontal, dashed, gray line is set at QFIT = 0.05.The middle panels provide analogous trends but for the RADXS.The horizontal, dashed, gray line is set at RADXS = 0.The blue and green squares highlight the position of a representative star and galaxy in each plot, respectively.The rightmost panels show these two objects in each image.Structures in the lowermost right panels are due to the high noise in the image.White pixels are those flagged by the pipeline as not to use.

Figure 6 .
Figure 6.In the top panel, we show the rms of the wave front error (WFE) of the observatory as a function of time (see alsoMcElwain et al. 2023).The black dots mark the mirror sensing visits.The red arrows highlights when the primary mirror segments assemblies were corrected, because the rms WFE was above the correction threshold (light-blue, dashed, horizontal line).The variations in the Encircled Energy (EE) of the F560W MIRI WebbPSF PSFs retrieved from the OPD files is presented in the middle panel.The red, dashed, horizontal line is set at 0 as a reference.Finally, the change of the fraction of light that falls in the centermost pixel of our F560W ePSF at the center of the FoV is plotted in bottom panel (again, a reference red, dashed, horizontal line is set at 0).

Figure 7 .
Figure7.QFIT as a function of instrumental magnitude for all PID-1464 images processed using the first 10 groups (top panel), 50 groups (middle panel) and all 125 groups (bottom panel).The red, horizontal, dashed line is set at 0.05, while the vertical one is set at instrumental magnitude −14.These two lines are set for reference.

Figure 8 .
Figure 8. Overview of the BF effect in the fit of our ePSF.The left panels shows the residuals from the ePSF subtraction in the centermost 5×5 pixel 2 for all stars in the faint regime, whereas the right panels focus on very bright sources that are likely impacted by the BF effect.From top to bottom, each row corresponds to the result obtained using images processed using only the first 10, 50 or all 125 groups in the ramp fit, respectively.

Figure 9 .
Figure 9. Impact of the BF effect on PSF photometry.The top panel shows the magnitude difference between the instrumental magnitude measured from the master frame made with images with 10 and 125 groups used in the fit, while the bottom panel refers to the comparison between the 10and 50-group cases.Error bars are defined as the sum in quadrature of the magnitude rms of the stars in each sample.The red lines are a weighted least-square straight-line fit to the point (the equation is reported in red).
x−x ref x ref ỹ = y−y ref y ref .14See the corresponding JDox page and references therein.

Figure 10 .
Figure 10.GD maps before (left panel) and after (right panel) applying the polynomial correction for the F560W data.Vectors are magnified by a factor of 50 (left) and 5 000 (right) to enhance the details.The positional x and y positional residuals as a function of x and y raw MIRI positions are shown in the side panels.

Figure 12 .
Figure 12.The x and y positional residuals as a function of x and y raw MIRI positions for each of the three filters used to compute the table of residuals correction.The red, dashed horizontal line is set to 0 as a reference.The median and the 68.7-th percentile about the median value of the residuals in 250-pixel-wide bins are shown in red.

Figure 13 .
Figure13.One-dimensional RMS (the sum in quadrature of the positional RMS along the x and y axes divided by √ 2) as a function of instrumental magnitude for the F560W, F770W and F1000W filters from top to bottom, respectively.The first column from the left refers to the results obtained using data from PID 1024, the second from PID 1040 and the third from PID 1521.The red horizontal line is set to the median value of the one-dimensional RMS of bright, unsaturated stars with a QFIT lower than 0.05.The pixel scale used for the pixel-to-mas conversion is the nominal pixel scale of MIRI (110 mas pixel −1 ).

Figure 14 .
Figure14.Pixel scale as a function of time.Blue, azure and red points refer to the results from the analysis of F560W, F770W and F1000W filters data, respectively.The dashed, horizontal lines, color-coded as for the points, mark the average pixel scale of PID 1521 data.

Figure 15 .
Figure 15.Calibrated mF560W versus (mF560W − mF770W) CMD for the stars in the LMC field imaged by PID 1171.The left panel shows the result obtained with aperture photometry on the LVL-3 mosaic image with the current JWST pipeline, while the right panel refers to the ePSF-based photometry on the LVL-2 images.

Figure 16 .
Figure 16.Overview of our PM test using the MIRI data of PID 1171 and the Gaia catalog.From left to right, we show G versus (G − mF560W) color-magnitude diagram (CMD); the PM error as a function of G; the vector-point diagram (VPD) of the stars with our PMs; the VPD of the same stars with Gaia's PMs.In all but the leftmost panel, black points refer to our PM catalog, whereas the red points marks the analog quantities in the Gaia catalog.

Table 1 .
Average pixel scale for three MIRI filters (corrected for velocity aberration).