# The James Webb Interferometer: Space-based interferometric detections of PDS
70 b and c at 4.8 $\mu$m

Dori Blakely Department of Physics and Astronomy, University of Victoria, 3800
Finnerty Road, Elliot Building, Victoria, BC, V8P 5C2, Canada  NRC Herzberg
Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC, V9E 2E7,
Canada  Doug Johnstone NRC Herzberg Astronomy and Astrophysics, 5071 West
Saanich Road, Victoria, BC, V9E 2E7, Canada  Department of Physics and
Astronomy, University of Victoria, 3800 Finnerty Road, Elliot Building,
Victoria, BC, V8P 5C2, Canada  Gabriele Cugno Department of Astronomy,
University of Michigan, Ann Arbor, MI 48109, USA  Anand Sivaramakrishnan Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Astrophysics Department, American Museum of Natural History, 79th Street at
Central Park West, New York, NY 10024  Department of Physics and Astronomy,
Johns Hopkins University, 3701 San Martin Drive, Baltimore, MD 21218, USA
Peter Tuthill Sydney Institute for Astronomy, School of Physics, University of
Sydney, NSW 2006, Australia  Ruobing Dong Department of Physics and Astronomy,
University of Victoria, 3800 Finnerty Road, Elliot Building, Victoria, BC, V8P
5C2, Canada  Benjamin J. S. Pope School of Mathematics and Physics, University
of Queensland, St Lucia, QLD 4072, Australia  University of Southern
Queensland, Centre for Astrophysics, Toowoomba, Queensland, Australia  Loïc
Albert Trottier Institute for Research on Exoplanets, Département de Physique,
Université de Montréal, 1375 Ave Thérèse-Lavoie-Roux, Montréal, QC, H2V 0B3,
Canada  Max Charles Sydney Institute for Astronomy, School of Physics,
University of Sydney, NSW 2006, Australia  Rachel A. Cooper Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Matthew De
Furio Department of Astronomy, The University of Texas at Austin, 2515
Speedway Stop C1400, Austin, TX 78712, USA  Louis Desdoigts Sydney Institute
for Astronomy, School of Physics, University of Sydney, NSW 2006, Australia
René Doyon Trottier Institute for Research on Exoplanets, Département de
Physique, Université de Montréal, 1375 Ave Thérèse-Lavoie-Roux, Montréal, QC,
H2V 0B3, Canada  Observatoire du Mont-Mégantic, Université de Montréal,
Montréal H3C 3J7, Canada  Logan Francis Leiden Observatory, Leiden University,
PO Box 9513, 2300 RA Leiden, The Netherlands  Alexandra Z. Greenbaum IPAC,
Caltech, 1200 E. California Blvd., Pasadena, CA 91125, USA  David Lafrenière
Trottier Institute for Research on Exoplanets, Département de Physique,
Université de Montréal, 1375 Ave Thérèse-Lavoie-Roux, Montréal, QC, H2V 0B3,
Canada  James P. Lloyd Carl Sagan Institute/Department of Astronomy, Cornell
University, Ithaca NY 14853  Michael R. Meyer Department of Astronomy,
University of Michigan, Ann Arbor, MI 48109, USA  Laurent Pueyo Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Shrishmoy Ray School of Mathematics and Physics, University of Queensland, St
Lucia, QLD 4072, Australia  Joel Sánchez-Bermúdez Instituto de Astronomía,
Universidad Nacional Autónoma de México, Apdo. Postal 70264, Ciudad de México,
04510, México  Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117
Heidelberg, Germany  Anthony Soulain Univ. Grenoble Alpes, CNRS, IPAG, 38000
Grenoble, France  Deepashri Thatte Space Telescope Science Institute, 3700 San
Martin Drive, Baltimore, MD 21218, USA  Thomas Vandal Trottier Institute for
Research on Exoplanets, Département de Physique, Université de Montréal, 1375
Ave Thérèse-Lavoie-Roux, Montréal, QC, H2V 0B3, Canada

###### Abstract

We observed the planet-hosting system PDS 70 with the James Webb
Interferometer, JWST's Aperture Masking Interferometric (AMI) mode within
NIRISS. Observing with the F480M filter centered at 4.8 $\mu$m, we
simultaneously fit a geometric model to the outer disk and the two known
planetary companions. We re-detect the protoplanets PDS 70 b and c at an SNR
of 21 and 11, respectively. Our photometry of both PDS 70 b and c provide
evidence for circumplanetary disk emission through fitting SED models to these
new measurements and those found in the literature. We also newly detect
emission within the disk gap at an SNR of $\sim$4, at a position angle of
$207^{+11}_{-10}$ degrees, and an unconstrained separation within $\sim$200
mas. Follow-up observations will be needed to determine the nature of this
emission. We place a 5$\sigma$ upper limit of $\Delta$mag = 7.56 on the
contrast of the candidate PDS 70 d at 4.8 $\mu$m, which indicates that if the
previously observed emission at shorter wavelengths is due to a planet, this
putative planet has a different atmospheric composition than PDS 70 b or c.
Finally, we place upper limits on emission from any additional planets in the
disk gap. We find an azimuthally averaged 5$\sigma$ upper limit of $\Delta$mag
$\approx$ 7.5 at separations greater than 125 mas. These are the deepest
limits to date within $\sim$250 mas at 4.8 $\mu$m and the first space-based
interferometric observations of this system.

††software: JAX (Bradbury et al., 2018), NumPy (Harris et al., 2020), dynesty
(Speagle, 2020; Koposov et al., 2023), Astropy (Astropy Collaboration et al.,
2013, 2018), Matplotlib (Hunter, 2007), corner (Foreman-Mackey, 2016),
pymultinest (Buchner et al., 2014).

##  1 Introduction

PDS 70 is one of the most extensively studied young stellar systems. It is the
only known multi-planet protoplanetary disk system, where two or more planets
have been robustly detected within the disks from which they formed (Keppler
et al., 2018; Haffert et al., 2019). The disk consists of a large outer
component and a smaller inner component (Keppler et al., 2018), separated by a
wide cavity (Hashimoto et al., 2012; Facchini et al., 2021) spanning between
$\lesssim$17 au and $\sim$54 au (Keppler et al., 2018). The inner disk was
first inferred via spectral energy distribution (SED) fitting (Hashimoto et
al., 2012; Dong et al., 2012), and was subsequently isolated in the near-
infrared (Keppler et al., 2018) and sub-millimeter (Keppler et al., 2019).
Recently observations with JWST/MIRI detected water vapour in the inner disk
(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 (Wagner et al., 2023) or AB
Aurigae (Currie et al., 2022). The PDS 70 outer disk appears as a close to
symmetric ring at sub-mm wavelengths (Keppler et al., 2019; Benisty et al.,
2021), with one significant azimuthal asymmetry observed to the north-west of
the star that is close to coincident with an arm-like structure seen in the
near-infrared (Wang et al., 2020; Juillard et al., 2022).

The planets PDS 70 b and c have been extensively imaged in the near-infrared
from $\sim$1 to 3.8 $\mu$m (e.g., Keppler et al., 2018; Müller et al., 2018;
Mesa et al., 2019; Christiaens et al., 2019; Wang et al., 2020, 2021b; Cugno
et al., 2021), with a few detections beyond $\sim$4 $\mu$m to 4.8 $\mu$m
(Stolker et al., 2020a; Christiaens et al., 2024), as well as detections in Hα
(Wagner et al., 2018; Haffert et al., 2019). Circumplanetary disk emission
from PDS 70 c has been directly detected in the sub-mm with ALMA, coincident
with detections in the near-infrared (Isella et al., 2019; Benisty et al.,
2021). Sub-mm emission has also been observed that is seemingly associated
with PDS 70 b, however it is offset from the near-infrared detections (Isella
et al., 2019; Benisty et al., 2021). Additionally, emission at the L5
Lagrangian point of PDS 70 b has tentatively been detected in the sub-mm with
ALMA (Balsalobre-Ruza et al., 2023). Finally, a third point source has been
tentatively detected in the PDS 70 disk cavity, using VLT/SPHERE over multiple
epochs (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 $\sim$1-1.6 $\mu$m, implying it
may be scattered light (Mesa et al., 2019). This source was tentatively re-
detected at 1.9 $\mu$m with JWST/NIRCam (Christiaens et al., 2024).

Many atmospheric models of massive giant planets predict a drop in flux in the
mid-infrared beyond $\sim$4.5 $\mu$m compared to between 2-2.5 $\mu$m and 3-4
$\mu$m, so observations at 4.8 $\mu$m may allow for excess emission from warm
circumplanetary disk (CPD) dust to be inferred. Wang et al. (2021b) showed
that the VLT/NACO M′ detection of PDS 70 b from Stolker et al. (2020a)
provided weak support for excess blackbody-like emission, so follow up
observations at $\sim$4.8 $\mu$m and beyond will allow for this tentative
detection of CPD emission in the mid-infrared to be confirmed. Recently,
Christiaens et al. (2024) re-detected PDS 70 b and c with the JWST/NIRCam (4.8
$\mu$m) F480M filter. Their detection of PDS 70 b was in agreement with the M′
detection, but with a larger uncertainty.

In our work, we re-detect PDS 70 b and c at F480M using the James Webb
Interferometer, an independent method to Christiaens et al.'s JWST/NIRCam
results, and make the most precise measurement of the flux of PDS 70 b at 4.8
$\mu$m, using the JWST/NIRISS F480M filter. We also detect PDS 70 c using the
F480M filter at a similar precision to Christiaens et al. (2024), after
including our conservative systematic noise estimation in our error budget.
Wang et al. (2021b) made the first interferometric detection of PDS 70 b and c
from the ground, using VLTI/GRAVITY. Here, using the power of the James Webb
Interferometer, 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: § 2 outlines the data acquisition,
reduction and cleaning procedure; § 3 details the disk plus two planet model
fitting and the conversion of the measured contrasts to fluxes; § 4 outlines
the derived planet parameters and the analysis of the data for any signal
beyond the two known planets and the disk; § 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 § 6.

##  2 Observations and data reduction

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 February 24 2023
using the James Webb Interferometer, the aperture masking interferometry mode
of the NIRISS instrument (Sivaramakrishnan et al., 2012, 2023; Doyon et al.,
2023) henceforth referred to as AMI, in the F480M filter ($\lambda=4.815$
$\mu$m, $\Delta\lambda=0.298$ $\mu$m)111https://jwst-docs.stsci.edu/jwst-near-
infrared-imager-and-slitless-spectrograph/niriss-instrumentation/niriss-
filters. The SUB80 sub-array 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.

Figure 1: Mean PDS 70 calints file, after bad pixel correction, using only the
first 28 groups.

The data were reduced using version 1.11.3 of the jwst
pipeline222https://github.com/spacetelescope/jwst, from the uncal data format.
We use the default calwebb_detector1 stage 1 pipeline and calwebb_image2 stage
2 pipeline. Note that we do not use the charge_migration333named
undersampling_correction in the 1.11.3 version of the pipeline that we used
step (Goudfrooij et al., 2024), that 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., 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 criteria because the BFE depends on the contrasts between pixels
(e.g., Plazas et al., 2017). Due to the fine structure in the AMI PSF, that is
somewhat undersampled at 4.8 $\mu$m with the NIRISS plate scale of $\sim$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
(Sivaramakrishnan et al., 2023; Kammerer et al., 2023; Ray et al., 2023;
Sallum et al., 2023b). 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$\times$3 pixels of the PSF is
above 30,000 data number (DN), or $\sim$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 $\sim$1%, and minimal bleeding into surrounding
pixels (Appendix B, Figure 10). We note that this loss in signal biases our
measured contrasts by significantly less than 1%, due to the central pixel
accounting for $\sim$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 $\sim$5800 DN. With fewer groups, the
data is likely to be dominated by the 1/f noise seen in the NIRISS detector
(e.g., Sallum et al., 2023b). 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.

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
(Ireland, 2013). We also perform sub-pixel 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 (Kammerer et al.,
2019).

The background star near PDS 70 that was originally observed by Riaud et al.
(2006), and confirmed to be a background star by Hashimoto et al. (2012), is
at the edge of the SUB80 array, at a separation of $\sim$2.5 arc seconds (as)
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 $\sigma$ = 30 pixels, is
applied to both the PDS 70 data and the HD 123991 data.

We calculate the squared visibilities, $V^{2}$, and closure phases, of PDS 70
and calibrator, HD 123991, using amical (Soulain et al., 2020; Soulain &
Robert C. M. T., 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^{2}$ of PDS 70 by the $V^{2}$ of HD 123991. Similarly, the
closure phases are calibrated by subtracting the closure phases of HD 123991
from those of PDS 70.

##  3 Methods

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

$V_{*,b,c}(u,v)=I_{*}+I_{b}e^{-i2\pi(ux_{b}+vy_{b})}+I_{c}e^{-i2\pi(ux_{c}+vy_{c})},$
(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 was
used to succesfully reconstruct the extended emission of LkCa 15, in the near
infrared, with ground based AMI (Blakely et al., 2022). 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 optically thick emission from the disk surface, in a similar
manner to the procedure used by Stolker et al. (2016) and Sissa et al. (2018).
We describe the height of the surface of the disk using a power law profile,
given by

$z=H_{100}\left(\frac{\sqrt{x^{2}+y^{2}}}{100\mathrm{\>au}}\right)^{\beta},$
(2)

where x and y are the potentially rotated and/or inclined coordinates of the
disk mid-plane, in units of au, $H_{100}$ is the aspect ratio at 100 au
($\sim$0.9 arc seconds), and $\beta$ 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,

$r=\sqrt{x^{2}+(y+z\sin(i))^{2}+z^{2}}.$ (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

$I_{d_{1}}(r)=\exp\left(-\frac{(r-r_{0})^{2}}{2\sigma_{r}^{2}}\right),$ (4)

and $I_{d_{2}}(r)$ is

$I_{d_{2}}(r)=\frac{1}{2}\left(1+\mathrm{erf}\left(\alpha\frac{(r-r_{0})}{\sqrt{2}\sigma_{r}}\right)\right),$
(5)

where $r_{0}$ describes the location of the peak of the disk brightness,
$\sigma_{r}$ describes the radial extent of the disk emission and $\alpha$
describes the degree to which the inner edge of the disk is truncated.

The azimuthal brightness distribution is a sum of an axisymmetric component
and an asymmetric component described by the Henyey-Greenstein scattering
phase function (Henyey & Greenstein, 1941)

$I_{d_{3}}(\theta)=\frac{1-g^{2}}{4\pi(1+g^{2}-2g\cos(\theta-\theta_{0}))^{3/2}},$
(6)

where $g$ is the scattering parameter, a free parameter constrained to be
between 0 (isotropic scattering) and 1 (forward scattering), and $\theta_{0}$
is the azimuthal location of the scattering peak, which is fixed to the near
side of the disk.

Combining Equations 2 through 6, the full disk model is given by

$I_{d}(r,\theta)=\left(A_{a}I_{d_{3}}(\theta)+A_{s}\right)I_{d_{1}}(r)I_{d_{2}}(r),$
(7)

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

Finally, any over-resolved 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 7, normalized by the total flux of the model. The complex
visibilities of this 17 parameter model are given by

$V_{m}=\frac{\mathcal{F}\\{I_{d}\\}+V_{*,b,c}}{\Sigma
I_{d}+I_{*}+I_{b}+I_{c}+I_{o}}.$ (8)

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.

We follow the general procedure used by Blakely et al. (2022) and Sallum et
al. (2023a) of fitting geometrical models to the squared visibilities and
closure phases calculated from AMI data. Our log-likelihood is the sum of an
unnormalized Gaussian log-likelihood for both the closure phases and squared
visibilities. For the closure phases, the uncertainties are simply the
statistical uncertainties, calculated across integrations. However, for the
squared visibilities, to account for a measurable drift at each baseline,
observed as the signal level increases, even at low flux levels, an estimated
systematic error budget of 1% is added in quadrature with the statistical
uncertainties. We estimate this 1% systematic on the squared visibilities by
comparing the calculated squared visibilities for PDS 70 and HD 123991 for all
signal levels (we show the evolution of pixel intensities in Appendix B,
Figure 10). We observe a drift from group 1 up to our cutoff signal level
(groups 28 and 10 for PDS 70 and HD 123991, respectively) that is inconsistent
with the calculated statistical uncertainties. We find that the difference in
the evolution of the squared visibilities between PDS 70 and HD 123991 is
$\sim$1% at the cutoff signal. We do not include a similar systematic for the
closure phases because the evolution we see for increasing signal level is
largely consistent with the calculated statistical uncertainities.

We use a Bayesian modelling approach to measure the locations and contrasts of
PDS 70 b and c, as well as our model disk parameters. We estimate our model
posterior using dynamic nested sampling (Higson et al., 2018), with dynesty
(Speagle, 2020; Koposov et al., 2023). We use uniform priors on the contrasts
of both planets, between 0 and 10 $\Delta$mag. We also use uniform priors on
the locations of PDS 70 b and c and then we test this result against a model
that uses tightly constrained Gaussian priors, henceforth referred to as
GRAVITY priors, based on the planets predicted locations from Wang et al.
(2021b), made using VLTI/GRAVITY astrometry. The planet location priors are
given in Table 3, in Appendix A. We use priors on the the disk geometry from
Keppler et al. (2018) and Keppler et al. (2019), that are described in
Appendix A.

We note that there are limitations to our joint model fitting approach due to
its necessarily simple nature. This means that, despite the ability of our
joint-modelling approach to measure the correlations between the planet and
disk parameters, the planet parameter uncertainties derived from our nested
sampling fitting do not capture the full uncertainties in the physical
parameters that we are trying to measure. To provide conservative estimates on
the uncertainties of our derived parameters for PDS 70 b and c, we estimate a
systematic uncertainty on their separations, position angles and contrasts. We
use an injection/recovery approach similar to what is used by Stolker et al.
(2024). Because of the non-linearity of closure phases (especially in the
moderate contrast regime, where the visibility amplitudes are not $\sim$1), we
mimic the injection of point sources into the data by injecting a companion
with a negative flux into the model complex visibilities before we calculate
the model closure phase and squared visibilities. We inject the companion at
the median separation and (negative) contrast of PDS 70 b and c found from our
nested sampling fits at position angles between 0 and 360 degrees, in
increments of 5 degrees. We recover the parameters using nested sampling for
each case, while fixing the disk parameters to the median parameters from the
original nested sampling fit. To estimate the systematic uncertainty on the
planet parameters, we use the root-mean-square deviation of the median
estimated parameters from the true parameters, at each position angle.

The mean parameters for the locations and contrasts of PDS 70 b and c, along
with the standard deviation of the distribution calculated with nested
sampling with the estimated systematic uncertainty added in quadrature are
shown in Table 1. For PDS 70 b, we estimate a systematic in the contrast of
0.08 mags, 3.4 mas in separation and 1.2 degrees in position angle. Only the
contrast systematic has a significant contribution to our total error budget,
compared to the uncertainties estimated with nested sampling (Appendix A,
Table 3). For PDS 70 c, we find a contrast systematic of 0.08 mags, a
separation systematic of 7.0 mas, and a position angle systematic of 0.8
degrees. Both the contrast and separation systematic uncertainties have a non-
negligible contribution to the total error budget, indicating that there is
some level of degeneracy between the disk signal and a companion with the
contrast that we estimate for PDS 70 c that we are able to capture using this
injection/recovery procedure.

Figure 2: Geometrical model fit to the F480M closure phase and squared
visibility data of PDS 70 (images), using uniform priors on the positions of
PDS 70 b and c. With respect to the central source (yellow star), PDS 70 b is
to the south-east while PDS 70 c is directly west. The white contours denote
the 1, 2, 3, 4 and 5 $\sigma$ contours of the marginalized posterior of the
positions of PDS 70 b and c. The two grey dashed circles are centered on the
predicted locations of the planets at the time of the observations (Wang et
al., 2021b) from whereistheplanet (Wang et al., 2021a). We mask the central
region to denote the inner working angle of $\sim$0.5$\lambda/B=$ 94 mas, the
diffraction limit of the data.

To convert the measured contrasts to fluxes, we model the spectral energy
distribution of the unresolved stellar plus inner disk emission using a power
law, assuming that any spectral features over the range of the F480M filter
($\sim$4.6-5.0 $\mu$m) are negligible, using the same approach as is done by
Stolker et al. (2020a). We fit the power law, in log-log space, to 2MASS
(Skrutskie et al., 2006) J, H and Ks as well as WISE (Wright et al., 2010) W1
and W2 magnitudes that we convert to fluxes using species (Stolker et al.,
2020b). We employ nested sampling (Skilling, 2004), using dynesty, to estimate
the posterior of the power-law model. To account for the outer disk flux,
which we measure in our fits to be $\sim$4% of the stellar flux, we scale the
model spectrum by 0.96, finding a stellar magnitude of 7.56 $\pm$ 0.02, in
F480M. To convert the measured contrasts of the planets to physical units we
integrate the model spectrum across the F480M filter transmission profile. To
estimate the uncertainties on the derived fluxes we do this for all
combinations of the planet contrast posterior samples and 500 randomly drawn
samples from the power law model posterior, and calculate confidence intervals
from this distribution.

Table 1: PDS 70 b and c parameters

| Separation | Position Angle | Contrast | Flux | SNR
---|---|---|---|---|---
| (mas) | (∘) | ($\Delta$mag) | (W m-2 $\mu$m-1) |
PDS 70 b | $147.8\pm 8.2$ | $131.5\pm 1.6$ | $5.83\pm 0.10$ | $8.75\pm 0.82\times 10^{-17}$ | 21.1
PDS 70 c | $221.1\pm 8.7$ | $269.6\pm 1.4$ | $6.48\pm 0.13$ | $4.83\pm 0.58\times 10^{-17}$ | 10.6

Notes: We report the mean and the standard deviation of the marginalized
posteriors, with the systematic uncertainties estimated from the
injection/recovery added in quadrature with the nested sampling uncertainties.
The SNR is defined here as the mean of the contrast posterior samples divided
by its standard deviation from nested sampling only (displayed in Appendix A,
Table 3).

##  4 Results

###  4.1 Derived Planet Parameters

Figure 2 shows the two planet plus disk model generated from the median
parameters from the uniform planet location model posterior. The 1, 2, 3, 4
and 5 $\sigma$ contours from the marginalized posterior of the planet
locations are shown by the white contours. A comparison of the model squared
visibilities and closure phases to the data is given in Figure 3, 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 Wang et
al. (2021b), which are shown by the dashed grey circles in Figure 2.

Figure 3: Squared visibility data, black, with the model squared visibilities
plotted in red (top) and the closure phase data, black, with the model closure
phases shown in red (bottom), for the two point source (with uniform priors on
the planet locations) plus disk geometrical model fit to the F480M data. For
both panels, the large red circles denote the values calculated with the
median parameter model, and 100 random posterior samples are denoted by the
small, transparent red points. The noise normalized residuals are shown below
both panels.

Table 1 summarizes the derived planet parameters, from both of the 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 (SNR) ratio, in locations that are
consistent between the uniform prior and tightly constrained Gaussian prior
models. We also find that the measured contrasts of the planets are consistent
for the uniform prior and GRAVITY prior models, with nearly identical
contrasts derived for PDS 70 b, and a small difference for PDS 70 c, at the
$\sim$$1\sigma$ level.

Figure 4 shows a corner plot of the planet parameters from the uniform prior
model, along with a correlation matrix of all of the model parameters. We note
that there are correlations between the parameters of both planets. Notably,
there is a moderate anti-correlation between the separation and the contrast
of PDS 70 c, which may explain the slight discrepancy between the derived
contrast of PDS 70 c between the uniform prior model and the GRAVITY prior
model. From the correlation matrix, 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.

###  4.2 Residual Signal and Contrast Limits

We compute an SNR map, shown in Figure 5, made by dividing the best fit point
source binary contrast to the residual closure phase signal by its associated
uncertainty fitting, using our new python package
drpangloss444https://github.com/benjaminpope/drpangloss. drpangloss provides
similar capabilities to the packages CANDID (Gallenne et al., 2015) and
fouriever (Kammerer et al., 2023) but is written using JAX to provide just-in-
time (JIT) compilation, GPU compatibility and auto-differentiation.

Because closure phases are non-linear, we calculate the complex visibilities
of the uniform prior model, using the median of the marginalized posterior for
all parameters, and add this to our binary model visibilities before
calculating the model closure phases. Each cell in the grid shows the best fit
contrast divided by the uncertainty in the contrast calculated using the
Laplace approximation (Tierney & Kadane, 1986), assuming the uncertainties are
Gaussian. We clearly see that there is residual emission at an SNR of $\sim$4,
to the south-west 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 closure phases and find a contrast of
$7.6^{+0.6}_{-1.3}$ magnitudes, a separation of $100^{+45}_{-42}$ mas and a
position angle of $207^{+11}_{-10}$ degrees. The position angle of the signal
is well constrained, however, there is a strong anti-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$\sigma$ contrast curve, from the
residual closure phase signal of the uniform prior model. We use the approach
described above to fit to the residual closure phase signal, and calculate the
contrast curve using the Absil method (Absil et al., 2011; Gallenne et al.,
2015), as is implemented in drpangloss. We note that the Absil method is a
standard approach to measuring contrast limits for interferometric data,
however, it is based on a model comparison between a single star model and a
binary model, and calculates the contrast at which the binary model deviates
by a certain probability from the single star model (Kammerer et al., 2023).
In this work, we wish to calculate an upper limit on the flux, at each
location, in our field-of-view. To measure Bayesian upper limits, we adapt the
method presented by Ruffio et al. (2018) to interferometric data. We first,
calculate a grid of likelihoods for a binary model at fixed R.A. and Dec.
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
(Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; 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 Dec. point. Finally, we plug the
calculated contrasts, uncertainties and a 99.999971% cutoff probability into
Equation 8 in Ruffio et al. (2018), to calculate the 5$\sigma_{eq}$ upper
limit at each point in our field of view, which azimuthally average over to
generate a contrast curve. We henceforth refer to this method as the Ruffio
method.

The derived 5$\sigma$ contrast/upper limit curves, and their $1\sigma$
uncertainties, calculated using the Absil and Ruffio methods, are shown in
Figure 6. We reach a contrast upper limit of $\sim$7.5 $\Delta$mag, which is
reasonably consistent with the expected photon noise limit of $\sim$7.8, given
by Equation 10 by Sivaramakrishnan et al. (2023). The contours show the 1, 2
and 3 $\sigma$ confidence intervals on the contrast and the separation for PDS
70 b and c calculated from the uniform prior model posterior.

We also calculate a 5$\sigma_{eq}$ upper limit on any flux at the reported
position of the planet candidate PDS 70 d from Christiaens et al. (2024) (sep.
= 115.2 mas, P.A. = 291.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
$-5.3\times 10^{-4}\pm 1.9\times 10^{-4}$. To convert this measured contrast
to an upper limit, we use Equation 8 by Ruffio et al. (2018), with the
measured mean, standard deviation and a cutoff probability of 99.999971%,
which is equivalent to a 5$\sigma$ upper limit, with a positivity prior on the
contrast. From this, we calculate a 5$\sigma_{eq}$ upper limit of 7.56 mag,
corresponding to a flux of $1.76\times 10^{-17}$ W m-2 $\mu$m-1. We also
calculate the 5$\sigma_{eq}$ upper limit using the Laplace approximation and
find an upper limit of 7.57 mag, consistent with the nested sampling result
within $<$1%.

###  4.3 SED Fitting

We fit the the available data points of both protoplanets including the new
NIRISS/AMI F480M measurement. For planet b, we use the SPHERE/IFS spectrum and
SPHERE/IRDIS photometries from Müller et al. (2018), the VLTI/GRAVITY spectrum
from Wang et al. (2021b), the VLT/NaCo photometries at $3-5~{}\mu$m from
Stolker et al. (2020a), and the F187N and F480M NIRCam measurements from
Christiaens et al. (2024). For planet c, we consider the SPHERE/IRDIS and
VLT/NaCo photometries from Stolker et al. (2020a), the SPHERE/IFS and
VLTI/GRAVITY from Wang et al. (2021b), and the F187N and F480M NIRCam
measurements from Christiaens et al. (2024).

Given the best-fit presented by Wang et al. (2021b), we use the Drift-PHOENIX
models (Woitke & Helling, 2003, 2004; Helling & Woitke, 2006; Helling et al.,
2008) to describe the atmospheric emission from the protoplanets. Following
Wang et al. (2021b), we do not apply any correction for extinction. The
parameter space is explored using pymultinest (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 circumplanetary disk
emission of GQ Lup B up to $11.7~{}\mu$m (Cugno et al., 2024). Although a
blackbody is a simple model and it 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.4.

For planet b, the best fit without a CPD component is reported by the green
solid line in Figure 7. The spectrum reproduces well the IFS spectrum, the
SPHERE $H$ band photometries, 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 $\lambda>3~{}\mu$m, almost
none of the datapoints appear to agree within $1\sigma$ 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
7). In terms of planet model parameters, the presence of MIR excess emission
in the modeling suggests a hotter and smaller planet. The reduced chi square
is 0.9 for the extincted Drift-PHOENIX model alone and 0.53 for the extincted
Drift-PHOENIX+blackbody model. Additionally, the Bayes factor comparing the
blackbody model to the no blackbody model is $1.6\times 10^{8}$. Hence, the
presence of the blackbody provides a better description of the observations.

The results for planet c are displayed in Figure 8. As reported by Christiaens
et al. (2024), the F187N shows an excess emission, possibly due to Pa$\alpha$
emission. The blue region of the GRAVITY spectrum does not agree with the fit,
and the $L^{\prime}$ measurement is off by almost $2\sigma$. 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 reduced
chi squared for these fits are 0.55 and 0.52 respectively, and the Bayes
factor in favour of the model including the blackbody is 5.6, suggesting
moderate support for including the blackbody component. Furthermore, assuming
that the NIRISS/AMI long-wavelength measurement is less biased, weighting it
more strongly would obtain a warmer/larger CPD contribution.

Figure 4: Corner plot showing the 1D and 2D marginalized posteriors of the
companion parameters from the joint two point source plus geometrical model
fit to the F480M squared visibilities and closure phases, with uniform priors
on the planet locations. The top right figure shows the correlation matrix of
all of the model parameters.

##  5 Discussion

###  5.1 PDS b and c F480M Photometry

For PDS 70 b, and the uniform prior model, we find a contrast of $4.67\pm
0.43\times 10^{-3}$ ($5.83\pm 0.10$ mags). Similarly, for PDS 70 c, and the
uniform prior model, we find a contrast of $2.58\pm 0.31\times 10^{-3}$
($6.48\pm 0.13$ mags). While our result for PDS 70 b is consistent with the
NIRCam F480M results of $3.44\pm 1.50\times 10^{-3}$ published by Christiaens
et al. (2024), and reasonably consistent with the $M^{\prime}$ contrast of
6.15 $\pm$ 0.27 mags published by Stolker et al. (2020a), the PDS 70 c value
is inconsistent with the NIRCam F480M results of $1.47\pm 0.30\times 10^{-3}$,
at roughly the 3$\sigma$ level.

There are several possible reasons for the PDS 70 c discrepancy between NIRCam
and NIRISS/AMI. It is possible that the contrasts published by Christiaens et
al. are systematically underestimated due to self/over subtraction, as a
result of the small roll angle of $\sim$5 degrees between their two
observations. Another possibility for an underestimated brightness of PDS 70 c
could be due to over subtraction of the spatially coincident disk. Christiaens
et al. use a radiative transfer model based on models published by Portilla-
Revelo et al. (2022) and Portilla-Revelo et al. (2023), that were originally
fit to ALMA data and 1.25 $\mu$m polarized VLT/SPHERE data of PDS 70. These
models were iterated on by fitting for the minimum grain size and settling
parameters to less than half of the forward scattering side of the disk, in
F480M. This tightly constrained fitting approach could lead to biased results
that do not adequately reproduce the disk emission at the location of PDS 70
c. Specifically, an incorrect estimation of the scattering phase function will
significantly impact the PDS 70 c measured contrast.

Separating out the disk emission is less of an issue for the interferometric
data set despite the simple disk model employed, using a symmetric, optically
thick disk component plus the asymmetric Henyey-Greenstein scattering phase
function. The signals of the disk and the planet in Fourier space are well
separated and thus we can disentangle the asymmetric contributions from the
forward scattering peak of the disk, and PDS 70 c, offset by $\sim$20 degrees,
due to the majority of our likelihood being from the closure phases (the
squared visibilities are nearly flat at a signal level of $\sim$0.90, Figure
3). Furthermore, our joint modelling of the planet and disk parameters allows
us to measure the correlations between the parameters of both planets and the
disk, which are shown by the correlation matrix in the top right of Figure 4.
Notably for the contrast of PDS 70 c, there is a moderate correlation with g,
the scattering parameter (more forward scattering leads to less emission), a
strong correlation with the separation of PDS 70 b (larger separation leads to
less emission), a moderate anti-correlation with $\alpha$, the disk inner-edge
truncation parameter and a moderate anti-correlation with its own separation
from the star (plus unresolved inner disk emission).

We note that these correlations, between the contrast of PDS 70 c and the disk
model parameters, are small compared to the correlation between the contrast
of PDS 70 c and the separation of PDS 70 b, which is evident in Figure 4. This
means a larger separation of PDS 70 b would correspond to a fainter contrast
of PDS 70 c. To asses the impact of this correlation, we calculate the
contrast of PDS 70 c using only posterior samples within 3$\sigma$ of the
predicted location of PDS 70 c from the orbit fitting presented in Wang et al.
(2021b). We find a contrast of 6.53 $\pm$ 0.12 mags, consistent with our value
of 6.48 $\pm$ 0.13 mags, calculated using the full posterior. As an additional
test to explore the impact of the scattering phase function profile, we
briefly explored other simple parametric scattering phase functions. We tested
an azimuthal Gaussian profile and the sum of an azimuthal Gaussian profile and
the Henyey-Greenstein scattering phase function. We repeated the model fitting
procedure described in Section 3, using uniform priors on the planet
locations, and found contrasts for PDS 70 b and c that were nearly identical
to the results reported in Table 1 (and Table 3, in Appendix A) for both
models.

As another test, we briefly explored including an inner disk component in our
model, with its geometry fixed to that of the outer disk. We tested symmetric
models as well as models that included skewed emission in the forward
scattering direction. We found that symmetric models did not noticeably
improve the fit. For the skewed models, we found that the fit did somewhat
improve, for a large, bright, extended inner disk model, however, we saw
significant correlations between the inner disk and the planet location
parameters. Such a bright and large inner disk is inconsistent with the size
and mass of the inner disk that has been inferred from SED fitting (e.g.,
Gaidos et al., 2024). Additionally, the planet contrasts that we found from
these tests were all consistent with the planets plus outer disk only model,
well within $1\sigma$ of the contrasts that we report in Table 3. Given these
points, as well as the fact that the inner disk is marginally resolved, if at
all, and that it possibly has a complicated morphology (Mesa et al., 2019;
Casassus & Cárcamo, 2022), not well captured by a geometrical model, we
decided not to further explore an inner disk model component.

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

Figure 5: Point source model SNR map calculated from the closure phase
residuals from the two point sources plus geometrical model fit that used
uniform priors on the location of PDS 70 b and c. The white contours denote
the 1, 2, 3, 4 and 5 $\sigma$ contours from the posterior calculated with
nested sampling. The grey dashed lines denote the predicted locations of PDS
70 b and c at the time of the observations (Wang et al., 2021b, a). The black
star denotes the position of the star. The black contours denote the forward
scattering side of the outer disk in our model. As in Figure 2, the white
dashed circle denotes the diffraction limit of the data. Figure 6: Mean
contrast curves calculated from the residuals of the geometrical model plus
two point source fit, that used uniform priors on the planet locations. The
solid line denotes the 5$\sigma$ contrast curve calculated using the Absil
method. The dashed line denotes the 5$\sigma_{eq}$ upper limit calculated
using the Ruffio method. The coloured regions denote one standard deviation
from the mean. The 1, 2, and 3 $\sigma$ confidence intervals of the planet
contrast/separation, calculated from the nested sampling posterior of the
model fit using uniform priors on the planet locations, are denoted by the
black contours. We note that, for the signal depth of our data, the expected
contrast limit for photon noise dominated data is $\sim$7.8, from by Equation
10 in Sivaramakrishnan et al. (2023).

The residual signal shown in Figure 5, is in a direction significantly offset
from the expected forward scattering direction of the inner disk. This
indicates that what we observe is not due to a simple inner disk structure,
and may hint at a complex inner disk morphology such as a spiral or clumpy
features, as has been suggested by Casassus & Cárcamo (2022), from ALMA
observations in the mm, and by Mesa et al. (2019) from near-infrared
observations with VLT/SPHERE. It is also possible that the asymmetry we
observe is related to the larger separation, $\sim$180 mas, ALMA compact
emission seen by Balsalobre-Ruza et al. (2023), as it is at very similar
position angle. It might even highlight extended emission from an accretion
stream between PDS 70 b and c, tentatively detected by Christiaens et al.
(2024). Another scenario is that the signal we observe is due to an additional
planet interior to the orbit of PDS 70 b. To distinguish between these many
scenarios, follow-up observations, at similar wavelengths, will be necessary
with the James Webb Interferometer. Analysing follow-up data using an approach
similar to Thompson et al. (2023a), using
Octofitter555https://sefffal.github.io/Octofitter.jl/dev/ (Thompson et al.,
2023b), would allow for a robust detection of any orbital motion, even with
low significance detections or non-detections at individual epochs.

Additionally, it is clear from our SNR maps (Figure 5), and our contrast upper
limit of 7.56 mags derived in Section 4, that we do not detect a signal (to
the north-west of PDS 70 A, at a P.A. $\approx$ 291.8∘) consistent with the
point-like feature (PLF) seen by Mesa et al. (2019)/PDS 70 d seen by
Christiaens et al. (2024), detected at shorter wavelengths. This indicates
that if the signal is a planet, it likely traces bright scattered light
emission from a potential planetary envelope instead of indicating another
warm source like both PDS 70 b and c. This scenario would be similar to the
protoplanet candidate HD 169142 b (Hammond et al., 2023).

Our contrast upper limits, shown in Figure 6, are the deepest limits on
additional point source emission within $\sim$250 mas in the disk gap of PDS
70, at 4.8 $\mu$m. This clearly demonstrates the power of the James Webb
Interferometer at probing small angular scales compared with direct imaging
that achieves significantly better contrasts at larger separations, as is seen
by Christiaens et al. (2024).

Figure 7: Spectral fit of PDS70 b using Drift-PHOENIX models with (dashed
orange line) and without (solid green line) contribution from a CPD. Blue
circles represent the IFS and GRAVITY spectra, while the photometric
datapoints are reported with a blue square. Horizontal errorbars represent the
effective width of the filters. The red star shows the new NIRISS F480M
measurement. The thick lines are the spectra obtained from the set of
parameters providing the maximum likelihood, while the thin lines report 100
samples randomly drawn from the posterior distribution. Figure 8: Same as
Figure 7, for PDS 70 c.

##  6 Conclusions

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 re-detect 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 $\sim$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
$\sim$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 $\mu$m. We
demonstrate a NIRISS/AMI observing strategy for targets faint enough to
acquire greater than $\sim$10 groups up the ramp before the signal in the
central 9 pixels reaches a total (linearized) value of 30,000 DN ($\sim$48,000
electrons). We show that by using this 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 modelling of the full optical system and
detector systematics (e.g., the brighter-fatter effect, 1/f noise) as is made
possible with $\partial$Lux (Desdoigts et al., 2023).

##  7 acknowledgements

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). 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ḱwŋn (Songhees and Esquimalt) Peoples on whose
territory the University of Victoria stands, and the Lḱwŋn and WSÁNEĆ Peoples
whose historical relationships with the land continue to this day.

BJSP acknowledges the traditional owners of the land on which the University
of Queensland is situated, and PGT and LD would like to acknowledge the
Gadigal People of the Eora nation: upon whose unceded, sovereign, ancestral
lands they works, and pays respect to their Ancestors and descendants, who
continue cultural and spiritual connections to Country.

## References

  * Absil et al. (2011) Absil, O., Le Bouquin, J. B., Berger, J. P., et al. 2011, A&A, 535, A68, doi: 10.1051/0004-6361/201117719
  * Argyriou et al. (2023) Argyriou, I., Lage, C., Rieke, G. H., et al. 2023, A&A, 680, A96, doi: 10.1051/0004-6361/202346490
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Balsalobre-Ruza et al. (2023) Balsalobre-Ruza, O., de Gregorio-Monsalvo, I., Lillo-Box, J., et al. 2023, A&A, 675, A172, doi: 10.1051/0004-6361/202346493
  * Benisty et al. (2021) Benisty, M., Bae, J., Facchini, S., et al. 2021, ApJ, 916, L2, doi: 10.3847/2041-8213/ac0f83
  * Blakely et al. (2022) Blakely, D., Francis, L., Johnstone, D., et al. 2022, ApJ, 931, 3, doi: 10.3847/1538-4357/ac6586
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs, 0.3.13.  http://github.com/google/jax
  * Broyden (1970) Broyden, C. G. 1970, IMA Journal of Applied Mathematics, 6, 76, doi: 10.1093/imamat/6.1.76
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125, doi: 10.1051/0004-6361/201322971
  * Casassus & Cárcamo (2022) Casassus, S., & Cárcamo, M. 2022, MNRAS, 513, 5790, doi: 10.1093/mnras/stac1285
  * Christiaens et al. (2019) Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJ, 877, L33, doi: 10.3847/2041-8213/ab212b
  * Christiaens et al. (2024) Christiaens, V., Samland, M., Henning, T., et al. 2024, arXiv e-prints, arXiv:2403.04855, doi: 10.48550/arXiv.2403.04855
  * Cugno et al. (2024) Cugno, G., Patapis, P., Banzatti, A., et al. 2024, arXiv e-prints, arXiv:2404.07086, doi: 10.48550/arXiv.2404.07086
  * Cugno et al. (2021) Cugno, G., Patapis, P., Stolker, T., et al. 2021, A&A, 653, A12, doi: 10.1051/0004-6361/202140632
  * Currie et al. (2022) Currie, T., Lawson, K., Schneider, G., et al. 2022, Nature Astronomy, 6, 751, doi: 10.1038/s41550-022-01634-x
  * Desdoigts et al. (2023) Desdoigts, L., Pope, B. J. S., Dennis, J., & Tuthill, P. G. 2023, Journal of Astronomical Telescopes, Instruments, and Systems, 9, 028007 .  https://api.semanticscholar.org/CorpusID:259297527
  * Dong et al. (2012) Dong, R., Hashimoto, J., Rafikov, R., et al. 2012, ApJ, 760, 111, doi: 10.1088/0004-637X/760/2/111
  * Doyon et al. (2023) Doyon, R., Willott, C. J., Hutchings, J. B., et al. 2023, PASP, 135, 098001, doi: 10.1088/1538-3873/acd41b
  * Facchini et al. (2021) Facchini, S., Teague, R., Bae, J., et al. 2021, AJ, 162, 99, doi: 10.3847/1538-3881/abf0a4
  * Fletcher (1970) Fletcher, R. 1970, The Computer Journal, 13, 317, doi: 10.1093/comjnl/13.3.317
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaidos et al. (2024) Gaidos, E., Thanathibodee, T., Hoffman, A., et al. 2024, arXiv e-prints, arXiv:2403.09970, doi: 10.48550/arXiv.2403.09970
  * Gallenne et al. (2015) Gallenne, A., Mérand, A., Kervella, P., et al. 2015, A&A, 579, A68, doi: 10.1051/0004-6361/201525917
  * Goldfarb (1970) Goldfarb, D. 1970, Mathematics of Computation, 24, 23.  http://www.jstor.org/stable/2004873
  * Goudfrooij et al. (2024) Goudfrooij, P., Grumm, D., Volk, K., & Bushouse, H. 2024, PASP, 136, 014503, doi: 10.1088/1538-3873/ad1c98
  * Haffert et al. (2019) Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nature Astronomy, 3, 749, doi: 10.1038/s41550-019-0780-5
  * Hammond et al. (2023) Hammond, I., Christiaens, V., Price, D. J., et al. 2023, MNRAS, 522, L51, doi: 10.1093/mnrasl/slad027
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hashimoto et al. (2012) Hashimoto, J., Dong, R., Kudo, T., et al. 2012, ApJ, 758, L19, doi: 10.1088/2041-8205/758/1/L19
  * Helling et al. (2008) Helling, C., Dehn, M., Woitke, P., & Hauschildt, P. H. 2008, ApJ, 675, L105, doi: 10.1086/533462
  * Helling & Woitke (2006) Helling, C., & Woitke, P. 2006, A&A, 455, 325, doi: 10.1051/0004-6361:20054598
  * Henyey & Greenstein (1941) Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70, doi: 10.1086/144246
  * Higson et al. (2018) Higson, E., Handley, W., Hobson, M., & Lasenby, A. 2018, Statistics and Computing, 29, 891–913, doi: 10.1007/s11222-018-9844-0
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ireland (2013) Ireland, M. J. 2013, MNRAS, 433, 1718, doi: 10.1093/mnras/stt859
  * Isella et al. (2019) Isella, A., Benisty, M., Teague, R., et al. 2019, ApJ, 879, L25, doi: 10.3847/2041-8213/ab2a12
  * Juillard et al. (2022) Juillard, S., Christiaens, V., & Absil, O. 2022, A&A, 668, A125, doi: 10.1051/0004-6361/202244402
  * Kammerer et al. (2019) Kammerer, J., Ireland, M. J., Martinache, F., & Girard, J. H. 2019, MNRAS, 486, 639, doi: 10.1093/mnras/stz882
  * Kammerer et al. (2023) Kammerer, J., Cooper, R. A., Vandal, T., et al. 2023, PASP, 135, 014502, doi: 10.1088/1538-3873/ac9a74
  * Keppler et al. (2018) Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44, doi: 10.1051/0004-6361/201832957
  * Keppler et al. (2019) Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118, doi: 10.1051/0004-6361/201935034
  * Koposov et al. (2023) Koposov, S., Speagle, J., Barbary, K., et al. 2023, joshspeagle/dynesty: v2.1.0, v2.1.0, Zenodo, doi: 10.5281/zenodo.7600689
  * Lucy (1974) Lucy, L. B. 1974, AJ, 79, 745, doi: 10.1086/111605
  * Martinache (2010) Martinache, F. 2010, ApJ, 724, 464, doi: 10.1088/0004-637X/724/1/464
  * Martinache (2013) —. 2013, PASP, 125, 422, doi: 10.1086/670670
  * Martinache et al. (2020) Martinache, F., Ceau, A., Laugier, R., et al. 2020, A&A, 636, A72, doi: 10.1051/0004-6361/201936981
  * Mesa et al. (2019) Mesa, D., Keppler, M., Cantalloube, F., et al. 2019, A&A, 632, A25, doi: 10.1051/0004-6361/201936764
  * Müller et al. (2018) Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2, doi: 10.1051/0004-6361/201833584
  * Perotti et al. (2023) Perotti, G., Christiaens, V., Henning, T., et al. 2023, Nature, 620, 516, doi: 10.1038/s41586-023-06317-9
  * Plazas et al. (2017) Plazas, A. A., Shapiro, C., Smith, R., Rhodes, J., & Huff, E. 2017, Journal of Instrumentation, 12, C04009, doi: 10.1088/1748-0221/12/04/C04009
  * Portilla-Revelo et al. (2022) Portilla-Revelo, B., Kamp, I., Rab, C., et al. 2022, A&A, 658, A89, doi: 10.1051/0004-6361/202141764
  * Portilla-Revelo et al. (2023) Portilla-Revelo, B., Kamp, I., Facchini, S., et al. 2023, A&A, 677, A76, doi: 10.1051/0004-6361/202346607
  * Ray et al. (2023) Ray, S., Sallum, S., Hinkley, S., et al. 2023, arXiv e-prints, arXiv:2310.11508, doi: 10.48550/arXiv.2310.11508
  * Riaud et al. (2006) Riaud, P., Mawet, D., Absil, O., et al. 2006, A&A, 458, 317, doi: 10.1051/0004-6361:20065232
  * Richardson (1972) Richardson, W. H. 1972, Journal of the Optical Society of America (1917-1983), 62, 55
  * Ruffio et al. (2018) Ruffio, J.-B., Mawet, D., Czekala, I., et al. 2018, AJ, 156, 196, doi: 10.3847/1538-3881/aade95
  * Sallum et al. (2023a) Sallum, S., Eisner, J., Skemer, A., & Murray-Clay, R. 2023a, ApJ, 953, 55, doi: 10.3847/1538-4357/ace16c
  * Sallum et al. (2023b) Sallum, S., Ray, S., Kammerer, J., et al. 2023b, arXiv e-prints, arXiv:2310.11499, doi: 10.48550/arXiv.2310.11499
  * Shanno (1970) Shanno, D. F. 1970, Mathematics of Computation, 24, 647.  http://www.jstor.org/stable/2004840
  * Sissa et al. (2018) Sissa, E., Gratton, R., Garufi, A., et al. 2018, A&A, 619, A160, doi: 10.1051/0004-6361/201732332
  * Sivaramakrishnan et al. (2012) Sivaramakrishnan, A., Lafrenière, D., Ford, K. E. S., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, ed. M. C. Clampin, G. G. Fazio, H. A. MacEwen, & J. Oschmann, Jacobus M., 84422S, doi: 10.1117/12.925565
  * Sivaramakrishnan et al. (2023) Sivaramakrishnan, A., Tuthill, P., Lloyd, J. P., et al. 2023, PASP, 135, 015003, doi: 10.1088/1538-3873/acaebd
  * Skilling (2004) Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405, doi: 10.1063/1.1835238
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, The Astronomical Journal, 131, 1163, doi: 10.1086/498708
  * Soulain & Robert C. M. T. (2023) Soulain, A., & Robert C. M. T. 2023, AMICAL: Aperture Masking Interferometry Calibration and Analysis Library, Astrophysics Source Code Library, record ascl:2302.021
  * Soulain et al. (2020) Soulain, A., Sivaramakrishnan, A., Tuthill, P., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11446, Optical and Infrared Interferometry and Imaging VII, ed. P. G. Tuthill, A. Mérand, & S. Sallum, 1144611, doi: 10.1117/12.2560804
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Stolker et al. (2016) Stolker, T., Dominik, C., Min, M., et al. 2016, A&A, 596, A70, doi: 10.1051/0004-6361/201629098
  * Stolker et al. (2020a) Stolker, T., Marleau, G. D., Cugno, G., et al. 2020a, A&A, 644, A13, doi: 10.1051/0004-6361/202038878
  * Stolker et al. (2020b) Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020b, A&A, 635, A182, doi: 10.1051/0004-6361/201937159
  * Stolker et al. (2024) Stolker, T., Kammerer, J., Benisty, M., et al. 2024, A&A, 682, A101, doi: 10.1051/0004-6361/202347291
  * Thompson et al. (2023a) Thompson, W., Marois, C., Do Ó, C. R., et al. 2023a, AJ, 165, 29, doi: 10.3847/1538-3881/aca1af
  * Thompson et al. (2023b) Thompson, W., Lawrence, J., Blakely, D., et al. 2023b, AJ, 166, 164, doi: 10.3847/1538-3881/acf5cc
  * Tierney & Kadane (1986) Tierney, L., & Kadane, J. B. 1986, Journal of the American Statistical Association, 81, 82, doi: 10.1080/01621459.1986.10478240
  * Wagner et al. (2018) Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJ, 863, L8, doi: 10.3847/2041-8213/aad695
  * Wagner et al. (2023) Wagner, K., Stone, J., Skemer, A., et al. 2023, Nature Astronomy, 7, 1208, doi: 10.1038/s41550-023-02028-3
  * Wang et al. (2021a) Wang, J. J., Kulikauskas, M., & Blunt, S. 2021a, whereistheplanet: Predicting positions of directly imaged companions, Astrophysics Source Code Library, record ascl:2101.003
  * Wang et al. (2020) Wang, J. J., Ginzburg, S., Ren, B., et al. 2020, AJ, 159, 263, doi: 10.3847/1538-3881/ab8aef
  * Wang et al. (2021b) Wang, J. J., Vigan, A., Lacour, S., et al. 2021b, AJ, 161, 148, doi: 10.3847/1538-3881/abdb2d
  * Woitke & Helling (2003) Woitke, P., & Helling, C. 2003, A&A, 399, 297, doi: 10.1051/0004-6361:20021734
  * Woitke & Helling (2004) —. 2004, A&A, 414, 335, doi: 10.1051/0004-6361:20031605
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, The Astronomical Journal, 140, 1868, doi: 10.1088/0004-6256/140/6/1868

##  Appendix A Model Fitting Results

###  A.1 Disk Parameters

Table 2 shows the derived disk parameters for the model that used uniform
priors on the planet locations, as well as the model that used priors from
Wang et al. (2021b). In both cases, we used Gaussian priors on the disk
geometry, with the mean values taken from Keppler et al. (2018) and Keppler et
al. (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 using the
two sets of planet location priors. This confirms what we see in the
correlation plot in Figure 4, that there are no significant correlations
between the disk parameters and the planet parameters that are biasing our
results. If strong correlations were present we would see significant
differences between the model disk parameters between the two cases.

###  A.2 Planet Parameters

In Table 3 we show the planet parameters, calculated using nested sampling,
for both the model that used uniform priors on the locations of PDS 70 b and
c, and the model that used the GRAVITY priors on the planet locations. The
estimated systematic uncertainties in the contrast, separation and position
angle, calculated in Section 3, are not included in the errors shown in Table
3.

###  A.3 Image Plane Residuals

In Figure 9 we plot the deconvolved image plane residuals from subtracting off
different components of our PDS 70, disk, and planets model fit to the squared
visibilities and closure phases. The deconvolved image plane residuals are
primarily shown to verify the fitting to the interferometric observables and
are not directly comparable to PSF subtracted direct images, due to the
possibility of complex residual features introduced by the highly structured
AMI PSF and the deconvolution process.

As also noted in the main text, we tentatively observe compact emission at a
similar position to the mm emission seen by Balsalobre-Ruza et al. (2023) (top
left panel of Figure 9.). The tentative compact emission seen in the image
plane residuals is at the same position angle and is consistent with the
signal that we see in the residual closure phases, shown in Figure 5. We note
that the deconvolved image plane residuals are a poor diagnostic for measuring
the noise in AMI data, because much of the residual signal is caused by
differential piston errors (the dominant phase error term on the time scale of
hours) to which the closure phase is immune. Thus, the SNR plots shown in
Figure 5 give a much clearer picture of the level of signal left in the data
after model removal of the disk and two planets, which gives us confidence in
the detection of the compact emission to the south of PDS 70 A, as this cannot
be caused by differential piston errors.

###  A.4 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,
B12, comparing each model with a blackbody component to the model without a
blackbody component.

Figure 9: Deconvolved image plane residuals from subtracting the uniform
planet prior model from the mean PDS 70 calints file. The data, and model
image are smoothed by convolving it by a Gaussian kernel with $\sigma$ = 70
mas to eliminate shifting artifacts and minimize photon noise. We then
deconvolve the residuals with the smoothed calibrator PSF using 150 iterations
of the Richardson-Lucy deconvolution algorithm (Richardson, 1972; Lucy, 1974).
Each panel is shown in arbitrary linear units, with the top two panels and the
bottom two panels each shown on the same scale. The top left shows the full
residuals, the top right shows the residuals from subtracting off the star and
disk, the bottom left shows the residuals from subtracting off the star and
both planets and the bottom right shows the residuals from subtracting off the
star.  Table 2: PDS 70 disk model parameters

Parameters | Prior | Uniform prior | GRAVITY prior
---|---|---|---
$\log{A_{a}}$ (arb.) | $\mathcal{U}(-10,-4)$ | $-4.34^{+0.06}_{-0.09}$ | $-4.34^{+0.07}_{-0.12}$
$r_{0}$ (as) | $\mathcal{U}(0.35,0.65)$ | $0.50^{+0.03}_{-0.03}$ | $0.49^{+0.03}_{-0.03}$
FWHMr (as) | $\mathcal{U}(0.03,1.00)$ | $0.22^{+0.04}_{-0.03}$ | $0.22^{+0.06}_{-0.03}$
$i$ (∘) | $\mathcal{N}(51.7,1.0)^{1}$ | $53.9^{+0.7}_{-0.8}$ | $53.9^{+0.8}_{-0.8}$
P.A. (∘) | $\mathcal{N}(160.4,1.0)^{1}$ | $160.2^{+0.7}_{-0.7}$ | $160.3^{+0.7}_{-0.6}$
$\log{A_{s}}$ (arb.) | $\mathcal{U}(-10,-4)$ | $-8.4^{+1.1}_{-1.1}$ | $-8.4^{+1.0}_{-1.0}$
$H_{100}$ (au) | $\mathcal{N}(13,5)^{2,3}$ | $11.8^{+3.3}_{-3.3}$ | $12.2^{+2.9}_{-3.2}$
$\beta$ | $\mathcal{N}(1.25,0.01)^{2}$ | $1.25^{+0.01}_{-0.01}$ | $1.25^{+0.01}_{-0.01}$
$\log{I_{o}}$ (arb.) | $\mathcal{U}(-10,-2)$ | $-5.6^{+2.3}_{-2.8}$ | $-6.2^{+2.6}_{-2.6}$
$g$ | $\mathcal{U}(0,1)$ | $0.29^{+0.02}_{-0.01}$ | $0.30^{+0.01}_{-0.01}$
$\alpha$ | $\mathcal{U}(0,10)$ | $0.7^{+0.9}_{-0.5}$ | $0.8^{+1.6}_{-0.6}$

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, respectively.
$\mathcal{N}(\mu,\sigma)$ denotes a Gaussian prior, with mean $\mu$ and
standard deviation $\sigma$.

1 Reference: Keppler et al. (2019).

2 Reference: Keppler et al. (2018), using the same distance as in the original
paper of 113.43 pc (Gaia Collaboration et al., 2018).

3 Truncated to be positive.

Table 3: PDS 70 b and c parameters from nested sampling

| Sep. Prior | P.A. Prior | Separation | Position Angle | Contrast | Flux | SNR |  $\log{\mathcal{L}}^{*}$2
---|---|---|---|---|---|---|---|---
| (mas) | (∘) | (mas) | (∘) | ($\Delta$mag) | (W m-2 $\mu$m-1) |  |
PDS 70 b | $\mathcal{U}(100,250)$ | $\mathcal{U}(70,160)$ | $148.2^{+6.9}_{-7.8}$ | $131.5^{+1.1}_{-1.1}$ | $5.83^{+0.05}_{-0.05}$ | $8.73^{+0.44}_{-0.40}\times 10^{-17}$ | 21.1 | -30.0
| $\mathcal{N}(155.5,1.4)^{1}$ | $\mathcal{N}(132.7,0.4)^{1}$ | $155.0^{+1.4}_{-1.2}$ | $132.3^{+0.4}_{-0.4}$ | $5.82^{+0.04}_{-0.04}$ | $8.85^{+0.35}_{-0.34}\times 10^{-17}$ | 27.9 | -32.9
PDS 70 c | $\mathcal{U}(100,250)$ | $\mathcal{U}(230,320)$ | $221.3^{+4.8}_{-5.3}$ | $269.7^{+1.2}_{-1.2}$ | $6.47^{+0.10}_{-0.10}$ | $4.82^{+0.46}_{-0.44}\times 10^{-17}$ | 10.6 | $-$
| $\mathcal{N}(209.8,1.0)^{1}$ | $\mathcal{N}(270.1,0.3)^{1}$ | $210.4^{+0.9}_{-0.9}$ | $270.1^{+0.3}_{-0.3}$ | $6.61^{+0.09}_{-0.09}$ | $4.24^{+0.38}_{-0.34}\times 10^{-17}$ | 11.9 | $-$

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

1 From the orbit posterior published in (Wang et al., 2021b), calculated using
whereistheplanet (Wang et al., 2021a).

2 $\log{\mathcal{L}}^{*}$ is the maximum likelihood of the full model
including the star, both planets and the disk.

Table 4: Drift-PHOENIX SED model parameters for PDS 70 b and c

| Teff | $\log g$ | [M/H] | R | Tdisk | Rdisk | B12
---|---|---|---|---|---|---|---
| (K) |  |  | (RJ) | (K) | (RJ) |
PDS 70 b | $1339^{+26}_{-14}$ | $3.83^{+0.31}_{-0.78}$ | $0.03^{+0.18}_{-0.23}$ | $2.16^{+0.08}_{-0.14}$ | $-$ | $-$ | $-$
| $1404^{+26}_{-23}$ | $5.36^{+0.10}_{-0.18}$ | $0.04^{+0.19}_{-0.20}$ | $1.77^{+0.08}_{-0.07}$ | $207^{+19}_{-8}$ | $1267^{+449}_{-568}$ | $1.6\times 10^{8}$
PDS 70 c | $1031^{+26}_{-18}$ | $3.23^{+0.22}_{-0.16}$ | $0.13^{+0.09}_{-0.15}$ | $2.41^{+0.18}_{-0.17}$ | $-$ | $-$ | $-$
| $1050^{+43}_{-32}$ | $3.26^{+0.30}_{-0.18}$ | $-0.00^{+0.18}_{-0.18}$ | $2.19^{+0.18}_{-0.19}$ | $244^{+60}_{-22}$ | $179^{+149}_{-128}$ | 5.6

Notes: The median along with the 16th and 84th percentiles of the marginalized
posteriors are reported. B12 is the Bayes factor between the blackbody CPD
model and the no blackbody model.

##  Appendix B The Brighter-Fatter Effect

In Figure 10, 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 evolves (blue), 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.

Figure 10: Evolution of the central 3$\times$3 pixels as a function of the
intensity in each pixel (in the linearized-ramp level data). PDS 70 is shown
on the left and HD 123991 is shown on the right. The solid red line shows the
evolution of the mean of the calints frames, after data cleaning. Here we are
plotting the calculated rate minus the rate calculated from the first 28
groups and the first 10 groups, for PDS 70 and HD 123991, respectively. These
rates are divided by the rate of the ramp fit using all of the groups. The
dashed blue line shows the evolution of the linearized ramp-level data (ramp-
level immediately after it has undergone the linearity step, in stage 1 of the
jwst pipeline). To clearly illustrate the change in rate in each pixel as a
function of well-depth, we are plotting the difference between groups
separated by 15 groups and 5 groups for PDS 70 and HD 123991, respectively,
for increasing signal. From this, we additionally subtract the difference
between the 25th and 10th group for PDS 70 and the difference between the 10th
group and the 5th group for HD 123991 and then normalize this quantity by the
difference between the final group and one 15 and 5 groups earlier for PDS 70
and HD 123991, respectively. The vertical dashed red line shows the extent of
the data that were used in our analysis.