The JWST Early Release Science Program for Direct Observations of Exoplanetary Systems I: High-contrast Imaging of the Exoplanet HIP 65426 b from 2 to 16 μm

We present JWST Early Release Science coronagraphic observations of the super-Jupiter exoplanet, HIP 65426b, with the Near-Infrared Camera (NIRCam) from 2 to 5 μm, and with the Mid-Infrared Instrument (MIRI) from 11 to 16 μm. At a separation of ∼0.″82 (87 −31+108 au), HIP 65426b is clearly detected in all seven of our observational filters, representing the first images of an exoplanet to be obtained by JWST, and the first-ever direct detection of an exoplanet beyond 5 μm. These observations demonstrate that JWST is exceeding its nominal predicted performance by up to a factor of 10, depending on separation and subtraction method, with measured 5σ contrast limits of ∼1 × 10−5 and ∼2 × 10−4 at 1″ for NIRCam at 4.4 μm and MIRI at 11.3 μm, respectively. These contrast limits provide sensitivity to sub-Jupiter companions with masses as low as 0.3M Jup beyond separations of ∼100 au. Together with existing ground-based near-infrared data, the JWST photometry are fit well by a BT-SETTL atmospheric model from 1 to 16 μm, and they span ∼97% of HIP 65426b's luminous range. Independent of the choice of model atmosphere, we measure an empirical bolometric luminosity that is tightly constrained between logLbol/L⊙ = −4.31 and −4.14, which in turn provides a robust mass constraint of 7.1 ± 1.2 M Jup. In totality, these observations confirm that JWST presents a powerful and exciting opportunity to characterize the population of exoplanets amenable to high-contrast imaging in greater detail.


INTRODUCTION
Across the last twenty-five years a variety of observational techniques have been employed to uncover and characterise the current population of over 5000 confirmed exoplanets (NASA Exoplanet Science Institute 2020).Of these techniques, the direct detection of photons from an exoplanetary atmosphere -direct imaging -remains one of the most challenging due to the substantial contrast in flux between host stars and their exoplanetary companions.The emitted flux of an exoplanet can be many magnitudes fainter than the stellar diffraction halo at its angular separation, and bespoke instrumentation (e.g., Macintosh et al. 2014;Beuzit et al. 2019) and/or image post-processing (e.g., Chauvin et al. 2004;Marois et al. 2008) is needed to isolate the exoplanet emission.Even with state-of-the-art instruments and data analysis procedures, only ∼20 planetary mass companions (PMCs) have been detected and characterised through these "high-contrast" observations (Currie et al. 2022), and all exoplanets directly imaged to date have estimated or dynamically-measured masses 2 M Jup .Despite these drawbacks, high-contrast observations offer considerable advantages to other techniques.At present, high contrast imaging is the most viable technique for characterising the population of exoplanets at orbital separations greater than ∼10 au.Furthermore, beyond the large population of irradiated, close-in planets with atmospheric measurements obtained via exoplanet transit observations, direct observations of exoplanet emission remain the most readily accessible path towards the characterisation of exoplanet atmospheres.Constraints on atmospheric composition may improve our understanding of exoplanet formation and evolution (e.g., Öberg et al. 2011;Morley et al. 2019;Zhang et al. 2021;Mollière et al. 2022), although these determinations can be highly dependent on the post-formation accretion of solid material.Compared to close-in transiting exoplanets, directly imaged planets present a distinct advantage in this regard, as they are easier to detect at younger ages where they are less likely to have experienced significant migration and/or accretion.Additionally, at young ages bulk properties such as temperature, radius, and bolometric luminosity provide independent constraints on formation conditions (Marley et al. 2007;Marleau & Cumming 2014) that can be contrasted to atmosphere-driven conclusions on formation.Finally, the study of exoplanet atmospheres continues to advance towards smaller, and more Earth-like exoplanets, and could ultimately lead to the discovery of life outside our Solar System (Schwieterman et al. 2018).
Nevertheless, to fully realise the advantages of highcontrast imaging, upgraded or new observatories and instruments (e.g., Gardner et al. 2006;Males et al. 2018;Chilcote et al. 2020) will be necessary so that we can expand this technique across a broader wavelength range, and to a wider diversity of closer separation and/or lower mass exoplanets.

High-Contrast Observations with JWST
With a primary mirror diameter of 6.5 m, an operational wavelength range from 0.6−28.1 µm, and a diverse range of instrumental modes, JWST (Gardner et al. 2006) presents a revolutionary opportunity for scientific exploration and discovery across many branches of astronomy and astrophysics.Within this remit, highcontrast observations of exoplanets and exoplanetary systems are no exception.
JWST is located at the second Sun−Earth Lagrange point, far from the thermal background, telluric contamination, and wavefront aberrations generated by Earth's atmosphere.The combination of excellent optical performance (∼75−130 nm RMS wavefront error, depending on the instrument), highly stable wavefront (<2 nm drift over a few hours), and large telescope aperture, enables JWST to reach better photometric and spectroscopic limiting sensitivities than both past or current ground-(i.e., 8−10 m class telescopes) and space-based (e.g., Hubble, Spitzer ) observatories (Rigby et al. 2022).It is not only this increased sensitivity that improves our ability to detect and characterise faint objects such as exoplanets, but its combination with JWST 's access to the near-and mid-infrared.At these wavelengths, the flux emitted from a hotter host star steadily decreases as a function of increasing wavelength, whereas the flux emitted from cooler exoplanetary companions reaches a peak.Hence, the natural star-planet contrast is minimised.To realise these advantages, JWST offers a selection of instrumental modes designed for, or that can be applied to, high-contrast observations.Specifically, both NIRCam (Rieke et al. 2005) and MIRI (Rieke et al. 2015) are equipped with coronagraphic masks (Krist et al. 2009;Boccaletti et al. 2015), NIRISS (Doyon et al. 2012) is equipped with a non-redundant mask which enables aperture masking interferometry (AMI; Sivaramakrishnan et al. 2012), and although lacking any hardware for starlight suppression, both NIRSpec (Bagnasco et al. 2007) and MIRI are equipped with integral field units (IFUs; Wells et al. 2015;Böker et al. 2022).
In anticipation of the range of capabilities that JWST would provide, a similar range of predictions and simulations were constructed in an effort to forecast its potential for exoplanet imaging science.With JWST 's extraordinary sensitivity across 4−15µm (where cooler planets are more luminous, Morley et al. 2014), the first direct detections of sub-Jupiter, and even sub-Saturn, mass planets will be within reach (Brande et al. 2020;Beichman et al. 2020;Carter et al. 2021a;Ray et al. 2023).Already discovered companions will also be readily detectable across this broad wavelength range, allowing for deeper atmospheric characterisations that may result in the detection of a range of atmospheric species (Danielski et al. 2018;Patapis et al. 2022), or tighter constraints on bulk planetary properties.
Nevertheless, these predictions are based on groundbased testing and observatory simulations, whereas the true capabilities of JWST hinge on its on-sky performance.Preliminarily, the performance of both the NIR-Cam and MIRI coronagraphic modes exceeded expectations during observatory commissioning (Girard et al. 2022;Kammerer et al. 2022;Boccaletti et al. 2022), but the first scientific demonstrations of JWST 's capabilities are being conducted as part of the Director's Discretionary Early Release Science (ERS) Programs1 .
Our ERS program "High Contrast Imaging of Exoplanets and Exoplanetary Systems with JWST " (ERS-01386, Hinkley et al. 2022a) is the only ERS program that has tested the high-contrast exoplanet imaging modes of JWST and includes: coronagraphic imaging from 2−16 µm of the known exoplanet HIP 65426 b (Chauvin et al. 2017, this work) and circumstellar disk HD 141569 A (Weinberger et al. 1999, Millar-Blanchaer et al. in preparation, Choquet et al. in preparation), spectroscopy from 1−28 µm of the PMC VHS J125601.92-125723.9AB b (VHS 1256 b, Gauza et al. 2015;Miles et al. 2023), and AMI observations of HIP 65426 at 3.8 µm (Sallum et al. in preparation, Ray et al. in preparation).This program is rapidly disseminating these crucial initial data, and demonstrating the true capabilities of JWST for highcontrast imaging and spectroscopy for the first time.Furthermore, we will provide a range of science enabling products (e.g., data analysis pipelines, recommendations for best observing practices) to the community to support their own proposals and investigations in Cycle 2 and beyond (Hinkley et al. 2022a).
In this work we focus exclusively on the coronagraphic imaging observations of the HIP 65426 system within this ERS program, and their context within a broader understanding of JWST as a tool for high-contrast imaging.

HIP 65426 b
Discovered by Chauvin et al. (2017), HIP 65426 b is a super-Jupiter mass exoplanet at a wide physical separation of 110 +90 −30 au (Cheetham et al. 2019) to the star HIP 65426 (A2V, 2MASS J=6.826, J−K=0.055,M =1.96±0.04).HIP 65426 is located at a distance of 107.49±0.40pc (Gaia Collaboration 2020), has no signs of binarity from radial velocity and sparse aperture masking observations (Chauvin et al. 2017;Cheetham et al. 2019), and is a fast rotator (vsin(i)=299±9 km s −1 , Chauvin et al. 2017).Furthermore, HIP 65426 is a likely member of the Lower Centaurus-Crux association as derived from its proper motion and radial velocity measurements (89% probability, Gagné et al. 2018), constraining its age to 14±4 Myr.The interest in this association has grown over time with the increasing number of directly imaged exoplanet discoveries within it (e.g., HD 95086 b Rameau et al. 2013, PDS 70 b/c Keppler et al. 2018, and TYC 8998 b/c Bohn et al. 2020. Although HIP 65426 b was initially observed with a combination of low-resolution spectroscopy and photometry from ∼1−2 µm (Chauvin et al. 2017), follow up observations expanded this coverage to ∼1−5 µm (Cheetham et al. 2019;Stolker et al. 2020b), including a medium-resolution spectrum from ∼2−2.5 µm (R 5500, Petrus et al. 2021).Photometric analysis has demonstrated that HIP 65426 b is similarly located to mid-to-late L-dwarfs in colour-magnitude diagrams (Figure 1), and lies between already discovered early L-dwarf exoplanet companions (e.g., β Pic b, HD 106906 b) and those at the L/T transition (e.g., HR 8799 c/d/e).Using combined photometric and spectroscopic observations, Petrus et al. (2021) performed an atmospheric forward model analysis of HIP 65426 b, indicating that it has T eff =1560±100 K, log(g)≤4.40dex, [M/H]=0.05+0.24  −0.22 dex, and the atmospheric carbon-tooxygen ratio has an upper limit of C/O≤0.55.Furthermore, Petrus et al. (2021) also detect the presence of H 2 O and CO in the atmosphere of HIP 65426 b using cross-correlation molecular mapping, in addition to non-detections of CH 4 and NH 3 .Finally, from evolutionary model analyses to the data reported in Chauvin et al. (2017), Marleau et al. (2019) estimate the mass of HIP 65426 b to be 9.9 +1.1 −1.8 M Jup or 10.9 +1.4 −2.0 M Jup for hot-or cold-start initial entropy conditions, respectively (see e.g., Marley et al. 2007;Spiegel & Burrows 2012).
In Section 2 we describe our JWST observations of HIP 65426 b and all necessary data reduction.In Section 3 we describe the analysis steps taken to produce residual starlight subtracted images and measurements of contrast performance.We present a discussion of these observations in the context of both the overall performance of JWST in Section 4, and our understanding of HIP 65426 b in Section 5. Finally, we summarise our conclusions in Section 6.

OBSERVATIONS & DATA REDUCTION
The NIRCam and MIRI coronagraphic imaging observations of HIP 65426 b presented here were taken as part of program ERS-01386 (Hinkley et al. 2022a) and exist as a subset of a broad range of observations to assess the performance of JWST 's high-contrast imaging and spectroscopic modes with respect to the study of exoplanetary systems (Hinkley et al. 2022a).

Observational Structure
The observational strategies used for this program were adopted following the recommended best practices as known prior to launch and described in the JWST user documentation2,3 .All observations of HIP 65426 b are repeated at two independent roll angles separated by ∼10 • to enable subtraction of the residual stellar point spread function (PSF) through angular differential imaging (ADI, Müller & Weigelt 1985;Liu 2004;Marois et al. 2006).Although a large number of rolls across a larger angular range would be desirable for an optimal subtraction using this technique, the combina-tion of lengthy exposure times, increased overheads, and spacecraft orientation constraints prohibit an observing strategy more complex than described.Given the maximum possible roll offset for JWST at any given epoch is 14 • , a larger roll offset than we have adopted would also require multi-epoch observations.
We also perform similar observations of a bright reference star, HIP 68245 (B2IV, 2MASS J=4.628, J−K=0.137), to additionally enable subtraction through reference differential imaging (RDI).This star was selected as it is: a) bright, therefore reducing the exposure time required to attain a similar signal-to-noise as the science target, b) a similar spectral type to the science target, therefore reducing the impact of spectral type mismatch, c) is relatively close (∼10 • ) to the science target, therefore reducing slew overheads and minimising position dependent wavefront drift between science and reference observations, and d) has no evidence of binarity as determined by VLT /SPHERE AMI observations (Proposal ID: 108.22CD).For information on selecting a suitable reference star, see the JWST User Documentation4 .All exposure settings for the science and reference observations are shown in Table 1.
An RDI based subtraction is likely to reach superior contrast limits to ADI from pre-launch predictions (Lajoie et al. 2016;Perrin et al. 2018;Carter et al. 2021b), and is therefore also more representative of the optimal performance of JWST coronagraphy (also see Section 4.1).The exposure settings for these reference observations were chosen to reach an approximately equivalent fraction of full well saturation per integration as the corresponding target observations.Additionally, each reference observation was repeated at nine separate dither positions following small-grid dither patterns 9-POINT-CIRCLE and 9-POINT-SMALL-GRID for the NIRCam and MIRI observations, respectively (Soummer et al. 2014;Lajoie et al. 2016).The goal of this strategy is to produce a small library of reference PSFs for each science exposure which captures different misalignments between the star and the center of the coronagraphic mask, and can in turn facilitate more advanced PSF subtraction techniques (e.g., KLIP, Soummer et al. 2012).Further discussion on the relative benefits between ADI and RDI subtraction strategies, or a combination of the two, with respect to these JWST observations can be found in Section 4.1.
For MIRI, we also add background observations to both our science and reference observations in both filters of a nearby "empty" region of sky (as identified in WISE images, Wright et al. 2010) separated ∼1.5 away the target star to measure the stray light "glow stick" that is inherent to MIRI coronagraphic observations (Boccaletti et al. 2022).Specifically, this position corresponds to a right ascension and declination of α=13 • 24 44.2915 , δ=−51 • 29 31.54, respectively (ICRS J2000 coordinates).To best match the science and reference observations, the exposure parameters for each background observation exactly match the parameters for a single roll/dither of the associated science or reference target.These observations were intended to be performed at two separate dither positions to identify astrophysical sources that might impact the background subtraction, however, due to a previously unresolved issue they were instead repeated at an identical pointing (Dean Hines, private communication).
The NIRCam and MIRI observations were executed as two separate non-interruptible sequences, ensuring that observations between rolls, and also between science, reference, and background targets, are minimally separated in time.This reduces the extent to which the wavefront can vary across observations, due to variations in the telescope mirror alignment, the thermal evolution of the telescope, or both (Perrin et al. 2018;Carter et al. 2021b).Changes in the wavefront will lead to variations in the residual PSF between exposures, hinder our ability to perform an optimal subtraction of these residual PSFs, and suppress the overall achievable contrast.

spaceKLIP
For all observations we perform data reduction using the newly developed and publicly available python package, spaceKLIP5 (Kammerer et al. 2022).Briefly, spaceKLIP takes a collection of data products from the jwst pipeline6,7 (Bushouse et al. 2022) as inputs, and generates PSF subtracted images, contrast curves, and measurements of companion photometry and astrometry.The majority of this functionality is provided by the underlying pyKLIP (Wang et al. 2015) package, with spaceKLIP providing a user friendly interface, streamlined code execution, custom JWST data reduction routines, and built-in plotting procedures.

NIRCam Coronagraphy
The NIRCam observational sequence was executed from 23:00 July 29th to 05:16 July 30th 2022 UTC, with exposures taken using the MASK335R round coronagraphic mask (Krist et al. 2010) in the F250M, F300M, F410M, F356W, and F444W filters for the reference star (PA=110.2• ), then HIP 65426 (PA=110.0• ), and then finally HIP 65426 at a second roll angle (PA=120.4• ).This sequence structure significantly reduces overheads, as once the target acquisition has been performed it is not necessary to reacquire the target to switch the observational filter.We begin data reduction using the Stage 0 (*uncal.fits)files as generated by the jwst pipeline.These products are then processed to Stage 1 (*rateints.fits)files using spaceKLIP, which follows a slightly modified version of the jwst pipeline.Where possible the JWST detectors have reference pixels that can be used to track and correct drifts in the measured pixel counts due to readout electronics.In the absence of such a reference, these drifts may instead be misinterpreted as "jumps"8 from cosmic ray events during the up-the-ramp (MUL-TIACCUM) detector readout.The NIRCam coronagraphic subarrays do not have any embedded reference pixels and default pipeline processing leads to multiple erroneous jump detections and greatly increased noise in processed images.Therefore, we manually define all pixels within a four pixel border of each image as refer-ence pixels within the pipeline to mitigate these effects.Additionally, we identify a significant improvement in image quality by skipping the dark current subtraction step that is turned on as a default in the jwst pipeline.At present, the dark current calibration data exhibit a large number of hot pixels, persistence, and cosmic rays which cannot be averaged out or corrected due to limited number of available integrations.Attempts to perform a dark current subtraction using this data result in a variety of negative flux residuals in the reduced images and an overall reduced sensitivity of the final pipeline product.Once this calibration file is improved through further calibration observations, we do not anticipate a need to skip this step for NIRCam coronagraphic data.
During the jump detection step, the jwst pipeline will make use of a detection threshold value based on the estimated signal and noise to assess whether a deviation between groups is significant enough to be considered a jump (see 6 for further detail).The default value for this threshold is 4, but we repeat an early version of our F444W analysis across thresholds of 4 to 16 to search for potential improvements.We find that the contrast is slightly improved from a threshold of 4 to 5 by ∼5×10 −8 at 1 , but at larger thresholds it does not vary (deviations <1×10 −8 at 1 ).Hence, we adopt a detection threshold of 5 for all of our NIRCam analyses.
The Stage 1 products are processed further to Stage 2 (*calints.fits)files using spaceKLIP, with some addi-tional pixel cleaning procedures as follows.Firstly, every pixel with a data quality flag (e.g., indicating hot or warm pixels, unreliable data processing) is replaced by the median of its orthogonal and diagonal neighbours, with the notable exclusion of pixels with a jump flag which are typically grouped in clusters.We also inspect each pixel for temporal flux variations and identify situations where: a) the pixel is bright (MJy/sr > 1) for at least one integration, and b) the pixel is relatively faint (a >80% decrease in flux compared to the brightest integration) for at least one integration.The pixel values for the integrations that are not marked as faint are then replaced by the median value of the integrations that are marked as faint.We note here that an outlier identification process based on variations from the standard deviation of each pixel in time may be preferable, but is difficult to incorporate given the small number of integrations in these exposures (see Table 1).Despite the above corrections, ∼25 static hot pixels remain across all images.Although these pixels do not impact our ability to recover HIP 65426 b, they introduce residuals in the PSF subtraction process (see Section 2.6) that bias our measurements of the contrast performance.Therefore, we provide the locations of these pixels to spaceKLIP manually and correct them in an identical manner as the pixels marked with data quality flags.
As the NIRCam PSF in the F250M filter is undersampled, ringing artifacts are generated by the interpolation methods in pyKLIP used for image registration and spatial shifting of input PSFs.These artifacts bias the image registration process, and limit our ability to accurately inject synthetic PSFs for contrast curve calibration and companion astrometry and photometry.To overcome this issue, we smooth all of the F250M images by a Gaussian filter, as implemented by scipy, with σ=1.3 pixels.This value of σ was chosen as it was the minimum possible value that removed the observed artifacts across a test of ten equally spaced values from σ=1 to σ=2.We note that this factor of 1.3 is consistent with the ratio of the detector pixel scale and the theoretical Nyquist sampling at 2.5 µm (assuming a reduced primary mirror diameter of 5.2 m due to the NIRCam Lyot stop).Smoothing these images may lead to reduced precision in our astrometric analysis, however, it should not influence the accuracy of our retrieved photometry.

MIRI Coronagraphy
The MIRI observational sequence was executed from 21:05 July 17th 2022 to 05:19 July 18th 2022 UTC.Exposures were taken for HIP 65426 in the F1140C filter (PA=117.4• ), then once again at a second roll angle (PA=108.0• ), and then for the reference star (PA=109.2• ).This sequence was then repeated in reverse for the F1550C filter, except with slightly different target roll angles of 108.2 • and 117.4 • .This structure is different to that of the NIRCam observations as each filter is tied to a specific four-quadrant phase mask coronagraph (FQPM, Boccaletti et al. 2015), and target acquisition must be repeated when switching between them.Inserting the reference observations between the observations of HIP 65426 b minimises the time separation between science and reference exposures for each filter, and therefore the extent of the wavefront evolution between them.After all science and reference observations are complete, we perform the dedicated background observations that are used to subtract the dominant "glow stick" stray light feature (Boccaletti et al. 2022) as described in Section 2.1.
We begin data reduction using the Stage 0 (*uncal.fits)files as generated by the jwst pipeline.These products are then processed to Stage 1 (*rateints.fits)files using spaceKLIP.Similarly to NIRCam (see Section 2.3) we explore the impact of the jump detection threshold on our analysis.For these data in particular, we found that the default jump detection threshold value of 4 is too low and leads to a number of pixels being erroneously flagged as containing a jump.Flagged pixels are interpreted differently to unflagged pixels in the ramp fitting procedure and the resulting Stage 1 files contain a large number of pixels with unrealistic (negative) flux values as a result.After repeating an early version of our F1140C analysis across thresholds of 4 to 16, we observe an improvement in contrast between a threshold value of 4 and 5 of ∼2×10 −5 at 1 , and ∼1×10 −5 at 2 .Beyond this value there is only slight improvement, and the obtained contrast varies by less than ∼5×10 −6 at 1 .For our final analyses we select a threshold value of 8, as it has the best contrast at 2 and fewer pixels with unrealistic flux values than lower thresholds (as determined from visual inspection).
Following ramp fitting we found that the first integration of each exposure contained a significantly increased level of non-uniform noise indicative of a detector reset charge decay anomaly, which is driven by differences in how the MIRI field effect transistors reset the detector charge prior to an exposure versus between integrations 9 .Ideally this anomaly can be corrected by calibration dark exposures, however, the currently available calibration darks were acquired in quick succession, and do not have a similar amount of dead time (e.g., due to telescope slews/dithers) as our science exposures.As a result, during our observations the detector electronics were given a much longer time to settle, and our integrations exhibit an entirely different reset anomaly.This effect is most dominant for the first integration of each exposure, but appears to persist throughout the entire exposure as well.Quantitatively, the median flux of the first integration following background subtraction (see below) for the target/reference observations are 6.2/3.0σ and 7.3/3.9σdeviant from the average median flux across all integrations for the F1140C and F1550C observations, respectively.In comparison, the median flux of all other target/reference integrations deviates within ranges of 0.003−0.59/0.26−0.55σand 0.004−1.29/0.06−0.86σfor the F1140C and F1550C observations, respectively.As this increased noise is not accurately captured in the available calibration files and differs significantly from all other integrations, we opt to exclude the first integration of each exposure from all further analysis.This exclusion corresponds to a 2.5%/10% and 1.7%/5.3%cut of the target/reference data for the F1140C and F1550C observations, respectively.In the future better dark exposures will be taken as part of observatory calibrations that may well remove the need to exclude the first integration, and we encourage future observers to carefully evaluate the calibration status of their data before adopting similar cuts.
The Stage 1 products are processed further to Stage 2 (*calints.fits)files using spaceKLIP, with some additional pixel cleaning procedures as follows.First, every pixel with a data quality flag (e.g., indicating hot or warm pixels, unreliable data processing) is replaced by the median of its orthogonal and diagonal neighbours, with the notable exclusion of pixels with a jump flag which are typically grouped in clusters.Following this correction, ∼30 static hot pixels remain in our images and are corrected following an identical procedure by manually providing the pixel locations to spaceKLIP.Similarly to NIRCam, these final pixels primarily impact the measured contrast performance, and not our ability to recover HIP 65426 b.
As shown in Boccaletti et al. (2022), the MIRI coronagraphic fields of view are subjected to a stray light "glow stick" feature along the horizontal edges of the FQPMs that dwarfs the residual stellar flux (see Figure 2).We subtract this feature from our processed Stage 2 products for each filter using a median background image of every 4−5 integrations from the dedicated background observations to the corresponding 4−5 science or reference integrations.The value of 4−5 was selected as it provided a slightly improved contrast at 1 compared to other tested numbers of integrations per median, rang- ing from 1 (∼1×10 −4 improvement) to all available integrations (∼1×10 −5 improvement).Grouping the median subtraction in this manner better captures the diffuse reset anomaly noise between integrations mentioned above.Additionally, it may help capture variations in the stray light feature, which varies with a standard deviation of ∼0.5% (as identified by variations in the median pixel flux for pixels above 50% of the peak pixel flux across integrations).
Following the median background subtraction we find that the residual stellar flux is easily recovered (see Figure 2).We also attempted to model the stray light glow stick as a principal component within both KLIP (Soummer et al. 2014) and LOCI (Lafrenière et al. 2007) based subtractions, but found that they were susceptible to oversubtraction of the residual stellar flux and/or could not additionally account for background variations between integrations.We anticipate that with careful masking and optimisation of these algorithms it may be possible to overcome these issues, however, the median frame background subtraction is already highly effective and improvements to the achieved contrast are likely to be minimal.

Image Alignment
The NIRCam and MIRI coronagraphic modes adopt independent target acquisition procedures to correctly center a star behind each focal plane mask.The in-flight positions of the NIRCam mask centers are known to better than ∼10 mas but the distortion model is still being refined.At present the target acquisition error for the MASK335R is as large as ∼12−30 mas, or 0.2−0.5 pixels (Girard et al. 2022).This error is dominated by the precision of the centering algorithm, which is not well adapted to the PSF shape with coronagraphic optics (wider in x−axis), and not by the small angle maneuver (SAM) that places a target at its desired position behind the mask (which for NIRCam is repeatable to ∼6 mas).For MIRI the mask center positions have been measured to ∼5−10 mas, or ∼0.1 pixels (Boccaletti et al. 2022).However, the SAM has a typical uncertainty of ∼ 10 − 20 mas, leading to positional offsets between different rolls or targets.Finally, between integrations the pointing stability of JWST (∼1 mas) and the accuracy of the small-grid dither manoeuvres (∼2−4 mas) will lead to further positional shifts for both the NIRCam and MIRI coronagraphs (Rigby et al. 2022;Boccaletti et al. 2022).
For NIRCam the absolute star position is only explicitly measured for the first science image in each filter.This position is measured using a cross correlation of a model coronagraphic PSF as obtained from webbpsf ext to the science PSF using the scikit-image package (van der Walt et al. 2014).To best match the science PSF, the model PSF is generated using the telescope optical path difference (OPD) map as obtained on July 29th 2022.To identify the accuracy of this process we repeat the procedure, except comparing the model PSF to itself, and to a second model PSF that was generated using a pre-launch measurement of the telescope OPD (which exhibits comparable differences to the model PSF as our data).In each case, we manually shift the comparison PSFs across a range of 0.01 to 0.5 pixels and attempt to recover these offsets using the cross correlation process.For the self comparison, we recover the injected shift to at least ∼0.03 pixels, or 1.9 mas, whereas for the second model PSF comparison, we recover the injected shift to ∼0.1 pixels, or 6.3 mas.Given the comparable difference between the model PSF for both the second model PSF and our data, we adopt the latter uncertainty as a systematic uncertainty in our astrometric measurements (see 5.1).The shifts of the other science and reference images are obtained relative to the first science image through a similar cross correlation procedure.However, we instead use a box of 11×11 pixels around the coronagraph center position, where the flux is dominated by the central core of the coronagraphic PSF.All measured relative shifts match expectations for the JWST pointing precision of ∼1 mas (1σ, radial) (Rigby et al. 2022).
To estimate the absolute star position for MIRI, we use the current measurements for the centers of the coronagraphic masks and assume the star is perfectly centered behind the coronagraphic mask.In zero-indexed subarray x-y coordinates these values are (119.749, 112.236) and (119.746, 113.289) for the FQPM1140 and FQPM1550, respectively (Jonathan Aguilar, private communication).As a result, the stellar position is not known better than a minimum of ∼10 mas (Boccaletti et al. 2022).Similarly, the relative alignment between images may be discrepant by ∼1−4 mas (Rigby et al. 2022), or <0.05 pixels.Attempts to estimate the absolute star position in a similar manner to the NIRCam data were unsuccessful, and this process led to significantly larger estimated shifts than the known pointing stability of JWST (∼50 mas, vs ∼1 mas).This is most likely due to the variable spatial structure of the residual PSF, which significantly changes with small pointing shifts.An effort was also made to fit each MIRI image with model PSFs across all feasible pointings using webbpsf ext10 within an emcee (Foreman-Mackey et al. 2013) Markov chain Monte Carlo (MCMC) framework.However, these fits were unable to converge and upon visual inspection the spatial structure of the empirical and model PSFs were different, despite modelling the PSFs based on measurements of the telescope OPD within ∼1 day of our observations.We do not believe that using the coronagraphic mask centers as a proxy for the absolute star position has significantly affected our results, but it is certainly an area of improvement for future studies using MIRI coronagraphy.
For NIRCam all images are aligned to a common center based on the measured shifts, however, for MIRI we opt to not perform any realignment.This decision is made under the assumption that because a pointing shift does not primarily cause a translation of the residual PSF in the MIRI images (in contrast to NIRCam), the unshifted reference images are more descriptive of variations in the science images.When comparing a RDI subtraction using unshifted reference images, and a separate RDI subtraction with a realignment based on the ideal small-grid dither positions, the measured contrasts are in agreement and this choice does not significantly impact our results.

PSF Subtraction
We perform a subtraction of the residual stellar PSF in each filter following three different principal component analysis (PCA) based methods as implemented in spaceKLIP.First, we take the two independent rolls of HIP 65426 b and perform an ADI subtraction.Second, we perform an RDI subtraction by using the corresponding observations of the reference star, HIP 68245, as a PSF library.Finally, we perform an ADI+RDI subtraction, which is identical to the RDI subtraction except that images at the opposite roll angle are also included ) is also marked.Additionally, the intensity of all images for a given filter are identically scaled.The exoplanet, HIP 65426 b, can be easily identified at a position angle of ∼150 • in the subtracted images.We note that the distinct "hamburger" shaped central core and six-lobed structure of the companion PSF in the NIRCam images is an expected feature that is related to the Lyot stop design, and not indicative of discrete astrophysical sources.
in the PSF library.In each case the subtraction is performed on each integration from both science rolls individually, before being rotated to a common orientation as marked in Figure 3 and summed together.Although the number of annuli and subsections the PSF subtraction is performed across can be adjusted, we find that this does not improve the observed contrast.Hence, we perform all subtractions using a single annulus and a single subsection (i.e., the entire image).We leave future optimisation of these parameters for future analysis, but note that any improvements to the measured contrast are likely to be small as the noise in our images is close to azimuthally symmetric.The number of KLIP PCA modes can also be adjusted to tune the aggressiveness of the PSF subtraction.Hence, we perform the PSF subtraction across the full range of possible PCA modes to investigate the impact on our measured contrast and companion fitting.The maximum number of PCA modes is dependent on the exposure settings for each filter and corresponds to: the number of integrations in a single roll for ADI, the total number of integrations across all 9 dithers for RDI, and the sum of the two for ADI+RDI (see Table 1 for precise values).
Pre-and post-subtraction images for the NIRCam F444W and MIRI F1140C filters are shown in Figure 3, and images for all filters are shown in Appendix A. We note that the distinct "hamburger" shaped central core and six-lobed structure of the companion PSF in the NIRCam images is an expected feature that is related to the Lyot stop design, and not indicative of discrete astrophysical sources.

Contrast Calibration
All proceeding contrast measurements are determined relative to a synthetic spectrum of HIP 65426 in each of the JWST filters, as estimated from fitting stellar and disk models to existing photometry following Yelverton et al. (2019), see Figure 4. We use data of HIP 65426 from Hipparcos/Tycho-2 (Høg et al. 2000), Gaia DR2 (Gaia Collaboration 2018), 2MASS (Cutri et al. 2003), ALLWISE (Wright et al. 2010), AKARI IRC (Ishihara et al. 2010), and Spitzer MIPS (Chen et al. 2012).The fitting procedure compares synthetic photometry of models to the data to compute a χ 2 value, and posterior distributions are found using MultiNest (Feroz et al. 2009;Buchner et al. 2014).We derive our own zero points using the CALSPEC Vega spectrum (Bohlin et al. 2014).We use PHOENIX models (Allard et al. 2012) for the stellar photosphere, and a Planck function for the disk model.There is a small excess at 24 µm that was previously reported at 3.5σ by Chen et al. (2012), though not considered significant in that paper.The best fit model has an effective temperature of 8600±200 K and luminosity 16±1 L .The dust temperature and luminosity are very poorly constrained (T dust = 300 +200 −100 K, and L dust /L ∼ 2 × 10 −5 ), though this uncertainty does not significantly impact the flux estimation in JWST bandpasses because the excess is small.The flux excess at 24 µm is not high, but if real the dust would reside relatively close to the star, probably less than 1 au (hence we use the total model flux to compute the contrast, as any dust component that contributes IR flux would remain unresolved).We use the posterior distribution of model parameters and synthetic photometry to generate a distribution of fluxes in the JWST bands, and adopt the maximum likelihood solution for the stellar flux in each bandpass.

Contrast Curves
Following PSF subtraction, we obtain metrics of the sensitivity as a function of angular separation (i.e., contrast curves) for all observational filters using spaceKLIP.To avoid biasing the contrast measurement we mask regions of the subtracted images near HIP 65426 b, background sources, and the FQPM edges.These "5σ" contrast curves report the flux level corresponding to a 5σ-equivalent false alarm probability of 2.9 × 10 −7 after correcting for small sample statistics at small separations (Mawet et al. 2014).We call them 5σ contrast curves for brevity.
To obtain a more accurate measurement of the contrast performance, the throughput of the coronagraph and the intrinsic throughput of the KLIP subtraction must be accounted for by injecting and then recovering the flux of artificial sources (Adams & Wang 2020).All artificial sources are generated using webbpsf11 (Perrin et al. 2012(Perrin et al. , 2014) ) at an initial intensity equivalent to a signal-to-noise ratio of 25.Immediately prior to injection and based on a desired injection location, each source is modulated by the coronagraphic throughput using a synthetic throughput map.
Synthetic coronagraphic throughput maps are provided in the calibration reference files12 , however, both the provided NIRCam and MIRI FQPM maps are inaccurate.For the NIRCam MASK335R, the position of the occulting mask within the throughput map does not correspond with its actual location.Therefore, we modify the throughput map by extracting the pixels impacted by the occulting mask and repositioning them at the true mask center location of (149.9, 174.4) in zero-indexed subarray x-y coordinates (Jarron Leisenring, private communication).In the case of MIRI, all of the FQPM maps are rudimentary and do not accurately capture the spatial throughput variations.Therefore, we instead use custom simulated maps of the FQPM throughput produced using WebbPSF.In brief, for each coronagraphic mode we generate 1681 positiondependent PSFs across a 25.6" by 25.6" grid spanning the subarray field of view (FOV).The PSFs are generated with logarithmic spacing such that they are most densely sampled along the FQPM axes.For each position in the FOV, we calculate two PSFs, one occulted and one unocculted, and take the ratio of the integrated PSF fluxes in each case to provide a throughput estimate at each position.From this 2D sampling of throughput, we use scipy.interpolate.griddatato linearly interpolate the throughput estimates across the FOV to produce a smooth 2D map matching the subarray dimensions for each MIRI FQPM.
Due to the target acquisition errors described in Section 2.5, the measured star centers are offset from the the center of the coronagraphic occulter.However, in addition to these pointing errors, for NIRCam coronagraphy we expect wavelength-dependent spatial shifts of the entire image as viewed on the detector focal plane due to refraction through the coronagraphic mask sap-phire substrate along with deflections through filter optics further downstream in the optical train.Therefore, the star position in an image is not in isolation a suitable proxy for its position behind the coronagraphic mask.The wavelength-dependent shifts must be accounted for to more accurately determine the impact of the coronagraphic throughput across the detector focal plane and apply the correct throughput scaling to injected PSFs.
The aforementioned location for the MASK335R occulter at the detector focal plane is set based on the F335M filter, which is used for target acquisition and astrometric confirmation observations, and its image shift is defined to be zero.To measure the relative shifts of the remaining observational filters, we first determine the center of the projected stellar image for each filter observation of the science target roll positions and the reference target by cross-correlating an observed image in the Stage 1 *rate.fitsfile with a perfectly centered synthetic PSF generated with webbpsf.These synthetic PSFs are created using the on-sky OPD map from July 29, 2022, and were recentered to remove pre-flight model shifts, which do not fully capture contributions from all optical elements, such as the different filter optics.The measured sub-pixel locations for each observed filter are then compared to the similarly measured F335M astrometric confirmation image to determine the filterdependent offsets.The two independent target rolls and reference exposures provide three independent measures that are averaged together and presented in Table 2; uncertainties are on the order of 1-2 mas.Finally, to apply the correct coronagraphic throughput to an injected synthetic planet in a given filter at a given position, we simply realign the throughput map according to the measured offset for that filter.
Once scaled by the coronagraphic throughput, PSFs are injected into multiple copies of the unsubtracted science images across a range of separations extending to 4 , and for a range of position angles from 0 to 360 • .Sources are not injected within 2λ/D of each other or a masked region.These images then undergo KLIP subtraction in an identical manner to the science images, and the relative flux of an initial source PSF and the KLIP processed PSF as estimated within pyKLIP describes the overall coronagraphic mask plus KLIP throughput at each location.Finally, the basic 5σ contrast is divided by an interpolation of the median throughput across all position angles to obtain the calibrated contrast.
For the ADI subtraction, the measured contrast does not vary significantly when using more than 1 PCA mode for the NIRCam filters and more than 2 PCA modes for the MIRI filters.For RDI, the measured con- trast does not vary beyond ∼6 modes for both NIRCam and MIRI.Finally, for ADI+RDI we find that there is a transition to improved contrast at P max − P ADI modes, where P max is the maximum possible number of PCA modes possible, and P ADI is the maximum number of modes for the ADI subtraction.This is likely a result of the much larger number of reference images weighting the calculation of the PCA modes to be mostly RDI-like, until a sufficient limit is reached where the influence of the opposing roll images appears in the principal components.Beyond this transition value, the measured contrast does not vary significantly.
Example calibrated contrast curves for the NIRCam F444W and MIRI F1140C filters using the maximum number of PCA modes are shown in Figure 5, and contrast curves for all filters are shown in Appendix B.

HIP 65426 b PSF Fitting
To analyse the properties of HIP 65426 b in greater detail, we make use of the forward model PSF fitting routine provided by pyKLIP and implemented in spaceKLIP.Briefly, this routine takes a model of the companion PSF and uses a forward model of the KLIP subtraction to apply PSF distortions that arise naturally from the KLIP process.This resultant PSF can then be scaled/shifted to best match the observed companion PSF and obtain a measurement of its location and intensity (Pueyo 2016;Wang et al. 2016).
For our analysis we adopt an independent model PSF for each filter using webbpsf ext functionality as implemented in spaceKLIP.Specifically, we use webbpsf ext to generate an offset PSF at the predicted location of HIP 65426 b as adopted from whereistheplanet (Wang et al. 2021) and at the appropriate position for each science roll image.Each NIRCam and MIRI PSF is generated using a measured OPD map as determined from wavefront sensing and control observations on July 29th 2022 and July 17th 2022, respectively.Particularly for MIRI, the PSF of an off-axis source is still sensitive to its location relative to the FQPM edges, and explicitly generating a model PSF close to this location is important for obtaining a close match to the true companion  2019).We note that as the initial model PSF is normalised to unit intensity, the primary purpose of selecting a model SED is not to accurately predict the absolute flux of HIP 65426 b in each bandpass, but instead to capture the relative variation in flux as a function of wavelength across each bandpass (which is largest for the broader NIRCam wide filters).Finally, the impact of the coronagraphic throughput on the received flux is applied by multiplying the normalised PSF by a scale factor equal to the relative difference between the integrated flux of the model PSF prior to normalisation and a matching webbpsf ext model PSF which excludes the coronagraphic elements.
The input model PSF is converted to physical units by scaling it to the flux of HIP 65426 as estimated by the stellar model described in Section 3.1.Therefore, all derived photometry is anchored relative to our assumption of the stellar flux.Furthermore, any comparisons between the intensity of this PSF and the observed PSF is subject to an implicit assumption that the absolute flux calibration of JWST is perfect.In reality, the absolute flux calibration accuracy requirement for NIRCam and MIRI coronagraphy is 5% and 15%, respectively 13 .
13 https://jwst-docs.stsci.edu/jwst-data-calibration-considerationsBoth these fundamental uncertainties on the absolute flux calibration and the uncertainty on the stellar model flux as derived from the distributions described in Section 3.1 are propagated as an increased uncertainty on all of our flux measurements.
For each filter we fit the location and intensity of the model PSF to the true PSF of HIP 65426 b from the ADI+RDI subtraction using 6 PCA modes.The fitting procedure is executed in an MCMC framework as implemented by emcee (Foreman-Mackey et al. 2013), with 50 walkers for 2200 steps, of which the first 200 steps are discarded as burn-in.spaceKLIP allows for the spatial scale of noise in an image to be fit with a variety of different kernels as implemented in a Gaussian process framework.For NIRCam the noise in our images appears to be correlated at the separation of HIP 65426 b, so, following Wang et al. (2016), we initially adopt a Matérn 3/2 kernel to better capture this spatially correlated noise structure.However, the generated model PSF is slightly mismatched to the observed data and a positive flux residual was present at the PSF core.Future improvements to model PSF generation may alleviate this issue, but in this work we instead assume uncorrelated noise (using the "diagonal kernel" option), which is able to better capture the flux of the companion at the expense of underestimating the obtained error bars.Therefore, more realistic error bars for the NIRCam photometric measurements are determined through a process of companion injection and recovery.For each fil-ter, the best fit model PSF is used to subtract away the companion flux, and is then injected at 20 different position angles spanning 0−105 • and 195−360 • across an equivalent number of duplicate science images.The HIP 65426 b PSF fitting procedure is then repeated on these synthetic PSFs and the standard deviation in the measured flux across all 20 position angles is adopted as the estimated error.In the case of MIRI the observed noise is visually consistent with uncorrelated noise, so we adopt a diagonal kernel here as well.
Results from the PSF fitting procedure are discussed further in Section 5, and images of the data, model, and residuals to each PSF fit are displayed in Appendix C. At the MASK335R nominal IWA of 0.63 , the F444W data achieve a contrast of ∼4×10 −5 , sloping down to ∼1×10 −5 at 1 and then ∼2×10 −6 beyond 3 .In comparison, at the FQPM1140C nominal IWA of 0.36 , the MIRI F1140C data achieve a contrast ∼1×10 −2 , sloping down to ∼2×10 −4 at 1 and then ∼5×10 −5 beyond 3 .For brevity we do not describe the achieved contrast in the other filters, and instead refer the reader to Appendix B.
In the background limited regime beyond 2 , the measured contrasts for NIRCam approximately match the predicted sensitivity of the ADI and ADI+RDI subtractions, and are up to 2 times more sensitive than the prediction for the RDI subtraction.For MIRI in this regime, all three subtraction methods outperform their predicted sensitivity by a factor of 1.5 − 2. In the contrast limited regime below ∼2 we observe further improvements upon the predicted contrast, with the ADI and RDI subtractions demonstrating a factor of up to 5−10 times deeper contrast, and the ADI+RDI subtraction improving by a factor of ∼2.In some instances at the shortest separations below ∼0.6 the contrast does underperform by up to a factor of ∼2 compared to predictions for RDI subtractions.The pri-mary driver for these improvements is likely the improvement in the overall optical and stability performance of JWST compared to expectations.The total throughput is ∼10−20% larger than predictions, driving analogous improvements in signal-to-noise; the overall telescope wavefront error is ∼ 50% smaller than requirements (75/110 nm vs 150/200 nm for NIRCam/MIRI), improving the raw contrast; and the pointing stability of 1 mas is ∼6−7 times smaller than predictions, meaning smaller drifts in alignment behind the coronagraphic mask throughout an observation (Rigby et al. 2022).
The different contrasts as achieved by the ADI, RDI, and ADI+RDI subtractions allow for more concrete recommendations in observing structure for future programs.Most significantly, the improvement between RDI and ADI+RDI subtractions is negligible.Future observers will may be able to achieve their science goals with less telescope time by focusing on purely ADI or RDI subtraction strategies, however, a broader range of observations will be required to rule out ADI+RDI strategies entirely.Within 1−2 the ADI subtraction for MIRI F1550C, and to some extent F1140C, struggles to fully subtract the residual stellar PSF and an RDI based subtraction strategy should be preferred.Although RDI subtractions can provide improvements of up 5−10 times deeper contrast below 1 for some filters (see Figure B), this technique requires a larger amount of observing time due to the cost-intensive nature of dithered reference observations.These time costs can be reduced by performing a smaller number of reference dithers, however, this will in turn reduce the achieved contrast.For example, by comparing subtractions using individual reference dithers from a nine-point dither strategy, Girard et al. (2022) demonstrate a range of 2−10 times worse contrast at 1 compared to a subtraction combining all dithers in a PCA framework similar to that adopted for this work.A precise assessment of the trade space between contrast, observing time, dither pattern, and subtraction strategy for JWST coronagraphy is beyond the scope of this work.Nevertheless, future observers should use simulation tools such as PanCAKE (Carter et al. 2021b) to identify the contrast performance necessary to meet their science goals, and adjust their PSF subtraction strategy accordingly to improve observational efficiency.
Similar NIRCam and MIRI observations have been performed as part of observatory commissioning and are detailed in Girard et al. (2022), Kammerer et al. (2022), andBoccaletti et al. (2022).The observed contrast between this work and these commissioning efforts are in broad agreement, with the regime within ∼2 being contrast limited, and the regime be-yond ∼2 being background limited.For the NIR-Cam F356W observation of HIP 65426 (0.5 Jy) in the background limited regime we reach sensitivities of ∼0.4 µJy (∼20 minutes exposure time), compared to ∼1 µJy in Kammerer et al. (2022) for HD 114174 (G3IV, 2MASS K s =5.202) in the F335M filter (2.6 Jy, ∼55 minute exposure).For MIRI F1140C/F1150C observations of HIP 65426 (0.06/0.03Jy) in the background limited regime we reach sensitivities of ∼2.7/3 µJy (∼30/120 minutes exposure time), compared to ∼9/30 µJy in Boccaletti et al. (2022) for HD 158165(K5V, 2MASS K s =4.704)/HD 163113(K0V, 2MASS K s =2.749) in the F1140C/F1150C filters (0.45/1.45 Jy, ∼75/150 minutes exposure time).Differences in these measured sensitivities are likely driven by the complex interplay of source spectral type, magnitude, and exposure time, and will be better understood in the context of a broader range of future JWST coronagraphic observations.Beyond this work, there is significant potential for the contrast performance for both NIRCam and MIRI to improve.Kammerer et al. (2022) have already demonstrated that the NIRCam bar masks can provide deeper contrasts at shorter separations, if a full 360 • field-ofview is not required.Similarly, in later cycles it may be possible to position a star at the "NARROW" bar mask position in an attempt to reduce the effective inner working angle (IWA) even further.The efficacy of RDI subtractions may improve with the use of a larger PSF library populated by on-sky observations across multiple programs.It may also be possible to perform an effective PSF subtraction using a much larger PSF library composed either entirely or partially of model PSFs.A particular opportunity for improvement will also come from JWST program GO-02627 (Ygouf et al. 2021), which aims to estimate the on-sky instrumental aberrations that drive variations in the observed PSF structure with a model based phase retrieval algorithm.

Mass Sensitivity
Using the obtained contrast curves we also determine the detectable mass limits for our observations following the approach of Carter et al. (2021a).Briefly, we convert our contrast curves to magnitude sensitivity curves using the stellar magnitudes as described in Section 3.1, and then convert these to a mass sensitivity using an interpolation of the evolutionary models of Linder et al. (2019) (BEX) and Phillips et al. (2020) (ATMO) assuming an age of 14±4 Myr (Chauvin et al. 2017) and distance of 107.49±0.40pc (Gaia Collaboration 2020).As in Carter et al. (2021a) we select the chemical equilibrium, non-cloudy models to maintain model consis-tency across mass ranges.Clouds and disequilibrium chemistry likely play a significant role in sculpting the emission of sub-stellar atmospheres, however, an investigation into these effects is beyond the scope of this work.
Following the calculation of mass sensitivity curves we use the Exoplanet Detection Map Calculator (Exo-DMC, Bonavita et al. 2012Bonavita et al. , 2013;;Bonavita 2020) 14 to estimate detection probability maps.In this case, we produce a population of synthetic companions with masses and semi-major axes from 0.1 M Jup to 100 M Jup and 0.1 au to 10,000 au respectively.The inclination is uniformly distributed in sin (i), the eccentricity is distributed using a Gaussian with µ = 0 and σ = 0.3 (excluding negative eccentricities; Hogg et al. 2010), and all other orbital parameters are uniformly distributed.Implicit in this is an assumption that these synthetic companions do not necessarily have a similar inclination to HIP 65426 b.The resulting map takes into account the effects of projection when estimating the detection probability, and the probability for a potential companion to truly lie in the instrumental field of view.
For these particular observations of HIP 65426, we identify the F444W and F1140C filters as the most sensitive to the lowest mass companions for NIRCam and MIRI, respectively, and display their detection probability maps in Figure 6.Of the two, the F444W filter is the most sensitive, reaching a sub-Jupiter mass sensitivity from ∼150−2000 au at a 50% probability, and a minimum sensitivity of ∼0.4 M Jup (∼350 K) ∼300−1500 au at a 50% probability.In contrast, the F1140C is unable to reach sub-Jupiter mass sensitivity, and has a minimum sensitivity of 1.5 M Jup (∼550 K) from ∼150−2000 au.As we detect no sources within the observed field-of-view that have colours consistent with a planetary mass companion, we set equivalent limits on the presence of additional companions in the HIP 65426 system.
Despite the F444W observation being at a shorter wavelength, corresponding to a finer angular resolution, it does not probe significantly closer separations than the F1140C observation.This is driven by the competing influence of the larger IWA for the NIR-Cam MASK335R of ∼0.63 , compared to ∼0.36 for the MIRI FQPM1140C (both assumed at 50% transmission radius).At the wider separations, the F444W observation is sensitive to the lower mass companions than the F1140C observation due to the lower thermal back- ground, which is ∼100 times fainter at 4.5µm compared to 11.3µm (Rigby et al. 2022) 15 .
As discussed in Carter et al. (2021a), at a given distance A stars are generally poor targets for detecting the lowest mass planets in terms of detection sensitivity, whereas M stars are among the most favourable.Given the improved performance of JWST, it is likely that for nearby targets within (or outside of) the M star sample from Carter et al. (2021a) it will be possible to detect Uranus and Neptune mass objects beyond ∼100−200 au, and Saturn mass objects beyond ∼10 au.Initial searches across a small sample of stars for sub-Jupiter mass objects will be performed as part of guaranteed time observations (Schlieder et al. 2017).Furthermore, for the nearest targets we may be able to push these limits even further, with planned observations of α Cen A aiming to be sensitive to 5 R ⊕ companions from 0.5−2.5 au (Beichman et al. 2020(Beichman et al. , 2021)).

HIP 65426 b IN CONTEXT
The known companion HIP 65426 b is clearly detected in all seven of the observational filters using RDI, and all filters except the MIRI F1550C using ADI.These observations represent the first images of an exoplanet to be obtained with JWST, and the first ever direct detection of an exoplanet beyond 5 µm.

Astrometry
The measured astrometry in each of the observed filters is obtained from the ADI+RDI reduction and shown in Table 3.Each of the measured uncertainties for the NIRCam astrometry are propagated with an additional uncertainty of 6.3 mas (0.1 pixels), and the uncertainties of the MIRI astrometry are propagated with an additional uncertainty of 10 mas (0.1 pixels) to account for the assumed precision of the absolute star centering as described in 2.5.All NIRCam and MIRI astrometry are consistent within 1σ, and in combination the NIR-Cam and MIRI astrometry provide a measurement of the separation, ρ=819±6 mas, and the position angle, θ=149.8±0.4 • .To compute these average values, we do not treat the NIRCam filters as independent and instead average both the quantities and their uncertainties.The absolute position of the star on the detector does not change significantly (<1 pixel) between NIRCam filters, and it is feasible that the measured position has a similar systematic offset in each filter (see Section 2.5).As the dominant noise source for the NIRCam alignment is this systematic offset, it is not appropriate to propagate the uncertainties of the NIRCam astrometry in a typical fashion.
We combine these new measurements with the existing astrometry from Chauvin et al. (2017) and Cheetham et al. (2019) to determine updated orbital parameters using the orbitize!package (Blunt et al. 2020).As in Cheetham et al. (2019), we exclude the NaCo epochs due to the inconsistency with the SPHERE epochs.Additionally, we exclude the MIRI The same random sample of 100 orbits is also displayed (grey lines).Bottom: Posterior distributions for the semi-major axis (a), eccentricity (e), inclination (i), argument of periastron (ω), and position angle of nodes (Ω).The 50% quantile from these distributions (green dot-dashed line) are also indicated.In this particular case, these additional JWST observations do not significantly increase the constraints on the measured orbital parameters.
epoch as the observations were taken just two weeks after the NIRCam observations and have larger uncertainties.orbitize! is initialised assuming one companion to the primary (HIP 65426 b), a total system mass of 1.97±0.046M (Chauvin et al. 2004), and a parallax of 9.3031±0.0346mas (Gaia Collaboration 2020).Orbit generation is performed using the Orbits for the Impatient (OFTI) algorithm (Blunt et al. 2020) until 100,000 possible orbits are identified.A random sample of 100 of the possible orbits along with posterior distributions for the entire sample are shown in Figure 7.We are able to constrain the semi-major axis, a = 86 +116 − 31 au, and the inclination, i = 99 +14 − 6 degrees (relative to the equatorial plane).Additionally, as the motion of HIP 65426 b is primarily radial, solutions that place the line of nodes close to its position angle are preferred, and the posi-tion angle of nodes is also constrained to two possible solutions.
The addition of the NIRCam astrometry does not significantly improve the orbital constraints for HIP 65426 b, and all retrieved properties are consistent with those from Cheetham et al. (2019) and Bowler et al. (2020).Although our eccentricity distribution more strongly favours higher eccentricities compared to Cheetham et al. (2019), it remains essentially unconstrained and should not be interpreted as evidence for a highly eccentric orbit for HIP 65426 b.However, if this high eccentricity is real, it would give credence to the scenario proposed in Marleau et al. (2019), where HIP 65426 b initially formed through core accretion before being scattered to a wider separation by an additional companion.The ability of JWST to provide high precision astrometry may improve with improvements to the measurement of the absolute star position in the NIRCam images, which in this case is limited by the precision with which we can locate the star center behind the coronagraphic mask (see Section 2.5).
Separate to JWST, HIP 65426 b has been observed using VLT interferometry as part of the ExoGRAVITY program (Lacour et al. 2020, Program ID: 1104.C-0651), which has routinely demonstrated sub−mas astrometric precision (e.g., Gravity Collaboration et al. 2019;Lacour et al. 2021;Hinkley et al. 2022b) and has an even greater potential to improve upon our reported constraints.

Photometry
As with the astrometry, the measured photometry in all of the observed filters is obtained from the ADI+RDI reduction and shown in Table 3, the subtracted images of all of these filters is shown in Figure 8, and the photometric data points themselves are shown in Figure 9.The measured contrast of the planet relative to the star ranges from 10.132 mag in the F250M filter to 8.029 mag in the F1550C filter.Images for the ADI and RDI subtractions can be found in Appendix A. Additional literature photometric measurements as shown in Figure 9 are provided in Appendix D.

Bolometric Luminosity
With the addition of JWST NIRCam and MIRI photometric observations, the SED of HIP 65426 b is measured across 1 µm to 15 µm.The measurements span the majority of its luminous wavelength range and enable a tight constraint on the bolometric luminosity of the planet.
To calculate the luminosity, a full SED was created by distributing the flux-density from photometric measurements over the effective bandwidth for each filter and using a model atmosphere to extrapolate beyond and interpolate between measured bands.Luminosity is then determined by integrating this semi-empirical SED over wavelength.
Since all of the flux measured in the NIRCam/F410M photometry is accounted for in the F444W measurement, we used only the wider band for our analysis.We also added measurements from the literature at shorter wavelengths, including the SPHERE-IFS YH-band spectrum (Cheetham et al. 2019), and SPHERE-IRDIS H3, K1, and K2-band photometry (Chauvin et al. 2017;Cheetham et al. 2019; also see Appendix D).
To explore the dependence on the details of atmospheric model assumptions, we calculated the bolometric luminosity multiple times, using a broad range of models for interpolation and extrapolation of the SED.
Atmospheric models spanned T eff from 1200 K to 1900 K and log(g) spanning 3.5 to 5.5.These models were drawn from three different grids, including two with different cloud implementations -BT-SETTL (Baraffe et al. 2015), and DRIFT-Phoenix (Witte et al. 2009) -and the cloud-free Sonora-bobcat models (Marley et al. 2021).No matter which model we used to fill in the gaps between the measured portions of the SED, log(L bol /L ) is always between −4.14 and −4.31.Consequently, the luminosity is constrained at the 0.17 dex level and the result is robust across all considered model atmospheres.In totality, the measured luminosity fraction ranges from 61−89% dependant on the adopted model atmosphere.

Estimates of Mass and other Companion
Properties from Hot-Start Evolutionary Models The mass of HIP 65426 b is estimated using a method similar to that described in Dupuy & Liu (2017).We first built an interpolated grid of model luminosities as a function of age and mass, with 10,000 equallyspaced age values spanning from 5 to 30 Myr and 10,000 equally-spaced mass values spanning from 0.3 to 21 M Jup using scipy.interpolate.griddatawith cubic interpolation in python.We adopted an age for HIP 65426 b of 14±4 Myr based on the Lower Centaurus-Crux age given in Chauvin et al. (2017) and a measured bolometric luminosity uniformly distributed between log(L bol /L ) = −4.14 and −4.31 from the previous section.
We then generated 1×10 6 samples of age and mass from a Gaussian distribution in age around 14 Myr, with σ=4 Myr, and a uniform distribution in mass from 0.3 to 21 M Jup .For each age, mass sample, we then look up the corresponding model luminosity from the interpolated grid of model luminosities.For each age, mass sample, we accept the sample if the corresponding model luminosity is within the measured range of uniformlydistributed bolometric luminosities and reject the sample if it lies outside this range.
We implemented this procedure using the hybrid cloud grid from Saumon & Marley (2008).Given that this is a dusty, young red object, we expect the Saumon & Marley (2008) models, which take clouds into account, to provide the most reliable estimates of the properties of these objects among the model choices available.To sample the corresponding effective temperatures, surface gravities, and radii corresponding to our accepted mass values, we built interpolated grids of model effective temperatures, surface gravities, and radii with the same spacing in age and mass as for the interpolated grid of luminosities, then looked up the corresponding prop-    10.The best value of each parameter was taken as the median of the accepted distribution, with error bars given by the 68% confidence interval as calculated from the histogram of each distribution.We find a mass of 7.1±1.2M Jup , radius of 1.44±0.03R Jup (which therefore imply a surface gravity log(g) = 3.93 +0.07 −0.09 ), and effective temperature T eff = 1283 +25 −31 K.

Atmospheric Forward Model Comparison
To explore the atmospheric properties of HIP 65426 b we performed a forward modelling analysis using the tool ForMoSA (Petrus et al. 2020), which compares spectroscopic and/or photometric data with grids of precomputed synthetic spectra.The code is based on the nested sampling algorithm (Skilling 2004), a Bayesian inversion method, that allows a global exploration of the parameter space provided by the grid.In this work, we limited our analysis to the BT-SETTL grid (CIFIST version, Allard et al. 2012;Baraffe et al. 2015) that accounts for convection using mixing-length theory, and works at hydrostatic, radiative-convective, and chemical equilibrium.
For the fit, we used a data set composed of the lowresolution spectra (R λ ∼54) between 1.00 and 1.65 µm provided by: VLT /SPHERE-IFS (Chauvin et al. 2017), VLT /SPHERE-IRDIS H 2 (1.58 µm), H 3 (1.66 µm), K 1 (2.11 µm), and K 2 (2.25 µm) photometry (Cheetham et al. 2019), NaCo L (3.77 µm), NB4.05 (4.06µm), and M (4.76 µm) photometry, and our new JWST NIRCam and MIRI photometry.We first adapted the BT-SETTL synthetic spectra to our data, by reducing their spectral resolution to that of SPHERE-IFS and calculating the synthetic photometric flux at each bandpass using throughputs as obtained from spaceKLIP for the JWST data, and from the SVO filter service for all other data 16 16 http://svo2.cab.inta-csic.es/theory/fps/(Rodrigo et al. 2012;Rodrigo & Solano 2020).We then defined flat priors on the T eff and the log(g) according to the limits of the grid, and applied nested sampling to estimate the posterior distributions of these two parameters.We also add the radius, R, to the list of the parameters explored by the nested sampling.At each iteration, a radius is picked randomly (uniform prior), and a dilution factor C K = (R/d) 2 is calculated and multiplied to the model, where d is the distance to the object (107.49pc).
The best fit models to our data combined with the existing SPHERE and NaCo data are displayed in   Best fit values and 1σ uncertainties are indicated, however, these should be interpreted as the model phase space that fits these data, and not the precision to which these properties can be empirically measured.
By comparing the integrated flux of the best fit model across all wavelengths, and the integrated flux between the shortest and longest wavelength observations, we determine that these observations span ∼97% of HIP 65426 b's luminous range.These results are also in agreement with a similar BT-SETTL model fitting procedure to VLT /SINFONI data of HIP 65426 b (see Petrus et al. 2021).The uncertainties of all parameters are given at 3σ, however, we emphasise that they do not necessarily describe our confidence in the true planetary properties and are better considered as the model phase space that best fits our data.The precision on these measurements is primarily driven by the SPHERE/IFS data, however, we do note some differences that result from the addition of the JWST data.Specifically, when fitting just the SPHERE and NaCo data in isolation, we obtain T eff = 1619 +18 −17 K, log(g) = 3.85 +0.10  −0.09 dex, R = 1.10 +0.09 −0.09R Jup and log(L bol /L ) = −4.12+0.06 −0.06 , again with uncertainties given at 3σ.Therefore, the JWST data improve the precision of the radius and bolometric luminosity by a factor of ∼2, but do not significantly improve the precision on the temperature and surface gravity.
The atmospheric forward model fit yields a luminosity within the bolometric luminosity range of −4.14 to −4.31 that was found from combining SED measurements with models in the regions not covered by the SED.However, the effective temperatures and radii found using the atmospheric forward model are considerably in tension with predictions from evolutionary model fits to the measured bolometric luminosity range (see Section 5.4).In particular, to obtain similar bolometric luminosity values, the forward models favour higher effective temperatures (∼1600 K) and smaller radii (∼1.06 R Jup ) compared to the evolutionary models (∼1300 K, ∼1.44 R Jup ).In fact, the atmospheric forward model used here corresponds to an unphysically small radius for an exoplanet which may still be contracting at these young ages.Thus, we consider the effective temperature and radius predictions from the evolutionary models to be more robust here.
The tension we find between the atmospheric models and evolutionary models is well-documented in the community; atmospheric models have a long history of requiring unphysically small radii and high effective temperatures to fit spectroscopy (see e.g., Marois et al. 2008;Patience et al. 2012).This stems from the different approaches and fundamental parameters underpinning evolutionary and atmospheric modelling techniques.
Atmospheric models produce model spectra as a function of effective temperature, T eff , surface gravity, g, and composition, irrespective of the mass or age of the ob-ject modelled.When fitting atmospheric model spectra to observed spectra, the value of the radius necessary to produce the observed luminosity of the object can be derived if the distance to the object is known.From the radius and the best-fit model surface gravity, a mass value can be derived as well.However, these are not fundamental parameters of the model, but rather extrapolated quantities.
In contrast, evolutionary models couple a similar atmospheric model with a stellar-like interior model, solving for hydrostatic equilibrium, mass conservation, and assuming conservation of energy in shells within the planet (Saumon & Marley 2008).The interior model has fundamental parameters of mass and age; running the model until radiative-convective equilibrium is reached yields a bolometric luminosity and radius.Mass and radius here are then properties of the model, as opposed to the case of the atmospheric models where they are derived quantities dependent in particular on the fit for T eff .
When fitting evolutionary models to observations, we fit directly to the bolometric luminosity.In this work the bolometric luminosity is measured to very high accuracy, as we are integrating over many wavelength bins.Thus, we find a robust fit to other fundamental properties such as mass, surface gravity, and T eff from evolutionary models.However, if we considered the "best-fit" model spectrum corresponding to the parameters (T eff , log(g)) of the best evolutionary model, it would poorly fit the observed spectrum.In contrast, the atmospheric model fits to the spectra involves directly fitting over many wavelength bins.Uncertainties in cloud parameterisation in current models pushes these fits to higher temperatures to explain the observed spectral features; as a result, to balance out to the measured luminosity of the object at its artificially high T eff , the radius derived from the spectral fit is pushed to implausibly low values.As a consequence of the small radius, the mass estimate from the atmospheric models is also unphysically low.
Further JWST observations across a broad diversity of exoplanets and/or brown dwarfs will be able to empirically constrain these model discrepancies as a function of properties such as temperature, mass, and age, and may in turn uncover their precise origin and extent.Until then, the competing benefits and drawbacks between atmospheric and evolutionary mean that they are best considered in tandem as opposed to in isolation.

Future Work
There is a range of additional investigation that can be performed on the data presented here that is worth highlighting, but ultimately falls outside the scope of this work.
Most importantly, it is possible to perform the atmospheric forward model fitting procedure shown here in Section 5.5 with a wide range of state-of-the-art models (e.g., ATMO, Tremblin et al. (2015); Phillips et al. (2020); Exo-REM, Charnay et al. (2018); Sonora, Karalidi et al. ( 2021)), each with their own treatment for the effects of clouds and atmospheric chemistry.Additionally, atmospheric fitting can also be performed using retrieval techniques (e.g., Mollière et al. 2020;Gonzales et al. 2021).Divergences in the measured planetary properties between these models are expected, and a more complete analysis in the context of the relative assumptions of each model will greatly improve our understanding of the true properties of HIP 65426 b.
The precision of the 3−5µm data may be sufficient to provide constraints on the relative atmospheric abundances of CH 4 and CO, which can be impacted by disequilibrium chemistry (Zahnle & Marley 2014;Miles et al. 2020).The F1140C photometry falls slightly under the best fit model, albeit only at ∼1σ, and may be indicative of absorption by small silicate dust grains (Cushing et al. 2006;Suárez & Metchev 2022;Miles et al. 2023).Similarly, the F1550C photometry falls ∼1σ above the best fit model, and may be sensitive to circumplanetary disk emission (Sterzik et al. 2004;Stolker et al. 2020a).

CONCLUSION
In this work we present the first ever scientific observations using the JWST high-contrast imaging modes of both NIRCam from 2−5 µm, and MIRI from 11−16 µm.The known exoplanet companion, HIP 65426 b, is clearly detected in all seven observational filters, representing the first ever direct detection of an exoplanet beyond 5 µm.These observations provide a variety of insights into: a) the performance and best practices of JWST high-contrast imaging, and b) the properties of the HIP 65426 b system, which we summarise below: • Contrast: JWST is exceeding its anticipated contrast performance for both NIRCam and MIRI coronagraphy by up to a factor of 10 in the contrast limited regime (see Section 4.1).For the contrasts achieved, we are sensitive to sub-Jupiter companions with masses as small as 0.3M Jup beyond separations of ∼100 au.Furthermore, for more optimal targets such as young, nearby M stars it is highly likely that both NIRCam and MIRI will be sensitive to sub-Saturn mass objects beyond ∼10 au (Carter et al. 2021a).
• Subtraction Strategy: For these data at small separations <2 , the best contrast is obtained using a small-grid dither RDI subtraction strategy for both NIRCam and MIRI.Additionally, an ADI+RDI subtraction does not significantly improve the measured contrast compared to RDI.
For the MIRI F1550C observations in particular, we were unable to recover HIP 65426 b using ADI alone.At wider separations however, the observational efficiency of ADI may make it preferable to RDI.These conclusions may aid future observers in selecting their PSF subtraction strategy, although we emphasise that a clearer understanding of whether they apply under all circumstances will require the analysis of a broader range of JWST coronagraphic observations.
• Photometry: These photometric observations of HIP 65426 b provide exquisite sensitivity at a precision of ∼7% for NIRCam and ∼16% for MIRI.Furthermore, prior to propagation of the uncertainty in the stellar flux (∼5%), and the current absolute flux calibration accuracy (5/15% for NIR-Cam/MIRI), the uncertainty in the measured relative flux is even smaller at ∼2% for both NIRCam and MIRI.These measurements are a significant step forward from ground-based observations from 3−5 µm, which have comparative uncertainties of ∼13−32% for HIP 65426 b, and are restricted to particular wavelength regions due to telluric contamination.With this improved precision we will be able to constrain directly imaged exoplanet atmospheres in much greater detail, in addition to more complex effects such as variability, disequilibrium chemistry, and the emission of circumplanetary material.
• Atmospheric Model Fitting: Using a BT-SETTL atmospheric forward model we are able to fit all data, in addition to the majority of ground-based observations to within 2σ.
This agreement provides precise constraints on the T eff = 1624 +16 −15 K, log(g) = 3.88 +0.08 −0.08 dex, R = 1.06 +0.05 −0.05 R Jup , and log(L bol /L )=−4.15+0.03 −0.03 .Compared to a fit excluding the JWST data, this corresponds to a factor of ∼2 improvement in the precision of the radius and bolometric luminosity.Despite the excellent model agreement, both the temperature and unphysically small radius are in disagreement with the values obtained from the evolutionary models, further emphasising a long standing tension for this class of objects.
• Empirical Bolometric Luminosity: As JWST offers a uniquely broad spectral coverage in comparison to ground-based instruments, we are able to obtain a very precise measurement of the bolometric luminosity of HIP 65426 b that is constrained between a log(L bol /L )=−4.14 to −4.31, irrespective of the model atmosphere adopted for the wavelengths not covered by observations.In combination with evolutionary models, this provides tight constraints on the properties of HIP 65426 b with M =7.2±1.1 M Jup , T eff =1283 +25 −31 K, and R=1.44±0.03R Jup .Given the achieved sensitivity, similar JWST observations will facilitate this analysis for a broader range of PMCs than ever before and provide comparable constraints on their bolometric luminosities, and therefore mass.These measurements will in turn be valuable for investigating discrepancies between atmospheric and evolutionary models of exoplanets.

Figure 1 .
Figure 1.Colour-magnitude diagram showing the position of HIP 65426 b (Chauvin et al. 2017) relative to both the population of brown dwarf objects (circles) and a selection of directly imaged planetary mass companions (PMCs, hexagons), as obtained from Best et al. (2020).

Figure 2 .
Figure 2. Left: A single integration in the Stage 2 (*calints.fits)file for MIRI coronagraphy of HIP 68245 in the F1550C filter.Right: As on the left, except following subtraction of a median background frame of an "empty" region of the sky.Both images are identically scaled.Before subtraction the residual stellar flux is completely obscured by the stray light "glow stick" (Boccaletti et al. 2022), but can be easily recovered.

Figure 3 .
Figure3.Unsubtracted and KLIP subtracted image stamps for the NIRCam F444W (top row) and MIRI F1140C (bottom row) filters.The leftmost column displays the median unsubtracted image for a single science roll, and all other columns display the KLIP subtracted images for ADI, RDI, and ADI+RDI subtraction methods using the maximum number of KLIP PCA modes.All images are oriented as shown by the directional arrow in unsubtracted image column, the position of the star (white star) is also marked.Additionally, the intensity of all images for a given filter are identically scaled.The exoplanet, HIP 65426 b, can be easily identified at a position angle of ∼150 • in the subtracted images.We note that the distinct "hamburger" shaped central core and six-lobed structure of the companion PSF in the NIRCam images is an expected feature that is related to the Lyot stop design, and not indicative of discrete astrophysical sources.

Figure 4 .
Figure 4.The best fit stellar model (purple) to existing: Hipparcos/Tycho-2, Gaia DR2, 2MASS, ALLWISE, AKARI IRC, and Spitzer MIPS data for HIP 65426 (black circles).The stellar (blue) and disk (red) components of the model are also shown, along with the equivalent model fluxes in each photometric band (solid circles).Error bars are smaller than the circle diameter, except for the 70 µm MIPS point, which is instead marked as an upper limit (black triangle).

Figure 5 .
Figure 5. Contrast curves for observations in the F444W and F1140C filters using both an ADI (dotted lines), RDI (dashed lines) and ADI+RDI (solid lines) subtraction using 20 PCA modes.Both the measured contrast of the true on-sky observations (black lines) and predicted contrasts as generated from PanCAKE (Carter et al. 2021b) (light blue lines) are displayed.The effective inner working angles (IWA) corresponding to the separation at 50% transmission (red dashed lines), and the companion HIP 65426 b (orange hexagon), are also marked.Contrast curves for all other filters are displayed in Figure 14.PSF.As the spatial intensity of the model PSF depends on an assumed spectral energy distribution (SED), we use an existing best fit BT-SETTL model for HIP 65426 b from Cheetham et al. (2019).We note that as the initial model PSF is normalised to unit intensity, the primary purpose of selecting a model SED is not to accurately predict the absolute flux of HIP 65426 b in each bandpass, but instead to capture the relative variation in flux as a function of wavelength across each bandpass (which is largest for the broader NIRCam wide filters).Finally, the impact of the coronagraphic throughput on the received flux is applied by multiplying the normalised PSF by a scale factor equal to the relative difference between the integrated flux of the model PSF prior to normalisation and a matching webbpsf ext model PSF which excludes the coronagraphic elements.The input model PSF is converted to physical units by scaling it to the flux of HIP 65426 as estimated by the stellar model described in Section 3.1.Therefore, all derived photometry is anchored relative to our assumption of the stellar flux.Furthermore, any comparisons between the intensity of this PSF and the observed PSF is subject to an implicit assumption that the absolute flux calibration of JWST is perfect.In reality, the absolute flux calibration accuracy requirement for NIRCam and MIRI coronagraphy is 5% and 15%, respectively 13 .
Achieved ContrastThese observations provide a first look at the achieved contrast for JWST 's high-contrast imaging modes following the completion of observatory commissioning and demonstrate that, in addition to the existing work ofGirard et al. (2022), Kammerer et al. (2022), and Boccaletti et al. (2022), both NIRCam and MIRI coronagraphic imaging are exceeding their predicted contrast performance.Examples of this improvement in comparison to pre-launch contrast estimates as obtained by PanCAKE (Carter et al. 2021b) are demonstrated in Figure 5, and contrast curves for all seven filters are displayed in Appendix B.

Figure 6 .
Figure 6.Detection probability maps as generated by Exo-DMC for the most sensitive NIRCam (F444W, top) and MIRI (F1140C, bottom) filters.Solid black contours signify the 10%, 50%, 80%, and 95% detection thresholds.The true location of HIP 65426 b is marked (hexagon) with errors taken from the astrometric fitting in Section 5.1, and evolutionary model analysis in Section 5.4.

Figure 7 .
Figure 7. Orbital fitting of both JWST NIRCam (this work) and SPHERE (Chauvin et al. 2017; Cheetham et al. 2019) astrometric measurements of HIP 65426 b using orbitize!(Blunt et al. 2017).Top Left: A random sample of 100 orbit models from the retrieved posterior (purple curves).The position of the planet (white circle) and the star (white star) are marked, and the epoch at a given position in an orbit is indicated by the colour bar.Top Right: Separation and position angle versus epoch for both the SPHERE (squares) and JWST NIRCam astrometry (circles).The same random sample of 100 orbits is also displayed (grey lines).Bottom: Posterior distributions for the semi-major axis (a), eccentricity (e), inclination (i), argument of periastron (ω), and position angle of nodes (Ω).The 50% quantile from these distributions (green dot-dashed line) are also indicated.In this particular case, these additional JWST observations do not significantly increase the constraints on the measured orbital parameters.

Figure 8 .Filter
Figure 8. Images of the exoplanetary companion, HIP 65426 b, in all seven NIRCam and MIRI filters used in our observations.Each image is produced following an ADI+RDI KLIP subtraction of the residual stellar PSF.The measured position of the star is marked (white stars), and the orientation and pixel scales of all images are marked in the top left panel.

Figure 9 .
Figure 9.All existing spectroscopic and photometric observations of HIP 65426 b as obtained from SPHERE/IFS (triangles), SPHERE/IRDIS (squares), NaCo (diamonds), and JWST (circles).Top: Data are plotted alongside the 1, 2, and 3σ confidence intervals obtained from fitting to a collection of BT-SETTL atmospheric forward models (blue shaded regions), and the model values in the photometric bandpasses (small blue circles).At 3σ, the best fit models occupy parameter ranges of T eff = 1624 +16 −15 K, log(g) = 3.88 +0.08 −0.08 dex, and R = 1.06 +0.05 −0.05 RJup.Also plotted are the normalised filter throughput profiles for all photometric observations, with the NaCo throughputs scaled by a factor of 2 to improve clarity.Bottom: Residuals of each data point relative to the best fit model in addition to 1, 2, and 3σ regions (grey shading).

Figure 10 .
Figure 10.Histograms of the final sets of accepted model properties for the hybrid cloud grid from Saumon & Marley (2008).The median value for each property is shown as a solid black line, with the 68% confidence region falling between the two dotted black lines.

Figure 11 .
Figure 11.Posterior distributions for the BT−Settl atmospheric model fitting to both JWST and V LT /SPHERE observations of HIP 65426 b.Best fit values and 1σ uncertainties are indicated, however, these should be interpreted as the model phase space that fits these data, and not the precision to which these properties can be empirically measured.
demonstrate that JWST provides an transformative opportunity to study exoplanets through high-contrast imaging.Beyond this work, we also highlight existing and future publications from our program of: 3−16 µm NIR-Cam and MIRI coronagraphy of the circumstellar disk HD 141569 A (Millar-Blanchaer et al. in preparation, Choquet et al. in preparation); NIRSpec and MIRI spectroscopy from 1−28 µm of the PMC VHS J1256 b (Miles et al. 2023) and NIRISS AMI observations of HIP 65426 at 3.8 µm (Ray et al. in preparation, Sallum et al. in preparation).
grateful for the countless hours that thousands of people have devoted to the design, construction, and commissioning of JWST.ALC acknowledges the significant harm caused to members of the LGBTQIA+ community in the Department of State and NASA, while under the leadership of James Webb as Under Secretary of State and NASA Administrator, respectively.This project was supported by a grant from STScI (JWST -ERS-01386) under NASA contract NAS5-03127.This work benefited from the 2022 Exoplanet Summer Program in the Other Worlds Laboratory (OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Foundation.ALC and this work have greatly benefited from ExoExplorers, which is sponsored by the Exoplanets Program Analysis Group (ExoPAG) and NASA's Exoplanet Exploration Program Office (ExEP).This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTH-EAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon/CRAL, and is supported by a grant from Labex OSUG@2020 (Investissements d'avenir -ANR10 LABX56).SP acknowledges the support of ANID, -Millennium Science Initiative Program -NCN19 171.MBo acknowledges support in France from the French National Research Agency (ANR) through project grant ANR-20-CE31-0012.This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (COBREX; grant agreement n • 885593; EPIC, grant agreement n • 819155)).All the JWST presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute and can be accessed via: 10.17909/2bdf-3p61.This research has also made use of the NASA Astrophysics Data System Bibliographic services and the python (Van Rossum & Drake 2009) modules listed below.Software: NumPy (Harris et al. 2020), matplotlib (Hunter 2007), Astropy (Astropy Collaboration et al. 2013, 2018), SciPy (Virtanen et al. 2020), scikit-image (Van der Walt et al. 2014), synphot (STScI Development Team 2018), emcee (Foreman-Mackey et al. 2013), MultiNest

Figure 13 .
Figure13.As in Figure3, except for the MIRI F1140C and F1550C filters.Here we show subtractions using the maximum number of PCA modes for ADI, RDI, and ADI+RDI, respectively.To aid visual clarity, the subtracted F1550C images are shown with a peak image intensity five times smaller than the unsubtracted image.
Figure 15.The data (left column), model PSF (middle column), and residuals (right column) for the spaceKLIP PSF fitting of the NIRCam observations of HIP 65426 b.Pixel counts are in MJy/sr and are indicated by the colour bar on the right hand side, images are oriented with North upwards and North through East counterclockwise, and the image size is 30×30 pixels (1.9×1.9 ).

Table 1 .
Target and reference exposure settings.Background observations were also performed for the MIRI F1140C and F1550C filters with parameters identical to two exposures of a single roll or dither of the target and reference observations, respectively (see Section 2.1).Filter mean wavelengths (λmean) and bandwidths (W eff ) are taken from spaceKLIP (see Section 2.2).See https://jwst-docs.stsci.edu/understanding-exposure-timesfor further detail on JWST exposure settings.

Table 2 .
Filter dependent PSF offsets relative to the center of the NIRCam MASK335R as defined by the F335M filter.

Table 3 .
JWST astrometry and photometry of HIP 65426 b, m * corresponds to the stellar magnitude in each filter, and ∆corr corresponds to the relative magnitude following the propagation in uncertainties of a 5% or 15% absolute flux calibration accuracy for NIRCam and MIRI, respectively.The position angle (θ) is provided from North through East, and all apparent magnitudes are relative to Vega.erty in the appropriate table for each accepted age, mass sample.Histograms of the final set of accepted masses, effective temperatures, surface gravities, and radii for each model are shown in Fig.

Table 4 .
Additional photometric measurements of HIP 65426 b considered in this work.