A Kinematical Detection of Two Embedded Jupiter-mass Planets in HD 163296

Despite the detection of close to 4000 fully formed planets, there are only a handful of planets detected during their formative stages. Characterizing the formation environment and witnessing the delivery of volatiles to the atmosphere of a young planet is an essential step in understanding the planet formation process. Toward this end, the Atacama Large Millimetre Array (ALMA) has revealed a variety of substructure in continuum emission from protoplanetary disks suggestive of planet–disk interactions sculpting the dust density distribution (ALMA Partnership et al. 2015; Andrews et al. 2016; Isella et al. 2016; Zhang et al. 2016).

Inferences of embedded planets and estimates of their mass are typically made from these observations by extrapolating a surface brightness profile of the thermal mm continuum to a gas surface density which, in turn, can be compared with either analytical prescriptions (e.g., Duffell 2015) or directly to hydrodynamic simulations (Fedele et al. 2018). Similar perturbations observed in molecular line emission have been used as additional constraints on the gas surface density profile (Isella et al. 2016; Fedele et al. 2017; Teague et al. 2017).

It is well known from attempts to measure the total disk gas mass that the extrapolation of continuum and molecular line emission to a gas surface density is fraught with uncertainty (Bergin et al. 2013; Miotello et al. 2017). Ill-constrained gas-to-dust ratios, radially varying grain opacities, complex grain evolution dictated by a coupling of the solid particles to the gas (Birnstiel et al. 2012), and local chemical and excitation effects (Öberg et al. 2015; Cleeves 2016) all conspire to limit the accuracy of the recovered surface density profile.

It has been shown that a gap is not a unique signature of a planet. Massive planets are able to excite spiral waves that open up secondary and tertiary gaps (Bae et al. 2017; Fedele et al. 2018), while grain growth around ice lines (Zhang et al. 2015) and the shepherding of dust by (magneto-)hydrodynamic instabilities have also been shown to produce ring-like structure in the continuum (Birnstiel et al. 2015; Flock et al. 2015; Okuzumi et al. 2016). In total, these effects can result in a significant uncertainty in the derived planet mass and, in some cases, lead to an incorrect inference of a planet.

An alternative approach is to use the rotation of the gas to probe the local pressure gradient. As the disk is in both radial and vertical hydrostatic equilibrium, the rotation velocity is given by

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\mathrm{rot}}^{2}}{r}=\displaystyle \frac{{{GM}}_{\star }r}{{({r}^{2}+{z}^{2})}^{3/2}}+\displaystyle \frac{1}{{\rho }_{\mathrm{gas}}}\displaystyle \frac{\partial P}{\partial r},\end{eqnarray} \tag{ 1 }$$

where M_⋆ is the mass of the star and ∂P/∂r is the radial pressure gradient (Rosenfeld et al. 2013). We have neglected the impact of self-gravity as this term requires knowledge of the disk mass. Rosenfeld et al. (2013) demonstrated that this term will introduce a slight hastening of the rotation at large radii for only the most massive disks, and will not introduce small-scale perturbations. When a planet perturbs the local gas density there will be a change in the local pressure gradient that manifests as a change in the rotation velocity as shown in Figure 1 (see also Kanagawa et al. 2015). As this method traces the pressure gradient the deviations in rotation velocity, δv_rot = (v_rot − v_kep)/v_kep, will inform us on the shape of the perturbation. By directly tracing the gas pressure, this technique is free from the uncertainties discussed above, and to which traditional methods are prone.

Zoom In Zoom Out Reset image size

In this Letter we present a new method to precisely measure v_rot, which is used to make the first kinematical detection of embedded planets in the disk of HD 163296, an A1 star at 101.5 pc (Bailer-Jones et al. 2018).⁵ It has been suggested that this source plays host to multiple planets due to the substructure observed in its continuum and CO line emission (Isella et al. 2016; Liu et al. 2018, hereafter I16 and L18, respectively).

In Section 2 we briefly describe the archival data used before presenting the method and results in Sections 3 and 4. We discuss the robustness of the method in Section 5 and provide a summary in Section 6.

We use archival data of the J = 2−1 transitions of ¹²CO, ¹³CO, and C¹⁸O from the disk HD 163296 (2013.1.00601.S; PI: A. Isella). These data were originally presented by I16, who reported dips in the emission profile coincident with ringed substructure in the mm continuum. Modeling of this emission suggested the presence of three planets of 0.46, 0.46, and 0.58 M_Jup at 59, 105, and 160 au, respectively (L18).

Our reduction closely follows the procedure outlined in I16. The data were initially calibrated with CASA v4.4.0 before transferring to CASA v4.7.2 for self-calibration and imaging. Phase self-calibration was performed on the continuum by combining the three continuum spectral windows before being applied to the whole data set. No amplitude self-calibration was performed. Continuum emission was subtracted from the line data using the task uvcontsub, then the lines imaged using Briggs weighting and a robust value of 0.5, resulting in a beam size of 026 × 018 for all three lines. The channel spacing was averaged down to 50 m s⁻¹ in order to improve the signal-to-noise ratio (S/N). We refer the reader to I16 for a discussion about the morphology of the emission.

As the emission arises from a vertically elevated region with a lower limit set by the freeze-out of the CO molecules (Rosenfeld et al. 2013; Schwarz et al. 2016), we first measure the emission height to (a) accurately deproject pixels into radial bins and (b) account for any contributions to $\delta {v}_{\mathrm{rot}}$ from changes in z (see Equation (1)). For this we follow the method presented in Pinte et al. (2018), using a position angle of 132° that we obtained from fitting a Keplerian rotation pattern to the first moment map and is consistent with previous determinations (Flaherty et al. 2015, 2017). The emission surface is modeled as a Gaussian Process (GP), which assumes that the observations are drawn from a smoothly varying function. This approach removes the need to assume an analytical form for the function, providing more flexibility in the model in addition to not requiring the binning of data (Foreman-Mackey et al. 2017). The resulting GP model for the C¹⁸O emission is presented in Figure 2(d), and shows slight depressions in the emission surface over the gap locations. This is indicative of a reduction in the local gas scale height or a perturbation in the CO column density.

Zoom In Zoom Out Reset image size

To derive v_rot, we exploit the azimuthal symmetry of the system and build upon the spectral deprojection and azimuthally averaging technique used in previous studies of protoplanetary and debris disks (Teague et al. 2016; Yen et al. 2016; Matrà et al. 2017). At a location (r, θ), where r is the radius and θ is polar angle relative to the major axis of the disk, the center of an emission line will be Doppler shifted by an amount v_rot · cosθ. Conversely, other properties of the line profile, such as amplitude and width, vary only a function of radius and are constant around θ. Thus, correcting for this offset allows for the lines to be stacked in order to provide a significant boost in the S/N. Rather than assuming v_rot a priori for this deprojection, we are able to infer its value from the data. Any error in the assumed v_rot value will result in slightly offset line profiles, resulting in a broadened stacked profile. We assert that the correct v_rot is the value that minimizes this width.

This approach was validated on a suite of forward models with known v_rot. Although the chosen bin width is below the spatial resolution of the data, wider bins were found to sample a sufficiently large range of ${v}_{\mathrm{kep}}$ that would overwhelm any signal from changes in the pressure gradient. For a S/N of ≳10, as with the observations, we could recover v_rot accurately to a precision of 2 m s⁻¹.

To apply this to the observations, we first binned the data into annuli roughly two pixels wide (∼9 au). To derive v_rot we used the L-BFGS-B method implemented in the scipy.optimize package. Panel (e) in Figure 2 shows the derived v_rot profile for C¹⁸O in blue. The uncertainties on individual bins are roughly 2 m s⁻¹ between 50 and 250 m s⁻¹, while the uncertainty on the GP models, shown in the top row of Figure 2, is larger at 8 m s⁻¹. We find an uncertainty on the GP model that is larger than the individual bins because the extreme flexibility of the GP model makes it harder to distinguish between small-scale fluctuations and noise, unlike when fitting an analytical function. In essence, precision in the resulting model is sacrificed to allow for a larger population of smooth models to be considered. Relative residuals from a Keplerian curve assuming 2.3 M_Sun and i = 477, appropriate for HD 163296 (Flaherty et al. 2015, 2017) are shown for the three transitions in the top row of Figure 2. Clear deviations are observed centered on the gaps in the continuum (shown by the vertical dotted lines; I16), consistent with the predicted profile shown in Figure 1. The significant deviation at smaller radii is an effect of the finite resolution of the observations. Testing with forward models show that this method is significantly biased within two beam FWHM from the disk center, equal to ≈52 au for the current observations.

In the following, we focus on the C¹⁸O emission as this arises from a region closer to the disk midplane so that 2D hydrodynamic models are sufficient to model the velocity perturbations. Conversely, ¹²CO and ¹³CO emission arises from much higher altitudes and thus requires the use of considerably more computationally expensive 3D hydrodynamic models. As such, the interpretation of these lines is left for future work focusing on the vertical structure of the perturbation.

Using Equation (1) we are able to quantify the relative importance of the height, gas temperature, and density. Although all three factors likely act in concert, isolating the main sources of perturbations helps in constraining the physical processes at play. Changes in the emission height can account for deviations in v_rot of ≲0.5% at 165 au; however, no significant change in z is observed at 100 au or closer in. If the deviations were solely due to temperature perturbations then we would require changes in temperature of 45% and 75% at 100 and 165 au, respectively. These correspond to changes in linewidth of 19% and 50%, and changes in the integrated intensity of 48% and 87%. Such large temperature changes can be ruled out as the linewidth that shows deviations of ≲10% across the gaps (Figure 2(f)) and the radial C¹⁸O intensity profile (Figure 1 of I16) that shows deviations of ≈50% and ≲25% at 100 and 165 au. This suggests that while changes in temperature are likely present, the dominant source of the perturbation is changes in the density profile, consistent with the previously proposed embedded planet scenario (I16; L18). Therefore, in the following we focus on changes in the density structure and leave the effects of temperature for future work.

To test the hypothesis of planetary causes for $\delta {v}_{\mathrm{rot}}$, we ran hydrodynamical simulations of embedded protoplanets using FARGO-3D (Benítez-Llambay & Masset 2016) including the orbital advection algorithm FARGO (Masset 2000). These used a polar grid with 1024 logarithmically spaced radial cells between 16 and 480 au, and 1920 uniformly spaced azimuthal cells covering 2π radians. The underlying physical structure was taken from the best-fit model of Flaherty et al. (2017).

We opted to run 2D simulations in order to efficiently explore the parameter space. An isothermal equation of state was assumed and a uniform viscosity of α = 10⁻³ was adopted, consistent with the constraints on the turbulence level in HD 163296 (Flaherty et al. 2015, 2017). The calculations run for 1000 orbital times for the outermost planet, corresponding to 4000 and 2100 for the innermost and middle planet, respectively. In all simulations, we found that the gas densities and velocities had reached steady state before the end of the simulation.

To compare these simulations to the observations we generate C¹⁸O image cubes, simulate the observations, and apply the methodology outline in Section 3 to the simulated observations. This approach will fully account for changes in the emission height (and subsequent changes in the temperature) that can arise from perturbed surface densities. To generate the C¹⁸O image cubes we inflated the azimuthally averaged surface density from the simulations assuming hydrostatic equilibrium. Above the 19 K isotherm and below a shielding column of $1.2\times {10}^{21}\,{{\rm{H}}}_{2}\,{\mathrm{cm}}^{-2}$ (Visser et al. 2009), a molecular abundance of 8.67 × 10⁻⁸ was assumed for C¹⁸O. Figure 3 shows the azimuthally averaged physical structure of the best-fit model and the resulting $\delta {v}_{\mathrm{rot}}$ profile, which is consistent with the observed profile in Figure 2. Figure 4 shows the simulations before azimuthal averaging, the density structure in (a), $\delta {v}_{\mathrm{rot}}$ in (b), and a zoom-in around the planet at 165 au.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Radiative transfer was performed with the non-LTE code LIME (Brinch & Hogerheijde 2010) with image properties matching the observations. As we do not expect significant spatial filtering from the data (I16), the images were convolved with a 2D Gaussian beam consistent with the C¹⁸O observations to provide a fair comparison. We find an excellent match to the observations with a 1 M_Jup planet at 100 au, and a 1.3 M_Jup planet at 165 au as shown in Figure 4; however, were unable to find a convincing fit to the feature at <70 au. A planet of 0.6 M_Jup at 65 au provided the best fit, but other potential scenarios are discussed in Section 5. We were able to constrain the mass of the two outer planets to within 50% and their locations to within 10%. Figure 5 compares the best-fit hydrodynamic model in blue with the observed $\delta {v}_{\mathrm{rot}}$ profile. The orange line shows the initial, unperturbed model demonstrating that a smooth model is a poor match to the data.

Zoom In Zoom Out Reset image size

The results derived from this method are dependent on relative flux measurements rather than absolute intensities, as in previous works. To demonstrate the robustness of this method to the model uncertainties, we present a model where we have increased the mass of the disk by a factor of five through a rescaling of the surface density (and thus an increase in C¹⁸O column density). This model is plotted in red in Figure 5. The derived $\delta {v}_{\mathrm{rot}}$ profile is minimally affected, other than a slight shift to lower values as the emission arises from a higher (and therefore slower rotating) region. Conclusions based on the presence and mass of the planets are therefore insensitive to the disk mass, thus these constraints are exceptionally complimentary to current methods.

For this pilot study we have not included any explicit heating or cooling across the gaps as this introduces additional uncertainty. The level of heating depends upon the accretion assumed onto the planet, which can vary by orders of magnitude for a specific planetary mass (Mordasini et al. 2017), while changes in the thermal structure due to changes in the dust distribution require computationally expensive modeling of the 2D grain distribution (Facchini et al. 2017; Teague et al. 2017). As we have not included cooling of the gas within the gap, our planet masses should be seen as an upper limit and the inclusion of a decrease in temperature in the gaps could account for the difference in derived masses between this work and that in L18.

Although in no way optimized to reproduce emission profiles, our best-fit model is able to broadly recover the observed deviations in the C¹⁸O normalized integrated intensity profiles (Figure 1, I16). Deriving a model consistent with both the velocity signatures and the radial emission profiles is left for future work. Such an approach will result in unparalleled constraints on the physical conditions around planet-opened gaps.

As with I16, no reasonable fit was found for the $\delta {v}_{\mathrm{rot}}$ perturbation at 80 au with a single planet; even the continuum ring is too wide to be well described by multiple embedded planets. An alternative solution for this perturbation is the pressure confinement of grains at 80 au due to the edge of the deadzone of the magneto-rotational instability (Flock et al. 2015; Pinilla et al. 2016). Such pressure confinement would require a pressure maximum at the center of the bright continuum ring at 80 au, and thus a local maximum in $\delta {v}_{\mathrm{rot}}$ slightly inward of this location. Recovery of the rotation profile is limited here by spatial resolution, and higher observations would be sufficient to see whether a local minimum in $\delta {v}_{\mathrm{rot}}$ could be resolved at ≈45 au, thereby distinguishing between these scenarios.

As we have unparalleled constraints on the gas surface density profile and the size and shape of the perturbations, we are able to rule out other scenarios that may produce similar features in the gas surface density. For example, the vertical shear instability has been shown to also produce concentric rings and gaps with a width that is much narrower than the gas scale height (Flock et al. 2017). This would result in gaps far narrower than those required to match the v_rot observations at 100 and 165 au. Similarly, ambipolar diffusion-driven self-organization has been shown to perturb the surface density (Béthune et al. 2017); however, this results in rings of constant width across the disk, which is again inconsistent with the observations. Finally, reconnection of magnetic field lines can produce ringed substructure (Suriano et al. 2018), however the strength of this in the outer disk, r ≳ 10 au, has yet to be demonstrated. For the outer two perturbations, only the embedded planet scenario is able to succinctly account for all of the observations. Nonetheless, without the detection of a point source, we cannot unambiguously dismiss scenarios without a planet.

We have presented a new method that enables the direct measurement of the gas pressure profile. This allows for significantly tighter and more accurate constraints on the gas surface density profile than traditional methods. Furthermore, as this method is sensitive to the gap profile, it provides essential information about the gap width in the gas, which is typically poorly constrained from brightness profiles.

Application of this method to CO isotopologue emission from HD 163296 revealed the predicted deviations for a significant perturbation in the local pressure. By accurately measuring the height of emission and using the linewidth as a proxy of local gas temperature, we were able to isolate changes in the local gas density as the primary drivers of these perturbations for C¹⁸O.

Comparisons with hydrodynamic models found excellent agreement with a 1 M_Jup planet at 100 au and a 1.3 M_Jup at 165 au (83 and 137 au at the new distance of 101.5 au) and allowed for a constraint on the mass and radius of 50% and 10%, respectively. A fit to the inner perturbation was less successful and requires high-resolution observations to distinguish between scenarios.

This method represents a new approach to searching for planets still embedded in their parental protoplanetary disk. By tracing the total gas component directly via the pressure, this method is free from the numerous uncertainties involved with more traditional approaches of mapping flux measurements to gas surface densities.

We would like to thank the referee for helpful comments leading to a much clearer presentation of the results. This paper makes use of the following ALMA data: JAO.ALMA#2013.1.00601.S. ALMA is a partnership of European Southern Observatory (ESO) (representing its member states), National Science Foundation (USA), and National Institutes of Natural Sciences (Japan), together with National Research Council (Canada), National Science Council and Academia Sinica Institute of Astronomy and Astrophysics (Taiwan), and Korea Astronomy and Space Science Institute (Korea), in cooperation with Chile. The Joint ALMA Observatory is operated by ESO, Associated Universities, Inc/National Radio Astronomy Observatory (NRAO), and National Astronomical Observatory of Japan. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. T.B. acknowledges funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program under grant agreement No. 714769. J.B. acknowledges support from NASA grant NNX17AE31G and computing resources provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. This work was supported by funding from NSF grants AST-1514670 and NASA NNX16AB48G.

Facility: ALMA. -

Software: emcee (Foreman-Mackey et al. 2013), celerite (Foreman-Mackey et al. 2017), CASA (McMullin et al. 2007), LIME (Brinch & Hogerheijde 2010), FARGO-3D (Benítez-Llambay & Masset 2016), scipy (Jones et al. 2001).