The James Webb Interferometer: Space-based Interferometric Detections of PDS 70 b and c at 4.8 μm

PDS 70 is one of the most extensively studied young stellar systems. It is the only known multiplanet protoplanetary disk system, where two or more planets have been robustly detected within the disks from which they formed (M. Keppler et al. 2018; S. Y. Haffert et al. 2019). The disk consists of a large outer component and a smaller inner component (M. Keppler et al. 2018) separated by a wide cavity (J. Hashimoto et al. 2012; S. Facchini et al. 2021) spanning between ≲17 au and ~54 au (M. Keppler et al. 2018). The inner disk was first inferred via spectral energy distribution (SED) fitting (R. Dong et al. 2012; J. Hashimoto et al. 2012) and was subsequently isolated in the near-infrared (M. Keppler et al. 2018) and submillimeter (M. Keppler et al. 2019). Recently, observations with JWST/MIRI detected water vapor in the inner disk (G. Perotti et al. 2023). The outer disk of PDS 70 is not significantly perturbed, in contrast to some other protoplanetary disk systems where planet candidates have been detected, such as MWC 758 (K. Wagner et al. 2023) or AB Aurigae (T. Currie et al. 2022). The PDS 70 outer disk appears as a close-to-symmetric ring at submillimeter wavelengths (M. Keppler et al. 2019; M. Benisty et al. 2021), with one significant azimuthal asymmetry observed to the northwest of the star that is close to coincident with an armlike structure seen in the near-infrared (J. J. Wang et al. 2020; S. Juillard et al. 2022).

The planets PDS 70 b and c have been extensively imaged in the near-infrared from ~1 to 3.8 μm (e.g., M. Keppler et al. 2018; A. Mülleret al. 2018; V. Christiaens et al. 2019; D. Mesa et al. 2019; J. J. Wang et al. 2020, 2021b; G. Cugno et al. 2021), with a few detections beyond ~4–4.8 μm (T. Stolker et al. 2020; V. Christiaens et al. 2024), as well as detections in H_α (K. Wagner et al. 2018; S. Y. Haffert et al. 2019). Circumplanetary disk (CPD) emission from PDS 70 c has been directly detected in the submillimeter with the Atacama Large Millimeter/submillimeter Array (ALMA), coincident with detections in the near-infrared (A. Isella et al. 2019; M. Benisty et al. 2021). Submillimeter emission has also been observed that is seemingly associated with PDS 70 b; however, it is offset from the near-infrared detections (A. Isella et al. 2019; M. Benisty et al. 2021). Additionally, emission at the L5 Lagrangian point of PDS 70 b has tentatively been detected in the submillimeter with ALMA (O. Balsalobre-Ruza et al. 2023). Finally, a third point source has been tentatively detected in the PDS 70 disk cavity, using the Very Large Telescope (VLT) SPHERE over multiple epochs (D. Mesa et al. 2019). The nature of the source is still not well understood, as it has a distinct spectrum compared to both PDS 70 b and c, more similar to the spectrum of PDS 70 A between ~1 and 1.6 μm, implying that it may be scattered light (D. Mesa et al. 2019). This source was tentatively redetected at 1.9 μm with JWST/NIRCam (V. Christiaens et al. 2024).

Beyond ~4 μm, emission from warm (≳100 K) CPD material is expected to contribute comparably (and greater at longer wavelengths) to atmospheric emission (e.g., J. J. Wang et al. 2021b; G. Cugno et al. 2024). So, observations into the mid-infrared may allow for the presence of CPD emission to be inferred by looking for excess emission compared to atmospheric models. J. J. Wang et al. (2021b) showed that the VLT/NACO M$^{\prime}$ detection of PDS 70 b from T. Stolker et al. (2020) provided weak support for excess blackbody-like emission, so follow-up observations at ~4.8 μm and beyond will allow for this tentative detection of CPD emission in the mid-infrared to be confirmed. Recently, V. Christiaens et al. (2024) redetected PDS 70 b and c with the JWST/NIRCam (4.8 μm) F480M filter. Their detection of PDS 70 b was in agreement with the M$^{\prime}$ detection but with a larger uncertainty.

In our work, we redetect PDS 70 b and c at F480M using the aperture masking interferometry (AMI) mode on JWST/NIRISS, an independent method to the JWST/NIRCam results, presented by V. Christiaens et al. (2024), and make the most precise measurement of the flux of PDS 70 b at 4.8 μm using the JWST/NIRISS F480M filter. We also detect PDS 70 c using the F480M filter at a similar precision to V. Christiaens et al. (2024). Here, using the power of NIRISS AMI, as introduced in previous works using AMI (S. Ray et al. 2023; R. M. Lau et al. 2024; S. Sallum et al. 2024), we detect PDS 70, its outer disk, and its two protoplanets, b and c. These are the first planets detected with space-based interferometry.

This paper is structured as follows. Section 2 outlines the data acquisition, reduction, and cleaning procedure; Section 3 details the disk plus two-planet model fitting and the conversion of the measured contrasts to fluxes; Section 4 outlines the derived planet parameters and the analysis of the data for any signal beyond the two known planets and the disk; and Section 5 discusses the implications of the measured contrasts and the nature of the tentative detection of residual emission in the data. We close with a summary in Section 6.

We observed PDS 70 as a part of the NIRISS GTO program (PID 1242; PI: D. Johnstone). The data presented in this paper were acquired on 2023 February 24 using the James Webb Interferometer, the AMI mode of the NIRISS instrument (A. Sivaramakrishnan et al. 2012, 2023; R. Doyon et al. 2023), in the F480M filter (λ = 4.815 μm, Δλ = 0.298 μm).²⁰ The SUB80 subarray was used along with the NISRAPID readout pattern. The point-spread function (PSF) calibrator star HD 123991 was also observed using the same configuration. Both targets were observed using no dithers and no rolls. The PDS 70 data set consists of 96 groups and 418 integrations, while the brighter HD 123991 data set consists of 29 groups and 418 integrations. This corresponds to exposure times of 50.45 minutes for PDS 70 and 15.24 minutes for HD 123991.

The data were reduced using version 1.11.3 of the jwst pipeline²¹ (H. Bushouse et al. 2023), with Calibration Reference Data System context jwst_1110.pmap, from the uncal data format. We use the default calwebb_detector1 stage 1 pipeline and calwebb_image2 stage 2 pipeline, skipping only the IPC, photom, and resample steps. Note that we do not use the charge_migration²² step (P. Goudfrooij et al. 2024), which was designed to minimize the impact of the detector brighter-fatter effect (BFE). The charge_migration algorithm discards all data in each pixel above a set signal threshold, along with the rest of the data in the four directly adjacent neighbors. This strategy is not optimal for AMI data because, as we show in Appendix B and has been shown in other works (e.g., I. Argyriou et al. 2023), charge bleeding from a bright pixel into its neighbors affects not only the four adjacent pixels but also the four pixels along the diagonals, albeit somewhat less. Additionally, considering only the brightness of individual pixels is not the optimal threshold criterion because the BFE depends on the contrasts between pixels (e.g., A. A. Plazas et al. 2017). Due to the fine structure in the AMI PSF, which is somewhat undersampled at 4.8 μm with the NIRISS plate scale of ~66 mas, there are many parts of the PSF with large contrasts between neighboring pixels. The current pipeline BFE mitigation algorithm, therefore, is not useful for correcting AMI data. A modification to the existing pipeline, beyond the scope of this work, could be to base the threshold on the contrast between adjacent and diagonal neighboring pixels.

We investigate the effect of the BFE on the data (see Appendix B for details), as it has been shown to be the limiting noise floor for AMI observations (J. Kammerer et al. 2023; S. Ray et al. 2023; A. Sivaramakrishnan et al. 2023; S. Sallum et al. 2024). Based on our analysis, and to minimize the effect of the BFE on the calculated squared visibilities and closure phases, we discard all groups beyond where the sum of the central 3 × 3 pixels of the PSF is above 30,000 data number (DN), or ~48,000 electrons, in the linearized-ramp-level data. We use this cutoff, as opposed to the intensity in the central pixel, because it is more robust against changes in PSF centering. This constraint corresponds to an apparent systematic change in the calculated rate in the central pixel of ~1% and minimal bleeding into surrounding pixels (Appendix B, Figure 11). We note that this loss in signal biases our measured contrasts by significantly less than 1% due to the central pixel accounting for ~1% of the total signal in the PSF. This can be seen in Figure 1, which presents the mean PDS 70 image-plane data using only the first 28 groups. We also note that this BFE mitigation method is only possible because PDS 70 and HD 123991 are relatively faint, such that we have at least 10 groups before a peak signal level of ~5800 DN. With fewer groups, the data are likely to be dominated by the 1/f noise seen in the NIRISS detector (e.g., S. Sallum et al. 2024). For PDS 70, we are therefore limited to 28 out of the total 96 groups, while for HD 123991, we are limited to 10 out of the total 29 groups.

Zoom In Zoom Out Reset image size

We correct bad pixels that were flagged by the pipeline by finding the pixel value that minimizes the power outside of the Fourier support of the AMI mask (M. J. Ireland 2013). We also perform subpixel centering of the data using a Fourier-based shift by finding the position that minimizes the absolute value of the phase calculated from the mean of the cleaned data (J. Kammerer et al. 2019).

The background star near PDS 70 that was originally observed by P. Riaud et al. (2006) and confirmed to be a background star by J. Hashimoto et al. (2012) is at the edge of the SUB80 array, at a separation of ~2'' from PDS 70 A. This is outside of the interferometric field of view, set by the shortest baseline, but the signal will still contaminate all baselines. To minimize the contribution to the calculated interferometric observables, a super-Gaussian window (${e}^{-{(r/\sigma )}^{4}}$), with σ = 30 pixels, is applied to both the PDS 70 data and the HD 123991 data.

We calculate the squared visibilities, V², and closure phases of PDS 70 and calibrator HD 123991 using amical (A. Soulain et al. 2020; A. Soulain & C. M. T. Robert 2023). We also compute the statistical uncertainties (i.e., the standard error of the mean) of these quantities across integrations. The interferometric observables are calibrated to remove the instrumental bias by dividing the V² of PDS 70 by the V² of HD 123991. Similarly, the closure phases are calibrated by subtracting the closure phases of HD 123991 from those of PDS 70. The calibrated closure phases and squared visibilities are shown in Figure 2.

Zoom In Zoom Out Reset image size

Due to the interferometric nature of the data, we construct a model to jointly fit for flux from the star, extended disk emission, and the two known planets. This is required because the signal from all components is entangled and thus cannot be measured independently. We model the star and the two planets as delta functions, with their analytic Fourier transform given by

$$\begin{eqnarray}\begin{array}{rcl}{V}_{* ,b,c}(u,v) & = & {I}_{* }+{I}_{b}{e}^{-i2\pi (u{x}_{b}+v{y}_{b})}\\ & & +{I}_{c}{e}^{-i2\pi (u{x}_{c}+v{y}_{c})},\end{array}\end{eqnarray} \tag{ 1 }$$

where I_* is the brightness of the star, I_b is the brightness of planet b, and I_c is the brightness of planet c. Additionally, x_b, y_b, x_c, and y_c denote the x and y offset of planets b and c relative to the star (which is fixed at the phase center), and u and v are the baseline coordinates in the Fourier domain in units of wavelength.

To model the extended disk emission, we use a simple geometrical model, as has been used to successfully reconstruct the extended emission of LkCa 15, in the near-infrared with ground-based AMI (D. Blakely et al. 2022; S. Sallum et al. 2023). A simple model is sufficient due to the sparsity of the data and the fact that the dominant signal is expected to be from the forward-scattering peak of the disk. To construct this model, we assume that all of the disk signal that we are observing is either scattered light or optically thick emission from the disk surface, in a similar manner to the procedure used by T. Stolker et al. (2016) and E. Sissa et al. (2018). We describe the height of the surface of the disk using a power-law profile, given by

$$\begin{eqnarray}&&z={H}_{100}{\left(\displaystyle \frac{\sqrt{{x}^{2}+{y}^{2}}}{100\,\mathrm{au}}\right)}^{\beta },\end{eqnarray} \tag{ 2 }$$

where x and y are the potentially rotated and/or inclined coordinates of the disk midplane in units of au, H₁₀₀ is the aspect ratio at 100 au (~0$\mathop{.}\limits^{\unicode{x02033}}$9), and β is the flaring angle. To project this geometry into the plane of the sky, we define r to be the apparent distance from the star to the disk surface, given an observed inclination angle, i, relative to the plane of the sky,

$$\begin{eqnarray}&&r=\sqrt{{x}^{2}+{(y+z\sin (i))}^{2}+{z}^{2}}.\end{eqnarray} \tag{ 3 }$$

We describe the radial brightness distribution as a skewed Gaussian projected onto the flared disk geometry, given by ${I}_{{d}_{1}}(r)\cdot {I}_{{d}_{2}}(r)$, where ${I}_{{d}_{1}}(r)$ is defined as

$$\begin{eqnarray}&&{I}_{{d}_{1}}(r)=\exp \left(-\frac{{(r-{r}_{0})}^{2}}{2{\sigma }_{r}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

and ${I}_{{d}_{2}}(r)$ is

$$\begin{eqnarray}&&{I}_{{d}_{2}}(r)=\frac{1}{2}\left(1+{\rm{erf}}\left(\alpha \frac{(r-{r}_{0})}{\sqrt{2}{\sigma }_{r}}\right)\right),\end{eqnarray} \tag{ 5 }$$

where r₀ describes the location of the peak of the disk brightness, σ_r describes the radial extent of the disk emission, and α describes the degree to which the inner edge of the disk is truncated.

The azimuthal brightness distribution is the sum of an axisymmetric component and an asymmetric component. For the asymmetric component, we explore three independent parameterizations: a Henyey–Greenstein model (L. G. Henyey & J. L. Greenstein 1941), a Gaussian model, and a power-law model. The Henyey–Greenstein model is given by

$$\begin{eqnarray}&&{I}_{{d}_{3}}(\theta )=\frac{1-{g}^{2}}{4\pi {(1+{g}^{2}-2g\cos (\theta -{\theta }_{0}))}^{3/2}},\end{eqnarray} \tag{ 6 }$$

where g is the scattering parameter, a free parameter constrained to be between 0 (isotropic scattering) and 1 (forward scattering); θ is the azimuthal angle in the disk midplane (defined by the coordinates x and y); and θ₀ is the azimuthal location of the scattering peak, which is fixed to the minor axis on the near side of the disk. The Gaussian model is given by

$$\begin{eqnarray}&&{I}_{{d}_{3}}(\theta )=\exp \left(-\frac{{(\theta -{\theta }_{0})}^{2}}{2{\sigma }_{\theta }^{2}}\right),\end{eqnarray} \tag{ 7 }$$

where σ_θ describes the azimuthal extent of the disk emission. The power-law model is given by

$$\begin{eqnarray}&&{I}_{{d}_{3}}(\theta )=\cos {\left(\frac{1}{2}(\theta -{\theta }_{0})\right)}^{N},\end{eqnarray} \tag{ 8 }$$

where N describes the azimuthal extent of the disk emission.

Combining Equations (2)–(5) and any one of the scattering phase functions given by Equations (6), (7), or (8), the full disk model is given by

$$\begin{eqnarray}&&{I}_{d}(r,\theta )=\left({A}_{a}{I}_{{d}_{3}}(\theta )+{A}_{s}\right){I}_{{d}_{1}}(r){I}_{{d}_{2}}(r),\end{eqnarray} \tag{ 9 }$$

where A_a and A_s control the brightness of the asymmetric and symmetric components, respectively. The Fourier transform of Equation (9) is calculated using the one-sided discrete Fourier transform (DFT) as implemented in XARA (F. Martinache 2010, 2013; F. Martinache et al. 2020) to exactly compute the DFT at the baseline coordinates of the NIRISS aperture mask baselines.

Finally, any overresolved emission is included in the model using a single free parameter I_o that contributes additional flux to the normalization of the visibilities. The full model consists of the sum of Equation (1) and the DFT of Equation (9), normalized by the total flux of the model. The complex visibilities of these 17 parameter models are given by

$$\begin{eqnarray}&&{V}_{m}=\frac{{ \mathcal F }\{{I}_{d}\}+{V}_{* ,b,c}}{{\rm{\Sigma }}{I}_{d}+{I}_{* }+{I}_{b}+{I}_{c}+{I}_{o}}.\end{eqnarray} \tag{ 10 }$$

We calculate the closure phase of the model by summing the phase measured between baselines that form a (closure) triangle for all 35 triplets of holes in the aperture mask, corresponding to 15 independent measurements. Squared visibilities are the squared amplitude of the complex visibilities, corresponding to each of the 21 baselines formed by pairs of mask holes.

To fit our model to the data, we further construct a more robust set of observables from the squared visibilities and closure phases, following a similar approach to the Bayesian analysis in J. W. Xuan et al. (2024; and described in more detail in W. Thompson et al. 2024, in preparation). We calculate self-calibrating log-closure amplitudes, which calibrate out hole-based gain terms affecting the squared visibilities. We calculate the log-closure amplitude uncertainties using a Monte Carlo approach by calculating the standard deviation of 1,000,000 noise realizations from the statistical uncertainties of the squared visibilities. We also consider correlations between shared baselines, using analytic correlation matrices for the log-closure amplitudes and closure phases. For the closure phases, we use the closure phase correlation matrix, C_Ψ, described by J. Kammerer et al. (2020; with no spectral component). We construct a similar log-closure amplitude correlation matrix, C_a, following L. Blackburn et al. (2020). The closure phase and log-closure amplitude correlation matrices can be constructed using Equations (B5) and (B14) derived by L. Blackburn et al. (2020), respectively, setting $\sigma =1/\sqrt{3}$ for closure phases and $\sigma =1/\sqrt{4}$ for log-closure amplitudes. Finally, we construct linearly independent closure phase and log-closure amplitude observables. We calculate projection matrices by taking the singular value decomposition of TT^T, where T is the design matrix of either the closure phases or log-closure amplitudes, as defined by L. Blackburn et al. (2020), and only taking the nonzero singular values. This is similar to what is described by M. J. Ireland (2013) for closure phases. We henceforth refer to these new observables as kernel phases and log kernel amplitudes. The advantage of using the "kernel" observables is twofold. By using closure phases, we are "double counting" the data in the likelihood calculation, leading to overly confident results and underestimated parameter uncertainties. The kernel approach also allows for the covariance matrices, described above, to be used in the likelihood, as the matrices are rank-deficient, whereas this is not the case in the new basis.

We use a Bayesian modeling approach to measure the locations and contrasts of PDS 70 b and c, as well as our disk model parameters. We estimate our model posterior using dynamic nested sampling (E. Higson et al. 2018) with dynesty (J. S. Speagle 2020; S. Koposov et al. 2023). We use uniform priors on the contrasts, separations, and position angles of both planets. Our contrast prior ranges from 0 to 10 mag. Our location priors for PDS 70 b are from 100 to 250 mas in separation and 70°–160° in position angle. Similarly, for the position of PDS 70 c, we use a range of 100–250 mas in separation and 230°–320° in position angle. We use priors on the disk geometry from M. Keppler et al. (2018, 2019) that are described in Appendix A. To account for an unknown level of miscalibration seen in the squared visibilities, which is also present in the log-closure amplitudes, we include a systematic uncertainty, σ_a,sys, term in our model, added in quadrature with the statistical uncertainty, that is recovered along with our star plus disk plus planet model parameters.

We estimate parameters for both PDS 70 b and c using Bayesian model averaging over the three asymmetric scattering models we explore. We follow the approach used by E. Nasedkin et al. (2024) and M. C. Nixon et al. (2024) of assuming that each model is equal likely a priori and thus combine the posterior distributions by weighting each model by their Bayesian evidence (listed in Table 2 in Appendix A).

To convert the measured contrasts to flux measurements, we use aperture photometry. We follow the procedure outlined in Section 2 of reducing and cleaning the data. The only difference is that we now apply the photom step, so that the data are in units of surface brightness. We crop the data to a size of 67 × 67 pixels, with PDS 70 A at the center. We calculate the total flux within a 67 pixel diameter circular aperture and recover a total flux of 128.6 mJy. We note that this value includes resolved disk flux, which can be clearly seen by looking at the squared visibilities decrement from 1, seen at all baselines, in Figure 2. From the visibility amplitudes, we calculate a resolved flux component of 5.2% ± 0.6%, independent of baseline. Additionally, our photometry is missing a nonnegligible amount of flux that is outside of the field of view. We calculate a correction factor by calculating PSF models using WebbPSF (M. D. Perrin et al. 2012). We calculate a 67 × 67 pixel model and 2000 × 2000 pixel model, both with the native NIRISS pixel size. We calculate the correction factor by dividing the total flux within a 2000 pixel diameter circular aperture for the full PSF and a 67 pixel diameter for the cutoff PSF, giving a correction factor of 1.12. After applying these correction factors, we find a flux of PDS 70 A (and the unresolved inner disk emission) of 136.4 ± 6.8 mJy, assuming a 5% absolute flux calibration uncertainty.²³

4.1. Derived Planet Parameters

Figure 3 shows the two-planet plus disk model generated from the median parameters from the power-law scattering model (Equation (8)) model posterior. We focus on the power-law model in particular, because it had the highest Bayesian evidence out of the three models we explored (Table 2 in Appendix A). The 1σ, 2σ, and 3σ contours from the marginalized posterior of the planet locations are shown by the white contours. A comparison of the model log kernel amplitudes and kernel phases to the data is given in Figure 4, showing that the model provides an excellent fit. Furthermore, the measured positions of both PDS 70 b and c are consistent with VLTI/GRAVITY predictions from J. J. Wang et al. (2021b), which are shown by the dashed gray circles in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1 summarizes the derived planet parameters, calculated from the Bayesian average of all three models. The disk model parameters are listed in Table 2 in Appendix A. We find that both planets are found at a high signal-to-noise ratio (SNR), 14.7 for PDS 70 b and 7.0 for PDS 70 c (where we define the SNR as the contrast divided by the uncertainty in the contrast marginalized over all of the model parameters), in locations that are consistent between the three models (planet parameters for the individual models are shown in Table 3 in Appendix A). We also find that the measured contrasts of the planets are consistent between the three models.

Table 1. PDS 70 b and c Parameters

	PDS 70 b	PDS 70 c
Separation (mas)	$150.{5}_{-10.9}^{+9.5}$	$218.{4}_{-8.2}^{+7.2}$
Position angle (deg)	$131.{2}_{-1.6}^{+1.5}$	$270.{0}_{-1.7}^{+1.7}$
Contrast (Δmag)	$5.8{4}_{-0.07}^{+0.07}$	$6.4{8}_{-0.14}^{+0.17}$
Fluxa (μJy)	$629.{7}_{-39.0}^{+42.2}$	$349.{5}_{-49.6}^{+48.5}$

Note. We report the median along with the 16th and 84th percentiles of the Bayesian average of the marginalized posteriors over the three disk (Equations (6), (7), and (8)) plus star plus planets models we explored.

^aWe convert the planet contrasts to flux using our measured star plus unresolved inner disk flux of 136.4 ± 6.8 mJy.

Download table as:  ASCIITypeset image

Figure 5 shows a corner plot of the planet parameters for each of the three models, as well as the Bayesian average of the three models. We note that there are correlations between the parameters of both planets. Notably, there is a moderate anticorrelation between the separation and the contrast of PDS 70 c. From the correlation matrices of the model parameters, shown in Figure 6, there are no strong correlations between the disk parameters and the planet contrasts, indicating that the disk model has a minimal effect on the measured contrasts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Residual Signal and Contrast Limits

We compute an SNR map, shown in Figure 7, made by dividing the best-fit point-source binary contrast to the residual kernel phase signal by its associated uncertainty fitting.

Zoom In Zoom Out Reset image size

Because closure phases are nonlinear, we calculate the complex visibilities of the power-law model using the median of the marginalized posterior for all parameters of the power-law model and add this to our binary model visibilities before calculating the model closure phases and then kernel phases. Each cell in the grid shows the best-fit contrast divided by the uncertainty in the contrast calculated using the Laplace approximation (L. Tierney & J. B. Kadane 1986), assuming the uncertainties are Gaussian. We clearly see that there is residual emission at an SNR of ~4 to the southwest of PDS 70 A and very close to the inner working angle of 94 mas that is not fully captured by our two-planet plus disk model. To assess the nature of the signal that we observe, we use nested sampling to fit a star plus planet model to the residual kernel phases and find a contrast of $7.{6}_{-1.3}^{+0.8}$ mag, a separation of $11{8}_{-64}^{+52}$ mas, and a position angle of $22{0}_{-15}^{+10}$@. The position angle of the signal is reasonably well constrained; however, there is a strong correlation between the separation and contrast, which is expected for a high-contrast source near the diffraction limit. Interestingly, the position angle of the observed signal is not consistent with forward scattering from an inner disk that has the same geometry as the outer disk. However, due to the unconstrained nature of the separation of the emission, it is unclear whether the signal is due to compact emission or an inner disk feature.

To assess the limits that we are able to place on the presence of additional planets in the system, we calculate a 5σ contrast curve from the residual kernel phase signal of the power-law model. We use the approach described above to fit to the residual kernel phase signal and calculate Bayesian upper limits. To do this, we adapt the method presented by J. B. Ruffio et al. (2018) to interferometric data. We first calculate a grid of likelihoods for a binary model at fixed R.A. and decl. locations for a range of contrasts. We next use the maximum-likelihood point in each grid cell to initialize gradient descent, using the BFGS algorithm (C. G. Broyden 1970; R. Fletcher 1970; D. Goldfarb 1970; D. F. Shanno 1970), to find the contrast that maximizes the likelihood. We use the calculated contrast in each cell along with the Laplace approximation to calculate the uncertainty in the calculated contrast at each R.A. and decl. point. Finally, combining the calculated contrasts, uncertainties, and a 99.999971% (equivalent to a 5σ) cutoff probability and using Equation (8) in J. B. Ruffio et al. (2018), we calculate the 5σ upper limit at each point within a radius of 400 mas, which we azimuthally average over to generate a contrast curve. We henceforth refer to this method as the Ruffio method.

The derived 5σ contrast/upper-limit curves, along with the 1σ and 3σ curves, calculated using the Ruffio method are shown in Figure 8. Beyond 110 mas, we reach a 5σ contrast upper limit of >7 mag and a 3σ limit of ~7.6–7.8 mag, which is consistent with the expected photon noise limit of ~7.8, given by Equation (10) in A. Sivaramakrishnan et al. (2023). The contours show the 1σ, 2σ, and 3σ confidence intervals on the contrast and the separation for PDS 70 b and c calculated from the Bayesian average of the three model posteriors.

Zoom In Zoom Out Reset image size

We also calculate a 5σ upper limit on any flux at the reported position of the planet candidate PDS 70 d from V. Christiaens et al. (2024; sep. = 115.2 mas, P.A. = 291$\mathop{.}\limits^{\unicode{x000B0}}$8), as our observations are only separated by 12 days. We first calculate the contrast posterior distribution using dynesty, with uniform priors on the contrast, from −1 to 1. We find a mean contrast of −0.1(±3.0) × 10⁻⁴. To convert this measured contrast to an upper limit, we use Equation (8) from J. B. Ruffio et al. (2018) with the measured mean and standard deviation and a cutoff probability of 99.999971%, which is equivalent to a 5σ upper limit, with a positivity prior on the contrast. From this, we calculate a 5σ upper limit of 7.04 mag, corresponding to a flux of 208 ± 10 μJy. We also calculate the 5σ upper limit using the Laplace approximation, again retrieving an upper limit of 7.04 mag, in line with our assumption of a Gaussian posterior.

4.3. SED Fitting

We fit the available data points of both protoplanets including the new NIRISS/AMI F480M measurement. For planet b, we use the SPHERE/integral field spectrograph (IFS) spectrum and SPHERE/IRDIS photometry from A. Mülleret al. (2018), the VLTI/GRAVITY spectrum from J. J. Wang et al. (2021b), the VLT/NaCo photometry at 3–5 μm from T. Stolker et al. (2020), and the F187N and F480M NIRCam measurements from V. Christiaens et al. (2024). For planet c, we consider the SPHERE/IRDIS and VLT/NaCo photometry from T. Stolker et al. (2020), the SPHERE/IFS and VLTI/GRAVITY from J. J. Wang et al. (2021b), and the F187N and F480M NIRCam measurements from V. Christiaens et al. (2024).

Given the best fit presented by J. J. Wang et al. (2021b), we use the Drift-PHOENIX models (P. Woitke & C. Helling 2003, 2004; C. Helling & P. Woitke 2006; C. Helling et al. 2008) to describe the atmospheric emission from the protoplanets. Following J. J. Wang et al. (2021b), we do not apply any correction for extinction. The parameter space is explored using pymultinest (J. Buchner et al. 2014), employing 1000 live points. For both planets, we first fit the available data without including any CPD contribution. We then include a CPD contribution in the form of a single-temperature blackbody. Blackbody emission has been shown to be the simplest model able to describe well the CPD emission of GQ Lup B up to 11.7 μm (G. Cugno et al. 2024). Although a blackbody is a simple model and is not the only possible solution to characterize CPD emission, more complex models would not provide useful information, as the CPD contribution is mostly determined by the F480M photometry. The model parameters for both PDS 70 b and c are shown in Section A.3.

For planet b, the best fit without a CPD component is reported by the green solid line in Figure 9. The spectrum reproduces well the IFS spectrum, the SPHERE H-band photometry, and the F187N measurement obtained with NIRCam. However, the first half of the GRAVITY spectrum is underestimated, while the second appears to be overestimated. Finally, at λ > 3 μm, almost none of the data points appear to agree within 1σ with the best-fit model, with the exception of the NIRCam F480M photometry, which has a very large uncertainty. When including the CPD contribution, the second half of the GRAVIITY spectrum is better reproduced, while at longer wavelengths, the predicted flux is closer to the measurements (see orange dashed line in Figure 9). In terms of planet model parameters, the presence of mid-IR excess emission in the modeling suggests a hotter and smaller planet. The Bayes factor comparing the blackbody model to the no-blackbody model is 2.8 × 10¹⁰. Hence, the presence of the additional blackbody emission provides a better description of the observations.

Zoom In Zoom Out Reset image size

The results for planet c are displayed in Figure 10. As reported by V. Christiaens et al. (2024), the F187N shows an excess emission, possibly due to Paα emission. The blue region of the GRAVITY spectrum does not agree with the fit, and the $L^{\prime}$ measurement is off by almost 2σ. Also, our NIRISS/AMI F480M measurement is clearly at odds with respect to the green solid line model fit, which does not include a CPD component. The spectrum of the best model once including the CPD presents a lower $L^{\prime}$ flux, consistent with the NaCo measurement. Furthermore, the CPD contribution passes right in between the F480M measurements from NIRCam and NIRISS. The Bayes factor of the model including the blackbody compared with the blackbody-free model is 0.88, indicating that the model without the additional blackbody component is preferred. However, we note that our measurement of PDS 70 c hints that there is in fact excess emission that cannot be explained by the atmosphere-only model, and the fit is being constrained by the fainter measurement of PDS 70 c made by V. Christiaens et al. (2024). Follow-up observations at 4.8 μm with NIRISS and NIRCam and beyond 5 μm with MIRI could be used to confirm the presence of warm CPD emission and constrain the CPD properties.

Zoom In Zoom Out Reset image size

5.1. PDS b and c F480M Photometry

5.2. Contrast Limits and the Nature of the Signal Seen in the Residuals

In this work, we present James Webb Interferometer observations of PDS 70 with the NIRISS F480M filter, the first space-based interferometric observations of this system. Using a joint model fitting approach to simultaneously fit for star, disk, and planet emission, we redetect the protoplanets PDS 70 b and c and derive fluxes of both planets in this filter (Table 1). Additionally, we place the deepest constraints on additional planets within the disk gap of PDS 70 inside ~250 mas in F480M and calculate an F480M upper limit on the flux of the candidate PDS 70 d. We also detect a new feature at an SNR of ~4 to the south of PDS 70 A whose nature is uncertain and will require follow-up observations to confirm.

Furthermore, our results show that NIRISS/AMI can reliably measure relative astrometry and contrasts of young planets in a part of parameter space (small separations and moderate-to-high contrasts) that is unique to this observing mode and inaccessible to all other present facilities at 4.8 μm. We demonstrate a NIRISS/AMI observing strategy for targets faint enough to acquire greater than ~10 groups up the ramp before the signal in the central 9 pixels reaches a total (linearized) value of 30,000 DN (~48,000 electrons). We show that by using these stringent data selection criteria, we achieve nearly photon-noise-limited performance. For observing brighter targets and to overcome the limitation of the sparsity of the NIRISS AMI uv-coverage, methods that are able to analyze the data directly in the image plane (or using the full extent of the Fourier "splodges") will be necessary. A promising approach will be forward modeling of the full optical system and detector systematics (e.g., the BFE, 1/f noise), as is made possible with ∂Lux (L. Desdoigts et al. 2023).

The authors would like to thank the referee for their comments and suggestions that significantly strengthened the paper. The authors also thank Gregory Herczeg for his helpful comments. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program GTO 1242 and can be accessed via doi:10.17909/6qvy-zr60.

D.B. and D.J. acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC). G.C. thanks the Swiss National Science Foundation for financial support under grant numbers P500PT\_206785 and P5R5PT\_225479. M.D.F. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2303911. J.S.-B. acknowledges the support received from the UNAM PAPIIT project IA 105023.

We acknowledge and respect the Lək̓ʷəŋən n (Songhees and Esquimalt) peoples, on whose territory the University of Victoria stands, and the Lək̓ʷəŋən n and W̱SÁNEĆ peoples, whose historical relationships with the land continue to this day.

B.J.S.P. acknowledges the traditional owners of the land on which the University of Queensland is situated, and P.G.T. and L.D. would like to acknowledge the Gadigal People of the Eora nation, upon whose unceded, sovereign, ancestral lands they work, and pay respect to their ancestors and descendants, who continue cultural and spiritual connections to the country.

Software: JAX (J. Bradbury et al. 2018), NumPy (C. R. Harris et al. 2020), dynesty (J. S. Speagle 2020; S. Koposov et al. 2023), Astropy (Astropy Collaboration et al. 2013, 2018, 2022), Matplotlib (J. D. Hunter 2007), corner (D. Foreman-Mackey 2016), pymultinest (J. Buchner et al. 2014), jwst (H. Bushouse et al. 2023).

A.1. Disk Parameters

Table 2 shows the derived disk parameters for each of the Henyey–Greenstein, Gaussian, and power-law models. In all cases, we used Gaussian priors on the disk geometry, with the mean values taken from M. Keppler et al. (2018, 2019) and sensibly small standard deviation values chosen so as to break degeneracies between the disk model parameters using the previously measured geometry of the disk. This was done because we are not primarily interested in independently measuring the parameters of the disk but wish to produce a sensible model of the disk emission so as to measure the emission from the known planets.

The derived disk parameters are nearly identical between the models. This confirms what we see in the correlation plot in Figure 5, that there are no significant correlations between the disk parameters and the planet parameters that could be biasing our results. If strong correlations were present, we would see significant differences between the model disk parameters between the two cases.

Table 2. PDS 70 Disk Model Parameters

Parameters	Prior	Henyey–Greenstein	Gaussian	Power Law
$\mathrm{log}{A}_{a}$ (arb.)	${ \mathcal U }(-10,-4)$	$-4.5{3}_{-0.09}^{+0.09}$	$-5.1{1}_{-0.10}^{+0.12}$	$-5.1{4}_{-0.10}^{+0.11}$
r0 (arcsec)	${ \mathcal U }(0.35,0.65)$	$0.4{4}_{-0.02}^{+0.02}$	$0.4{6}_{-0.03}^{+0.04}$	$0.4{6}_{-0.03}^{+0.03}$
FWHMr (arcsec)	${ \mathcal U }(0.03,1.00)$	$0.2{9}_{-0.06}^{+0.08}$	$0.2{2}_{-0.05}^{+0.06}$	$0.2{2}_{-0.05}^{+0.06}$
i (deg)	${ \mathcal N }(51.7,1.0)$a	$52.{7}_{-1.0}^{+0.9}$	$52.{2}_{-0.9}^{+0.9}$	$52.{3}_{-0.9}^{+0.9}$
P.A. (deg)	${ \mathcal N }(160.4,1.0)$a	$160.{0}_{-0.9}^{+0.8}$	$160.{2}_{-0.8}^{+0.9}$	$160.{4}_{-0.9}^{+0.9}$
$\mathrm{log}{A}_{s}$ (arb.)	${ \mathcal U }(-10,-4)$	$-8.{6}_{-1.0}^{+1.1}$	$-8.{2}_{-1.2}^{+1.2}$	$-8.{2}_{-1.2}^{+1.1}$
H100 (au)	${ \mathcal N }(13,5)$b,c	$9.{9}_{-4.7}^{+5.3}$	$13.{7}_{-4.7}^{+4.9}$	$14.{7}_{-4.8}^{+4.5}$
β	${ \mathcal N }(1.25,0.01)$b	$1.2{5}_{-0.01}^{+0.01}$	$1.2{5}_{-0.01}^{+0.01}$	$1.2{5}_{-0.01}^{+0.01}$
$\mathrm{log}{I}_{o}$ (arb.)	${ \mathcal U }(-10,-2)$	$-6.{7}_{-2.2}^{+3.0}$	$-6.{2}_{-2.6}^{+2.8}$	$-5.{9}_{-2.7}^{+2.6}$
α	${ \mathcal U }(0,10)$	$5.{4}_{-3.5}^{+3.2}$	$2.{5}_{-1.9}^{+4.9}$	$2.{7}_{-2.0}^{+4.4}$
σa,sys	log${ \mathcal U }(0.0001,0.1)$	$0.00{9}_{-0.002}^{+0.003}$	$0.00{6}_{-0.001}^{+0.002}$	$0.00{6}_{-0.001}^{+0.002}$
g	${ \mathcal U }(0,1)$	$0.3{3}_{-0.04}^{+0.04}$	⋯	⋯
σθ (deg)	${ \mathcal U }(60,360)$	⋯	$10{6}_{-9}^{+12}$
N	${ \mathcal U }(0,100)$	⋯	⋯	$6.{6}_{-1.2}^{+1.2}$
$\mathrm{log}{ \mathcal Z }$d	⋯	−31.3 ± 0.1	−25.0 ± 0.1	−24.7 ± 0.1

Notes. The median along with the 16th and 84th percentiles of the marginalized posteriors calculated with nested sampling are reported. ${ \mathcal U }(a,b)$ denotes a uniform prior with lower and upper bounds a and b. ${ \mathcal N }(\mu ,\sigma )$ denotes a Gaussian prior with mean μ and standard deviation σ.

^aReference: M. Keppler et al. (2019). ^bReference: M. Keppler et al. (2018), using the same distance as in the original paper of 113.43 pc (Gaia Collaboration et al. 2018). ^cTruncated to be positive. ^dThe logarithm of the Bayesian evidence for each model.

Download table as:  ASCIITypeset image

A.2. Planet Parameters

In Table 3, we show the planet parameters, calculated using nested sampling, from each of the Henyey–Greenstein, Gaussian, and power-law models.

Table 3. PDS 70 b and c Parameters from Individual Models

	PDS 70 b	PDS 70 b
Henyey–Greenstein
Separation (mas)	$152.{4}_{-10.2}^{+9.0}$	$215.{7}_{-8.9}^{+7.8}$
Position angle (deg)	$130.{2}_{-1.8}^{+1.6}$	$270.{6}_{-1.7}^{+1.8}$
Contrast (Δmag)	$5.8{5}_{-0.07}^{+0.07}$	$6.4{6}_{-0.14}^{+0.16}$
Gaussian
Separation (mas)	$150.{7}_{-10.8}^{+9.4}$	$218.{7}_{-8.0}^{+7.1}$
Position angle (deg)	$131.{1}_{-1.6}^{+1.5}$	$270.{0}_{-1.7}^{+1.7}$
Contrast (Δmag)	$5.8{4}_{-0.07}^{+0.07}$	$6.4{7}_{-0.14}^{+0.17}$
Power Law
Separation (mas)	$150.{3}_{-10.8}^{+9.6}$	$218.{3}_{-8.3}^{+7.2}$
Position angle (deg)	$131.{3}_{-1.5}^{+1.5}$	$270.{0}_{-1.7}^{+1.8}$
Contrast (Δmag)	$5.8{4}_{-0.07}^{+0.07}$	$6.4{8}_{-0.14}^{+0.17}$

Note. We report the median along with the 16th and 84th percentiles of the marginalized posteriors of each model.

Download table as:  ASCIITypeset image

A.3. SED Fit Parameters

Table 4 shows the Drift-PHOENIX model parameters, with and without a blackbody component, to model any excess CPD emission. We also report the Bayes factor, B₁₂, comparing each model with a blackbody component to the model without a blackbody component.

Table 4. Drift-PHOENIX SED Model Parameters for PDS 70 b and c

	Teff	$\mathrm{log}g$	[M/H]	R	Tdisk	Rdisk	B12
	(K)			(RJ)	(K)	(RJ)
PDS 70 b	$133{1}_{-10}^{+15}$	$3.1{4}_{-0.11}^{+0.81}$	$0.0{2}_{-0.16}^{+0.13}$	$2.2{2}_{-0.07}^{+0.06}$	⋯	⋯	⋯
	$140{2}_{-22}^{+27}$	$5.3{7}_{-0.18}^{+0.09}$	$0.0{4}_{-0.20}^{+0.18}$	$1.7{7}_{-0.08}^{+0.08}$	$20{7}_{-9}^{+20}$	$121{5}_{-569}^{+492}$	2.8 × 1010
PDS 70 c	$103{2}_{-21}^{+32}$	$3.2{3}_{-0.16}^{+0.22}$	$0.0{9}_{-0.14}^{+0.12}$	$2.3{5}_{-0.17}^{+0.19}$	⋯	⋯	⋯
	$104{6}_{-28}^{+42}$	$3.2{6}_{-0.18}^{+0.28}$	$-0.0{4}_{-0.20}^{+0.14}$	$2.2{2}_{-0.18}^{+0.20}$	$23{0}_{-62}^{+55}$	$16{7}_{-114}^{+146}$	0.88

Note. The median along with the 16th and 84th percentiles of the marginalized posteriors are reported. B₁₂ is the Bayes factor between the blackbody CPD model and the no-blackbody model.

Download table as:  ASCIITypeset image

In Figure 11, we show the evolution of the pixel intensities for the central 9 pixels of the PSF for PDS 70 and HD 123991. We show how the linearized ramp data evolve (orange) as each pixel accumulates signal, along with how this affects the calibrated rate data (red) output by the jwst pipeline (the calints data). We note that for the central pixel, the detected signal decreases, apparently at all signal levels, with the signal in most of the adjacent pixels either rising or falling, with noticeable differences between PDS 70 and HD 123991. The chosen cutoff value, where we discard all data above this group, is denoted by the vertical dashed red line. For all pixels other than the central pixel, there is no significant signal deficit/excess up to this cutoff value.

Zoom In Zoom Out Reset image size