MIRI MRS Observations of β Pictoris. I. The Inner Dust, the Planet, and the Gas

We present JWST MIRI Medium Resolution Spectrograph (MRS) observations of the β Pictoris system. We detect an infrared excess from the central unresolved point source from 5 to 7.5 μm which is indicative of dust within the inner ∼7 au of the system. We perform point-spread function (PSF) subtraction on the MRS data cubes and detect a spatially resolved dust population emitting at 5 μm. This spatially resolved hot dust population is best explained if the dust grains are in the small grain limit (2πa ≪ λ). The combination of unresolved and resolved dust at 5 μm could suggest that dust grains are being produced in the inner few astronomical units of the system and are then radiatively driven outwards, where the particles could accrete onto the known planets in the system, β Pictoris b and c. We also report the detection of an emission line at 6.986 μm that we attribute to [Ar ii]. We find that the [Ar ii] emission is spatially resolved with JWST and appears to be aligned with the dust disk. Through PSF-subtraction techniques, we detect β Pictoris b at the 5σ level in our MRS data cubes and present the first mid-infrared spectrum of the planet from 5 to 7 μm. The planet’s spectrum is consistent with having absorption from water vapor between 5 and 6.5 μm. We perform atmosphere model grid fitting of the spectra and photometry of β Pictoris b and find that the planet’s atmosphere likely has a substellar C/O ratio.


Introduction
Debris disks are planetary systems that consist of dust, gas, planetesimals, and planets that typically correspond to the late stages of planetary system formation (Wyatt 2008;Hughes et al. 2018).They provide a unique laboratory to study the processes involved in the later stages of planet formation and evolution.Unlike protoplanetary disks, the dust seen in debris disks is thought to be constantly replenished through collisional processes between minor bodies in the system (Hughes et al. 2018).This is because the dust grains in debris disks are subject to radiation pressure and Poynting-Robertson drag that remove dust from orbits around the central star on timescales that are short compared to the age of the system (Guess 1962;Krivov et al. 2006;Wyatt 2008).The detection of these dust grains in debris disks points to ongoing stochastic and steady-state collisions between planetesimals that actively replenish the small dust grains (Backman & Paresce 1993).
The particles in debris disks range in size from the parent planetesimals to the collisionally produced dust, but it is the dust that is observable in both thermal emission at mid-infrared to millimeter wavelengths (e.g., Holland et al. 1998;Koerner et al. 1998;Telesco et al. 2005;MacGregor et al. 2018) and scattered light at optical and near-infrared wavelengths (e.g., Smith & Terrile 1984;Kalas et al. 2006;Esposito et al. 2020;Ren et al. 2023).Infrared spectra and the spectral energy distributions (SEDs) of debris disks can reveal the temperature of the dust (e.g., Ballering et al. 2013), which depends on the location, size, and composition of the dust particles.SED modeling as well as multiwavelength imaging have revealed some debris disks to have multiple populations of dust at various stellocentric distances (e.g., Chen et al. 2014;Jang-Condell et al. 2015;Gáspár et al. 2023).Knowledge of the innermost regions of debris disks has been mostly limited to infrared spectroscopic and photometric analyses, mainly because of the low spatial resolution of previous space-based infrared observatories like Spitzer.JWST provides a unique opportunity to study the spectra and structure of dust in debris disks at high angular resolution with the MIRI Medium Resolution Spectrograph (MRS).
β Pictoris (β Pic hereafter) is a ∼23 Myr (Mamajek & Bell 2014) A6V star that is host to the first ever imaged debris disk (Smith & Terrile 1984).The β Pic disk is oriented close to edge on from our line of sight and has been studied in scattered light as well as thermal emission at mid-infrared and millimeter wavelengths (e.g., Telesco et al. 2005;Golimowski et al. 2006;Matrà et al. 2019;Rebollido et al. 2024).At a distance of 19.6 pc (Gaia Collaboration et al. 2023), β Pic provides a great laboratory for studying the spatial structure of its debris disk.Scattered light imaging with the Hubble Space Telescope revealed a warp in the inner disk (Golimowski et al. 2006) that was later attributed to interactions with the now confirmed Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
giant planet β Pic b (∼9 M J , a = 9.9 au; Lagrange et al. 2009;Dawson et al. 2011;Lagrange et al. 2012;GRAVITY Collaboration et al. 2020;Nowak et al. 2020).Radial velocities provided evidence for the presence of a second planet in the system, β Pic c (Lagrange et al. 2019; ∼8 M J , a = 2.7 au), which was later directly confirmed by interferometric observations with the Very Large Telescope (VLT)/GRAV-ITY (Nowak et al. 2020).
As a young nearby system with both giant planets and a debris disk, β Pic provides a unique opportunity to study the interactions between dust in the disk and the giant planets in the system.Giant planets present in debris disks are thought to impart structures on the dust in debris disks (e.g., Dawson et al. 2011;Crotts et al. 2021), however, it is unclear if and to what extent the dust and minor bodies in debris disks can affect giant planets via accretion.Giant planets likely form within ∼10 Myr, before the debris disk evolutionary stage (Williams & Cieza 2011;Li & Xiao 2016), although it is still possible for material from debris disks to accrete onto giant planets (e.g., Frantseva et al. 2020;Kral et al. 2020).In our solar system for instance, the impact of comet Shoemaker-Levy 9 in 1994 delivered refractory material to Jupiter's atmosphere (Harrington et al. 2004;Fletcher et al. 2010).Marley et al. (2012) suggested that micron-and submicron-sized silicate grains can remain at higher altitudes in the atmospheres of young, low surface gravity giant planets and brown dwarfs and thus affect their observed spectra.If young giant planets are accreting dust from their debris disks, it is then possible for the dust grains to remain in their atmosphere and affect the observed spectrum of the planet.For young giant planets that are in systems with debris disks, the amount of dust that is accreted from the disks onto the planets is not well observationally constrained.
Mid-infrared spectroscopy of β Pic has revealed the presence of small submicron silicate grains in the system (Okamoto et al. 2004;Li et al. 2012;Lu et al. 2022).By modeling silicate features at 10, 18, and 23 μm detected with Spitzer IRS, Lu et al. (2022) found that the features are best produced by dust grains in the Rayleigh limit (2πa = λ), indicating the presence of submicron-sized silicate grains that are subject to blowout from radiation forces.Depending on the location of the dust, these subblowout-sized grains have the potential to interact with the known planets in the system as they are being radiatively driven outwards.Lu et al. (2022) also found tentative evidence for an infrared excess at 5 μm and the presence of ∼600 K dust in the system.Because of the limited angular resolution of Spitzer, the location of this hot dust was not clear.
In this paper, we present JWST MIRI MRS observations of the β Pic system, which provide higher spatial and spectral resolution space-based data in the mid-infrared than previously obtained with Spitzer.The observations and data reduction are described in Section 2. In Section 3, we present our main findings including (1) an infrared excess from 5 to 7.5 μm, (2) the discovery of a spatially resolved hot dust population, (3) the first detection of [Ar II] at 6.986 μm in the β Pic disk, and (4) the extraction of a low-resolution MRS spectrum of β Pic b.In Section 4, we discuss the implications of the new spatially resolved hot dust population and how, along with the 5 μm excess, it provides evidence for an outflowing wind of small dust grains that could be blown onto the known planets where they may be accreted.In Section 5, we state our conclusions and summarize.

Data Acquisition
As a part of GTO program 1294, we observed β Pic (K = 3.48; Ducati 2002) with MIRI MRS (Wells et al. 2015;Argyriou et al. 2023) on 2023, January 11, preceded by dedicated background observations and followed immediately by observations of the nearby star N Car (A0II, K = 4.218 ;Houk 1978;Cutri et al. 2003), which we used as a point-spread function (PSF) reference star.N Car is offset from β Pic by 7°.6.All observations were taken in all three grating settings (short, medium, and long) in all four channels with the FASTR1 readout pattern to cover the wavelength range 4.9-28 μm.We used a four-point point-source dither pattern with the negative dither orientation10 to observe β Pic and N Car.We observed both stars with target acquisition to ensure that both stars were well aligned on the detector to optimize reference PSF subtraction.Both β Pic and N Car were used themselves as target acquisition stars.The aperture position angle of JWST was 23°.8 for the observations of β Pic to make sure that the MRS image slicers were orthogonal to the edge-on disk.This was to ensure that the MRS cross artifact (Argyriou et al. 2023) would not mimic disk signatures in the integral field unit (IFU) data cubes.
For the dedicated background observations, we used a twopoint dither pattern.The position of the dedicated background field was R.A., decl.= 05 47 6.3034, −51 02 55.85, which is 1 7 offset from β Pic.The exposure time for N Car was longer than for β Pic to ensure that the N Car observations had a similar signal-to-noise ratio (S/N) as the β Pic observations.For β Pic, the number of groups and integrations for each dither position for channels 1 and 2 (4.9-11.71μm) was 5 and 14, respectively (exposure time = 921 s).The number of groups and integrations for each dither position for channels 3 and 4 (11.55-28.1 μm) was 15 and 5, respectively (exposure time = 877 s).For N Car, the number of groups and integrations for each dither position for channels 1 and 2 was 5 and 28, respectively (exposure time = 1854 s).The number of groups and integrations for each dither position for channels 3 and 4 was 15 and 10, respectively (exposure time = 1765 s).For the dedicated background observation, the number of groups and integrations for each dither position for channels 1 and 2 was 5 and 16, respectively (exposure time = 264 s).For channels 3 and 4, the number of groups and integrations for the background observations was 15 and 6, respectively (exposure time = 264 s).11

Data Reduction
We processed the raw detector files for β Pic and N Car using v1.11.0 of the JWST calibration pipeline (Bushouse et al. 2023) using Calibrated Reference Data System context "jwst_1094.pmap."The pipeline consists of three stages: Detector1, Spec2, and Spec3.We processed the β Pic and N Car raw files with the exact same pipeline setup.
The Detector1 stage converts the raw ramp images to uncalibrated slope images and also includes a jump detection step that flags jumps in the ramp between consecutive groups.This step mitigates ramp jumps that are often caused by cosmic-ray hits.We changed the three group rejection threshold to be 100σ in the jump detection step because the default setting in the pipeline overflags jumps in the raw data and creates artifacts in the final spectrum.The default three group rejection threshold in the pipeline was 6σ.The Spec2 pipeline takes the uncalibrated slope images and applies calibrations and instrumental corrections, including a fringe correction.The MRS is known to suffer from fringing that can have effects of up to 30% of the spectral baseline (e.g., Argyriou et al. 2020).In the Spec2 pipeline, we applied both fringe corrections available in the pipeline, the fringe flat and the residual fringe correction step, where the residual fringe correction is not turned on by default.In Spec2, we also performed the stray light subtraction step, which is necessary to subtract the cross artifact.
We used the Spec3 pipeline to combine all four dither positions into spectral cubes for each of the 12 subbands.We created the spectral cubes with the IFUalign coordinate system so that the data cubes were aligned with the MRS image slicers, and the PSFs of β Pic and N Car were well aligned for optimal PSF subtraction.The background subtraction from the dedicated background observations was done in the Spec3 pipeline along with the outlier rejection step.We built the spectral cubes using the drizzle algorithm in the pipeline (Law et al. 2023).We left the pixel size to be the pipeline default, which is 0 13 for channel 1, 0 17 for channel 2, 0 20 for channel 3, and 0 35 for channel 4. The output of the Spec3 pipeline was calibrated spectral cubes, which we used for the analysis in the following sections.

Point-spread Function Subtraction
To search for resolved disk structure as well as β Pic b in our MRS data, we used N Car to perform classical reference differential imaging PSF subtraction on the calibrated spectral cubes.Because of position repeatability issues of the dichroic grating wheel assembly (DGA) of the MRS (Patapis et al. 2023), the science star β Pic and reference star N Car were not exactly centered in the same location on the detector.The position repeatability of the DGA wheel is ∼30 mas radially (Patapis et al. 2023).To achieve the best contrast performance from PSF subtraction, we aligned the PSF reference N Car with the PSF of β Pic prior to subtracting.
We first measured the position of the center of the PSF of β Pic and N Car in the spectral cubes by fitting a 2D Gaussian to each wavelength slice and taking the median of the best-fit center of all the wavelength slices.We did this for each subband separately (1A, 1B, etc.) because not all subbands were observed simultaneously, leading to position changes of the PSF in the spectral cube between the different subbands.We then calculated the difference between the position of the center of the PSF of β Pic and N Car and shifted and interpolated the N Car PSF to the center location of the β Pic PSF using the SciPy ndimage.shiftfunction.After aligning the N Car PSF to the β Pic PSF, we scaled the N Car PSF to have the same spectrum of the unresolved point source of β Pic.We scaled to the spectrum of β Pic rather than a photosphere model because there is an unresolved infrared excess in the central point source as shown in Figure 1, and scaling to just the photosphere model does not completely subtract the β Pic PSF.
Once the N Car PSF was aligned and scaled to the flux level of β Pic, we subtracted the N Car PSF from β Pic in each wavelength slice of the spectral cube, producing a PSFsubtracted spectral cube.We were left with large PSFsubtraction residuals at the center location of the star; so to search for faint resolved disk structure as well as β Pic b, we masked out spaxels within a 3 spaxel radius of the center of the PSF in channel 1 where the spaxel size is 0 13.We also binned every 100 image slices in the data cube to increase the S/N of the individual image slices.

Spectrum of the Central Point Source
To extract the spectrum of the unresolved point source of β Pic, we performed aperture photometry on each wavelength slice of the spectral cubes.To find the center of the PSF, we collapsed each spectral cube along the wavelength axis and fit a 2D Gaussian to the collapsed PSF.We placed the center of the circular aperture for each wavelength slice at the center of the PSF determined from the Gaussian fit.We then used a wavelength-dependent aperture radius of 1.5 times the FWHM of the PSF radius at each wavelength where we determined the FWHM of the PSF by taking a line cut through the center of the PSF in the direction orthogonal to the disk.We fit a Gaussian to the line cut of the PSF to measure the FWHM.At 5 μm, the FWHM of the PSF from the Gaussian fit was 0 325.This is larger than a diffraction-limited PSF at 5 μm, which would have a FWHM of 0 2. Law et al. (2023) similarly found that the size of the MRS PSF at 5 μm is larger than what is expected from a diffraction-limited PSF.After we extracted the flux within the circular aperture in each wavelength slice of the spectral cube, we applied an aperture correction to account for flux missed outside the extraction aperture.We applied an aperture correction for an extraction radius of 1.5 times the PSF FWHM derived in Argyriou et al. (2023) to the extracted spectra.We also applied a correction to the 1D spectrum, as described in Gasman et al. (2023), for the spectral leak at 12.2 μm.
After performing aperture photometry on the spectral cubes, there were still noise sources, like fringing, present in the extracted spectra that were not completely removed from the pipeline processing.To refine the calibration systematics, mitigate noise sources, and obtain the highest S/N spectrum of the unresolved point source of β Pic, we used the observations of N Car to create an RSRF and applied it to the spectrum of the unresolved point source of β Pic.The RSRF was also used to align the different subbands.
We extracted the spectrum of N Car using the exact same method as for β Pic described above.We then fit the UBVGJHK s (Cutri et al. 2003;Reed 2003;Gaia Collaboration 2018) photometric points for N Car with a T = 8800 K, log (g) = 4.0, and [Fe/H] = 0 BT-NextGen stellar photosphere model (Allard et al. 2011) using synthetic photometry and then divided this photosphere model by the N Car spectrum.This removed the stellar continuum and photospheric absorption lines from the N Car spectrum and provided an RSRF (consisting of detector effects present in the spectra that are present after pipeline processing).We then multiplied the RSRF into the β Pic spectrum to remove the noise sources captured in the RSRF.In channel 1 (∼5-7.5 μm), applying the RSRF increased the S/N of the spectrum from ∼130 to ∼230.
The final MRS spectrum of the unresolved point source of β Pic is shown in Figure 1.A surprising difference between the MIRI MRS spectrum shown here and the Spitzer IRS spectrum of β Pic is that the sharp 18 μm silicate feature seen with Spitzer (Lu et al. 2022), is not present in our MRS spectrum.Here, we show the MRS spectrum and leave the analysis of the properties of the silicate feature and the discussion of the disappearance of the 18 μm silicate feature for a future paper (C.H. Chen et al. 2024, in preparation).

5 μm Infrared Excess
We compared our extracted MRS spectrum to a model of the β Pic photosphere from Lu et al. (2022) to search for an IR excess at the shortest MRS wavelengths.This photosphere model is a BT-NextGen model fit to UBVRJHK s photometry as well as an Infrared Telescope Facility (IRTF) SPEX spectrum of β Pic from 0.7 to 2.5 μm.This photosphere model is used throughout the rest of the analysis.
At 5 μm, we find an infrared excess that is 4% above the stellar photosphere model (see Figure 1), potentially indicating the presence of dust within a stellocentric radius ∼ 6.5 au (from the size of the extraction aperture) in the system.To determine the temperature of the dust emitting from 5 to 7.5 μm, we subtracted the photosphere model from the MRS spectrum, smoothed it with a boxcar smoothing filter with a kernel of three data points, and then fit the channel 1 spectrum (4.9 to 7.5 μm) with a blackbody.We also included the L-band photometric point of β Pic from Bonnefoy et al. (2013) of 3.454 ± 0.003 mag.To subtract the photosphere from the Lband point, we calculated synthetic photometry using the stellar photosphere model of β Pic and subtracted this value from the measured L-band photometry.Figure 2 shows the best-fit blackbody to the 5-7.5 μm excess.This blackbody fit gives a dust temperature of 500 ± 20 K.There appears to be some residual structure in the spectrum shown in Figure 2.This structure is likely due to an imperfect subtraction of the absorption lines in the stellar photosphere model.We think this because the structure in the spectrum (such as the dips at ∼5.6 and ∼7.5 μm) are at the same wavelength as the stellar photosphere lines shown in Figure 1.
Even if the absolute flux calibration of the MRS spectrum was incorrect by 4%, the shape of the spectrum still indicates an excess between 5.5 and 7.5 μm.If we pin the shortest wavelength of the MRS spectrum (4.9 μm) to the stellar photosphere model, there is still an ∼3% infrared excess at 5.5 μm.Performing the same photosphere subtraction and blackbody fitting analysis, but with the spectrum pinned to the photosphere model at 4.9 μm, gives a dust temperature of 370 K instead of 500 K.This still indicates the presence of dust in the inner few astronomical units of the system, just at a lower temperature.However, given that dust emission is spatially resolved at 5 μm (see Section 3.3), the excess we see in the spectrum at 5 μm is likely real and not due to uncertainties in the flux calibration.We are mainly interested in the temperature of the hot dust, so we only fit a blackbody to the 5-7.5 μm excess and exclude longer wavelengths, including the silicate feature from ∼8 to 12 μm.Fitting the MRS spectrum at wavelengths longer than 12 μm would require a twocomponent, two-temperature model, and here we are only interested in the hot component.Adding a cold dust component does not affect the blackbody fit from 5 to 7.5 μm.A complete modeling of the entire MRS spectrum including the silicate features will be done in an upcoming paper (C.H. Chen et al. 2024, in preparation).
Assuming that the dust grains are blackbodies in radiative equilibrium, and using the luminosity of β Pic of 8.13 L e (Bonnefoy et al. 2013), we calculated the blackbody location of  2022) that was best fit to a near-infrared spectrum.This shows an infrared excess from 5 to 7 μm which can be explained by hot dust in the system.The vertical dashed gray line shows the detection of [Ar II] emission at 6.986 μm.The other spikes in the spectrum do not appear to be real emission lines.This is because they only cover one wavelength bin, unlike the [Ar II] line which covers five.Furthermore, these spikes also only appear after applying the relative spectral response function (RSRF) and are well aligned to noise features seen in the spectrum of N Car.They then likely result from noise that is injected by performing the RSRF with N Car.The [Ar II] line, however, is present both before and after applying the RSRF.
the 500 K dust around β Pic to be 0.9 au.This blackbody distance is a lower limit to the stellocentric distance of the grains, as smaller dust grains can be at higher temperatures further away from the star than blackbody grains.This is further discussed in Section 3.3.

Spatially Resolved Hot Dust
With Spitzer, Lu et al. (2022) found a tentative infrared excess at 5 μm that could indicate the presence of hot dust within the inner parts of the system.Interferometric observations in H band found evidence for an infrared excess that could either be from a population of hot dust (1500 K) at a distance of less than 4 au from the star, or scattered light from the outer part of the edge-on disk (Defrère et al. 2012).Here, we found additional evidence for the 5-7.5 μm infrared excess seen by Lu et al. (2022) in the spectrum of the unresolved point source (see Figure 1).We confirmed that this population of hot dust emitting at 5 μm is spatially extended in our PSFsubtracted image cubes.
The binned and PSF-subtracted image slice at 5.2 μm is shown in Figure 3 with the beam of Spitzer and the JWST extraction aperture overlaid.The purple line shows a 5σ contour of the detected dust emission, indicating that we detect spatially resolved hot dust emitting at 5 μm at the 5σ level.Figure 3 illustrates that the excess thermal dust emission at 5 μm was spatially unresolved given the angular resolution of Spitzer, but is now spatially resolved with JWST out to ∼20 au, where the disk flux density drops below the 5σ level.By applying PSF subtraction to the MRS cubes, we have discovered a new population of spatially extended hot dust in the β Pic debris disk emitting at 5 μm.
For dust at 10-20 au to emit at 5 μm, the dust grains are likely to be submicron to micron sized.In radiative equilibrium, the energy absorbed by a dust grain at a certain stellocentric distance is equal to the energy emitted from that dust grain, which can be written as where Q abs is the absorption coefficient at each wavelength for a specific material, J λ (T * ) is the mean intensity of the star, and B λ (T d ) is the Planck function representing emission from a dust grain of a certain temperature T d .In the small-grain limit, Q abs can be approximated as Q abs ∝ (1/λ), where λ is the wavelength of light.Using this small-grain approximation and solving for the dust grain temperature as a function of stellocentric distance as in Jura et al. (1998), we get where T d is the temperature of the dust, R * is the stellar radius, D is the stellocentric distance of the dust, and T * is the effective temperature of the star.In the blackbody limit for larger dust grains, where Q abs = 1 at all wavelengths, the equation for dust temperature as a function of distance has the same form, however, the exponent changes from 1/5 to 1/4.We calculated the temperature as a function of stellocentric distance from β Pic for dust grains in the blackbody and small-grain limits assuming stellar parameters of β Pic of T * = 8000 K (Lu et al. 2022) and R * = 1.73 R e (Kervella et al. 2004).The temperature as a function of distance for dust grains around β Pic in both the small-grain and blackbody approximations is shown in Figure 4.
To test whether the dust is better explained in the small-grain limit (∼300-400 K, see Figure 4) versus the blackbody approximation (∼100-160 K), we first summed the flux from all the pixels in the extended disk within the 5σ contour shown in Figure 3.This gives a total flux density from the extended emission at 5.2 μm of 0.03 ± 0.01 Jy.This is about 10 times less than the flux density of the spatially unresolved excess at the same wavelength shown in Figure 2. We then estimated the number of dust particles that would be required to produce an observed flux density of 0.03 Jy at 5.2 μm if the dust was in the blackbody approximation or the small-grain limit.At a stellocentric distance of 10 au, the dust temperature predicted by the blackbody approximation is 150 K.The observed flux from a single dust grain at a given temperature is given by the equation where B(T, λ) is the Planck function, a is the radius of the dust grain, and d is the distance of the dust to Earth.The number of dust particles required to produce the observed flux density from the extended emission at 5.2 μm is then where F obs = 0.03 Jy.For a dust temperature of 150 K in the blackbody limit and assuming a grain size equal to the wavelength of light of 5 μm, we found N dust = 1 × 10 34 .We then used Equation (3) to calculate the predicted flux density at 20 μm from this number of dust particles at 150 K.We found that the predicted 20 μm flux is 450 Jy, which is much greater than, and thus inconsistent with, the observed flux at 20 μm shown in Figure 1.We performed the same calculation, but with dust in the small-grain limit (a = 1 μm) and at a temperature of 350 K (corresponding to a distance of 10 au, as shown in Figure 4).We found that the number of dust particles required to produce the 0.03 Jy flux density was  (2013).The best-fit blackbody temperature is 500 ± 20 K.
N dust = 8 × 10 30 and that the predicted 20 μm flux from this number of dust grains at a temperature of 350 K was 0.2 Jy.This predicted flux density is less than the observed excess at 20 μm and thus consistent with the observations.Because the number of dust particles in the blackbody approximation required to produce the observed flux density from the spatially extended dust at 5 μm is inconsistent with the observed flux density at longer wavelengths, the spatially extended 5 μm excess is best explained by hotter submicron-sized grains in the small-grain limit.This calculation yields the same conclusion regardless of the assumed grain size and dust temperature for all temperatures in the small-grain and blackbody approximations within the red colored area in Figure 4.
We thus cannot explain the spatially extended thermal emission at 5 μm out to ∼20 au with blackbody grains, because at temperatures of 100-160 K, they do not emit significant flux at 5 μm.The spatially extended dust seen at 5 μm is better explained by dust grains in the small-grain approximation.In the small-grain limit (2πa = λ) at 5 μm, the grains would have submicron radii and could be below the blowout size for the system (see Section 4.2 for a calculation of the blowout size).

The Detection of [Ar II] Emission
We searched for emission from atomic and molecular gas in our MRS spectrum of the unresolved point source.We did not have any clear detections of molecular gas emission (we search for H 2 O, CO, CO 2 , and CH 4 ).We do, however, detect an emission line at 6.986 μm that we believe to be due to [Ar II] (rest wavelength of 6.9853 μm; Yamada et al. 1985).To verify this emission line further, we reduced each dither position independently and found that the emission line is present in all four dither positions.We also did not detect the emission line in N Car's spectrum, suggesting that the emission line is not due to a instrument systematic or pipeline artifact.We fit the spatially unresolved component of the [Ar II] line with a Gaussian profile to determine its line flux, line width, and radial velocity.The argon emission line with the best-fit Gaussian is shown in Figure 5.
From the best-fit Gaussian, we measured the FWHM of the line to be 86 ± 7 km s −1 and the line center to be at a velocity of 21.0 ± 2.4 km s −1 .At the end of channel 1C (7.65 μm), the spectral resolution of the MRS is ∼83 km s −1 (Argyriou et al. 2023), so the [Ar II] line is consistent with being spectrally unresolved.The measured barycentric radial velocity of β Pic is 20.0 ± 0.7 km s −1 (Gontcharov 2006), which is consistent with the center velocity of the [Ar II] line.We calculated a line flux of 2.4 × 10 −14 erg s −1 cm −2 .
We checked to see if the emission from [Ar II] is spatially resolved to determine if it is from circumstellar gas by subtracting a slice in the spectral cube outside of the argon line (6.9824μm) from the slice of the spectral cube at the peak of the [Ar II] line (6.9856μm).Doing this subtracted off the PSF of the unresolved point source as well as the continuum emission from the dust and planet (β Pic b) so only emission from [Ar II] gas remains.The continuum-subtracted cube slice   of the [Ar II] emission is shown in Figure 6.In the 6.9856 μm slice, there is a ∼5-10σ detection of spatially resolved [Ar II] emission, where σ is the standard deviation of the pixels in the annulus, at this wavelength slice, centered on the star with an inner radius of 1″ and an outer radius of 2 5.By visual inspection, the argon appears to have a similar spatial distribution as the dust (see Figure 6), indicating that the argon is indeed a part of the β Pic disk.To the best of our knowledge, this is the first detection of argon in the β Pic disk as well as in any debris disk.We used the line flux from the unresolved and resolved where M [Ar II] is the mass of [Ar II], d is the distance to the star, F is the line flux, h is Planck's constant, ν is the frequency of the light, A ul is the Einstein A coefficient, χ u is the fraction of atoms in the upper state, and m s is the mass of argon.We used an Einstein A coefficient of 5.3 × 10 −2 s −1 from Yamada et al. (1985).We calculated χ u assuming the argon is in local thermodynamic equilibrium (LTE) and that it has an excitation temperature equal to the radiative equilibrium temperature profile shown for blackbody grains in Figure 4.These assumptions of the excitation temperature and LTE make the estimated [Ar II] mass highly uncertain because, with only one [Ar II] line, we have no knowledge of the excitation temperature of the gas or if it is in LTE.
Since we expect the gas temperature to change as a function of disk radius, we calculated the line flux separately for the spatially resolved [Ar II] component.We did this by using rectangular extraction apertures of 2 pixels in height and the width was set by the 3σ contour of the [Ar II] emission shown in Figure 6.We placed these apertures centered at stellocentric distances of 10 and 15 au on both the northeast and southwest sides of the disk.We then fit the spectra from each aperture with a Gaussian and integrated the best-fit Gaussian to obtain a line flux as described above.We then summed together the line fluxes from the two sides of the disk at the same stellocentric distance, giving a total [Ar II] line flux at 10 and 15 au.We then used Equation (4) and excitation temperatures at 10 and 15 au equal to the radiative equilibrium temperatures shown in Figure 4 (150 and 125 K, respectively) to calculate the [Ar II] mass at each stellocentric distance.
For the unresolved [Ar II] emission, the location and LTE temperature of the gas is unknown because the spectral resolution of the MRS is not high enough to resolve gas kinematics.We then assumed that the unresolved [Ar II] has a temperature (180 K) corresponding to a stellocentric distance equal to the PSF FWHM at 6.986 μm, which is ∼7 au.This assumption is likely incorrect; however, we made it because it then gives an upper limit to the total gas mass in the unresolved point source.We then computed the total mass by summing the calculated mass at each stellocentric distance.This gives the total [Ar II] mass of 1 × 10 −3 M ⊕ .Given the uncertainty on the inner edge of the argon disk and the excitation of the gas, there is likely at least 1 to 3 orders of magnitude of uncertainty on this [Ar II] mass.

β Pictoris b
We also searched for β Pic b in our binned PSF-subtracted spectral cubes as it was predicted to be at an angular separation of 0 54 on 2023, January 11, the day of our MRS observations (Lacour et al. 2021;Wang et al. 2021).We detect a point source in channels 1A, 1B, and 1C in the PSF-subtracted cubes.Slices from the PSF-subtracted cubes from all three subbands of channel 1 showing the detection of a point source are shown in Figure 7.We detect this point source throughout channels 1A, 1B, and about half of channel 1C (up to ∼6.8 μm).At longer wavelengths, thermal emission from the edge-on disk becomes too bright and we cannot recover the point source from the emission of the cospatial disk (see Figure 8).At the three different wavelength slices shown in Figure 7, we fit a 2D Gaussian to the point source and measured its separation and position angle from the PSF center of the star in each subband.We computed the median of the position angle and angular separation of the point source for the three image slices from the three different subbands and took their standard deviations to be the uncertainty.We measured a separation of the point source from the star of 0 55 ± 0 02 and a position angle of 31°± 1°.The predicted location of β Pic b on the date of our observations, based on high-precision GRAVITY astrometry measurements, was a separation of 0 540 ± 0.″003 and a position angle of 31°.571 ± 0°.006 (Lacour et al. 2021), which are consistent with the location of the detected point source in our PSF-subtracted spectral cubes.
We computed contrast curves by masking out the planet and then calculating the standard deviation (1σ) in concentric annuli of a 2 pixel width in the binned PSF-subtracted image slices as a function of angular separation from the center of the PSF.We then divided these standard deviations as a function of angular separation by the spectrum of the unresolved point source to estimate the contrast performance.We did this both in the direction of the disk, to determine how well we recover β Pic b, and also outside of the disk in order to show the general contrast performance of our PSF subtraction with the MRS.To compute the contrast curves in the direction of the disk, we included all columns of pixels that included disk emission at the 3σ level.The contrast curves in the direction of the disk and outside of the disk for the three wavelength slices in Figure 7 are shown in Figure 9.The contrast curves computed in the direction of the disk show that we are achieving a 5σ contrast∼ 1 × 10 −4 at a separation of 1 0 at 5.3 and 5.7 μm.At 6.5 μm, we achieve a 5σ contrast of 3 × 10 −4 .The decrease in the contrast performance from 5.7 to 6.5 μm between separations of 0 6-2 5 seen in the left plot of Figure 9 is likely because the thermal emission from the disk is brighter and extended out to larger stellocentric distances at longer wavelengths.We also computed contrast curves with the extended disk emission masked out.This was done to estimate the general contrast performance of our PSF-subtraction technique with the MRS without contamination from the extended disk emission.These contrast curves are shown in the right panel of Figure 9.With the disk emission masked, our PSF subtraction and binning method achieves 5σ contrast levels ∼ 1 × 10 −4 at all three wavelengths at separations 1″.
We estimated the contrast of β Pic b by dividing its spectrum, extracted from aperture photometry in Section 3.7, by that of the central unresolved point source (star + disk) from Section 3.1.We then compared this contrast to the calculated contrast curves shown in Figure 9 to estimate the S/N of β Pic b.The S/N of β Pic b at 5.3, 5.7, and 6.5 μm in our binned MRS image cubes is 5.1, 6.3, and 5.5 respectively.

Cross-correlation
We also searched for β Pic b using cross-correlation methods with a template planetary atmosphere spectrum.This method involves cross-correlating the spectrum of each spaxel of the data cube with a planetary atmosphere model with various molecular absorption lines to search for molecules present in the atmosphere of the planet.This creates a cross-correlation map, where each spaxel in the IFU cube has one computed cross-correlation value.At the location of the planet in the spectral cube, where there is molecular absorption lines from its atmosphere, the cross-correlation function will peak, while the spaxels in the cube that do not contain a planet with the absorption features of the specific molecule being tested will appear as a featureless background.
For the atmosphere template spectrum, we generated petitRADTRANS (Mollière et al. 2019)  ).For each model we generated, we only included absorption from one molecular species to search specifically for each individual species in the atmosphere of  β Pic b.We searched for H 2 O, CO, CH 4 , CO 2 , and NH 3 in all wavelength channels.For each molecule, we input a mass fraction of 1 × 10 −4 into the atmosphere model.We tried a range of mass fractions from 1 × 10 −3 to 1 × 10 −6 and found that the result of the cross-correlation detections did not change for any molecule.
We first removed the continuum from the atmosphere model by fitting the model with a B-spline of 10 breakpoints and subtracting off this continuum fit, leaving only the molecular absorption lines.The same was done for the spectrum of each spaxel in the spectral cube, to remove the continuum spectrum.Once we removed the continuum of the atmosphere model and each spaxel, we then cross-correlated the atmosphere model templates with the continuum-subtracted spectra of each spaxel to calculate a cross-correlation value.We then estimated the S/N per spaxel of any potential cross-correlation detection by dividing by the standard deviation of cross-correlation values of all the spaxels, which is similar to what is done by Mâlin et al. (2023).
This produces a S/N detection map where each spaxel in the image cube has a S/N value.Figure 10 shows the crosscorrelation map of the full channel 1 cube with a petitRADTRANS atmosphere model containing H 2 O as the only absorbing molecular species.We find a bright pixel in the cross-correlation map at the same location as the brightest pixel of the point source recovered in the PSF-subtracted images slices.Out of all the molecules and MRS channels we tested, H 2 O in channel 1 was the only one that resulted in at least one pixel above a S/N value of 5.The bright pixel in our detection map has a S/N value of 5. Previous studies that use cross-correlation techniques with near-infrared ground-based data have detections of planets that cover multiple IFU spaxels (Hoeijmakers et al. 2018), while here we only detect one spaxel with a significant crosscorrelation value, which lines up with the position of the planet from the PSF-subtracted image slices.Doing the same crosscorrelation method with all the molecular species listed above did not result in any spaxel in the cubes having a cross-correlation value greater than five in any of the MRS channels (1-4), so finding a spaxel with a S/N value of 5.7 using H 2 O does not appear to be random chance.Furthermore, with 1498 spaxels in the channel 1 data cube, it is unlikely that the crosscorrelation method would have randomly selected the same spaxel as the point source in the PSF-subtracted cubes to have a significant cross-correlation value.
The detection of a point source through PSF subtraction at the predicted location of β Pic b based on orbit fits from Lacour et al. (2021) in three MRS subbands as well as the tentative detection of water through cross-correlation at the same location as the point source provides enough confidence for us to conclude that we are in fact detecting β Pic b in the MRS data.

Spectrum of β Pictoris b
We performed aperture photometry on β Pic b in each slice of the binned PSF-subtracted spectral cubes to extract a spectrum of the planet.We found that the colocated spatially extended dust emission contaminates the spectrum of the planet.In order to figure out the best strategy for subtracting the contribution of the dust, we performed an injection-andrecovery test using the N Car PSF.We first scaled the PSF of N Car in each wavelength slice such that the brightest pixel of the N Car PSF matches the brightest pixel of β Pic b in the first wavelength slice in each subband.In the last wavelength slice of each subband, the brightest pixel in the injected point source was within 5% of that of β Pic b.The spectrum of N Car was not changed aside from the scaling.We then injected the scaled N Car PSF on the opposite side of the disk at the same separation as β Pic b before performing PSF subtraction.We PSF subtracted each spectral cube with N Car injected the same way as described in Section 2.3 and produced a binned and PSF-subtracted spectral cube, which is shown in Figure 11.
We performed aperture photometry on the injected planet and subtracted off the disk background by placing a background aperture of the exact same size at a further separation of 0 8 on the same side of the disk as the injected planet.This location of the background aperture was chosen because it is the closest location to the injected planet without the extraction and background apertures overlapping.We used Figure 9. Left: contrast curves in the direction of the disk from the binned PSF-subtracted image cubes at wavelength slices in channels 1A, 1B, and 1C (same slices shown in Figure 7).The red point shows the contrast of β Pic b calculated at 5.7 μm from our data.We do not show the contrast of β Pic b at the other two wavelengths for clarity.Right: contrast curves from the binned PSF-subtracted image cubes at the same wavelength slices as in the left plot but over the entire field of view with the extended emission from the disk masked.This shows the contrast performance of our PSF subtraction with the MRS without contamination from the extended disk emission.
an aperture radius of 0.5 times the FWHM of the PSF because we found that using a larger aperture includes more thermal emission from the disk, which results in a worse recovery of the injected spectrum.We created our own aperture correction for this aperture size using N Car.We binned the spectral cube of N Car the same way as for the PSF-subtracted cubes, and then extracted the spectrum of N Car the same way as described above for the unresolved point source of β Pic.We also extracted the binned spectrum of N Car using an aperture of 0.5 times the FWHM and divided the two to derive an aperture correction for an aperture size of 0.5 times the FWHM.Using this background-subtraction method and extraction aperture size, we were able to recover the general shape of the spectrum of the injected point source; however, we were unable to recover the absolute flux density values.The injected and recovered spectra for channel 1B are shown in Figure 12.At the same wavelength that we were unable to detect β Pic b (∼7 μm), we were unable to recover the injected point source as well (see Figure 11).This is potentially because beyond this wavelength, the thermal emission from the disk is bright enough that it dominates over the planet flux, and we cannot recover the planet at longer wavelengths.We calculated the difference between the injected and recovered spectra and then subtracted this off of the extracted spectrum of β Pic b.This was to remove the excess flux from the cospatial disk that was not subtracted with the background aperture, as indicated by the injection-and-recovery test.
We used the same background aperture, same aperture size, and same aperture corrections for extracting the spectrum of β Pic b as we did for the injected fake planet.We placed the extraction aperture on the center of β Pic b measured by the 2D Gaussian fit for each subband.We estimated the uncertainty on the spectrum of β Pic b from the contrast curves shown in Figure 9.That is, the standard deviation of the annulus containing the separation of β Pic b in each wavelength slice is used as the uncertainty.The final spectrum of β Pic b, along with an ExoREM (Charnay et al. 2018;Blain et al. 2021) planetary atmosphere model, is shown in Figure 13.

Atmospheric Model Fitting
We used the Species package (Stolker et al. 2020) along with three different sets of model grids of synthetic spectra to fit the spectra of β Pic b.The first was the ATMO model grid (Tremblin et al. 2015) presented in Petrus et al. (2023).The second model grid we used was the DRIFT-PHOENIX model grid (Helling et al. 2008) and the third was the ExoREM grid (Charnay et al. 2018;Blain et al. 2021).The model parameters and the prior ranges for all three model grids are shown in Table 1.We used uniform priors for each parameter in each of the model grids.In the modeling fitting, we included our MRS spectrum, the GRAVITY K-band spectrum from GRAVITY Collaboration et al. (2020), and the GPI YJH spectra from Chilcote et al. (2017).We also included photometric measurements from Magellan and Gemini NICI (Males et al. 2014) and VLT/NACO (Bonnefoy et al. 2013).Since we were unable to recover the absolute flux density of the injected planet in the injection-and-recovery test because of the cospatial thermal dust emission, we left the overall flux scaling of the MRS spectrum as a free parameter in the model fitting, similar to what was done in Kammerer et al. (2021) to align spectra from GPI and GRAVITY for the substellar object HD 206893 B. We weighted each spectrum and photometric point in the log-likelihood function in the Species model fitting such that each data set, spectral and photometric, had the same weighting.This was to ensure that the model fitting is not dominated by the GRAVITY spectrum, which contains the most data points.We then infered the posterior distribution of the parameters using nested sampling with PyMultiNest (Feroz et al. 2009;Buchner et al. 2014;Feroz et al. 2019).The MRS spectrum of β Pic b is shown in Figure 13 with the best-fit ExoREM model.The posterior distributions for the model fits are shown in the Appendix.
Figure 14 shows the best-fit DRIFT-PHOENIX, ATMO, and ExoREM models to all of the spectra and photometry points.Table 1 shows the best-fit model parameters for each of the model grids.Both the ATMO and DRIFT-PHOENIX model grids give consistent measurements within the uncertainties for log(g) and effective temperature.The ExoREM grid gives a lower effective temperature by ∼200 K than the other two model grids.The temperatures from the ATMO and DRIFT-PHOENIX grids are within the uncertianties of each other.Fits with ExoREM models have given lower effective temperatures than other model grids in other fits in the literature as well (e.g., GRAVITY Collaboration et al. 2020;Kammerer et al. 2021).All three model grids produce inconsistent metallicities.The ATMO model grid gives a best-fit metallicity that is subsolar, while the DRIFT-PHOENIX model grid gives a supersolar metallicity and the ExoREM grid gives a larger supersolar metallicity (see Table 1).The uncertainty in the metallicity from the ATMO grid does include solar metallicity, although the uncertainties of the metallicity between the three model grids do not overlap, suggesting that we do not constrain the metallicity of the planet well with this method of grid model fitting.
To test how the new MRS data affect the constraints on the atmosphere of β Pic b from our modeling, we repeated the grid model fitting excluding the MRS spectrum.We find that the addition of the MRS spectrum does not significantly change the results from the grid model fitting.All of the atmospheric parameters from fitting each of the three model grids are consistent with each other (within the uncertainties) when we performed the model fits with and without the MRS spectrum.The C/O ratio from the ExoREM and ATMO model fits is similarly not significantly changed.Excluding the MRS spectrum with the ExoREM grid yielded a C/O ratio of -+ 0.38 0.02 0.10 , which is consistent with the C/O ratio obtained when including the MRS spectrum (See Table 1).Excluding the MRS spectrum from the ATMO model grid fit gives a C/O ratio of -+ 0.38 0.06 0.10 , which is consistent with the ATMO fit that included the MRS spectrum.The shape of the MRS spectrum of β Pic b is therefore consistent with what is predicted by the atmosphere model fits to the near-infrared spectra and photometry.

β Pictoris b's Atmosphere Composition
The composition of a giant planet's atmosphere, usually the C/O ratio, can potentially be used as a tracer to infer formation mechanisms and location in the parent protoplanetary disk (Öberg et al. 2011).There are two main hypothesis for forming a giant planet.The first is through gravitational collapse, which is similar to star formation (Bodenheimer 1974).In this scenario, a region of the protoplanetary disk becomes gravitationally unstable and quickly collapses to form a planet that then cools over time.The second is core accretion, where a solid core is formed that slowly accretes gas from the surrounding disk until the mass of the accreted gas is similar to that of the core.Then, a phase of runaway gas accretion occurs and the planet accretes a significant amount of gas over Figure 12.Channel 1B spectrum of the injected point source (orange) and the spectrum we recover from the injected point source (blue).We generally recover the overall slope of the injected spectrum, however, we are unable to recover the absolute flux density of the spectrum.Even with subtracting off the flux density in a background aperture in an attempt to remove contamination from the spatially extended disk, the recovered spectrum still has a greater flux density than the injected spectrum by a factor ∼ 1.3-1.4.
Figure 13.MIRI MRS spectrum of β Pic b (black) with the best-fit ExoREM model overplotted (blue) and a best-fit blackbody (orange).The broad absorption bands from 5.1 to 6.0 μm and from 6.4 to 6.8 μm are due to absorption from water vapor.The local maximum in the spectrum at 6.3 μm is due to a relative absence of water vapor lines near this wavelength.The spectrum is better fit by the atmosphere model that contains water absorption rather than the blackbody with no molecular features since the blackbody cannot reproduce the shape of the spectrum between 6 and 6.5 μm.The spectrum showing signs of water absorption is consistent with the tentative detection of water from the cross-correlation technique.a short period (Lissauer & Stevenson 2007).As discussed by GRAVITY Collaboration et al. (2020), in the gravitational collapse scenario, the C/O ratio of the planet is expected to match that of the host star because the formation happens quickly and all solid and gaseous material in the disk that collapses into the planetary atmosphere have a combined stellar C/O ratio.The core accretion scenario takes longer than the gravitational collapse scenario, which provides more time for the accretion of solid materials, including planetesimals.In the core accretion scenario, without enrichment from solids before Middle: same as the top panel but now with the best-fit model from the ATMO grid.Bottom: same as the top panel but with the ExoREM grid.The gray spectral models are sampled from the posterior distribution.The residuals from the best-fit model are shown in the bottom panel of all plots with the two dotted lines representing ±5σ residuals.All three models appear to be well fit, but the ExoREM grid yields a systematically lower temperature than the other two.
Note.The c red 2 value reported here is for all the photometry and spectra shown in Figure 14, not just the MRS spectra.
the runaway gas accretion phase, the atmosphere of the planet is not comprised of a combination of gas and solids like in the gravitational collapse scenario, but is only composed of the gas.Therefore, without the enrichment of solids, the C/O ratio of the atmosphere in this scenario is expected to match the C/O ratio of the gas in the disk, which is superstellar.However, with enrichment from solids before the runaway gas accretion stage, the C/O ratio of the atmosphere can be lowered to substellar values as there is more time to accrete solid material in the core accretion scenario than the gravitational collapse scenario (GRAVITY Collaboration et al. 2020).
The wavelength range of our MRS spectrum of β Pic b contains a water absorption band, which, along with the CO band in the GRAVITY spectrum of β Pic b, could potentially help constrain the C/O ratio of the planet's atmosphere.From our modeling with the ATMO grid of all the spectra and photometry in Figure 14 Using the C/O ratio of a planetary atmosphere to infer its formation history requires knowledge of the C/O ratio of the host star; however, there is not a published value for the C/O ratio of β Pic in the literature.Because of this, we compared to the C/O ratio measured from the β Pic moving group member HD 181327 as in Reggiani et al. (2024).The C/O ratio of HD 181327 was used as a proxy for β Pic since they likely formed out of the same molecular cloud.The C/O ratio of HD 181327 is 0.62 ± 0.08 (Reggiani et al. 2024).The C/O ratio for β Pic b we got from the ATMO grid fitting is substellar and does not contain the stellar C/O in its uncertainty.Similarly, the best-fit C/O ratio from the ExoREM grid does not contain the stellar C/O.Both C/O ratios from the two model grids are within the uncertainties of each other and suggest a substellar C/O ratio for the planet.
The two C/O ratios we inferred from the grid model fitting are both consistent within the uncertainties with the C/O ratio measured from a free retrieval of the GRAVITY and GPI spectra of

Dust Accretion onto β Pictoris b and c
The infrared excess seen at 5 μm suggests that there is dust produced in the inner few astronomical units of the system.The dust emitting at 5 μm that we see to be spatially extended out to ∼20 au in the PSF-subtracted image slices is best explained by grains in the small-grain limit (2πa = λ).This spatially extended dust seen at 5 μm is then likely to be below the blowout size for the system (see Figure 15 for the blowout size) and can be driven outwards over time by radiation pressure.The blackbody stellocentric distance for the 500 K dust (from our best-fit blackbody to the 5-7.5 μm excess) is 0.9 au.If small dust grains are being produced at 0.9 au, and are then radiatively driven outwards, some will likely collide with the planets β Pic b and β Pic c, which have semimajor axes of 9.9 ± 0.05 au and 2.72 ± 0.02 au (Nowak et al. 2020), respectively.We used our MIRI MRS data to perform an order-of-magnitude estimate for the dust accretion rate from the inner hot dust onto the two known planets in the β Pic system.
To estimate the dust accretion rate onto β Pic b and β Pic c, we first calculated the dust mass from the 5 μm unresolved excess using the following equation from Lisse et al. (2009) where F λ is the flux density from the dust at a given wavelength, D is the distance to the star (19.6 pc; Gaia Collaboration et al. 2023), B λ (T) is the blackbody function at a given temperature T, Q abs (λ) is the absorption coefficient of the material at a given wavelength, S is a scaling factor for the particle size distribution, a is the grain radius, and dn da is the particle size distribution.We assumed the particle size distribution that is expected for collisional equilibrium where µ - dn da a 3.5 (Dohnanyi 1969).We assumed a dust temperature of 500 K based on the fit to the 5-7 μm excess and for the absorption coefficient, we calculated Q abs (λ) assuming Mie theory and the optical constants of amorphous olivine (Mg 2 SiO 4 ) from Jäger et al. (2003) for the different dust grain sizes spanning the region = a 0.03 min μm to = a 10 max μm.Because we do not know the composition of the dust grains producing the 5 μm excess, we assume they are amorphous silicate grains, since these were detected by Spitzer around β Pic (Lu et al. 2022).Because of this uncertainty in composition, we also repeat this calculation assuming optical The black dashed horizontal line shows where β = 0.5.Grains sizes with β > 0.5 become unbound and are blown out of the system.The vertical black dashed dotted line shows a = λ/(2π) for λ = 5 μm, which is a = 0.8 μm.The small-grain limit is defined as 2πa = λ.
At 5 μm, the flux density from the excess dust emission is 0.2 Jy (see Figure 2).We determined a dust mass estimate by scaling the particle size distribution (using the scaling factor S) such that the dust flux density calculated on the right-hand side of Equation ( 5) is equal to the observed dust flux density at 5 μm.We then calculated the dust mass by integrating over the scaled particle size distribution, using the equation where M is the dust mass, S is the scaling factor from Equation (5), ρ is the bulk density, and a is the grain radius.The p a 4 3 3 term is the volume of a spherical grain and a −3.5 is from the assumed particle size distribution.We assumed a bulk density of ρ = 3.3 g cm −3 .Integrating Equation (6) gives a dust mass from the 5 μm excess of 7 × 10 21 g, assuming that all of the excess is produced from olivine.Doing the same calculation for pyroxene grains gives a dust mass of 1 × 10 22 g.
We next calculated the blowout size for olivine and pyroxene silicate grains around β Pic.The ratio of the radiation force to the gravitational force can be written as where L * is the stellar luminosity, M * is the mass of the star, c is the speed of light, G is the gravitational constant, a is the radius of the dust grains, ρ is the bulk density of the material, and 〈Q pr 〉 is the average radiation pressure coupling coefficient for a given material.We used a mass of β Pic of 1.75 M e (Crifo et al. 1997).〈Q pr 〉 is given as pr pr where F λ is the flux density at a given wavelength of the star.
For F λ , we used the photosphere model from Lu et al. (2022) and we calculated Q pr assuming Mie theory and amorphous olivine and pyroxene grains.Dust grains with β-values >0.5 become unbound and are blown outwards, setting the blowout size for the system (Krivov et al. 2006).Figure 15 shows the calculated β-values as a function of grain size along with a horizontal dashed black line indicating where β = 0.5.Grain sizes with β-values above the horizontal black line are blown out of the system, which sets the blowout radius for the system to be between 0.03 and 1.14 μm, assuming olivine grains.For pyroxene grains, the blowout radius is between 0.03 and 1.24 μm.
Using the β-values of the dust grains, we calculated their terminal velocity and the time it takes for the dust grains to reach the planets assuming that their initial location is at the blackbody distance of the 5 μm excess (0.9 au, since we do not know the physical location of the grains producing the unresolved excess).We also assumed that the dust has a distribution that is azimuthally symmetric, and that the dust grains travel at their terminal velocity until they accrete onto the planets.We calculated the terminal radial velocity of the dust grains using this expression from Su et al. (2005) as where v r is the terminal radial velocity and r init is the stellocentric distance where the dust grains are released (we assumed the minimum blackbody distance of r init = 0.9 au).
To calculate the dust accretion rate onto β Pic b and β Pic c, we binned the scaled particle size distribution we determined above into radius bins of length 0.1 μm from 0.03 to 1.14 μm (the blowout size range for olivine) and calculated a dust mass in each radius bin by integrating Equation (6) over each radius bin.We then calculated the terminal velocity of the dust grains in each radius bin by calculating β for each bin center radius using Equation (7).Assuming that the dust grains start at 0.9 au and using the semimajor axes of β Pic b and β Pic c of 9.9 and 2.72 au, respectively, we then calculated the time it takes for the dust to travel from their initial locations to the distance of the planets assuming that they travel the entire distance at the terminal velocity.We checked this assumption by integrating the equation of motion of a dust grain in orbit around β Pic in two dimensions, assuming that the dust grain starts out on a circular orbit with Keplerian velocity at 0.9 au.The force in the radial direction for a dust grain under gravitational and radiation forces, as shown in Lisse et al. (1998) and Krivov et al. (2006), is where r is the distance from the dust grain to the central star.
By numerically integrating Equation (10) in two dimensions with a β-value of β = 0.5, we found that the dust grain will reach its terminal velocity in about a tenth of the time it takes to reach the semimajor axis of β Pic b and 40% of the time it takes to get to β Pic c.For dust grains of larger β-values, they reach their terminal velocity quicker.Thus, the assumption that the dust grains travel the entire distance at the terminal velocity likely does not have a large effect on the estimated dust accretion rate.
We estimated the dust accretion rate in each radius bin of the particle size distribution using the equation where  M is the accretion rate, M d is the dust mass in each size bin, R p is the planet radius, t is the time it takes the dust in each radius bin to reach the planet from its origin location, d is the distance between the dust origin radius and the planet semimajor axis, and Ω is the solid angle subtended by the disk.We used planetary radii of 1.36 R J for β Pic b (GRAVITY Collaboration et al. 2020) and 1.2 R J for β Pic c (Nowak et al. 2020).We summed this equation across all radius bins in the particle size distribution to get a total dust accretion rate.
To estimate Ω, we first measured the vertical profile of the spatially resolved dust at 5 μm by collapsing the disk along the midplane to create a single integrated vertical profile of the disk with a high S/N.Then, we fit the vertical profile with a Lorentzian to measure the FWHM of the disk.The disk vertical structure and the best-fit Lorenztian is shown in Figure 16.To measure the FWHM of the disk, we first subtracted in quadrature the FWHM of the PSF from the best-fit FWHM of the Lorentzian.We used the measured FWHM of the PSF described above.The FWHM of the disk measured from the best-fit Lorentzian after deconvolving with the PSF FWHM is 11.5 ± 0.3 au.
We calculated the solid angle Ω using the vertical profile of the 5 μm spatially extended disk determined above.We calculated the opening angle of the 5 μm disk by using half the FWHM of the vertical profile (11.5 au) and the outer radius of the disk (20 au) and then solving ( ) q = h r tan , where h is 5.75 au and r is 20 au.This gives an opening angle of θ = 16°.We then integrated the differential for solid angle dΩ = ( ) q q f d d sin over the opening angle and over 2π radians azimuthally.This gives a solid angle of Ω = 1.92 sr.Then, using Equation (11), we calculated a dust accretion rate for β Pic b and β Pic c of  = ´- M 2 10 17 M J yr −1 and  = ´-M 2 10 15 M J yr −1 , respectively, assuming olivine dust grains.Repeating the calculation with pyroxene dust grains gives dust accretion rates on β Pic b and β Pic c of  = ´-M 3 10 17 M J yr −1 and  = ´- M 3 10 15 M J yr −1 , respectively.
We also performed this calculation for an Fe/Mg ratio of 1:1 for olivine (MgFeSiO 4 ) and pyroxene (Mg 0.5 Fe 0.5 SiO 3 ), and we found that the β-values increase for these more iron-rich grains compared to the pure magnesium grains.For pyroxene, using a composition of Mg 0.5 Fe 0.5 SiO 3 , the β-values increased by about a factor of 1.4 from those shown in Figure 15.The blowout size for Mg 0.5 Fe 0.5 SiO 3 is a 1.8 μm.For this composition of 50% Fe and Mg, the value for β does not come below 0.5 at the small particle size like it does for MgSiO 3 and Mg 2 SiO 4 , as Figure 15 shows.Instead, all particles below 1.8 μm are blown outwards, so we have to estimate the smallest size grains in the particle size distribution to get a limit on the smallest particle sizes that are present in the disk.Krijt & Kama (2014) estimated the smallest particle size in a debris disk and relate it to the blowout size with the equation where a is the semimajor axis of the dust, L * is the luminosity of the star, f is the ratio of the relative velocity of the colliding bodies and their Keplerian velocity, η is the fraction of the precollision kinetic energy of the colliding bodies that is converted into making a new surface, and γ is the surface energy per unit surface area of the material.As done in Krijt & Kama (2014), we assumed that the collisions are below the hypervelocity regime and f = 10 −2 and η = 10 −2 .We used γ = 0.05 J m −2 for silicate grains (Krijt & Kama 2014) and we assumed the dust grains are at the blackbody distance inferred from the blackbody fit above of 0.9 au.We find that for silicate grains, the minimum dust particle radius in the disk is 0.028 μm.We then used this size as the minimum particle size for the more iron-rich grains (MgFeSiO 4 and Mg 0.5 Fe 0.5 SiO 3 ) and calculated the dust accretion rate using the same method as above.We estimated a dust accretion rate onto β Pic b for Mg 0.5 Fe 0.5 SiO 3 of 6 × 10 −17 M J yr −1 and 5 × 10 −17 M J yr −1 for MgFeSiO 4 .For β Pic c we obtained dust accretion rates of 6 × 10 −15 M J yr −1 and 5 × 10 −15 M J yr −1 for Mg 0.5 Fe 0.5 SiO 3 and MgFeSiO 4 respectively.In this calculation, we did not assume that the grains have any porosity; however, Arnold et al. (2019) showed that for silicate grains with 97.5% porosity around β Pic, the blowout size increases to ∼10 μm, meaning that more dust particles are blown outwards, increasing the dust accretion rates onto the planets.
Given that the masses of these planets are 9.0 ± 1.6 M J for β Pic b and 8.2 ± 0.8 M J for β Pic c (Nowak et al. 2020), these dust accretion rates do not contribute significantly to the total masses of these planets, even over the ~23 Myr age of the system.If, however, the small dust grains that are accreted can remain aloft in the planet's atmosphere, they could potentially impact the observed spectra of the planet at near-infrared and mid-infrared wavelengths.Based on near-infrared photometry of β Pic b, Bonnefoy et al. (2013) found that the atmosphere of β Pic b is likely dusty, which suggests that there could be an absorption feature in the spectrum of β Pic b due to silicates at 10 μm.However, since the thermal emission of the disk at 10 μm dominates over the flux of the planet in our MRS data (see Figure 8), we were unable to search for silicate absorption in the atmosphere of β Pic b.Based on these dust accretion rate estimates, we expect the spectrum of β Pic c to be more affected by dust than β Pic b, although we cannot test this here because we cannot recover β Pic c with the MRS since it is too close to the host star.

[Ar II]
The detection of circumstellar [Ar II] emission around β Pic with the MRS is surprising, so here, we speculate on possible explanations for this detection of argon.Similar to the dust in the disk, the gas around β Pic is also subject to blowout from radiation pressure (Beust et al. 1989).The radiation pressure for the atomic gas around β Pic depends on the number and strength of ground state transitions for a given species (see Fernández et al. 2006).Fernández et al. (2006) found that singly ionized argon experiences zero radiation force and is not subject to radiative blowout because it does not have any ground state transitions in the 0.1-5 μm range.It is then possible for Ar II to accumulate over time in the β Pic disk if it is continuously produced through the destruction of minor bodies in the system.Other atomic species detected around β Pic are C, O, Na, Mg, Al, Si, S, Ca, Cr, Mn, Fe, and Ni (Brandeker et al. 2004;Roberge et al. 2006), which are thought to be of secondary origin (Fernández et al. 2006), meaning not from the protoplanetary disk, but rather created from the destruction of minor bodies.Similarly, the molecular CO in the β Pic disk is thought to be of secondary origin (Matrà et al. 2017).
If the argon is also produced from the destruction of minor bodies, then it could persist in the disk since it does not feel radiation pressure.If we assume that the argon production rate is constant over the age of the system (∼20 Myr), we estimate an argon production rate of 3 × 10 14 kg of argon per year to produce the total argon mass seen with the MRS.For reference, this is roughly equivalent to producing the mass of Halley's Comet (2.2 × 10 14 kg; Hughes 1985) of Ar II each year in the β Pic system.This estimate is highly uncertain given the uncertainties in the argon mass and the likely incorrect assumption that the gas can persist over millions of years in the disk since it does not feel radiation pressure.
In the solar system, minor bodies such as comets and meteorites have been found to contain argon.The argon found in solar system meteorites is thought to be primordial presolar gas that has become "trapped" within the matrices of meteorites (Huss et al. 1996;Patzer & Schultz 2002;Vogel et al. 2003).Argon has also been detected in the coma of comets 67P/ Churyumov-Gerasimenko (Balsiger et al. 2015) and Hale-Bopp (Stern et al. 2000).Argon sublimates at a temperature of 40 K, so the detection of argon in comet 67P/Churyumov-Gerasimenko was interpreted to mean that the comet formed in a cold outer region of the protosolar nebula (Balsiger et al. 2015).If the minor bodies in the β Pic system also contain trapped primordial argon like those in the solar system, it is possible that the argon detected around β Pic is from the destruction of minor bodies.
Although the production of argon in the β Pic system could potentially be explained by the destruction of minor bodies, the mechanism of ionization is unclear.Using a photoionization code, Fernández et al. (2006) found that argon should remain mostly neutral throughout the β Pic disk with an ionization fraction on the order of 10 −6 .Their calculation considered UV photons from the star, the interstellar radiation field, and ionization from cosmic rays, and they found that the UV radiation field dominated the ionization.They did not consider, however, coronal X-ray emission from β Pic, which was detected by Hempel et al. (2005) and Günther et al. (2012).The ionization fraction of argon predicted by Fernández et al. (2006) would require an unphysically large amount (thousands of Earth masses) of neutral argon in the disk to produce the detected [Ar II] emission seen with the MRS.The ionization of neutral argon could then potentially be due to X-ray emission from the corona of β Pic.However, confirmation of this would require a photoionization model that considers coronal emission from the star, which is beyond the scope of this paper.

Conclusion
As a part of JWST GTO program 1294, we present MIRI MRS observations of the β Pic exoplanetary system.Our main findings are: 1. We detect an infrared excess from the unresolved point source from 5 to 7 μm that is best fit by a 500 K blackbody, indicating the presence hot dust in the inner few astronomical units of the system.
2. Through PSF subtraction, we detect a spatially resolved dust population emitting at 5 μm extending out to stellocentric distances ∼ 20 au.This dust population is best explained by small subblowout-sized grains because larger blackbody grains in radiative equilibrium are not hot enough between 10 and 20 au to emit significant radiation at 5 μm. 3. We detect and spatially resolve circumstellar [Ar II]  emission for the first time around β Pic.This is also the first argon detection in a debris disk.4. The unresolved and resolved hot dust population suggests that dust is produced in the inner few astronomical units of the system where the subblowout-sized grains are radiatively driven outward where they could accrete onto β Pic b and β Pic c.We use our MRS data to estimate a dust accretion rate onto β Pic b and β Pic c of  = M -10 17 M J yr −1 and  = - M 10 15 M J yr −1 , respectively.5. We detect β Pic b with both PSF subtraction and crosscorrelation with atmosphere models.We find that water is the only molecule we detect in cross-correlation.We also present the first mid-infrared spectrum of the planet from ∼5 to 7 μm, which includes a water absorption band.6. Grid model fitting for β Pic b reveals that the metallicity of the planet is entirely model dependent and that the C/O ratio is likely between 0.3 and 0.5.The atmospheric modeling with both ExoREM and ATMO mode grids appears to favor a substellar C/O ratio for β Pic b.

Figure 1 .
Figure1.Left: MRS spectrum of the β Pic unresolved point source of channels 1-4.The feature from ∼8 to 12 μm is from emission from silicate grains.Right: channel 1 MRS spectrum of the β Pic central point compared with a stellar photosphere model of the β Pic star fromLu et al. (2022) that was best fit to a near-infrared spectrum.This shows an infrared excess from 5 to 7 μm which can be explained by hot dust in the system.The vertical dashed gray line shows the detection of [Ar II] emission at 6.986 μm.The other spikes in the spectrum do not appear to be real emission lines.This is because they only cover one wavelength bin, unlike the [Ar II] line which covers five.Furthermore, these spikes also only appear after applying the relative spectral response function (RSRF) and are well aligned to noise features seen in the spectrum of N Car.They then likely result from noise that is injected by performing the RSRF with N Car.The [Ar II] line, however, is present both before and after applying the RSRF.

Figure 3 .
Figure 3. Binned 5.2 μm slice of the PSF-subtracted spectral cube showing the extended emission from hot dust.We apply an image normalization where =v 200 min

Figure 4 .
Figure 4. Dust temperature as a function of stellocentric distance from β Pic for dust grains in the small-grain approximation (black solid line) and for dust grains in the blackbody approximation (black dashed line).The red colored area shows the approximate location of the spatially resolved dust seen at 5 μm (based on where the flux drops below a 5σ detection).

Figure 5 .
Figure 5. [Ar II] emission line at 6.986 μm from the β Pic unresolved point source with the best-fit Gaussian overlaid.The best-fit blackbody, as shown in Figure 2, was subtracted from the spectrum in this figure.
[Ar II] emission to estimate the total mass of [Ar II].Under the assumption that the emission is optically thin, the total mass of [Ar II] is related to the line flux by the equation

Figure 6 .
Figure 6.Left: [Ar II] image created by subtracting the 6.9824 μm slice (outside the argon line) from the 6.9856 μm slice (the peak of the argon line).The red contours shown are 3, 5, and 10σ contours.The black circle shows the FWHM of the JWST PSF at this wavelength, indicating that the [Ar II] emission is spatially resolved out to ∼20 au.The red cross shows the position of the central star.The north-east orientation of the images are the same as in Figure 3. Right: binned PSF-subtracted cube slice at 6.5 μm with 3, 5, and 10σ contours of the continuum-subtracted [Ar II] emission overlaid.The point source on the northeast side is β Pic b.
models at their full spectral resolution with absorption from various molecules.We used the physical parameters of β Pic b determined by GRAVITY Collaboration et al. (2020; T eff = 1742 K, ( ) = log g 4.34

Figure 7 .
Figure 7. Binned PSF-subtracted image slices showing the detection of a point source at the predicted location of β Pic b (separation = 0 54, PA = 31°.57; Lacour et al. 2021; Wang et al. 2021).The left panel shows a binned slice from channel 1A, the middle panel shows one from 1B, and the right panel a binned slice from 1C.The red star indicates the position of the central star determined from the 2D Gaussian fit to the PSF.The pixels within a 3 pixel radius of that center position are masked for clarity.The S/N of the brightest pixel in the point source, going from left to right, is 5.1, 6.3, and 5.5.

Figure 8 .
Figure 8. PSF-subtracted and binned slice of the channel 2B data cube at 10 μm.The red star indicates the position of the central star and the red cross indicates the position of β Pic b at the position angle and separation we found from the channel 1 data cubes.The pixels within a 2 pixel radius of the PSF center have been masked for clarity.The image slice is dominated by thermal emission from the disk and we do not detect a clear point source on one side of the disk as we did in channel 1.The orientation of this image relative to north is the same as those shown in Figure 7.
7. Mâlin et al. (2023) simulated cross-correlation techniques using simulated MRS observations of the β Pic system and similarly found that H 2 O was the only molecule that yielded a significant S/N value for β Pic b (see Figure C.2 of Mâlin et al. 2023).

Figure 10 .
Figure10.Cross-correlation map of the channel 1 IFU cube with an atmosphere template spectrum with water being the only molecular absorber.We detect a bright pixel with a cross-correlation S/N value of 5.7 at the same pixel location as the brightest pixel of the point source shown in Figure7.The red cross shows the position of the central star.

Figure 11 .
Figure 11.Left: same as the 5.7 μm slice shown in Figure 7 but now with N Car injected as a fake planet on the opposite side of the disk.The point source on the northeast side of the disk is β Pic b.Right: same as the left panel but for the 7 μm slice, where there is no clear detection of the injected point source or of β Pic b.

Figure 14 .
Figure14.Top: best-fit DRIFT-PHOENIX model (black) with spectra and photometry of β Pic b.Middle: same as the top panel but now with the best-fit model from the ATMO grid.Bottom: same as the top panel but with the ExoREM grid.The gray spectral models are sampled from the posterior distribution.The residuals from the best-fit model are shown in the bottom panel of all plots with the two dotted lines representing ±5σ residuals.All three models appear to be well fit, but the ExoREM grid yields a systematically lower temperature than the other two.
, we get a C/O ratio for β Pic b of -+ 0.39 0.06 0.10 .This is consistent with the C/O ratio determined by GRAVITY Collaboration et al. (2020) from a petitRAD-TRANS retrieval of the GRAVITY and GPI spectra of β Pic b.This retrieval yielded a C/O ratio of -+ 0.43 0.03 0.04 .The ExoREM model grid gives a C/O ratio of -+ 0.36 0.05 0.13 , which also contains the C/O ratio from GRAVITY Collaboration et al. (2020) within the uncertainty.
(GRAVITY Collaboration et al. 2020).This substellar C/O ratio was argued to favor the core accretion scenario with planetesimal enrichment to lower the C/O ratio of the planetary atmosphere.Our C/O ratios from the grid fitting are consistent with what was found by GRAVITY Collaboration et al. (2020).The MIRI MRS data and grid model fitting presented, however, are not enough to constrain the C/O ratio, and thus the formation history of β Pic b, better than what was presented by GRAVITY Collaboration et al. (2020); so for a complete discussion on the formation history of β Pic b, see GRAVITY Collaboration et al. (2020).

Figure 15 .
Figure15.β (ratio of radiation pressure to gravitational force) as a function of particle radius for olivine (Mg 2 SiO 4 ) and pyroxene (MgSiO 3 ) silicate grains.The black dashed horizontal line shows where β = 0.5.Grains sizes with β > 0.5 become unbound and are blown out of the system.The vertical black dashed dotted line shows a = λ/(2π) for λ = 5 μm, which is a = 0.8 μm.The small-grain limit is defined as 2πa = λ.

Figure 16 .
Figure 16.Vertical cut of the disk at 5.2 μm with the best-fit Lorentzian profile overlaid.

Figure 17 .
Figure 17.Posterior distribution for the model fitting of the β Pic b spectra and photometry with the DRIFT-PHOENIX model grid.The values show the 68% confidence intervals around the median.

Figure 18 .
Figure 18.Same as Figure 17 but for the ExoRem model grid.

Table 1
Best-fit Model Parameters for β Pictoris b