Highly structured disk around the planet host PDS 70 revealed by high-angular resolution observations with ALMA

Context. Imaged in the gap of a transition disk and found at a separation of about 195 mas ( ∼ 22 au) from its host star at a position angle of about 155 ◦ , PDS70b is the most robustly detected young planet to date. This system is therefore a unique laboratory for characterizing the properties of young planetary systems at the stage of their formation. Aims. We aim to trace direct and indirect imprints of PDS70b on the gas and dust emission of the circumstellar disk in order to study the properties of this ∼ 5 Myr young planetary system. Methods. We obtained ALMA band 7 observations of PDS70 in dust continuum and 12 CO(3–2) and combined them with archival data. This resulted in an unprecedented angular resolution of about 70mas ( ∼ 8au). Results. We derive an upper limit on circumplanetary material at the location of PDS70b of ∼ 0.01 M ⊕ and ﬁnd a highly structured circumstellar disk in both dust and gas. The outer dust ring peaks at 0.65 (cid:48)(cid:48) (74 au) and reveals a possible second unresolved peak at about 0.53 (cid:48)(cid:48) (60 au). The integrated intensity of CO also shows evidence of a depletion of emission at ∼ 0.2 (cid:48)(cid:48) (23au) with a width of ∼ 0.1 (cid:48)(cid:48) (11au). The gas kinematics show evidence of a deviation from Keplerian rotation inside (cid:46) 0.8 (cid:48)(cid:48) (91 au). This implies a pressure gradient that can account for the location of the dust ring well beyond the location of PDS70b. Farther in, we detect an inner disk that appears to be connected to the outer disk by a possible bridge feature in the northwest region in both gas and dust. We compare the observations to hydrodynamical simulations that include a planet with different masses that cover the estimated mass range that was previously derived from near-infrared photometry ( ∼ 5–9 M Jup ). We ﬁnd that even a planet with a mass of 10 M Jup may not be sufﬁcient to explain the extent of the wide gap, and an additional low-mass companion may be needed to account for the observed disk morphology.


Introduction
In recent years, high angular resolution observations of protoplanetary disks have revolutionized our view of disk evolution and showed that small scale structures, such as concentric rings and spiral arms are ubiquitous (e.g., Andrews et al. 2018;Long et al. 2018a), suggesting that planet formation might occur very early in the young stellar system history.Although these substructures are often interpreted as direct imprints of planet-disk interactions, it is still challenging to understand and constrain the architectures of planetary systems needed to account for them (e.g., Bae et al. 2018), or to rule out alternative scenarios (e.g., magneto-hydrodynamical instabilities, Ruge et al. 2016;Flock et al. 2017).In addition, an accurate determination of the properties of young planets (e.g., luminosity and mass at a given age) is needed to constrain the formation mechanisms at work (e.g., Mordasini et al. 2017).
The theory of interactions of embedded planets with their natal environment, the protoplanetary disk, and their relation to the observational signatures have been studied by many authors (e.g., Paardekooper & Mellema 2004;Jin et al. 2016;Dipierro et al. 2018;Liu et al. 2018;Zhang et al. 2018).Currently the most promising methods to understand the interaction of a young planet with its environment and its further evolution are to detect it through direct imaging (e.g., Keppler et al. 2018) or through the perturbations that it induces in the disk velocity field (e.g., Pérez et al. 2015).
Current direct imaging infrared surveys reach detection limits of a few Jupiter masses (e.g., Maire et al. 2017;Uyama et al. 2017) but are often limited by the presence of bright and complex disk features.Indeed, numerous claims of companion candidates in disks that show asymmetric features are still debated (e.g., HD100546, HD169142, MWC758, LkCa15; see Quanz et al. 2015;Follette et al. 2017;Rameau et al. 2017;Biller et al. 2014;Reggiani et al. 2014;Ligi et al. 2018;Reggiani et al. 2018;Kraus & Ireland 2012;Sallum et al. 2015;Mendigutía et al. 2018) and require confirmation through additional observations at e.g.different filter bands.
The presence of three different planets in the disk around HD 163296 was claimed by two teams with a complementary method based on perturbations in the Keplerian velocity field of the disk.Pinte et al. (2018) detected a localized (both in space and velocity) deformation of the isovelocity curves in 12 CO transitions, consistent with the spiral wake induced by the presence of a 2 M Jup planet at 260 au.Teague et al. (2018a) measured the rotation velocity curves of CO isotopologues as a function of distance to the star, and found local pressure gradients consistent with gaps carved by two ∼1 M Jup planets at 83 au and 137 au.
Using the VLT/SPHERE instrument and complementary datasets covering multiple epochs and various near-infrared (NIR) wavelengths, we recently discovered a companion to the 5.4±1.0Myr old (Müller et al. 2018) and 113.4±0.5 pc distant (Gaia Collaboration et al. 2016, 2018) T Tauri star PDS 70 (Keppler et al. 2018;Müller et al. 2018).Comparison of the NIR photometry to evolutionary models implies that the companion is in the planetary mass regime (∼5-9 M Jup , Keppler et al. 2018) consistent with the mass range inferred from atmospheric modeling (∼2-17 M Jup , Müller et al. 2018).PDS 70 b is located at a projected separation of about 22 au from the central star, within the large gap of its host transition disk, between an inner disk and a well-resolved outer disk (Hashimoto et al. 2012(Hashimoto et al. , 2015;;Long et al. 2018b;Keppler et al. 2018).Follow-up direct imaging observations with MagAO in the Hα line enabled a 2-3σ detection of the companion at two different epochs, and imply that it is likely still accreting gas from the disk (Wagner et al. 2018).This object is therefore a unique case of a directly imaged planet still shaping its natal environment.
In this paper, we present new ALMA band 7 observations of PDS 70 obtained in Cycle 5. We combined the data with archival observations presented by Long et al. (2018b) obtaining an unprecedented angular resolution of ∼0.07 .In Sect. 2 we describe the observing setup and data reduction, Sect. 3 presents our results which are discussed and compared to hydrodynamical simulations in Sect. 4.

Observations and data reduction
We obtained ALMA Cycle 5 DDT observations (Project ID: 2017.A.00006.S, PI: M. Keppler) of PDS 70 in band 7 on December 2nd, 3rd and 6th, 2017 under very good weather conditions (mean pwv ≤ 0.9 mm).For three of the four spectral win-dows, the correlator was tuned to a center frequency of 357.2, 355.3 and 344.3 GHz for continuum observations in dual polarization mode with a bandwidth of 2.0 GHz.The fourth spectral window was centered around the 12 CO(3-2) transition at 345.3 GHz with a bandwidth of 0.938 GHz.The quasars J1427-4206, J1337-1257 and J1517-2422 were used as bandpass, phase and flux calibrators.The calibration was performed using the Common Astronomy Software Package (CASA), version 5.1.1.The total on-source integration time was 1.9 hours.
Since the extended antenna configuration is filtering out the largest spatial scales in the disk, we made use of the archival Cycle 3 data taken in a similar spectral setup and presented by Long et al. (2018b) to recover the short baselines.Details regarding the observing strategy and setup are described in Long et al. (2018b).We transferred both Cycle 3 and Cycle 5 data to CASA v.5.3.0 and subtracted the continuum emission from the line data using the task UVCONTSUB.We corrected the phase center of the Cycle 3 data for the shift due to the proper motion of the star ((-29.7, -23.8) mas/yr, Gaia Collaboration et al. 2016, 2018) with respect to to the Cycle 5 data set.We then combined the two data sets and shifted the phase center by an amount of (0.509 , 0.490 ), which was found to be the center of the disk by fitting a two-dimensional Gaussian to the continuum Cycle 5 emission using the UVMODELFIT tool.We finally used the task TCLEAN for imaging, applying Briggs weighting with a robust parameter of 0.5.Since self-calibration of both continuum and CO data did not significantly improve the images, we will base our analysis on the non-self-calibrated data.The resulting beam size for the dust continuum at a mean frequency of 350.6 GHz (855 µm) is 74 × 57 mas (8.4×6.5 au) with a position angle (PA) of 63 o .We measure an rms noise level of 0.026 mJy beam −1 from emission free regions.For the CO, we obtained a beam size of 76 × 61 mas (8.6×6.9 au) with a PA of 60 o and a channel width of 425 m/s.The noise level per channel is determined to be 1.26 mJy beam −1 .

855 µm dust continuum
Figure 1 (right column) shows the continuum image at 350.6 GHz (855 µm).The disk is detected at a high signal-tonoise ratio (∼65 at the peak).The integrated flux density of the disk inside 1.3 after applying 2σ clipping is 230±23 mJy, where the error bar is dominated by the ∼10% uncertainty of the absolute amplitude calibration of ALMA in band 7 1 .This is consistent with the value found by Long et al. (2018b).The dust continuum shows evidence of a large cavity, a dust ring with a brightness distribution that is slightly asymmetric, both in radial and azimuthal direction, an inner disk, as well as a possible bridge feature, all of which we will describe in the following paragraphs.
By fitting a two-dimensional Gaussian to both datasets using the task UVMODELFIT we find a disk inclination of 51.7±0.1 o and 52.1±0.1 o , and a PA of 156.7±0.1 o and 159.7±0.1 o , for the Cycle 5 and Cycle 3 datasets, respectively.We verified the inclination by using only short baselines (<150 kλ, which correspond The bottom row provides a closer view of the observations including annotations where the color scaling has been stretched to bring out detail.The contours for the 12 CO are starting at 20% of the peak value to the peak in steps of 10%.For the continuum, the gray dashed contour is 5σ with black contours starting at 10σ and increasing in steps of 10σ where σ = 26 µJy beam −1 .The synthesized beams are shown in the bottom left of each panel. to the location of the null in the real part of the visibilities, see Fig . A.4) for the Gaussian fit, which ensures that the cavity is not resolved, as well as by using a disk model.These efforts yielded in all cases similarly good fits with consistent values for the inclination within 3 o .We note however that all these models assume axial symmetry and therefore none of them reproduces the real morphology of the disk.Considering the complexity of the continuum emission that appears to be highly structured, such simple modeling appears limited and we adopt a final value of 51.7 o , as it corresponds to the model with the least assumptions.

Disk radial and azimuthal morphology
Figure 2 (uppermost, gray line) shows the azimuthally averaged and deprojected radial profile of the dust continuum which clearly reveals a large gap and a ring component.The emission strongly decreases inside the ring where the flux is reduced by more than 90%.
The radial profile of the ring is asymmetric, best seen in the cuts along the major and minor axes (Fig. 2, colored lines).The inner edge of the continuum ring reveals the presence of a second peak located at a deprojected distance of about 0.53 (60 au).The feature is most pronounced along the major axes, which can be explained by the projection effect as well as by the beam whose major axis is oriented roughly along the disk minor axis.Observations at even higher angular resolution are required to quantify this structure in greater detail.
To quantify the radial brightness distribution of the dust ring we use the same approach as Pinilla et al. (2018).We first deproject the data assuming an inclination of 51.7 o , and fit the real part of the deprojected visibilities with a radially asymmetric Gaussian ring using a Markov Chain Monte Carlo (MCMC) method using emcee (Foreman-Mackey et al. 2013).The best fit model has a peak radial position of 73.7±0.1 au, and an inner and outer width of 14.8±0.1 au and 13.4±0.1 au.The ring is therefore radially resolved by our observations.The best fit model is overplotted in Fig. 2   Fig. 2. Radial profiles of the deprojected dust continuum image, along the semi-major (red, orange) and the semi-minor (green, blue) axes, as well as averaged over the entire azimuth (grey).The black line in the uppermost plot corresponds to the best-fit model of the radial profile found in Sect.3.1.1.The deprojection assumes that the continuum is geometrically flat.Radial samples are taken every ∼1/4 beam (20 mas) and the cuts along the minor/major axes are azimuthally averaged in a cone of ±10 o around the corresponding axes.The black arrow highlights a bump in the profile close to the location of PDS 70 b, and the dotted circles mark the location of the second peak.
We confirm an azimuthal brightness enhancement of the ring reported by Long et al. (2018b) on the North-West side of the disk, which peaks at a PA of ∼327 • and which is roughly 13% brighter than the opposite disk side. 2 If the dust is optically thin, the asymmetry could trace the presence of an overdensity.As we will argue below, the dust is likely close to optically thick.The brightness enhancement is therefore likely a combination of both, differences in mass density and temperature.Observations at longer wavelengths are required to break the degeneracies of temperature and density effects, and to conclude on the origin driving the azimuthal brightness asymmetry.

Inner disk
Our image also confirms the detection of a compact signal towards the location of the star, which was already detected and attributed to be a possible inner disk component by Long et al. (2018b) and evidenced by the NIR excess detected in the SED.Our observations marginally resolve the emission inside the innermost ∼80 mas (9 au) at a 5-σ level.Observations at longer wavelengths will enable us to establish the spectral index of this central emission which is required to exclude the possible contribution from free-free emission.
2 Value found by comparing the peak pixel value of the NW side with the peak pixel value of the SE side.

Possible bridge feature
We detect a spur that projects from the dust ring into the gap into the direction of the inner disk at a PA of about 285 o (referred to as 'spur' in Fig. 1 and best seen in panel (d)).This signal is even more clearly detected in the DDT data alone which has a slightly higher resolution (71×56 mas, see Fig. A.5).It is possible that the signal forms a bridge feature connecting the outer and inner disks.Whereas the spur is detected at high confidence (> 5σ), the presence of a continuous connection to the inner disk in the dust continuum remains to be confirmed with deeper observations.Interestingly, this feature is cospatial with an extended feature found in scattered light (Keppler et al. 2018;Müller et al. 2018, see Fig.A.5).Furthermore, the CO shows evidence of a feature at that same location that seems indeed to connect the outer and inner disk (see Sect. 3.2).

Upper limits on CPD dust mass
Figure 2 shows that the radial profile along the SE semi-major axis presents a marginally (SNR ∼3) enhanced signal at ∼0.2 .This corresponds roughly to the expected location of PDS 70 b.We note however that flux density variations of similar amplitude are present at several other position angles as well, and the persistence of this signal is therefore to be tested with deeper observations.
Circumplanetary disks (CPD) are expected to have outer radii R out of a fraction (∼30-70%) of the Hill radius R H (e.g., Quillen & Trilling 1998;D'Angelo et al. 2003;Ayliffe & Bate 2009;Szulágyi et al. 2014), with R H = a P (M P /3M ) 1/3 and a P the planetary distance to the star.For a 5 M Jup companion at 22 au this corresponds to ∼0.8-1.9 au and the disk is therefore expected to be unresolved.Our measured noise level of 0.026 mJy beam −1 translates into a 5-σ upper limit on the flux density of an unresolved CPD around PDS 70 b of 0.130 mJy beam −1 .
We compare this value to the theoretically expected emission from a CPD in order to derive an upper limit on the dust mass.For this aim, we follow the approach presented by Isella et al. (2014), where the dust temperature T d in the CPD at a given radius r from the planet is described as: where T irr, is the temperature of the surrounding circumstellar disk heated by the central star at the distance of the planet to the star, T irr,p is the temperature due to the heating by the planet itself, and T acc denotes the contribution from viscous accretion within the CPD.
For T irr, we adopt a value of 19 K at a distance of 22 au from the star which is estimated from our RT models (Keppler et al. 2018).The irradiation by the planet, T irr,p , can be estimated (assuming a CPD aspect ratio of 0.1; Zhu et al. 2018) as: where we use L p ∼1.5×10 −4 L as the luminosity of PDS 70 b (Müller et al. 2018).Finally, the heating due to accreting material is given by: with Ṁacc the mass accretion rate onto the planet and r p the planetary radius.From Wagner et al. (2018), we assume Ṁacc ∼ 10 −8 M Jup yr −1 , and r P ∼ 3 R Jup (Müller et al. 2018).
As Isella et al. (2014), we assume a power-law surface density Σ(r) = C × r −3/4 , where C is the normalisation constant for the total CPD dust mass M d = r out r in Σ(r)2πrdr.We therefore can compute the expected mm-flux F d for a given M d by integrating the flux density contribution from each radius over the entire CPD: Here, κ denotes the dust opacity which we assume to be 3.5 cm 2 g −1 at 855 µm, linearly scaled from Andrews et al. (2012), B ν the Planck function evaluated at T d , and i the CPD inclination, which we assume to be equal to the inclination of the circumstellar disk (51.7 o ).
We compute the expected flux densities for different CPD dust masses considering outer CPD radii of 0.3-0.7 r H and assuming that the CPD touches the planetary surface (e.g., r in = r p , but note that regions in which the temperature exceeds the sublimation temperature of silicates (∼1500 K) are taken out of the integral).The result is compared to our noise level of 0.026 mJy beam −1 and shown in Fig. 3.With the given choice of parameters, we find a 5-σ upper dust mass limit of ∼0.01 M ⊕ (∼0.8 lunar masses).This value is roughly independent from the outer CPD radius, which means that the emission is likely optically thin.As shown in Appendix A.1, this detection limit holds for the entire estimated mass range of PDS 70 b.3.2.12 CO J = 3 − 2 Figures 1 (a) and (c) show the 12 CO J = 3 − 2 integrated intensity (zeroth moment) map, with the latter including annotations of the major features.There is a clear asymmetry with respect to the disk major axis.This is due to the significantly elevated τ ∼ 1 surface of the 12 CO, which is typically assumed to trace disk layers where z /r ∼ 0.25 (Rosenfeld et al. 2013).In addition, several other features are seen, including two gaps (a prominent one at ∼ 0.2 and a faint one at ∼ 0.6 ), a bridge-like feature similar to the one observed in the continuum, and apparent shadowing along the major and minor axes, previously reported by Long et al. (2018b).
Towards the center of the image, the inner disk component is clearly detected, extending out to approximately 15 au, which is consistent with estimates from scattered light (Keppler et al. 2018).For disks shaped by planets, the presence of a bright gaseous inner disk (implying a gas gap rather than a cavity) is in agreement with the predictions from hydrodynamical models, even for the cases where the planet's mass is as high as 10 M Jup (Facchini et al. 2018).
Just at about the same location as the spur found in the continuum, the zeroth moment map shows evidence of extended signal connecting the inner disk and the outer ring in the North-West region.This signal may be connected to the extended feature detected in the NIR (Keppler et al. 2018;Müller et al. 2018, and Fig. A.5, right panel, of this paper), and may also possibly be related to the features seen in CO and HCO + by Long et al. (2018b) at similar locations.In case this feature is indeed connecting the outer and inner disks, it may be tracing gas flow through the gap from the outer to the inner disk (e.g., Tang et al. 2017;Casassus et al. 2015;Price et al. 2018).Such a hypothesis could be confirmed through the detection of localised velocity changes in the given region, which we however do not detect with the present spectral resolution.The nature of this feature needs therefore to be tested with observations at higher spectral and angular resolution.
The inner gap at ∼0.2 is likely due to a gap opened by PDS 70 b and will be discussed further in Section 4.1.The outer gap at ∼0.6 can be explained by continuum absorption of the bottom side of the disk: as shown in Fig. 4, the contours of equal projected velocity at the top and bottom sides of the disk in regions between the disk major and minor axes are spatially offset.While emission from the bottom side travels through the midplane towards the observer, it is absorbed by the dust, reducing the integrated flux at that location (e.g., Isella et al. 2018).As emission from the bottom side of the disk is nearly totally absorbed, we conclude that the dust ring is likely optically thick at ν = 345 GHz, a result which is found at mm-wavelengths for other disks as well (e.g., Pinilla et al. 2017).
Along the disk major and minor axes, on the other hand, the iso-velocity contours do overlap.Because the 12 CO is optically thick, emission from the bottom side of the disk is self-absorbed and only the top side is visible.This causes the apparent shadowing along the major and minor axes of the disk (and the shadowing observed in the HCO + data presented by Long et al. 2018b).A more elevated emission layer results in a larger azimuthal variance, as the two sides become more spatially resolved.The difference between the value along an inter-axis region and along an axis will peak at roughly a factor of two, a feature which is commonly seen in the integrated intensity maps of high spatial resolution observations of 12 CO (e.g., Rosenfeld et al. 2013).

Deriving a 12 CO Emission Surface
Since the 12 CO emission is coming from an elevated layer above the midplane, we need to deproject the data in order to precisely analyse the emission and velocity structure as a function of the radius.For this aim, we want to derive constraints on the emission height of the 12 CO.Following Teague et al. (2018b), we generate a map of the rotation velocity using the method presented in Teague & Foreman-Mackey (2018), which is robust against confusion from the near and far sides of the disk 3 .We then fit a Keplerian rotation pattern to the data, including a flared emission surface parameterized as z(r) = z 0 × (r / 1 ) ϕ and fix v LSR ± 8 ∆v chan Fig. 4. The iso-velocity contours for the upper (blue) and lower (red) sides of the disk at different velocities with a flared emission surface.Along the major and minor axes, shown by the black dotted lines, the iso-velocity contours overlap, as in the left most and right most panels, and thus only emission from the upper side of the disk is visible.Conversely, in inter-axis regions the iso-velocity contours are spatially separated, as in the central panels, so that emission from both sides of the disk reaches the observer.Based on Fig. 4 from Rosenfeld et al. (2013).the inclination at i = 51.7 • to break the degeneracy with the stellar mass.We note that our modeling of the surface height is limited to a generic model of a flared surface due to the limited resolution of our data.To perform more detailed modeling of the emission surface under consideration of spatial variations of the underlying gas density structure higher resolution is required.Our modeling results in a tight constraint on the emission surface of: (5) with the additional parameters of M star = 0.875 ± 0.03 M sun , PA = 160.4• ± 0.1 • and v LSR = 5505 ± 2 m s −1 .These uncertainties describe the 16th to 84th percentile range of the posterior distributions for each parameter which are symmetric about the median.We note that these uncertainties correspond to the statistical uncertainties and do not take into account the systematic uncertainties that may be significantly larger.Figure 5 shows the best fit emission surface overlaid on the rotation map.Using this emission surface, the data is deprojected into bins of constant radius and azimuthally averaged with the resulting integrated intensity profiles shown in Fig. 6.The radial profile of the integrated flux density, top panel, shows a clear gap at 0.2 (∼23 au), consistent with the orbit of PDS 70b (Keppler et al. 2018;Müller et al. 2018) and a gap width of ∼ 0.1 .Due to the very high optical depth of 12 CO, any visible gap feature requires a significant depletion of gas or considerable change in gas temperature (e.g., Facchini et al. 2018).
Using the brightness temperature, T B , presented in Fig. 6 (lower panel) as a proxy of the gas temperature, we infer a drop in the local gas temperature across the gap.This is consistent with a surface density depletion of the gas, which would move the τ = 1 surface of the 12 CO deeper within the disk, closer to the cooler midplane, therefore dropping the temperature.One possibility to clearly disentangle the effects of temperature and density on the brightness temperature is to use the CO line width as a tracer for temperature variations (Teague et al. 2018a), for which however higher spectral resolution would be required than given by our data.
From the integrated flux density profile, we find that the gap extends from about 0.1 to 0.3 (∼11 to 34 au).It is spatially resolved, and does not seem to extend out to the location of the dust continuum ring, although it is not possible to measure 12 CO depletion accurately due to its large optical depth.This preferential depletion of grains compared to gas within a cavity is a common feature for transition disks (van der Marel et al. 2015Marel et al. , 2016)).3.2.2. 12 CO Rotation Curve Radial gas pressure gradients perturb the gas rotation velocity and are used as tracers for planet-induced perturbations (Pérez et al. 2015;Pinte et al. 2018;Teague et al. 2018a).Velocity distortions by the planet at the close-in location of 22 au are small in size, such that their detection in single channel maps as done by Pinte et al. (2018) is hampered by our limited angular and spectral resolution (see also Sect.4.1.2),and further by the relatively low signal-to-noise (SNR) of the CO emission at the location of the planet (well within the CO integrated flux density gap).To improve the SNR of potential kinematic perturbations we therefore make use of an azimuthally averaged rotation curve of the 12 CO data to probe the underlying gas density structure (Teague et al. 2018a).This is even possible in cases where the line emission is optically thick.Whereas a negative pressure gradient induces sub-Keplerian rotation, a positive pressure gradient would cause super-Keplerian rotation.
Following the method described in Teague et al. (2018b), we infer the rotation profile by finding the rotation velocity for each radius which allows for all spectra in an annulus to be shifted back to the same systemic velocity 4 .We run ten different realizations of this, randomizing the pixels taken from each annulus (making sure they are separated by at least one FWHM of the beam), and randomizing the radial locations of the annuli while maintaining a radial bin width of a quarter beam width.The resulting rotation curve and the residual relative to the best-fit Keplerian profile are plotted in Fig. 7.
The absolute scale of the deviation from Keplerian rotation depends on the reference Keplerian velocity and therefore on the assumed stellar mass.The systematic uncertainties on the dynamical determination of the stellar mass as well as the parametrization of the surface together with the fact that our fiducial model for the rotation velocity does not take into account the overall pressure gradient in the disk may cause that the uncertainty of the absolute scaling is as large as 10%.Figure 7 (bottom panel) shows the residuals of the rotation curve (blue), where the green hatched area marks the uncertainty of the zero point of δv rot inferred by the 3-σ statistical uncertainties of the stellar mass.Within these uncertainties, the peak of the continuum ring (∼0.65 ) lies close to the location where δv rot recovers Keplerian rotation and therefore where pressure reaches its maximum.
A significant deviation of up to ∼ 12% at ∼ 0.2 is observed, suggestive of significant changes in the gas pressure at this location, consistent with the structure observed in the rotation map in Fig. 5.The rotation curve clearly demonstrates a positive pressure gradient between ∼0.4 and 0.8 , reaching a maximum at about 0.55 .This implies that the gas density is likely depleted beyond ∼0.4 , and suggests therefore that the gap is in reality larger than what is observed in integrated emission: if the gap was only as wide as the gap in the 12 CO integrated emission then one would expect the peak residual of the rotation curve to fall at the edge of the gap at ∼ 0.3 (see Fig. 1 in Teague et al. 2018a, for example), however the peak is found closer to 0.55 .The shape of the residual curve in the inner disk, r < 0.3 , is dominated by the steep gradients in the intensity profile due to both the inner disk and the gap, challenging a direct analysis.A more thorough discussion of this effect and the impact of beam smearing is discussed in Sect.4.1.2in the context of hydrodynamical models and in Appendix A.2.

Potential Point Source
We tentatively detect a point source in the 12 CO emission maps at a projected separation of ∼0.39 and a PA of ∼260 o .This corresponds to a deprojected radius of ∼71 au, if coming from the midplane.The peak is detected at a ∼ 6σ level and spatially offset from the Keplerian emission pattern.Figure 8 shows the spectrum extracted at the location of the source, and three channel maps showing the offset nature of the emission.The signal appears at a velocity of around 6.45 km s −1 , corresponding to a redshift of roughly 1 km/s with respect to the line center of the Keplerian profile.The spectrum also shows a blueshifted peak which emission, however, may be biased from the bottom side of the disk.Interestingly, if located in the midplane, the source would be located well within the dust continuum ring, close to the dip between the main and the tentative second peak detected in the continuum profiles (see Sect. 3.1.1).Spatially offset emission has been shown to be potentially a signature of a CPD (Pérez et al. 2015), as the additional rotation of the CPD would shift the emission from the Keplerian pattern.If the signal were indeed connected to a forming, embedded planet this may explain the azimuthal gap found in the HCO + emission at a similar location (Long et al. 2018b), since chemical changes due to heating from the planet may locally deplete HCO (Cleeves et al. 2015).Additional observations are required to confirm the potential point source.

Discussion
As shown by theoretical studies, the interaction of a massive body with the disk opens a gap in the gas (e.g., Lin & Papaloizou 1986).The perturbation of the local gas density causes a change in the local pressure gradient, which manifests itself in two ways.First, it generates a pressure bump outside the planetary orbit, trapping large dust particles (while small ones, well coupled to the gas, may still enter the gap).This leads to a spatial segregation of large and small grains (e.g., Pinilla et al. 2012).Second, the change in pressure gradient manifests itself in a local deviation from Keplerian rotation the amplitude of which is sensitive to the planet's mass (Teague et al. 2018a).
Our aim is to investigate the impact of PDS 70 b on the observed disk morphology.For this purpose we carried out hydrodynamic and radiative transfer simulations that we will present in the next section.

Model setup
To simulate interaction between PDS 70 b and the circumstellar disk of PDS 70, we carry out three-dimensional hydrodynamic calculations using FARGO3D (Benítez-Llambay & Masset 2016; Masset 2000).We adopt the disk density and aspect ratio profiles used in Keppler et al. (2018): and where R c = 40 au, R p = 22 au is the distance of PDS 70 b assuming a circular orbit, (H/R) p = 0.089 and f = 0.25.Σ c = 2.87 g cm −2 is chosen such that the total gas mass in the disk is 0.003 M , consistent with the model presented in Keppler et al. (2018).The surface density profiles are shown in Fig. 9 (a).We assume a vertically isothermal disk temperature structure and use an isothermal equation of state.
The simulation domain extends from r = 0.2 R p to 9 R p in the radial direction, from π/2 − 0.4 to π/2 in the meridional direction, and from 0 to 2π in the azimuthal direction.We adopt 256 logarithmically space grid cells in the radial direction, 48 uniformly spaced grid cells in the meridional direction, and 420 uniformly spaced grid cells in the azimuthal direction.A disk viscosity of α = 10 −3 is added to the simulations.This value of turbulence is consistent with the level of turbulence constrained for the protoplanetary disks around TW Hya (Teague et al. 2016(Teague et al. , 2018c;;Flaherty et al. 2018) and HD 163296 (Flaherty et al. 2015(Flaherty et al. , 2017)).
We test three planet masses: 2, 5, and 10 M Jup , covering the range of potential planet masses proposed by Keppler et al. (2018), assuming a 0.85 solar-mass star.The simulations ran for 1000 orbits, after which we find the gap width and depth reach a quasi-steady state.This is in agreement with other planetdisk interaction simulations from the literature (e.g., Duffell & MacFadyen 2013;Fung et al. 2014;Kanagawa et al. 2015).The radial profile of the deviations from Keplerian rotation after 1000 orbits are shown in Fig. 10 (b).
We generate 12 CO image cubes using the radiative transfer code RADMC3D version 0.415 .We first compute the disk thermal structure by running a thermal Monte Carlo calculation.To do so, we place a 0.85 solar-mass star at the center, which has an effective temperature of 3972 K and a radius of 1.26 R (Pecaut & Mamajek 2016;Keppler et al. 2018), emitting 10 8 photon packages.As in Keppler et al. (2018), we consider two grain size distributions whose number density follows a power-law as a function of the grain size a with n(a) ∝ a −3.5 : small grains range from 0.001 to 0.15 µm and large grains range from 0.15 to 1000 µm in size.The relative mass fraction of small to large grains is 1/31, implying that about 3% of the total dust mass is confined within the small grain population.This is consistent with previous RT models of PDS 70 (Dong et al. 2012;Keppler et al. 2018).We assume that the grains are composed of 70 % astronomical silicates (Draine 2003) and 30 % amorphous carbon grains (Zubko et al. 1996).The grain opacity is computed according to the Mie theory using the BHMIE code (Bohren & Huffman 1983).
CO line radiative transfer is done under LTE assumptions, assuming a constant 12 CO to H 2 ratio of 10 −4 (e.g., Lacy et al. 1994;Williams & Best 2014).A local, spatially unresolved microturbulence is added at a constant level of 30 m s −1 .This choice is equivalent to α of a few ×10 −3 .We simulate the ALMA observations using the SIMOBSERVE task in CASA version 5.1.2.using the same velocity resolution, synthesized beam, and onsource integration time to those used in the observations.Thermal noise from the atmosphere and from the antenna receivers is added by setting the thermalnoise option in the simobserve task to tsys-atm.Using the same tools as for the observations, we derived the velocity integrated flux density, as well as the rotation profiles for each simulation (Figures 9 and 10).

Comparison with observations
The disk density distribution from the hydrodynamic model and a simulated 12 CO zeroth moment map are presented in Fig. 11.The 5 M Jup planet opens a gap around its orbit, which is clearly seen in the simulated zeroth moment map.We find that velocity kinks associated with the planet-driven spiral arms are present in raw simulated channel maps, similar to what is found in HD 163296 (Pinte et al. 2018).However, the velocity distortions are too small in size and thus smeared out after convolution with the ALMA beam.
We compare the radial profiles of the simulated and observed integrated flux densities in Figure 9 (b).The profiles show evidence of a depletion in integrated flux density at the location of the planet, which is stronger for larger planet masses.Both the width and depth of the depleted flux density in the observations are reasonably well reproduced by a 5 M Jup planet.We note that the models appear to overestimate the increase of flux density towards the inner disk.Since CO is optically thick, this is likely caused by a different temperature structure of the inner disk region, the reason for which could be a different density profile than assumed (e.g., overestimation of the actual density in the inner part of the disk, or a different gap shape), but needs further investigation with higher angular resolution that is able to better resolve these inner regions.
Figure 10 (panel a) presents the absolute rotation profiles while the residual δv rot profiles before and after radiative transfer and ALMA simulations are shown in panels b and c, respectively.We note two points: first, comparison of the residual model profiles before and after convolution alters the overall shape of the rotation curves.Second, the residual curve of the PDS 70 disk follows the general shape of the modelled curves, but differs with respect to the location of the maximum, as well as the velocity gradient towards the inner disk.
The change in the shape of the rotation curve when simulating the observations is due to beam convolution effects in the presence of strong radial gradients in intensity and velocity.This is described in detail in Appendix A.2.In brief, sharp edges in the flux density profile induce a distortion in the measurement of the rotation curve because the velocities measured within one beam are biased towards those at the highest line intensity.This causes the velocity to be overestimated in the inner region of the gap, and underestimated at the outer edge of the gap.The resulting rotation curve is of characteristic shape which is asymmetric with respect to the gap center (see Fig. A.2).It shows evidence of strong super-Keplerian rotation in the inner gap region, and a less strongly pronounced region of sub-Keplerian rotation at the outer gap edge.This effect is now superimposed on the effect of the planet-induced pressure gradient on the rotation profile (sub/super-Keplerian rotation inside/outside the planet's orbit).This effect can be fully accounted for when performing forward modelling.As Fig. 10 (c) shows, all convolved model profiles show this characteristic shape with the amplitudes of their minima and maxima depending on the planet's mass.
The observed rotation curve of PDS 70 shows the same characteristic transition from sub-Keplerian to super-Keplerian transition as the models.While we found that the width and depth of the integrated flux density profiles seem consistent with the effect of a 5 M Jup planet, we find that the radial location and the amplitude of the minimum δv rot of the rotation curve of the PDS 70 disk is best matched by the perturbations created by a 10 M Jup planet.We note however that our hydrodynamic models consider a vertically isothermal temperature structure, whereas in a more realistic approach (introducing a more physical prescription for the vertical temperature structure) the deviation from Keplerian rotation may actually be higher in the disk surface than in the midplane, implying that the δv rot in the current models may be underestimated (Bae et al in prep.; see also Fig. 3 of Teague et al. 2018a).Relaxing the isothermal assumption and introducing a more physical prescription of the vertical temperature structure may be able to solve this discrepancy, but is beyond the scope of this study.
Towards the inner region the observed rotation curve is less steep, which, again may be due to a slightly different gap shape (i.e., a less steep inner edge).The most conspicuous difference to the models is the region of super-Keplerian rotation beyond the planet which is extending far more out than in the models.As seen in Sect.3.2.2,within the uncertainties, the observed rotation curve returns to Keplerian rotation close to the location of maximum emission in the continuum ring (∼0.65 or 74 au) (see Fig. 7).This is consistent with the interpretation of large grains being trapped in the region of maximum pressure (e.g., Pinilla et al. 2012).While we have shown that the observed integrated flux density profile can be reproduced well by one planet of 5 M Jup , the large extension of super-Keplerian rotation and the concomitant far-out location of the continuum ring imply that the gap is in reality wider than predicted by all the models.It therefore appears -within our model assumptions -that only one planet located at the orbit of PDS 70 b may not be sufficient to generate a kinematic signature in the disk with the inferred width, nor maintain a continuum ring at ∼74 au, a scenario that needs to be probed by future observations at higher spectral resolution.This is consistent with gap width considerations in the literature.As an example, hydrodynamical and dust simulations suggest that the accumulation of large dust grains is expected to be found at roughly 10 R H outwards of the planet's orbit (Pinilla et al. 2012;Rosotti et al. 2016).For a 10 M Jup planet at the location of PDS 70 b's orbit, the dust ring would therefore be expected at about 46-56 au, assuming a stellar mass of 0.88 M .This suggests that an additional low-mass planet located beyond PDS 70 or the combination with other physical mecha-nisms such as photoevaporation or dead zones may be needed to explain the outwards shifted location of the pressure bump.Indeed, models predict that large gaps in transitional disks can be reproduced by introducing multiple planets (Dodson-Robinson & Salyk 2011;Zhu et al. 2011;Duffell & Dong 2015).Detailed modeling of the system introducing multiple planets, as well as deep observations are required to constrain the planetary architecture responsible for the observed features, which is beyond the scope of this study.
An alternative scenario to explain the distant location of the ring compared to the position of PDS 70 b is to consider that the ring is tracing a secondary pressure bump.Indeed, single planets can open multiple gaps in a disk with low viscosity or alternatively, vortices generated at the edge of a gap can lead to a secondary ring (Lobo Gomes et al. 2015).In that latter scenario, the primary ring, located at ∼50 au (corresponding to ∼10 R H from a 5 M Jup planet at 22 au), would be depleted.The secondary ring would be located at ∼1.5× the location of the primary ring (Lobo Gomes et al. 2015), corresponding to ∼75 au, which is where the dust ring is found in the PDS 70 system.Furthermore, secondary vortices may be generated at the edge of the secondary ring.If this is the case, this may also explain the azimuthal asymmetry observed in the dust continuum.A detailed exploration of this scenario will be the subject of a follow-up study.

Upper limit on CPD dust mass
The detection of Hα emission at the location PDS 70 implies that PDS 70 b is actively accreting (Wagner et al. 2018) and therefore likely posseses an accretion disk.Still, we can only derive upper limits on the circumplanetary disk with our data.Models of planet formation predict the presence of circumplanetary dust around young planets, implying that CPDs should be frequent.However, searches for circumplanetary material in the submillimeter/millimeter continuum around other young substellar companions have been unsuccessful though the detection of active accretion through the Hα and/or Paβ lines was seen in some of these cases (e.g., Isella et al. 2014;Bowler et al. 2015;MacGregor et al. 2017;Wolff et al. 2017;Ricci et al. 2017;Pineda et al. 2019).Our upper limit on the CPD dust content of ∼0.01 M ⊕ is similar to that derived for other systems (Pineda et al. 2019).
There are several reasons for which the detection of CPD's in the (sub-)millimeter regime may be challenging.First, CPD's are expected to be very small, substantially reducing the emitting area and therefore the expected signal.Second, since the large grains are substantially trapped in the outer dust ring, the replenishement of large grains within the gap is expected to be inefficient.Even if small grains go trough the gap, and replenish the CPD, once they grow the radial drift is expected to be extremely efficient, depleting the large grains very fast (Pinilla et al. 2013;Zhu et al. 2018).The search for the CPD using gas kinematics as a tracer or NIR observations might therefore be more promising.

Summary and conclusions
PDS 70 b is the most robust case of a directly imaged young, forming planet in the gap of a transition disk.We obtained ALMA band 7 DDT observations in Cycle 5 and combined them with previous Cycle 3 data (Long et al. 2018b) to study the planet's natal environment at high angular resolution (∼0.07 ) in dust continuum and at the 12 CO J=3-2 transition.
-We detect the emission from the dust continuum as a highly structured ring.Its radial distribution peaks at ∼74 au.The inner edge of the ring shows evidence of a marginally resolved second ring component, peaking at around 60 au.We also detect a spur projecting into the gap at a PA of about 285 o and confirm an azimuthal brightness asymmetry with an brightness enhancement of about 13% in the North-West part of the ring.-We can derive upper limits on the circumplanetary disk.
From the noise level of the image we infer a 5-σ upper dust limit of less than ∼0.01 M ⊕ .-The CO integrated intensity shows evidence of two radial intensity depression, with the inner depression of flux density at ∼0.2 (corresponding to the location of PDS 70 b) and a second gap at about 0.6 .The inner gap is most likely carved by PDS 70 b.Comparison of the flux density profile to hydrodynamical simulations show that the gap width and depth is best reproduced by a 5 M Jup body.The outer gap can be explained by the dust being optically thick.Further, we find evidence for an azimuthal intensity modulation which is due to self-absorption by optically thick CO.We also detect a bridge-like feature in the CO at the location of the spur seen in the continuum, as well as the inner disk which extends out to ∼15 au.Finally, we report on the tentative detection of a possible point source in the 12 CO emission maps, the existence of which needs to be confirmed with additional observations.-We detect significant deviation from Keplerian rotation inside ∼0.8 .The width of the δv rot feature is consistent with the far-out location of the dust ring.Comparison to hydrodynamical simulations implies that the depth of the kinematic signature is best matched by a ∼10 M Jup object (within our model assumptions of an isothermal disk), but the width of the feature suggests that only one planet located at the orbit of PDS 70 b may not be sufficient to generate a gap with the inferred extension.An additional physical mechanisms or a second low-mass body may be required to explain the disk morphology.Future observations at higher angular and spectral resolution will allow us to put tighter constraints on the planetary system architecture that can account for all of the observed features in the PDS 70 disk morphology.across the beam is steep, this will cause the sampled velocity to be strongly biased towards the region of highest intensity, rather than the beam center.This effect is illustrated in Fig. A.2 (left).The figure shows, assuming a smooth radial disk intensity profile, for each distance the sampled radius, i.e., the radius within the beam at which the velocity receives the highest weighting and which therefore corresponds to the effective radius at which the velocity is observed.This is shown for different power-law exponents of the disk radial intensity profile.We note that the steeper the intensity profile, the larger the sampled radius is biased towards smaller radii, and the more the measured velocity will be overestimated.This is even more complex if the intensity profile deviates from a simple power law, as in the presence of a gap structure.The additional steep gradients on the gap edges cause regions closer to the inner gap edge being even more biased towards smaller distances and therefore higher velocities, whereas regions close to the outer gap edge are biased towards larger distances, and smaller velocities.Figure A.2 (bottom right) shows the deviation from Keplerian rotation, assuming an intensity profile with a gap structure centered around 0.2 and a beam size of 76 mas (top right).The resulting δv rot profile is asymmetric with respect to the gap center, with super-Keplerian rotation in the inner regions going into sub-Keplerian rotation beyond ∼ 0.3 , and the strength of the deviation is sensitive on the gap depth.This beam smearing effect is additonal to the deviation from Keplerian rotation due to the planet-induced pressure gradient.
Figure A.3 demonstrates the effects of this bias using the radial intensity profile from Fig. 6 and compares this effect to the functional form from the pressure gradient (shown in light blue).The resulting profile will be the combination of both factors with their relative amplitudes dependent on the gap shape.While this limits interpretation, these effects will be full accounted for with forward modelling, as presented in Section 4.1.

Fig. 1 .
Fig.1.Observations of the 12 CO (left column) and the 350.6 GHz continuum (right column).The bottom row provides a closer view of the observations including annotations where the color scaling has been stretched to bring out detail.The contours for the 12 CO are starting at 20% of the peak value to the peak in steps of 10%.For the continuum, the gray dashed contour is 5σ with black contours starting at 10σ and increasing in steps of 10σ where σ = 26 µJy beam −1 .The synthesized beams are shown in the bottom left of each panel.
(black line) and is shown in Fig. A.4.

Fig. 3 .
Fig. 3. Theoretically expected flux densities from a CPD around a 5 M Jup planet at the location of PDS 70 b with different dust masses and outer disk radii, following the prescription from Isella et al. (2014).The contours mark the 2, 3, and 5-σ detection limits from the observations.

Fig. 5 .
Fig. 5. Rotation profile of the 12 CO emission, left, using the method presented in Teague & Foreman-Mackey (2018), with the best fit surface overlaid, right.The solid lines show the top surface, while the dotted lines show the far surface.

Fig. 6 .
Fig. 6.Radial profiles of the 12 CO integrated intensity, top, and brightness temperature, bottom.Radial samples are taken every 1/4 beam and the error bar shows the standard deviation in the azimuthal bin.The vertical dotted line shows the orbit of PDS 70b, while the gray shaded region shows the extent of the continuum ring.The beam size is shown in the top right of each panel.

Fig. 7 .
Fig. 7. Top: The measured rotation curve with 1σ uncertainties.The blue line and blue shadowed area show the running mean and its standard deviation.The dashed grey lines show the Keplerian rotation curve assuming the best-fit stellar mass (0.88 M , thick) and including the 3σ uncertainties on the stellar mass (corresponding to 0.79 and 0.97 M respectively, thin) derived from the rotation map fitting.Note that the uncertainties of the stellar mass correspond to the statistical uncertainties and do not include the systematics.Bottom: relative residuals (blue solid) and uncertainties (blue shaded area) between a smooth Keplerian curve and the inferred rotation curve.The green hatched area highlights the uncertainty of the absolute scaling of δv rot inferred by the 3-σ statistical uncertainties on the stellar mass.In both panels, the gray shaded region shows the extent of the continuum ring.The vertical dotted line show the orbit of PDS 70 b and the shaded vertical gray region traces the location of the continuum emission.

Fig. 8 .
Fig. 8. Top panel: CO line profile extracted at the location of the point source.Bottom panel: CO channel maps around 6.45 km/s.The white circle indicates the location of the point source.

Fig. 9 .
Fig. 9. Comparison of hydrodynamical models including a 2 M Jup (yellow), 5 M Jup (green), and 10 M Jup (red) planet located at 0.2 with the observations (blue).(a): Azimuthally averaged surface density profiles of hydrodynamical simulations.The dotted line corresponds to the initial, unperturbed surface density profile.(b): Integrated, azimuthally averaged CO flux density of observations and ALMA-simulated models, after applying 2-σ clipping.In each panel, the grey shaded area indicates the extension of the continuum ring, and the vertical dotted line corresponds to the approximate location of PDS 70 b.The black bar in the second panel indicates the major axis of the beam (0.076 ).

Fig. 10 .
Fig. 10.Comparison of hydrodynamical models including a 2 M Jup (yellow), 5 M Jup (green), and 10 M Jup (red) planet located at 0.2 with the observations (blue).(a): Rotation velocity as a function of deprojected distance.The grey dash-dotted line indicates the unperturbed Keplerian profile around a 0.88 M star.(b): Deviation from Keplerian rotation of the hydrodynamical simulations at the τ=1 surface.(c): Deviation from Keplerian rotation of ALMA-simulated models and observations.The plot shows the running mean and standard deviations.The inner region up to 160 mas is affected by beam confusion effects and is therefore blocked out.In each panel, the grey shaded area indicates the extension of the continuum ring, and the vertical dotted line corresponds to the approximate location of PDS 70 b.The black bar in the first and third panel indicates the major axis of the beam (0.076 ).

Fig. 11 .
Fig. 11.Left: 3D volume rendering of the gas density (in a normalized unit with logarithmic scaling) after evolution of 1000 orbits in the inner 100 au of the model disk with a 5 M Jup planet at 22 au.Ticks on the axes mark every 25 au.Right: A simulated 12 CO zeroth moment map based on the hydrodynamic model presented in the left panel.