JWST-TST Proper Motions: I. High-Precision NIRISS Calibration and Large Magellanic Cloud Kinematics

TL;DR

This paper presents high-precision calibration methods for JWST NIRISS, enabling detailed proper motion studies of stars in the Large Magellanic Cloud with accuracy surpassing Gaia for faint stars, and demonstrates JWST's potential for advanced astrometry.

Contribution

We developed and validated calibration techniques for JWST NIRISS that enable high-precision astrometry, and applied these to measure proper motions in the LMC, including the first absolute PM of a star cluster.

Findings

Achieved median PM uncertainty of ~13 μas/yr for G~20 stars.

First measurement of the absolute proper motion of the star cluster OGLE-CL LMC 407.

Detected kinematic substructure and measured velocity dispersions consistent with or refining Gaia results.

Abstract

We develop and disseminate effective point-spread functions and geometric-distortion solutions for high-precision astrometry and photometry with the JWST NIRISS instrument. We correct field dependencies and detector effects, and assess the quality and the temporal stability of the calibrations. As a scientific application and validation, we study the proper motion (PM) kinematics of stars in the JWST calibration field near the Large Magellanic Cloud (LMC) center, comparing to a first-epoch Hubble Space Telescope (HST) archival catalog with a 16-yr baseline. For stars with G~20, the median PM uncertainty is ~13 $\mu$as yr$^{-1}$ (3.1 km s$^{-1}$), better than Gaia DR3 typically achieves for its very best-measured stars. We kinematically detect the known star cluster OGLE-CL LMC 407, measure its absolute PM for the first time, and show how this differs from other LMC populations. The inferred cluster dispersion sets an upper limit of 24 $\mu$as yr$^{-1}$ (5.6 km s$^{-1}$) on systematic uncertainties. Red-giant-branch stars have a velocity dispersion of 33.8 $\pm$ 0.6 km s$^{-1}$, while younger blue populations have a narrower velocity distribution, but with a significant kinematical substructure. We discuss how this relates to the larger velocity dispersions inferred from Gaia DR3. These results establish JWST as capable of state-of-the-art astrometry, building on the extensive legacy of HST. This is the first paper in a series by our JWST Telescope Scientist Team (TST), in which we will use Guaranteed Time Observations to study the PM kinematics of various stellar systems in the Local Group.

Figures

Tables

Table 1. Table 1: Average pixel scale of each filter corrected for velocity aberration. The filter-averaged scale of NIRISS that we measure is 65.562 ± plus-or-minus \pm 0.003 mas pixel -1 .

Filter	Pixel scale [mas pixel-1]	Error [mas pixel-1]
F090W	65.56196	0.00003
F115W	65.55825	0.00004
F140M	65.56039	0.00004
F150W	65.54495	0.00004
F158M	65.53851	0.00003
F200W	65.55927	0.00004
F277W	65.57951	0.00004
F356W	65.56238	0.00004
F380M	65.57509	0.00003
F430M	65.57361	0.00003
F444W	65.56745	0.00003
F480M	65.55983	0.00003

Table 2. Table A1: Information about the processing of the \jwst data used in this work.

Program ID	CRDS_CTX	CAL_VER	SDP_VER	PRD_VER
1086	jwst_0937.pmap	1.5.3	2022_2a	PRDOPSSOC-055
1334	jwst_0874.pmap	1.5.2	2022_2	PRDOPSSOC-052
1501	jwst_0977.pmap	1.6.2	2022_3	PRDOPSSOC-057

Equations

$$P_{i,j}=z_{\ast }\cdot {\psi }_{\mathrm{E}}(i-x_{\ast },j-y_{\ast })+s_{\ast },$$

$${\psi }_{\mathrm{E}}(i-x_{\ast },j-y_{\ast })=\frac{P_{i,j}-s_{\ast }}{z_{\ast }}.$$

$$\left\{\begin{array}[]{c}\tilde{x}=\frac{x-x_{\textrm{ref}}}{x_{\textrm{ref}}}\\ \tilde{y}=\frac{y-y_{\textrm{ref}}}{y_{\textrm{ref}}}\end{array}\right.$$

$$\{\begin{array}{r}\delta x=\sum _{i=1,5}\sum _{j=1,5-i}{a}_{ij}{\tilde{x}}^i{\tilde{y}}^j\\ \delta y=\sum _{i=1,5}\sum _{j=1,5-i}{b}_{ij}{\tilde{x}}^i{\tilde{y}}^j\end{array}$$

$$({\mu }_{\alpha }\cos \delta ,{\mu }_{\delta })=(1.873\pm 0.003,0.321\pm 0.002){\text{ mas yr}}^{-1}.$$

Full text

astrometry – photometry – proper motions – stellar clusters – Large Magellanic Cloud

††journal: ApJ††facilities: \jwst(NIRISS), HST, Gaia††software: FORTRAN, python, astropy, photutils

\suppressAffiliations

1 Introduction

\jwst

’s primary objectives (Gardner et al., 2023) are to study the first galaxies in the Universe, shed light on the birth of stars and planets, and study/explore other worlds, but its cutting-edge design and state-of-the-art detectors are poised to revolutionize astronomy over a much broader context.

From simple source centroiding to large-area image mosaics, a multifaceted variety of scientific applications with \jwst, from stars and exoplanets in the Solar neighborhood to structures in the Local Group, requires accurate and/or precise astrometry and photometry that may (and often do) rely on better calibrations than those generally offered by the official calibration pipeline (e.g., distortion uncertainty not to exceed 5 mas per coordinate in all detectors; see Anderson 2016111Document JWST-STScI-005361, “Verification of Plan to Solve for the Distortion Solution”.).

Here we begin to extend the techniques that made the Hubble Space Telescope (HST) and its cameras among the best resources for high-precision photometry and astrometry (e.g., Anderson & King, 2000, 2004, 2006) to \jwst’s Near Infrared Imager and Slitless Spectrograph (NIRISS) camera (Doyon et al., 2012). NIRISS provides various observing modes, including imaging, and—because of its detector size and filter properties—it can be considered a complementary channel to the Near InfaRed Camera (NIRCam). In Cycle 1, NIRISS imaging was offered only for parallel observations but, starting from Cycle 2, it can also be used as the primary instrument. So far, most of the attention from the imaging community has been on NIRCam (Nardiello et al., 2022; Griggio et al., 2022), while many of the astrometric and photometric capabilities of NIRISS remain to be demonstrated.

We present careful procedures with which we modeled the effective point-spread function (ePSF222We refer to the PSF discussed and solved for here as the “effective” PSF, in that it represents the “instrumental” PSF that impinges on the detector convolved with an image pixel., Anderson & King, 2000) and geometric-distortion (GD) correction of the NIRISS detector for all 12 of its imaging filters. We also briefly discuss the quality and stability of our models over the limited temporal baseline available ($\sim$2 months). Finally, we demonstrate the scientific potential of our methods and calibrations by computing proper motions (PMs) of stars in a field in the Large Magellanic Cloud (LMC).

This paper is the first of many to be written by the JWST Telescope Scientist Team333https://www.stsci.edu/$\sim$marel/jwsttelsciteam.html.. The team, led by M. Mountain, was convened in 2002, following a competitive NASA selection process. In addition to providing scientific support for observatory development through launch and commissioning (e.g., Gersh-Range et al., 2012; Perrin et al., 2014; Gersh-Range & Perrin, 2015; Perrin et al., 2016; Laginja et al., 2018), the team was awarded 210 hours of Guaranteed Time Observer (GTO) time. This time is being used over the first three \jwst observing cycles for studies in three different subject areas: (a) Transiting Exoplanet Spectroscopy (lead: N. Lewis); (b) High Contrast Imaging of Exoplanetary Systems (lead: M. Perrin); and (c) Local Group Proper Motion Science (lead: R. van der Marel). A common theme of these investigations is the desire to pursue and demonstrate science for the astronomical community at the limits of what is made possible by the exquisite optics and stability of \jwst. The present paper is part of our work on Local Group Proper Motion Science, which focuses on the Galactic Center and selected satellite companions of both the Milky Way and Andromeda galaxies. We present here methods, calibrations, and scientific validations using archival data outside our GTO program that are essential for our subsequent analysis of new GTO data obtained with the NIRISS instrument.

The paper is organized as follows. Section 2 describes the data sets used in our analyses. Sections 3 and 4 describe the methodology used to obtain accurate ePSF models and GD corrections, respectively, and include analyses of their quality and temporal stability. Section 5 presents a scientific application focused on measuring the PM kinematics of various populations in the \jwstcalibration field in the Large Magellanic Cloud (LMC). Finally, Section 6 presents the ePSF models and GD corrections that we release to the community with this paper. Section 7 summarizes our main conclusions.

2 Data sets

We mostly made use of data from Commissioning program 1086 (PI: Martel) that observed the \jwstcalibration field in the LMC centered at $(\alpha,\delta)=(80.5,-69.5)$ deg (Anderson et al., 2021) on 2022 May 1. The data set consists of 9 pointings over a 3$\times$3 dither pattern and two images per pointing for each of the 12 NIRISS filters. Each exposure was taken with the NISRAPID444https://jwst-docs.stsci.edu/jwst-near-infrared-imager-and-slitless-spectrograph/niriss-instrumentation/niriss-detector-overview/niriss-detector-readout-patterns. readout pattern, one integration, and between five and 15 groups, depending on the filter.

Our external validation and stability analyses made use of two additional data sets: the Early Release Science (ERS) program 1334 (PI: Weisz) and the Cycle-1 calibration program 1501 (PI: Sohn). The ERS program observed the globular cluster M92 through the F090W and F150W filters using the NIS readout pattern, one integration, and seven groups. These ERS NIRISS images were taken between 2022 June 20-21, as parallel observations to NIRCam primary exposures, and as such the dither pattern was driven by NIRCam’s needs. The Cycle-1 images were taken on 2022 September 18 and were focused on the LMC astrometric field for calibration purposes, and consist of four exposures per filter with the same readout pattern of Commissioning program 1086. The present study requires the use of unresampled images, which we downloaded from the MAST555https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html. archive as level-2, _cal products. The three data sets666All of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via https://doi.org/10.17909/hv07-tg13 (catalog DOI: 10.17909/hv07-tg13). were downloaded from the archive at different times, and they have been calibrated using different reference calibration files given the fast update rate of calibration reference files in the Calibration Reference Data System (CRDS)777https://jwst-crds.stsci.edu/. during Commissioning and in the early days of Cycle-1 (see Table A1 in Appendix A).

3 ePSF modeling

The first step to achieve precise measurements of positions and fluxes for point-like sources is to obtain accurate ePSF models. The NIRISS ePSFs are undersampled below 3.8 $\mu$m, i.e., for 8 out of its 12 filters, and extra care must be taken in modeling the ePSFs to avoid systematic effects.

The construction of accurate PSFs involves an inherently circular problem: to extract a PSF from an image, we need accurate positions of stars, but we cannot measure accurate positions of stars without an accurate model of the PSF. The use of non-optimal ePSFs tends to result in positions that are affected by significant pixel-phase biases, i.e., positions that are systematically measured at a specific location with respect to the pixel boundaries regardless of where sources truly are. Anderson & King (2000) showed that well-dithered images can help to break this degeneracy, since dithering allows a given star to land at different locations with respect to the pixel boundaries. The use of a common reference frame, in which positions of sources measured in different dithered exposures are averaged together, also helps in mitigating pixel-phase systematics, as we can transform these averaged positions back into the raw coordinate system of each individual _cal image and use them as improved positions from which to construct the ePSF models. The process is iterated to break the circularity problem mentioned above: once improved ePSF models are made, we can measure improved stellar positions in each exposure, resulting in a more precise averaged position on the common frame that can be inverted to construct even better ePSF models. The cycle is repeated until convergence is reached. We closely follow the prescriptions provided in Anderson & King (2000), Anderson (2016), and Libralato et al. (2016, 2021). We refer the interested readers to these papers for a more detailed description of the method, and we provide a brief overview below.

The value $P_{i,j}$ of a pixel $(i,j)$ in the vicinity of a star located at $(x_{\ast},y_{\ast})$ can be defined as:

[TABLE]

where $z_{\ast}$ is the stellar flux, $s_{\ast}$ is the local sky background, and the fraction of star’s light that should fall on that pixel according to the ePSF model is given by $\psi_{\rm E}(i-x_{\ast},j-y_{\ast})$. The ePSF fit allows us to obtain the stellar parameters $(x_{\ast},y_{\ast},z_{\ast})$. Since we are interested in modeling the ePSF, we invert the equation above to obtain:

[TABLE]

The knowledge of $(x_{\ast},y_{\ast},z_{\ast})$ for a star makes each pixel in the star’s vicinity a potential sampling of the ePSF. Using many samplings from many stars in many images that uniformly map the pixel space, we can build reliable ePSF models.

The iterative process described above starts by measuring positions and fluxes of bright, isolated stars in each individual NIRISS image. Positions were initially obtained by simply estimating the light-center of the flux distribution in the innermost pixels of each source, and fluxes were initially measured via aperture photometry. As soon as ePSF models are available, these stellar parameters are re-measured via ePSF fitting.

The same stars were then cross-identified in multiple catalogs, and their positions and fluxes were averaged together after being transformed on to a common reference-frame system (hereafter, the master frame888Orientation and scale for the master frame was set up using the HST catalog of Anderson et al. (2021), while the photometry was registered to that of a NIRISS image at the center of the dither patter; see also Sect. 5.) by means of six-parameter linear transformations. These averaged positions and fluxes, once transformed back into the frames of each exposure, represented the (current) best estimates of $(x_{\ast},y_{\ast},z_{\ast})$, which were used to collect the samplings for the ePSF reconstruction. It is at this point that a GD correction is often necessary to transform positions back and forth between different frames999If the dither is small and the distortion is small, then sometimes the differential distortion among the exposures can be ignored, but the dithers here were large.. Without GD correction, the ePSF samplings would be systematically misplaced, thus preventing us from mitigating pixel-phase biases. For this reason, we made use of a GD correction derived following the procedure highlighted in Sect. 4. This correction was itself improved during the iteration process.

Thanks to the availability of many thousands of stars in several dithered images, we were able to model the ePSFs using an oversampling factor of 4. We initially constrained one single ePSF model for the entire detector, but later on during the iteration process we let the ePSF model vary spatially. We found that a 5$\times$5 array of ePSF models is able to fully capture the spatial variability of the ePSF across the field of view (FoV) of NIRISS. At the beginning of each iteration, newly-derived ePSF models were used to measure improved positions and fluxes, which in turn defined an improved master frame. Convergence of the ePSF models was reached in 20 iterations, with subsequent iterations providing no significant improvements. During the last few iterations, we allowed the ePSF fitting procedure to also fit source positions instead of imposing positions from the master frame.

As discussed in Anderson & King (2000) and Anderson (2016), some constraints should be applied to the ePSF models to make sure they are smooth and continuous. We achieved this by convolving the ePSFs with a series of low-pass filters that allowed less high-frequency structure at larger radii. These smoothing kernels were empirically chosen by trial and error (see also Anderson, 2016) and are a compromise between the need to have smooth and continuous ePSF models while preserving true sharp ePSF features. The ePSF models were also normalized to sum up to 1 within a radius of 5 NIRISS pixels (chosen arbitrarily).

Figures 1, 3 and 3 illustrate the main pixel-phase issues that affect undersampled detectors and that can be minimized using a careful modeling of the ePSFs. The Figures refer to the F090W case, which is the most undersampled among the NIRISS filters. At the beginning of our iterative process: (i) stars are preferentially measured to be at the corners of the pixels, even though their true position is random within the pixels (left panel of Fig. 1); (ii) the core of the ePSF is too sharp and not well constrained (left panels of Fig. 3); and (iii) the positional residuals (defined as the difference between the raw positions in each exposure and those predicted by the inverted averaged master frame) show a clear trend as a function of pixel phase (left panels of Fig. 3). All these features are caused by the undersampled nature of the ePSF: a sharp PSF does not transfer flux linearly from pixel to pixel, and therefore a linear center of light model mis-measures positions. After 20 iterations in which the true shape of the PSF is teased out: (i) stellar positions measured by fitting our precise ePSF models do not show any preferential location with respect to the pixel boundaries (right panel of Fig. 1); (ii) the ePSF samplings are smooth, with less scatter than before, and the shape of the ePSF core is well constrained (right panels of Fig. 3); and (iii) pixel-phase errors have been significantly reduced (right panels of Fig. 3).

To summarize, our final ePSF models are spatially-variable 5$\times$5 arrays of ePSFs covering the full NIRISS detector. Each ePSF has a size of 101$\times$101 oversampled pixels, with the ePSF centered at location (51,51), extending out 12.5 real NIRISS pixels in each direction, and it is normalized to have unity flux within a radius of 5 NIRISS pixels. The core of the ePSF models becomes progressively less sharp towards redder wavelengths (i.e., the less undersampled the ePSFs are), with the fraction of the flux in the centermost pixel of the ePSF model decreasing from $\sim$46% (F090W filter) to $\sim$11% (F480M filter).

The left panel of Figure 4 shows the full 5$\times$5 ePSF array for the F090W filter. The differences between each spatially-varying ePSF and the average ePSF are shown in the right panel. The spatial variation of the F090W models causes about $6\%$ of the central-pixel flux to be redistributed into the wings. The ePSFs of redder filters are found to be less affected by spatial variations, with as low as only $0.3\%$ of the central-pixel flux being redistributed in the reddest F480M filter.

In the following, we refer to these final ePSF models as the library ePSFs. We will discuss strategies to further improve these models on an image-by-image basis in Sect. 3.2.

3.1 Quality of our ePSFs

We used two diagnostics to test our library ePSF models. The first one is the “Quality of PSF fit” or QFIT (Libralato et al., 2014), which is defined as the absolute fractional error in the fit of a source. Well-measured stars have a QFIT close to 0, whereas increasingly poorly-measured objects have increasingly larger QFIT values. The result of this first test for one indicative image per filter from the Commissioning data set is shown in Fig. 5. The red lines are set at the median QFIT values of the bright stars in a region where the QFIT trend is flat (shaded pink regions). These medians, reported on the left side of each panel, are always below 0.05, a value generally used as a proxy for adequate ePSF models. The QFIT increases towards fainter magnitudes due to the lower signal-to-noise ratio. For the bluest, more undersampled filters, the brightest stars also show an increase of the QFIT trend. This is likely due to a combination of saturation, non-linearity and/or brighter-fatter effects. We also do not exclude that quantum yield101010Quantum yield describes the possibility, for the most energetic photons, to generate more than one electron-hole pair per detection. could play a role for the bluest filters. We investigated the possible presence of chromatic effects in our ePSFs by comparing the QFIT values of the bluest and reddest stars with high signal to noise, and found no systematic differences. This is somewhat expected since we did not restrict our ePSF modeling to stars of specific colors.

The second diagnostic parameter is RADXS: the excess/deficiency of flux outside the core of the sources with respect to the ePSF prediction (Bedin et al., 2008). RADXS is a powerful tool to discriminate between point-like sources (RADXS $\sim 0$), cosmic rays or hot pixels (negative RADXS), and galaxies (positive RADXS). If our ePSF models are a good representation of the flux distribution of point-like sources outside their cores, the RADXS of stars should be close to zero. Figure 6 presents RADXS trends as a function of instrumental magnitude111111The instrumental magnitude is defined as $-2.5\log$(total surface brightness), where “total surface brightness” is the total surface brightness (level-2 _cal FITS files are in units of MJy sr*-1*) under the fitted ePSF. for the same exposures used in Fig. 5. While the general RADXS trends are all around zero, small gradients are visible in all but the reddest filters, with RADXS preferentially positive for bright sources and preferentially negative for faint sources. We suspect this behavior is due to brighter-fatter effects (e.g., Antilogus et al., 2014; Plazas et al., 2018), which are known to be present in the NIRISS detector121212https://jwst-docs.stsci.edu/jwst-near-infrared-imager-and-slitless-spectrograph/niriss-instrumentation/niriss-detector-overview/niriss-detector-performance#NIRISSDetectorPerformance-Detectorfullwellcapacity.. The brighter-fatter effects causes photons that should land on pixels that are closer to saturation to end up generating a photo-electron in a neighboring pixel, resulting in bright point-like sources appearing “fatter”, i.e., broader, than fainter ones. Since our library ePSF models were constructed using stars within specific luminosity ranges, fitting sources outside these ranges might provide poorer results due to brighter-fatter effects. The systematic trends in Fig. 6 are more pronounced in bluer (more undersampled) filters. We find that the impact of the brighter-fatter effect on astrometry with the NIRISS data is very small, and is only relevant for the bluest, most undersampled filters. For example, Fig. 3 shows that after 20 iterations in our ePSF-modeling procedure, the pixel-phase errors are within $\pm$0.01 pixel along the $x$ direction. The amplitude of these residuals seems not to depend on the magnitude of the star, but the patterns for bright and faint stars have a different phase. This is another piece of evidence that the ePSF of bright and faint stars is slightly different because the flux distributions for the two groups of stars differ, as one would expect from the brighter-fatter effect. A detailed analysis of this phenomenon is beyond the scope of this paper, but we want to make users aware of its presence. Flux-dependent ePSF models, similarly to what is done in the Gaia DR4 data reduction (Rowell et al., 2021), might be able to produce a more accurate description of how the flux of bright and faint point-like sources is actually distributed on the detector.

3.2 Temporal stability

The true PSF of \jwstis expected to change over time, primarily due to small thermal variations, mirror movements including from occasional structural relaxations (“tilt events”), and mirror corrective moves from time to time. We monitored the temporal stability of our library ePSF models in two ways.

First, we studied the variation of flux in the centermost pixel of the library ePSFs (see Fig. 7 for variations at the center of the NIRISS detector) using all data sets at our disposal, which cover a temporal baseline of $\sim$51 days. The library ePSF models were perturbed on an image-by-image basis following the prescriptions given in, e.g., Bellini et al. (2017a) and Libralato et al. (2018a, 2022). The ePSF perturbation is an iterative process consisting of increasingly adjusting the models using residuals of the fit of well-measured sources. We report in the figure these perturbed central pixel values as a function of time for all 12 filters. The P2SF of \jwstappears to be relatively stable over the available temporal baseline, with deviations below $\sim$5% except for the F090W and F150W ERS data, for which we measured a central pixel variation of about $\sim$7% and 6%, respectively. This appears to be consistent with the stability of the Optical Telescope Element of \jwstdiscussed in McElwain et al. (2023). Temporal variations appear to be less prominent in redder filters.

The second test consisted of comparing the QFIT of well-measured, bright stars obtained by fitting the library ePSFs and that coming from perturbed models. We narrowed the investigation to just the ERS data where the largest discrepancies were found in the previous test. The top two panels of Fig. 8 show the results using the library ePSFs, while the bottom two panels are for perturbed ePSFs. It is clear that the perturbed models provide a better representation of point-like sources (lower QFIT values) in the ERS data, pointing to a real change in the telescope PSF during this time.

It is clear from these tests that perturbing the library ePSFs to better match the PSF in an individual exposure can result in significantly better fits to stars. Since the perturbation process does not include information from multiple dithers, it should not be expected to significantly improve astrometry, but it can improve photometry (particularly if the perturbation is not constant across the detector) and can significantly improve star-galaxy discrimination and multiple-star fitting procedures. One must be careful, however, to ensure a minimum number of high signal-to-noise stars in each perturbation zone.

4 Geometric-distortion correction

In this section, we describe the steps that lead us to achieve a very precise GD correction. In general, there are two main ways to solve for the GD (see an extensive discussion in, e.g., WFIRST Astrometry Working Group et al., 2019). The self-calibration approach makes use of observations characterized by large dither patterns and different roll angles. In this way, a star is imaged in different parts of the detector, ideally from corner to corner of the FoV. This helps to randomize the direction of systematic GD residuals, so that the average stellar position on a common reference frame provides a somewhat “distortion-free” representation of its true position that can be used to solve for the GD (e.g. Anderson & King, 2006; Libralato et al., 2014, 2015; Häberle et al., 2021). This approach works best when many observations at different dither positions are available, and hence requires a considerable amount of observing time, thus not a viable option for \jwstgiven that sufficiently diverse data are not available at this early phase of the mission.

The second approach, which is the one that was officially adopted by the \jwstinstrument teams during Commissioning, calibrates the GD through an external catalog.131313A combination of these two approaches will likely be used for the Nancy Grace Roman Space Telescope (WFIRST Astrometry Working Group et al., 2019). To this end, Anderson et al. (2021) used two HST data sets of the same LMC field taken 11 years apart to construct a high-precision, astro-photometric catalog for the GD correction of \jwst’s imagers. This reference catalog offers the ability to calibrate \jwst’s astrometry to a precision level similar to that of HST.

We also decided to follow the second approach based on the Anderson et al. (2021) catalog as a reference, with the goal of achieving an astrometric precision better than 0.01 pixel. Unaccounted-for internal motions of LMC stars ($\lesssim$50 km s*-1*, which at a distance of 49.5 kpc, correspond to 2.1 mas or $\sim$0.03 NIRISS pixel over 10 years) would prevent us from reliably detecting and correcting fine details in the GD of NIRISS. To mitigate this issue, we epoch-matched stellar positions of Anderson et al.’s catalog to the average epoch of the NIRISS Commissioning observations using the PM information provided by the catalog. We additionally selected only stars that are bright and relatively isolated in NIRISS images.

We started by measuring stellar positions in each NIRISS exposure by fitting the corresponding perturbed ePSF models, and kept only those with QFIT lower than 0.1 for the next stages. In the following, we describe each step of our GD correction using the F277W-filter data as guidance, but the same methodology and conclusions extend to all filters.

4.1 Polynomial correction

The bulk of the non-linear terms of the GD were fit with two fifth-order polynomial functions relative to the center of the NIRISS detector: $(x_{\textrm{ref}},y_{\textrm{ref}})=(1024,1024)$ pixels. In addition, the use of normalized positions:

[TABLE]

guarantees that the values of the coefficients of each polynomial term tell us the size of the distortion, center-to-corner, of that term in pixel units. The distortion correction for each star is defined as:

[TABLE]

The values of the $a_{ij}$ and $b_{ij}$ coefficients were determined through an iterative procedure as described in Bellini & Bedin (2010). First, we cross-identified stars between each NIRISS catalog and the HST reference catalog, and transformed the epoch-matched HST positions of common stars on to the raw coordinate frame of each NIRISS catalog by means of four-parameter linear transformations (two offsets, one rotation and one change of scale). The use of four parameters here instead of the usual six allows us to include the off-axis linear “skew” terms of the GD in the polynomial solution.

Once in the same raw coordinate frame, each cross-identification generates a pair of positional residuals $(\Delta x,\Delta y)$ defined as the difference between the transformed HST positions and the raw NIRISS positions. These residuals were iteratively 3$\sigma$-clipped and averaged into a look-up table with 20$\times$20 square elements covering the full NIRISS detector, to minimize single-measurement noise. The averaged residuals were least-squares fit to obtain the coefficients of the two fifth-order polynomial functions. To ensure a smooth convergence of the solution, only half the values of the computed coefficients were added to the tally, and the process was repeated until the values in the look-up tables were smaller than $10^{-5}$ pixel.

Figure 9 shows the F277W, two-dimensional distortion maps and the ($\Delta x$,$\Delta y$) positional residuals as a function of the $x$ and $y$ raw NIRISS coordinates before (left panels) and after (right panels) our polynomial correction. The NIRISS distortion has a magnitude of about 10 pixels from corner to corner of the FoV, with an RMS of 2.6 pixels. After the polynomial correction, the impact of systematic residuals is significantly lower, reaching a corner-to-corner amplitude of $\sim$0.1 NIRISS pixel and an RMS of 0.03 pixel. To show detail, the size of the vectors in the two-dimensional maps is magnified by factors of 50 and 5 000 for the before and after correction, respectively. Note that the orientation of the two-dimensional residuals after the correction is applied is significantly not random; this is even more evident by looking at the wavy residual patterns in the $\Delta x$ versus $x$ and $\Delta y$ versus $y$ side panels on the right.

4.2 Look-up table of residual corrections

It is clear from Fig. 9 that there are higher-frequency variations in the distortion that cannot be accounted for with our fifth-order polynomial corrections. A better view of the wavy pattern present in the side panels of Fig. 9 is provided in Fig. 10, where $\Delta x$ residuals are shown as a function of $x$ position in four different 512-pixel-tall bands. The NIRISS readout electronics consist of four amplifiers that read 512 rows of pixels each, suggesting readout electronics could, at least in part, be responsible for the presence of these wavy patterns, in a similar fashion to the hysteresis effects in the readout electronics of the HAWK-I@VLT camera (see, e.g., Libralato et al., 2014).

Regardless of the true nature of these patterns –the investigation of which is outside the scope of the present work– their correction can be empirically modeled using a look-up table of residuals, as done in, e.g., Bellini et al. (2011). Polynomial-correction residuals were collected and averaged (with iterative 3$\sigma$ clipping) over a 32$\times$32 look-up table (square elements of 64 pixels per side), which allowed us to remove the vast majority of what was left after the polynomial correction. The look-up table is populated in such a way that elements adjacent to the edges of the detectors are applied as if they were measured exactly on the edges, to allow the use of bi-linear interpolation of the correction across the entire table (see the scheme in Libralato et al., 2014). As for the polynomial-correction derivation, we iterated the look-up table part of the distortion solution and only applied 50% of the correction values each time until convergence, which was reached after 250 iterations.

The two-dimensional distortion map and the one-dimensional positional residuals after both the polynomial and the look-up table residual corrections were applied are shown in Fig. 11. The figure shows no significant residual trends in either the $x$ or $y$ direction. The RMS of the positional residuals is $\lesssim$0.01 NIRISS pixel.

Another quantitative assessment of the goodness of our distortion solution was obtained as follows. We constructed a common master frame of averaged positions by combining all 18 F277W single-exposure catalogs corrected with our polynomial and look-up table corrections. We then inverted the averaged positions of about 10 000 well-measured stars with instrumental magnitude between $-6.5$ and $-5$ back into the raw coordinate systems of each exposure, and collected the differences $(\Delta x,\Delta y)$ between them and the corresponding single-exposure measurements. These residuals are shown in Fig. 12. We also show the one-dimensional histogram of the residuals along each axis and fit them with Gaussian functions (red lines in the side panels). It is clear that: (i) the two-dimensional distribution of the positional residuals is “round”, which is what we would expect when no systematic errors are present; (ii) the $\sigma$ of the fitted Gaussian functions are of the order of 0.01 NIRISS pixel or less, which is comparable to the state-of-the-art distortion solutions of HST’s imagers; and (iii) the one-dimensional distributions are almost-perfectly fit by Gaussian functions, without any statistically-significant deviation in the wings. These are all indicators of the high quality of our GD correction.

4.3 Epoxy voids

The results shown in Figs. 11 and 12 clearly demonstrate the adequacy of our GD solutions, and we could have stopped here without further tests and checks. Nonetheless, we decided to look at the F277W residuals one more time, scanning through the $x$ and $y$ raw-coordinate space in 10-pixel-wide bands. It is during this final test that we noticed a suspicious, possibly significant and localized $\Delta x$ residual trend for $y$ coordinates between 800 and 900 pixels. A negative $\Delta x$ peak, of about 0.05 pixel, appeared around $x=1550$ pixels, and progressively shifted to $x=1250$ while becoming positive by approximately the same amount as we progressed the $x$-coordinate scan from $y=800$ to $y=900$. We also noticed a few other possible smaller systematic trends in different locations on the detector, but their significance was at best marginal.

We then looked at the residuals in the other filters, and found very similar trends at around the same location in raw-coordinate space and with about the same amplitudes in all filters. These findings suggest that what we saw while scanning the F277W residuals could be a real distortion feature, rather than something related to ePSF models and/or a particular filter, and strongly favors a detector origin.

If indeed the origin is related to the detector, by collating the distortion residuals from all images of all filters we should be able to maximize its signal. The two panels of Fig. 13 show the results coming from this collation of $\sim$1 000 000 residuals. A distinct structure, mainly oval in shape and as large as 1/6 of the detector, stands out in the $\Delta x$ residual map on the left, and traces of it are also present in the $\Delta y$ residual map on the right. An ellipse on the left panel marks the location of the trend initially identified during the F277W scanning. Also worth noting is what seems to be a vertical-band pattern in the left panel, and a similarly-spaced, horizontal-band pattern in the right panel.

It is usually the case that pixels affected by detector-related residual patterns also stand out in flat-field images (see, e.g., Bellini et al., 2011). Figure 14 shows the F277W flat-field calibration reference file obtained from the \jwstCRDS. The large pattern identified in Fig. 13, as well as other smaller substructures related to secondary residual patterns are clearly visible.

Most mercury cadmium telluride (HgCdTe) diodes are grown on a substrate that is then removed to increase sensitivity at short wavelengths below 800 nm and to lower dark-current levels (e.g., Murakami et al., 2003). To improve mechanical stability, the gap between the HgCdTe diode and the multiplexer is then filled with epoxy (see also Brown et al., 2006). The features we identified in the distortion residuals and in the flat-field reference file are found in correspondence to voids in the epoxy underfill. While epoxy voids have been studied with respect to dark current and photometric performance141414For the specific case of NIRISS, see: https://jwst-docs.stsci.edu/jwst-near-infrared-imager-and-slitless-spectrograph/niriss-instrumentation/niriss-detector-overview/niriss-detector-performance. (e.g., Regan & Bergeron, 2020), to our knowledge no one was expecting them to also have a measurable impact on astrometry, as we have discovered here for the NIRISS detector.

To account for these fine structures in our distortion solution, we opted for a pixel-based correction rather than an additional look-up table of residuals. The correction is made up of two parts, one for the $x$ direction and one for the $y$ direction, and is derived as follows: we assigned to each pixel of the 2048$\times$2048 pixel array the 2.5$\sigma$-clipped median of the closest 50 residuals coming from all filters. This was done as a compromise between smoothing out most of the noise coming from the individual measurements and preserving the residual patterns at both small and large scales. After including this new correction, the systematic trends in the positional residuals have been significantly reduced (Fig. 15).

4.4 Astrometric precision

The distortion maps and residuals presented in the previous Sections are clear qualitative indicators of the accuracy of our GD corrections. Quantitatively, state-of-the-art GD corrections for HST have residuals $\lesssim$0.01 pixel (e.g., Anderson & King, 2006; Bellini et al., 2011), measured as the RMS around the mean of averaged dithered positions for the same sources. Note that these residuals are the sum in quadrature of ePSF-fitting errors and the GD-correction errors themselves. Therefore, they are more an indicator of how precise the whole astrometry package is than just the GD correction alone.

Our quantitative assessment of the astrometric precision reached with the NIRISS Commissioning data was done as follows. For each filter, we averaged together positions of sources measured in all NIRISS single-exposure catalogs, using the full, 3-part GD solutions and six-parameter linear transformations. The process was iterated three times to improve both the averaged positions on the master frame and the transformations. Only stars measured in at least three images were considered in the analysis. The one-dimensional master-frame positional RMS as a function of instrumental magnitude in all filters are presented in Fig. 16. The red lines in each panel are set at the median values of the RMS of the best-measured stars that lie within the shaded pink regions. Pixel values were transformed to mas using the average plate scale of all pixels, which is listed at the top of the Figure (see also Sect. 4.6). For all filters, the median values of the positional RMS are well below 0.01 NIRISS pixel ($<$0.66 mas). Pixel-wise, the precision reached is superior to that achieved with HST’s imagers, which is expected given the more stable environment of \jwst.

4.5 Filter dependence

Any optical element traversed by light upstream of the detector can contribute to the GD. So, we investigated whether, and to what extent, the filter components affect the shape and amplitude of the GD.

We averaged together the total two-dimensional GD correction of all filters to obtain a proxy for a filter-independent solution. Then, we subtracted it from the GD correction of each filter to isolate their filter-dependent component. The result is shown in Fig. 17. The different patterns in each panel point to significant (up to 0.5-pixel) filter dependencies in the GD correction. Residuals $\lesssim$0.5 pixel are small enough for some scientific applications (e.g., catalog cross-matching), but are definitely too large for high-precision astrometric studies. This validates our choice of deriving filter-dependent GD corrections.

4.6 Pixel scale

The HST catalog of Anderson et al. (2021) provides us with an opportunity to measure the absolute pixel scale of the NIRISS detector. For each filter, we cross-identified bright, well-measured stars between each single-exposure NIRISS catalog of the Commissioning data and the HST catalog of Anderson et al. (2021) and used six-parameter linear transformations to relate the two sets of positions. One of these linear parameters is the relative scale between NIRISS and the HST catalog, which was defined to be exactly 50 mas pixel*-1*. We use this information to estimate the absolute pixel scale of NIRISS, allowing for a dependence on filter.

Velocity aberration causes cyclical changes in a detector’s pixel scale as \jwstorbits the Earth-Sun system (e.g., for HST see Cox & Gilliland, 2003). Although small, these scale changes need to be accounted for in the computation of the absolute scale of a detector. The velocity-aberration correction factor is provided in the header of every NIRISS image with the keyword VA_SCALE. We used it to correct the values of the NIRISS pixel scales that we derived from the HST catalog cross-matching. It is worth mentioning that this correction is quite small for the specific observations used, because–by design–the LMC calibration field targeted for Commissioning is in the Continuous Viewing Zone of \jwst, and the pointing of the telescope there is almost always orthogonal to its orbit around the Sun.

Table 1 collects our measurements of the absolute pixel scales for each filter. We find significant scale differences among the filters, which is consistent with what is stated in the official documentation of NIRISS151515https://jwst-docs.stsci.edu/jwst-near-infrared-imager-and-slitless-spectrograph/niriss-instrumentation/niriss-detector-overview.. We measure the filter-averaged scale of NIRISS to be 65.562 $\pm$ 0.003 mas pixel*-1*.

We also investigated the temporal variation of the pixel scale using all available data sets (Fig. 18). For the ERS data, which are focused on the globular cluster M 92, we used the Gaia Data Release 3 (DR3) catalog (Gaia Collaboration et al., 2016, 2021a, projected on to a tangent plane centered at the center of one ERS NIRISS image, with a fixed scale of 66 mas pixel*-1*) as a reference to compute the pixel scale of each image. Unfortunately, the only stars in common between the ERS data and the Gaia catalog are at the faint end of Gaia (${\it G}>19$), which resulted in less precise transformations in Fig. 18.

We find the scale of Cycle-1 images to be significantly different (at the 8-$\sigma$ level) from that of Commissioning data for the F090W filter. The scale difference becomes progressively less noticeable the redder the filters, down to less than 0.7$\sigma$ for the F480M filter.

Figure 19 shows the relative pixel scale variation in Commissioning and Cycle 1 images. There is a clear apparent trend as a function of time during the $\sim$12-hour duration of the Commissioning observations. While many potentially contributing factors were changing during Commissioning, there were no commanded adjustments in telescope mirror alignment or instrument focus (McElwain et al., 2023) over the time period of these measurements and, further, the observatory had already reached thermal stabilization weeks earlier, such that there were no substantial ongoing thermal drifts (Menzel et al. 2023, PASP in press). We also computed the scale against the HST-based reference frame (Anderson et al., 2021), against a pure Gaia DR3 frame (for the reddest filter that has more stars with $G<19$ in common with the Gaia catalog), and against the central-pointing NIRISS image, and in all cases saw a similar trend. If we hypothesize that the change in pixel scale is due to a drift in the Optical Telescope Element (OTE) secondary mirror position, the necessary amount of secondary-mirror defocus to change the pixel scale by $5\times 10^{-6}$ according to optical modeling is significantly inconsistent with the measured stability in the wavefront sensing observations over this time period. The origin of these apparent changes in pixel scale remains unclear. Continued monitoring of the NIRISS pixel scale over time may help understand the cause of this subtle effect.

5 Sample Scientific Application: Population-dependent kinematics of the LMC

An immediate scientific application demonstrating the quality of our ePSF models and GD solutions is provided by the Commissioning data. We computed PMs of sources in the LMC calibration field by combining HST and \jwstdata. As a first-epoch data set, we took advantage of the catalog made by Anderson et al. (2021) using HST data taken at epoch 2006.39. This catalog defines our first-epoch master frame. It is oriented with North up and East to the left, and it has scale of 50 mas pixel*-1*. We created filter-dependent NIRISS master frames by averaging together positions and fluxes of common sources measured in the individual images. Orientation and scale for the NIRISS master frames were set up using the HST catalog, while the photometry was registered to that of a NIRISS image at the center of the dither pattern. For PM purposes, we used the NIRISS F090W as our second-epoch master frame, since it provides the highest number of sources in common with the HST catalog. The master frames of the remaining NIRISS filters are used for additional photometric information in our analysis.

We transformed positions as measured in the NIRISS master frame into to the reference frame of the 2006.39 HST catalog by means of general, six parameters linear transformations. The coefficients of these transformations were obtained using a set of bright, well-measured, old red-giant branch (RGB) stars –selected in the instrumental F090W versus (F090W$-$F150W) color-magnitude diagram (CMD)– as reference stars, as shown in the left panel of Fig. 21. RGB stars define a narrow sequence in the CMD, as opposed to the very broad blue sequence typical of rich centrally-located LMC fields, which is a superposition of populations with different ages, chemical compositions and, most importantly, kinematics (Gaia Collaboration et al., 2021b).

Once in the same reference system, relative PMs were computed as the difference between the original HST and transformed \jwstpositions, divided by the temporal baseline ($\sim$15.94 yr). We multiplied the result by the first-epoch pixel scale of 50 mas pixel*-1* to obtain PMs in units of mas yr*-1*. PM errors were estimated as the sum in quadrature of the positional errors of the two catalogs divided by the temporal baseline. Finally, our relative PMs were registered onto an absolute reference frame by means of the absolute PM zero-point provided by the comparison between our relative PMs and the absolute PMs in the Anderson et al. (2021) catalog.

The right panel in Fig. 21 presents the one-dimensional PM error as a function of F090W instrumental magnitude. For bright, well-measured stars in the instrumental magnitude range between $-6.5$ and $-5$, the median PM error is $\sim$13 $\mu$as yr*-1*. This corresponds to 3.1 km s*-1* at the distance of the LMC (1 mas yr*-1* $\sim$ 235 km s*-1*; van der Marel et al., 2002). This is exquisite for studies of the internal kinematics of the LMC disk. Stars at these instrumental magnitudes have apparent magnitudes around ${\it G}\sim 20$ (given the ${\it G}-{\rm F090W}$ colors in the right panel of Fig. 21). By contrast, even for its best-measured stars ${\it G}\sim 9$–12, Gaia DR3 only achieves uncertainties of $\sim 16$–23 $\mu$as yr*-1* (Table 4 of Lindegren et al., 2021). Around the LMC tip of the RGB with ${\it G}\sim 16$ (Gaia Collaboration et al., 2021b), Gaia DR3 median errors are $\sim 38$–65 $\mu$as yr*-1*, some 3–5 times larger than what we report here for stars that are some 4 magnitudes fainter. The uncertainties reported here are comparable to those obtained in state-of-the-art PM studies combining HST-only data over a similar temporal baseline (e.g., Bellini et al., 2014; Libralato et al., 2018b, 2022). This is a testament of the overall quality of our \jwstNIRISS ePSF and GD models.

We compared our PMs with those in the Gaia DR3 catalog. Unsaturated stars in common with the Gaia catalog are between ${\it G}=17$ and ${\it G}=20$. Figure 21 summarizes the results of the comparison. The vector-point diagrams (VPDs) of stars in three magnitude bins – each width ${\it G}=1$ mag – are shown for our HST$+$\jwstPMs (left panels) and Gaia’s PMs (middle panels). We report the observed PM dispersion in red in each panel. While Gaia’s PM errors are larger than the intrinsic dispersion of stars in the field, our PMs are clearly able to reveal substructures in the internal motion of LMC stars.

These substructures are population dependent, as can be seen from our samples illustrated in Fig. 22. We selected the bluest and the reddest stars in the F090W versus (F090W$-$F150W) CMD (left panel in Fig. 22) within the instrumental magnitude range $-6.5<\rm{F090W}<-5$. This is the range where stars have the best PM precision according to Fig. 21. Stars in the red selection are mostly old RGB stars (hereafter, the old population), while those in the blue selection are mostly a mix of young and intermediate-age sources with ages between 50 Myr and 2 Gyr (hereafter, the “young populations”, e.g., Gaia Collaboration et al., 2021b).

The distribution of stars in these two groups in the VPD is clearly different. The old population is characterized by a broad distribution centered around ($\mu_{\alpha}\cos\delta,\mu_{\delta})=(1.90,0.38)$ mas yr*-1*, whereas the young populations seem to be made up of three main groups (highlighted in the rightmost VPD in Fig. 22), with hints of some of them being further fragmented into additional substructures. Group 1 is the least populated and most isolated of the three; Group 2 seems to be the most concentrated, and it overlaps substantially with a broader Group 3.

Perhaps a clearer way to highlight substructures in the motion of the young populations is through the analysis of their histogram distribution. To this end, we used the methodology described in Bellini et al. (2017b) and summarized in Fig. 23. The distributions of the Group 2 and 3 stars in the VPD appear elongated along a specific direction (blue dashed line in the left panel). We projected the motions onto this line and then constructed a histogram (right panel). To remove any dependence on the bin width and the starting point of the histogram, we followed the approach of Libralato et al. (2019) and used 1 000 individual samplings of the data, each with a different amount of added PM noise according to the stars’ PM-error distribution. We simultaneously fit three Gaussian functions to the resulting histogram that are shown in gold, magenta and green in the right panel of the figure, for Groups 1, 2 and 3, respectively. The areas of these Gaussians are indicative of the relative abundance of young stars of the three groups in our FoV. We find that about 12% of the stars belong to Group 1, 20% to Group 2 and the remaining 68% to Group 3.

Among the young population substructures in the VPD, Group 2 stands out with the tightest distribution. It has a dispersion of $\sim$0.026 mas yr*-1* (i.e., 6.2 km/s), measured as the $\sigma$ of the Gaussian fitted in Fig. 23. and it is centered at $(\mu_{\alpha}\cos\delta,\mu_{\delta})\sim(1.8,0.3)$ mas yr*-1*. Stars within $\pm$3$\sigma$ from the center of the magenta Gaussian are preferentially located in proximity of the South-East edge of the FoV, at odds with members of the other two groups that have a more uniform spatial distribution.

Interestingly, Group-2 stars are found to be in the same region where the 120-Myr-old stellar cluster OGLE-CL LMC 407 is located (Narloch et al., 2022), and are likely members of this cluster. We then further restricted Group-2 objects to only those within a 25*′′* radius from the center of the overdensity of stars (defined by hand) seemingly associated with the OGLE-CL LMC 407 cluster, which is visible in the stacked HST image from the epoch-one data (see Fig. 24). Using this refined sample, we were able to measure the absolute PM of the cluster as the 3$\sigma$-clipped median of the distribution in the VPD:

[TABLE]

Figure 25 presents the CMD of cluster members from our refined sample and with a further constraint on a magnitude-dependent PM consistent with the cluster’s bulk motion. The selection radii in the VPD as a function of magnitude bin were chosen as a compromise between excluding as many of the field stars with cluster-like PMs as possible, while retaining a statistical sample of bona-fide members. Apart from a few intermediate-age and old contaminants, the main sequence of the cluster appears to be well defined in the brighter part of the CMD, while field-star contamination becomes important for $\textrm{F090W}>-3.5$, due to larger PM errors.

Finally, we computed the intrinsic velocity dispersion of the cluster using the maximum-likelihood approach described in Sect. 4 of Libralato et al. (2022). We considered only bright, well-measured stars with $-6.5<{\rm F090W}<-5$ within 25 arcsec from the center of the cluster, and with a PM within 0.091 mas yr*-1* (3.5 times the $\sigma$ of the magenta Gaussian in Fig. 23) from the bulk motion of the cluster. We find a velocity dispersion of $23.7\pm 1.8$ $\mu$as yr*-1*, i.e., $5.6\pm 0.4$ km s*-1* assuming the cluster is at the same distance of the LMC (49.5 kpc; Pietrzyński et al., 2019). LMC clusters with masses similar to that of OGLE-CL LMC 407 ($\sim$4000 $M_{\odot}$; Popescu et al., 2012) are expected to have velocity dispersions of the order of 1 km s*-1* (McLaughlin & van der Marel, 2005). It is thus possible that the true dispersion of OGLE-CL LMC 407 is in fact lower, in which case the value of the velocity dispersion we computed can be interpreted as an upper limit on the possible presence of uncorrected systematic errors in either the first-epoch HST or the second-epoch \jwstcatalogs, or in both.

The intrinsic velocity dispersion of RGB stars implied by Fig. 22 is $33.8\pm 0.6$ km s*-1* ($0.144\pm 0.003$ mas yr*-1*). The distribution for the young populations is significantly more compact than this (but with significant substructure, per the preceding discussion), implying lower velocity dispersions. This is consistent with our general theoretical understanding of disk-heating mechanisms, which imply that old stars tend to have higher velocity dispersions than young stars. The values we report here are comparable to what is implied by line-of-sight velocities of LMC tracers (Zhao et al., 2003; van der Marel et al., 2009; Olsen & Massey, 2007; van der Marel & Kallivayalil, 2014). Our findings are also consistent with those inferred from HST-only PMs computed by Anderson et al. (2021, see also their Fig. 15). As a result of asymmetric drift, older populations have lower disk rotation than young populations, as has been verified observationally (van der Marel & Kallivayalil, 2014; Gaia Collaboration et al., 2021b). The LMC is a ${\sim}10$ times less massive galaxy than the Milky Way (Bland-Hawthorn & Gerhard, 2016; Erkal et al., 2019). So it also makes physical sense that the LMC velocity dispersion is several times lower than that of the Milky Way bulge, which has a central velocity dispersion $\sim$130–140 km s*-1* (Sanders et al., 2019; Valenti et al., 2018).

It is thus intriguing that our results conflict with the large central velocity dispersion values reported for the LMC disk by Gaia Collaboration et al. (2021b). Our field is at 0.4 degrees from their adopted LMC center. At that distance, they report velocity dispersions for RGB stars and several young populations of order 70–80 km s*-1*, well in excess of what is implied by Fig. 22. We do not know the cause of this discrepancy. Our samples have significantly smaller PM uncertainties than the Gaia DR3 PMs, so we may just be better resolving the kinematic structure near the crowded LMC center. It is possible that blending effects may have artificially inflated the scatter in the Gaia measurements near the dense LMC center. Gaia Collaboration et al. (2021b) include a systemic error term of 0.01 mas yr*-1* in their analysis, but do not explicitly account for additional uncertainties introduced by blending and crowding. Rejection of stars based on the use of flags in the Gaia catalog will remove some of these, but maybe not all of them. Studies of runaway stars from the LMC 30 Dor region based on Gaia DR2 found that the PM scatter in the catalog in dense regions was significantly inflated (i.e., implying unrealistic numbers of runaway stars; Lennon et al., 2018). It will be of interest to reassess the velocity dispersions in the LMC disk based on future Gaia data releases that may include more sophisticated treatment of blended sources.

Overall, our astrometric studies clearly demonstrate the exquisite quality of our ePSF models and GD solutions, as well as of the overall astrometric and photometric potential of NIRISS, and that of \jwstin general. Moreover, the PMs for the science application presented here were obtained, for simplicity, by only using one twelfth of the available Commissioning data. A larger data set and more sophisticated techniques (e.g., Bellini et al., 2014; Libralato et al., 2021, 2022) instead of simple global transformations between two frames could potentially improve the obtainable PM precision from Commissioning data. With further analyses it would not be out of the question to obtain PM uncertainties for the LMC approaching 1–2 km s*-1*.

6 Public release of our ePSF and GD models

We make our library ePSF models and GD corrections publicly available to the community161616The ePSF models and GD corrections in format #1 are available at STDPSFs and STDGDCs, respectively. The GD corrections in format #2 are available at OTHERGDCs. as data-cube FITS files. The format of the ePSF models is described at the end of Sect. 3. The GD-correction files are provided in two formats. The FITS files in format #1 contain five layers. The value of each pixel in layer (1) gives the $x$ coordinate of each pixel in the 2048$\times$2048-pixel2 chip into the GD-corrected reference frame. This value is provided only for the center of the pixel; the GD correction in any subpixel position should be derived using a bi-linear interpolation. Layer (2) is similar but for the $y$ coordinate. Layer (3) contains the pixel-area photometric corrections to be added to the magnitudes measured in the level-2 _cal images. Layers (4) and (5) are the analogs of layers (1) and (2) but for the inverse GD correction, to map pixels from the GD-corrected frame into the raw NIRISS frame. We refer to Appendix G of Anderson (2022a, b) for a detailed description of the GD-correction files in this format. The GD-correction FITS files in format #2 are similar to those in format #1, but they instead provide the actual values of the GD correction in layers (1) and (2), and again the pixel-area corrections in layer (3).

Recently, Anderson (2022a, b) released a FORTRAN routine, hst1pass, designed for ePSF-fitting of undersampled HST images. hst1pass can provide coordinates in the raw system of a given HST image, but it can also output GD-corrected positions. A version of hst1pass adapted to handle \jwstimages, jwst1pass, has been used in the analysis described in this study and a preliminary version of this code is released along with this paper171717The FORTRAN code is available at CODE. and will be described to the community in a future paper (Anderson et al., in preparation). Our ePSFs and GD corrections in format #1 are designed to be used directly by jwst1pass, and users should be able to use them in their analysis of NIRISS images following the prescriptions in Anderson (2022a, b).

Our ePSFs can potentially be used with other software packages, e.g., DOLPHOT (Dolphin, 2000, 2016) or photutils (Bradley et al., 2022), so long as these packages follow the same ePSF treatment and conventions described at the end of Sect. 3 or convert the ePSFs into a format that they can use. Similarly, our GD-correction FITS files in format #2 can be easily used with any program that can handle the FITS standard, including python.

7 Discussion and Conclusions

Proper motions have long been critical to our understanding of stellar systems and structures in the solar neighborhood. Starting in the early 2000’s, HST greatly advanced this field thanks to improved imagers, observations covering increasingly longer temporal baselines and new techniques. This made it possible to address for the first time the transverse kinematics of a wide range of stellar systems and structures throughout the entire Local Group (van der Marel, 2015), while also enabling new studies of stellar and exoplanetary systems in the Solar neighborhood (Bedin et al., 2017; Fontanive et al., 2021). Gaia further revolutionized this field starting in 2016, by combining high astrometric accuracy and well-calibrated systematics with all-sky coverage (Brown, 2021). Despite the advances of Gaia, space observatories like HST continue to be unique for a wide range of studies, even given their relatively small field of view. Thanks to the larger mirror size and the ability to integrate for long times on specific fields, the same astrometric accuracy can be achieved at much fainter magnitudes. This offers many advantages, including the ability to accurately determine PMs further out in the Local Group (e.g., Sohn et al., 2020). Moreover, the high spatial resolution makes it possible to resolve, e.g., the centers of star clusters that are too crowded for Gaia (Anderson & van der Marel, 2010; Bellini et al., 2014, 2017a; Watkins et al., 2015; Libralato et al., 2018a, 2019; Watkins et al., 2022; Libralato et al., 2022, 2023; Sabbi et al., 2022). Finally, the combination of data from HST and Gaia offers additional advantages (Massari et al., 2020; del Pino et al., 2022), in that it allows the intrinsically differential nature of small-field astrometry to be calibrated to an absolute system.

This is the first in a series of papers by our \jwstTelescope Scientist Team to extend the astrometric legacy and lessons learned from HST to \jwst. The \jwst observatory is uniquely suited for this, due to its high image quality and temporal stability (McElwain et al., 2023). This offers several exciting new future opportunities, including: (a) the ability to study much fainter objects, owing to large mirror size of \jwst compared to HST. In turn, this yields better statistics and accuracy in measurements of the mean and dispersion in kinematic properties, while also enabling PM studies of more distant objects; (b) the ability for PM studies of regions that are too extincted for HST or Gaia, due to \jwst’s infrared sensitivity; and (c) the ability to extend the time baseline of imagery with space-based resolution to an era when HST may no longer be available. All this will provide the opportunity to further advance the types of studies in which HST and Gaia have already excelled, while also opening the door for entirely new studies and discoveries.

Here we have focused on developing the tools necessary to obtain high-precision astrometry and photometry with the NIRISS instrument of \jwst. We described in great detail the careful analyses we performed to construct accurate ePSF models and GD corrections, and the rigorous testing we put our products through. This has also allowed us to obtain new insights on the NIRISS detector performance and on \jwstin general. We make our ePSF models and GD corrections publicly available, together with a preliminary version of the dedicated software we are developing to analyze data taken with all \jwst’s imagers, which will be described to the community in a future paper (Anderson et al., in preparation). Our combined efforts will provide to the community the tools necessary to achieve with \jwst the same outstanding results obtained for point sources that have been possible with HST data.

As a scientific application and validation of our work, we studied the population-dependent PM kinematics of the stars in the standard \jwstcalibration field near the center of the LMC disk. We determined the stellar PMs by combining second-epoch \jwstCommissioning data with a first-epoch HST archival catalog based on data taken 16 years prior. For stars with apparent magnitudes around ${\it G}\sim 20$, the median PM uncertainty is $\sim$13 $\mu$as yr*-1*, corresponding to 3.1 km s*-1* at the distance of the LMC. This is 1.2–1.8 times better than Gaia DR3 achieves for its best-measured stars (${\it G}\sim 9$–12), and 3–5 times better than what it achieves at $G\sim 16$, around the tip of the LMC RGB.

These stunning accuracies are further validated by our ability to kinematically resolve tight substructures in the VPD of young and intermediate-age blue stars. We kinematically detect a known star cluster, OGLE-CL LMC 407. We measure for the first time its absolute PM and show how this differs from other blue LMC disk populations in the field. We also measure the cluster velocity dispersion to be $23.7\pm 1.8$ $\mu$as yr*-1*, corresponding to $5.6\pm 0.4$ km s*-1*. This provides an upper limit on the possible presence of any remaining uncorrected systematic errors. These results provide a testament to the overall quality of our \jwstNIRISS ePSF and GD-correction models.

The RGB stars in the LMC field show no obvious kinematic substructures, and have a velocity dispersion of $33.8\pm 0.6$ km s*-1* ($0.144\pm 0.003$ mas yr*-1*). This is comparable to results obtained from a large body of line-of-sight velocity studies, and reasonable by comparison to the central kinematics of such stars in the more massive Milky Way. However, this dispersion is 2–2.5 times less than what has been reported from studies of the Gaia DR3. This is the case for most of the younger blue populations in the field as well. We hypothesize that these discrepancies could be due to uncorrected increases in Gaia PM scatter due to crowding and blending near the dense LMC center.

In summary, our study successfully validates that \jwstis an extremely capable observatory for high-accuracy astrometry and PM studies in the Local Group. Unlike with HST, where it took more than a decade for its full astrometric potential to be realized, this capability is available and accessible for the astronomical community with \jwstright from the start. We therefore anticipate that PM studies will become a key new component of \jwst’s already vast arsenal of observational capabilities to further revolutionize astronomy over the coming decade.

Acknowledgments

The authors thank the anonymous referee for the thoughtful suggestions that improved the quality of our paper. This paper reports work carried out in the context of the JWST Telescope Scientist Team181818https://www.stsci.edu/$\sim$marel/jwsttelsciteam.html. (PI: M.Mountain). Funding is provided to the team by NASA through grant 80NSSC20K0586.

NHA acknowledges support by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1746891. DRL acknowledges research support by an appointment to the NASA Postdoctoral Program at the NASA Goddard Space Flight center, administered by Oak Ridge Associated Universities (ORAU) under contract with NASA. AG acknowledges support from the Robert R Shrock Graduate Fellowship. KKWH acknowledges funding from the Giacconi Fellowship at the Space Telescope Science Institute. JG acknowledges support from SERB through SERB Startup Research Grant SRG/2022/000727. RJM acknowledges support from NASA through the NASA Hubble Fellowship grant HST-HF2-51513.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

This research was pursued in collaboration with the HSTPROMO (High-resolution Space Telescope PROper MOtion) collaboration191919https://www.stsci.edu/$\sim$marel/hstpromo.html., a set of projects aimed at improving our dynamical understanding of stars, clusters and galaxies in the nearby Universe through measurement and interpretation of proper motions from HST, JWST, Gaia, and other space observatories. We thank the collaboration members for the sharing of their ideas and software.

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. Also based on observations with the NASA/ESA HST, obtained at the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. 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).

Appendix A \jwstdata processing information

Table A1 provides information about the processing of the \jwstdata used in this work. We report the CRSD context (CRDS_CTX), the version of the calibration pipeline (CAL_VER), the version number of the data processing software (SDP_VER) and the Science and opertaion (S&OC) Project Reference Database (S&OC PRD) version number used in data processing (PRD_VER).