JWST-TST High Contrast: JWST/NIRCam Observations of the Young Giant Planet β Pic b

We present the first JWST/NIRCam observations of the directly imaged gas giant exoplanet β Pic b. Observations in six filters using NIRCam's round coronagraphic masks provide a high-signal-to-noise-ratio detection of β Pic b and the archetypal debris disk around β Pic over a wavelength range of ∼1.7–5 μm. This paper focuses on the detection of β Pic b and other potential point sources in the NIRCam data, following a paper by Rebollido et al. that presented the NIRCam and MIRI view of the debris disk around β Pic. We develop and validate approaches to obtaining accurate photometry of planets in the presence of bright, complex circumstellar backgrounds. By simultaneously fitting the planet’s point-spread function and a geometric model for the disk, we obtain planet photometry that is in good agreement with previous measurements from the ground. The NIRCam data support the cloudy nature of β Pic b’s atmosphere and the discrepancy between its mass as inferred from evolutionary models and the dynamical mass reported in the literature. We further identify five additional localized sources in the data, but all of them are found to be background stars or galaxies based on their color or spatial extent. We can rule out additional planets in the disk midplane above 1 M Jup outward of 2″ (∼40 au) and away from the disk midplane above 0.05 M Jup outward of 4″ (∼80 au). The inner giant planet β Pic c remains undetected behind the coronagraphic masks of NIRCam in our observations.


INTRODUCTION
The archetypal member of the nearby β Pic moving group (Zuckerman et al. 2001), the A-type star β Pic itself, was the first star around which a circumstellar disk was directly imaged (Smith & Terrile 1984), and one of the first to have a planet directly imaged orbiting it (Lagrange et al. 2009).Already in the discovery paper by Smith & Terrile (1984), the disk was found to be extended to more than 400 au from the star and to be viewed from the Earth in an edge-on configuration.The shape of the disk was suspected to be influenced by planet formation, a hypothesis that was later supported by a variety of observational and theoretical studies (Burrows et al. 1995;Kalas & Jewitt 1995;Mouillet et al. 1997;Heap et al. 2000;Augereau et al. 2001;Wahhaj et al. 2003;Golimowski et al. 2006).In 2009, a gas giant exoplanet was imaged in the inner gap of the debris disk.Discovered by Lagrange et al. (2009) using VLT/NACO observations in the L'-band (∼ 3.8 µm), β Pic b was found to orbit its young (∼ 18.5 Myr, Miret-Roig et al. 2020) host star in a highly inclined ∼ 10 au orbit, coplanar with the debris disk.At this orbital separation, which is similar to that of Saturn, β Pic b remains one of the closest-in directly imaged exoplanets to date (e.g., Currie et al. 2022) and represents a prime testbed for studying the formation and early evolution of gas giant planets on Solar System-like scales.
Thanks to its proximity (distance ∼ 19.44 pc, van Leeuwen 2007) and apparent brightness, β Pic b quickly became one of the most extensively studied exoplanets.Observations from ∼ 1-5 µm probing the thermal emission of the planet are discussed in numerous publications (Lagrange et al. 2009;Quanz et al. 2010;Bonnefoy et al. 2011;Currie et al. 2011;Bonnefoy et al. 2013;Currie et al. 2013;Absil et al. 2013;Males et al. 2014;Baudino et al. 2015;Morzinski et al. 2015;Stolker et al. 2020), and the first low-resolution J-band spectroscopy (obtained with Gemini/GPI) was presented in Bonnefoy et al. (2014b), later complemented by GPI H-and Kband spectroscopy (Morzinski et al. 2015;Chilcote et al. 2017).Based on the aforementioned studies, β Pic b was found to have an effective temperature of ∼ 1700 K, a radius of ∼ 1.4-1.5 R Jup , and a surface gravity of log g ∼ 4.2 dex estimated from the bolometric luminosity of the planet and the age of the system using warmor hot-start evolutionary models with moderate to high formation entropies ≳ 9.75 k B /baryon (e.g., Spiegel & Burrows 2012).Its spectral type lies somewhere between L0-T0.The red L-M color of β Pic b further suggests a cloudy atmosphere (Currie et al. 2013;Bonnefoy et al. 2013;Stolker et al. 2020) with possibly even thicker clouds or stronger water absorption than predicted by common cloudy atmosphere models such as AMES-Dusty (Baraffe et al. 2003) or DRIFT-PHOENIX (Helling et al. 2008).
The age of β Pic and the β Pic moving group has been the subject of numerous publications.Couture et al. (2023) compile an age range of ∼ 11-26 Myr from the literature, with their own dynamical traceback age estimate of 20.4 ± 2.5 Myr.A more recent study from Lee & Song (2024) also finds a kinematic traceback age (16.3 +3.4  −2.1 Myr) that is younger than the age determined from lithium depletion boundary (∼ 24 Myr).As noted in Couture et al. (2023), the dynamical traceback age is the time since when the members of a group became gravitationally unbound, which might be systematically different from the time of their formation if young stars remain bound by interstellar gas and dust for a few million years after their formation.However, studying potential systematics in different stellar age dating methods shall not be the subject of this work.Here, we as-sume the dynamical traceback age of 18.5 Myr inferred by Miret-Roig et al. (2020), which lies somewhere in between the estimates from Lee & Song (2024), Crundall et al. (2019), and Couture et al. (2023), but we note that this age is slightly younger than the overall reported age average.We recognize that assuming a slightly younger age might yield a slightly smaller planet mass and a slightly higher expected planet luminosity based on evolutionary models which will be discussed in Section 4.
Using R ∼ 5000 VLT/SINFONI integral field spectroscopy and cross-correlation techniques, Hoeijmakers et al. (2018) directly detected water and CO in the atmosphere of β Pic b, albeit without the possibility to derive abundance measurements.Two years later and taking advantage of the spatial resolution of long-baseline interferometry to separate the light of β Pic b from contamination by its host star and the debris disk, GRAVITY Collaboration et al. (2020) obtained high signal-to-noise ratio (SNR) R ∼ 500 K-band spectroscopy of the young gas giant planet, clearly showing CO absorption features at ∼ 2.3 µm.With the help of atmospheric forward models and spectral retrievals, they were able to derive a subsolar C/O ratio of ∼ 0.43 for the atmosphere of β Pic b, suggesting strong volatile enrichment during its formation process if the host star β Pic is assumed to be of similar composition as the Sun.We note that measuring the C/O abundance ratio for β Pic is difficult, but it can be done for other stars in the β Pic moving group which can serve as a proxy for the chemical composition of β Pic (e.g., Reggiani et al. 2024).While the C/O abundance ratio of the host star β Pic can hence be expected to be ∼solar or supersolar, the subsolar C/O ratio of β Pic b has recently been challenged by Kiefer et al. (2024) who found that with a more careful treatment of telluric and stellar features, the SINFONI spectrum of β Pic b also supports a ∼solar C/O ratio for the planet.Landman et al. (2024), however, also detect water and CO in β Pic b's atmosphere and confirm its subsolar C/O ratio using high-resolution (R ∼ 100 ′ 000) spectroscopy with the upgraded VLT/CRIRES+ instrument, updating the result that was earlier obtained by Snellen et al. (2014) with the original CRIRES instrument.
More recently, a second gas giant planet was detected in the β Pic system based on extensive radial velocity monitoring of the host star with the HARPS spectrograph at ESO's 3.6 m Telescope (Lagrange et al. 2019).Using interferometric observations with the VLTI/GRAVITY instrument, Nowak et al. (2020) were able to directly confirm this inner planet β Pic c at an orbital separation of only ∼ 3 au.
Combined together, the precise astrometry of the two planets measured with GRAVITY, the long-term radial velocity monitoring with HARPS, and the Hipparcos-Gaia proper motion anomaly of the system (Brandt 2021;Kervella et al. 2022) yield dynamical mass constraints of ∼ 9 ± 2 M Jup and ∼ 8 ± 1 M Jup for β Pic b and β Pic c, respectively (Grandjean et al. 2019;Nielsen et al. 2020;Nowak et al. 2020;Brandt et al. 2021).Lacour et al. (2021) were even able to measure the dynamical mass of β Pic c from the astrometric deflection that it causes on the orbit of β Pic b, similar to the discovery of Neptune from the orbital motion of Uranus by Le Verrier in 1846.The availability of model-independent dynamical masses makes the β Pic planets benchmark objects for probing gas giant planet formation and evolution.

Programmatic context
This paper is part of a series by the JWST Telescope Scientist Team (JWST-TST) 1 , which uses Guaranteed Time Observer (GTO) time awarded by NASA in 2003 (PI: M. Mountain) for studies in three different subject areas: (a) Transiting Exoplanet Spectroscopy (lead: N. Lewis); (b) Exoplanet and Debris Disk High-Contrast Imaging (lead: M. Perrin); and (c) Local Group Proper Motion Science (lead: R. van der Marel).Here, we present new results in the second area; previously reported results across these areas include Grant et al. (2023) and Gressier et al. (in prep.);Ruffio et al. (2023) and Rebollido et al. (2024);and Libralato et al. (2023), respectively.A common theme of these investigations is the desire to pursue and demonstrate science for the astronomical community at the limits of what is made possible by the exquisite optics and stability of JWST.
In this paper, we continue to present the first coronagraphic observations of the β Pic system with the James Webb Space Telescope's (JWST, Gardner et al. 2006Gardner et al. , 2023) ) Near-InfraRed Camera (NIRCam, Rieke et al. 2003Rieke et al. , 2023)), spanning a wavelength range from ∼ 1.7-5 µm.We previously reported in Rebollido et al. (2024) on the edge-on debris disk around β Pic as seen in these same NIRCam data, along with coronagraphic observations at 15.5 and 23 µm obtained with the Mid-InfraRed Instrument (MIRI, Rieke et al. 2015;Wright et al. 2015Wright et al. , 2023)), which led to the discovery of a dramatic curved "tail" of dust that appears to extend sharply away from the plane of the disk, and may be unbound debris from a recent major collision of large planetesimals.
1 https://www.stsci.edu/∼ marel/jwsttelsciteam.html In this work, we turn our attention to the planets within the β Pic system.The outer giant planet β Pic b is detected in all six observed NIRCam filters close to its greatest orbital elongation.The inner giant planet β Pic c remains undetected behind the coronagraphic masks of NIRCam.Both known planets β Pic b and β Pic c remain undetected in the MIRI coronagraphy data due to the tremendously bright debris disk in the mid-infrared, so that we choose to not further discuss those MIRI data here, except for one brief usage in Section ?? for evaluating the nature of additional candidate companions detected in the NIRCam data.However, there are also MIRI Medium Resolution Spectroscopy (MRS, Wells et al. 2015;Argyriou et al. 2023) observations of the β Pic system (program 1294, PI: C. Chen) which detect the outer planet β Pic b at longer wavelengths (∼ 5-7 µm, Worthen et al. 2024); we make use of their spectrum in joint analyses below.
The remainder of this paper is organized as follows.Section 2 presents the observations and data reduction, including two independent approaches to PSF subtraction.Section 3 presents the measured photometry and astrometry of β Pic b, and detection limits on lowermass planets in the outer system.In Section 4, we discuss the implications for β Pic b's atmospheric and bulk properties and its evolutionary status using the NIRCam measurements, prior literature values, and exoplanetary model grids.We summarize our findings and conclude in Section 5.

OBSERVATIONS AND DATA REDUCTION
The JWST data presented here were taken as part of GTO program 14112 ("Coronagraphy of the Debris Disk Archetype Beta Pictoris", PI: C. Stark) on 18 March 20233 .The program was mainly designed to characterize the disk around β Pic using NIRCam filters sensitive to the presence of water, CO ices, and organic tholins and MIRI filters probing the warm inner asteroid belt analogue and the cooler outer main disk.However, given the timing of JWST 's launch, the well-known gas giant planet β Pic b happened to be close to its greatest orbital elongation during the observations, enabling us to confidently detect and characterize its atmosphere from ∼ 1.7-5 µm.

Observational structure
The NIRCam observations targeted β Pic in six different filters from ∼ 1.7-5 µm, two of which were observed using NIRCam's short wavelength (SW) channel with a pixel scale of ∼ 31 mas (F182M, F210M) and four of which were observed using its long wavelength (LW) channel with a pixel scale of ∼ 63 mas (F250M, F300M, F335M, F444W).The F444W wide-band filter was included because it is the most sensitive to additional, yetundiscovered planets in the system (e.g., Carter et al. 2023).The SW observations were taken using NIR-Cam's 210R round coronagraphic mask; these were the first-ever science observations taken with SW channel coronagraphy.The LW observations were taken using the 335R mask.These masks were chosen because they are the smallest available round coronagraphic masks providing the best performance at small angular separations while retaining access to the entire 360 deg fieldof-view.We note that the observations were taken sequentially in the two channels, as parallel operation of the SW and LW channels was not yet enabled at the time.
The observations were conducted using the standard high-contrast imaging strategy of two rolls on the science target plus a PSF reference star.The science target β Pic was observed at two roll angles offset by a position angle (PA) of ∼ 10 deg from one another.This enables point-spread function (PSF) subtraction using angular differential imaging (ADI) techniques (Liu 2004;Marois et al. 2006), which are frequently used to detect faint companions (e.g., Bowler 2016;Currie et al. 2022).In addition, a PSF reference star (α Pic) was observed in sequence with the science target using the 5-POINT-BOX4 small grid dither (SGD) pattern.This enables reference differential imaging (RDI), including classical PSF subtraction methods which are typically less prone to self-or oversubtraction of extended circumstellar structure such as the debris disk around β Pic (Lafrenière et al. 2007;Soummer et al. 2012).The SGD can be used to enhance the quality of the PSF subtraction (i.e., to suppress residual stellar speckles in the PSF-subtracted images) by sampling changes in the shape of the PSF as a function of the relative misalignment between the star and the coronagraphic mask (Soummer et al. 2014).This is necessary because the target acquisition (TA) procedure for NIRCam coronagraphic imaging can have uncertainties of typically up to ∼ 20 mas as measured in flight (Girard et al. 2022;Rigby et al. 2023).Table 1 shows a summary of the NIRCam coronagraphic observations.

Data reduction
All data were reduced with the spaceKLIP5 community pipeline (Kammerer et al. 2022b;Carter et al. 2023).This pipeline was developed during JWST commissioning and Early Release Science (ERS) and has since been updated and improved steadily.The data reduction process is briefly summarized in the following, with the individual steps described in more detail in the subsequent subsections.
1. Reduce the raw ("uncal") files with the jwst6 stage 1 and 2 pipelines (Bushouse et al. 2023) first to "rateints" and then to flux-calibrated "calints" files, with several custom adaptions specific to NIRCam coronagraphy data that are implemented in spaceKLIP.
2. Use the spaceKLIP image processing library to clean bad pixels and align the PSFs in preparation for PSF subtractions.
3. Perform PSF subtractions.Several complementary methods were employed, as described below.In particular, we used the pyKLIP7 community pipeline (Wang et al. 2015) invoked through spaceKLIP to perform PSF subtractions using projections on Karhunen-Loève eigenimages (KLIP, Soummer et al. 2012).We also applied Model-Constrained Reference Differential Imaging (MCRDI, Lawson et al. 2022) as an alternative approach to help separate disk and planetary fluxes.
4. For the KLIP subtractions, we use the pyKLIP forward-modeling routines (Pueyo 2016) with NIRCam coronagraphy PSF models generated through spaceKLIP to extract companion astrometry and photometry and compute contrast curves.

Calibration of individual images
In the first stage, we downloaded the "uncal" files from the MAST archive and processed them with the Coron1Pipeline and Coron2Pipeline within spaceKLIP.These two pipelines are customized implementations of the jwst Detector1Pipeline and Image2Pipeline.We used version 1.12.1 of the jwst pipeline and calibration context jwst 1174.pmap for reducing the NIRCam data.As in Carter et al. (2023), we disabled the flagging of pixels that are diagonal neighbors to saturated pixels, skipped the pipeline's dark subtraction step, flagged a four pixels wide border around the edge of the subarrays to be used as "pseudo" reference pixels, and set the jump rejection threshold to 4. In addition, and as already discussed in Rebollido et al. (2024), we also employed a custom 1/f stripe noise correction step that is now implemented in spaceKLIP and yields a mostly cosmetic improvement.The output of this stage are flux-calibrated "calints" files in units of MJy/steradian.This reduction made use of the first inflight photometric calibrations for NIRCam coronagraphy delivered in fall 20238 .That on-sky flux calibration accounts for the reduced throughput of the NIRCam coronagraphic Lyot stops (pupil plane masks), but not for the attenuation of the focal plane masks for sources located at small angular separation from the observed target9 .It also accounts for the wavelength-dependent transmission of the coronagraphic mask (COM) substrate (Krist et al. 2010).The attenuation of the coronagraphic (focal plane) masks is accounted for later in the companion fitting step, when the model PSF of the companion is computed (see Section 3.1).
In the second stage, we prepared the reduced images for processing with pyKLIP.First, we median-subtracted each image to remove sky background flux, and then cleaned bad pixels.Due to the extended debris disk visible in all NIRCam images, the automatic bad pixel identification routine in spaceKLIP did not work satisfyingly well and we used a custom bad pixel map in addition to the DO NOT USE pixels flagged by the jwst pipeline.This custom map was obtained by inspecting the cleaned images by eye and manually flagging bad pixels that the jwst pipeline had missed.While we found no additional bad pixels in the SW images, we identified and manually flagged 28 bad pixels in the LW images.The routines that were used to clean the flagged bad pixels are the same as the ones in Carter et al. (2023).
We include additional steps to mitigate the flux from the debris disk and to deal with spatial undersampling.The shortest filters in each channel are spatially undersampled, which can lead to numerical artifacts when shifting or rotating images.Before recentering and aligning the images, we first applied a Gaussian highpass filter to subtract spatially extended and smooth flux from the debris disk.We explored a variety of filter sizes ranging from standard deviations of 1-9 pixels, as well as no high-pass filtering at all.To address the undersampling, we then blurred all images with a Gaussian kernel with a full width half maximum (FWHM) of Here, λ min is the minimum wavelength of the filter bandpass, D is the effective telescope diameter (5.2 m for NIRCam coronagraphy due to the undersized Lyot stops), s is the detector pixel scale, and the desired FWHM is 2.3 • 1.5 = 3.45 pixels 10 .The blurring helps to avoid "Fourier ripples" when numerically shifting or rotating undersampled images (due to the Gibbs phenomenon, e.g., Gottlieb & Shu 1997).

Image alignment and target acquisition performance
Image alignment to precisely register the science images and reference PSFs is critical for achieving highquality PSF subtractions.For each filter independently, the first science target image of each observation was recentered on the position of the target star.This is a difficult task because the star is located behind the coronagraphic mask so that its center cannot be easily determined.We followed the same procedure that was already used by Greenbaum et al. (2023) to center the star based on its speckle pattern with the help of a simulated PSF computed from a contemporaneous wavefront map of the telescope using WebbPSF ext 11 (Perrin et al. 2014;Girard et al. 2022).Greenbaum et al. (2023) reported a centroiding error of ∼ 7 mas for this procedure.However, given the bright debris disk around β Pic, the centroiding error is likely larger in our case.Once the first image of each observation was recentered, we then used the same image registration routine as in Kammerer et al. (2022b) to align all subsequent science and PSF reference target images to the first one.As discussed in Rigby et al. (2023), the pointing stability of JWST is ∼ 1 mas, and we are able to recover a ∼ 1 mas root-mean-square (RMS) jitter between individual images.We are also able to recover the injected 5-POINT-BOX SGD for the PSF reference target and to measure the relative TA offset between the science and the PSF reference target.
Figure 1 shows this alignment for one of the SW and one of the LW filters, revealing good TA performance with an offset of only ∼ 4 mas for the SW channel, but poor performance with an offset of ∼ 50 mas for the LW channel PSF reference observation.This is atypical, in and the LW (bottom) observations.The dots are color-coded by dither position and show the measured offsets between the reference star and the science target observation (first roll), while the crosses show the commanded and expected offsets in case of a perfect TA.The TA performance is good for the SW channel, but worse than expected for the LW channel (see Section 2.2.2 for details).The plots show the F210M and F335M datasets; all other filters in the same NIRCam detector channels show similar performance.
fact we believe one of the largest-observed TA offsets in the ensemble of NIRCam coronagraphy data yet obtained.We subsequently determined that the reason for this outlier TA performance is that the SIMBAD coor- dinates and proper motion of α Pic are not sufficiently accurate12 , so that after the initial telescope slew and fine guidance sensor acquisition, α Pic was offset by ∼ 0.7 arcsec from the expected position.As a result, its PSF did not entirely fit within the NIRCam LW TA subarray.This resulted in an inaccurate centroid measurement.We note that the same systematic position offset of α Pic was also seen in the NIRCam SW and MIRI TA images, but for those modes the relative size of the PSF with respect to the TA subarray is smaller so that the PSFs did still fit within the TA subarrays and the centroid measurements were accurate.We expect the poor NIRCam LW TA performance to impact the quality of the RDI PSF subtraction at small separations from the target for the LW datasets, specifically to reduce the contrast performance by a factor of ∼ 3-4 at a separation of one arcsecond with respect to what would have been achievable with a good TA performance (see Figure 16 in Girard et al. 2022).
After image registration, we then averaged all aligned integrations in an individual exposure (each dither position is an individual exposure) to speed up the subsequent processing with pyKLIP.This is beneficial for rerunning the reduction several times with different parameters later on in the paper.We note that averaging the images has a negligible impact on the KLIP reduction because the line-of-sight pointing and PSF of JWST are extremely stable.For comparison, we also ran the pyKLIP processing with our baseline parameters once without averaging the images and found that ∼ 99.9% of the signal is contained in the first KL mode of each individual exposure which represents the noise-weighted average of the individual images in that exposure.In other words, there is negligible information loss in averaging together all integrations within an exposure.

PSF subtractions (KLIP)
We tested and compared several methods for PSF subtractions, including reference differential imaging (RDI, Lafrenière et al. 2007), angular differential imaging (ADI, Marois et al. 2006), and model-constrained refer-ence differential imaging (MCRDI, Lawson et al. 2022).The former two were employed using the implementation of the Karhunen-Loève Image Processing (KLIP) algorithm (Soummer et al. 2012) in pyKLIP (Wang et al. 2015) through the spaceKLIP community pipeline and are described in this Section, and MCRDI is described in the following Section.
To remove post-coronagraph residual stellar speckles from the NIRCam images of β Pic, we employed ADI, RDI, and ADI+RDI techniques through pyKLIP.The objective of the KLIP algorithm is to project each science target image onto a covariance-weighted orthogonal basis of eigenimages of the PSF reference library, which consists of the PSF reference images of all five dither positions in the case of RDI, or the science target images from the other telescope roll in the case of ADI.The reference library constructed from the SGD observations of the PSF reference star will hence capture changes in the PSF shape as a function of the TA offset in order to provide a good representation of the science target image which is similarly affected by a TA uncertainty of typically up to ∼ 20 mas (Girard et al. 2022).Then, the projection of the science target image onto the reference library is subtracted from the science target image itself, ideally leaving behind flux from off-axis companions and circumstellar structure.
For β Pic, RDI subtraction alone does not work well due to the extended debris disk contributing significant flux in all NIRCam images, and additionally the poor PSF reference star TA performance in the LW channel.The KLIP routine tends to oversubtract the residual stellar speckles because it tries to minimize the total residual flux, including the flux of the debris disk.ADI alone yields a much better speckle subtraction, but its performance at small angular separations is limited due to the small telescope roll angle of only 10 deg.
Combining ADI and RDI yields the visually best performance.The KLIP-subtracted NIRCam images are shown in Figure 2. We did not split the images into multiple annuli or subsections because we did not find any improvement in fake companion injection and recovery tests from doing that (see Appendix A).We note that these PSF subtractions are optimized for measurements of the planet β Pic b; different choices are optimal for the study of the debris disk, as discussed in Rebollido et al. (2024).

PSF subtractions (MCRDI)
In Model-Constrained RDI (MCRDI), the effects of RDI oversubtraction (e.g., Pueyo 2016) are avoided by explicitly assuming the presence of circumstellar flux while the stellar PSF model is constructed (Lawson et al. 2022).For this purpose, a synthetic model of the circumstellar scene is optimized using standard forward modeling techniques for an initial unconstrained RDI reduction.A final MCRDI reduction is then carried out using the resulting best-fit model (Lawson et al. 2022).
The MCRDI procedure predominantly follows the approach outlined in Rebollido et al. (2024) but uses slightly different input data resulting from the preprocessing details described in Sections 2.2.1 and 2.2.2 (e.g., the larger blurring kernel to mitigate the effects of the Gibbs phenomenon).The circumstellar model prescription effectively follows the one in Rebollido et al. (2024) -assuming a circumstellar scene that is the superposition of a simple ring-like disk and a point source.We point out that this simplified disk model is only used for computing the PSF subtraction coefficients in the RDI framework, the PSF-subtracted scene images do however still comprise the disk in its full complexity.For convolution, we use a more finely spatially sampled grid of synthetic PSFs from WebbPSF (Perrin et al. 2014), now sampling the origin along with 12 logarithmicallyspaced radial positions at each of 8 linearly-spaced azimuthal positions (for 97 spatial samples in total).The finer spatial sampling has no discernible impact on the final MCRDI image, but does slightly change the diskmodel-subtracted residuals within the IWA.Otherwise, the MCRDI procedure is as described in Rebollido et al. (2024).
The MCRDI-reduced NIRCam images are shown in Figure 3, zoomed in to show just the area centered on β Pic b, after the subtraction of the disk model for each filter.Fairly good PSF subtractions are achieved for the two SW filters, but a small blob can be seen in the residuals of the F250M data towards the North-East of the host star position.We note that a similarly shaped blob can also be seen in other of the LW filters if the color stretch is optimized for the residuals instead of the planet β Pic b and the blob's brightness scales with the brightness of the other residual stellar speckles.We hence attribute this blob to the subpar TA performance for the LW observations.

Impact of the disk on the β Pic b photometry
The debris disk around β Pic is clearly detected in all JWST /NIRCam images and thus affects the extracted photometry of the giant planet β Pic b.The impact of this disk on the extracted planet flux in the KLIP subtractions can be complex.Residual disk flux may add on top of the planet flux, leading to an overestimation of the planet flux, or the disk may be oversubtracted and part of the planet flux may be removed, leading to an underestimation of the planet flux.Self-subtraction artifacts from running ADI on extended sources may further introduce biases in the planet flux that differ from these expectations.An important parameter that controls the impact of the disk on the extracted planet flux is the size of the high-pass filter that is used to remove (part of) the disk flux from the NIRCam images (see Section 2.2.1).Fake companion injection and recovery tests (see Appendix A) reveal that the best performance in the two SW filters is achieved using ADI alone (with a high-pass filter size of 7 pixels) whereas the best performance in the four LW filters is achieved by combining ADI and RDI (with a high-pass filter size of 3 pixels), so that we adopt these as the baseline approach for extracting planet photometry in the respective cases.We note that even with these optimized parameters, the retrieved flux in the F250M dataset is still much higher than the injected one in these tests, showing that the coarse detector sampling at 2.5 µm together with the extended debris disk prohibits an accurate companion flux measurement using KLIP techniques.
In the MCRDI approach, rather than having parameters for spatial filtering, instead there are additional free parameters for the disk model which are fit as part of the MCRDI process.In effect, these disk model parameters become nuisance parameters which must be solved for simultaneously, increasing the overall dimensionality of the fitting process, but providing a loosely physical approach for subtracting the disk's light around the location of the planet.As in Rebollido et al. (2024), we note that the disk model fit in MCRDI is intentionally simplified and not intended as a physically-correct model of the β Pic debris disk, but it does suffice to largely remove the disk light while avoiding biases from oversubtraction (Lawson et al. 2022).
Given the substantial impact of β Pic's bright disk on the PSF subtractions, much effort and iteration went into optimizing these methods and their parameters until the two independent PSF subtraction approaches yielded consistent measurements for the photometry of β Pic b.This cross-validation was essential for establishing confidence in the removal of subtraction systematics.Potential sources of systematic errors were taken into account by inflating the error bars on the planet photometry as described in the third paragraph of Section 3.1.

RESULTS
β Pic b is robustly detected in all six observed filters, from which we extract photometry and astrometry measurements.We also use these data to set deep limits on the presence of additional outer planets in the system.As expected, the inner giant planet β Pic c remains un-detected behind the coronagraphic masks of NIRCam.Its angular separation at the time of these observations was only 95 ± 10 mas, resulting in a 7.3 +0.5 −0.3 mag attenuation due to the coronagraphic masks and leading to an effective contrast of ∼ 18.1 mag or ∼ 6e−8 in the K-band which is beyond the capabilities of NIRCam at such small separations.

JWST photometry
We extracted new photometry for the giant planet β Pic b in all six observed filters, spanning a wavelength range of ∼ 1.7-5 µm.This new photometry is especially interesting since L-and M-band observations from the ground are typically affected by a high thermal background from the Earth's atmosphere, leading to large uncertainties in the measured photometry (e.g., Stolker et al. 2020).Our new measurements from JWST are not affected by this issue.Furthermore, they also open up bands such as F182M and F300M which are inaccessible from the ground due to telluric absorption from water vapor in the Earth's atmosphere.
The photometry was obtained by forward-modeling the planet PSF in all of the six observed NIRCam bands and fitting it to the KLIP-subtracted images using Forward Model Matched Filtering (FMMF, Pueyo 2016).The methods used follow the ones in Kammerer et al. (2022b) and Carter et al. (2023), but now using new inflight photometric calibrations for NIRCam coronagraphy instead of the pre-flight ones.Briefly, a model PSF of β Pic b was computed at its expected distance from the coronagraphic mask center with the WebbPSF ext package using the closest in time available JWST wavefront measurement from MAST.This PSF model accounts for the attenuation by the coronagraphic mask.Due to different TA and therefore coronagraphic mask offsets, this needed to be done separately for each of the two telescope rolls.The model PSFs were then fed into the pyKLIP BKA routine (Wang et al. 2016) which forward-models the KLIP-processed PSF to account for companion self-subtraction (due to ADI) and KLIP algorithm throughput losses.Finally, the forwardmodeled PSFs were fitted to the KLIP-subtracted images using Markov Chain Monte-Carlo (MCMC) sampling with emcee (Foreman-Mackey et al. 2013).The planet photometry (and astrometry, see next Section), including uncertainties, were obtained from the converged MCMC chains.An advantage of JWST over ground-based imagers is that the planet photometry can directly be measured in the images since they were fluxcalibrated by the jwst pipeline, except for the coronagraphic mask throughput which is accounted for separately in our PSF forward models.We note that this en-tire procedure is integrated within a single user-friendly function within spaceKLIP.Figure 4 shows our best fit forward model PSFs compared to the data.
The new β Pic b NIRCam photometry is shown in Table 2 and Figure 5.The KLIP photometry was obtained using the baseline extraction parameters defined in Section 2.3.Besides the statistical uncertainties from the MCMC fit of the forward-modeled PSF to the data, we estimate systematic uncertainties of ∼ 2% from the uncertainty of the absolute flux calibration of NIRCam13 , ∼ 1% from numerical inaccuracies in the forward-modeled PSFs from WebbPSF, and ∼ 2% from the uncertainty on the coronagraphic mask throughput (due to the uncertainty on the companion/mask position) that need to be added in quadrature to the statistical uncertainties.The companion flux uncertainty from the uncertainty on the companion/mask position was derived by shifting companion and mask closer together and further away from each other by their respective position uncertainties and evaluating the resulting change in flux using WebbPSF.The uncertainty in the WebbPSF models is motivated by Weisz et al. (2024) who find sub-1% photometric errors for WebbPSF models in regular NIRCam imaging, and we budget 1% for this effect to be slightly conservative for the more complex optical situation of coronagraphy.In addition, and based on the results from the fake companion injection and recovery tests (Appendix A), we added another systematic error of 5% to the KLIP photometry which accounts for the typical offset that we find between the injected and recovered fake companion flux.
Figure 5 shows the KLIP and the MCRDI photometry together with existing photometry and spectra and a DRIFT-PHOENIX model atmosphere of β Pic b from the literature.The model atmosphere was fitted to the literature data only (and not to the new NIRCam data) using the species14 toolkit (Stolker et al. 2020).The purpose of this Figure is to show whether the new NIRCam photometry agrees with the existing literature data and qualitatively assess whether the NIRCam data points at wavelengths inaccessible from the ground are broadly consistent with models used to interpret the ground-based data; model atmosphere fitting including the new NIRCam photometry will be more rigorously conducted in Section 4.1.Figure 5   ing literature data of β Pic b.Both the KLIP and the MCRDI data points fall within 3σ from the prediction of the model atmosphere.
Notably, all NIRCam data points fall below the prediction of the model atmosphere.In the L-and M-band, this makes the new NIRCam photometry consistent with the existing VLT/NACO photometry which also falls slightly below the prediction of the model atmosphere.
In the K-band, however, there is a ∼ 10-20% discrepancy between the new NIRCam photometry and the existing VLTI/GRAVITY spectrum.We note that the blurring, although being accounted for in the forward modeling, might be responsible for blending some of the planet flux with the disk flux and leading to an underestimation of the planet flux.The disk being brighter at shorter wavelengths makes this especially problematic in the K-band.We further note that the fake companion injection and recovery tests (Figure 14) show that the retrieved planet flux declines with increasing high-pass filter size.In these tests, a rather large FWHM was ideal to retrieve the injected flux for the fake companion, but it could well be that for the real planet β Pic b, a smaller FWHM might be ideal so that we are underestimating its flux.This is extremely difficult to quantify though, since the PSF-and disk-subtraction residuals are not exactly the same on the North-East and the South-West side of the debris disk.Despite these uncertainties, it is worth mentioning that our NIRCam photometry is in good agreement with the NACO data point at ∼ 2.1 µm though.Besides issues with the data reduction, astrophysical variability might also play a role for the measured discrepancies.While variability in the planet itself is expected to be below ∼ 5% based on the observed variability of low-mass brown dwarfs and planetary-mass objects (Metchev et al. 2015;Vos et al. 2020Vos et al. , 2022)), the findings from Rebollido et al. (2024) suggest that the debris disk around β Pic is a highly dynamic environment.Moreover, Worthen et al. (in prep.)find significant variability between their MIRI/MRS and archival Spitzer spectra (Lu et al. 2022) of the debris disk supporting this hypothesis.Nevertheless, the impact that this disk variability would have on the flux of the planet β Pic b is expected to be only a fraction of the disk variability itself and hence smaller than the observed discrepancy between GRAVITY and NIRCam, NACO, and GPI (see also Section 4.1), so that systematics in the flux calibration remain the more likely hypothesis.
Due to the complications of the debris disk affecting the planet flux measured in the KLIP-subtracted images and the advantages of using MCRDI in such cases as discussed in Section 2.3, we decided to use the planet flux measurements from the MCRDI reductions as our baseline for the atmospheric characterization of β Pic b in Section 4.

JWST astrometry
From the KLIP-subtracted JWST /NIRCam images of the β Pic system, we also extracted relative astrometry of the well-known outer giant planet β Pic b.This astrometry is shown in Table 3 and Figure 6.It was obtained using KLIP ADI+RDI and a high-pass filter size of 3 pixels for all observed NIRCam bands15 .While the photometry of β Pic b was taken from the MCRDI reduction, we found that the planet astrometry from MCRDI is rather sensitive to the exact choice of disk model parameters (e.g., one-component vs twocomponent disk model).On the other hand, the KLIP astrometry is much more consistent for different highpass filter size choices, so that we decided to use this more robust measurement.
In the NIRCam images, β Pic b is detected close to its maximum orbital elongation on the North-East side of the disk.The astrometric measurements from NIR-Cam are not yet competitive with high-contrast imagers from the ground, mainly due to the difficulties in finding the position of the star behind the coronagraphic mask and the missing distortion correction in spaceKLIP.Together, they lead to systematic uncertainties on the order of ∼ 10 mas in the NIRCam astrometry.We further expect that the numerous residual stellar speckles caused by the poor reference star TA performance in the LW channel negatively affect the precision of the NIRCam astrometry, so that its systematic offset with respect to the GRAVITY orbit shown in Figure 6 is not a surprise.This expectation is supported by the LW channel astrometry generally agreeing worse with the GRAVITY orbit than the SW channel astrometry.Given these systematic errors and the precision of the existing GRAVITY astrometry of only ∼ 0.1-0.3mas (Lacour et al. 2021), we did not attempt any orbital fits for the β Pic system with our new NIRCam measurements.Instead, we only show the orbits of both planets inferred by Lacour et al. (2021) from the GRAVITY data and overplot the new β Pic b NIRCam astrometry to demonstrate that we indeed detected the planet at its expected location.Improving the NIRCam astrometry and making it competitive with ground-based imagers is left for future work.Filter ∆RA [mas]  regime in the wide-band filters (e.g., Carter et al. 2023).

Limits for other companions in the β Pic system
Here, we use the new NIRCam data of β Pic to search for previously undiscovered companions.
By visually inspecting the KLIP-subtracted images, we identify five other sources at separations > 5 arcsec from β Pic in addition to the known outer giant planet β Pic b (see Table 4 and Figure 7).We note again that NIRCam coronagraphy does not achieve the IWA re-quired to detect the inner known giant planet β Pic c at a separation of ∼ 100 mas at the time of our observations.With the exception of one source (BG4), the additional sources are only detected in the NIRCam LW channel, and are visible at a low number of KL modes (≤ 4) in all four LW filters.BG4 is also visible in the MIRI F1550C four quadrant phase mask data from Rebollido et al. (2024).Three of the five sources (BG1, BG2, BG3) are located to the East of the edge-on debris disk and the other two sources (BG4, BG5) are located to the West of the disk.The two sources to the West are clearly resolved and appear diffuse, so that we classify them as background galaxies.Two of the three sources to the East (BG2, BG3) also appear resolved so that we count four background galaxies in total.Only one of the sources to the East (BG1) appears point-like and unresolved.Its approximate position with respect to β Pic is ∆RA ≈ 6 ′′ and ∆Dec ≈ 8 ′′ .Approximate astrometry of all five sources can be found in Table 4.
To determine the nature of the point-like source (BG1), we extract its photometry by fitting forwardmodeled PSFs to the KLIP-subtracted images, similar to how we measured the photometry of β Pic b.   (2019).This was necessary to cover the entire planet mass range from ∼ 0.1-13 M Jup to which JWST is sensitive.We first interpolated the evolutionary models at an age of 18.5 Myr, and then stitched together the ATMO and BEX models by interpolating linearly between them in log(mass) space.The transition between the two models happens between 0.5 and 1 M Jup where both models are defined.We focused on chemical equilibrium models without clouds to maintain consistency among the different companion mass ranges.We further assumed a distance of 19.44 pc (van Leeuwen 2007) for β Pic.The point-like source is far too blue in order to be a young and bound planetary-mass companion.Instead, its colors are more consistent with a background star.We note that we did also inspect archival HST /STIS and WFC3 data (program ID 12551, PI: D. Apai, and 11150, PI: J. Graham, respectively), but did not identify any background star whose position and proper motion would agree with the point-like source seen in the NIR-Cam data.Given the position of the source near the edge-on debris disk and the proper motion of β Pic, it seems possible that the source might have been right behind the disk in the archival HST images.In any case, its colors indicate it is almost certainly a star.
Due to its unparalleled sensitivity, the NIRCam data also allow us to put deep constraints on the presence of additional companions in the β Pic system.First, contrast curves were computed from the KLIP-subtracted images assuming a noise distribution following a student's t-distribution and correcting for small sample statistics following Mawet et al. (2014).To avoid be- Source ing impacted by residual flux from the debris disk when computing the contrast curves, we here consider the high-pass filtered data (with a high-pass filter size of 3 pixels).Nevertheless, since some amount of residual noise from the debris disk remains in the high-passfiltered images, we always show two kinds of detection limits: (1) those in a region far away from the disk and (2) those in a small region around the disk midplane.Then, using the ATMO + BEX hot start evolutionary models introduced in the previous paragraph, the companion mass sensitivity of our observations was derived from the contrast curves.Figure 9 shows the contrast and companion mass limits obtained for the different filters as a function of the angular separation from the host star β Pic.The contrast with respect to β Pic was obtained by comparing the flux measured in the NIR-Cam images to the flux of β Pic according to a BOSZ A6V-type stellar atmosphere model (Bohlin et al. 2017) fitted to the archival 1-5 µm photometry of β Pic obtained from VizieR16 .The contrast limits reach down to better than 1e−7 beyond ∼ 4 arcsec (∼ 80 au).In the background-limited regime, this is about 10 times deeper in the L-band and 30 times deeper in the M-band if compared to the NACO detection limits reported in Stolker et al. (2019).The ADI+RDI and ADI reductions perform fairly similar at close-in and wide separations and there is an intermediate regime (∼ 0.5-2 arcsec) where ADI+RDI performs better than ADI alone.We note that this regime is where the coronagraphic optics create a strong quasi-static speckle field in the images.
The RDI only reduction performs the worst.As discussed in Section 2.2.2, this is due to the poor reference star subtraction which resulted from the suboptimal TA performance for the reference star observations of α Pic in the LW channel.
Looking at the companion mass limits, the F444W wide-band filter outperforms all other filters by a factor of ≳ 10 at wide separations and reaches a 5σ detection limit of ∼ 0.05 M Jup (approximately a Neptune/Uranus mass) beyond ∼ 4 arcsec (∼ 80 au) away from the disk midplane.The far deeper sensitivity of the F444W filter has multiple reasons.Firstly, the wide-band filters have a significantly increased throughput compared to the medium-band filters, resulting in higher sensitivity.Secondly, the contrast between a young giant planet and an A-type host star is smaller at longer wavelengths (e.g., Beiler et al. 2023;Leggett & Tremblin 2023), which benefits the F444W filter the most.We can hence rule out the presence of companions away from the disk midplane with masses above ∼ 0.05 M Jup between ∼ 4 and 10 arcsec (∼ 80-200 au) from β Pic.In the disk midplane, the companion mass limits are worse due to residual photon noise and flux from the disk.We reach a 5σ sensitivity of ∼ 1 M Jup beyond 2 arcsec (∼ 40 au) from β Pic.(Allard et al. 2001) and Sonora Bobcat (Marley et al. 2021) isochrones.For reference, two other substellar companions that were already observed with JWST /NIRSpec spectroscopy (VHS 1256B, Miles et al. 2023, and TWA 27 B, Luhman et al. 2023;Manjavacas et al. 2024) are also shown.These diagrams can inform the evolution of substellar companions from the correlation of their luminosity and spectral type, and from correlations between the spectral slopes of two different wavelength regimes, respectively.Both diagrams were made with the species toolkit (Stolker et al. 2020).Given the dynamical mass of ∼ 9 ± 2 M Jup of β Pic b (Nowak et al. 2020;Brandt et al. 2021), the young giant planet agrees better with the AMES-Dusty isochrone, which is not surprising given that previous works have already identified β Pic b to have a cloudy atmosphere (e.g., Currie et al. 2013;Bonnefoy et al. 2013;Morzinski et al. 2015;Chilcote et al. 2017).Both isochrones assume equilibrium chemistry which plays a role for the depth of the methane feature discussed below.While Sonora Bobcat describes a cloudless atmosphere with only rainout chemistry included, AMES-Dusty completely neglects the gravitational settling of the grains so that layers of clouds automatically build up through condensation.The dust clouds absorb light at shorter wavelengths and re-emit it in the near-IR, leading to a redder color and thus an increased brightness in the near-IR if compared to a cloudless atmosphere at the same surface gravity (i.e., planet mass).In the colorcolor diagram, both isochrones predict that objects will become bluer in F300M − F444W with increasing mass, because the effective temperature also increases, moving the blackbody peak to shorter wavelengths.These colors are hence probing mostly the pseudo-continuum of the atmosphere.The combination of the F300M and the F335M filters provides a tracer of the methane absorption feature at ∼ 3.3 µm.The F300M − F335M color first becomes redder at small masses, because the depth of the methane absorption feature in the F335M filter decreases.Then, at a mass of about ∼ 7 M Jup for AMES-Dusty and ∼ 14 M Jup for Sonora Bobcat, the methane feature completely disappears and objects become bluer with increasing mass as expected from a pure blackbody.It can be seen that the giant planet β Pic b is already in the regime where the methane feature has disappeared.We note, however, that the evolutionary models used here do not include disequilibrium chemistry so that methane absorption already appears at higher effective temperatures compared to what has been observed (Yamamura et al. 2010;Zahnle & Marley 2014;Moses et al. 2016;Stolker et al. 2020).For a low-gravity L-type object such as β Pic b, methane absorption is not expected.
To derive a set of physical parameters of β Pic b based on our JWST /NIRCam photometry, we fitted model atmospheres to its observed spectrophotometry using the species toolkit (Stolker et al. 2020).Besides the new JWST /NIRCam photometry, we included the Gemini/GPI spectra from Chilcote et al. (2017), the VLTI/GRAVITY spectrum from GRAVITY Collaboration et al. ( 2020), the recent JWST /MIRI spectrum from Worthen et al. (2024), the Magellan/VisAO and Gemini/NICI photometry from Males et al. (2014), and the VLT/NACO photometry from Bonnefoy et al. (2011), Currie et al. (2013), Stolker et al. (2019), andStolker et al. (2020).We explored three different model atmosphere grids: the DRIFT-PHOENIX (Helling et al. 2008), the Exo-REM (Baudino et al. 2015;Charnay et al. 2018), and the BT-Settl (Allard et al. 2012) grid.All of them assume radiative-convective equilibrium to calculate the temperature structure of the atmosphere and include photospheric absorption by dust clouds.Cloudy model atmospheres have been found to provide good fits to the observed spectrophotometry of β Pic b by many previous works (e.g., Currie et al. 2013;Bonnefoy et al. 2013;Stolker et al. 2020;GRAVITY Collaboration et al. 2020).A cloudy atmosphere is also expected for β Pic b based on its effective temperature and the condensation temperatures of common refractory species in its atmosphere (e.g., Marley & Robinson 2015).To be able to fit the different model grids to the data, we first bin them to the spectral resolution of the respective instrument and calculate synthetic photometry for all considered filters.Then, we interpolate the binned spectra and synthetic photometry linearly between the different grid points and infer the posterior distribution of the grid parameters using nested sampling with the PyMultiNest17 package (Buchner et al.

2014
). Uniform distributions were chosen for the priors in the parameter ranges shown in Table 5.
To investigate the characterization power of NIRCam by itself, we perform two types of fits: (1) using only the new NIRCam photometry and (2) including the existing literature data in addition to the new NIRCam photometry.Given that β Pic b is such a well-studied object, these tests are ideally suited to check which atmospheric parameters our chosen set of NIRCam filters can constrain by themselves and which not.This knowledge will be especially interesting for future NIRCam observations of previously unknown exoplanetary companions.The best fit parameter values for the three considered model atmosphere grids and the two types of fits are shown in Table 6.We note that fits using only the existing literature data without the new NIRCam photometry do not significantly differ from the ones including the new NIRCam photometry and are hence not shown or further discussed here.There are multiple reasons for this.Firstly, the GPI and GRAVITY spectra have so many data points that the weight of the NIRCam photometry is very small.Secondly, there is already NACO photometry in the L-and M-bands and the error bars of the NIRCam photometry are too large to further constrain the best fit models in this wavelength regime.
Figure 11 shows the best fit model atmospheres for β Pic b for all six considered cases.The fits that consider only the new NIRCam photometry favor an object which emits less flux over the shown spectral range and has a smaller bolometric luminosity.Notably, the NACO photometry and the GPI H-band spectrum appear to fit better to the model atmosphere derived from the NIRCam photometry alone.This means that the NACO and GPI H-band data seem to be more con-sistent with the new NIRCam photometry than with the existing GRAVITY spectrum, which pulls the entire near-IR spectrum of β Pic b up to higher fluxes in the fits that include the literature data.We note that these fits use an equal weight for each data point, so that the GRAVITY spectrum, although correlations between spectral channels are taken into account, dominates the fit.Investigating the origin of the systematic differences between the GRAVITY spectrum and the other instruments is beyond the scope of this paper, but we note in Section 3.1 that variability in the scattered light from the debris disk could play a role.
Comparing, for each model, the fits with the NIR-Cam data alone to the ones including existing literature data, Figure 12 shows that the retrieved posteriors for the effective temperature, the surface gravity, and the radius are mostly consistent with one another within 1σ.As observed frequently in the literature, there are however systematic differences between the different model grids.In general, the retrieved effective temperature ranges from ∼ 1450-1750 K, which is consistent with previous studies (e.g., Bonnefoy et al. 2014b;GRAV-ITY Collaboration et al. 2020).However, it becomes evident that at the current photometric precision our specific set of NIRCam filters is not suitable to constrain the planet's metallicity or C/O abundance ratio.Figure 12 illustrates that the posterior distributions of these quantities remain unconstrained regardless of the chosen model atmosphere grid.When including literature data in the fits, the metallicity and C/O abundance ratio posterior distributions converge sharply.With the Exo-REM grid, we retrieve a C/O abundance ratio of 0.37 +0.01  −0.01 which is consistent with the value retrieved in GRAVITY Collaboration et al. (2020) with the GRAV-  ITY K-band spectrum alone.The retrieved metallicity is strongly model-dependent.We note that we can obtain much tighter constraints on the metallicity and C/O abundance ratio from our NIRCam observations alone when artificially increasing the photometric precision to ∼ 1%.When doing so, the uncertainty on the retrieved metallicity drops to ∼ 0.05 dex and the uncertainty on the retrieved C/O ratio drops to ∼ 0.1.We suspect that the six photometric points from NIRCam alone can constrain these quantities quite tightly because they cover such a broad wavelength range (∼ 1.7-5 µm).However, similar as above, there remain systematic differences between the different model grids so that the retrieved parameter values need to be treated with caution.For all three model grids, though, we find that the luminosity retrieved from the NIRCam photometry alone is significantly smaller than the one retrieved from the fit including literature data.This is not completely surprising given that the new NIRCam photometry falls consistently below the best fit DRIFT-PHOENIX model atmosphere fitted to literature data shown in Figure 5.Our specific set of NIRCam filters covers both the water  6).The right panel shows the radius, with the constraints from the three model atmosphere grids highlighted separately.The samples from NGPPS were selected from NG73 (i.e., one planetary embryo per simulated system) at an age of 20 Myr.The gray area shows the dynamical mass of β Pic b, indicating a discrepancy with the bulk constraints from the SED.
bands at ∼ 2 and ∼ 2.5 µm and the CO/CO 2 bands in the ∼ 4-5 µm regime, which are the main carbon-and oxygen bearing species in the atmosphere of β Pic b.However, we are lacking medium band observations between 4-5 µm which could provide a better handle on the relative strength of the CO and CO 2 features.

The formation of β Pic b
The formation history of β Pic b has been investigated in several studies, and most of them concluded that either its luminosity-age relationship (e.g., Bonnefoy et al. 2014b;Marleau & Cumming 2014;Nowak et al. 2020) or its atmospheric composition (e.g., GRAVITY Collaboration et al. 2020) are most consistent with a moderate formation entropy (≳ 9.75 k B /baryon) and suggesting a formation by warm-start core-accretion (Pollack et al. 1996;Spiegel & Burrows 2012;Mordasini 2013) followed by strong planetesimal enrichment in the inner regions of the β Pic system (i.e., within the CO 2 iceline, e.g., Öberg et al. 2011;Eistrup et al. 2018).Based on the new JWST /NIRCam photometry which improves the constraints on β Pic b's thermal emission especially in the ∼ 3-5 µm regime, we obtain slightly updated constraints on the planet's bolometric luminosity and radius.We note that while the ∼ 5-7 µm MIRI MRS spectrum of β Pic b is included in the fits, it does not contribute significantly to the estimated bolometric luminosity of the planet because we let the spectrum scale up or down freely in these fits, following the same procedure as in Worthen et al. (2024).Moreover, as in Stolker et al. (2020), we keep the GRAVITY and GPI Y-J-band spectra fixed while also fitting a scale factor for the GPI H-band spectrum which does otherwise not agree with the GRAVITY spectrum.A similar discrep-ancy between the GPI H-band spectrum and the GRAV-ITY and SPHERE spectra of the brown dwarf companion HD 206893 B has been reported by Kammerer et al. (2021).
We compare these updated values to evolutionary model predictions from the New Generation Planetary Population Synthesis (NGPPS) simulation (Emsenhuber et al. 2021) in Figure 13.The NGPPS simulation uses a global end-to-end planetary formation and evolution model based on the core-accretion process to predict all common observational quantities of exoplanets and exoplanetary systems.It can be seen that the dynamical mass from Nowak et al. (2020) and Brandt et al. (2021) is in tension with the luminosity measured from the SED fits which suggests that β Pic b has a mass closer to the deuterium burning limit.This was already noted by Brandt et al. (2021).The same conclusion can be drawn from the inferred planet radius, which (regardless of the used model atmosphere grid) is also in tension with the one inferred from the NGPPS and the dynamical mass of β Pic b.Investigating the reason for this discrepancy is beyond the scope of this paper.However, we note that Lacour et al. ( 2021) find a more consistent dynamical mass of 11.9 ± 3 M Jup for β Pic b when only considering the radial velocities and the astrometry, ignoring the Gaia-Hipparcos proper motion anomaly.
Our inferred bolometric luminosity is consistent with the one reported in Morzinski et al. (2015) and Chilcote et al. (2017) to within ±0.05 dex, however, the comparison of dynamical mass and NGPPS predicts a bolometric luminosity that is ∼ 0.5 dex lower.We further note that the radius inferred by the Exo-REM grid is very large, almost unphysical if compared to the NG-GPS population.However, it is roughly consistent with the radius of 1.9 R Jup inferred by GRAVITY Collaboration et al. (2020) also using the Exo-REM grid.With the free retrieval code petitRADTRANS, the same authors obtain a radius of 1.36 ± 0.02 R Jup which is much closer to the one that we infer using the DRIFT-PHOENIX model atmosphere grid.The large discrepancies between the inferred radii can originate from the different cloud prescriptions in the different model atmosphere grids and retrieval codes.Clouds can have a similar effect on a planet's SED as a change in effective temperature and radius and are therefore often correlated with these quantities.In the literature, the DRIFT-PHOENIX grid and the cloud prescription therein has been found to provide good fits to young and dusty lowgravity objects (Patience et al. 2012;Bonnefoy et al. 2014a;Lachapelle et al. 2015) such as β Pic b (Stolker et al. 2020), HD 95086 b (De Rosa et al. 2016), and HD 206893 B (Kammerer et al. 2021).

SUMMARY AND CONCLUSIONS
• We observed the β Pic exoplanetary system with JWST /NIRCam coronagraphy.We detected the directly-imaged giant planet β Pic b in all six observed NIRCam filters.
• We measure new photometry for β Pic b which mostly agrees with existing data from the ground.We observe a ∼ 10-20% discrepancy with respect to the GRAVITY spectrum in the K-band.We note that a similar discrepancy can be seen between the NACO and the GRAVITY data, and that our NIRCam photometry is consistent with NACO within the uncertainties.As it was found in previous studies, our new NIRCam photometry supports that β Pic b has an atmosphere composed of a thick layer of clouds.
• We measure new astrometry for β Pic b which is still affected by significant systematic errors.
Given the unprecedented accuracy of the existing VLTI/GRAVITY astrometry, our new NIRCam astrometry does not provide updated constraints on the orbital parameters or the dynamical mass of β Pic b.Making the JWST astrometry competitive with measurements from single-dish groundbased telescopes requires further work.While the implementation of an updated distortion correction is already being worked on by the instrument and spaceKLIP teams, the largest systematic error is introduced by the unknown position of the host star behind the coronagraphic mask.Observers who care about accurate companion astrometry are hence encouraged to take astrometric confir-mation full frame images and use nearby field stars to interpolate the position of the occulted host star behind the coronagraphic mask.
• Thanks to the unparalleled infrared sensitivity of JWST especially in the NIRCam F444W filter, we are able to rule out additional companions in the disk midplane above ∼ 1 M Jup beyond 40 au from the star.Away from the disk midplane, the companion mass limits reach down to ∼ 0.05 M Jup beyond 80 au from the star.In particular, we can rule out companions more massive than ∼ 1 M Jup in the disk midplane and ∼ 0.1 M Jup away from the disk midplane outward of the location of the southwestern clump in the debris disk (e.g., Han et al. 2023;Skaf et al. 2023;Rebollido et al. 2024).We further identify five additional localized sources in the NIRCam LW data, but all of them are found to be background stars or galaxies based on their colors or spatial extent.
• By fitting model atmosphere grids to the new NIR-Cam photometry alone and combined with existing data from the literature, we determine the physical and atmospheric parameters of β Pic b.We found that our NIRCam observations are not able to constrain the metallicity and C/O abundance ratio of the planet.However, the effective temperature, surface gravity, and radius inferred from the NIRCam only fits are typically consistent with the ones inferred from existing literature data within the uncertainties.The bolometric luminosity inferred from the NIRCam only fits is slightly smaller than the one inferred from existing literature data.
• The presence of the extended debris disk around β Pic makes the extraction of accurate planet photometry for β Pic b challenging.We used fake companion injection and recovery tests to address this issue, but these assume that the debris disk and residual speckles are symmetric with respect to the star.Ultimately, we used MCRDI techniques to simultaneously fit for the disk and the planet, but also the accuracy of this approach depends on how accurately the MCRDI disk model is able to resemble the data.These complications illustrate that measuring accurate planet photometry in the presence of an exozodiacal dust disk in the era of the Habitable Worlds Observatory is a potential problem that deserves further attention (e.g., Kammerer et al. 2022a;Currie et al. 2023).
• The known inner planet β Pic c remains undetected behind the occulting spot of the NIRCam coronagraphs.The direct detection of β Pic c with JWST would require the use of high-resolution imaging techniques such as NIRISS Aperture Masking Interferometry (AMI, Sivaramakrishnan et al. 2023), although even for AMI this planet is challenging as its contrast is expected to be close to the reported limit of ∼ 1e−4 in the L-and Mband.We note that NIRISS Kernel Phase Interferometry (KPI, Kammerer et al. 2023) is not feasible for this target because β Pic would highly saturate the detector.Injected fake companion (F444W) to recover it using the same routines that have also been used to extract the astrometry and photometry of β Pic b.Since we forward model the companion PSFs, they already take into account KLIP algorithm throughput losses and companion self-subtraction due to ADI.However, residual flux from the debris disk and remaining stellar speckles may still affect the measured companion properties.Since the debris disk around β Pic is highly symmetric in the NIRCam bands (see Rebollido et al. 2024), the ideal position to inject the fake companion is at the same angular separation as β Pic b, but on the opposite side of the star.This allows us to study the impact of the disk on the measured fake companion properties in an environment that is very similar to that of the true companion β Pic b.The flux at which we inject the fake companion is always the flux that was measured for β Pic b from the MCRDI reductions.
Figure 14 shows the results from these fake companion injection and recovery tests using KLIP ADI+RDI, ADI, and RDI techniques.With RDI alone, the recovered companion flux without high-pass filtering is typically higher than the injected flux, because flux from the debris disk adds on top of the companion flux and positively biases the companion flux.With ADI alone, it is the other way around, because disk self-subtraction adds on top of the planet flux and negatively biases the companion flux.If high-pass filtering is included, the recovered companion flux often depends strongly on the high-pass filter size.RDI alone yields either too high or too low flux.ADI alone yields reasonable flux in the SW channel, but too low flux in the LW channel.ADI+RDI yields too high flux in the SW channel, but reasonable flux in the LW channel.Altogether, we adopt an ideal high-pass filter size of 7 pixels for ADI in the SW channel and of 3 pixels for ADI+RDI in the LW channel.While we could in principle choose a different high-pass filter size for each NIRCam filter, this would unnecessarily complicate the analysis given that the recovered fluxes for our adopted parameters are consistent with the recovered fluxes for the ideal high-pass filter sizes within the uncertainties.

B. LITERATURE PHOTOMETRY OF β Pic b
Table 7 summarizes the literature photometry and spectra of β Pic b used for model atmosphere fitting in this work.We also quote the epochs at which the data has been observed.This might be relevant if variability in the debris disk affects the measured planet fluxes.

Figure 1 .
Figure1.Reference star TA performance for the SW (top) and the LW (bottom) observations.The dots are color-coded by dither position and show the measured offsets between the reference star and the science target observation (first roll), while the crosses show the commanded and expected offsets in case of a perfect TA.The TA performance is good for the SW channel, but worse than expected for the LW channel (see Section 2.2.2 for details).The plots show the F210M and F335M datasets; all other filters in the same NIRCam detector channels show similar performance.

Figure 2 .
Figure 2. JWST /NIRCam coronagraphy images of the β Pic system.The six rows show the six observed filters and the three columns show the pre-KLIP images (non-PSF-subtracted, left), the post-KLIP images (PSF-subtracted, middle), and the post-KLIP and high-pass filtered images (high-pass filter size of 3 pixels, right).The host star is located at the origin.All images are shown in the same linear color stretch.The small insets in the middle and right columns show a zoom on the giant planet β Pic b in a different color stretch that aims to highlight the planet PSF.

Figure 3 .
Figure 3. MCRDI-modeling of β Pic b.The six rows show the six observed filters and the three columns show the host star PSF-and disk model-subtracted images (left), the best fit planet PSF models (middle), and the residuals between the two (right).The position of the host star β Pic is indicated with a black star.F250M and F300M are particularly affected by spatial undersampling and the PSF subtraction has higher residuals for these filters.
shows that the new NIRCam photometry agrees well with a DRIFT-PHOENIX model atmosphere fitted to previously exist-

Figure 4 .
Figure 4. Forward-modeling of β Pic b's PSF in the KLIP approach.The six rows show the six observed filters and the three columns show the KLIP-subtracted images (left), the best fit forward model planet PSFs (middle), and the residuals between the two (right).The position of the host star β Pic is indicated with a black star.
JWST /NIRCam photometry of the giant planet β Pic b.While the flux values were directly measured from the flux-calibrated NIRCam images, the contrast (∆mag) was computed relative to a stellar model atmosphere of β Pic and thus contains an additional systematic uncertainty.TPMSK and TPCOM denote the throughput of the coronagraphic mask at the position of the planet and the throughput of the coronagraphic mask substrate, respectively.

Figure 5 .Figure 6 .
Figure 5. JWST /NIRCam KLIP (top) and MCRDI (bottom) photometry of the giant planet β Pic b in red plotted together with other photometry and spectra and a DRIFT-PHOENIX model atmosphere of β Pic b from the literature.The filled symbols show the data and the open symbols show the values predicted by the best fitting model atmosphere.The model atmosphere was fitted to the literature data only (and not the new NIRCam data) and its parameters are printed at the top of the middle panel.The residuals between the data and the model atmosphere are shown in the bottom panel and the filter transmission curves are shown in the top panel.JWST /NIRCam coronagraphy has demonstrated exquisite sensitivity to new companions at wide separa-

Figure 7 .
Figure 7. PSF-subtracted JWST /NIRCam (left) and MIRI (right) images of the β Pic system (arbitrary color stretch).Five background sources were identified based on their NIRCam colors or morphology whose positions are highlighted in both images with red circles.The orange circles in the NIRCam image show the source positions rotated by the telescope roll angle, where negative ADI residuals are visible.BG4 appears to be the only source which is also visible in the MIRI image.The MIRI image was adapted from Rebollido et al. (2024).

Figure 8 .
Figure 8. JWST /NIRCam color-magnitude diagrams showing an ATMO + BEX hot start evolutionary track for an 18.5 Myr old bound companion in blue.The point-like background source identified in the NIRCam data of β Pic and the giant planet β Pic b are shown by an orange and a red star.The dashed and transparent blue line shows an ATMO + BEX cold start evolutionary track for reference.arytrack for an 18.5 Myr old bound companion in blue.For reference, β Pic b is also shown by a red star.The evolutionary track was obtained by combining ATMO equilibrium chemistry models fromPhillips et al. (2020) with BEX hot and cold start models fromLinder et al. (2019).This was necessary to cover the entire planet

Figure 9 .
Figure9.5σ contrast limits (top) and companion mass limits (bottom) measured from the PSF-subtracted JWST /NIRCam coronagraphy images of the β Pic system.The left panels show the ADI+RDI limits for all six observed filters; the right panels compare the three employed KLIP PSF subtraction techniques (ADI+RDI, ADI, RDI) for the F444W wide-band filter.All curves were computed from high-pass-filtered images to remove the flux of the edge-on debris disk around β Pic as much as possible.The solid lines show the limits away from the disk midplane and the transparent lines show the limits in the disk midplane.The companion mass limits were computed using ATMO + BEX hot start evolutionary models and assuming an age of 18.5 Myr and a distance of 19.44 pc for β Pic.The masses of the gas and ice giants in the Solar System are shown for reference (planet symbol sizes are not to scale, and planet projected separations are five times the true separations from the Sun).

Figure 10 .
Figure 10.JWST /NIRCam color-magnitude diagram (left) and color-color diagram (right) showing Sonora Bobcat and AMES-Dusty isochrones for an 18.5 Myr old companion in blue and orange, respectively.The giant planet β Pic b is shown by a red star and two other substellar companions that were already observed with JWST spectroscopy are shown by pink stars.The dashed gray line shows a blackbody and the dashed blue and orange lines show the isochrones for a 100 Myr old companion for reference.
Figure10shows the position of β Pic b in a NIR-Cam color-magnitude and color-color diagram together with AMES-Dusty(Allard et al. 2001) and Sonora Bobcat(Marley et al. 2021) isochrones.For reference, two other substellar companions that were already observed with JWST /NIRSpec spectroscopy(VHS 1256B, Miles et al. 2023, and TWA 27 B, Luhman et al. 2023;Manjavacas et al. 2024) are also shown.These diagrams can inform the evolution of substellar companions from the correlation of their luminosity and spectral type, and from correlations between the spectral slopes of two different wavelength regimes, respectively.Both diagrams were made with the species toolkit(Stolker et al. 2020).Given the dynamical mass of ∼ 9 ± 2 M Jup of β Pic b(Nowak et al. 2020;Brandt et al. 2021), the young giant planet agrees better with the AMES-Dusty isochrone, which is not surprising given that previous works have already identified β Pic b to have a cloudy atmosphere (e.g.,Currie et al. 2013;Bonnefoy et al. 2013;Morzinski et al. 2015;Chilcote et al. 2017).Both isochrones assume equilibrium chemistry which plays a role for the depth of the methane feature discussed below.While Sonora Bobcat describes a cloudless atmosphere with only rainout chemistry included, AMES-Dusty completely neglects the gravitational settling of the grains so that layers of clouds automatically build up through condensation.The dust clouds absorb light at shorter wavelengths and re-emit it in the near-IR, leading to a redder color and thus an increased brightness in the near-IR if compared to a cloudless atmosphere at the same surface gravity (i.e., planet mass).In the colorcolor diagram, both isochrones predict that objects will become bluer in F300M − F444W with increasing mass, because the effective temperature also increases, moving the blackbody peak to shorter wavelengths.These colors are hence probing mostly the pseudo-continuum of the atmosphere.The combination of the F300M and the F335M filters provides a tracer of the methane absorption feature at ∼ 3.3 µm.The F300M − F335M color first becomes redder at small masses, because the depth of the methane absorption feature in the F335M

TTTFigure 11 .
Figure 11.Best fit model atmospheres from the BT-Settl (top), DRIFT-PHOENIX (middle), and Exo-REM (bottom) grids fitted to the JWST /NIRCam photometry of the giant planet β Pic b alone (light gray) and together with existing data from the literature (dark gray).

Figure 12 .
Figure 12.Model atmosphere parameter posterior distributions for the BT-Settl (top), DRIFT-PHOENIX (middle), and Exo-REM (bottom) grids fitted to the JWST /NIRCam photometry of the giant planet β Pic b alone (orange) and together with existing data from the literature (blue).The markers with horizontal error bars above the histograms show the 16th, 50th, and 84th percentile of the distribution.

Figure 13 .
Figure 13.Comparison of β Pic b's inferred luminosity and radius with the New Generation Planetary Population Synthesis (NGPPS; Emsenhuber et al. 2021).The left panel shows the bolometric luminosity, with the orange area highlighting the combined constraints from the SED fits (all six models in Table6).The right panel shows the radius, with the constraints from the three model atmosphere grids highlighted separately.The samples from NGPPS were selected from NG73 (i.e., one planetary embryo per simulated system) at an age of 20 Myr.The gray area shows the dynamical mass of β Pic b, indicating a discrepancy with the bulk constraints from the SED.

Figure 14 .
Figure 14.Recovered companion fluxes as a function of the high-pass filter size from fake companion injection and recovery tests using ADI+RDI (solid lines), ADI (dashed lines), and RDI (dotted lines) techniques.The injected fluxes (ground truth) are shown by a horizontal solid line.The fake companion was injected at the same angular separation as β Pic b, but on the opposite side of the star.

Table 1 .
Target Filter Readout pattern Dither pattern Nints/Ngroups/N frames Texp [s] PA [deg] Summary of the JWST /NIRCam observations of the young giant planet β Pic b presented in this paper (JWST program ID 1411).The science target is β Pic and the PSF reference target is α Pic.The observations were taken on 18 March 2023.The total exposure time Texp sums up all five dither positions for the reference target.

Table 4 .
Astrometry and photometry of the detected background sources.Only the values for BG1 are based on an accurate fit with a model PSF.The four other sources appear extended and their astrometric positions are the approximate centroids.

Table 5 .
Prior ranges adopted for the considered model atmosphere grids.The priors were distributed uniformly within the quoted ranges.

Table 6 .
Best fit bulk parameters inferred for the young giant planet β Pic b using different model atmosphere grids.For each model grid, we separately report the results of fits to all available data, and to the NIRCam data only.