key: cord-0318041-04celmw2
authors: Benisty, Myriam; Bae, Jaehan; Facchini, Stefano; Keppler, Miriam; Teague, Richard; Isella, Andrea; Kurtovic, Nicolas T.; Perez, Laura M.; Sierra, Anibal; Andrews, Sean M.; Carpenter, John; Czekala, Ian; Dominik, Carsten; Henning, Thomas; Menard, Francois; Pinilla, Paola; Zurlo, Alice
title: A Circumplanetary Disk Around PDS70c
date: 2021-08-16
journal: nan
DOI: 10.3847/2041-8213/ac0f83
sha: 04a66a4657f2aa8e6394187c1177c836932adbc7
doc_id: 318041
cord_uid: 04celmw2

PDS70 is a unique system in which two protoplanets, PDS70b and c, have been discovered within the dust-depleted cavity of their disk, at $sim$22 and 34au respectively, by direct imaging at infrared wavelengths. Subsequent detection of the planets in the H$alpha$ line indicates that they are still accreting material through circumplanetary disks. In this Letter, we present new Atacama Large Millimeter/submillimeter Array (ALMA) observations of the dust continuum emission at 855$mu$m at high angular resolution ($sim$20mas, 2.3au) that aim to resolve the circumplanetary disks and constrain their dust masses. Our observations confirm the presence of a compact source of emission co-located with PDS70c, spatially separated from the circumstellar disk and less extended than $sim$1.2au in radius, a value close to the expected truncation radius of the cicumplanetary disk at a third of the Hill radius. The emission around PDS70c has a peak intensity of $sim$86$pm$16 $mu mathrm{Jy}~mathrm{beam}^{-1}$ which corresponds to a dust mass of $sim$0.031M$_{oplus}$ or $sim$0.007M$_{oplus}$, assuming that it is only constituted of 1 $mu$m or 1 mm sized grains, respectively. We also detect extended, low surface brightness continuum emission within the cavity near PDS70b. We observe an optically thin inner disk within 18au of the star with an emission that could result from small micron-sized grains transported from the outer disk through the orbits of b and c. In addition, we find that the outer disk resolves into a narrow and bright ring with a faint inner shoulder.

Recent surveys have revealed that almost ubiquitously, protoplanetary disks appear highly structured with rings and gaps, spiral arms and asymmetries (e.g., Garufi et al. 2018; Andrews 2020) . While other scenarios are discussed, these features are often interpreted as resulting from the presence of planets embedded in disks (e.g., Dong et al. 2015; Bae & Zhu 2018a; Lodato et al. 2019) . Additional observational support for such a scenario can be found in the form of local perturbation of the gas velocity field from Keplerian rotation (Pinte et al. 2018; Teague et al. 2019; . The quest to detect protoplanets embedded in their host disk through direct imaging has been challenging, with current detection limits on the order of a few Jupiter masses (M Jup ) at large radii (e.g., Huélamo et al. 2018; Asensio-Torres et al. 2021) . A few protoplanet candidates have been claimed in the infrared (IR) and in the H α line (e.g., Sallum et al. 2015; Reggiani et al. 2018) but remain controversial (Mendigutía et al. 2018 ).

The first robust detection through direct imaging techniques of a protoplanet still embedded in its natal disk was obtained in the young system PDS 70 (spectral type K7; M∼0.8 M ; age∼5.4 Myr old; Müller et al. 2018 ) located at ∼112.4 pc (Gaia Collaboration et al. 2020) in the Upper Centaurus Lupus association (Pecaut & Mamajek 2016) . PDS 70 b was discovered with an orbital radius of ∼22 au, and imaged at multiple IR wavelengths Müller et al. 2018) as well as in a filter centered on the Hα line (Wagner et al. 2018a ). Subsequently, PDS 70 c was discovered in Hα imaging at the outer edge of the cavity with an orbital radius of ∼34 au ). These two planets carve a large cavity in the disk, evidenced by a cavity in dust (e.g., Hashimoto et al. 2012; Dong et al. 2012) and a gap in the 12 CO gas emission along the orbit of PDS 70 b ) that indicates significant gas depletion. Observations and hydrodynamic simulations indicate that the planets' orbital configuration is stable, close to a 2:1 mean motion resonance, with PDS 70 b in a slightly eccentric orbit (e∼0.2; Bae et al. 2019; Toci et al. 2020; Wang et al. 2021) . The masses of the two planets are still uncertain, although both planets are likely lighter than 10 M Jup to ensure dynamical stability ) and a non-eccentric outer disk . Spectro-photometric analyses, limited to the IR regime (1-5 µm) remain inconclusive, but suggest planet masses between 1 and a few M Jup * NASA Hubble Fellowship Program Sagan Fellow (e.g., Müller et al. 2018; Mesa et al. 2019; Stolker et al. 2020 ) as well as a clear contribution from dust grains in clouds and/or circumplanetary disks (CPDs) (Christiaens et al. 2019; Stolker et al. 2020; Wang et al. 2020) .

CPDs play a fundamental role in planet formation, as they regulate the gas accretion onto the planet and determine the conditions for satellite formation. As gas enters the planet's sphere of influence, it falls at supersonic velocities onto the surface of the CPD (Tanigawa et al. 2012; Szulágyi & Mordasini 2017) , possibly episodically (Gressel et al. 2013) , leading to shocks that can ionize hydrogen and be traced in the Hα line. From observations of the Hα line, PDS 70 b and PDS 70 c are found to be accreting material from their host disk at a rate of ∼ 10 −8 M Jup per year (Wagner et al. 2018b; Haffert et al. 2019; Thanathibodee et al. 2019; Aoyama & Ikoma 2019; Hashimoto et al. 2020 ). Using Atacama Large Millimeter/ submillimeter Array (ALMA) observations at ∼67 mas×50 mas resolution, Isella et al. (2019) showed evidence for sub-millimeter continuum emission co-located with PDS 70 c, interpreted as tracing a dusty CPD, and for another compact continuum emission source located at ∼74 mas offset in a South West direction from b. The emission around c however was not spatially separated from the outer ring. In this Letter, we present new ALMA observations with 20 mas resolution that provide an independent detection of a compact source of emission colocated with PDS 70 c and of low surface brightness emission within the cavity close to PDS 70 b. The Letter is organized as follows. Section 2 presents the observations and the procedure to calibrate the data. Section 3 presents our new images and analysis. Finally, we discuss our findings in Section 4.

This Letter presents new ALMA observations, hereafter referred to as LB19 (for 'Long Baselines 2019'), obtained in Band 7 (λ = 855 µm), under a Director's Discretionary Time (DDT) program with ID 2018.A.00030.S. PDS 70 was observed during 4 execution blocks (EB) with the C-8 configuration on 2019 July 27, 28, and 30, for a total on-source time of 43 minutes per execution. An observing log including the precipitable water vapor (PWV) levels and calibrator names is given in Appendix A.1. The spectral set-up was tuned to optimize continuum detection, but includes the 12 CO J=3-2 line at 345.8 GHz and the HCO + J=4-3 line at 356.7 GHz, which will be presented in forthcoming papers. The raw data calibration was done with the CASA v.5.6.1 pipeline (McMullin et al. 2007 ) and the self-calibration and post-processing imaging were done using CASA v.5.4.0. We first flagged the channels that included the 12 CO and the HCO + lines and spectrally averaged the remaining channels to produce a continuum dataset. We imaged the resulting visibilities with the tclean task using the multi-scale CLEAN algorithm with scales of 0, 1, 3 and 6 times the beam FWHM, and an elliptic CLEAN mask encompassing the disk emission.

To reduce the size of the data, we time averaged it to 6.06 seconds, i.e., 3 times the original integration time.

After imaging, one EB image appeared of much lower SNR and therefore the corresponding visibilities were rejected. The individual images of the three remaining execution blocks (EBs 0,1,3) did not appear astrometrically offset with respect to each other, which is as expected because they were taken very close in time. As the fluxes of all EBs match within 2%, we concatenated the three EBs and self-calibrated them all together. To determine a good initial model for the self-calibration, we used multi-scale cleaning with the tclean task using a threshold of ∼7 times the rms noise level of the image. Using the tasks gaincal and applycal, we corrected for phase offsets between spectral windows, and between polarizations considering a solution interval of the scan length (solint=inf). Another iteration of phase selfcalibration was done with a solution interval of 30s. We reached an overall improvement in peak SNR of 34% after self-calibrating the LB19 data. The LB19 data were combined with archival observations previously published in Isella et al. (2019) Isella et al. (2019) where the procedure for the self-calibration of SB16 and IB17 data is described in detail. For all datasets, we used the statwt task to weight the visibilities according to their scatter. Before combining the LB19 data with the previously published data, we fitted an elliptical ring to the maximum of the outer ring in the image plane, for all datasets separately, to derive the center of the image and then used the fixvis task to shift the image to the phase center, and assign it to a common phase center using the fixplanets task on the center coordinate derived by Isella et al. (2019) . The fluxes of the executions in LB19 differed by ∼3% from the archival datasets (IB17+SB16; Isella et al. 2019 ) and were rescaled using the DSHARP rescale flux function 1 . After concatenation of the data, we followed the same procedure as explained above, with three rounds of phase selfcalibration.

We proceeded with imaging of the final data using CLEAN. In a normal CLEAN workflow, after the CLEAN iterations terminate when the peak value of the residual image drops below a threshold value (4× rms noise level in the observations considered here), a restored CLEAN model is combined with the residual image to form the CLEANed image. As discussed in Czekala et al. in press, however, the units of these two quantities differ: the units of the restored CLEAN model are Jy {CLEAN beam} −1 while the units of the residual image are Jy {dirty beam} −1 , since it originated as the dirty image. When the CLEAN beam (typically chosen to be an elliptical Gaussian) poorly approximates the dirty beam (as is common with multi-configuration ALMA datasets), the normal CLEAN workflow produces a CLEANed image with an incorrect flux scale and compromised image fidelity, especially for faint emission. This phenomenon was first described in Jorsater & van Moorsel (1995) , and so we term it the "JvM effect". To correct for the unit mismatch, before combining the residual image with the restored CLEAN model, we first rescaled the residual image by the ratio of the CLEAN beam / dirty beam "volumes" (see "JvM correction", Czekala et al. in press).

To test the effect of the angular resolution on the image features and assess their robustness, we performed a grid of CLEANed and JvM-corrected images, using Briggs weighting (Briggs & Cornwell 1992) with different robust parameters. A gallery of continuum images (and corresponding fluxes), synthesized from the new dataset alone (LB19) and from dataset combinations includ- (Table 4 ). We note that while the uv coverage and sensitivity are maximized when all datasets are combined (LB19+IB17+SB16), such a combination does not take into account the intrinsic changes of the emission that are due to the rotation of the system, and the change in the location of the dust surrounding the planets. Based on the orbital solutions of Wang et al. (2021) , we expect a motion of ∼14 mas for both planets between December 2017 and July 2019.

3. RESULTS Figure 1 presents a selection of images of the continuum emission of PDS 70 at 855 µm, synthesized from the new ALMA observations combined with short baseline data (LB19+SB16). The disk is well detected with a spatially integrated flux density of ∼176±18 mJy (all images give similar values). After deprojecting the image with an inclination of ∼51.7 • and a position angle of ∼160.4 • , we computed an azimuthally averaged radial profile and found that the outer disk resolves in a ring extending radially from ∼0.4 (45 au) and ∼0.9 (100 au). The outer disk is not radially symmetric and shows a clear azimuthal asymmetric feature in the North West (∼27% brighter at peak compared to the mean ring value), as already discussed by Long et al. (2018) and Keppler et al. (2019) . When imaged at high resolution, the outer disk resolves into a narrow and bright ring with a faint inner shoulder detected in the image at the 3-4σ level (Appendix A.2). To better assess the presence of such substructures, we model the azimuthally averaged radial visibility profile using the frank package (Jennings et al. 2020 ). Our analysis, presented in Appendix B recovers a double peaked profile for the outer disk. Such a substructure was already hinted in the data presented in Keppler et al. (2019) . Inward of the outer disk, the dust-depleted cavity includes an inner disk that radially extends up to 0.16 (18 au) and presents faint additional emission in the West and in the South of the inner disk that will be discussed in the next subsection.

Within the cavity, the inner disk appears well resolved with an integrated flux ranging between 727±27 µJy and 888±59 µJy depending on the dataset (Table 5) . When imaged at high angular resolution (e.g., Figure 1 , left), it appears irregular and the emission is discontinuous in the North.

Continuum emission is also detected near the locations of the planets, confirming the findings of Isella et al. (2019) . We use the same nomenclature as Isella et al. (2019) and label the continuum emission located close to planet b and c, b smm and c smm , respectively. The continuum emission around PDS 70 c, c smm , is recovered in all images, and in particular in the standalone new, high resolution, dataset (LB19), where it appears as a 5.4 to 16σ feature depending on the robust parameter. c smm clearly separates from the outer disk when imaged at resolutions finer than ∼40 mas. It appears unresolved even at our best angular resolution (∼20 mas;∼2.3 au). We find that its peak intensity is similar in all the images that spatially resolve it from the outer disk (see Appendix A.3), confirming its point-source nature. Depending on the dataset (IB17+SB16 or LB19+SB16) and the robust parameter, its peak intensity ranges between 80±6 and 107±15 µJybeam −1 . In the following, we will consider 86±16 µJybeam −1 as a reference for further discussion. The emission located near PDS 70 b, b smm , is on the other hand, only recovered when the new high resolution data is combined with short baselines, and when the beam is larger than ∼50 mas. This indicates that it is low surface brightness, extended emission. Its peak intensity and morphology vary greatly between images of different datasets (Table 5) , which makes its morphology and properties difficult to recover accurately.

In order to assess whether the signal within the cavity could result from imaging artifacts, following Andrews et al. (2018) , we subtracted the Fourier transform of the CLEAN model of the outer disk, after blanking out the pixels within the cavity (using an elliptical mask of 0.25 × 0.4 ), and image and model the visibilities carrying the residual signal from within the cavity. Figure 2 show two residual images, hereafter called 'cavity images', for LB19+SB16 and IB17+SB16, that clearly show that the inner disk emission and c smm are recovered in both epochs, the latter with a significance up to 18σ. On the other hand, b smm is detected at a 3σ level only in some cavity images obtained from combined datasets. A gallery of cavity images are given in the Appendix A.2.

As an additional test, we perform a model fit of the cavity visibilities using the dataset LB19+SB16, obtained after subtracting the Fourier transform of the CLEAN model of the outer disk using a robust parameter of 1. We consider a simple model for all three sources of emission within the cavity, namely the inner disk, b smm and c smm , compute the Fourier transform using galario (Tazzari et al. 2018 ) and explore the parameter space using the Monte Carlo Markov chains implementation in emcee (Foreman-Mackey et al. 2013). Our model consists in a Gaussian ring for the inner disk, that enables to model an additional structure within the inner disk, a point source for c smm (between PA=250 • and 280 • ), and a circular Gaussian for b smm , located in the South (between PA=70 • and 250 • ). A uniform prior was used over the allowed range for each parameter. Our best-fit model and residual maps are shown in Figure 3 , and cor- Table 2 . Best-fit parameters for the model to the cavity data for the datasets LB19+SB16 and IB17+SB16, with the 1σ error. The flux, radial peak position, and width of the Gaussian for the inner disk are finn, rinn, σinn, respectively. The total flux and polar coordinates in the disk plane of bsmm and csmm are f b , r b , θ b and fc, σc, rc, θc, respectively. The relative apparent astrometry ∆(RA, Dec) is also provided. From the orbital fits of Wang et al. (2021) , the expected motions of the planets between the epoch of the long baselines observations (December 2017 and July 2019) is similar for both, ∼14 mas, smaller than the angular resolution of our observations. To search for possible motion of c smm between the two epochs, we performed the same modeling as above on the IB17+SB16 dataset. b smm was not recovered in this fit, but the inner disk and c smm were. Using the best-fit positions for c smm at the two epochs, and considering a 2 mas error in the centering of the two datasets, we find marginal evidence for a movement of the peak position of 10.9±6.9 mas. We note that the nominal positional accuracy is defined as beam FWHM /SNR/0.9 (Thompson et al. 2017 , and ALMA Cycle 8 2021 Technical Handbook), with 0.9 a factor to account for a nominal 10% signal decorrelation. We consider two images in which c smm is imaged at a decent SNR and separated from the outer disk, LB19+SB16 (r=0.5) and IB17+SB16 (r=-0.3). With corresponding SNR of 8.9σ and 7.1σ on the peak intensity of c smm respectively, and a beam FWHM of 36 and 60 mas, respectively, the positional accuracies are ∼4.5 mas and 9.4 mas, respectively, comparable to the uncertainty that we derived for the apparent displacement of c smm . Additional observations with ALMA in the coming years, providing a longer time baseline, are needed to confirm such a movement.

A circumplanetary disk around PDS 70 c -Isella et al. (2019) reported the detection of c smm using ∼67 mas resolution observations. We confirm this detection with higher angular resolution observations that enable us to separate the emission from the outer disk. Given that the location of c smm is very close to the existing Hα and NIR measurements of PDS 70 c (Isella et al. 2019) , and to the expected positions of PDS 70 c at the time of our observations (Figure 3 ), we interpret it as tracing the millimeter emission of dust grains located in a CPD. Assuming that c smm is optically thin, its flux density can be converted into a dust mass estimate, for a given dust opacity and temperature. We note that if the emission is optically thick, such an assumption would provide a lower limit in the dust mass. The CPD temperature is also uncertain. It is determined by the sum of various sources of heating, namely viscous heating due to accretion of material through the CPD, accretion shocks, and external irradiation from both the planet and the star (Isella et al. 2014 (Isella et al. , 2019 Andrews et al. 2021 ). Using 2 M Jup , 2 R Jup , and 1055 K as the mass, radius, and temperature of PDS 70 c ), a mass accretion rate of 10 −8 M Jup /year ), we find that at a radial distance of 1 au from the planet, T vis = 3 K, and T p,irr = 18 K. Considering a stellar-irradiation temperature of T s,irr = 24 K at the location of PDS 70 c (obtained from the radiative transfer model of Keppler et al. (2019) ), the CPD temperature at 1 au is T 4 CP D = T 4 vis + T 4 p,irr + T 4 s,irr , that is T CP D ∼ 26 K. Considering a typical dust opacity for 1 mm sized grains of 3.63 cm 2 g −1 (Birnstiel et al. 2018 ) and a temperature of 26 K, we estimate a CPD dust mass of ∼0.007 M ⊕ . A lower dust mass would be inferred if the dust temperature is higher than considered here (Schulik et al. 2020) .

However, PDS 70 c is massive enough to carve a gap, and, as a consequence, large grains are trapped in a pressure maximum in the outer disk while small grains, well coupled to the gas, can flow inward. This is confirmed by the different cavity outer radii measured in scattered light compared to mm wavelengths (probing small and large grains, respectively; Keppler et al. 2019 ). The CPD is therefore only replenished with small dust particles that leak into the cavity ) through meridional flows from the upper protoplanetary disk layers (e.g., Kley et al. 2001; Ayliffe & Bate 2009 ). If the CPD contains only small 1 µm sized grains (with an opacity of 0.79 cm 2 g −1 ; Birnstiel et al. 2018 ) the CPD dust mass increases to ∼0.031 M ⊕ . It is of course possible that the CPD hosts a range of particle sizes if the grains can grow. Bae et al. (2019) find that, if a steady state is achieved between the mass inflow to the CPD and the mass accretion rate onto the planet, the amount of sub-micron grains in the CPD would largely underestimate the observed mm flux and that accumulation of grains beyond the steady-state amount and/or in-situ grain growth is needed to account for it. In Appendix D, we show the range of dust masses that the CPD would have for various dust grain size distributions, as a function of the maximum grain size. With these mass estimates, the ratio between the CPD dust mass and the planet mass, considering 2 M Jup , ranges between 1 and 5×10 −5 .

If small grains can grow to mm sizes within the CPD, they could rapidly be lost as they efficiently drift toward the planet and it only takes 100-1000 years for an accreting CPD to lose all its mm dust ). However, as in protoplanetary disks, local gas pressure maxima can act as particle traps, and prevent these grains from drifting. Interestingly, this can naturally occur in CPDs. Most of the gas that is feeding the CPD through meridional flows is then radially flowing outward in a decretion disk. The balance between the sub-Keplerian headwind and viscous outflow associated with a decretion flow leads to a global dust trap (Batygin & Morbidelli 2020) . As a consequence, dust grains with sizes 0.1-10 mm may be trapped in the CPD and as the dust-to-gas ratio increases, streaming instabilities might be triggered (Drażkowska & Szulágyi 2018) , or gravitational fragmentation in the outer regions of the CPD (Batygin & Morbidelli 2020) that will eventually lead to the formation of satellitesimals. At the same time, dust particles can accrete via pebble-accretion onto the satellitesimals formed in situ or captured from the disk edge (e.g., Ronnet & Johansen 2020) .

Our observations also put a strong constraint on the spatial extent of the CPD as seen in the dust emission at mm wavelengths. The emission c smm is unresolved even at our highest angular resolution, and its peak intensity is similar over a range of beam sizes, until ∼40 mas, beyond which the CPD does not separate from the outer disk anymore (Appendix A.3). This indicates that it is more compact than 1.2 au in radius. On the other hand, there is a lower limit to the CPD extent needed to account for the observed flux. Assuming that it is a uniform disk with an optical depth of 1, and considering a temperature of 26 K, we find that it has a radius of 0.58 au. These two values (0.58 and 1.2 au) are therefore the lower and upper limits on the CPD radial extent constrained from our observations. The CPD is expected to be truncated (in gas) at a third of the Hill radius, which for PDS 70 c, assuming a planet mass of 2 M Jup at 34 au, is 1/3 ×3.1 ∼1 au. 3D simulations show that isothermal CPD are bound within 10% of the Bondi radius (Fung et al. 2019) , that is 1/10 ×11 ∼1.1 au for PDS 70 c assuming a local temperature of 26 K. Both estimates are therefore consistent with our constraints. However, we cannot rule out the possibility that the gas component of the CPD extends beyond the dust component, in particular if some dust grains in the CPD drift inward. (Hashimoto et al. 2020) . It is interesting to understand why PDS 70 b, at the sensitivity of our observations, does not seem to host a compact, dusty, circumplanetary disk as PDS 70 c does. A possibility would be that PDS 70 b has a much smaller Hill radius than PDS 70 c, as it orbits at smaller separation. Another natural explanation could be that PDS 70 b is starved of dust grains, as only the small grains that leak through the orbit of PDS 70 c and are transported through a streamer from the outer to the inner planet would enter the region of influence of PDS 70 b. Finally, it could be that the nature of the CPD is different around the two planets, with a decretion disk around PDS 70 c, and an accretion disk around PDS 70 b that is fed through a streamer coming from PDS 70 c rather than through meridional flows. More theoretical work looking at formation of CPDs in systems hosting two giant planets is needed to assess the potential differences between CPD formation in the inner and outer planet.

Inner disk. -An inner dusty disk, evidenced in the IR spectral energy distribution and scattered light images is also clearly detected in our images up to ∼0.16 (∼18 au) (see also, Long et al. 2018; Keppler et al. 2019) . Considering that the planets are filtering material from the outer disk such that only small dust particles can flow in the cavity, as for the CPD, it is unclear whether the inner disk mm emission is due to a population of small or large dust grains. To address this question, we computed the dust surface density and optical depth radial profiles of the continuum emission, using the combined dataset (SB16+IB17+LB19) imaged with robust=1. We consider 4 models for the dust grain population, that follow a size distribution n(a)da ∝ a −3.5 da with a maximum grain size a max of 10 µm, 100 µm, 1 mm and 1 cm, and a minimum size of 0.05 µm. We use the DSHARP opacities (Birnstiel et al. 2018 ) and the temperature profile output of the radiative transfer model of Keppler et al. (2018) . The dust surface density as well as the total optical depth τ ν is numerically computed, considering scattering and absorption opacities (Sierra & Lizano 2020; Sierra & MAPS team 2021) . Figure 4 , left, shows the total optical depth τ ν for all 4 models. The right panels show the dust surface density profiles (top) and corresponding cumulative masses (bottom). The dust surface density is maximum at the outer disk that is obviously the disk region that contributes to most of the dust mass (∼0.24×10 −3 M for a max = 1 mm). We note that without the inclusion of scattering, the optical depth would follow the curve of the dust population with a max =10 µm, as the albedo at mm wavelengths is negligible for these small grains. In all these models, the inner disk is optically thin, with a total dust mass of ∼ 2 × 10 −7 − 10 −6 M (i.e., 0.08-0.36M ⊕ ).

It therefore appears that the emission at 855 µm from the inner disk regions located within the orbit of PDS70 b could be accounted for by a population of small grains. Interestingly, we note that the near infrared excess apparent in the spectral energy distribution of PDS70 is very low ). This emission is mostly due to the thermal emission of small grains located within the innermost au and such a low excess could indicate a low small-dust mass content in the inner disk, and therefore suggest the additional presence of larger dust grains in order to account for the measured flux at 855 µm. However, the inner disk emission in the infrared could still be optically thick ), making it difficult to directly relate to our sub-millimeter observations and multiple wavelength observations in the millimeter regime are needed to constrain the grain size population in the inner disk. We note that the brightness temperature might be under- estimated near the star because of our limited angular resolution and that it is possible that the innermost disk regions are optically thick also at sub-millimeter wavelengths. The longevity of the inner disk remains unclear; the replenishment flow is controlled by the planets, if it is so strongly depleted (in gas) it may not allow grains to grow efficiently. It is possible that some of the dust in the inner disk is of second generation produced by collision of larger bodies, perhaps stirred up by PDS 70 b. The star exhibits a small, but non negligible, mass accretion rate, for which an additional mass reservoir in the inner disk, such as a dead zone, was recently suggested (Thanathibodee et al. 2020) . Determining the physical conditions there-in, in particular the dust to gas ratio, would be crucial for understanding whether such an inner disk can still grow terrestrial planets within a system hosting two outer giant planets. The current dust mass estimates are so low that it is unlikely that planets could form through pebble accretion (Lambrechts et al. 2019) .

Outer disk structure. -Our observations at high angular resolution indicate that the outer disk hosts substructures. In addition to an 'arc' in the North-West, already seen at lower resolution images (Long et al. 2018; Keppler et al. 2019) , it resolves into two components, that can be either a double-ring structure with a dip at ∼0.55 or a bright ring with an inner shoulder. Interestingly, Huang et al. (2020) also find with high resolution observations, a two-component structure in GM Aur, with a bright ring and an outer shoulder. It is unclear if such two-component structure in PDS 70 could be due to a secondary gap induced by PDS 70 c as an outer sec-ondary gap opens only when the disk is sufficiently cold (Bae & Zhu 2018b) , with (h/r) p 0.06 where (h/r) p is the disk aspect ratio at the location of the planet ((h/r) p 0.08 at PDS 70 c's location; Bae et al. 2019 ). On the other hand, recent three-dimensional planet-disk interaction simulations including both gas and dust components showed that dust grains at the gap edge can have radial structures (Bi et al. 2021) , potentially induced by corrugated vertical flows driven by the spiral wave instability (Bae et al. 2016a,b) or meridional flows (Fung & Chiang 2016) . Alternatively, such a substructure could be due to the presence of an additional, yetundetected low-mass planet embedded within the outer disk. Similar multiple-ring substructures were also observed in other transition disks, such as HD 169142 in which three narrow rings were found and interpreted as tracing a migrating 10 M ⊕ in a low viscosity disk . However, hydrodynamical simulations show that thermodynamics can dramatically affect the structure of gas and dust, with different disk cooling timescales leading to different planet-induced substructures (Facchini et al. 2020) . Further chemical surveys will help to constrain the density and temperature structures (Facchini et al. 2021) , enabling to test the possibility that an additional, low-mass planet is responsible for the structured outer ring and constrain the mass and radial location of that planet. We note that it is unlikely that an additional planet within the outer continuum ring disrupt the planetary system. In a two-planet system neglecting the eccentricity damping from the protoplanetary disk gas, the plan-ets can avoid close encounters and are Hill-stable when their orbital separation is greater than 3.46 R H , where R H = a 1 [(M 1 + M 2 )/3M * ] 1/3 is the mutual Hill radius (Gladman 1993; Barnes & Greenberg 2006) . The addition of a third planet generally makes the stability criteria more stringent because the conservation of the total angular momentum and energy can no longer guarantee the avoidance of close encounters even for initially large separations beyond the Hill-stability criteria (Tamayo et al. 2015) . However, provided that the protoplanetary disk gas provides sufficient eccentricity damping, Tamayo et al. (2015) argued that the two-planet criteria can still be used in three-planet cases. Assuming a range of 1 − 10 M Jup for PDS 70 c and a Saturn mass for the hypothesized additional planet, this criteria is met when the latter is located beyond 44 − 53 au. Therefore, the system would be dynamically stable if the additional planet is located within the dip in the outer continuum ring at ∼ 60 au. Future numerical simulations will allow to further test our conclusions.

In this Letter, we report new ALMA observations obtained at high angular resolution (∼20 mas) at 855 µm of the PDS 70 system. We confirm the tentative detection by Isella et al. (2019) of a compact source colocated with the position of PDS 70 c with an independent dataset at higher angular resolution. These new observations provide the most compelling evidence of the presence of a CPD around an accreting planet to date. Future molecular line infrared observations at very high angular resolution may be able to detect rotating gas around PDS 70 c, providing conclusive results on the nature of the continuum mm emission. The detec-tion of unresolved (r < 1.2 au) emission around planet c confirms that circumplanetary material is able to retain dust for long timescales, as required in satellite formation models.

These ALMA observations shed new light on the origin of the mm emission close to planet b. The emission is diffuse with a low surface brightness and is suggestive of a streamer of material connecting the planets to the inner disk, providing insights into the transport of material through a cavity generated by two massive planets. The non-detection of a point source around PDS 70 b indicates a smaller and/or less massive CPD around planet b as compared to planet c, due to the filtering of dust grains by planet c preventing large amounts of dust to leak through the cavity, or that the nature of the two CPDs differ. We also detect a faint inner disk emission that could be reproduced with small 1 µm dust grains, and resolve the outer disk into two substructures (a bright ring and an inner shoulder).

PDS 70 is the best system to date to study and characterize circumplanetary disks, but also planet-disk interactions and disk cavity clearing by massive planets. The two massive planets, likely migrating outward in a grand tack-like scenario , are reminiscent of the Jupiter-Saturn pair, at larger distances from the star. Detailed studies of the circumplanetary disks, and of the leakage of material through the cavity, will provide strong constraints on the formation of satellites around gas giants, and on the ability to provide the mass reservoir needed to form terrestrial planets in the inner regions of the disk. Upcoming studies of the gas kinematics and chemistry of PDS 70 will complement the view provided by this work, serving as a benchmark for models of satellite formation, planet-disk interactions and delivery of chemically enriched material to planetary atmospheres. Table 3 provides the observing log of the new ALMA observations presented in this Letter. EB2, indicated in italic, was not included in the images.

To test the effect of the angular resolution on the image features and assess whether they are recovered in various images, we performed a grid of CLEANed and JvM-corrected images using different datasets, and different Briggs robust parameters. Figure 5 presents the resulting images. Corresponding image properties and fluxes are reported in Tables 4 and 5. Figure 6 shows the residual images (called 'cavity images') obtained after subtracting the Fourier transform of the CLEAN model of the outer disk, for robust values of 0.5, 1 and 2.0. A.3. Peak intensity of the continuum emission associated with the planets Figure 7 shows the peak intensity of c smm as a function of angular resolution. Depending on the dataset, and the robust parameter, its peak intensity ranges between 80±6 and 107±15 µJybeam −1 when it is well separated from the outer ring. At larger resolution than ∼60 mas, the peak intensity increases because the beam contains contribution from the outer disk. The grey area reports the estimate of Isella et al. (2019) . In contrast the peak intensity of b smm varies between 46±5, 56±6 and 27±4 µJy beam −1 for three different datasets (LB19+SB16; IB17+SB16; LB19+IB17+SB16, respectively) imaged at resolutions of 51 mas×44 mas, 87 mas×69 mas, 63 mas×54 mas, respectively (see Table 5 ). Table 4 . Summary of disk and CPD properties for various datasets. The csmm peak intensities were computed with the CASA task imstat in an aperture centered on the CPD, with major/minor axis twice the beam major/minor axis (col. 7) and with a Gaussian fit when possible (col. 8). The rms is computed considering an annulus between 2.4 and 6 . We considered 10% calibration uncertainty as the flux uncertainty. Figure 7 . Peak intensity of csmm as a function of angular resolution. The peak intensity of the CPD is constant around ∼80-105 µJy beam −1 when it is well separated from the outer ring. At larger resolution than ∼60 mas, the peak intensity increases because the beam contains contribution from the outer disk. The grey area reports the estimate of Isella et al. (2019) .

To better assess the presence of substructures in the outer disk, we fit azimuthally averaged deprojected visibilities using the frank package that models an axisymmetric surface density profile using a flexible Gaussian process (Jennings et al. 2020) . To do so, we considered the combined dataset LB19+IB17+SB16, which has the best uv coverage, assuming that the outer disk brightness distribution has not changed between these observations. The data was averaged into intervals of 30 s and 1 channel per spectral window to reduce data volume.

As the disk presents an asymmetric arc-like feature in the North-West that can lead to overestimate the disk radial intensity profile when fitting with an axisymmetric model, we followed the procedure developed in Andrews et al. (2021) that modifies the visibilities by removing a model for the asymmetry before fitting with frank. We mask the emission between position angles of -85 • and 40 • . The asymmetry is selected in the CLEAN model image, and the mean radial profile of the CLEAN model from the disk outside the asymmetric region is subtracted from the model image, allowing us to obtain a model for the asymmetry only. The Fourier transform of the asymmetry model was then subtracted from the original visibilities and the resulting set of visibilities are further modelled. frank fits deprojected visibilities, and inaccurate estimates of the geometric parameters for the deprojection, (∆RA, ∆Dec, inc, PA), can lead to significant residuals. Automatic procedures performed poorly to find the best parameter set and we therefore optimized those parameters by hand, exploring different values of spatial offset and geometry in steps of 1 mas and 0.5 deg. The final values adopted for the frank fit are (12 mas, 15 mas) for (dRA, dDec), and (49.5, 161.0) for (inc, pa). We tested the sensitivity of the fit to the hyperparameters α (varied between 1.05 and 2.00) and w smooth (varied between 10 −4 to 10 −1 ) and found no significant difference on the residuals. We therefore fixed for standard values w smooth = 10 −4 and α=1.30.

The fit to the data and the corresponding profile are shown in Figure 8 (top panels). The best fit model indicates a radial profile with two local maxima for the emission of the outer disk, confirming the findings of Keppler et al. (2019) with lower resolution observations. At the angular resolution of our observations, the two peaks are separated by ∼7 beams. It is however unclear whether the outer disk is composed of two separated, broad, rings, or of a ring with an inner shoulder. No clear gap or ring is evidenced in the inner disk. We note the presence of a possible additional shoulder at 0.85 . The model and residuals are imaged exactly as the observations, with a robust parameter of 0.5, and are presented in Figure 8 , bottom panels. The residuals indicate that the axisymmetric model does not account for the full complexity of the emission, in particular the disk can be vertically thick and flared. A dedicated 2D modeling of the visibilities, that is beyond the scope of this Letter, is needed to properly assess the morphology of the disk and will be presented in a forthcoming study. Table 6 provides the published astrometry of PDS 70 b and PDS 70 c, as well as the predicted locations at the time of our observations based on the best orbital fits by Wang et al. (2021) .

We consider 3 models for the dust grain population in the CPD around PDS 70 c, that follow different size distribution n(a)da ∝ a −3.5 da,∝ a −3.0 da, and ∝ a −2.5 da and show the predicted dust mass as a function of the maximum grain size a max in Figure 9 . We consider a minimum grain size of 0.05 µm, a temperature of 26 K and use the DSHARP opacities (Birnstiel et al. 2018 ). 

Astronomical Society of the Pacific Conference Series

Astronomical Society of the Pacific Conference Series

& MAPS team

Interferometry and Synthesis in Radio Astronomy

We are very grateful to Marco Tazzari and Antonella Natta for their support during the Covid-19 pandemic. We thank the North American ALMA ARC for their help. We acknowledge Ryan Loomis, Jason Wang and Rens Waters for insightful discussions and the reviewer for helpful suggestions that improved the manuscript. This Letter makes use of the following ALMA data: ADS/JAO. (Jennings et al. 2020)