Debris Disks Can Contaminate Mid-infrared Exoplanet Spectra: Evidence for a Circumstellar Debris Disk around Exoplanet Host WASP-39

The signal from a transiting planet can be diluted by astrophysical contamination. In the case of circumstellar debris disks, this contamination could start in the mid-infrared and vary as a function of wavelength, which would then change the observed transmission spectrum for any planet in the system. The MIRI/Low Resolution Spectrometer WASP-39b transmission spectrum shows an unexplained dip starting at ∼10 μm that could be caused by astrophysical contamination. The spectral energy distribution displays excess flux at similar levels to that which are needed to create the dip in the transmission spectrum. In this Letter, we show that this dip is consistent with the presence of a bright circumstellar debris disk, at a distance of >2 au. We discuss how a circumstellar debris disk like that could affect the atmosphere of WASP-39b. We also show that even faint debris disks can be a source of contamination in MIRI exoplanet spectra.


INTRODUCTION
Transit spectroscopy is one of the best ways to study atmospheres of exoplanets.By measuring how much light from the star the planet blocks as a function of wavelength, we can learn important information about planets, such as their atmospheric composition (e.g.Kreidberg 2018).
However, transit spectroscopy has the same weaknesses as single-band transit data.One major weakness is that blended targets can result in a shallower transit depth.The most common blends are stellar companions (e.g.Damiano et al. 2017;Edwards et al. 2020), and emission from the planet itself (Kipping & Tinetti 2010;Morello et al. 2021).If the contaminant spectrum differs from that of the target star, it may imprint features on the transit spectrum and bias our interpretation of the data (e.g.Edwards et al. 2020).
Another common source of flux from planetary systems are circumstellar debris disks.Debris disks are the remains from planet formation, akin to our own Asteroid Belt or Kuiper Belt.Collisions of the planetesimals in these belts produce dust that emits in the infrared (IR) and millimeter wavelengths.This is evident from the spectral energy distribu- * NHFP Sagan Fellow tion (SED) covering those wavelengths (Wyatt 2008, and sources therein).While dust levels from the Asteroid Belt or Kuiper Belt are not yet detectable in other Solar Systems (Hughes et al. 2018), larger amounts of dust are frequently seen in other stellar and planetary systems.While often quite faint, disks can be the dominant source of flux at mid and far IR wavelengths (Chen et al. 2014).In particular, small silicate particles cause a well-studied emission feature at 10 µm (Henning 2010, and sources therein).At high enough flux ratios, which are plausible for debris disks, this extra emission could contaminate transit spectra of orbiting exoplanets.
WASP-39b is a hot Saturn around a mainsequence G8 star and a target of the JWST Transiting Exoplanet Early Release Science Program (JWST-DD-1366) and an ancillary DDT program (JWST-DD-2783).As such, its atmosphere has been extensively characterized through transmission spectroscopy with 4 different instruments over a wavelength range between 0.6 and 12 µm (Rustamkulov et al. 2023;Alderson et al. 2023;Feinstein et al. 2023;Ahrer et al. 2023;Powell et al. 2024;Welbanks & et al. 2024).Analysis from the spectra show the planet has a low C/O ratio and a high metallicity.The transmission spectrum also shows a mysterious (and unexplained) dip longward of 10 µm (Powell et al. 2024).
In this paper, we look at the case of the MIRI transit spectrum of WASP-39b and show that this dip could plausibly be caused by a circumstellar debris disk around the host star diluting the transit at wavelengths longer than 10 µm.We show that the system's SED is also consistent with a star surrounded by a debris disk.We model the debris disk and discuss the range of parameters a debris disk could have.We also discuss how a debris disk could affect the atmosphere of WASP-39b.

WASP-39 MIRI-LRS DATA
The 5-12 µm transmission spectrum of WASP-39b (Director's Discretionary Time, PID: 2783) was observed using JWST MIRI/LRS (Kendrew et al. 2015) on 2023-02-14 to confirm the presence of atmospheric SO 2 The details of the observations and analysis were presented in Powell et al. (2024).

Modelling the Stellar SED
To model the stellar flux in the IR, we used photometry from WISE (Wright et al. 2010), in addition to the out-of-transit stellar spectra taken for the transit observations.For the purposes of evaluating the stellar flux, we only used the post-transit observations in order to minimize the effect of the ramp seen at the beginning of MIRI time series observations.We used the rateints files created from the Eureka pipeline for Powell et al. (2024); those files were produced using version 0.9 of the Eureka!(Bell et al. 2022) pipeline, CRDS version 11.16.16 and context 1045, jwst package version 1.8.3 (Bushouse et al. 2022).During the Stage 1 processing, the jump rejection threshold was increased to 7 and the lastframe step was also applied to remove the excessively noisy last frame from each integration.
From there, we ran the JWST pipeline, version 1.11.0 (Bushouse et al. 2023), with the following modifications: we included the flux calibration and the pixel replace step to deal with bad pixels, and for the extraction, we used a tapered profile that is the width of 3 times the full width at half maximum.Our reference file is at https://doi.org/10.5281/zenodo.8423535.We derive uncertainties by calculating the scatter between integrations using the JWST pipeline.The data was then binned to match the resolution of the transmission spectrum from Powell et al. (2024).The resulting spectrum is plotted in Figure 1, top.We then compared the observed flux from MIRI and WISE/W3 to what we would expect from the photosphere of a star alone, using a PHOENIX/BT-Settl photospheric model at 5400 K, log(g)=4.5 at solar metallicity (Allard et al. 2003(Allard et al. , 2012) ) 1 .We scaled the model to a distance of 213.3 pc (Bailer-Jones et al. 2021), assuming a radius of 0.9 R ⊙ , and then by a factor of 0.956 to match the flux measured by WISE's W2 filter (centered at ∼4.6 µm) using synphot.The photospheric model matches the MIRI data well between 7 and 10 mum validating our choice to scale the model to the W2 flux and showing that the flux calibration of the MIRI data at those wavelengths is reasonable.Given we are in the Rayleigh-Jeans tail of the spectrum, we assume the photospheric model is accurate and do not propagate any potential uncertainties in the model flux from assumed stellar properties or the W2 flux when comparing the model flux to the observed MIRI data.The W3 bandpass goes from ∼8 to 15 µm, while the MIRI data ends at 12 µm, so we cannot directly compare the MIRI data to the W3 photometry.In the top panel of Figure 1, we see hints of excess flux past 10 µm, with the caveats that currently (i.e. using the pipeline and calibration files available in August 2023 with jwst 1112.pmap) the systematic uncertainties on the flux-calibrated stellar extraction for MIRI-LRS are likely on the order of 10% (private communication, JWST Help Desk).To verify our data reduction was not inducing the excess, we repeated our analysis on three calibrator stars, none of which showed excess.

WASP-39b Transmission Spectrum
To assess the excess flux in the transmission spectrum, we used the Eureka!transmission spectrum from Powell et al. (2024).The full reduction details are fully described in Powell et al. (2024) and closely followed those of Bell et al. (2023), and the key steps are briefly summarized here.A column-by-column background subtraction was run on each integration, the spectrum was extracted using variance-weighted optimal spectral extraction methods (Horne 1986), and lightcurves were binned to a constant 0.25 µm resolution.A starry (Luger et al. 2019) transit model was fitted to each channel assuming the orbital parameters of Carter et al. (2023) and placing a Gaussian prior on the stellar limb-darkening coefficients using ExoTETHyS (Morello et al. 2020b,a) models based on the Stagger grid (Chiavassa et al. 2018).The systematic noise model consisted of a linear trend in time, a linear decorrelation against changes in the spatial position and PSF-width, an exponential ramp in time, and a white noise multiplier.Lightcurves were fit using the No U-Turn Sampler from PyMC3 (Salvatier et al. 2016).

Atmospheric Modelling of WASP-39b
Here we aim to estimate the atmospheric spectrum of WASP-39b at longer wavelengths (e.g., ≳ 10µm) as predicted by models that have been informed by the MIRI LRS data at shorter wavelengths (e.g.,≲ 10µm).To do this, we perform an atmospheric retrieval using Au- We also plot the excess stellar flux from three G star calibrators.Note: the uncertainties are the 1σ, random uncertainties calculated from the scatter in the data.Systematic uncertainties are not included.In the case of the flux-calibrated stellar spectrum, the systematic uncertainties may be on the order of 10%.
rora (Welbanks & Madhusudhan 2021) on just those wavelengths of the MIRI data.Aurora computes the spectrum for a parallel-plane at-mosphere in transmission geometry assuming hydrostatic equilibrium.The chemical abundances are assumed to be constant with height and parameterized with a free parameter for their individual volume mixing ratios.We use the same atmospheric model setup as in Powell 2023, with the same priors and sources for absorption cross-sections.Briefly summarized here, we use an isothermal model with inhomogeneous gray clouds and power-law hazes following the single-sector prescription in Welbanks & Madhusudhan (2021), and we fit for the volume mixing ratios of H 2 O (Rothman et al. 2010) and SO 2 (Underwood et al. 2016).Additionally, we fit for the reference pressure for an assumed planetary radius of R p = 1.279R J .
The retrieved atmospheric properties are weakly constrained using this model and data, in agreement with the results from Powell et al.
After performing the retrieval, we randomly sample 200 equally weighted samples from the retrieved posterior distribution from fitting the model to the observations below 10µm and compute their associated spectrum from 5 to 12µm at a constant resolution of 10,000.These 200 spectra are then used to compute a median spectrum and 1σ and 2σ confidence intervals, shown alongside the observations in Figure 1, middle.Then, the inferred median spectrum was binned to the resolution of the data assuming a top-hat response function.This binned spectrum is used as the reference model for the remainder of this work.

Calculating Excess Flux from the Transmission Spectrum
We created an empirical excess flux model assuming that the lack of agreement between the transmission spectrum and the model is due entirely to excess flux i.e., the retrieved spectrum (Figure 1, middle) is the true, geometric transit spectrum for the planet.If there is a background source contaminating the transit, the transit will be diluted by a dilution factor as a function of wavelength, D(λ), where D can be calculated as: where λ is wavelength, F * (λ) is the stellar flux, F d (λ) is the excess flux or the disk flux, d obs (λ) is the observed transit depth and d geo (λ) is the real, geometric transit depth.Using the retrieved model as d geo (λ), the calculated transit spectrum as d obs (λ), we can easily calculate D; using the photospheric model from Section 3.1 as F * , we can also simply solve for F d (λ).In Figure 1 (bottom), we plot the empirical excess flux calculated from the dilution factor, which is 1/D, in green.We also plot the excess flux calculated from comparing the observed stellar flux to the photospheric model in blue.This gives us a spectrum of the contamination source as a function of wavelength.In both cases, at wavelengths beyond ∼10 µm, we calculate excess on the order of a few percent

Modeling the Excess Flux as a Debris Disk
To find the basic characteristics of a debris disk that could explain the excess flux, we generate simple models of a narrow circular ring of optically thin dust grains in radiative equilibrium with the star.The flux density emitted by the disk is entirely determined by the semimajor axis of the ring (and therefore its temperature) and the amount of dust it contains.Using these two free parameters, we try models of two types: 1) dust grains that behave like blackbodies and 2) small grains composed of crystalline olivine that have a pronounced emission peak at ∼11µm.Although we could generate much more complex models that distribute the dust in a broad ring and have multiple compositional components, even a very simple model can explain the data.For both model types, we assume a stellar luminosity of 0.63 L ⊙ and T eff =5400 K, and we test fits to both methods of calculating the excess flux (from the stellar and transmission spectra as described above) while also constraining the model to be consistent with the W4 upper limit.WISE W3-band shows a small excess over the photosphere and over the JWST spectra (1.9σ) that we do not attempt to model.We use the least squares method to find best fit models by minimizing the chi-squared goodness-of-fit metric.Because the W4 point is a 95% confidence upper limit, we fit models constrained by 30 different W4 values spanning the range of possibilities from the 99% confidence upper limit (1.94 mJy) down to a value equal to an excess equal to the stellar photosphere (0.69 mJy).Overall, we used 30 different adopted W4 points, evenly spaced in flux.We calculate formal uncertainties on the two parameters, using formulae for the variation in chi-squared as described e.g. in Bevington & Robinson (2003), but because of the large uncertainties and short wavelength coverage and existence of only a W4 upper limit, the 3σ ranges are quite large.

RESULTS
The excess derived from the transmission spectrum has larger uncertainties and is compatible with a larger range of disk models including those that best fit the data from the stellar spectrum method.For blackbody grains, the best fits (shown in orange solid lines in Figure 2, parameters in Table 1) have distances from the star of ∼3.4 and 4 au (dust temperatures ∼135 and 124 K), for the transmission Figure 2. The best fits to the excess flux as calculated from the transmission spectrum (top) or from the stellar spectrum (bottom).We fit a crystalline silicate feature and a blackbody curve to both spectra separately.These plots assume the 3σ upperlimit to the WISE/W4 flux.
and stellar spectrum respectively.The uncertainties on the flux densities imply that a wide range of dust temperatures are actually allowed by the data; a 3σ lower limit of about 1.3 au is set by the lack of excess shorter than 10µm.The effect of including the W4 upper limit is to set an outer radius for the disk, because cooler grains further from the star produce more 25 µm emission.The maximum distance allowed by the W4 upper limit is about 8 au.
For crystalline silicate grains, the best fit temperatures are somewhat cooler, with best fit disk radii of 13 or 16 au for the transit and stellar spectra, respectively, and a minimum ring size of 4 au and a maximum ring size of about 30 au.
One important difference between the models is what they predict for the 12-25 µm emission, because the blackbody models continue to rise monotonically past 12 µm.Cooler dust is allowed in the silicate models because of the emission peaks are fairly narrow.
All best fit disk models are relatively bright, with L IR /L * > 3× 10 −4 .L IR /L * is the modeled disk luminosity in the IR integrated over wavelengths out to 1 mm divided by the total stellar luminosity.For debris disks detected in other systems, L IR /L * ranges from ∼10 −5 to ∼10 −2 over a range of temperatures usually <100 K. Below that systems are too faint too detect; above that are protoplanetary disks.In comparison, the Asteroid Belt has L IR /L * ∼10 −7 (Hughes et al. 2018).While we do not measure the mass directly, as measurements out to 12 µm are only sensitive to small dust grains, this implies the disk is fairly massive.Micronsized dust in mature debris disks come from the collisions of planetesimals, so a detection of dust requires the presence of larger planetesimals in the debris disk (Wyatt 2008, and sources therein).The mass in dust can be esimated from L IR /L * , given the typical dust size and density and the simple assumption that the grains absorb in proportion to their geometric cross sections (Chen & Jura 2001).For grains to act like blackbodies, they must be larger than the observed wavelength, so we can assume a size of 50 µm and a density typical of slightly porous silicates of 2.7 g cm −3 .For a ring radius of 4 au, the mass is ∼8×10 20 kg, or similar to the most massive asteroid in the solar system.The original mass in parent bodies would presumably be larger than this, as there would actually be dust of a range of sizes from 50 µm on up, which if described as a power law puts most of the mass in the largest bodies.The crystalline silicate model requires small grains, much smaller than the wavelength of observations, say 0.5 µm.The mass in small grains is an order of magnitude lower, ∼4×10 19 kg for these models.
Despite our simple models, we found that we could successfully model the data as excess flux arising in a disk.All of our models result in a better model fit to the data at wavelengths >9 µm (as measured by χ 2 ν ) than assuming no excess (Table 1).This would propagate to the exoplanet transmission spectrum, too; allowing for even these simple debris disk models results in a better fit.These disk models result in the transmission spectra and residuals shown in Figure 3.Given the systematic uncertainties, we cannot conclusively prove that the dip in transmission spectrum or the excess flux in the stellar spectrum is caused by a debris disk, but we can clearly show that a debris disk is a plausibleand as of now, the only -explanation.Photometric observations at longer wavelengths could confirm or reject this hypothesis.C/O ratio and metallicity can, in theory, be used to trace the formation location of an exoplanet ( Öberg et al. 2011;Espinoza et al. 2017).However, a low C/O ratio cannot uniquely trace a formation location, largely because a low C/O ratio can be reached via a number of paths, including forming within the H 2 O ice line or significant planetesimal impacts.Characterizing the planetesimal population in a debris disk would allow us to better understand planet formation, because it will help break this degeneracy between formation location and planetesimal impacts.
Analysis of near-IR JWST data for WASP-39b (Rustamkulov et al. 2023;Alderson et al. 2023;Feinstein et al. 2023;Ahrer et al. 2023;Welbanks & et al. 2024) showed that the system likely has a sub-stellar C/O ratio and supersolar metallicity.Feinstein et al. (2023) speculated that one way these could have occurred was by significant planetesimal accretion after the planet formed with planetestimals from 2-10 au.The potential debris disk described by  Note-The parameters for our four best fit models (shown in Figure 2) assuming the 3σ upperlimit to the WISE/W4 flux.The χ 2 ν for the No Excess Flux Fits were calculated with a flat line at 0 Jy.several of our best fit models (i.e.large and between 4-10 au) would be consistent with the source of planetesimals needed to lower WASP-39b's C/O ratio and raise its metallicity.

Potential Contamination of MIRI Data for Other Exoplanets
Regardless of whether WASP-39 has a debris disk, debris disk contamination is a serious potential issue for JWST/MIRI observations of exoplanets.In Figure 4, we plot the magnitude of the dip induced by varying amounts of excess flux from a disk based on Equation 1.
The decrease in transit depth due to a faint debris disk, even fainter than the potential disk around WASP-39, can be on the order of several hundred parts-per-million (ppm).This is larger than the amplitude of many potentially detectable spectral features produced by exoplanets in this wavelength range.It is also comparable to or larger than typical uncertainties in transmission spectra of exoplanets.
Since JWST/MIRI is sensitive to excesses induced by debris disks on the order of 1%, we need to be able to measure fluxes to that level to determine whether debris disk contamination is likely to be an issue for any given system.Unfortunately, WISE alone is not sensitive enough.For exoplanet targets observed thus far with JWST/MIRI (Nikolov et al. 2022), the median uncertainty in WISE/W3 for the system is ∼2.2%, with ∼85% of targets having a W3 uncertainty less than 4.5%, but no target had an uncertainty under 1%.
We can also consider population statatistics to estimate what percentage of systems will have excess flux at the 1% level.The most sensitive population surveys for warm dust use Spitzer Space Telescope or WISE data.WISE/W3 constrains this to some extent for particularly bright and warm disks.Several studies have looked at the percentage of systems with excess from debris disks.Absolute photometric calibration is the main source of uncertainty in population studies.The largest surveys, which use WISE W3 or W4 (12 or 22 µm) photometry, such as Kennedy & Wyatt (2013) and Patel et al. (2014), the typical excesses detected are >15% of photospheric flux in at least one band.We then need to extrapolate to estimate the probability of disks with excess of ∼1%.Kennedy & Wyatt (2013) estimate that 3% of Sun-like stars will have excess greater than 1% at 12 µm from circumstellar disks.The LBTI is more sensitive to low levels of dust; Ertel et al. (2020) showed that 20% of systems are "significantly dusty" based on their flux levels at 12 µm.

CONCLUSION
We have shown that the emitted flux from circumstellar debris disks could contaminate MIRI transmission spectra of exoplanets, and may be a potential explanation for the the dip in the transit spectra seen in Powell et al. (2024).We showed that the stellar spectrum also has excess flux at wavelengths longer than 10 µm, again consistent with a debris disk.If there is a circumstellar debris disk around WASP-39b caus- ing these features in the spectra, it is relatively large (L IR /L * >10 −4 ), further out than 2 au, and could be the reason why WASP-39b has a sub-stellar C/O ratio and a super-solar metallicity.Data at longer wavelengths could confirm the presence of the disk and better constrain its properties.
ence Institute.The specific observations analyzed can be accessed via DOI.
Software: SpectRes (Carnall 2017), NumPy (Oliphant 2006;Van Der Walt et al. 2011), Mat-plotlib (Hunter 2007), astropy (Collaboration et al. 2013) A star whose flux partly falls into the aperture used for spectral extraction is a stellar blend.It could be either a gravitationally bound companion or a chance aligned background star.WASP-39 is a single star (Faedi et al. 2011;Mancini et al. 2018).We used tpfplotter (Aller et al. 2020) to search for other potential blends from the GAIA DR3 catalog, finding only two significantly fainter sources (∆mag =4.25 and 5.32) within ∼1 ′ .We will show that they cannot be the cause of the observed spectroscopic feature.
Figure 5 reports the inverse of the dilution factor (or the excess flux), normalized to the 5-5.25 µm bin, for a variety of stellar contaminants.We adopted the MPS-Atlas of stellar spectra (Kostogryz et al. 2023).In order to rise by a few percent at ≳10 µm, the contaminant should be an M dwarf or cooler star at roughly the same distance of WASP-39 b.We note that such a companion would critically affect the transmission spectrum even at shorter wavelengths.

A.2. Planetary emission
The emission from the planet itself may cause a similar dilution effect to that of a stellar blend, but typically smaller (Kipping & Tinetti 2010;Martin-Lagarde et al. 2020).However, the socalled planet self-blend effect can be significant in the infrared, where planets have their peak emission.Morello et al. (2021) estimated the possible self-contamination bias in the JWST transit spectra for a list of exoplanets, finding effects below 50 ppm for WASP-39 b.We conclude that planetary emission cannot cause the observed variation in transit depth at wavelengths >10 µm.

Figure 1 .
Figure 1.(top) A portion of the observed flux from the WASP-39 system along with a BT-Settl photospheric model, using WISE-W2 as an anchor point.Both the MIRI data, as well as the W3 point from WISE hint at a slight amount of excess flux.(middle) The WASP-39b MIRI transmission spectrum along with the best-fit model retrieved to the data short of 10 µm, along with 1σ uncertainties.(bottom) The excess flux calculated using the stellar flux compared with a photospheric model flux from the top plot (in blue) and the excess flux calculated from the transmission spectrum (in green).We also plot the excess stellar flux from three G star calibrators.Note: the uncertainties are the 1σ, random uncertainties calculated from the scatter in the data.Systematic uncertainties are not included.In the case of the flux-calibrated stellar spectrum, the systematic uncertainties may be on the order of 10%.
The Effect of a Debris Disk on the C/O ratio and Metallicity of WASP-39b

Figure 3 .
Figure 3. (top) WASP-39b MIRI observed transmission spectrum along with its model retrieval and the model retrievals adjusted by the dilution factor induced by our best debris disk models.(bottom) The residuals between the data and the models.The inclusion of the debris disk decreases the residuals.

Figure 4 .
Figure 4.The decrease in transit depth due to a relatively faint debris disks can be on the order of several hundred ppm.

Figure 5 .
Figure 5. Dilution factor on the MIRI transit spectrum for a hypothetical stellar blend with various temperatures.The labels in the legend indicate the effective temperature and flux fraction from the blended star (1.0 is the maximum possible fraction).

Table 1 .
Best Fit Debris Disk Model Parameters