The Disk Substructures at High Angular Resolution Project (DSHARP) - IX. A high definition study of the HD 163296 planet forming disk

TL;DR

This high-resolution ALMA study of the HD 163296 protoplanetary disk reveals detailed dust and gas structures, including rings, gaps, and asymmetries, providing insights into planet formation processes and disk-planet interactions.

Contribution

The paper presents the first high-resolution imaging of dust and gas in HD 163296, revealing new small-scale structures and constraining disk properties, advancing understanding of planet formation environments.

Findings

Resolved dust rings at 67 and 100 au.

Discovered new structures: a gap at 10 au, a ring at 15 au, and a crescent at 55 au.

Observed asymmetries consistent with planet-disk interactions.

Abstract

ALMA observations of protoplanetary disks acquired by the Disk Substructure at High Angular Resolution Project (DSHARP) resolve the dust and gas emission on angular scales as small as 3 astronomical units, offering an unprecedented detailed view of the environment where planets form. In this article, we present and discuss observations of the HD 163296 protoplanetary disk that imaged the 1.25 mm dust continuum and $^{12}$CO J=2-1 rotational line emission at a spatial resolution of 4 and 10 au, respectively. The continuum observations resolve and allow us to characterize the previously discovered dust rings at radii of 67 and 100 au. They also reveal new small scale structures, such as a dark gap at 10 au, a bright ring at 15 au, a dust crescent at a radius of 55 au, and several fainter azimuthal asymmetries. The observations of the CO and dust emission inform about the vertical structure of the disk and allow us to directly constrain the dust extinction optical depth at the dust rings. Furthermore, the observed asymmetries in the dust continuum emission corroborate to the hypothesis that the complex structure of the HD 163296 disk is the result of the gravitational interaction with yet unseen planets.

Tables

Table 1. Table 1: Results of elliptical fitting to the bright and dark rings observed in the 1.25 mm continuum map.

Feature	$\Delta x_0$ (mas)	$\Delta y_0$(mas)	$r$(”)	$r$(au)	$i$(°)	PA(°)
D10	-6.2$\pm$0.8	10.5$\pm$0.7	0.0987$\pm$0.0009	9.96$\pm$0.07	37.9$\pm$1.2	127.90$\pm$1.9
B14	-4.8$\pm$0.7	11.4$\pm$0.7	0.1430$\pm$0.0010	14.44$\pm$0.07	47.24$\pm$0.64	131.10$\pm$0.88
D45	-5.4$\pm$2.0	5.2$\pm$1.9	0.4433$\pm$0.0026	44.77$\pm$0.19	42.22$\pm$0.67	133.67$\pm$0.95
B67	-5.3$\pm$1.0	7.3$\pm$1.1	0.6633$\pm$0.0014	66.99$\pm$0.11	46.78$\pm$0.21	133.13$\pm$0.29
D86	6.7$\pm$2.3	-0.3$\pm$2.3	0.8575$\pm$0.0028	86.61$\pm$0.22	47.34$\pm$0.32	132.78$\pm$0.45
B100	-2.3$\pm$0.9	8.6$\pm$0.9	0.9870$\pm$0.0011	99.69$\pm$0.08	46.59$\pm$0.11	133.46$\pm$0.15
D141	-0.2$\pm$9.9	6.4$\pm$9.5	1.3923$\pm$0.0126	140.62$\pm$0.96	47.2$\pm$0.9	131.3$\pm$1.2
B159	-16.2$\pm$12.0	1.1$\pm$12.0	1.572$\pm$0.016	158.7$\pm$1.2	45.7$\pm$1.0	132.0$\pm$1.4
mean*	-3.5$\pm$0.5	9.1$\pm$0.5	-	-	46.7$\pm$0.1	133.3$\pm$0.1

Table 2. Table 2: Intensity ratios between bright rings and dark gaps along different radial directions

Pair	+90°	-90°	0°	180°	Mean
	(Maj E)	(Maj W)	(Min N)	(Min S)
B14/D10	1.2	1.1	1.0	0.9	1.1
B67/D45	13.3	$>$37	6.6	4.7	$>$17
B100/D86	11.7	7.1	4.5	2.5	7.0
B159/D141	2.1	1.1	1.0	0.5	1.3

Table 3. Table 3: Results of Gaussian model fitting of the dust continuum emission

	—— Radial intensity profile ——			—— Visibilities ——
Feature	$r$ (au)	$w_d$ (au)	I (Jy/as2)	$r$ (au)	$w_d$ (au)	I (Jy/as2)
-	-	-	-	3.8$\pm$0.1	2.3$\pm$0.1	$1.7\pm 0.2$
B14	-	-	-	$15.5\pm 0.2$	$8.7\pm 0.2$	1.1$\pm$0.1
-	-	-	-	31.8$\pm$0.3	4.4$\pm$0.1	0.62$\pm$0.06
B67	67.7	6.84	0.38	67.08$\pm$0.03	6.56$\pm$0.05	0.38$\pm$0.02
B100	100.0	4.67	0.24	101.16$\pm$0.04	5.8$\pm$0.1	0.20$\pm$0.01

Equations

$$T_b^{pl}=\frac{h\nu }{k_b}{\left[\ln \left(\frac{2h{\nu }^3\Omega }{c^2{F}_{\nu }}+1\right)\right]}^{-1},$$

$$I(r)=\sum _{i=1}^n{I}_i{e}^{-(r-r_i{)}^2/2w_{d,i}^2}.$$

$${\tau }_d^{abs}=ϵ\ln \frac{I_{co}}{I_{co}^{obs}}$$

$$T_d=\frac{h\nu }{k_b}{\left[\ln \left(\frac{2h{\nu }^3\Omega (1-e^{-{\tau }_d^{abs}})}{c^2{F}_{\nu }}+1\right)\right]}^{-1},$$

$$T_d=\frac{h\nu }{k_b}{\left[\ln \left(\frac{2h{\nu }^3\Omega (1-I_{co}^{obs}/I_{co})}{c^2{F}_{\nu }}+1\right)\right]}^{-1}.$$

$$\delta \mathrm{RA}$$

$$\delta \mathrm{Dec}$$

$$x^{\prime }$$

$$y^{\prime }$$

Full text

protoplanetary disks, planet-disk interactions, submillimeter: planetary systems

1 Introduction

In recent years, millimeter-wave interferometers and near infrared high contrast cameras have imaged nearby protoplanetary systems at unprecendented angular resolution in both continuum and line emission (e.g., Andrews et al., 2016; ALMA Partnership et al., 2015; Isella et al., 2016; Keppler et al., 2018; Sallum et al., 2015). Although still limited to a few objects, these observations inform about the processes responsible for the formation of planets.

Several circumstellar disks observed at a spatial resolution better than 50 au reveal ring-like features in the emission of small and large dust particles, as well as of the molecular gas (Andrews et al., 2011b, a; Boehler et al., 2018; Dong et al., 2017, 2018; Fedele et al., 2017, 2018; Isella et al., 2013, 2014; Tang et al., 2017; van der Marel et al., 2015, 2016, 2018b; Zhang et al., 2014, 2016). In particular, the homogeneous survey performed by the Disk Substructure at High Angular Resolution Project (DSHARP) has revealed that multiple-ring systems are ubiquitous among the most massive circumstellar disks (Andrews et al. 2018). The number of rings and their structure vary substantially from object to object, and in some cases, even within the same disk (Huang et al. 2018).

A myriad of theoretical models have been proposed to explain the formation of rings in circumstellar disks. These include the interaction between the disk and yet-unseen giant planets (e.g. Bryden et al., 1999; Zhu et al., 2014; Dong et al., 2018), sharp opacity variations at gas-solid phase transitions (Banzatti et al., 2015; Zhang et al., 2015; Okuzumi et al., 2016), dust accumulations at the edge of low viscosity regions (Flock et al., 2015; Miranda et al., 2017), and zonal flows via spontaneous concentration of net vertical flux (Bai & Stone, 2014; Béthune et al., 2017; Suriano et al., 2018). However, to date, planet-disk interaction models have been the most successful in explaining the observed structures (Jin et al., 2016; Liu et al., 2018; Dong et al., 2018). Taken at face value, these results tend to suggest the existence of a population of young gas giant planets orbiting at several tens of au from the central star which challenge current planet formation models. Furthermore, the link between planets and rings is supported by the direct detection of possible young planets in the PDS 70 and LkCa 15 systems (Keppler et al., 2018; Sallum et al., 2015), although the nature of some candidates is debated (Mendigutía et al., 2018; Thalmann et al., 2016).

In this paper, we present new ALMA observations of the 1.25 mm dust continuum and 12CO J=2-1 line emission of the circumstellar disk around the Herbig Ae star HD 163296. The observations achieve a spatial resolution of 4 au and deliver the sharpest images of this source obtained to date. At a distance of 101 pc (Gaia Collaboration et al., 2018), HD 163296 is in several respects the poster child of disks thought to be perturbed by planets. The star is surrounded by a huge (more than 1000 au in diameter) Keplerian disk whose mass has been estimated to range between 0.01 and 0.15 M⊙ (Isella et al., 2007; Muro-Arena et al., 2018; Tilling et al., 2012). Previous ALMA observations that resolved the gas and dust emission on spatial scales of 20 au, revealed the presence of three circular gaps in the disk density with radii of 45, 87, and 140 au, and, correspondingly, three dense rings centered at 68, 100, and 160 au (Isella et al., 2016; Zhang et al., 2016). The recent detection of deviations from Keplerian rotation at the location of the gaps and rings by Teague et al. (2018) confirms that these structures correspond to variations in the gas pressure (i.e., variations in the gas density and/or temperature), as opposed to radial changes in the dust opacity. The comparison between the observations and planet-disk interaction models indicates that the density gaps might have been carved by planets with masses between 0.5 and 1 MJ orbiting at 48, 86, and 131 au from the central star (Liu et al., 2018; van der Marel et al., 2018a). However, models that assume very low viscosity ($\alpha<10^{-4}$) suggest that the observed multiple ring structure may also result from the gravitational interaction with a single planet less massive than Saturn at 100 au orbital radius (Dong et al., 2018). Furthermore, the presence of an additional planet has been proposed based on the detection of local deviations from Keplerian velocity. Current models suggest that this planet might be located at 260 au from the star and has a mass of about 2 MJ (Pinte et al., 2018).

The HD 163296 system has also been the target of a number of optical and infrared high-contrast imaging campaigns aimed at characterizing the morphology of the disk in scattered light emission and detecting low mass companions (Benisty et al., 2010; Garufi et al., 2014; Grady et al., 2000; Muro-Arena et al., 2018; Wisniewski et al., 2008). Whereas these observations revealed rings similar to those observed at millimeter wavelengths, they did not detect any stellar or planetary mass object at orbital radii larger than 25 au and down to a mass sensitivity of a few MJ (Guidi et al., 2018). Thanks to the unprecedented angular resolution, our new ALMA observations provide detailed information on the morphology of the HD 163296 disk that help understanding the origin of the observed structures.

The structure of the paper is the following. The observations and data reduction are discussed in Section 2. The map of the dust continuum emission and its analysis are presented in Section 3, while CO data are discussed in Section 4. The implications of our results in the context of constraining the temperature of the disk and the origin of the observed structures are presented in Section 5. Finally, the main results of our investigation are summarized in Section 6.

2 Observations

The observations used in this study were obtained during the DSHARP ALMA Large Program (2016.1.00484.L) and are discussed in detail in Andrews et al. (2018). In brief, the HD 163296 disk was observed in Band 6 ($\lambda\sim 1.3$ mm) in September 2017 in configuration C40-8, which delivered baselines $B$ between 41 m and 5.8 km, and a theoretical angular resolution $\lambda/B_{max}\sim 0.04\arcsec$. Archival data from project 2013.1.00366.S (Flaherty et al., 2015) and project 2013.1.00601.S (Isella et al., 2016) were used to improve sensitivity and uv-coverage on shorter baselines.

An initial calibration of each dataset was produced by the ALMA pipeline. A visual inspection of the data obtained from project 2013.1.00366.S revealed a number of visibilities characterized by very noisy amplitudes that required manual flagging. Short baseline observations were independently self-calibrated and re-centered to account for shifts in the position of the disk due to the proper motion of the source. We compared the flux calibration of each track and adjusted the flux scale so that the source visibility amplitude measured on overlapping baselines agrees within 5%. Finally, we combined all the tracks and performed phase and amplitude self-calibration on both short and long baseline data. The self-calibration procedure resulted in an improvement of 43% in the peak signal-to-noise ratio in the dust continuum image. A complete description of the DSHARP data calibration procedure is presented in Andrews et al. (2018), while the CASA script used to calibrate and image HD 163296 data, including the manual flagging and flux rescaling, is available online at https://almascience.org/alma-data/lp/DSHARP. The 1.25 mm dust continuum emission was imaged using the CASA task tclean and a robust parameter of -0.5 resulting in a synthesized beam FWHM of 0.038″$\times$0.048″, corresponding to a spatial resolution of 3.8$\times$4.8 au at the distance of the source. The rms noise is 0.023 mJy beam*-1*, and the peak signal to noise ratio is 185.

The complex gain solutions of the self-calibration of the continuum emission were then applied to the 12CO data. The line was observed at the same angular resolution of the continuum ($\sim$0.04″), but we imaged it using Briggs robust=0.5 (no uv-tapering) to achieve higher signal-to-noise ratio (see Andrews et al. 2018 for more details about CO imaging). The FWHM of the synthesized beam of the presented CO maps is 0.104″$\times$0.095″, which is about half the beam of the CO observations published in Isella et al. (2016), and about 7 times smaller than the resolution of the ALMA science verification data (Rosenfeld et al., 2013; de Gregorio-Monsalvo et al., 2013). Channels are spaced in velocity by 0.32 km s*-1*, but, due to Hanning smoothing, the velocity resolution is 0.64 km s*-1*. The velocity grid for HD 163296 slightly differs from the fiducial channel spacing of 0.35 km s*-1* adopted for the other DSHARP sources. The rms noise per channel is 0.84 mJy beam*-1*.

3 Dust continuum emission

The map of the 1.25 mm continuum emission (Figure 1) features two bright elliptical rings previously reported by Isella et al. (2016), as well as three new morphological features: an arc of emission inside the first ring (inset a), an inner dark gap and bright ring at about 10 au and 15 au, respectively, and a central azimuthal asymmetry (inset b). In this section we discuss the morphology of the dust rings, constrain the radial profile of the continuum emission, and investigate the presence of asymmetric structures in the distribution of the solid material.

3.1 Morphology of the dust rings

As a first step toward characterizing the structure of the dust rings, we fit the crests and troughs of the continuusm emission with circles defined by the position of the center ($\Delta x_{0}$, $\Delta y_{0}$), the radius $r$, the inclination $i$ at which they are observed, and the position angle (PA) of their projected major axis. The fit is performed as discussed in Huang et al. (2018). In brief, we locate the points (x,y) corresponding to the radial maxima (or minima) of each ring (or gap), and use the python emcee package to calculate the ellipses that best reproduce this set of points. The best fit parameters are listed in Table 1. Bright rings and dark gaps are labeled with the suffixes B and D, respectively, while the number appearing in the feature name corresponds to their radius in au. The uncertainties on the best fit values correspond to the 16th and 84th percentile of the marginalized posterior probability distribution.

The outermost features B67, D86, B100, D145, and B155, are well described by concentric circles with an average inclination of (46.7$\pm$0.1)°  and average position angle of (133.3$\pm$0.1)°. Taken at face value, the circular model fitting indicates that D10 and D45 have lower inclinations. We believe that the lower inclination of D45 is due to the presence of the dust crescent centered at 55 au, which partially fills the gap making it appear more circular (i.e. less inclined), while the lower inclination of D10 is compatible with the effect of beam smearing which is discussed in more details below.

Figure 2 shows the deprojected map and the azimuthally averaged radial profile of the continuum emission obtained using the values for the disk inclination and position angle derived from the outermost rings. In the deprojected map, the outer rings appear as vertical stripes, confirming that these structures have indeed similar inclination and position angle, and that they are intrinsically circular in shape. The intensity ratio between B14 and D10 varies with the azimuthal angle reaching a maximum of about 1.2 along the disk major axis ($\theta=\pm 90$°) and a minimum of about 1 along the disk minor axis ($\theta=0$°). This variation can be explained with the fact that the effective spatial resolution for an inclined disk is higher along the disk major axis compared to the disk minor axis. In practice, B14 and D10 are spatially resolved only along the disk major axis, while they blend together along the disk minor axis.

Azimuthal variations in the contrast between bright rings and dark gaps are also visible for the outer pairs. For example, along the westward semi-major axis ($\theta=-90$°), the intensity reaches a minimum consistent with the noise of the observations at the center of D45, and a maximum of 0.91 mJy beam*-1* at the center of the adjacent bright ring B67. This corresponds to an intensity ratio larger than 37. In the opposite direction, i.e. along the eastward disk semi-major axis ($\theta=90$°), the intensity ratio B67/D45 is 13. As comparison, the intensity ratio along the southward ($\theta=180$°) and northward ($\theta=0$°) direction of the disk minor axis is about 7 and 5, respectively (see Table 2 for a summary of the intensity ratio among rings pairs). The difference between intensity ratios measured along the disk major and minor axes is likely due to beam smearing, while the difference between the intensity ratios measured along the disk major axis suggests the presence of azimuthal asymmetries in the dust emission.

The center and right panels of Figure 2 show the azimuthally averaged profile of the deproject continuum emission expressed in units of brightness temperature $T^{pl}_{b}$ calculated using the full Planck equation as

[TABLE]

where the $F_{\nu}$ is the flux density integrated over a beam defined by the solid angle $\Omega=\pi\theta_{min}\theta_{max}/(4\ln 2)$, where $\theta_{min}$ and $\theta_{max}$ are the minimum and maximum FWHM of the synthesized clean beam. The brightness temperature varies between about 50 K (corresponding to an intensity of 3.8 mJy beam*-1*) in the innermost disk regions down to about 3 K (or a flux density of 0.0025 mJy beam*-1*) at about 180 au. The pairs B67/D45 and B100/D86 have brightness temperature ratios of 4.3 and 2.4, respectively. Note that the difference between brightness temperature and intensity ratios is due to the non-linearity of Equation 1.

3.2 Width of the continuum rings

Dullemond et al. (2018) finds that many of the continuum rings revealed by DSHARP observations, including the rings B67 and B100 of the HD 163296 disk, have a radial width narrower than the estimated pressure scale height of the gaseous disk, and conclude that this is strong evidence that these rings have formed by radial migration of dust particles toward local maxima of the gas pressure. Dullemond et al. (2018) defines the width of a ring as the dispersion of a Gaussian function that reproduces the deprojected and azimuthally averaged radial profile of the dust emission across the ring itself. This modeling approach is very fast, but it does not account for the observational noise, nor for the discrete sampling of the uv plane performed by an interferometer, or for the fact that the synthesized beam is not circular.

Here, we compare the widths of the HD 163296 rings measured by Dullemond et al. (2018) with those estimated by a more accurate, but much slower, modeling of the continuum visibilities in the uv domain. We start by assuming that the dust emission is axisymmetric and expressed by a sum of $n$ Gaussian functions

[TABLE]

We generate synthetic 2D images of the continuum emission that are inclined and rotated to simulate all possible disk viewing angles. The center of emission is allowed to vary with respect to the center of the image (i.e. the phase center). Synthetic visibilities ($V^{mod}$) are then calculated by taking the Fourier transform of the image at the spatial frequencies of the observed visibilities ($V^{obs}$). The goodness of the model fit is evaluated through a usual $\chi^{2}$ test, where $\chi^{2}=\sum{w(V^{obs}-V^{mod})^{2}}$, and $w$ is the weight of each visibility measurement provided by the ALMA pipeline. The calculation of $\chi^{2}$ was performed using the cpu-based version of python package galario (Tazzari et al., 2017). Finally, the $\chi^{2}$ is used to sample the posterior likelihood using the python package emcee (Foreman-Mackey et al., 2013).

An initial exploration of the parameter space indicates that a minimum of 5 Gaussian components are required to reproduce the observed continuum emission. In total, the model therefore has 19 free parameters: the coordinates ($x_{0}$, $y_{0}$) of the center of the emission, the disk inclination $i$ and position angle PA, and the radius $r$, peak intensity $I$, and width $w_{d}$ of five Gaussian rings. We run emcee using 200 walkers initialized around the parameter values obtained by a Gaussian fitting of the radial intensity profile as in Dullemond et al. (2018). Each walker is evolved for an initial burn-in run of $10^{4}$ steps, after which they all converge toward a stable set of parameters. Following Goodman & Weare (2010), we sample the posterior distribution by letting emcee run until convergence, which is established based on the mean autocorrelation time of all the model parameters. We find that convergence is achieved after about $4\times 10^{4}$ steps. The optimal model parameters are then calculated as the median of the marginalized probability distribution and are listed in Table 3, while the related uncertainties correspond to the 0.3th and 99.7th percentile of the marginalized posterior probability distribution. This would be equivalent to a 3$\sigma$ uncertainty for a Gaussian posterior distribution. An image constructed from the residual visibilities is shown in Figure 3.

For B67, we find a width of $\sim 6.6$ au, that is very close to the 6.8 au derived by Dullemond et al. (2018). Conversely, the width of B100 is about 25% larger than that measured by fitting the intensity profile. As noted in Figure 3 of Dullemond et al. (2018), the intensity radial profile of B100 has a peaked central region and bright wings that are not well reproduced by a Gaussian function. The value of the width in Dullemond et al. (2018) results from fitting only the innermost part of the ring between 94 au and 104 au, while our value resulted from modeling the global profile including the bright wings of the ring. The fact that the two modeling procedures achieve the same result for rings that are intrinsically Gaussian (i.e., B67) indicates that the results of Dullemond et al. (2018) are independent of the exact shape of the beam, or, more generally, from non linear effects caused by the discrete uv-sampling and image deconvolution. In practice, the uv-sampling of the observation is sufficiently uniform and the signal-to-noise ratio of the continuum map is sufficiently high to allow to perform the data analysis in the image domain without loss of accuracy. This is a great advantage since modeling the dust rings in the image domain takes only a few seconds, compared to the several days, or weeks, required to perform the model fitting in the Fourier domain.

Our analysis indicates that the ring B14 has a radius of about 15.5 au, slightly larger than the value calculated by fitting the crest of the emission as discussed in Section 3.1 (see also Table 1). The best fit model also includes a small ring of 4 au in radius, therefore suggesting the presence of a inner circular cavity like those observed in transitional disks. However, as discussed in the next section, our azimuthally symmetric best fit model poorly matches the dust emission from the innermost disk regions which instead seems to indicate the presence of a dust crescent.

Finally, the last component of the model is a ring with a radius of about 32 au that is required to reproduce the kink in the radial profile of the emission observed at about 30 au from the center. Similar kinks in the radial profiled of the dust continuum emission are observed in other DSHARP sources and are discussed in Huang et al. (2018).

3.3 Asymmetries of the continuum emission

The map obtained by imaging the residual visibilities obtained by subtracting the best fit Gaussian model from the observations (Figure 3) reveals several features with absolute intensities above 5 times the rms noise and peak intensities as high as 30 times the rms noise. The two most prominent structures have the shape of crescents and were already introduced at the beginning of Section 3 (see panels b and c of Figure 1). The innermost crescent is located at about 4 au from the center and it overlaps with the innermost ring component of the Gaussian model discussed above. Given its azimuthally asymmetric structure, and following the nomenclature adopted for the rings and gaps, we label this feature as C4 (C meaning “crescent”). The peak of C4 has a position angle of -44°  and an intensity is 0.55 mJy beam*-1*, or about 20% of the continuum intensity measured at the same position. C4 extends by about 180°  in azimuth and by about one resolution element in radius. C4 resembles the crescents observed in other Herbig Ae disks that have been attributed to the presence of anticyclonic vortices (e.g. Casassus et al., 2013; Isella et al., 2013; van der Marel et al., 2013; Boehler et al., 2018) or to an optically thick inner disk warped with respect to the outer disk (e.g. Marino et al., 2015; Benisty et al., 2017). However, higher angular resolution observations are required to investigate the morphology of HD 163296 innermost disk regions in greater details.

The second crescent, labelled as C55, is centered at an orbital radius of about 55 au and position angle of 99°. The peak intensity is 0.64 mJy beam*-1*, corresponding to a signal-to-noise ratio of 28. Modeling C55 intensity as a 2D Gaussian function in the polar plane $I(r,\theta)\propto e^{-[(r-r_{c})^{2}/2\sigma_{r}^{2}+(\theta-\theta_{c})^{2}/2\sigma_{\theta}^{2}]}$, we find $\sigma_{\theta}=17$ au and $\sigma_{r}=3$ au. Accounting for beam convolution (Equation 3), the intrinsic radial and azimuthal widths of C55 are $\sim 2.2$ au and $\sim 16.9$ au, respectively.

In addition to C4 and C55, the image of the residuals shows fainter asymmetries along the B67 and B100 dust rings. If the dust emission is optically thin, intensity variations along the rings might probe the clumpiness of the dust distribution, while if the emission is optically thick, they inform about variations in the dust temperature. The intensity profiles along B67 and B100 are shown in the bottom right panel of Figure 3. In the case of B67, the residual intensity varies between $\pm$0.15 mJy beam*-1*, to be compared with an azimuthally averaged intensity value measured at the same radius of 1.0 mJy beam*-1*. The azimuthal intensity profile of B67 shows small amplitude variations characterized by an angular scale of about 4°, and larger variations with angular scales of 30-50°. Whereas the small scale variations have the same size of the synthesized beam and are most likely caused by the noise of the observations, the larger scale structures might trace real asymmetries in the dust emission and, since the emission from B67 is estimated to be optically thin (see Section 5.1), in the dust density.

4 12CO J=2-1 line emission

Figure 4 shows continuum-subtracted channels maps of the 12CO (2-1) emission (hereafter referred simply as CO emission). Below we describe the morphology and kinematics of the CO emission, and use the observations to constrain disk properties such as the disk temperature and the optical depth of the dust continuum emission.

4.1 Morphology of the CO emission

The CO emission extends far beyond the outer edge of millimeter-wave continuum emission up to an angular distance of 5.6″  or about 560 au, from the central star. As discussed in Isella et al. (2016), the discrepancy between the radial extent of the millimeter wave dust and CO emission indicates a sharp drop in either the dust opacity or dust column density beyond about 200 au.

The CO emission observed at intermediate velocities (e.g., $v_{lsr}=3.40-5.00$ km s*-1* and $v_{lsr}=6.28-8.20$ km s*-1*) shows a characteristic “dragonfly” structure tracing both the front and rear warm CO emitting layers of the disk (Figure 5). The emission coming from the rear of the disk is fainter than that from the front due to a combination of lower temperature of the emitting gas and absorption from intervening material. As discussed in Rosenfeld et al. (2013), the fact that the front and rear sides of the disk appear as two distinct emitting regions is explained by the large optical depth of the CO line (whose intensity therefore mostly depend on the gas temperature) and a vertical gradient of the disk temperature characterized by warm surfaces and a cold disk midplane. Furthermore, the separation between the front and rear CO emission measures the vertical geometry of the CO emitting surfaces.

We compare the observed geometry of the CO emission with that of a parametric model where the CO line originates from a geometrically thin layers characterized by a distance from the midplane $z_{co}=\pm z_{co,0}(r/r_{0})^{q}$. As demonstration of the technique, we show in Figure 5 parametric isovelocity curves corresponding to the front side of the disk calculated as presented in the Appendix A. The geometry of the isovelocity curves depends on $z_{co,0}$, $q$, as well as on the stellar mass and disk inclination. Assuming the disk inclination derived from the continuum ($i=46.7$°) and a stellar mass of 2.0 M⋆, we obtain a good match between model and observations for $q=0.5$ and $z_{co}(100\mathrm{au})=30$ au.

Our observations reveal that the CO emission coming from the back side of the disk is dimmer at the location of the bright dust rings B67 and B100 compared to the emission observed at the location of the dark gaps D45 and D87 (Figure 6). This dimming, which is likely due to absorption of the line emission as it crosses the dusty rings, provides a tool to constrain the optical depth of continuum emission independently from the assumptions on the dust temperature. The CO emission coming from the back side of the disk can be expressed as $I_{co}^{obs}=I_{co}e^{-\tau_{d}^{ext}}$, where $\tau_{d}^{ext}$ is the line of sight dust extinction optical depth equal to the sum of the absorption and scattering optical depths ($\tau^{ext}_{d}=\tau^{abs}_{d}+\tau^{sca}_{d}$), while $I_{co}$ is the initial unattenuated CO emission. Defining the parameter $\epsilon$ as $\epsilon\equiv\tau_{d}^{abs}/\tau_{d}^{ext}$, the absorption optical depth of the dust continuum can be written as

[TABLE]

The observations of the CO emission directly measure $I_{co}^{obs}$ along the bright continuum rings, while the unattenuated emission $I_{co}$ can be estimated through the procedure illustrated in Figure 6. In practice, in the channels where CO absorption is observed, we measure the intensity of the CO emission along the isovelocity curves corresponding to the back side emission. We then use the values measured within the dust gaps, where, based on the results of the previous section, $\tau_{d}\ll 1$, to estimate $I_{co}$ at the position of the bright rings. We measure the ratio $I_{co}/I_{co}^{obs}$ for B67 and B100 along 18 different lines of sight (i.e., 18 different azimuthal angles along the dust rings) and find mean values of 1.9$\pm$0.1 and 2.1$\pm$0.1. We also perform this measurement for B168 finding $I_{co}/I_{co}^{obs}=1.1\pm 0.2$. The corresponding dust extinction optical depths are $\tau^{ext}_{d}(\mathrm{B}67)=0.64\pm 0.05$, $\tau^{ext}_{d}(\mathrm{B}100)=0.74\pm 0.05$, $\tau^{ext}_{d}(\mathrm{B}168)=0.1\pm 0.2$. It is worth noting that the ratio $I_{co}/I_{co}^{obs}$ could vary with the azimuthal angle if $\tau_{d}^{ext}$ is not constant along the dust rings. However, the sensitivity of our observations does not allow to perform such an analysis, and we must instead rely on the azimuthally averaged value of the CO intensity ratio to investigate the dust extinction properties.

As from Equation 3, the absorpion optical depth of the dust rings depends on the dust scattering properties through the parameter $\epsilon$. For dust grains much smaller than the wavelength of the observations (in our case, for particle sizes $a\ll\lambda/2\pi\sim 0.2$ mm), scattering is negligible compared to absorption, and $\epsilon\simeq 1$. For grain sizes comparable to the wavelength of the observations ($a\sim\lambda/2\pi\sim 0.2$ mm), scattering and absorption properties of spherical grains follow Mie theory, and $\epsilon$ can vary between about 0.1 and 1, depending on the dust composition and internal structure. Finally, if $a\gg 0.2$ mm, scattering and absorption are equal, and $\epsilon=0.5$ (see, e.g., Bohren & Huffman, 1998). Without any prior information about the size and composition of dust grains, the observed absorption of the CO emission constrains the absorption optical depth of B67 and B100 between about 0.07 and 0.7, while the lack of CO absorption at the location of B168 sets an upper limit to the absorption optical depth of this ring at about 0.3. A caveat of this approach is that the emission from the front side of the disk could partially contaminate the emission from the back side, particularly close to the star where the separation between the two surfaces is small. Also, the ratio $I_{co}/I_{co}^{obs}$ could be affected by beam smearing, implying that our measurements provide only a lower limit for the true contrast between the initial and absorbed CO emission, and consequently, the corresponding dust optical depth. However, the angular resolution of the CO observations is sufficiently high to believe that this effect is small.

If the absorption optical depth of the dust rings is measured from CO absorption, the corresponding continuum intensity $I_{\nu}$ can be used to estimate the temperature $T_{d}$ of the emitting dust. If dust scattering is negligible ($\tau_{d}^{ext}=\tau_{d}^{abs}$ and $\epsilon=1$), the conversion between flux density and dust temperature is

[TABLE]

or, using Equation 3, as

[TABLE]

From the measured continuum intensities we obtain dust temperatures of 24 K and 15 K for B67 and B100, respectively, while the upper limits for the optical depth of B168 translates in a minimum temperature of about 6 K.

If dust scattering is not negligible ($\tau_{d}^{ext}>\tau_{d}^{abs}$ and $\epsilon<1$), the dust temperature can be estimated using the formalism presented in Birnstiel et al. (2018). In this case, the intensity is written as $I_{\nu,d}\sim S(T_{d},\tau_{\nu,d}^{ext},\epsilon)(1-e^{-\tau_{\nu,d}^{ext}})$, where the source function $S_{\nu}$ depends on the dust temperature, and on the extinction and scattering optical depth along the line of sight. Adopting $\epsilon$ in the range between 0.1 and 1, we obtain dust temperatures between 130 K and 24 K for B67, and between 70 K and 15 K for B100. The fact that the dust temperature depends on the dust scattering properties hampers the application of this technique as an independent measure of disk temperature. However, if the dust temperature is measured through other means, then the dust scattering properties can by constrained (see Section 5).

4.2 Gas temperature

At typical densities and temperatures of protoplanetary disks, the peak of the 12CO J=2-1 line emission is optically thick and the CO rotational levels are in local thermodynamic equilibrium (see, e.g., Weaver et al., 2018). The intensity at the peak of the line therefore measures the temperature of the emitting gas. Following Weaver et al. (2018), we generate a CO temperature map from non-continuum subtracted line emission to avoid underestimating the line intensity along line of sights where the continuum is optically thick (Figure 7). The CO temperature reaches values as high as 120 K in the innermost disk regions and as low as 20 K at the disk outer edge. The peak intensity map is characterized by an evident asymmetry with respect to the apparent disk major axis: the south-west side of the disk is hotter than the north-east side. This is explained in first approximation by the fact that the CO emitting layer is a conic surface observed from a direction inclined with respect to the axis of the cone. The orientation of the asymmetry implies that the northern side of the disk is the closest to the observer (see, also, Rosenfeld et al., 2013).

Assuming, as in the previous section, that the peak of the line is emitted from a geometrically thin layer at a distance from the disk midplane given by $z_{co}=z_{{co},0}(r/r_{0})^{q}$, it is possible to deproject the observed CO temperature and estimate its profile as a function of the cylindrical radius $r$. Using the values for $z_{co}$ and $q$ presented above, we derive the CO temperature profile shown in the right panel of Figure 7. Between 30 and 500 au, the radial profile of the CO temperature scales almost linearly with the radius and follows the relation $T_{co}(\mathrm{K})\sim 87-0.14(r/\mathrm{au})$. Within 30 au, the CO brightness temperature drops. This is consistent with the expected drop in brightness temperature as the emitting area becomes smaller than the angular resolution of the observations (beam dilution).

4.3 Kinematics of the CO emission

The overall kinematics of the CO gas is well described by Keplerian rotation around a star with a mass of about 2.0 M⊙ (Figure 8). Both the integrated velocity map and the position-velocity diagram are in first approximation symmetric with respect to the apparent disk minor axis. The systemic velocity inferred from CO observations is 5.8 km s*-1* and is in agreement with previous results. The velocity-position diagram also reveal diffuse emission at a velocities between about 12.5 and 14 km s*-1* that might trace a redshifted stellar wind Klaassen et al. (2013) or perhaps emission from the left-over of the HD 163296 parent cloud.

The CO map registered at a velocity of 6.92 km s*-1* (Figure 4 and Figure 5) is characterized by a kink at $\delta$RA$\sim$-1″,$\delta$Dec$\sim 1.5$″. This feature was previously reported by Pinte et al. (2018) and has been attributed to perturbations in the gas kinematics caused by a planet with a mass of 2 MJ orbiting at 260 au from the star. Unfortunately, despite the better angular resolution, the observations obtained using the extended ALMA configuration achieve a spectral resolution about 10 times worse than that of those presented by Pinte et al. (2018), and do not allow us to investigate the origin of this feature. Similarly, the coarse velocity resolution of our data might not allow to investigate the small deviations from keplerian rotation observed within the dust gaps and reported by Teague et al. (2018).

5 Discussion

The two main results obtained so far are that (i) the combination of continuum and CO observations directly inform about the temperature and geometry of the HD 163296 disk and (ii) the analysis of the continuum emission constrains radial profile and level of asymmetry of the observed dust rings. In this section, we firstly compare the disk temperature inferred from the observations with current models for HD 163296 and, secondly, we discuss the possible origin of the morphology of the gas and dust emission.

5.1 Disk temperature and dust scattering

To date, the HD 163296 disk has been the subject of several studies aimed at constraining the density and thermal structure of the disk through forward modeling of the spatially resolved millimeter-wave molecular line and continuum emission (e.g., Isella et al., 2007, 2016; Rosenfeld et al., 2013; Flaherty et al., 2015; Tilling et al., 2012). This approach relies on a large number of assumptions including the choice of the parametric forms for the gas and dust surface density and temperature (when the latter is not calculated using radiative transfer models), the dust and molecular abundances, the dust opacities, the stellar properties (mass, luminosity, distance), and the gas kinematics (e.g. Keplerian motion, turbulent velocity). Because of the large computational time required to generate synthetic images at high spatial and spectral resolution, the comparison between models and observations is generally performed by either manually adjusting the model parameters until obtaining a satisfactory result (e.g. Isella et al., 2016), or by using automatic algorithms to search for best fit models among a subset of model parameters. Not surprisingly, forward modeling performed by different investigators lead to different results.

Figure 9 compares the temperature inferred from the 1.25 mm dust continuum intensity to the models discussed in Isella et al. (2016), Rosenfeld et al. (2013), and Flaherty et al. (2015). We choose these three studies because they assume a similar parametric form for the disk temperature and are based on high spectral and angular resolution observations of the 12CO line emission. If scattering from dust grains is negligible ($\tau^{ext}_{d}=\tau^{abs}_{d}$, $\epsilon=1$), the brightness temperature $T_{b}^{pl}$, the physical dust temperature $T_{d}$, and the optical depth $\tau^{abs}_{d}$ are linked by the ralation $B_{\nu}(T_{b}^{pl})=B_{\nu}(T_{d})(1-e^{-\tau_{d}^{abs}})$. Within 30 au, the model temperatures are bracketed between Rosenfeld et al. (2013) and Isella et al. (2016), and would imply optical depths of about 1.3 and 0.7, respectively. Beyond 30 au, $T_{b}^{pl}$ drops well below the model predictions, suggesting that the emission becomes more optically thin. In particular, the model temperatures give $\tau_{d}^{abs}<0.02$ and $\tau_{d}^{abs}\sim 0.04-0.05$ at the center of the dark gaps D45 and D86, respectively, and $\tau_{d}^{abs}\sim 0.5-0.7$ and $\tau_{d}^{abs}\sim 0.3-0.5$ at the bright rings B67 and B100.

These measurements rely on the assumption that scattering is negligible, but, as discussed in Section 4.1, this is likely not the case if dust grains have sizes comparable or larger than the wavelength of the observations. The inclusion of scattering significantly complicates the interpretation of the observed continuum emission and a simple relation between $T_{b}^{pl}$, $T_{d}$, and the absorption and scattering optical depths, $\tau^{abs}_{d}$ and $\tau^{sca}_{d}$ does not exist. Fortunately, however, it is possible to derive an approximated analytic solution to constrain $T_{d}$, $\tau^{abs}_{d}$ or $\tau^{sca}_{d}$ if the other two quantities are known (see Section 5 of Birnstiel et al. 2018).

In the case of HD 163296, the direct measurement of the extinction opacity at the bright rings B67 and B100 discussed in Section 4.1 provide a unique tool to constrain the scattering properties of the circumstellar dust. Using the formalism discussed in Birnstiel et al. (2018), we calculate the physical dust temperature required to reproduce the dust continuum emission measured at B67 (0.77 mJy beam*-1*) and B100 (0.45 mJy beam*-1*) for values of the scattering parameter $\epsilon$ varying from 0.1 ($\tau_{d}^{sca}=9\tau_{d}^{abs}$) to 1 (no scattering). The comparison between these temperatures, which are indicated with horizontal segments in the left panel of Figure 9, and those predicted by Flaherty et al. (2015), Isella et al. (2016), Rosenfeld et al. (2013) constrain $\epsilon$, and therefore the dust scattering properties. For B67, we find that $\epsilon\gtrsim 0.6$ leads to dust temperature consistent with the models. We note that the midplane temperature of Flaherty et al. (2015) is slightly below the minimum temperature consistent with the observations, which corresponds to $\epsilon=1$. Interestingly, for B100, the consistency between models and observations requires $\epsilon$ between 0.4 and 0.6, while $\epsilon=1$ would imply a physical dust temperature of only 15 K. The difference in scattering properties between B67 and B100 is intriguing and, taken at face value, might suggest a variation in the dust properties between B67 and B100. However, the difference might also result from the fact that, due to beam smearing, the measured value of $\tau_{d}^{ext}(\mathrm{B}67)$ might be a lower estimate of the extinction opacity, and consequently, of the scattering opacity. Future ALMA observations expressly designed to image the molecular gas emission from the HD 163296 at both high angular resolution and sensitivity are required to place more stringent constraints on the dust extinction optical depth.

We conclude this section by comparing the peak brightness temperature of the CO emission measured in Section 4.2 to the temperature of the CO emitting layer predicted by Flaherty et al. (2015), Isella et al. (2016), and Rosenfeld et al. (2013) models (see right panel of Figure 9). In Section 4.1 we found that the CO emitting layer corresponds to the surface defined by $z(r)=30\mathrm{au}(r/100\mathrm{au})^{0.5}$. Within about 100 au, the measured gas temperature is below the model predictions. This is likely due to the effect of beam dilution as discussed in Weaver et al. (2018). Between 100 and 400 au, the temperature of the CO emitting layer predicted by Rosenfeld et al. (2013) and Isella et al. (2016) is in good agreement with the observations, while the temperature predicted by Flaherty et al. (2015) is about 30% lower. Also, while the temperature of the CO emitting layer predicted by Rosenfeld et al. (2013) and Isella et al. (2016) is similar to that of the stellar irradiated disk surface (dashed line), as expected in the case of a very optically thick line, Flaherty et al. (2015) temperature is closer to that achieved in the disk midplane (dotted lines).

5.2 On the origin of the dust rings

In Section 3, we found that the radial profile of the continuum emission from the HD 163296 disk can be represented by a sum of azimuthally symmetric Gaussian rings. However, significant asymmetric structures appear once the symmetric component of the emission is removed. In this section, we discuss these results in the framework of the planet-disk interaction mechanism, while we refer to Huang et al. (2018) for a discussion of other processes capable of generating dust rings. We also refer the reader to Dullemond et al. (2018) for a discussion of the observed rings properties in the context of dust trap models.

Massive planets are expected to open gaps in the surface density of protoplanetary disks along their orbits through the strong mutual gravitational interaction with the circumstellar material (see, e.g., Bryden et al., 1999). While the planet minimum mass for this to happen varies as a function of the disk physical properties, most of the simulations agree that for typical disk conditions planets more massive than 20-40 Earth masses should produce gaps observable at the spatial scales probed by ALMA (e.g. Dong et al., 2012; Isella et al., 2018). In the limit of nearly inviscid disks, simulations have also shown that even planets with masses as low as few Earth masses might generate visible rings (Dong et al., 2018). Given the ubiquity of planets in the Universe, it seems therefore reasonable to associate the ring structures observed by ALMA with the presence of planetary systems in the act of forming.

In the case of HD 163296, the planet-disk interaction hypothesis is supported by the detection of gas depletion and deviation from Keplerian rotation within the dark gaps D87 and D143 (Isella et al., 2016; Teague et al., 2018), which imply that the emission gaps correspond to gaps in the gas density. The comparison with hydrodynamic simulations indicates that D45, D87 and D143 might results from the dust clearing operated by planets with masses between 0.5 and 2 MJ (Liu et al., 2018; Teague et al., 2018). However, the ring morphology of the HD 163296 disk might also be consistent with gravitational perturbations by a single planet with a mass as low as 65 ME orbiting at about 100 au from the star if the disk viscosity is very low ($\alpha<10^{-4}$, Dong et al., 2018). Overall, these results are consistent with the planet masses discussed in Zhang et al. (2018), which have been estimated from the azimuthally averaged profile of the continuum emission shown in Figure 2. In addition, Zhang et al. (2018) find that the newly discovered gap B10 might be generated by a planet with a mass between 0.2 and 1.5 MJ.

Here we want to point out that the detection of asymmetries in the dust emission discussed in Section 3.3 corroborates the hypothesis that the HD 163296 disk is shaped by the interaction with yet unseen planets. As shown in Figure 6 and 7 of Zhang et al. (2018), features similar to the crescent observed inside the D45 gap (inset b of Figure 1) are naturally produced by the disk-planet interaction whenever the disk viscosity is low ($\alpha<10^{-3}$) and the planet mass is above about 0.1 MJ. The crescents observed in numerical simulations have two different origins: they either trace dust particles trapped in anticyclonic vortices that form at the edge of the gap opened by the planet (Li et al., 2001), or probe material trapped at the Lagrangian points along the planet orbit. Without any pretense of accuracy, we compare in Figure 10 the map of the 1.25 mm dust continuum emission recorded toward HD 163296 with a synthetic image of a disk model characterized by the presence of a 0.15 MJ planet orbiting at 54 au from the central star. The model is part of the suite of models presented in Zhang et al. (2018) and was generated assuming a viscosity parameter $\alpha=10^{-4}$ and a disk pressure scale height $h=0.05r$. This comparison is purely qualitative but nevertheless shows the resemblance between the corotation features predicted by models and the circular arc observed in HD 163296.

6 Conclusion

The unprecedented resolution of the DSHARP/ALMA data of the HD 163296 circumstellar disk enabled us to dive into a detailed characterization of the morphology of the dust continuum emission and of some of the physical properties of the circumstellar material. The comparison between images of the 1.25 mm continuum and 12CO line emission and parametric disk models led us to the following conclusions:

• •

The new ALMA observations confirm the presence of two bright rings in the millimeter-wave dust continuum emission centered at 67 and 101 au, and reveal an additional ring with a radius of 15 au. The radial profile of the dust emission across the rings is well described by a Gaussian profiles with widths of 8.7, 6.6, and 5.8 au, respectively.

• •

The dust continuum emission is characterized by several asymmetric structures. The most prominent consist in crescents centered at about 4 and 55 au. We also find asymmetric structures along the ring B67 characterized by amplitude variations of $\pm 15\%$ relative to the ring mean intensity. We argue that this asymmetries trace local variations in the dust densities and support the hypothesis that the HD 163296 disk is shaped by the gravitational interaction with yet unseen planets.

• •

The observations of the 12CO line emission probe the temperature and geometrical structure of the emitting gas. We find that the peak of the line arises from a disk layer at a vertical distance from the disk midplane given by $z_{co}=30\mathrm{au}(r/100\mathrm{au})^{0.5}$. Between 30 au and 500 au, the temperature of the CO emitting layer scales linearly from the orbital radius and follows the relation $T_{co}(\mathrm{K})\sim 87-0.14(r/\mathrm{au})$.

• •

The CO maps show the line emission coming from both the front and rear side of the disk. The latter is attenuated as it passes through the dust rings, providing a tool to directly measure the extinction optical depth $\tau_{d}^{ext}$ of the dust. We measure $\tau_{d}^{ext}\sim 0.7$ at the position of the B67 and B100 rings, and $\tau_{d}^{ext}<0.3$ for B168. We point out that if the absorption optical depth of the dust rings $\tau_{d}^{abs}$ can be measured through other means, a direct measurement of $\tau_{d}^{ext}$ allows to constrain the dust scattering properties.

• •

Finally, we compare the CO temperature inferred from our ALMA observations to theoretical predictions based on previous forward modeling of millimeter wave observations of HD 163296. Overall, we find that models predict temperatures within 30% of the measured values. By adopting the predicted temperatures for the disk midplane, we find that $\tau_{d}^{abs}$ is between about 0.5 and 0.6 at B67, and between 0.3 and 0.5 at B100. This implies $\tau_{d}^{sca}$ between 0.1 and 0.2 at B67, and between 0.2 and 0.4 at B100.

We thank the anonymous referee for the comments that helped improving the manuscript. A.I. acknowledges support from the National Aeronautics and Space Administration under grant No. NNX15AB06G issued through the Origins of Solar Systems program, and from the National Science Foundation under grant No. AST-1715719. S.A. and J.H. acknowledge support from the National Aeronautics and Space Administration under grant No. 17-XRP17$\_$2-0012 issued through the Exoplanets Research Program. J.H. acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1144152. L.P. acknowledges support from CONICYT project Basal AFB-170002 and from FCFM/U. de Chile Fondo de Instalación Académica. C.P.D. acknowledges support by the German Science Foundation (DFG) Research Unit FOR 2634, grants DU 414/22-1 and DU 414/23-1. V.V.G. and J.C acknowledge support from the National Aeronautics and Space Administration under grant No. 15XRP15$\_$20140 issued through the Exoplanets Research Program. T.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 714769. M.B. acknowledges funding from ANR of France under contract number ANR-16-CE31-0013 (Planet Forming disks). Z.Z. and S.Z. acknowledge support from the National Aeronautics and Space Administration through the Astrophysics Theory Program with Grant No. NNX17AK40G and the Sloan Research Fellowship. L.R. acknowledges support from the ngVLA Community Studies program, coordinated by the National Radio Astronomy Observatory, which is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. M.B. acknowledges funding from ANR of France under contract number ANR-16-CE31-0013 (Planet Forming disks).

Appendix A: calculation of isovelocity lines

The isovelocity lines shown in Figure 5 are calculated as follow. Firstly, we assume that CO molecules rotate around the star at the Keplerian velocity $v_{k}(R)=\sqrt{GM_{\star}/R}$, where $R$ is the spherical radius. Secondly, we assume that the line emission arises from a geometrically thin layer characterized by the distance from the disk midplane $z_{co}(r)=z_{0}(r/r_{0})^{q}$, where $r$ is the cylindrical radius, $r=\sqrt{R^{2}+z_{co}^{2}}$. Finally, we assume that the disk midplane is inclined with respect to the line of sight by an angle $i$ and is rotated in the sky by the position angle PA. The component of the Keplerian velocity along the line of sight is $v_{l}(R,\theta)=v_{k}(R)\sin{i}\sin{\theta}$, where $\theta$ is the azimuthal angle. The isovelocity contours corresponding to the velocity $v$ in a spherical reference frame centered at the center of the disk is calculated by finding the radii that satisfy the relation $v-v_{l}(R,\theta)=0$, for $\theta$ between 0 and 2$\pi$. Finally, the conversion between the reference frame centered on the disk and the image plane is

[TABLE]

where

[TABLE]